symplectic-chp 0.1.0.0 → 0.1.0.1

raw patch · 6 files changed

+827/−3 lines, 6 filesdep +timedep ~stim-parsernew-component:exe:verify-large-tableauPVP: minor bump suggested

API additions: PVP suggests at least a minor version bump

Dependencies added: time

Dependency ranges changed: stim-parser

API changes (from Hackage documentation)

Files

• CHANGELOG.md +7/−0
• README.md +111/−0
• app/VerifyLargeTableau.hs +370/−0
• src/SymplecticCHP/BitVec.hs +104/−0
• src/SymplecticCHP/LargeTableau.hs +217/−0
• symplectic-chp.cabal +18/−3

@@ -1,5 +1,12 @@ # Revision history for symplectic-chp +## 0.1.0.1 -- 2026-07-17++* Support `stim-parser` 0.2.0.0: relax upper bound to `< 0.3`.  +  The `Float` → `Double` change in `stim-parser`'s AST is transparent to+  `symplectic-chp` because probabilities and coordinates are ignored during+  translation.+ ## 0.1.0.0 -- 2026-04-06  ### Major Features

@@ -92,6 +92,53 @@ - O(1) indexing with bounds guarantees - No out-of-bounds errors at runtime +## Arbitrary Qubit Support++In addition to the standard `Tableau n` (optimized for up to 64 qubits using `Word64`), we now provide `LargeTableau` for **arbitrary qubit counts** using chunked bit-vector storage.++### Standard API (≤ 64 qubits, Word64-optimized)++```haskell+import SymplecticCHP++-- Type-safe, compile-time sized (fast for small circuits)+bellState :: Tableau 2+bellState = +  applyGate (CNOT 0 1) $+  applyGate (Local (Hadamard 0)) $+  emptyTableau @2+```++### Large Tableau API (arbitrary qubits, BitVec-backed)++```haskell+import SymplecticCHP.BitVec+import SymplecticCHP.LargeTableau++-- Runtime-sized, works with any number of qubits+largeCircuit :: Int -> LargeTableau+largeCircuit n = +  largeApplyGate (LargeCNOT 0 1) $+  largeApplyGate (LargeLocal (LargeHadamard 0)) $+  largeEmpty n++-- Simulate 1000-qubit circuit+main = do+  let tab0 = largeEmpty 1000+  let tab1 = largeApplyGate (LargeLocal (LargeHadamard 0)) tab0+  print $ largeIsValid tab1  -- True+```++### Feature Comparison++| Feature | Standard `Tableau n` | `LargeTableau` |+|---------|---------------------|----------------|+| Max qubits | 64 (Word64) | Unlimited (BitVec) |+| Storage | Unboxed Word64 | Chunked Vector Word64 |+| Type safety | Compile-time n | Runtime n |+| Performance | Optimal | Good (slight overhead) |+| Best for | Small circuits, education | Large QEC codes |+ ## Quick Start  ### Library Usage@@ -214,6 +261,70 @@ - `.stim` - The circuit file - `.expected` - Expected results for automated testing - `.derive.md` - Mathematical derivation of the circuit's behavior++### Verifying LargeTableau Correctness++We provide a comprehensive verification suite to validate the `LargeTableau` implementation against known quantum states:++```bash+# Build the verification executable+cabal build verify-large-tableau++# Run with default settings (10,000 qubits for Bell pairs test)+cabal run verify-large-tableau++# Run with custom qubit counts+cabal run verify-large-tableau -- --bell-pairs=5000 --rep-code=1000 --random=50++# Quick test with smaller circuits+cabal run verify-large-tableau -- --bell-pairs=100 --rep-code=100 --random=10+```++#### Verification Tests++| Test | Description | Qubits | Verification Method |+|------|-------------|--------|---------------------|+| **Bell Pairs** | Creates N/2 independent \|Φ⁺⟩ states | Configurable (default 10,000) | Stabilizer validity, pair-wise commutation |+| **Repetition Code** | Creates \|+⋯+⟩ GHZ-like state | Configurable (default 1,000) | X₀Xᵢ stabilizer properties |+| **Phase Identity** | Verifies S² = Z algebra | 100 | Eigenvalue verification |+| **Random Circuits** | Property-based fuzzing | 100 | Tableau validity preservation |+| **Performance** | Benchmarks gate throughput | 100-10,000 | Timing measurements |++#### Example Output++```+========================================+  LargeTableau Verification Suite+========================================+Configuration:+  Bell pairs test: 10000 qubits+  Rep code test: 1000 qubits+  Random circuits: 100++=== Test 1: Pairwise Bell States ===+Creating 5000 Bell pairs with 10000 qubits...+Circuit creation time: 0.23s+Tableau valid: True+Sampled stabilizers commute: True+Stabilizer count: 10000 (expected: 10000)++=== Test 5: Performance Benchmark ===+Benchmarking 10000 qubits:+  Creation: 0.01s+  100 Hadamards: 0.15s+  100 CNOTs: 0.32s+  Tableau valid: True+  Estimated memory: 2500 KB++========================================+  Summary+========================================+Total time: 2.34s++✅ ALL TESTS PASSED+```++These tests ensure that `LargeTableau` produces correct results for arbitrary qubit counts by validating against analytically known quantum states.  ## Learn More 

