Shor's Algorithm: factoring 15 with quantum phase estimation

// Shor's Algorithm (N=15, a=7)
//
// Uses phase estimation to find the order r of 7 mod 15.
// where r is the order such that a^r ≡ 1 (mod N).
// From the measured phase (e.g., 1/4 or 3/4), we recover r = 4,
// which yields the factors 3 and 5.
//
// 010: 0.010 = 0.25 = 1/4; r=4.
// 110: 0.110 = 0.75 = 3/4; r=4.
// otherwise: reject and repeat.
//
// gcd(pow(a, r/2)-1, N) = gcd(pow(7, 4/2)-1, 15) = 3.
// gcd(pow(a, r/2)+1, N) = gcd(pow(7, 4/2)+1, 15) = 5.

OPENQASM 3.0;

gate x q { U(pi, 0, pi) q; }
gate h q { U(pi/2.0, 0, pi) q; }
gate cx c, t { ctrl @ U(pi, 0, pi) c, t; }
gate ccx c0, c1, t { ctrl(2) @ U(pi, 0, pi) c0, c1, t; }
gate cr(theta) c, t { ctrl @ U(0, 0, theta) c, t; }

def modexp(qubit[3] q, qubit[4] a) {
    // controlled-U^(2^0)
    cx q[0], a[1];
    cx q[0], a[2];

    // controlled-U^(2^1)
    cx        a[0], a[2];
    ccx q[1], a[2], a[0];
    cx        a[0], a[2];

    cx        a[3], a[1];
    ccx q[1], a[1], a[3];
    cx        a[3], a[1];
}

def inv_qft(qubit[3] q) {
    h q[2];

    cr(-pi/2) q[2], q[1];
    h q[1];

    cr(-pi/4) q[2], q[0];
    cr(-pi/2) q[1], q[0];
    h q[0];
}

// N=15, a=7
qubit[3] q;
qubit[4] a;
reset q;
reset a;

h q;
x a[3];

modexp(q, a);
inv_qft(q);

measure a;