python/qiskit_goldilocks_reduce.py

# qiskit_goldilocks_reduce.py
# Validation of Goldilocks Controlled Reduction Circuit (2^32 - 1)
# Requires: qiskit, qiskit-aer, numpy

import numpy as np
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit_aer import AerSimulator
from qiskit.compiler import transpile
from qiskit.visualization import plot_histogram
import matplotlib.pyplot as plt

# ---------------------------------------------------------
# 1. Circuit Construction (Matching OpenQASM Decomposition)
# ---------------------------------------------------------

def lower_step(qc: QuantumCircuit, x: int, c_in: int, c_out: int):
    """Gate sequence for lower 32 bits (Constant=1): x += 1 + c_in"""
    qc.x(x)                    # Constant 1
    qc.cx(c_in, x)             # Add carry
    qc.x(x)                    # Restore x_orig for Toffoli control
    qc.ccx(x, c_in, c_out)     # c_out = x_orig & c_in
    qc.x(x)                    # Restore x (now holds Sum)

def upper_step(qc: QuantumCircuit, x: int, c_in: int, c_out: int):
    """Gate sequence for upper 32 bits (Constant=0): x += 0 + c_in"""
    qc.ccx(x, c_in, c_out)     # c_out = x_orig & c_in (Compute carry FIRST)
    qc.cx(c_in, x)             # Sum = x_orig ^ c_in

def lower_step_dg(qc: QuantumCircuit, x: int, c_in: int, c_out: int):
    """Inverse of lower_step"""
    qc.x(x)                    # x_orig = NOT(x_current)
    qc.ccx(x, c_in, c_out)     # Uncompute Carry
    qc.x(x)                    # Restore x (Sum)
    qc.cx(c_in, x)             # Uncompute Sum (x ^= c_in)
    qc.x(x)                    # Uncompute Constant (x = NOT(x) -> Original)

def upper_step_dg(qc: QuantumCircuit, x: int, c_in: int, c_out: int):
    """Inverse of upper_step"""
    qc.cx(c_in, x)             # Uncompute Sum
    qc.ccx(x, c_in, c_out)     # Uncompute Carry

