Quantum Algorithms

"""
Quantum algorithm implementations using Qiskit and NumPy.
"""

import numpy as np
from math import gcd, pi, sqrt, floor

# ================================================================
# Manual state-vector simulation helpers
# ================================================================

def tensor(*ops):
    """Tensor product of matrices."""
    result = ops[0]
    for op in ops[1:]:
        result = np.kron(result, op)
    return result

H = np.array([[1, 1], [1, -1]]) / sqrt(2)
I = np.eye(2)
X = np.array([[0, 1], [1, 0]])
Z = np.array([[1, 0], [0, -1]])

def R(k):
    """Phase rotation gate R_k."""
    return np.array([[1, 0], [0, np.exp(2j * pi / 2**k)]])

# ================================================================
# Deutsch-Jozsa (statevector simulation)
# ================================================================

def deutsch_jozsa_sim(f_bitstring):
    """
    f_bitstring: list of 2^n values (0 or 1), defines f on {0,1}^n.
    Returns 'constant' or 'balanced'.
    """
    n = int(np.log2(len(f_bitstring)))
    N = 2**n
    # Build oracle matrix U_f for phase kickback
    U_f = np.eye(2 * N, dtype=complex)
    for x in range(N):
        if f_bitstring[x] == 1:
            U_f[2*x, 2*x]     =  0;  U_f[2*x, 2*x+1]   = 1
            U_f[2*x+1, 2*x]   =  1;  U_f[2*x+1, 2*x+1] = 0
    # Initial state |0...0>|1>
    state = np.zeros(2 * N, dtype=complex)
    state[1] = 1.0  # last qubit = |1>
    # Apply H to all n+1 qubits
    H_all = tensor(*([H] * (n + 1)))
    state = H_all @ state
    # Apply oracle
    state = U_f @ state
    # Apply H to first n qubits (acting on 2N dimensional space)
    H_n = tensor(*([H] * n), I)
    state = H_n @ state
    # Measure first n qubits (project onto |0...0> in first register)
    # Sum probabilities of |0...0>|0> and |0...0>|1>
    prob_zero = abs(state[0])**2 + abs(state[1])**2
    return 'constant' if prob_zero > 0.99 else 'balanced'

# Test: constant function
f_const = [0] * 8
print(deutsch_jozsa_sim(f_const))  # 'constant'

# Test: balanced function (first half 0, second half 1)
f_bal = [0]*4 + [1]*4
print(deutsch_jozsa_sim(f_bal))    # 'balanced'

# ================================================================
# Grover's search (statevector simulation)
# ================================================================

def grover_search(target, n):
    """
    Find 'target' in database of size 2^n using Grover's algorithm.
    Returns found index.
    """
    N = 2**n
    # Uniform superposition
    state = np.ones(N, dtype=complex) / sqrt(N)
    k_opt = floor(pi * sqrt(N) / 4)

    for _ in range(k_opt):
        # Oracle: flip sign of target
        state[target] *= -1
        # Diffusion: reflect about mean
        mean = np.mean(state)
        state = 2 * mean - state
        state[target] = -(state[target] - 2*mean) + 2*mean  # undo double
        # Redo correctly:
        state_before_diff = state.copy()
        state_before_diff[target] *= -1  # undo target flip for diffusion
        # Actually do it right:

    # Redo from scratch properly
    state = np.ones(N, dtype=complex) / sqrt(N)
    for _ in range(k_opt):
        # Step 1: Oracle
        state[target] *= -1
        # Step 2: Diffusion D = 2|s> 0.1)[0])
# Should show peaks at multiples of N/r = 2

# ================================================================
# Qiskit: Deutsch-Jozsa
# ================================================================

try:
    from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
    from qiskit_aer import AerSimulator

    def dj_circuit_balanced(n=3):
        """Build a balanced Deutsch-Jozsa circuit."""
        qr = QuantumRegister(n + 1, 'q')
        cr = ClassicalRegister(n, 'c')
        qc = QuantumCircuit(qr, cr)
        qc.x(qr[n])                 # ancilla to |1>
        qc.h(range(n + 1))          # H to all qubits
        # Balanced oracle: XOR each input qubit into ancilla
        for i in range(n):
            qc.cx(qr[i], qr[n])
        qc.h(range(n))              # H to first n qubits
        qc.measure(range(n), range(n))
        return qc

    qc = dj_circuit_balanced(n=3)
    sim = AerSimulator()
    result = sim.run(qc, shots=1024).result()
    print("DJ balanced result:", result.get_counts())
    # Should show only non-zero strings (balanced)
