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
  45   *k = n; // we will get it to be k=n-d by caling shrinkstar d times below
  46
  47   randomX = calloc(n*d,sizeof(float));
  48   randomXcopy = calloc(n*d,sizeof(float));
  49   randomXsingvals = calloc(d,sizeof(float));
  50   superb = calloc((d),sizeof(float));
  51   U = calloc((d*d),sizeof(float));
  52   VT = calloc((n*n),sizeof(float));
  53
  54   int info;
  55   int numsingvals = -1;
  56
  57   while(numsingvals != d) {
  58     for(i=0; i<d; i++) {
  59       for(j=0; j<n; j++) {
  60         randomX[i*n+j] = (float)(rand()%2);
  61         randomXcopy[i*n+j] = randomX[i*n+j];
  62         //printf("%1.0f ", randomX[i*n+j]);
  63       }
  64       //printf("\n");
  65     }
  66
  67     info = LAPACKE_sgesvd( LAPACK_ROW_MAJOR, 'N', 'N', d, n, randomXcopy, n, randomXsingvals, U, n, VT, n, superb );
  68
  69     numsingvals = 0;
  70
  71     for(i=0; i<d; i++) {
  72       if(fabs(randomXsingvals[i]) >= ZEROTHRESHOLD)
  73         numsingvals++;
  74     }
  75
  76     //printf("Number of singular values: %d\n", numsingvals);
  77   }
  78
  79   *G = calloc(n, sizeof(int*)); *GBar = calloc(n, sizeof(int*));
  80
  81   *h = calloc(n, sizeof(int));
  82
  83   for(i=0; i<n; i++) {
  84     (*G)[i] = calloc(n, sizeof(int));
  85     (*GBar)[i] = calloc(n, sizeof(int));
  86
  87     (*G)[i][i] = 1;
  88     (*GBar)[i][i] = 1;
  89   }
  90
  91   int* xi = calloc(n, sizeof(int));
  92   for(i=0; i<d; i++) {
  93     for(j=0; j<n; j++){
  94       xi[j] = (int)randomX[i*n+j];
  95       //printf("%1.0f ", randomX[i*n+j]);
  96     }
  97     //printf("\n");
  98     //printVector(xi, n);
  99     shrinkstar(n, k, *h, *G, *GBar, xi, 0);
 100   }
 101   free(xi);
 102   //printf("n-d=%d\n", n-d);
 103   //printf("k=%d\n", *k);
 104   //printMatrix(*G, n, n);
 105
 106   for(i=0; i<n; i++)
 107     (*h)[i] = rand()%2;
 108
 109   *Q = (rand()%8); // Q \in Z/8Z
 110
 111   *D = calloc(*k, sizeof(int));
 112   for(i=0; i<*k; i++)
 113     (*D)[i] = (rand()%4)*2; // D_a \in {0, 2, 4, 6}
 114
 115   *J = calloc(*k, sizeof(int*));
 116   for(i=0; i<*k; i++) {
 117     (*J)[i] = calloc(*k, sizeof(int));
 118     for(j=(i+1); j<*k; j++)
 119       (*J)[i][j] = (rand()%2)*4; // J_{a,b} \in {0, 4} for a!=b
 120     (*J)[i][i] = (2*(*D)[i])%8;
 121   }
 122   for(i=0; i<*k; i++) {
 123     for(j=(i+1); j<*k; j++) {
 124       (*J)[j][i] = (*J)[i][j];
 125     }
 126   }
 127
 128   free(randomX); free(randomXcopy); free(randomXsingvals); free(superb); free(U); free(VT);
 129   return;
 130
 131 }