--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <complex.h>
+#include <math.h>
+
+void deallocate_mem(complex double ***arr, int rows);
+void printMatrix(complex double** a, int rows, int cols);
+complex double** multMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2);
+complex double** outerMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2);
+complex double** transp(complex double **a, int rows, int cols);
+complex double** conjtransp(complex double **a, int rows, int cols);
+complex double trace(complex double **a, int rows, int cols);
+complex double** scalarmultMatrix(complex double scalar, complex double **a, int rows, int cols);
+complex double** addMatrix(complex double **A, complex double **B, int rows, int cols);
+int readPaulicoeffs(int* omega, int* alpha, int* beta, int* gamma, int* delta, int numqubits);
+
+// order of matrix elements is [row][column]!!!
+
+static complex double (*(I2[])) = { (double complex[]) {1.0+0.0*I, 0.0+0.0*I}, (double complex[]) {0.0+0.0*I, 1.0+0.0*I} };
+static complex double (*(X[])) = { (double complex[]) {0.0*I, 1.0+0.0*I}, (double complex[]) {1.0+0.0*I, 0.0*I} };
+static complex double (*(Y[])) = { (double complex[]) {0.0*I, 0.0-1.0*I}, (double complex[]) {0.0+1.0*I, 0.0*I} };
+static complex double (*(Z[])) = { (double complex[]) {1.0+0.0*I, 0.0+0.0*I}, (double complex[]) {0.0+0.0*I, -1.0+0.0*I} };
+
+int main()
+{
+
+ int i, j;
+
+ int N; // number of qubits
+ scanf("%d", &N);
+
+ int K; // number of T gate magic states (set to the first 'K' of the 'N' qubits -- the rest are set to the '0' computational basis state)
+ scanf("%d", &K);
+
+ complex double **IN; // N-qubit identity matrix
+ IN = I2;
+ for(i=1; i<N; i++) {
+ IN = outerMatrix(IN,I2,pow(2,i),pow(2,i),2,2);
+ }
+
+ complex double (*(psiT[])) = { (double complex[]) {1.0/sqrt(2.0)}, (double complex[]) {0.5*(1.0+1.0*I)}}; // T gate magic state
+ complex double (*(psi0[])) = { (double complex[]) {1.0+0.0*I}, (double complex[]) {0.0*I}}; // T gate magic state
+
+ complex double **psiN;
+ if(K > 0) {
+ psiN = psiT;
+ for(i=1; i<K; i++)
+ psiN = outerMatrix(psiN, psiT, pow(2,i), 1, 2, 1);
+ for(i=K; i<N; i++)
+ psiN = outerMatrix(psiN, psi0, pow(2,i), 1, 2, 1);
+ } else {
+ psiN = psi0;
+ for(i=1; i<N; i++)
+ psiN = outerMatrix(psiN, psi0, pow(2,i), 1, 2, 1);
+ }
+
+
+ int omega, alpha[N], beta[N], gamma[N], delta[N];
+
+ int Paulicounter = 0;
+
+ complex double **fullP; // full product: \prod_i 1/2*(1+P_i)
+ complex double **P; // P (P_i above) is made up of products of one-qubit Paulis, P1
+ complex double **P1[N]; // one-qubit Paulis
+
+ complex double tr;
+
+ printf("psiN:\n");
+ for(i=0; i<pow(2,N); i++) {
+ printf("%d: %lf+%lfI\n", i, creal(psiN[i][0]), cimag(psiN[i][0]));
+ }
+
+
+ while(readPaulicoeffs(&omega, alpha, beta, gamma, delta, N)) { // go over the product of 1/2*(I+Paulis) that makes up the full projector
+
+ Paulicounter++;
+ if(Paulicounter > N) {
+ printf("Error: Number of Paulis is greater than N!\n");
+ return 1;
+ }
+
+
+ printf("%d\n", omega);
+ for(i=0; i<N; i++) {
+ printf("%d %d %d %d\n", alpha[i], beta[i], gamma[i], delta[i]);
+ P1[i] = addMatrix(addMatrix(addMatrix(scalarmultMatrix(alpha[i],I2,2,2),scalarmultMatrix(beta[i],Z,2,2),2,2),scalarmultMatrix(gamma[i],X,2,2),2,2),scalarmultMatrix(delta[i],Y,2,2),2,2);
+ }
+
+ P = P1[0];
+ for(i=1; i<N; i++)
+ P = outerMatrix(P,P1[i],pow(2,i),pow(2,i),2,2);
+ P = scalarmultMatrix(cpow(I,omega),P,pow(2,N),pow(2,N));
+
+ if(Paulicounter == 1)
+ fullP = scalarmultMatrix(0.5,addMatrix(IN,P,pow(2,N),pow(2,N)),pow(2,N),pow(2,N));
+ else {
+ fullP = multMatrix(scalarmultMatrix(0.5,addMatrix(IN,P,pow(2,N),pow(2,N)),pow(2,N),pow(2,N)),fullP,pow(2,N),pow(2,N),pow(2,N),pow(2,N));
+ }
+ deallocate_mem(&P, pow(2,N));
+
+ }
+
+ complex double **psiNfinal = multMatrix(fullP,psiN,pow(2,N),pow(2,N),pow(2,N),1);
+
+ printf("fullP:\n");
+ for(i=0; i<pow(2,N); i++) {
+ for(j=0; j<pow(2,N); j++) {
+ printf("%lf+%lfI ", creal(fullP[i][j]), cimag(fullP[i][j]));
+ }
+ printf("\n");
+ }
+ printf("psiNfinal:\n");
+ for(i=0; i<pow(2,N); i++) {
+ printf("%d: %lf+%lfI\n", i, creal(psiNfinal[i][0]), cimag(psiNfinal[i][0]));
+ }
+
+ tr = 0.0 + 0.0*I;
+ //printf("tr:\n");
+ for(i=0; i<pow(2,N); i++) {
+ tr += conj(psiN[i][0])*psiNfinal[i][0];
+ //printf("%d: %lf+%lfI\n", i, creal(tr), cimag(tr));
+ }
+
+ if(creal(tr+0.00000001)>0)
+ printf("%.10lf %c %.10lf I\n", cabs(creal(tr)), cimag(tr+0.00000001)>0?'+':'-' , cabs(cimag(tr)));
+ else
+ printf("%.10lf %c %.10lf I\n", creal(tr), cimag(tr+0.00000001)>0?'+':'-' , cabs(cimag(tr)));
+ //printf("%lf %c %lf I\n", creal(tr), cimag(tr)>0?'+':'-' , cabs(cimag(tr)));
+ //printf("%lf\n", cabs(creal(tr)));
+ /* cabs the creal part because it's always positive, but sometimes the 0.0 gets a minus sign which is annoying to see when comparing outputs */
+
+ deallocate_mem(&psiNfinal, pow(2,N));
+
+
+ // deallocate_mem(&psiN, pow(2,N));
+
+ return 0;
+}
+
+
+complex double** addMatrix(complex double **A, complex double **B, int rows, int cols)
+{
+ int i, j;
+
+ complex double** C;
+
+ C = calloc(cols, sizeof(complex double*));
+ for(i=0; i<cols; i++)
+ C[i] = calloc(rows, sizeof(complex double));
+
+ for(i=0; i<rows; i++)
+ for(j=0; j<cols; j++)
+ C[i][j] = A[i][j] + B[i][j];
+
+ return C;
+}
+
+complex double** scalarmultMatrix(complex double scalar, complex double **a, int rows, int cols)
+{
+ int i, j;
+
+ complex double** C;
+
+ C = calloc(cols, sizeof(complex double*));
+ for(i=0; i<cols; i++)
+ C[i] = calloc(rows, sizeof(complex double));
+
+ for(i=0; i<rows; i++)
+ for(j=0; j<cols; j++)
+ C[i][j] = scalar*a[i][j];
+
+ return C;
+}
+
+complex double trace(complex double **a, int rows, int cols)
+{
+ int i;
+ complex double tr = 0.0*I;
+
+ for(i=0; i<rows; i++)
+ tr += a[i][i];
+
+ return tr;
+}
+
+complex double** transp(complex double **a, int rows, int cols)
+{
+ int i, j;
+ complex double** C;
+
+ C = calloc(cols, sizeof(complex double*));
+ for(i=0; i<cols; i++)
+ C[i] = calloc(rows, sizeof(complex double));
+
+ for(i=0; i<cols; i++)
+ for(j=0; j<rows; j++) {
+ C[i][j] = a[j][i];
+ }
+
+ return C;
+}
+
+complex double** conjtransp(complex double **a, int rows, int cols)
+{
+ int i, j;
+ complex double** C;
+
+ C = calloc(cols, sizeof(complex double*));
+ for(i=0; i<cols; i++)
+ C[i] = calloc(rows, sizeof(complex double));
+
+ for(i=0; i<cols; i++)
+ for(j=0; j<rows; j++) {
+ C[i][j] = conj(a[j][i]);
+ }
+
+ return C;
+}
+
+void printMatrix(complex double** a, int rows, int cols)
+{
+ int i, j;
+ printf("Matrix[%d][%d]\n", rows, cols);
+ for(i=0; i<rows; i++) {
+ for(j=0; j<cols; j++) {
+ printf("%lf+%lfI ", creal(a[i][j]), cimag(a[i][j]));
+ }
+ printf("\n");
+ }
+}
+
+complex double** multMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
+{
+ int i, j, k;
+ complex double **C;
+ C = calloc(ro1, sizeof(complex double*));
+ for(i=0; i<ro1; i++)
+ C[i] = calloc(co2, sizeof(complex double));
+
+ for(i=0; i<ro1; i++) {
+ for(j=0; j<co2; j++) {
+ C[i][j] = 0;
+ for(k=0; k<co1; k++)
+ C[i][j] += A[i][k] * B[k][j];
+ }
+ }
+
+ return C;
+}
+
+complex double** outerMatrix(complex double **A, complex double **B, int ro1, int co1, int ro2, int co2)
+{
+ int i, j, k, l;
+ complex double **C;
+ C = calloc(ro1*ro2, sizeof(complex double*));
+ for(i=0; i<ro1*ro2; i++)
+ C[i] = calloc(co1*co2, sizeof(complex double));
+
+ for(i=0; i<ro1; i++)
+ for(j=0; j<ro2; j++)
+ for(k=0; k<co1; k++)
+ for(l=0; l<co2; l++)
+ C[j+ro2*i][l+co2*k] = A[i][k]* B[j][l];
+
+ return C;
+}
+
+
+
+void deallocate_mem(complex double ***arr, int rows)
+{
+ int i;
+ for(i=0; i<rows; i++)
+ free((*arr)[i]);
+ free(*arr);
+}
+
+int readPaulicoeffs(int *omega, int *alpha, int *beta, int *gamma, int *delta, int numqubits)
+{
+
+ int i;
+
+ if(scanf("%d", omega) != EOF) {
+ for(i=0; i<numqubits; i++) {
+ if(scanf("%d %d %d %d", &alpha[i], &beta[i], &gamma[i], &delta[i]) == EOF) {
+ printf("Error: Too few input coeffs!\n");
+ exit(0);
+ }
+ if(alpha[i]+beta[i]+gamma[i]+delta[i] > 1) {
+ printf("Error: Too many coefficients are non-zero at Pauli %d!\n", i);
+ exit(0);
+ }
+ }
+ return 1;
+ } else
+ return 0;
+
+}
+
+
+
projects / weak_simulation_stab_extent.git / blobdiff