Let's Build Your Own Quantum Circuit Simulator in JavaScript!

Quantum Circuit Simulator in JS.js

// -----------------------------------------------------------------------------
// Content by rootEnginear | MIT License
// Read the Blog: https://rootenginear.gitbook.io/blog/qcom/qsim
// -----------------------------------------------------------------------------
import * as math from "mathjs";

// -----------------------------------------------------------------------------
// Utils
// -----------------------------------------------------------------------------
const toString = (statevector) => statevector.toString();

// -----------------------------------------------------------------------------
// Define Basis States
// -----------------------------------------------------------------------------
const ket_0 = math.matrix([[1], [0]]);
const ket_1 = math.matrix([[0], [1]]);

console.log("|0>", ket_0.toString());
console.log("|1>", ket_1.toString());

// -----------------------------------------------------------------------------
// Make Gates from Matrices
// -----------------------------------------------------------------------------
const i_matrix = math.matrix([
  [1, 0],
  [0, 1],
]);

const x_matrix = math.matrix([
  [0, 1],
  [1, 0],
]);

const applyGate = (gate) => (statevector) => math.multiply(gate, statevector);

const I = applyGate(i_matrix);
const X = applyGate(x_matrix);

const execute = (starting_statevector, ...instructions) =>
  instructions.reduce(
    (previous_statevector, current_instruction) =>
      current_instruction(previous_statevector),
    starting_statevector
  );

console.log("|0> ─I─");
console.log(execute(ket_0, I, toString));
console.log("|0> ─X─I─");
console.log(execute(ket_0, X, I, toString));

// -----------------------------------------------------------------------------
// Generalize Gates
// -----------------------------------------------------------------------------
const U = (theta, phi, lambda) =>
  math.matrix([
    [math.evaluate(`cos(${theta}/2)`), math.evaluate(`-e^(i*${lambda})*sin(${theta}/2)`)],
    [
      math.evaluate(`e^(i*${phi})*sin(${theta}/2)`),
      math.evaluate(`e^(i*${phi + lambda})*cos(${theta}/2)`),
    ],
  ]);

const z_matrix = U(0, 0, math.pi);
const Z = applyGate(z_matrix);

const h_matrix = U(math.pi / 2, 0, math.pi);
const H = applyGate(h_matrix);

console.log("|0> ─H─");
console.log(execute(ket_0, H, toString));

console.log("|0> ─H─Z─H─");
console.log(execute(ket_0, H, Z, H, toString)); // Should get |1>

// -----------------------------------------------------------------------------
// Measurement
// -----------------------------------------------------------------------------
const complexMatrixNorm = (statevector) =>
  statevector
    .map((state) => math.abs(state) ** 2)
    .toArray()
    .flat()
    .reduce((acc, cur) => acc + cur, 0);

const normalizeState = (state) => {
  const sqrt_norm = math.sqrt(complexMatrixNorm(state));
  // Norm will never < 0, but can = 0
  if (sqrt_norm === 0) throw new Error("Norm is 0!");
  return math.divide(state, math.sqrt(complexMatrixNorm(state)));
};

const postselect0_matrix = math.matrix([
  [1, 0],
  [0, 0],
]);

const postselect1_matrix = math.matrix([
  [0, 0],
  [0, 1],
]);

const M = (qubit) => (state) => {
  // 1.) Extract relevant probability
  const qubit_ket0_prob = state
    .toArray()
    .flat()
    .filter((_, i) => (i & (1 << qubit)) === 0)
    .map((c) => math.abs(c) ** 2)
    .reduce((a, b) => a + b, 0);

  // 2.) Categorize probability
  const qubit_count = math.log2(state.toArray().length);

  let postselect_set = new Array(qubit_count).fill(i_matrix);
  postselect_set[qubit] =
    Math.random() <= qubit_ket0_prob ? postselect0_matrix : postselect1_matrix;
  const postselector = postselect_set.reduce((prev, cur) => math.kron(cur, prev));

  // 3.) Apply postselect gate & normalize
  const postselected_state = math.multiply(postselector, state);
  try {
    return normalizeState(postselected_state);
  } catch (error) {
    // This can only happen if the postselected state is zero matrix.
    // In this case, we just return the `postselected_state` (because it's zero).
    return postselected_state;
  }
};

console.log("|0> ─H─M─");
new Array(10).fill().forEach(() => {
  // M(0) -> Measure at qubit 0
  console.log(execute(ket_0, H, M(0), toString));
});

// -----------------------------------------------------------------------------
// Repeat Execution & Aggregate Results
// -----------------------------------------------------------------------------
const repeatExecution =
  (starting_statevector, ...instructions) =>
  (shots) =>
    new Array(shots).fill().map(() => execute(starting_statevector, ...instructions));

