13. Example Programs

13.1. Hello World - Simple Bell State

#include <cudaq.h>

struct bell {
  int operator()(int num_iters) __qpu__ {
    cudaq::qreg q(2);
    int nCorrect = 0;
    for (int i = 0; i < num_iters; i++) {
      h(q[0]);
      x<cudaq::ctrl>(q[0], q[1]);
      auto results = mz(q);
      if (results[0] == results[1])
        nCorrect++;

      reset(q[0]);
      reset(q[1]);
    }
    return nCorrect;
  }
};

int main() { printf("N Correct = %d\n", bell{}(100)); }

13.2. GHZ State Preparation and Sampling

#include <cudaq.h>

struct ghz {
  void operator()(const int n_qubits) __qpu__ {
    cudaq::qreg q(n_qubits);
    h(q[0]);
    for (int i = 0; i < n_qubits - 1; ++i)
      // note use of ctrl modifier
      x<cudaq::ctrl>(q[i], q[i+1]);

    mz(q);
  }
};

int main() {
  // Sample the state produced by the ghz kernel
  auto counts = cudaq::sample(ghz{}, 10);
  for (auto [bits, count] : counts) {
    printf("Observed %s %lu times.\n", bits.c_str(), count);
  }
  return 0;
}

13.3. Quantum Phase Estimation

#include <cudaq.h>
#include <cmath>

__qpu__ void iqft (cudaq::qspan<> q) {
  int N = q.size();
  // Swap qubits
  for (int i = 0; i < N / 2; ++i) {
    swap(q[i], q[N - i - 1]);
  }

  for (int i = 0; i < N - 1; ++i) {
    h(q[i]);
    int j = i + 1;
    for (int y = i; y >= 0; --y) {
      const double theta = -M_PI / std::pow(2.0, j - y);
      r1<cudaq::ctrl>(theta, q[j], q[y]);
    }
  }

  h(q[N - 1]);
}

struct PhaseEstimation {
  double operator()(const int n_counting_qubits, const int n_state_qubits,
                   auto&& statePrep, auto&& unitary) __qpu__ {
    // Allocate a register of qubits
    cudaq::qreg q(n_counting_qubits + n_state_qubits);

    // Extract sub-registers, one for the counting qubits
    // another for the eigen state register
    auto counting_qubits = q.front(n_counting_qubits);
    auto state_register = q.back(n_state_qubits);

    // Prepare the eigenstate
    statePrep(state_register);

    // Put the counting register into uniform superposition
    h(counting_qubits);

    // Perform ctrl-U^j
    for (int i = 0; i < n_counting_qubits; ++i) {
      for (int j = 0; j < 1UL << i; ++j) {
        cudaq::control(unitary, counting_qubits[i], state_register);
      }
    }

    // Apply inverse quantum fourier transform
    iqft(counting_qubits);

    // Measure and compute the phase...
    auto bits = mz(counting_qubits);
    return cudaq::to_integer(bits) / std::pow(2, n_counting_qubits);
  }
};

int main() {
  auto statePrep = [](cudaq::qspan<> q) __qpu__ { x(q); };
  auto unitary = [](cudaq::qspan<> q) __qpu__ { t(q); };
  double phase = PhaseEstimation{}(3, 1, statePrep, unitary);
  printf("Phase = %lf\n", phase);
}

13.4. Deuteron Binding Energy Parameter Sweep

#include <cudaq.h>
#include <cudaq/algorithm.h>

struct deuteron_n2_ansatz {
  void operator()(double theta) __qpu__ {
    cudaq::qreg q(2);
    x(q[0]);
    ry(theta, q[1]);
    x<cudaq::ctrl>(q[1], q[0]);
  }
};

int main() {
  using namespace cudaq::spin;
  cudaq::spin_op h = 5.907 - 2.1433 * x(0) * x(1) - 2.1433 * y(0) * y(1) +
              .21829 * z(0) - 6.125 * z(1);

  // Perform parameter sweep for deuteron N=2 Hamiltonian
  const auto param_space = cudaq::linspace(-M_PI, M_PI, 25);
  for (const auto& param : param_space) {
    // KERNEL::observe(...) <==>
    // E(params...) = <psi(params...) | H | psi(params...)>
    double energy_at_param = cudaq::observe(deuteron_n2_ansatz{}, h, param);
    printf("<H>(%lf) = %lf\n", param, energy_at_param);
  }
}

13.5. Grover’s Algorithm

#include <cudaq.h>

__qpu__ void reflect_about_uniform(cudaq::qspan<> q) {
  auto ctrl_qubits = q.front(q.size() - 1);
  auto& last_qubit = q.back();
  // Compute (U) Action (V) produces
  // U V U::Adjoint
  cudaq::compute_action([&]() {
                         h(q);
                         x(q);
                       },
                       [&]() {
                         z<cudaq::ctrl>(ctrl_qubits, last_qubit);
                       }
                      );
}

struct run_grover {
  template <typename CallableKernel>
  auto operator()(const int n_qubits, const int n_iterations,
                        CallableKernel&& oracle) __qpu__ {
    cudaq::qreg q(n_qubits);
    h(q);
    for (int i = 0; i < n_iterations; i++) {
      oracle(q);
      reflect_about_uniform(q);
    }
    mz(q);
  }
};

struct oracle {
  void operator()(cudaq::qspan<> q) __qpu__ {
    cz(q[0], q[2]);
    cz(q[1], q[2]);
  }
};

int main() {
  auto counts = cudaq::sample(run_grover{}, 3, 1, oracle{});
  counts.dump();
  return 0;
}

13.6. Iterative Phase Estimation

 #include <cudaq.h>

struct iqpe {
  void operator()() __qpu__ {
    cudaq::qreg<2> q;
    h(q[0]);
    x(q[1]);
    for (int i = 0; i < 8; i++)
      r1<cudaq::ctrl>(-5 * M_PI / 8., q[0], q[1]);

    h(q[0]);
    auto cr0 = mz(q[0]);
    reset(q[0]);

    h(q[0]);
    for (int i = 0; i < 4; i++)
      r1<cudaq::ctrl>(-5 * M_PI / 8., q[0], q[1]);

    if (cr0)
      rz(-M_PI / 2., q[0]);

    h(q[0]);
    auto cr1 = mz(q[0]);
    reset(q[0]);

    h(q[0]);
    for (int i = 0; i < 2; i++)
      r1<cudaq::ctrl>(-5 * M_PI / 8., q[0], q[1]);

    if (cr0)
      rz(-M_PI / 4., q[0]);

    if (cr1)
      rz(-M_PI / 2., q[0]);

    h(q[0]);
    auto cr2 = mz(q[0]);
    reset(q[0]);
    h(q[0]);
    r1<cudaq::ctrl>(-5 * M_PI / 8., q[0], q[1]);

    if (cr0)
      rz(-M_PI / 8., q[0]);

    if (cr1)
      rz(-M_PI_4, q[0]);

    if (cr2)
      rz(-M_PI_2, q[0]);

    h(q[0]);
    mz(q[0]);
  }
};

int main() {
  auto counts = cudaq::sample(/*shots*/ 100, iqpe{});
  counts.dump();
  return 0;
}