hgp canonical codeword fixed

Author: mkangquantum

doc/intro.ipynb

@@ -98,15 +98,15 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Z0 :  [3 5 6 7 8 9] , X0 :  [ 36  60  72  84  96 108]\n",
-      "Z1 :  [ 0  1  5  6  8 10] , X1 :  [ 37  61  73  85  97 109]\n",
-      "Z2 :  [ 2  4  5  6  8 11] , X2 :  [ 38  62  74  86  98 110]\n",
-      "Z3 :  [15 17 18 19 20 21] , X3 :  [  0  12  60  72  96 120]\n",
-      "Z4 :  [12 13 17 18 20 22] , X4 :  [  1  13  61  73  97 121]\n",
-      "Z5 :  [14 16 17 18 20 23] , X5 :  [  2  14  62  74  98 122]\n",
-      "Z6 :  [27 29 30 31 32 33] , X6 :  [ 24  48  60  72  96 132]\n",
-      "Z7 :  [24 25 29 30 32 34] , X7 :  [ 25  49  61  73  97 133]\n",
-      "Z8 :  [26 28 29 30 32 35] , X8 :  [ 26  50  62  74  98 134]\n"
+      "Z0 :  [111 113 114 115 116 117] , X0 :  [ 45  69  81  93 105 117]\n",
+      "Z1 :  [108 109 113 114 116 118] , X1 :  [ 46  70  82  94 106 118]\n",
+      "Z2 :  [110 112 113 114 116 119] , X2 :  [ 47  71  83  95 107 119]\n",
+      "Z3 :  [123 125 126 127 128 129] , X3 :  [  9  21  69  81 105 129]\n",
+      "Z4 :  [120 121 125 126 128 130] , X4 :  [ 10  22  70  82 106 130]\n",
+      "Z5 :  [122 124 125 126 128 131] , X5 :  [ 11  23  71  83 107 131]\n",
+      "Z6 :  [135 137 138 139 140 141] , X6 :  [ 33  57  69  81 105 141]\n",
+      "Z7 :  [132 133 137 138 140 142] , X7 :  [ 34  58  70  82 106 142]\n",
+      "Z8 :  [134 136 137 138 140 143] , X8 :  [ 35  59  71  83 107 143]\n"
      ]
     }
    ],
@@ -292,7 +292,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 12,
    "id": "68975ae8",
    "metadata": {},
    "outputs": [],
@@ -304,7 +304,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 13,
    "id": "d1eb9f95",
    "metadata": {},
    "outputs": [
@@ -341,7 +341,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 14,
    "id": "56a210c1",
    "metadata": {},
    "outputs": [
@@ -367,7 +367,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 15,
    "id": "e5bdce14",
    "metadata": {},
    "outputs": [],
@@ -392,15 +392,15 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 16,
    "id": "59700b1d",
    "metadata": {},
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "100%|██████████| 100/100 [00:11<00:00,  8.75it/s]"
+      "100%|██████████| 100/100 [00:10<00:00,  9.45it/s]"
      ]
     },
     {
@@ -455,7 +455,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "100%|██████████| 100/100 [00:08<00:00, 11.63it/s]\n"
+      "100%|██████████| 100/100 [00:09<00:00, 10.29it/s]\n"
      ]
     }
    ],

src/quits/gf2_utility.py

