--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <string.h>
+#include <math.h>
+
+#include "matrix.h"
+
+#include "randomstabilizerstate_equatorial.h"
+
+#define ZEROTHRESHOLD (0.00000001)
+
+// Note: Make sure you seed the random number generator before calling this!
+// i.e. srand((unsigned)time(NULL));
+void randomstabilizerstate_equatorial(int n, int* k, int** h, int*** G, int*** GBar, int* Q, int** D, int*** J)
+{
+ // vector and matrix pointers should all be passed unallocated! (no freeing of memory performed)
+
+ int d;
+
+ int i, j;
+
+ d = 0;
+
+ *k = n; // we will get it to be k=n-d by calling shrinkstar d times below
+
+ int info;
+ int numsingvals = -1;
+
+ *G = calloc(n, sizeof(int*)); *GBar = calloc(n, sizeof(int*));
+ *h = calloc(n, sizeof(int));
+
+ for(i=0; i<n; i++) {
+ (*G)[i] = calloc(n, sizeof(int));
+ (*GBar)[i] = calloc(n, sizeof(int));
+
+ (*G)[i][i] = 1;
+ (*GBar)[i][i] = 1;
+ }
+
+ // keep the equatorial states with h=0 since it shouldn't matter if they have equal amplitude across all computational basis states
+ for(i=0; i<n; i++)
+ (*h)[i] = 0;//rand()%2;
+
+ *Q = (rand()%8); // Q \in Z/8Z
+
+ *D = calloc(*k, sizeof(int));
+ for(i=0; i<*k; i++)
+ (*D)[i] = (rand()%4)*2; // D_a \in {0, 2, 4, 6}
+
+ *J = calloc(*k, sizeof(int*));
+ for(i=0; i<*k; i++) {
+ (*J)[i] = calloc(*k, sizeof(int));
+ for(j=(i+1); j<*k; j++)
+ (*J)[i][j] = (rand()%2)*4; // J_{a,b} \in {0, 4} for a!=b
+ (*J)[i][i] = (2*(*D)[i])%8;
+ }
+ for(i=0; i<*k; i++) {
+ for(j=(i+1); j<*k; j++) {
+ (*J)[j][i] = (*J)[i][j];
+ }
+ }
+
+ return;
+
+}