except ImportError:
    print("Qiskit not installed; skipping circuit demo.")

/**
 * Manual state-vector simulation of Grover's search on 4 elements.
 * No external libraries. Runs in any JavaScript environment.
 */
(function() {
  const N = 4;         // database size = 2^2
  const TARGET = 2;    // we're searching for index 2

  // State vector: complex numbers as {re, im}
  let state = Array.from({length: N}, () => ({re: 1/Math.sqrt(N), im: 0}));

  function cMul(a, b) {
    return { re: a.re*b.re - a.im*b.im, im: a.re*b.im + a.im*b.re };
  }
  function cAdd(a, b) { return { re: a.re+b.re, im: a.im+b.im }; }
  function cScale(s, a) { return { re: s*a.re, im: s*a.im }; }
  function cNeg(a) { return { re: -a.re, im: -a.im }; }
  function cAbs2(a) { return a.re*a.re + a.im*a.im; }

  function oracle(state, target) {
    // Multiply target amplitude by -1
    const s = state.map((a, i) => i === target ? cNeg(a) : {...a});
    return s;
  }

  function diffusion(state) {
    // D = 2|s> = uniform superposition
    // New amplitude: 2 * mean - old amplitude
    const n = state.length;
    const meanRe = state.reduce((s, a) => s + a.re, 0) / n;
    const meanIm = state.reduce((s, a) => s + a.im, 0) / n;
    return state.map(a => ({
      re: 2*meanRe - a.re,
      im: 2*meanIm - a.im
    }));
  }

  const kOpt = Math.floor(Math.PI * Math.sqrt(N) / 4);
  console.log(`Searching ${N} items for target=${TARGET}`);
  console.log(`Optimal iterations: ${kOpt}`);

  // Initial probabilities
  console.log('Initial probabilities:', state.map(a => cAbs2(a).toFixed(4)));

  for (let k = 0; k < kOpt; k++) {
    state = oracle(state, TARGET);
    state = diffusion(state);
    const probs = state.map(a => cAbs2(a).toFixed(4));
    console.log(`After iter ${k+1}:`, probs);
  }

  const foundIdx = state.reduce((best, a, i) =>
    cAbs2(a) > cAbs2(state[best]) ? i : best, 0);

  console.log(`Measured: ${foundIdx} (target was ${TARGET})`);
  console.log(`Probability of finding target: ${(cAbs2(state[TARGET])*100).toFixed(1)}%`);
})();

/**
 * QFT on a 4-element state vector (n=2 qubits).
 * Manual implementation without library dependencies.
 */
(function() {
  const N = 4;
  const TWO_PI = 2 * Math.PI;

  function qft(state) {
    // state is array of {re, im} complex numbers
    const out = Array.from({length: N}, () => ({re: 0, im: 0}));
    for (let k = 0; k < N; k++) {
      for (let j = 0; j < N; j++) {
        const angle = TWO_PI * j * k / N;
        const re = Math.cos(angle);
        const im = Math.sin(angle);
        out[k].re += (state[j].re * re - state[j].im * im) / Math.sqrt(N);
        out[k].im += (state[j].re * im + state[j].im * re) / Math.sqrt(N);
      }
    }
    return out;
  }

  // Signal with period 2 (alternating 1, -1, 1, -1)
  const signal = [
    {re: 1/Math.sqrt(2), im: 0},
    {re: 0, im: 0},
    {re: 1/Math.sqrt(2), im: 0},
    {re: 0, im: 0}
  ];
  const spectrum = qft(signal);
  console.log('QFT spectrum |amplitudes|^2:',
    spectrum.map(a => (a.re*a.re + a.im*a.im).toFixed(4)));
  // Peak at k=0 and k=2, corresponding to period N/2 = 2.
})();

