Add global phase gate.

Author: sergeisakov

apps/qsim_base_custatevec.cu

@@ -133,7 +133,7 @@ int main(int argc, char* argv[]) {
     }
 
     Simulator CreateSimulator() const {
-      return Simulator(custatevec_handle);
+      return Simulator(cublas_handle, custatevec_handle);
     }
 
     cublasHandle_t cublas_handle;

lib/circuit_qsim_parser.h

@@ -167,13 +167,19 @@ class CircuitQsimParser final {
     IO::errorf("invalid gate in %s in line %u.\n", provider.c_str(), line);
   }
 
+  /**
+   * Checks formatting for a zero-qubit gate parsed from 'ss'.
+   * @param ss Input stream containing the gate specification.
+   */
+  static bool ValidateGate(std::stringstream& ss) {
+    return ss && ss.peek() == std::stringstream::traits_type::eof();
+  }
+
   /**
    * Checks formatting for a single-qubit gate parsed from 'ss'.
    * @param ss Input stream containing the gate specification.
    * @param num_qubits Number of qubits, as defined at the start of the file.
    * @param q0 Index of the affected qubit.
-   * @param provider Circuit source; only used for error reporting.
-   * @param line Line number of the parsed gate; only used for error reporting.
    */
   static bool ValidateGate(std::stringstream& ss,
                            unsigned num_qubits, unsigned q0) {
@@ -261,7 +267,11 @@ class CircuitQsimParser final {
     unsigned q0, q1;
     fp_type phi, theta;
 
-    if (gate_name == "id1") {
+    if (gate_name == "p") {
+      ss >> phi;
+      if (!ValidateGate(ss)) return false;
+      gates.push_back(GateGPh<fp_type>::Create(time, phi));
+    } else if (gate_name == "id1") {
       ss >> q0;
       if (!ValidateGate(ss, num_qubits, q0)) return false;
       gates.push_back(GateId1<fp_type>::Create(time, q0));

lib/fuser.h

@@ -122,6 +122,38 @@ class Fuser {
     return epochs;
   }
 
+  template <typename GateSeq0, typename Parent, typename GateFused>
+  static void FuseZeroQubitGates(const GateSeq0& gate_seq0,
+                                 Parent parent, std::size_t first,
+                                 std::vector<GateFused>& fused_gates) {
+    GateFused* fuse_to = nullptr;
+
+    for (std::size_t i = first; i < fused_gates.size(); ++i) {
+      auto& fgate = fused_gates[i];
+
+      if (fgate.kind != gate::kMeasurement && fgate.kind != gate::kDecomp
+          && fgate.parent->controlled_by.size() == 0
+          && !fgate.parent->unfusible) {
+        fuse_to = &fgate;
+        break;
+      }
+    }
+
+    if (fuse_to != nullptr) {
+      // Fuse zero-qubit gates with the first available fused gate.
+      for (const auto& g : gate_seq0) {
+        fuse_to->gates.push_back(parent(g));
+      }
+    } else {
+      auto g0 = parent(gate_seq0[0]);
+      fused_gates.push_back({g0->kind, g0->time, {}, g0, {g0}, {}});
+
+      for (std::size_t i = 1; i < gate_seq0.size(); ++i) {
+        fused_gates.back().gates.push_back(parent(gate_seq0[i]));
+      }
+    }
+  }
+
  private:
   static bool AddBoundary(unsigned time, unsigned max_time,
                           std::vector<unsigned>& boundaries) {
@@ -144,7 +176,9 @@ inline void CalculateFusedMatrix(FusedGate& gate) {
   MatrixIdentity(unsigned{1} << gate.qubits.size(), gate.matrix);
 
   for (auto pgate : gate.gates) {
-    if (gate.qubits.size() == pgate->qubits.size()) {
+    if (pgate->qubits.size() == 0) {
+      MatrixScalarMultiply(pgate->matrix[0], pgate->matrix[1], gate.matrix);
+    } else if (gate.qubits.size() == pgate->qubits.size()) {
       MatrixMultiply(gate.qubits.size(), pgate->matrix, gate.matrix);
     } else {
       unsigned mask = 0;

lib/fuser_basic.h

@@ -161,16 +161,24 @@ class BasicGateFuser final : public Fuser<IO, Gate> {
     // Sequence of top level gates the other gates get fused to.
     std::vector<const RGate*> gates_seq;
 
+    // Sequence of zero-qubit gates.
+    std::vector<const RGate*> gates_seq0;
+
     // Lattice of gates: qubits "hyperplane" and time direction.
     std::vector<std::vector<const RGate*>> gates_lat(max_qubit1);
 
     // Current unfused gate.
     auto gate_it = gfirst;
 
+    std::size_t last_fused_gate_index = 0;
+
     for (std::size_t l = 0; l < times.size(); ++l) {
       gates_seq.resize(0);
       gates_seq.reserve(num_gates);
 
+      gates_seq0.resize(0);
+      gates_seq0.reserve(num_gates);
+
       for (unsigned k = 0; k < max_qubit1; ++k) {
         gates_lat[k].resize(0);
         gates_lat[k].reserve(128);
@@ -212,6 +220,8 @@ class BasicGateFuser final : public Fuser<IO, Gate> {
           gates_lat[gate.qubits[0]].push_back(&gate);
           gates_lat[gate.qubits[1]].push_back(&gate);
           gates_seq.push_back(&gate);
+        } else {
+          gates_seq0.push_back(&gate);
         }
       }
 
@@ -306,7 +316,14 @@ class BasicGateFuser final : public Fuser<IO, Gate> {
         gates_fused.push_back(std::move(gate_f));
       }
 
+      if (gates_seq0.size() != 0) {
+        Base::FuseZeroQubitGates(gates_seq0, [](const RGate* g) { return g; },
+                                 last_fused_gate_index, gates_fused);
+      }
+
       if (gate_it == glast) break;
+
+      last_fused_gate_index = gates_fused.size();
     }
 
     if (fuse_matrix) {

lib/fuser_mqubit.h

@@ -282,6 +282,7 @@ class MultiQubitGateFuser final : public Fuser<IO, Gate> {
     unsigned max_fused_size = std::min(unsigned{6}, param.max_fused_size);
     max_fused_size = std::min(max_fused_size, max_qubit1);
 
+    std::size_t last_fused_gate_index = 0;
     auto gate_it = gfirst;
 
     // Iterate over epochs.
@@ -432,6 +433,14 @@ class MultiQubitGateFuser final : public Fuser<IO, Gate> {
           FuseOrphanedGates(max_fused_size, stat, orphaned_gates, fused_gates);
         }
       }
+
+      if (fgates[0].size() != 0) {
+        Base::FuseZeroQubitGates(fgates[0],
+                                 [](const GateF* g) { return g->parent; },
+                                 last_fused_gate_index, fused_gates);
+      }
+
+      last_fused_gate_index = fused_gates.size();
     }
 
     if (fuse_matrix) {

lib/gates_cirq.h

@@ -74,6 +74,7 @@ enum GateKind {
   kMatrixGate1,  // One-qubit matrix gate.
   kMatrixGate2,  // Two-qubit matrix gate.
   kMatrixGate,   // Multi-qubit matrix gate.
+  kGlobalPhaseGate,
   kDecomp = gate::kDecomp,
   kMeasurement = gate::kMeasurement,
 };
@@ -85,6 +86,31 @@ constexpr double h_double = 0.5;
 constexpr double pi_double = 3.14159265358979323846264338327950288;
 constexpr double is2_double = 0.7071067811865475;
 
+// Gates from cirq/ops/global_phase_op.py:
+
+/**
+ * The global phase gate.
+ */
+template <typename fp_type>
+struct GlobalPhaseGate {
+  static constexpr GateKind kind = kGlobalPhaseGate;
+  static constexpr char name[] = "GlobalPhaseGate";
+  static constexpr unsigned num_qubits = 1;
+  static constexpr bool symmetric = true;
+
+  static GateCirq<fp_type> Create(unsigned time, fp_type phi) {
+    return Create(time, std::cos(phi), std::sin(phi));
+  }
+
+  static GateCirq<fp_type> Create(unsigned time, fp_type cp, fp_type sp) {
+    return CreateGate<GateCirq<fp_type>, GlobalPhaseGate>(
+        time, {}, {cp, sp}, {cp, sp});
+  }
+};
+
+template <typename fp_type>
+using global_phase_operation = GlobalPhaseGate<fp_type>;
+
 // Gates from cirq/ops/identity.py:
 
 /**

lib/gates_qsim.h

@@ -46,8 +46,9 @@ enum GateKind {
   kGateIS,      // iSwap
   kGateFS,      // fSim
   kGateCP,      // control phase
-  kGateMatrix1, // One-qubit matrix gate.
-  kGateMatrix2, // Two-qubit matrix gate.
+  kGateMatrix1, // one-qubit matrix gate
+  kGateMatrix2, // two-qubit matrix gate
+  kGateGPh,     // global phase gate
   kDecomp = gate::kDecomp,
   kMeasurement = gate::kMeasurement,
 };
@@ -59,6 +60,28 @@ using GateQSim = Gate<fp_type, GateKind>;
 constexpr double h_double = 0.5;
 constexpr double is2_double = 0.7071067811865475;
 
+// Zero-qubit gates:
+
+/**
+ * The global phase gate.
+ */
+template <typename fp_type>
+struct GateGPh {
+  static constexpr GateKind kind = kGateGPh;
+  static constexpr char name[] = "p";
+  static constexpr unsigned num_qubits = 1;
+  static constexpr bool symmetric = true;
+
+  static GateQSim<fp_type> Create(unsigned time, fp_type phi) {
+    return Create(time, std::cos(phi), std::sin(phi));
+  }
+
+  static GateQSim<fp_type> Create(unsigned time, fp_type cp, fp_type sp) {
+    return CreateGate<GateQSim<fp_type>, GateGPh>(
+        time, {}, {cp, sp}, {cp, sp});
+  }
+};
+
 // One-qubit gates:
 
 /**

lib/matrix.h

@@ -172,7 +172,7 @@ inline void MatrixMultiply(unsigned mask1,
 }
 
 /**
- * Multiply a matrix by a scalar value.
+ * Multiply a matrix by a real scalar value.
  * @c Scalar value.
  * @m Input matrix to be multiplied. Output matrix.
  */
@@ -183,6 +183,23 @@ inline void MatrixScalarMultiply(fp_type1 c, Matrix<fp_type2>& m) {
   }
 }
 
+/**
+ * Multiply a matrix by a complex scalar value.
+ * @re Real part of scalar value.
+ * @im Imaginary part of scalar value.
+ * @m Input matrix to be multiplied. Output matrix.
+ */
+template <typename fp_type1, typename fp_type2>
+inline void MatrixScalarMultiply(
+    fp_type1 re, fp_type1 im, Matrix<fp_type2>& m) {
+  for (unsigned i = 0; i < m.size() / 2; ++i) {
+    fp_type2 re0 = m[2 * i + 0];
+    fp_type2 im0 = m[2 * i + 1];
+    m[2 * i + 0] = re * re0 - im * im0;
+    m[2 * i + 1] = re * im0 + im * re0;
+  }
+}
+
 /**
  * Daggers a matrix.
  * @n Number of matrix rows (columns).

lib/simulator.h

@@ -41,22 +41,27 @@ class SimulatorBase {
                           uint64_t* ms, uint64_t* xss) {
     constexpr unsigned hsize = 1 << H;
 
-    uint64_t xs[H];
-
-    xs[0] = uint64_t{1} << (qs[L] + 1);
-    ms[0] = (uint64_t{1} << qs[L]) - 1;
-    for (unsigned i = 1; i < H; ++i) {
-      xs[i] = uint64_t{1} << (qs[L + i] + 1);
-      ms[i] = ((uint64_t{1} << qs[L + i]) - 1) ^ (xs[i - 1] - 1);
-    }
-    ms[H] = ((uint64_t{1} << num_qubits) - 1) ^ (xs[H - 1] - 1);
+    if (H == 0) {
+      ms[0] = uint64_t(-1);
+      xss[0] = 0;
+    } else {
+      uint64_t xs[H + 1];
+
+      xs[0] = uint64_t{1} << (qs[L] + 1);
+      ms[0] = (uint64_t{1} << qs[L]) - 1;
+      for (unsigned i = 1; i < H; ++i) {
+        xs[i] = uint64_t{1} << (qs[L + i] + 1);
+        ms[i] = ((uint64_t{1} << qs[L + i]) - 1) ^ (xs[i - 1] - 1);
+      }
+      ms[H] = ((uint64_t{1} << num_qubits) - 1) ^ (xs[H - 1] - 1);
 
-    for (unsigned i = 0; i < hsize; ++i) {
-      uint64_t a = 0;
-      for (uint64_t k = 0; k < H; ++k) {
-        a += xs[k] * ((i >> k) & 1);
+      for (unsigned i = 0; i < hsize; ++i) {
+        uint64_t a = 0;
+        for (uint64_t k = 0; k < H; ++k) {
+          a += xs[k] * ((i >> k) & 1);
+        }
+        xss[i] = a;
       }
-      xss[i] = a;
     }
   }
 

lib/simulator_avx.h

@@ -51,6 +51,9 @@ class SimulatorAVX final : public SimulatorBase {
     // Assume qs[0] < qs[1] < qs[2] < ... .
 
     switch (qs.size()) {
+    case 0:
+      ApplyGateH<0>(qs, matrix, state);
+      break;
     case 1:
       if (qs[0] > 2) {
         ApplyGateH<1>(qs, matrix, state);
@@ -137,6 +140,13 @@ class SimulatorAVX final : public SimulatorBase {
     }
 
     switch (qs.size()) {
+    case 0:
+      if (cqs[0] > 2) {
+        ApplyControlledGateHH<0>(qs, cqs, cvals, matrix, state);
+      } else {
+        ApplyControlledGateHL<0>(qs, cqs, cvals, matrix, state);
+      }
+      break;
     case 1:
       if (qs[0] > 2) {
         if (cqs[0] > 2) {
@@ -322,7 +332,6 @@ class SimulatorAVX final : public SimulatorBase {
   }
 
  private:
-
 #ifdef __BMI2__
 
   template <unsigned H>

lib/simulator_avx512.h

@@ -51,6 +51,9 @@ class SimulatorAVX512 final : public SimulatorBase {
     // Assume qs[0] < qs[1] < qs[2] < ... .
 
     switch (qs.size()) {
+    case 0:
+      ApplyGateH<0>(qs, matrix, state);
+      break;
     case 1:
       if (qs[0] > 3) {
         ApplyGateH<1>(qs, matrix, state);
@@ -143,6 +146,13 @@ class SimulatorAVX512 final : public SimulatorBase {
     }
 
     switch (qs.size()) {
+    case 0:
+      if (cqs[0] > 3) {
+        ApplyControlledGateHH<0>(qs, cqs, cvals, matrix, state);
+      } else {
+        ApplyControlledGateHL<0>(qs, cqs, cvals, matrix, state);
+      }
+      break;
     case 1:
       if (qs[0] > 3) {
         if (cqs[0] > 3) {

lib/simulator_basic.h

@@ -49,6 +49,9 @@ class SimulatorBasic final : public SimulatorBase {
     // Assume qs[0] < qs[1] < qs[2] < ... .
 
     switch (qs.size()) {
+    case 0:
+      ApplyGateH<0>(qs, matrix, state);
+      break;
     case 1:
       ApplyGateH<1>(qs, matrix, state);
       break;
@@ -92,6 +95,9 @@ class SimulatorBasic final : public SimulatorBase {
     }
 
     switch (qs.size()) {
+    case 0:
+      ApplyControlledGateH<0>(qs, cqs, cvals, matrix, state);
+      break;
     case 1:
       ApplyControlledGateH<1>(qs, cqs, cvals, matrix, state);
       break;