goldilocks-controlled-reduce / python /qiskit_goldilocks_reduce.py

October 2026 main drop: mirror from GitHub

4f2752b verified 15 days ago

16.5 kB

# 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))

	# 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 = 264 - 232 + 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))