def build_goldilocks_reduce_circuit() -> QuantumCircuit:
    """
    Builds the full 68-qubit circuit:
    - 64 qubits for register x
    - 1 qubit for control (ctrl)
    - 2 qubits for ancillae (a0, a1)
    - 64 classical bits for measurement of x
    - 1 classical bit for measurement of ctrl (verification)
    - 2 classical bits for measurement of ancillae (cleanliness check)
    """
    qr_x = QuantumRegister(64, 'x')
    qr_ctrl = QuantumRegister(1, 'ctrl')
    qr_a = QuantumRegister(2, 'anc') # a0, a1
    cr_x = ClassicalRegister(64, 'meas_x')
    cr_ctrl = ClassicalRegister(1, 'meas_ctrl')
    cr_anc = ClassicalRegister(2, 'meas_anc')
    
    qc = QuantumCircuit(qr_x, qr_ctrl, qr_a, cr_x, cr_ctrl, cr_anc)
    
    # Indices
    x = [qr_x[i] for i in range(64)]
    ctrl = qr_ctrl[0]
    a0 = qr_a[0]
    a1 = qr_a[1]
    
    # ---------------------------------------------------------
    # FORWARD PASS
    # ---------------------------------------------------------
    
    # 1. LOAD CONTROL
    qc.cx(ctrl, a0)
    
    # 2. LOWER 32 BITS (Ping-Pong a0 <-> a1)
    # Bit 0: a0 -> a1
    lower_step(qc, x[0], a0, a1)
    # Bit 1: a1 -> a0
    lower_step(qc, x[1], a1, a0)
    # Bit 2: a0 -> a1
    lower_step(qc, x[2], a0, a1)
    # Bit 3: a1 -> a0
    lower_step(qc, x[3], a1, a0)
    # Bit 4: a0 -> a1
    lower_step(qc, x[4], a0, a1)
    # Bit 5: a1 -> a0
    lower_step(qc, x[5], a1, a0)
    # Bit 6: a0 -> a1
    lower_step(qc, x[6], a0, a1)
    # Bit 7: a1 -> a0
    lower_step(qc, x[7], a1, a0)
    # Bit 8: a0 -> a1
    lower_step(qc, x[8], a0, a1)
    # Bit 9: a1 -> a0
    lower_step(qc, x[9], a1, a0)
    # Bit 10: a0 -> a1
    lower_step(qc, x[10], a0, a1)
    # Bit 11: a1 -> a0
    lower_step(qc, x[11], a1, a0)
    # Bit 12: a0 -> a1
    lower_step(qc, x[12], a0, a1)
    # Bit 13: a1 -> a0
    lower_step(qc, x[13], a1, a0)
    # Bit 14: a0 -> a1
    lower_step(qc, x[14], a0, a1)
    # Bit 15: a1 -> a0
    lower_step(qc, x[15], a1, a0)
    # Bit 16: a0 -> a1
    lower_step(qc, x[16], a0, a1)
    # Bit 17: a1 -> a0
    lower_step(qc, x[17], a1, a0)
    # Bit 18: a0 -> a1
    lower_step(qc, x[18], a0, a1)
    # Bit 19: a1 -> a0
    lower_step(qc, x[19], a1, a0)
    # Bit 20: a0 -> a1
    lower_step(qc, x[20], a0, a1)
    # Bit 21: a1 -> a0
    lower_step(qc, x[21], a1, a0)
    # Bit 22: a0 -> a1
    lower_step(qc, x[22], a0, a1)
    # Bit 23: a1 -> a0
    lower_step(qc, x[23], a1, a0)
    # Bit 24: a0 -> a1
    lower_step(qc, x[24], a0, a1)
    # Bit 25: a1 -> a0
    lower_step(qc, x[25], a1, a0)
    # Bit 26: a0 -> a1
    lower_step(qc, x[26], a0, a1)
    # Bit 27: a1 -> a0
    lower_step(qc, x[27], a1, a0)
    # Bit 28: a0 -> a1
    lower_step(qc, x[28], a0, a1)
    # Bit 29: a1 -> a0
    lower_step(qc, x[29], a1, a0)
    # Bit 30: a0 -> a1
    lower_step(qc, x[30], a0, a1)
    # Bit 31: a1 -> a0 (Final Lower Carry in a0)
    lower_step(qc, x[31], a1, a0)
    
    # 3. UPPER 32 BITS (Ping-Pong a0 <-> a1)
    # Bit 32: a0 -> a1
    upper_step(qc, x[32], a0, a1)
    # Bit 33: a1 -> a0
    upper_step(qc, x[33], a1, a0)
    # Bit 34: a0 -> a1
    upper_step(qc, x[34], a0, a1)
    # Bit 35: a1 -> a0
    upper_step(qc, x[35], a1, a0)
    # Bit 36: a0 -> a1
    upper_step(qc, x[36], a0, a1)
    # Bit 37: a1 -> a0
    upper_step(qc, x[37], a1, a0)
    # Bit 38: a0 -> a1
    upper_step(qc, x[38], a0, a1)
    # Bit 39: a1 -> a0
    upper_step(qc, x[39], a1, a0)
    # Bit 40: a0 -> a1
    upper_step(qc, x[40], a0, a1)
    # Bit 41: a1 -> a0
    upper_step(qc, x[41], a1, a0)
    # Bit 42: a0 -> a1
    upper_step(qc, x[42], a0, a1)
    # Bit 43: a1 -> a0
    upper_step(qc, x[43], a1, a0)
    # Bit 44: a0 -> a1
    upper_step(qc, x[44], a0, a1)
    # Bit 45: a1 -> a0
    upper_step(qc, x[45], a1, a0)
    # Bit 46: a0 -> a1
    upper_step(qc, x[46], a0, a1)
    # Bit 47: a1 -> a0
    upper_step(qc, x[47], a1, a0)
    # Bit 48: a0 -> a1
    upper_step(qc, x[48], a0, a1)
    # Bit 49: a1 -> a0
    upper_step(qc, x[49], a1, a0)
    # Bit 50: a0 -> a1
    upper_step(qc, x[50], a0, a1)
    # Bit 51: a1 -> a0
    upper_step(qc, x[51], a1, a0)
    # Bit 52: a0 -> a1
    upper_step(qc, x[52], a0, a1)
    # Bit 53: a1 -> a0
    upper_step(qc, x[53], a1, a0)
    # Bit 54: a0 -> a1
    upper_step(qc, x[54], a0, a1)
    # Bit 55: a1 -> a0
    upper_step(qc, x[55], a1, a0)
    # Bit 56: a0 -> a1
    upper_step(qc, x[56], a0, a1)
    # Bit 57: a1 -> a0
    upper_step(qc, x[57], a1, a0)
    # Bit 58: a0 -> a1
    upper_step(qc, x[58], a0, a1)
    # Bit 59: a1 -> a0
    upper_step(qc, x[59], a1, a0)
    # Bit 60: a0 -> a1
    upper_step(qc, x[60], a0, a1)
    # Bit 61: a1 -> a0
    upper_step(qc, x[61], a1, a0)
    # Bit 62: a0 -> a1
    upper_step(qc, x[62], a0, a1)
    # Bit 63: a1 -> a0 (Final Carry Out in a0)
    upper_step(qc, x[63], a1, a0)
    
    # ---------------------------------------------------------
    # REVERSE PASS (Uncompute Carries)
    # ---------------------------------------------------------
    
    # REVERSE UPPER 32 (Bits 63 down to 32)
    upper_step_dg(qc, x[63], a0, a1)
    upper_step_dg(qc, x[62], a1, a0)
    upper_step_dg(qc, x[61], a0, a1)
    upper_step_dg(qc, x[60], a1, a0)
    upper_step_dg(qc, x[59], a0, a1)
    upper_step_dg(qc, x[58], a1, a0)
    upper_step_dg(qc, x[57], a0, a1)
    upper_step_dg(qc, x[56], a1, a0)
    upper_step_dg(qc, x[55], a0, a1)
    upper_step_dg(qc, x[54], a1, a0)
    upper_step_dg(qc, x[53], a0, a1)
    upper_step_dg(qc, x[52], a1, a0)
    upper_step_dg(qc, x[51], a0, a1)
    upper_step_dg(qc, x[50], a1, a0)
    upper_step_dg(qc, x[49], a0, a1)
    upper_step_dg(qc, x[48], a1, a0)
    upper_step_dg(qc, x[47], a0, a1)
    upper_step_dg(qc, x[46], a1, a0)
    upper_step_dg(qc, x[45], a0, a1)
    upper_step_dg(qc, x[44], a1, a0)
    upper_step_dg(qc, x[43], a0, a1)
    upper_step_dg(qc, x[42], a1, a0)
    upper_step_dg(qc, x[41], a0, a1)
    upper_step_dg(qc, x[40], a1, a0)
    upper_step_dg(qc, x[39], a0, a1)
    upper_step_dg(qc, x[38], a1, a0)
    upper_step_dg(qc, x[37], a0, a1)
    upper_step_dg(qc, x[36], a1, a0)
    upper_step_dg(qc, x[35], a0, a1)
    upper_step_dg(qc, x[34], a1, a0)
    upper_step_dg(qc, x[33], a0, a1)
    upper_step_dg(qc, x[32], a1, a0)
    
    # REVERSE LOWER 32 (Bits 31 down to 0)
    lower_step_dg(qc, x[31], a0, a1)
    lower_step_dg(qc, x[30], a1, a0)
    lower_step_dg(qc, x[29], a0, a1)
    lower_step_dg(qc, x[28], a1, a0)
    lower_step_dg(qc, x[27], a0, a1)
    lower_step_dg(qc, x[26], a1, a0)
    lower_step_dg(qc, x[25], a0, a1)
    lower_step_dg(qc, x[24], a1, a0)
    lower_step_dg(qc, x[23], a0, a1)
    lower_step_dg(qc, x[22], a1, a0)
    lower_step_dg(qc, x[21], a0, a1)
    lower_step_dg(qc, x[20], a1, a0)
    lower_step_dg(qc, x[19], a0, a1)
    lower_step_dg(qc, x[18], a1, a0)
    lower_step_dg(qc, x[17], a0, a1)
    lower_step_dg(qc, x[16], a1, a0)
    lower_step_dg(qc, x[15], a0, a1)
    lower_step_dg(qc, x[14], a1, a0)
    lower_step_dg(qc, x[13], a0, a1)
    lower_step_dg(qc, x[12], a1, a0)
    lower_step_dg(qc, x[11], a0, a1)
    lower_step_dg(qc, x[10], a1, a0)
    lower_step_dg(qc, x[9], a0, a1)
    lower_step_dg(qc, x[8], a1, a0)
    lower_step_dg(qc, x[7], a0, a1)
    lower_step_dg(qc, x[6], a1, a0)
    lower_step_dg(qc, x[5], a0, a1)
    lower_step_dg(qc, x[4], a1, a0)
    lower_step_dg(qc, x[3], a0, a1)
    lower_step_dg(qc, x[2], a1, a0)
    lower_step_dg(qc, x[1], a0, a1)
    lower_step_dg(qc, x[0], a1, a0)
    
    # 5. UNLOAD CONTROL
    qc.cx(ctrl, a0)
    
    return qc

