1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
/*
Input format:
N
Bx
Bz_0
Bz_1
...
Bz_N
J_{0,1}
J_{1,2}
...
J{n-2,n-1}
*/

#include <stdio.h>


const int N = 10;
const double Bx = 2;
const double Bz[] = {.3, -.9, -1.2, 1.1, 1.5, .6, -1.6, .3, -1.3, 1.9};
const double J[] = {-1.5, 1.3, -1.8, -1.3, -.6, 1.9, 1.1, -1.3, .4};

void CZ(qbit q1, qbit q2, double phi) {
    // printf("performing Controlled -2phi Z rotation\n");
    Rz(q2, -2.0*phi);
    CNOT(q1, q2);
    Rz(q2, phi);
    CNOT(q1, q2);
}

void ZcrossZ(qbit q1, qbit q2, double phi) {
    // printf("performing sigma_z cross sigma_z Hamiltonian\n");
    Rz(q1, phi);
    Rz(q2, -phi);
    CZ(q1, q2, -2.0*phi);
}


int main() {
    qbit reg[N];

    // Initialize all qubits to |+> state
    int n;
    for(n=0; n < N; n++) {
        PrepZ(reg[n], 0);
        H(reg[n]);
    }

    double T = 3;  // total duration of adiabatic evolution
    int M = 5;  // number of Trotter steps
    double dt = M / T;  // time per Trotter step

    int m;
    for(m = 1; m <= M; m++) {

        // Z-Z hamiltonian, executed in parallel
        for (n = 0; n < N-1; n += 2) {  // red pairs
            double phi = J[n] * (2.0*m - 1) / M;
            ZcrossZ(reg[n], reg[n + 1], phi);
        }
        for (n = 1; n < N-1; n += 2) {  // black pairs
            double phi = J[n] * (2.0*m - 1) / M;
            ZcrossZ(reg[n], reg[n + 1], phi);
        }

        double theta1 = (1.0 - (2.0*m-1)/M) * -2 * Bx * T / M;
        for (n = 0; n < N; n++) {
            H(reg[n]);
            Rz(reg[n], theta1);
            H(reg[n]);
        }

        // Bz hamiltonian
        for (n = 0; n < N; n++) {
            double theta2 = (1.0 - (2.0*m-1)/M) * -2 * Bz[n] * T / M;
            Rz(reg[n], theta2);
        }
    }


    for (n = 0; n < N; n++) {
        MeasZ(reg[0]);
    }

    return 0;
}

/*
int main() {
    // make sure to define these as const:

    
    FILE *fp;
    fp = fopen("input.txt", "r");
    int N = 0;
    double Bx = 0;
    fscanf(fp, "%d", &N);
    fscanf(fp, "%lf", &Bx);
    double Bz[N];
    int i;
    for (i = 0; i < N; i++) {
        fscanf(fp, "%lf", &Bz[i]);
    }
    double J[N-1];
    for (i = 0; i < N-1; i++) {
        fscanf(fp, "%lf", &J[i]);
    }
    
    
    IsingModel();
    return 0;
}
*/