@@ -0,0 +1,76 @@
+import numpy as np
+
+def gf2_rref_pivots(H: np.ndarray):
+    """Return pivot columns of H in RREF over GF(2)."""
+    A = (np.asarray(H) & 1).astype(np.uint8, copy=True)
+    m, n = A.shape
+    pivots = []
+    r = 0
+    for c in range(n):
+        if r >= m:
+            break
+        rows = np.where(A[r:, c] == 1)[0]
+        if rows.size == 0:
+            continue
+        p = r + int(rows[0])
+        if p != r:
+            A[[r, p], :] = A[[p, r], :]
+        # eliminate column c in all other rows
+        ones = np.where(A[:, c] == 1)[0]
+        ones = ones[ones != r]
+        if ones.size:
+            A[ones, :] ^= A[r, :]
+        pivots.append(c)
+        r += 1
+    return np.array(pivots, dtype=int)
+
+def gf2_coset_reps_rowspace(H: np.ndarray) -> np.ndarray:
+    """
+    Canonical reps for F2^n / rowspace(H):
+    pick standard basis e_j for the non-pivot columns of RREF(H).
+    """
+    H = (np.asarray(H) & 1).astype(np.uint8, copy=False)
+    n = H.shape[1]
+    piv = set(gf2_rref_pivots(H).tolist())
+    nonpiv = [c for c in range(n) if c not in piv]  # size = n - rank(H) = k (if full row rank)
+    E = np.zeros((len(nonpiv), n), dtype=np.uint8)
+    for t, c in enumerate(nonpiv):
+        E[t, c] = 1
+    return E
+
+def gf2_rank(H: np.ndarray) -> int:
+    return len(gf2_rref_pivots(H))
+
+def in_rowspace(v: np.ndarray, H: np.ndarray) -> bool:
+    """Check if v is in rowspace(H) over GF(2) by solving H^T a = v."""
+    H = (np.asarray(H) & 1).astype(np.uint8, copy=False)
+    v = (np.asarray(v).reshape(-1) & 1).astype(np.uint8, copy=False)
+    m, n = H.shape
+    # Solve H^T a = v for a in F2^m
+    A = H.T.copy()
+    b = v.copy()
+
+    # Gaussian elimination
+    Aug = np.concatenate([A, b[:, None]], axis=1)
+    r = 0
+    pivcol = -np.ones(n, dtype=int)
+    for c in range(m):
+        if r >= n:
+            break
+        rows = np.where(Aug[r:, c] == 1)[0]
+        if rows.size == 0:
+            continue
+        p = r + int(rows[0])
+        if p != r:
+            Aug[[r, p]] = Aug[[p, r]]
+        pivcol[r] = c
+        ones = np.where(Aug[:, c] == 1)[0]
+        ones = ones[ones != r]
+        if ones.size:
+            Aug[ones, c:] ^= Aug[r, c:]
+        r += 1
+
+    # infeasible row: 0...0 | 1
+    if np.any(np.all(Aug[:, :m] == 0, axis=1) & (Aug[:, m] == 1)):
+        return False
+    return True

src/quits/qldpc_code.py

@@ -8,6 +8,7 @@
 from scipy.linalg import circulant
 from ldpc.mod2.mod2_numpy import nullspace
 from .ldpc_utility import compute_lz_and_lx
+from .gf2_utility import *
 
 # Parent class 
 class QldpcCode:
@@ -196,25 +197,24 @@ def __init__(self, h1, h2):
             self.lz, self.lx = compute_lz_and_lx(self.hx, self.hz)
 
     def get_logicals(self):
-        '''
-        :return: Logical operators of the code as a list of tuples (logical_z, logical_x)
-                 where logical_z and logical_x are numpy arrays of shape (num_logicals, num_data_qubits)
-                 The logicals are written in the "canonical form" as described in arXiv:2204.10812
-        '''
-        lz = np.zeros((self.k1*self.k2, self.hz.shape[1]), dtype=int)
-        lx = np.zeros((self.k1*self.k2, self.hx.shape[1]), dtype=int)
+        """
+        Canonical logicals for your HGP convention, assuming k1^T=k2^T=0.
+        Returns:
+          lz, lx: shape (k1*k2, num_data_qubits) as uint8.
+        """
+        E1 = gf2_coset_reps_rowspace(self.h1)  # (k1, n1) if H1 full row rank
+        E2 = gf2_coset_reps_rowspace(self.h2)  # (k2, n2)
+
+        lz = np.zeros((self.k1 * self.k2, self.hz.shape[1]), dtype=np.uint8)
+        lx = np.zeros((self.k1 * self.k2, self.hx.shape[1]), dtype=np.uint8)
 
         cnt = 0
         for i in range(self.k2):
-            ei = np.zeros(self.h2.shape[1], dtype=int)
-            ei[i] = 1
             for j in range(self.k1):
-                ej = np.zeros(self.h1.shape[1], dtype=int)
-                ej[j] = 1
-
-                lz[cnt,:self.n1*self.n2] = np.kron(ei, self.l1[j,:])
-                lx[cnt,:self.n1*self.n2] = np.kron(self.l2[i,:], ej)
-
+                # Z: (E2_i ⊗ L1_j | 0)
+                lz[cnt, :self.n1 * self.n2] = np.kron(E2[i, :], self.l1[j, :]) & 1
+                # X: (L2_i ⊗ E1_j | 0)
+                lx[cnt, :self.n1 * self.n2] = np.kron(self.l2[i, :], E1[j, :]) & 1
                 cnt += 1
 
         return lz, lx