console.log("|0> ─H─M─ × 100");
console.log(repeatExecution(ket_0, H, M(0), toString)(100));

const aggregateResult = (result_list) => {
  const aggregated_result = result_list
    .map((result) => result.toString())
    .reduce((cnt, cur) => {
      return (cnt[cur] = cnt[cur] + 1 || 1), cnt;
    }, {});

  return aggregated_result;
};

const runAndAggregate =
  (starting_statevector, ...instructions) =>
  (shots) =>
    aggregateResult(repeatExecution(starting_statevector, ...instructions)(shots));

console.log("Σ(|0> ─H─M─ × 1000)");
console.log(runAndAggregate(ket_0, H, M(0))(1000));

// -----------------------------------------------------------------------------
// Expanding Qubits
// -----------------------------------------------------------------------------
const ket_00 = math.kron(ket_0, ket_0);

const HX = applyGate(math.kron(h_matrix, x_matrix));

console.log(`|0> ─X─\n|0> ─H─`);
console.log(execute(ket_00, HX, toString));

// -----------------------------------------------------------------------------
// 2-qubit Gates
// -----------------------------------------------------------------------------
const IH = applyGate(math.kron(i_matrix, h_matrix));

const cnot_matrix = math.matrix([
  [1, 0, 0, 0],
  [0, 0, 0, 1],
  [0, 0, 1, 0],
  [0, 1, 0, 0],
]);
const CNOT = applyGate(cnot_matrix);

console.log("|0> ─H─╭@─M─\n|0> ───╰X─M─");
console.log(runAndAggregate(ket_00, IH, CNOT, M(0), M(1))(1000));

console.log("|0> ─H─╭@─M─\n|0> ───╰X───");
console.log(runAndAggregate(ket_00, IH, CNOT, M(0))(1000));

const IX = applyGate(math.kron(i_matrix, x_matrix));

const swap_matrix = math.matrix([
  [1, 0, 0, 0],
  [0, 0, 1, 0],
  [0, 1, 0, 0],
  [0, 0, 0, 1],
]);
const SWAP = applyGate(swap_matrix);

console.log(`|0> ─X─╭X─\n|0> ───╰X─`);
console.log(execute(ket_00, IX, SWAP, toString));

// -----------------------------------------------------------------------------
// Conjugate Transpose
// -----------------------------------------------------------------------------
const s_matrix = U(0, 0, math.pi / 2);
const S = applyGate(s_matrix);

const dag = (matrix) => math.ctranspose(matrix);

const sdg_matrix = dag(s_matrix);
const SDG = applyGate(sdg_matrix);

console.log("|0> ─H─S─Sdg─H─M─");
console.log(execute(ket_0, H, S, SDG, H, toString));

const ket_000 = math.kron(ket_0, math.kron(ket_0, ket_0));
const HHH = applyGate(math.kron(h_matrix, math.kron(h_matrix, h_matrix)));

// -----------------------------------------------------------------------------
// 3-qubit (and beyond!) Circuit
// -----------------------------------------------------------------------------
console.log("|0> ─H─M─\n|0> ─H─M─\n|0> ─H─M─");
console.log(runAndAggregate(ket_000, HHH, M(0), M(1), M(2))(1000));

// -----------------------------------------------------------------------------
// Adding More Controls
// -----------------------------------------------------------------------------
const directSum = (a, b) => {
  const [a_col, a_row] = a.size();
  const [b_col, b_row] = b.size();

  const top_right_matrix = math.zeros(a_col, b_row);
  const bottom_left_matrix = math.zeros(b_col, a_row);

  return math.concat(
    math.concat(a, top_right_matrix),
    math.concat(bottom_left_matrix, b),
    0
  );
};

const cx21_matrix = directSum(i_matrix, x_matrix);

const ii_matrix = math.kron(i_matrix, i_matrix);
const ccx321_matrix = directSum(ii_matrix, cx21_matrix);
const CCX321 = applyGate(ccx321_matrix);

const composeGates = (...gates) =>
  gates.reduce((composted, current) => math.multiply(composted, current));

const i3sw12_matrix = math.kron(i_matrix, swap_matrix);
const sw23i1_matrix = math.kron(swap_matrix, i_matrix);

const sw13_matrix = composeGates(i3sw12_matrix, sw23i1_matrix, i3sw12_matrix);
const SW13 = applyGate(sw13_matrix);

const IHX = applyGate(math.kron(i_matrix, math.kron(h_matrix, x_matrix)));

console.log("|0> ─X─╭X─╭X─╭X─M─\n|0> ─H─│──├@─│──M─\n|0> ───╰X─╰@─╰X─M─\n");
console.log(runAndAggregate(ket_000, IHX, SW13, CCX321, SW13, M(0), M(1), M(2))(1000));