/*
 * QFT on a classical complex array (n=3, N=8).
 * Illustrates the circuit structure via explicit gate application.
 * Compile: gcc -O2 -lm quantum_sim.c -o quantum_sim
 */
#include <stdio.h>
#include <math.h>
#include <complex.h>

#define N_QUBITS 3
#define N_STATES 8

typedef double complex cplx;

/* Apply Hadamard to qubit 'qubit' of an n-qubit state vector */
void apply_H(cplx *state, int qubit, int n) {
    int stride = 1 << qubit;
    int size   = 1 << n;
    for (int i = 0; i < size; i++) {
        if ((i & stride) == 0) {
            cplx a = state[i];
            cplx b = state[i | stride];
            state[i]        = (a + b) / sqrt(2.0);
            state[i|stride] = (a - b) / sqrt(2.0);
        }
    }
}

/* Apply controlled-Rk (phase gate) with control=ctrl, target=tgt */
void apply_CR(cplx *state, int ctrl, int tgt, int k, int n) {
    cplx phase = cexp(2.0 * M_PI * I / (1 << k));
    int size   = 1 << n;
    for (int i = 0; i < size; i++) {
        /* Both ctrl and tgt must be |1> to apply phase */
        if ((i >> ctrl & 1) && (i >> tgt & 1)) {
            state[i] *= phase;
        }
    }
}

/* Swap qubits a and b */
void swap_qubits(cplx *state, int a, int b, int n) {
    int size = 1 << n;
    for (int i = 0; i < size; i++) {
        int bit_a = (i >> a) & 1;
        int bit_b = (i >> b) & 1;
        if (bit_a != bit_b) {
            int j = i ^ (1 << a) ^ (1 << b);
            if (i < j) {
                cplx tmp = state[i];
                state[i] = state[j];
                state[j] = tmp;
            }
        }
    }
}

/* n-qubit QFT */
void qft(cplx *state, int n) {
    for (int qubit = n - 1; qubit >= 0; qubit--) {
        apply_H(state, qubit, n);
        for (int ctrl = qubit - 1; ctrl >= 0; ctrl--) {
            int k = qubit - ctrl + 1;
            apply_CR(state, ctrl, qubit, k, n);
        }
    }
    /* Reverse bit order */
    for (int i = 0; i < n / 2; i++) {
        swap_qubits(state, i, n - 1 - i, n);
    }
}

int main(void) {
    /* Period-4 signal: 1/sqrt(2), 0, 1/sqrt(2), 0, 0, 0, 0, 0 */
    cplx state[N_STATES] = {0};
    state[0] = 1.0 / sqrt(2.0);
    state[4] = 1.0 / sqrt(2.0);

    printf("Input state:\n");
    for (int i = 0; i < N_STATES; i++)
        printf("  |%d>: %.4f + %.4fi\n", i, creal(state[i]), cimag(state[i]));

    qft(state, N_QUBITS);

    printf("\nAfter QFT (probabilities):\n");
    for (int i = 0; i < N_STATES; i++) {
        double prob = creal(state[i]) * creal(state[i]) +
                      cimag(state[i]) * cimag(state[i]);
        if (prob > 0.01)
            printf("  |%d>: %.4f\n", i, prob);
    }
    /* Expected peaks: |0>, |2>  (multiples of N/period = 8/4 = 2) */
    return 0;
}

/**
 * C++ template state-vector simulator with Grover's search.
 * Compile: g++ -O2 -std=c++17 grover.cpp -o grover
 */
#include <iostream>
#include <vector>
#include <complex>
#include <cmath>
#include <algorithm>
#include <numeric>
#include <cassert>

using Complex = std::complex<double>;
using State   = std::vector<Complex>;

constexpr double PI = M_PI;

// ----------------------------------------------------------------
// State-vector helpers
// ----------------------------------------------------------------

State uniform_superposition(int n) {
    int N = 1 << n;
    return State(N, 1.0 / std::sqrt(N));
}

double prob(const State& sv, int idx) {
    return std::norm(sv[idx]);  // |amplitude|^2
}

int measure(const State& sv) {
    // Deterministic argmax (no randomness for demo)
    return std::max_element(sv.begin(), sv.end(),
        [](const Complex& a, const Complex& b){ return std::norm(a) < std::norm(b); })
        - sv.begin();
}