# ---------------------------------------------------------
# 2. Classical Reference Implementation
# ---------------------------------------------------------

GOLDILOCKS_P = 2**64 - 2**32 + 1
CONST_C = 2**32 - 1

def classical_reduce(x_val: int, ctrl_val: int) -> int:
    """Reference: x + ctrl * (2^32 - 1) mod 2^64"""
    if ctrl_val == 1:
        return (x_val + CONST_C) & ((1 << 64) - 1)
    return x_val

def int_to_bits(val: int, n: int) -> list:
    """Little Endian bit list"""
    return [(val >> i) & 1 for i in range(n)]

def bits_to_int(bits: list) -> int:
    """Little Endian bits to int"""
    return sum(b << i for i, b in enumerate(bits))

# ---------------------------------------------------------
# 3. Test Harness
# ---------------------------------------------------------

def run_test_case(qc_base: QuantumCircuit, x_init: int, ctrl_init: int, shots: int = 1000):
    """Runs a single test case on AerSimulator."""
    qc = qc_base.copy()
    
    # Initialize Registers
    # Qiskit uses Big Endian for `initialize` but our circuit is Little Endian (q[0] = LSB).
    # We must initialize qubits individually or reverse the bitstring.
    # Easiest: Apply X gates for |1> states.
    
    x_bits = int_to_bits(x_init, 64)
    for i, bit in enumerate(x_bits):
        if bit == 1:
            qc.x(qc.qubits[i]) # qr_x[i]
            
    if ctrl_init == 1:
        qc.x(qc.qubits[64]) # qr_ctrl[0]
        
    # Ancillae start at |0> (default)
    
    # Measure
    qc.measure(qc.qubits[0:64], qc.clbits[0:64])      # x
    qc.measure(qc.qubits[64], qc.clbits[64])           # ctrl
    qc.measure(qc.qubits[65:67], qc.clbits[65:67])     # ancillae
    
    # Execute
    backend = AerSimulator(method='statevector') # Statevector is exact for verification
    # backend = AerSimulator() # For shot noise simulation
    
    tqc = transpile(qc, backend)
    job = backend.run(tqc, shots=shots)
    result = job.result()
    counts = result.get_counts()
    
    # Parse Results (Qiskit returns Big Endian hex strings typically, 
    # but with explicit measure mapping, it respects clbit order).
    # `counts` keys are strings like '00 0 00' (cr_x cr_ctrl cr_anc) separated by spaces if registers named.
    # Since we used separate ClassicalRegisters, the key format is 'cr_anc cr_ctrl cr_x' (reversed reg order) 
    # or similar. Safest: use `memory` or parse carefully.
    
    # Let's use `memory` for single-shot deterministic check (statevector)
    if shots == 1:
        mem = result.get_memory()
        # Memory string format matches clbit index order: [cr_x[0]...cr_x[63], cr_ctrl[0], cr_anc[0], cr_anc[1]]
        # Wait, Qiskit memory string is usually MSB first (clbit N-1 ... clbit 0).
        # We measure clbits 0..63 (x), 64 (ctrl), 65,66 (anc).
        # Memory string: bit 66 ... bit 0.
        # This is confusing. Let's just sample counts.
        pass
    
    # Parse Counts (Hex string usually)
    # Key format: 'anc_ctrl_x' (Register order reversed in string representation often)
    # We will just check the most frequent outcome.
    most_likely = max(counts, key=counts.get)
    
    # The bitstring in counts keys is ordered by *register declaration order* reversed? 
    # Actually: `qc.measure(qr_x, cr_x)` maps qr_x[i] -> cr_x[i].
    # `counts` keys are bitstrings where leftmost is cr_x[63] ... cr_x[0] cr_ctrl[0] cr_anc[1] cr_anc[0]?
    # Standard Qiskit: `counts` keys are big-endian integers representing the classical register state.
    # If we have registers [cr_x, cr_ctrl, cr_anc], the integer is (cr_anc << 65) | (cr_ctrl << 64) | cr_x.
    # The hex string represents this integer.
    
    # Let's decode the integer key.
    # We expect: x_final, ctrl_final, anc_final.
    
    # Since we only ran 1 shot (statevector), there is only 1 key.
    key = list(counts.keys())[0]
    # Remove spaces
    key = key.replace(" ", "")
    
    # Total clbits = 67. Hex length = ceil(67/4) = 17 chars.
    # Parse as int
    try:
        val = int(key, 16)
    except ValueError:
        # Handle binary string format if hex fails
        val = int(key, 2)
        
    # Extract fields (Little Endian clbit index -> Bit position in integer)
    # clbit 0 -> LSB.
    meas_x = val & ((1 << 64) - 1)
    meas_ctrl = (val >> 64) & 1
    meas_anc = (val >> 65) & 3
    
    expected_x = classical_reduce(x_init, ctrl_init)
    
    passed = (meas_x == expected_x) and (meas_ctrl == ctrl_init) and (meas_anc == 0)
    
    return {
        "passed": passed,
        "x_init": x_init,
        "ctrl_init": ctrl_init,
        "expected_x": expected_x,
        "measured_x": meas_x,
        "measured_ctrl": meas_ctrl,
        "measured_anc": meas_anc,
        "counts": counts
    }

