X-Git-Url: https://s3miclassical.com/gitweb/?p=weak_simulation_stab_extent.git;a=blobdiff_plain;f=innerproduct_equatorial.c;fp=innerproduct_equatorial.c;h=37f1e2c532472f23fbb1548da45c2ae385c5577a;hp=0000000000000000000000000000000000000000;hb=7bc25656311c391af6a4d5814cc82080981c63c4;hpb=d914e7c183347b1340988e0869a6e053dd573991

diff --git a/innerproduct_equatorial.c b/innerproduct_equatorial.c
new file mode 100644
index 0000000..37f1e2c
--- /dev/null
+++ b/innerproduct_equatorial.c
@@ -0,0 +1,190 @@
+#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);
+    
+}