// ----------------------------------------------------------------
// Grover oracle: negate the amplitude of 'target'
// ----------------------------------------------------------------

void grover_oracle(State& sv, int target) {
    sv[target] = -sv[target];
}

// ----------------------------------------------------------------
// Grover diffusion: reflect about the mean
// ----------------------------------------------------------------

void grover_diffusion(State& sv) {
    int N = sv.size();
    Complex mean = std::accumulate(sv.begin(), sv.end(), Complex{0, 0})
                   / static_cast<double>(N);
    for (auto& a : sv)
        a = 2.0 * mean - a;
}

// ----------------------------------------------------------------
// Full Grover search
// ----------------------------------------------------------------

int grover_search(int target, int n) {
    int N = 1 << n;
    State sv = uniform_superposition(n);
    int k_opt = static_cast<int>(std::floor(PI * std::sqrt(N) / 4.0));

    std::cout << "Searching N=" << N << " items for target=" << target
              << " (" << k_opt << " iterations)\n";

    for (int k = 0; k < k_opt; k++) {
        grover_oracle(sv, target);
        grover_diffusion(sv);
        std::cout << "  iter " << k+1
                  << "  P(target) = " << prob(sv, target) << "\n";
    }
    return measure(sv);
}

// ----------------------------------------------------------------
// QFT
// ----------------------------------------------------------------

void apply_H_gate(State& sv, int qubit, int n) {
    int stride = 1 << qubit;
    int size   = 1 << n;
    const double s = 1.0 / std::sqrt(2.0);
    for (int i = 0; i < size; i++) {
        if ((i & stride) == 0) {
            Complex a = sv[i], b = sv[i | stride];
            sv[i]         = s * (a + b);
            sv[i | stride] = s * (a - b);
        }
    }
}

void apply_CR_gate(State& sv, int ctrl, int tgt, int k, int n) {
    Complex phase = std::exp(Complex{0, 2.0 * PI / (1 << k)});
    int size = 1 << n;
    for (int i = 0; i < size; i++)
        if ((i >> ctrl & 1) && (i >> tgt & 1))
            sv[i] *= phase;
}

void qft(State& sv, int n) {
    for (int q = n - 1; q >= 0; q--) {
        apply_H_gate(sv, q, n);
        for (int c = q - 1; c >= 0; c--)
            apply_CR_gate(sv, c, q, q - c + 1, n);
    }
    // Bit reversal
    for (int i = 0; i < n / 2; i++) {
        int bit_a = 1 << i, bit_b = 1 << (n-1-i);
        for (int s = 0; s < (1 << n); s++)
            if ((s & bit_a) != (s & bit_b ? bit_a : 0))
                ;  // swap handled below
        // Proper bit-reversal permutation
        int size = 1 << n;
        for (int s = 0; s < size; s++) {
            int t = s ^ ((s ^ (s << (n-1-2*i))) & (bit_b | bit_a));
            if (t > s) std::swap(sv[s], sv[t]);
        }
    }
}

int main() {
    // Grover demo
    int found = grover_search(7, 4);  // target=7 in N=16
    std::cout << "Found: " << found << "\n\n";

    // QFT demo: period-4 signal in N=8
    int n = 3, N = 8, period = 4;
    State sv(N, 0);
    for (int k = 0; k < N / period; k++)
        sv[k * period] = 1.0 / std::sqrt(N / period);

    std::cout << "QFT of period-" << period << " signal (N=" << N << "):\n";
    qft(sv, n);
    for (int i = 0; i < N; i++) {
        double p = std::norm(sv[i]);
        if (p > 0.01)
            std::cout << "  |" << i << "> prob=" << p << "\n";
    }
    return 0;
}

/**
 * Java: Grover's search and Deutsch-Jozsa state-vector simulation.
 * Compile: javac QuantumSim.java && java QuantumSim
 */
public class QuantumSim {

    // Complex number arithmetic
    static double[] add(double[] a, double[] b) {
        return new double[]{a[0]+b[0], a[1]+b[1]};
    }
    static double[] scale(double s, double[] a) {
        return new double[]{s*a[0], s*a[1]};
    }
    static double[] negate(double[] a) {
        return new double[]{-a[0], -a[1]};
    }
    static double abs2(double[] a) { return a[0]*a[0] + a[1]*a[1]; }