@@ -0,0 +1,370 @@+{-# LANGUAGE TypeApplications #-}++-- | Verification tests for LargeTableau implementation.+-- Validates correctness against known quantum states.+module Main where++import Control.Monad (foldM, when, forM_, forM)+import Data.Time (diffUTCTime, getCurrentTime)+import System.Environment (getArgs)+import System.Exit (exitFailure, exitSuccess)+import System.Random (randomRIO)+import Text.Printf (printf)++import Data.Bits (testBit, setBit, clearBit, xor)+import Data.Word (Word64)+import SymplecticCHP.BitVec+import SymplecticCHP.LargeTableau+import qualified Data.Vector as V++-- ============================================================================+-- Test Configuration+-- ============================================================================++data Config = Config+  { bellPairsQubits :: Int    -- ^ Qubits for Bell pairs test (default 10000)+  , repCodeQubits :: Int      -- ^ Qubits for repetition code test (default 1000)+  , randomCircuits :: Int     -- ^ Number of random circuits (default 100)+  , verbose :: Bool           -- ^ Verbose output+  }++defaultConfig :: Config+defaultConfig = Config+  { bellPairsQubits = 10000+  , repCodeQubits = 1000+  , randomCircuits = 100+  , verbose = False+  }++-- ============================================================================+-- Test 1: Pairwise Bell States+-- ============================================================================++-- | Create N/2 independent Bell pairs: |Φ⁺⟩^⊗N/2+-- Circuit: For each pair (2k, 2k+1): H on 2k, then CNOT(2k, 2k+1)+createBellPairsLarge :: Int -> LargeTableau+createBellPairsLarge n+  | odd n = error "createBellPairsLarge: qubit count must be even"+  | otherwise = +      let tab0 = largeEmpty n+          -- Apply H to even qubits+          tab1 = foldl (\t i -> largeApplyGate (LargeLocal (LargeHadamard i)) t) +                       tab0 [0,2..n-2]+          -- Apply CNOT from even to odd+          tab2 = foldl (\t i -> largeApplyGate (LargeCNOT i (i+1)) t) +                       tab1 [0,2..n-2]+      in tab2++-- | Test 1: Pairwise Bell States+testBellPairs :: Config -> IO Bool+testBellPairs config = do+  let n = bellPairsQubits config+  putStrLn $ "\n=== Test 1: Pairwise Bell States ==="+  putStrLn $ "Creating " ++ show (n `div` 2) ++ " Bell pairs with " ++ show n ++ " qubits..."+  +  start <- getCurrentTime+  let tableau = createBellPairsLarge n+  mid <- getCurrentTime+  +  putStrLn $ "Circuit creation time: " ++ show (diffUTCTime mid start)+  +  -- Verify tableau validity+  let valid = largeIsValid tableau+  putStrLn $ "Tableau valid: " ++ show valid+  +  -- Verify each pair commutes with all others (independent pairs should commute)+  putStrLn "Checking pair-wise commutation..."+  let stabs = ltStabs tableau+      nStabs = V.length stabs+      -- Sample some pairs to check (checking all would be O(n²))+      sampleIndices = take 100 [0..nStabs-1]  -- Check first 100+      checkCommutation = all (\i -> +        all (\j -> +          if i == j then True+          else not (lpOmega (stabs V.! i) (stabs V.! j)))+          sampleIndices)+        sampleIndices+  +  putStrLn $ "Sampled stabilizers commute: " ++ show checkCommutation+  putStrLn $ "  (checked " ++ show (length sampleIndices) ++ " stabilizers)"+  +  -- Verify stabilizer count+  putStrLn $ "Stabilizer count: " ++ show nStabs ++ " (expected: " ++ show n ++ ")"+  let countCorrect = nStabs == n+  +  end <- getCurrentTime+  putStrLn $ "Total test time: " ++ show (diffUTCTime end start)+  +  return (valid && checkCommutation && countCorrect)++-- ============================================================================+-- Test 2: Repetition Code State+-- ============================================================================++-- | Create repetition code state: |+⋯+⟩ = (|0⋯0⟩ + |1⋯1⟩)/√2+-- Circuit: H on qubit 0, then CNOT(0, i) for all i > 0+createRepCodeLarge :: Int -> LargeTableau+createRepCodeLarge n+  | n < 2 = error "createRepCodeLarge: need at least 2 qubits"+  | otherwise =+      let tab0 = largeEmpty n+          tab1 = largeApplyGate (LargeLocal (LargeHadamard 0)) tab0+          tab2 = foldl (\t i -> largeApplyGate (LargeCNOT 0 i) t) tab1 [1..n-1]+      in tab2++-- | Test 2: Repetition Code State+testRepCode :: Config -> IO Bool+testRepCode config = do+  let n = repCodeQubits config+  putStrLn $ "\n=== Test 2: Repetition Code State ==="+  putStrLn $ "Creating repetition code with " ++ show n ++ " qubits..."+  +  start <- getCurrentTime+  let tableau = createRepCodeLarge n+  mid <- getCurrentTime+  +  putStrLn $ "Circuit creation time: " ++ show (diffUTCTime mid start)+  +  -- Verify tableau validity+  let valid = largeIsValid tableau+  putStrLn $ "Tableau valid: " ++ show valid+  +  -- For repetition code, check that X_0 X_i are stabilizers for sampled i+  putStrLn "Checking X_0 X_i stabilizers (sampled)..."+  let sampleIndices = take 50 [1..n-1]  -- Check first 50+      checkXX = all (\i ->+        let xVec = bvSetBit (bvSetBit (bvEmpty n) 0) i+            pauli = LargePauli xVec (bvEmpty n) 0 n+        in largeIsDeterminate tableau pauli) sampleIndices+  +  putStrLn $ "Sampled X_0 X_i are stabilizers: " ++ show checkXX+  putStrLn $ "  (checked " ++ show (length sampleIndices) ++ " indices)"+  +  end <- getCurrentTime+  putStrLn $ "Total test time: " ++ show (diffUTCTime end start)+  +  return (valid && checkXX)++-- ============================================================================+-- Test 3: Phase Gate Identity (S² = Z)+-- ============================================================================++-- | Test that S² = Z on a superposition state+-- Start with |+⟩, apply S twice, should get |−⟩ (Z eigenvalue -1)+testPhaseIdentity :: Config -> IO Bool+testPhaseIdentity _config = do+  putStrLn $ "\n=== Test 3: Phase Gate Identity (S² = Z) ==="+  +  let n = 100  -- Use 100 qubits+  putStrLn $ "Testing on " ++ show n ++ " qubits..."+  +  start <- getCurrentTime+  +  -- Create |+⟩ states on all qubits+  let tab0 = largeEmpty n+      tab1 = foldl (\t i -> largeApplyGate (LargeLocal (LargeHadamard i)) t) +                   tab0 [0..n-1]+  +  -- Apply S to all qubits (first time)+  let tab2 = foldl (\t i -> largeApplyGate (LargeLocal (LargePhase i)) t) +                   tab1 [0..n-1]+  +  -- Apply S to all qubits (second time)+  let tab3 = foldl (\t i -> largeApplyGate (LargeLocal (LargePhase i)) t) +                   tab2 [0..n-1]+  +  mid <- getCurrentTime+  putStrLn $ "Circuit time: " ++ show (diffUTCTime mid start)+  +  -- Verify tableau valid+  let valid = largeIsValid tab3+  putStrLn $ "Tableau valid: " ++ show valid+  +  -- After H then S twice, we should have -X stabilizers (|−⟩ state)+  -- Check first few qubits+  let sampleIndices = take 10 [0..n-1]+      checkNegX = all (\i ->+        let xVec = bvSetBit (bvEmpty n) i+            pauli = LargePauli xVec (bvEmpty n) 2 n  -- phase 2 = -1+        in largeIsDeterminate tab3 pauli) sampleIndices+  +  putStrLn $ "Sampled qubits have -X stabilizer: " ++ show checkNegX+  +  end <- getCurrentTime+  putStrLn $ "Total test time: " ++ show (diffUTCTime end start)+  +  return (valid && checkNegX)++-- ============================================================================+-- Test 4: Random Circuit Property Tests+-- ============================================================================++-- | Apply random Clifford gate to large tableau+applyRandomGate :: Int -> LargeTableau -> IO LargeTableau+applyRandomGate n tab = do+  gateType <- randomRIO (0, 2) :: IO Int+  case gateType of+    0 -> do -- Hadamard+      q <- randomRIO (0, n-1)+      return $ largeApplyGate (LargeLocal (LargeHadamard q)) tab+    1 -> do -- Phase+      q <- randomRIO (0, n-1)+      return $ largeApplyGate (LargeLocal (LargePhase q)) tab+    2 -> do -- CNOT+      c <- randomRIO (0, n-1)+      t <- randomRIO (0, n-1)+      if c == t +        then return tab+        else return $ largeApplyGate (LargeCNOT c t) tab+    _ -> return tab++-- | Test 4: Random circuits preserve validity+testRandomCircuits :: Config -> IO Bool+testRandomCircuits config = do+  let n = 100  -- Use 100 qubits for random tests (fast but large enough)+      numCircuits = randomCircuits config+  putStrLn $ "\n=== Test 4: Random Circuit Properties ==="+  putStrLn $ "Testing " ++ show numCircuits ++ " random circuits on " ++ show n ++ " qubits..."+  +  results <- forM [1..numCircuits] $ \i -> do+    when (i `mod` 10 == 0) $ putStrLn $ "  Circuit " ++ show i ++ "/" ++ show numCircuits+    +    -- Generate random circuit+    let initialTab = largeEmpty n+    numGates <- randomRIO (10, 100) :: IO Int+    finalTab <- foldM (\t _ -> applyRandomGate n t) initialTab [1..numGates]+    +    -- Check invariants+    let valid = largeIsValid finalTab+        stabs = ltStabs finalTab+        destabs = ltDestabs finalTab+        correctSize = V.length stabs == n && V.length destabs == n+    +    return (valid && correctSize)+  +  let allPass = and results+  putStrLn $ "All random circuits valid: " ++ show allPass+  putStrLn $ "Passed: " ++ show (length (filter id results)) ++ "/" ++ show numCircuits+  +  return allPass++-- ============================================================================+-- Test 5: Performance Benchmark+-- ============================================================================++-- | Benchmark gate application performance+benchmarkGates :: IO ()+benchmarkGates = do+  putStrLn $ "\n=== Test 5: Performance Benchmark ==="+  +  let sizes = [100, 500, 1000, 5000, 10000]+  +  forM_ sizes $ \n -> do+    putStrLn $ "\nBenchmarking " ++ show n ++ " qubits:"+    +    -- Create tableau+    start <- getCurrentTime+    let tab0 = largeEmpty n+    mid1 <- getCurrentTime+    putStrLn $ "  Creation: " ++ show (diffUTCTime mid1 start)+    +    -- Apply 100 Hadamards+    let tab1 = foldl (\t i -> largeApplyGate (LargeLocal (LargeHadamard (i `mod` n))) t) +                     tab0 [0..99]+    mid2 <- getCurrentTime+    putStrLn $ "  100 Hadamards: " ++ show (diffUTCTime mid2 mid1)+    +    -- Apply 100 CNOTs+    let tab2 = foldl (\t i -> largeApplyGate (LargeCNOT (i `mod` n) ((i+1) `mod` n)) t) +                     tab1 [0..99]+    mid3 <- getCurrentTime+    putStrLn $ "  100 CNOTs: " ++ show (diffUTCTime mid3 mid2)+    +    -- Check validity+    let valid = largeIsValid tab2+    putStrLn $ "  Tableau valid: " ++ show valid+    +    -- Memory estimate: 2 vectors (X and Z) × n qubits × 8 bytes per 64 qubits+    let chunksPerQubit = (n + 63) `div` 64+        bytesPerPauli = 2 * chunksPerQubit * 8  -- X and Z+        totalBytes = 2 * n * bytesPerPauli  -- Stabs and destabs+    putStrLn $ "  Estimated memory: " ++ show (totalBytes `div` 1024) ++ " KB"++-- ============================================================================+-- Main+-- ============================================================================++parseArgs :: [String] -> Config+parseArgs args = foldl parseArg defaultConfig args+  where+    parseArg cfg "--verbose" = cfg { verbose = True }+    parseArg cfg ('-':'-':'b':'e':'l':'l':'-':'p':'a':'i':'r':'s':'=':n) = +      cfg { bellPairsQubits = read n }+    parseArg cfg ('-':'-':'r':'e':'p':'-':'c':'o':'d':'e':'=':n) = +      cfg { repCodeQubits = read n }+    parseArg cfg ('-':'-':'r':'a':'n':'d':'o':'m':'=':n) = +      cfg { randomCircuits = read n }+    parseArg cfg _ = cfg++printUsage :: IO ()+printUsage = do+  putStrLn "Usage: verify-large-tableau [OPTIONS]"+  putStrLn ""+  putStrLn "Options:"+  putStrLn "  --verbose              Enable verbose output"+  putStrLn "  --bell-pairs=N         Qubits for Bell pairs test (default: 10000)"+  putStrLn "  --rep-code=N           Qubits for repetition code test (default: 1000)"+  putStrLn "  --random=N             Number of random circuits (default: 100)"+  putStrLn "  --help                 Show this help"+  putStrLn ""+  putStrLn "Examples:"+  putStrLn "  verify-large-tableau                    # Run all tests with defaults"+  putStrLn "  verify-large-tableau --bell-pairs=5000  # Test with 5000 qubits"+  putStrLn "  verify-large-tableau --verbose          # Verbose output"++main :: IO ()+main = do+  args <- getArgs+  +  when ("--help" `elem` args) $ do+    printUsage+    exitSuccess+  +  let config = parseArgs args+  +  putStrLn "========================================"+  putStrLn "  LargeTableau Verification Suite"+  putStrLn "========================================"+  putStrLn $ "Configuration:"+  putStrLn $ "  Bell pairs test: " ++ show (bellPairsQubits config) ++ " qubits"+  putStrLn $ "  Rep code test: " ++ show (repCodeQubits config) ++ " qubits"+  putStrLn $ "  Random circuits: " ++ show (randomCircuits config)+  +  startTime <- getCurrentTime+  +  -- Run all tests+  results <- sequence+    [ testBellPairs config+    , testRepCode config+    , testPhaseIdentity config+    , testRandomCircuits config+    ]+  +  -- Performance benchmark+  benchmarkGates+  +  endTime <- getCurrentTime+  +  putStrLn "\n========================================"+  putStrLn "  Summary"+  putStrLn "========================================"+  putStrLn $ "Total time: " ++ show (diffUTCTime endTime startTime)+  +  let allPass = and results+  if allPass+    then do+      putStrLn "\n✅ ALL TESTS PASSED"+      exitSuccess+    else do+      putStrLn "\n❌ SOME TESTS FAILED"+      exitFailure

@@ -0,0 +1,104 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE GADTs #-}++-- | Simple BitVec implementation for arbitrary-qubit support.+-- This is a pragmatic alternative to the full type-family approach.+module SymplecticCHP.BitVec+  ( BitVec+  , bvEmpty+  , bvFromWord64+  , bvToWord64+  , bvXor+  , bvAnd+  , bvTestBit+  , bvSetBit+  , bvClearBit+  , bvPopCount+  , bvShow+    -- * Re-export for compatibility+  , Word64+  ) where++import Data.Bits (Bits(..), popCount, xor)+import Data.Word (Word64)+import qualified Data.Vector.Unboxed as V++-- | Arbitrary-length bit vector+data BitVec = BitVec+  { bvBits :: !(V.Vector Word64)  -- ^ Chunked storage+  , bvSize :: !Int                -- ^ Number of bits+  }+  deriving (Eq)++instance Show BitVec where+  show = bvShow++-- | Create empty BitVec for n bits+bvEmpty :: Int -> BitVec+bvEmpty n+  | n <= 0 = error "bvEmpty: size must be positive"+  | otherwise = +      let chunks = (n + 63) `div` 64+      in BitVec (V.replicate chunks 0) n++-- | Create from Word64 (for n <= 64)+bvFromWord64 :: Word64 -> BitVec+bvFromWord64 w = BitVec (V.singleton w) 64++-- | Convert to Word64 (only valid for n <= 64)+bvToWord64 :: BitVec -> Word64+bvToWord64 (BitVec v n)+  | n <= 64 && V.length v >= 1 = V.head v+  | otherwise = error "bvToWord64: BitVec too large"++-- | XOR operation+bvXor :: BitVec -> BitVec -> BitVec+bvXor (BitVec a sa) (BitVec b sb)+  | sa /= sb = error "bvXor: size mismatch"+  | otherwise = BitVec (V.zipWith xor a b) sa++-- | AND operation  +bvAnd :: BitVec -> BitVec -> BitVec+bvAnd (BitVec a sa) (BitVec b sb)+  | sa /= sb = error "bvAnd: size mismatch"+  | otherwise = BitVec (V.zipWith (.&.) a b) sa++-- | Test bit+bvTestBit :: BitVec -> Int -> Bool+bvTestBit (BitVec v n) i+  | i < 0 || i >= n = False  -- Out of bounds+  | otherwise = +      let (word, bit) = i `divMod` 64+      in (V.unsafeIndex v word) `testBit` bit++-- | Set bit+bvSetBit :: BitVec -> Int -> BitVec+bvSetBit bv@(BitVec v n) i+  | i < 0 || i >= n = bv  -- Out of bounds: no change+  | otherwise =+      let (word, bit) = i `divMod` 64+          oldVal = V.unsafeIndex v word+          newVal = oldVal `setBit` bit+      in BitVec (V.unsafeUpd v [(word, newVal)]) n++-- | Clear bit+bvClearBit :: BitVec -> Int -> BitVec+bvClearBit bv@(BitVec v n) i+  | i < 0 || i >= n = bv+  | otherwise =+      let (word, bit) = i `divMod` 64+          oldVal = V.unsafeIndex v word+          newVal = oldVal `clearBit` bit+      in BitVec (V.unsafeUpd v [(word, newVal)]) n++-- | Population count+bvPopCount :: BitVec -> Int+bvPopCount (BitVec v _) = V.sum (V.map popCount v)++-- | Show as string of bits+bvShow :: BitVec -> String+bvShow (BitVec v n) = +  "BitVec[" ++ show n ++ "] " ++ +  concatMap showChunk (V.toList v)+  where+    showChunk w = [if testBit w i then '1' else '0' | i <- [0..63]]

@@ -0,0 +1,217 @@+{-# LANGUAGE GADTs #-}+{-# LANGUAGE RankNTypes #-}++-- | Large tableau support using BitVec for arbitrary qubit counts.+-- This module provides a complement to the standard Word64-based tableau.+module SymplecticCHP.LargeTableau+  ( -- * Large Pauli operators+    LargePauli(..)+  , lpPauliX+  , lpPauliZ+  , lpPauliY+  , lpOmega+  , lpMultiply+    -- * Gates (redefined to avoid circular imports)+  , LargeLocalSymplectic(..)+  , LargeSymplecticGate(..)+  , LargeMeasurementResult(..)+    -- * Large Tableau+  , LargeTableau(..)+  , largeEmpty+  , largeNQubits+  , largeApplyGate+  , largeMeasure+    -- * Validation+  , largeIsValid+  , largeIsDeterminate+  ) where++import Data.Bits (Bits(..), popCount, xor)+import Data.Word (Word64)+import SymplecticCHP.BitVec++import System.Random (randomRIO)+import qualified Data.Vector as V++-- ============================================================================+-- Gate Types (local copies to avoid circular imports)+-- ============================================================================++data LargeLocalSymplectic +  = LargeHadamard !Int+  | LargePhase !Int+  deriving (Show, Eq)++data LargeSymplecticGate+  = LargeLocal !LargeLocalSymplectic+  | LargeCNOT !Int !Int+  deriving (Show, Eq)++data LargeMeasurementResult = LargeDeterminate Bool | LargeRandom Bool+  deriving (Show, Eq)++-- ============================================================================+-- Large Pauli (using BitVec)+-- ============================================================================++-- | Pauli operator for arbitrary qubit count using BitVec+data LargePauli = LargePauli+  { lpX :: !BitVec+  , lpZ :: !BitVec+  , lpPhase :: !Int+  , lpNQubits :: !Int+  }+  deriving (Eq, Show)++-- | Create Pauli X on qubit i+lpPauliX :: Int -> Int -> LargePauli+lpPauliX n i = LargePauli (bvSetBit (bvEmpty n) i) (bvEmpty n) 0 n++-- | Create Pauli Z on qubit i  +lpPauliZ :: Int -> Int -> LargePauli+lpPauliZ n i = LargePauli (bvEmpty n) (bvSetBit (bvEmpty n) i) 0 n++-- | Create Pauli Y on qubit i+lpPauliY :: Int -> Int -> LargePauli+lpPauliY n i = LargePauli (bvSetBit (bvEmpty n) i) (bvSetBit (bvEmpty n) i) 1 n++-- | Symplectic inner product+lpOmega :: LargePauli -> LargePauli -> Bool+lpOmega (LargePauli x1 z1 _ n1) (LargePauli x2 z2 _ n2)+  | n1 /= n2 = error "lpOmega: qubit count mismatch"+  | otherwise = odd (bvPopCount (bvAnd x1 z2) + bvPopCount (bvAnd z1 x2))++-- | Pauli multiplication+lpMultiply :: LargePauli -> LargePauli -> LargePauli+lpMultiply (LargePauli x1 z1 r1 n1) (LargePauli x2 z2 r2 n2)+  | n1 /= n2 = error "lpMultiply: qubit count mismatch"+  | otherwise =+      let x = bvXor x1 x2+          z = bvXor z1 z2+          symPhase = bvPopCount (bvAnd x1 z2) - bvPopCount (bvAnd z1 x2)+          r = (r1 + r2 + symPhase) `mod` 4+      in LargePauli x z r n1++-- ============================================================================+-- Large Tableau+-- ============================================================================++-- | Tableau for arbitrary qubit count+data LargeTableau = LargeTableau+  { ltStabs :: !(V.Vector LargePauli)+  , ltDestabs :: !(V.Vector LargePauli)+  , ltN :: !Int+  }+  deriving (Show)++-- | Create empty tableau (|0...0⟩ state)+largeEmpty :: Int -> LargeTableau+largeEmpty n+  | n <= 0 = error "largeEmpty: n must be positive"+  | otherwise = LargeTableau stabs destabs n+  where+    stabs = V.fromList [lpPauliZ n i | i <- [0..n-1]]+    destabs = V.fromList [lpPauliX n i | i <- [0..n-1]]++-- | Get qubit count+largeNQubits :: LargeTableau -> Int+largeNQubits = ltN++-- | Apply gate to large tableau+largeApplyGate :: LargeSymplecticGate -> LargeTableau -> LargeTableau+largeApplyGate g (LargeTableau s d n) =+  LargeTableau (V.map (lpApplyGate g) s) (V.map (lpApplyGate g) d) n++-- | Apply gate to large Pauli+lpApplyGate :: LargeSymplecticGate -> LargePauli -> LargePauli+lpApplyGate (LargeLocal (LargeHadamard i)) (LargePauli x z r n) =+  let xi = bvTestBit x i+      zi = bvTestBit z i+      x' = if zi then bvSetBit (bvClearBit x i) i else bvClearBit x i+      z' = if xi then bvSetBit (bvClearBit z i) i else bvClearBit z i+      r' = (r + if xi && zi then 2 else 0) `mod` 4+  in LargePauli x' z' r' n++lpApplyGate (LargeLocal (LargePhase i)) (LargePauli x z r n) =+  let xi = bvTestBit x i+      zi = bvTestBit z i+      -- Z' = Z XOR X+      z' = if xi then bvXor z (bvSetBit (bvEmpty n) i) else z+      r' = (r + if xi && not zi then 1 else 0) `mod` 4+  in LargePauli x z' r' n++lpApplyGate (LargeCNOT c t) (LargePauli x z r n) =+  let xc = bvTestBit x c+      zc = bvTestBit z c+      xt = bvTestBit x t+      zt = bvTestBit z t+      -- X'[t] = X[t] XOR X[c]+      x' = if xc then bvXor x (bvSetBit (bvEmpty n) t) else x+      -- Z'[c] = Z[c] XOR Z[t]+      z' = if zt then bvXor z (bvSetBit (bvEmpty n) c) else z+      phaseTerm = if xc && zt then (if xt /= zc then 2 else 0) + 1 else 0+      r' = (r + phaseTerm) `mod` 4+  in LargePauli x' z' r' n++-- | Measurement on large tableau+largeMeasure :: LargeTableau -> LargePauli -> IO (LargeTableau, LargeMeasurementResult)+largeMeasure tab@(LargeTableau s d n) p+  | largeIsDeterminate tab p = do+      let outcome = largeComputePhase tab p+      return (tab, LargeDeterminate outcome)+  | otherwise = do+      case largeFindAntiCommuting tab p of+        Nothing -> error "Internal error in largeMeasure"+        Just j -> do+          let s_j = s V.! j+              -- Update other stabilizers+              newStabs = V.imap (\k s_k ->+                if k == j+                  then p+                  else if lpOmega p s_k+                       then lpMultiply s_k s_j+                       else s_k) s+              -- New destabilizer is old stabilizer+              newDestabs = d V.// [(j, s_j)]+          +          outcome <- randomRIO (0, 1) :: IO Int+          let p' = p { lpPhase = (lpPhase p + if outcome == 0 then 2 else 0) `mod` 4 }+              finalStabs = newStabs V.// [(j, p')]+          +          return (LargeTableau finalStabs newDestabs n, LargeRandom (outcome == 1))++-- | Check if measurement is deterministic+largeIsDeterminate :: LargeTableau -> LargePauli -> Bool+largeIsDeterminate (LargeTableau s _ _) p =+  V.all (\s_i -> not (lpOmega p s_i)) s++-- | Find anticommuting stabilizer+largeFindAntiCommuting :: LargeTableau -> LargePauli -> Maybe Int+largeFindAntiCommuting (LargeTableau s _ _) p =+  V.ifoldl' (\acc i s_i ->+    case acc of+      Just _ -> acc+      Nothing -> if lpOmega p s_i then Just i else Nothing) Nothing s++-- | Compute deterministic measurement outcome+largeComputePhase :: LargeTableau -> LargePauli -> Bool+largeComputePhase (LargeTableau s d _) p =+  let scratch = V.ifoldl' (\acc j d_j ->+        if lpOmega p d_j+          then lpMultiply acc (s V.! j)+          else acc) (LargePauli (bvEmpty (lpNQubits p)) (bvEmpty (lpNQubits p)) 0 (lpNQubits p)) d+      totalPhase = (lpPhase p - lpPhase scratch) `mod` 4+  in totalPhase == 0++-- | Check tableau validity (duality condition)+largeIsValid :: LargeTableau -> Bool+largeIsValid (LargeTableau s d n) =+  -- Check ω(Dᵢ, Sⱼ) = δᵢⱼ+  V.all (\i ->+    V.all (\j ->+      let di = d V.! i+          sj = s V.! j+          omega = if lpOmega di sj then 1 else 0+      in if i == j then omega == 1 else omega == 0+    ) (V.fromList [0..n-1])+  ) (V.fromList [0..n-1])

@@ -28,7 +28,7 @@ -- PVP summary:     +-+------- breaking API changes --                  | | +----- non-breaking API additions --                  | | | +--- code changes with no API change-version:            0.1.0.0+version:            0.1.0.1  -- A short (one-line) description of the package. synopsis:           CHP Clifford simulator using symplectic geometry@@ -88,6 +88,9 @@     import:           warnings     exposed-modules:         SymplecticCHP+        SymplecticCHP.BitVec+        SymplecticCHP.LargeTableau+    other-modules:     build-depends:         , base                >= 4.17    && < 4.22         , vector              >= 0.12    && < 0.14@@ -131,7 +134,7 @@         , filepath            >= 1.4     && < 1.6         , directory           >= 1.3     && < 1.4         , containers          >= 0.6     && < 0.8-        , stim-parser         >= 0.1     && < 0.2+        , stim-parser         >= 0.1     && < 0.3     build-tool-depends:         hspec-discover:hspec-discover  >= 2.10 && < 2.12     default-language: GHC2021@@ -156,7 +159,7 @@     -- Other library packages from which modules are imported.     build-depends:           , base                >= 4.17    && < 4.22-          , stim-parser         >= 0.1     && < 0.2+          , stim-parser         >= 0.1     && < 0.3           , containers          >= 0.6     && < 0.8           , symplectic-chp           , random              >= 1.2     && < 1.3@@ -165,4 +168,16 @@     hs-source-dirs:   app      -- Base language which the package is written in.+    default-language: GHC2021++executable verify-large-tableau+    import:           warnings+    main-is:          VerifyLargeTableau.hs+    build-depends:+          , base                >= 4.17    && < 4.22+          , symplectic-chp+          , vector              >= 0.12    && < 0.14+          , random              >= 1.2     && < 1.3+          , time                >= 1.11    && < 1.14+    hs-source-dirs:   app     default-language: GHC2021