# ---------------------------------------------------------
# 4. Main Execution
# ---------------------------------------------------------

if __name__ == "__main__":
    print("=" * 60)
    print("Goldilocks Controlled Reduction Circuit Validation")
    print("=" * 60)
    print(f"Prime p = 2^64 - 2^32 + 1 = {GOLDILOCKS_P}")
    print(f"Constant C = 2^32 - 1 = {CONST_C}")
    print(f"Circuit Qubits: 64 (x) + 1 (ctrl) + 2 (anc) = 67")
    print("Building circuit...")
    
    base_circuit = build_goldilocks_reduce_circuit()
    print(f"Circuit Depth: {base_circuit.depth()}")
    print(f"Circuit Size: {base_circuit.size()}")
    print(f"Toffoli Count: {base_circuit.count_ops().get('ccx', 0)}")
    print(f"CNOT Count: {base_circuit.count_ops().get('cx', 0)}")
    print(f"X Count: {base_circuit.count_ops().get('x', 0)}")
    
    # Test Vectors
    test_cases = [
        (0, 0, "Zero Input, Ctrl=0"),
        (0, 1, "Zero Input, Ctrl=1 (Add Const)"),
        (1, 0, "One Input, Ctrl=0"),
        (1, 1, "One Input, Ctrl=1"),
        (2**63, 1, "MSB Set, Ctrl=1 (Carry Propagation)"),
        (2**64 - 1, 1, "Max U64, Ctrl=1 (Wrap)"),
        (2**64 - CONST_C, 1, "Boundary: x + C = 2^64 -> 0"),
        (GOLDILOCKS_P - 1, 1, "Field Max (p-1), Ctrl=1 -> p-1 + C mod 2^64"),
        (0xFFFFFFFF00000000, 1, "Upper 32 bits set"),
        (0x00000000FFFFFFFF, 1, "Lower 32 bits set (Const)"),
    ]
    
    all_passed = True
    for x_val, ctrl_val, desc in test_cases:
        print(f"\nTest: {desc}")
        print(f"  Input:  x={x_val:#018x}, ctrl={ctrl_val}")
        
        res = run_test_case(base_circuit, x_val, ctrl_val, shots=1) # Statevector = 1 shot exact
        
        status = "PASS" if res["passed"] else "FAIL"
        print(f"  Output: x={res['measured_x']:#018x}, ctrl={res['measured_ctrl']}, anc={res['measured_anc']} [{status}]")
        print(f"  Expected x={res['expected_x']:#018x}")
        
        if not res["passed"]:
            all_passed = False
            print(f"  DETAILS: {res}")
            
    print("\n" + "=" * 60)
    if all_passed:
        print("ALL TESTS PASSED: Circuit correctly implements Goldilocks Reduction.")
        print("Ancillae cleanly returned to |0>. Control preserved.")
    else:
        print("SOME TESTS FAILED.")
    print("=" * 60)
    
    # Optional: Draw circuit (first 20 layers)
    # print(base_circuit.draw(output='text', fold=100, idle_wires=False))