    // ----------------------------------------------------------------
    // Grover's search
    // ----------------------------------------------------------------
    static double[][] groverSearch(int target, int n) {
        int N = 1 << n;
        double[][] sv = new double[N][2];
        double amp = 1.0 / Math.sqrt(N);
        for (int i = 0; i < N; i++) sv[i][0] = amp;

        int kOpt = (int)(Math.PI * Math.sqrt(N) / 4.0);
        System.out.printf("Grover: N=%d, target=%d, iterations=%d%n", N, target, kOpt);

        for (int k = 0; k < kOpt; k++) {
            // Oracle: negate target
            sv[target] = negate(sv[target]);

            // Diffusion: reflect about mean
            double meanRe = 0, meanIm = 0;
            for (double[] a : sv) { meanRe += a[0]; meanIm += a[1]; }
            meanRe /= N; meanIm /= N;
            for (int i = 0; i < N; i++) {
                sv[i][0] = 2*meanRe - sv[i][0];
                sv[i][1] = 2*meanIm - sv[i][1];
            }
            System.out.printf("  iter %d  P(target)=%.4f%n", k+1, abs2(sv[target]));
        }
        return sv;
    }

    // ----------------------------------------------------------------
    // Deutsch-Jozsa
    // ----------------------------------------------------------------
    static String deutschJozsa(int[] f) {
        // f is a truth table of size N = 2^n
        int N = f.length;
        int n = 31 - Integer.numberOfLeadingZeros(N);  // log2(N)

        // State: n+1 qubits, size = 2*N
        double[][] sv = new double[2 * N][2];
        sv[1][0] = 1.0;  // |0...0>|1>

        // Apply H to all n+1 qubits
        for (int q = 0; q <= n; q++) applyH(sv, q, n + 1);

        // Oracle (phase kickback)
        for (int x = 0; x < N; x++) {
            if (f[x] == 1) {
                // Swap |x,0> and |x,1> (XOR into ancilla)
                double[] tmp = sv[2*x].clone();
                sv[2*x]   = sv[2*x+1];
                sv[2*x+1] = tmp;
            }
        }

        // Apply H to first n qubits
        for (int q = 0; q < n; q++) applyH(sv, q + 1, n + 1);
        // Note: bit ordering — qubit 0 is ancilla (index bit 0)

        // Measure: probability of first n qubits being all-zero
        double probAllZero = abs2(sv[0]) + abs2(sv[1]);
        return probAllZero > 0.99 ? "constant" : "balanced";
    }

    static void applyH(double[][] sv, int qubit, int totalQubits) {
        int stride = 1 << qubit;
        int size   = 1 << totalQubits;
        double s   = 1.0 / Math.sqrt(2.0);
        for (int i = 0; i < size; i++) {
            if ((i & stride) == 0) {
                double[] a = sv[i].clone();
                double[] b = sv[i | stride].clone();
                sv[i]         = scale(s, add(a, b));
                sv[i | stride] = scale(s, add(a, negate(b)));
            }
        }
    }

    public static void main(String[] args) {
        // Grover demo
        int n = 4, target = 11;
        double[][] sv = groverSearch(target, n);
        int found = 0;
        for (int i = 1; i < sv.length; i++)
            if (abs2(sv[i]) > abs2(sv[found])) found = i;
        System.out.printf("Measured: %d (expected %d)%n%n", found, target);

        // Deutsch-Jozsa demo
        int[] fConst = new int[8];  // constant-0
        int[] fBal   = {0,1,0,1,0,1,0,1};  // balanced
        System.out.println("Constant-0 function: " + deutschJozsa(fConst));
        System.out.println("Balanced function:   " + deutschJozsa(fBal));
    }
}

// Quantum circuit simulation in Go: QFT and Grover's search.
// Run: go run quantum_sim.go
package main

import (
	"fmt"
	"math"
	"math/cmplx"
)

// ----------------------------------------------------------------
// State vector type
// ----------------------------------------------------------------

type State []complex128

