--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <complex.h>
+
+#include "matrix.h"
+#include "shrink.h"
+#include "extend.h"
+#include "exponentialsum.h"
+#include "innerproduct_equatorial.h"
+
+
+/****************************************************************************
+ * Note: Arguments are copied and not modified!
+ ****************************************************************************/
+
+double complex innerproduct_equatorial(int n1, int k1, int *h1, int **G1, int **GBar1, int Q1, int *D1, int **J1, int n2, int k2, int *h2, int **G2, int **GBar2, int Q2, int *D2, int **J2) {
+ // n1 should equal n2 (we assume this)
+
+ int a, b, c;
+
+ int n = n1;
+ int k = k1;
+ int *h = calloc(n, sizeof(int));
+ int **G = calloc(n, sizeof(int*));
+ int **GBar = calloc(n, sizeof(int*));
+ for(a=0; a<n; a++) {
+ h[a] = h1[a];
+ G[a] = calloc(n, sizeof(int));
+ GBar[a] = calloc(n, sizeof(int));
+ for(b=0; b<n; b++) {
+ G[a][b] = G1[a][b];
+ GBar[a][b] = GBar1[a][b];
+ }
+ }
+
+ int Q1copy = Q1;
+ int *D1copy = calloc(k1, sizeof(int));
+ int **J1copy = calloc(k1, sizeof(int*));
+ for(a=0; a<k1; a++) {
+ D1copy[a] = D1[a];
+ J1copy[a] = calloc(k1, sizeof(int));
+ for(b=0; b<k1; b++) {
+ J1copy[a][b] = J1[a][b];
+ }
+ }
+
+ int Q2copy = Q2;
+ int *D2copy = calloc(k2, sizeof(int));
+ int **J2copy = calloc(k2, sizeof(int*));
+ for(a=0; a<k2; a++) {
+ D2copy[a] = D2[a];
+ J2copy[a] = calloc(k2, sizeof(int));
+ for(b=0; b<k2; b++) {
+ J2copy[a][b] = J2[a][b];
+ }
+ }
+
+ int alpha;
+ int **R = calloc(k, sizeof(int*));
+ for(a=0; a<k; a++)
+ R[a] = calloc(k2, sizeof(int));
+ int **RT = calloc(k2, sizeof(int*));
+ for(a=0; a<k2; a++)
+ RT[a] = calloc(k, sizeof(int));
+ char shrinkResult;
+
+ /* for(a=k2; a<n2; a++) { */
+ /* alpha = dotProductMod(h2, GBar2[a], n2, 2); */
+
+ /* shrinkResult = shrink(n, &k, h, G, GBar, &Q1copy, &D1copy, &J1copy, GBar2[a], alpha); */
+ /* if(shrinkResult == 'E') {// if EMPTY */
+ /* free(h); */
+ /* deallocate_mem(&G, n); deallocate_mem(&GBar, n); */
+ /* free(D1copy); */
+ /* deallocate_mem(&J1copy, k); */
+ /* free(D2copy); */
+ /* deallocate_mem(&J2copy, k2); */
+ /* deallocate_mem(&R, k); */
+ /* deallocate_mem(&RT, k2); */
+ /* return(0.0+0.0*I); */
+ /* } */
+ /* } */
+
+ // Now K = K1 intersect K2 = (n1, k1, h1, G1, GBar1)
+ int *y = calloc(k2, sizeof(int));
+ int *tmpVec = calloc(n2, sizeof(int));
+ addVectorMod(h,h2,tmpVec,n2,2);
+ for(a=0; a<k2; a++) {
+ y[a] = dotProductMod(tmpVec,GBar2[a],n2,2);
+ for(b=0; b<k; b++) {
+ R[b][a] = dotProductMod(G[b],GBar2[a],n2,2);
+ }
+ }
+
+ // IF YOU UNCOMMENT THIS THEN YOU WILL BE CHANGING h2!!!
+ /* for(a=0; a<k2; a++) { */
+ /* for(b=0; b<n2; b++) { */
+ /* h2[b] = (h2[b] + y[a]*G2[a][b])%2; */
+ /* } */
+ /* } */
+ /* printf("h2 (should now match h):\n"); */
+ /* printVector(h2, n2); */
+ /* printf("h:\n"); */
+ /* printVector(h, n2); */
+
+ free(tmpVec);
+
+ // update Q2 and D2 using Eqs. 52 and 53 with y
+ for(a=0; a<k2; a++) {
+ Q2copy = Q2copy + D2copy[a]*y[a];
+ for(b=a+1; b<k2; b++) {
+ Q2copy = Q2copy + J2copy[a][b]*y[a]*y[b];
+ }
+ }
+ Q2copy = Q2copy%8;
+
+ for(a=0; a<k2; a++) {
+ for(b=0; b<k2; b++)
+ D2copy[a] = D2copy[a] + J2copy[a][b]*y[b];
+ D2copy[a] = D2copy[a]%8;
+ }
+
+ free(y);
+
+ // update D2copy and J2copy using Eqs. 49 and 50 with R
+ // (note that R is a kxk2 matrix so this may reduce the dimension of D2copy and J2copy)
+ y = calloc(k, sizeof(int));
+ for(c=0;c<k;c++) {
+ y[c] = 0;
+ for(a=0;a<k2;a++) {
+ y[c] = y[c] + R[c][a]*D2copy[a];
+ for(b=a+1;b<k2;b++)
+ y[c] = y[c] + J2copy[a][b]*R[c][a]*R[c][b];
+ }
+ y[c] = y[c]%8;
+ }
+ free(D2copy);
+ D2copy = calloc(k, sizeof(int));
+ memcpy(D2copy, y, sizeof(int)*k);
+ free(y);
+
+ transp(R,RT,k,k2);
+ int **tmpMatrix = calloc(k2, sizeof(int*));
+ for(a=0; a<k2; a++)
+ tmpMatrix[a] = calloc(k, sizeof(int));
+ multMatrixMod(J2copy,RT,tmpMatrix,k2,k2,k2,k,8);
+ deallocate_mem(&J2copy, k2);
+ J2copy = calloc(k, sizeof(int*));
+ for(a=0;a<k;a++)
+ J2copy[a] = calloc(k, sizeof(int));
+ multMatrixMod(R,tmpMatrix,J2copy,k,k2,k2,k,8);
+
+ deallocate_mem(&tmpMatrix, k2);
+
+ deallocate_mem(&R, k);
+ deallocate_mem(&RT, k2);
+
+ // Now (Q1, D1, J1) and (Q2, D2, J2) are defined in the same basis
+
+ int Q = (Q1copy - Q2copy)%8 < 0 ? (Q1copy - Q2copy)%8 + 8 : (Q1copy - Q2copy)%8;
+ int *D = calloc(k, sizeof(int));
+ int **J = calloc(k, sizeof(int*));
+
+ for(a=0; a<k; a++) { // again, if k=0 then no need to fill these
+ J[a] = calloc(k, sizeof(int));
+ D[a] = (D1copy[a] - D2copy[a])%8 < 0 ? (D1copy[a] - D2copy[a])%8 + 8 : (D1copy[a] - D2copy[a])%8;
+ for(b=0; b<k; b++)
+ J[a][b] = (J1copy[a][b] - J2copy[a][b])%8 < 0 ? (J1copy[a][b] - J2copy[a][b])%8 + 8 : (J1copy[a][b] - J2copy[a][b])%8;
+
+ }
+
+ /* free(h); */
+ /* deallocate_mem(&G, n); deallocate_mem(&GBar, n); */
+
+ double complex amplitude = pow(2,-0.5*(k1+k2))*exponentialsum(&k, &Q, D, J);
+
+ free(D);
+ deallocate_mem(&J, k);
+
+ free(D1copy);
+ deallocate_mem(&J1copy, k);
+
+ free(D2copy);
+ deallocate_mem(&J2copy, k);
+
+ return(amplitude);
+
+}
projects / weak_simulation_stab_extent.git / blobdiff