randomstabilizerstate.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <time.h>
   4 #include <string.h>
   5 #include <math.h>
   6
   7 #include "lapacke.h"
   8 #include "matrix.h"
   9 #include "shrinkstar.h"
  10
  11 #include "randomstabilizerstate.h"
  12
  13 #define ZEROTHRESHOLD (0.00000001)
  14
  15 // Note: Make sure you seed the random number generator before calling this!
  16 // i.e. srand((unsigned)time(NULL));
  17 void randomstabilizerstate(int n, int* k, int** h, int*** G, int*** GBar, int* Q, int** D, int*** J, double** Pd)
  18 {
  19   // vector and matrix pointers should all be passed unallocated! (no freeing of memory performed)
  20
  21   float *randomX, *randomXcopy; // d*n matrix that is in array-form for LAPACKE
  22   float *randomXsingvals; // singular values of randomX
  23   float *U, *VT, *superb;
  24   int d;
  25
  26   int i, j;
  27
  28   //srand((unsigned)time(NULL));
  29
  30   double randPd = (double)rand()/(double)(RAND_MAX);
  31
  32   double cumprob = 0.0;
  33   for(i=0;i<n+1;i++) {
  34     cumprob += Pd[n-1][i]; // cumprob goes up to 1.0 once you count all Pd[n][0:n+1]
  35     if(cumprob > randPd) {
  36       d = i;
  37       break;
  38     }
  39     d=i;
  40   }
  41   //d = rand()%(n+1);
  42   //printf("d=%d total=%f randPd=%f\n", d, total, randPd);
  43   //printf("!!!d=%d\n", d);
  44   //d = 0;
  45
  46   printf("!!\n");
  47   fflush(stdout);
  48
  49   *k = n; // we will get it to be k=n-d by calling shrinkstar d times below
  50
  51   randomX = calloc(n*d,sizeof(float));
  52   randomXcopy = calloc(n*d,sizeof(float));
  53   randomXsingvals = calloc(d,sizeof(float));
  54   superb = calloc((d),sizeof(float));
  55   U = calloc((d*d),sizeof(float));
  56   VT = calloc((n*n),sizeof(float));
  57
  58   int info;
  59   int numsingvals = -1;
  60
  61   while(numsingvals != d) {
  62     for(i=0; i<d; i++) {
  63       for(j=0; j<n; j++) {
  64         randomX[i*n+j] = (float)(rand()%2);
  65         randomXcopy[i*n+j] = randomX[i*n+j];
  66         //printf("%1.0f ", randomX[i*n+j]);
  67       }
  68       //printf("\n");
  69     }
  70
  71     info = LAPACKE_sgesvd( LAPACK_ROW_MAJOR, 'N', 'N', d, n, randomXcopy, n, randomXsingvals, U, n, VT, n, superb );
  72
  73     numsingvals = 0;
  74
  75     for(i=0; i<d; i++) {
  76       if(fabs(randomXsingvals[i]) >= ZEROTHRESHOLD)
  77         numsingvals++;
  78     }
  79
  80     //printf("Number of singular values: %d\n", numsingvals);
  81   }
  82
  83   *G = calloc(n, sizeof(int*)); *GBar = calloc(n, sizeof(int*));
  84   *h = calloc(n, sizeof(int));
  85
  86   for(i=0; i<n; i++) {
  87     (*G)[i] = calloc(n, sizeof(int));
  88     (*GBar)[i] = calloc(n, sizeof(int));
  89
  90     (*G)[i][i] = 1;
  91     (*GBar)[i][i] = 1;
  92   }
  93
  94   int* xi = calloc(n, sizeof(int));
  95   for(i=0; i<d; i++) {
  96     for(j=0; j<n; j++){
  97       xi[j] = (int)randomX[i*n+j];
  98       //printf("%1.0f ", randomX[i*n+j]);
  99     }
 100     //printf("\n");
 101     //printVector(xi, n);
 102     shrinkstar(n, k, *h, *G, *GBar, xi, 0);
 103   }
 104   free(xi);
 105   //printf("n-d=%d\n", n-d);
 106   //printf("k=%d\n", *k);
 107   //printMatrix(*G, n, n);
 108
 109   for(i=0; i<n; i++)
 110     (*h)[i] = rand()%2;
 111
 112   *Q = (rand()%8); // Q \in Z/8Z
 113
 114   *D = calloc(*k, sizeof(int));
 115   for(i=0; i<*k; i++)
 116     (*D)[i] = (rand()%4)*2; // D_a \in {0, 2, 4, 6}
 117
 118   *J = calloc(*k, sizeof(int*));
 119   for(i=0; i<*k; i++) {
 120     (*J)[i] = calloc(*k, sizeof(int));
 121     for(j=(i+1); j<*k; j++)
 122       (*J)[i][j] = (rand()%2)*4; // J_{a,b} \in {0, 4} for a!=b
 123     (*J)[i][i] = (2*(*D)[i])%8;
 124   }
 125   for(i=0; i<*k; i++) {
 126     for(j=(i+1); j<*k; j++) {
 127       (*J)[j][i] = (*J)[i][j];
 128     }
 129   }
 130
 131   free(randomX); free(randomXcopy); free(randomXsingvals); free(superb); free(U); free(VT);
 132   return;
 133
 134 }