func uniformSuperposition(n int) State {
	N := 1 << n
	sv := make(State, N)
	amp := complex(1.0/math.Sqrt(float64(N)), 0)
	for i := range sv {
		sv[i] = amp
	}
	return sv
}

func prob(sv State, idx int) float64 {
	return math.Pow(cmplx.Abs(sv[idx]), 2)
}

func argmax(sv State) int {
	best := 0
	for i := 1; i < len(sv); i++ {
		if prob(sv, i) > prob(sv, best) {
			best = i
		}
	}
	return best
}

// ----------------------------------------------------------------
// Gate applications
// ----------------------------------------------------------------

func applyH(sv State, qubit, nQubits int) {
	stride := 1 << qubit
	size := 1 << nQubits
	s := complex(1.0/math.Sqrt(2), 0)
	for i := 0; i < size; i++ {
		if i&stride == 0 {
			a, b := sv[i], sv[i|stride]
			sv[i] = s * (a + b)
			sv[i|stride] = s * (a - b)
		}
	}
}

func applyCR(sv State, ctrl, tgt, k, nQubits int) {
	phase := cmplx.Exp(complex(0, 2*math.Pi/float64(1<<k)))
	size := 1 << nQubits
	for i := 0; i < size; i++ {
		if (i>>ctrl)&1 == 1 && (i>>tgt)&1 == 1 {
			sv[i] *= phase
		}
	}
}

func swapQubits(sv State, a, b, nQubits int) {
	size := 1 << nQubits
	for i := 0; i < size; i++ {
		bitA := (i >> a) & 1
		bitB := (i >> b) & 1
		if bitA != bitB {
			j := i ^ (1 << a) ^ (1 << b)
			if i < j {
				sv[i], sv[j] = sv[j], sv[i]
			}
		}
	}
}

// ----------------------------------------------------------------
// QFT
// ----------------------------------------------------------------

func qft(sv State, n int) {
	for q := n - 1; q >= 0; q-- {
		applyH(sv, q, n)
		for c := q - 1; c >= 0; c-- {
			applyCR(sv, c, q, q-c+1, n)
		}
	}
	// Bit reversal
	for i := 0; i < n/2; i++ {
		swapQubits(sv, i, n-1-i, n)
	}
}

// ----------------------------------------------------------------
// Grover's search
// ----------------------------------------------------------------

func groverOracle(sv State, target int) {
	sv[target] = -sv[target]
}

func groverDiffusion(sv State) {
	N := len(sv)
	var mean complex128
	for _, a := range sv {
		mean += a
	}
	mean /= complex(float64(N), 0)
	for i := range sv {
		sv[i] = 2*mean - sv[i]
	}
}

func groverSearch(target, n int) int {
	N := 1 << n
	sv := uniformSuperposition(n)
	kOpt := int(math.Pi * math.Sqrt(float64(N)) / 4.0)

	fmt.Printf("Grover: N=%d, target=%d, kOpt=%d\n", N, target, kOpt)
	for k := 0; k < kOpt; k++ {
		groverOracle(sv, target)
		groverDiffusion(sv)
		fmt.Printf("  iter %d  P(target)=%.4f\n", k+1, prob(sv, target))
	}
	return argmax(sv)
}

func main() {
	// ---- Grover demo ----
	found := groverSearch(13, 4)
	fmt.Printf("Measured: %d\n\n", found)

	// ---- QFT demo: period-4 signal in N=8 ----
	n, period := 3, 4
	N := 1 << n
	sv := make(State, N)
	for k := 0; k < N/period; k++ {
		sv[k*period] = complex(1.0/math.Sqrt(float64(N/period)), 0)
	}
	fmt.Printf("Input (period-%d signal, N=%d):\n", period, N)
	for i, a := range sv {
		if cmplx.Abs(a) > 0.01 {
			fmt.Printf("  |%d>: %.4f+%.4fi\n", i, real(a), imag(a))
		}
	}
	qft(sv, n)
	fmt.Println("After QFT (significant probabilities):")
	for i, a := range sv {
		p := prob(State(sv), i)
		_ = a
		if p > 0.01 {
			fmt.Printf("  |%d>: %.4f\n", i, p)
		}
	}
	// Expected peaks at |0> and |2> (multiples of N/period = 2)
}