From 80337ef5ab4e05f7a469a9257e59c75cae3756bd Mon Sep 17 00:00:00 2001
From: InfoTCube <infocube.yout@gmail.com>
Date: Mon, 17 Aug 2026 20:03:12 +0200
Subject: [PATCH 1/9] fix: correctness bugs in unitarity check, simulator, copy
 ctor and init
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- IsUnitary compared floats with exact equality, so any gate built from
  1/sqrt(2) was rejected. qc.Custom(hadamard, 0) threw "not unitary" —
  the first thing anyone tries with Custom(). Now compares with tolerance.
- Simulator.Run mutated the caller's gate list via Gates.RemoveAt(0), so a
  second Run(qc, shots) replayed a circuit stripped of its own gates and
  reported 100% |0..0>. Now works on a local copy.
- The copy constructor documented a deep copy but assigned references, so
  gates applied to the copy mutated the original. Now clones.
- Initialize accepted unnormalized states silently, making GetProbabilities
  return probabilities summing to >1. Now validates |a|^2 + |b|^2 == 1.
- Toffoli rendered its second control as a target marker in circuit diagrams.
- Qubit limit was documented as 30, which cannot allocate (2^30 Complex =
  17 GB, and the CLR caps one array at 2 GB). Lowered to a real 26.
- GetStringResult threw on an all-zero count array.
- Renamed ApplayGate -> ApplyGate; dropped two unused usings.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Mhwc4XeGtzM2UQx9sNLb3r
---
 README.md                             | 11 ++++--
 src/Qubit.NET/Math/QuantumMath.cs     | 32 ++++++++---------
 src/Qubit.NET/QuantumCircuit.cs       | 51 ++++++++++++++++++---------
 src/Qubit.NET/Simulation/Simulator.cs | 45 +++++++++++++----------
 src/Qubit.NET/Utilities/Helpers.cs    |  3 +-
 5 files changed, 86 insertions(+), 56 deletions(-)

diff --git a/README.md b/README.md
index 881e63f..df88d7c 100644
--- a/README.md
+++ b/README.md
@@ -4,7 +4,9 @@
 
 # 🧠 C# Quantum Computing Simulation Library
 
-**Qubit.NET** is a lightweight quantum circuit simulation library written in C#. It allows users to simulate quantum circuits up to 30 qubits, initialize qubits, apply common quantum gates, and measure results — all using a classical computer. Perfect for learning, prototyping, or integrating quantum logic into .NET applications.
+**Qubit.NET** is a lightweight quantum circuit simulation library written in C#. It lets you build quantum circuits, initialize qubits, apply common quantum gates, and measure results — all on a classical computer. Perfect for learning, prototyping, or integrating quantum logic into .NET applications.
+
+The state vector holds 2ⁿ complex amplitudes, so memory is the limit: **20 qubits ≈ 16 MB, 24 ≈ 256 MB, 26 ≈ 1 GB** (the hard ceiling, set by the CLR's 2 GB single-array limit).
 
 ---
 
@@ -62,12 +64,15 @@ You can initialize any qubit to one of the predefined basis states:
 qc.Initialize(0, State.Minus);
 ```
 
-or in any custom state
+or in any custom state, given as amplitudes α and β:
 
 ```csharp
-qc.Initialize(0, new Complex(1, 1), new Complex(2, 2));
+// |ψ⟩ = (|0⟩ + i|1⟩) / √2
+qc.Initialize(0, new Complex(1 / Math.Sqrt(2), 0), new Complex(0, 1 / Math.Sqrt(2)));
 ```
 
+> ⚠️ The state must be normalized — `|α|² + |β|² = 1` — or an `ArgumentException` is thrown.
+
 > ⚠️ Initialization can only be done **before any gate is applied** to that qubit.  
 > This is internally tracked using a private `_isQubitModified` array.
 
diff --git a/src/Qubit.NET/Math/QuantumMath.cs b/src/Qubit.NET/Math/QuantumMath.cs
index 7c8084d..3805c95 100644
--- a/src/Qubit.NET/Math/QuantumMath.cs
+++ b/src/Qubit.NET/Math/QuantumMath.cs
@@ -8,6 +8,14 @@ namespace Qubit.NET.Math;
 /// </summary>
 internal static class QuantumMath
 {
+    /// <summary>
+    /// Absolute tolerance used when comparing floating-point results against exact
+    /// values (unitarity checks, state normalization). Loose enough to absorb the
+    /// rounding of irrational gate entries such as 1/sqrt(2), tight enough to reject
+    /// genuinely wrong matrices.
+    /// </summary>
+    internal const double Tolerance = 1e-9;
+
     /// <summary>
     /// Applies a single-qubit gate to a specific qubit in a multi-qubit state vector.
     /// </summary>
@@ -333,27 +341,17 @@ internal static bool IsUnitary(Complex[,] matrix)
             }
         }
 
-        // Step 3: Check if the result is the identity matrix
+        // Step 3: Check if the result is the identity matrix.
+        // Compared with a tolerance, not for exact equality: gates built from 1/sqrt(2)
+        // (Hadamard, CH, SX...) produce 1.0000000000000002 on the diagonal.
         for (int i = 0; i < rows; i++)
         {
             for (int j = 0; j < cols; j++)
             {
-                if (i == j)
-                {
-                    // Diagonal elements must be 1
-                    if (!result[i, j].Equals(Complex.One))
-                    {
-                        return false;
-                    }
-                }
-                else
-                {
-                    // Off-diagonal elements must be 0
-                    if (!result[i, j].Equals(Complex.Zero))
-                    {
-                        return false;
-                    }
-                }
+                Complex expected = i == j ? Complex.One : Complex.Zero;
+
+                if (Complex.Abs(result[i, j] - expected) > Tolerance)
+                    return false;
             }
         }
 
diff --git a/src/Qubit.NET/QuantumCircuit.cs b/src/Qubit.NET/QuantumCircuit.cs
index af640cc..92a6b45 100644
--- a/src/Qubit.NET/QuantumCircuit.cs
+++ b/src/Qubit.NET/QuantumCircuit.cs
@@ -1,7 +1,4 @@
-using System.Collections;
-using System.Collections.ObjectModel;
-using System.Numerics;
-using System.Security.AccessControl;
+using System.Numerics;
 using System.Text;
 using Qubit.NET.Gates;
 using Qubit.NET.Math;
@@ -49,17 +46,26 @@ public class QuantumCircuit
     /// </summary>
     private readonly bool[] _isQubitModified;
 
+    /// <summary>
+    /// Largest supported qubit count. The state vector holds 2^n <see cref="Complex"/>
+    /// values at 16 bytes each, and the CLR caps a single array at 2 GB, so 2^26 elements
+    /// (1 GB) is the hard ceiling. Circuits near this limit need
+    /// <c>gcAllowVeryLargeObjects</c> and several GB of free RAM.
+    /// </summary>
+    public const int MaxQubitCount = 26;
+
     /// <summary>
     /// Initializes quantum circuit with specified number of qubits in 0 state.
     /// </summary>
-    /// <param name="qubitCount">Number of qubits. (0-30] qubits are supported.</param>
+    /// <param name="qubitCount">Number of qubits, from 1 to <see cref="MaxQubitCount"/>.</param>
     public QuantumCircuit(int qubitCount)
     {
         if (qubitCount <= 0)
-            throw new ArgumentException("Qubit count must be positive.");
-        if (qubitCount > 30)
-            throw new AggregateException("Qubit count can be at most 30.");
-        
+            throw new ArgumentOutOfRangeException(nameof(qubitCount), qubitCount, "Qubit count must be positive.");
+        if (qubitCount > MaxQubitCount)
+            throw new ArgumentOutOfRangeException(nameof(qubitCount), qubitCount,
+                $"Qubit count can be at most {MaxQubitCount}.");
+
         QubitCount = qubitCount;
         StateVector = new Complex[1 << qubitCount];
         StateVector[0] = new Complex(1, 0);
@@ -69,16 +75,19 @@ public QuantumCircuit(int qubitCount)
     /// <summary>
     /// Creates a new instance of the <see cref="QuantumCircuit"/> class by copying the properties of the provided <paramref name="qc"/> object.
     /// Initializes the new quantum circuit with the same number of qubits and state vector as the original.
-    /// This constructor creates a deep copy, ensuring the new instance is independent of the original.
+    /// The state vector, gate list and initialization list are copied, so applying gates to
+    /// either circuit afterwards leaves the other untouched. The recorded <see cref="Gate"/>
+    /// entries themselves are shared, which is safe because they are never mutated after being appended.
     /// </summary>
     /// <param name="qc">The <see cref="QuantumCircuit"/> object to copy.</param>
     public QuantumCircuit(QuantumCircuit qc)
     {
         QubitCount = qc.QubitCount;
-        StateVector = qc.StateVector;
-        Gates = qc.Gates;
-        Initializations = qc.Initializations;
-        _isQubitModified = qc._isQubitModified;
+        StateVector = (Complex[])qc.StateVector.Clone();
+        Gates = new List<Gate>(qc.Gates);
+        Initializations = new List<InitialState>(qc.Initializations);
+        _isQubitModified = (bool[])qc._isQubitModified.Clone();
+        RandomSource = qc.RandomSource;
     }
 
     /// <summary>
@@ -164,17 +173,27 @@ public void Initialize(int qubit, State state)
     /// <param name="qubit">The index of the qubit to initialize.</param>
     /// <param name="alpha">Amplitude for the |0⟩ component of the qubit.</param>
     /// <param name="beta">Amplitude for the |1⟩ component of the qubit.</param>
-    /// /// <param name="state">Optional description of the state.</param>
+    /// <param name="state">Optional description of the state.</param>
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
+    /// <exception cref="ArgumentException">
+    /// Thrown if the state is not normalized, i.e. |α|² + |β|² ≠ 1.
+    /// </exception>
     public void Initialize(int qubit, Complex alpha, Complex beta, State state = State.Custom)
     {
         CheckQubit(qubit);
-        
+
         if (_isQubitModified[qubit])
             throw new InvalidOperationException($"Qubit {qubit} has already been modified and cannot be re-initialized.");
 
+        double normSquared = alpha.Real * alpha.Real + alpha.Imaginary * alpha.Imaginary
+                             + beta.Real * beta.Real + beta.Imaginary * beta.Imaginary;
+
+        if (System.Math.Abs(normSquared - 1.0) > QuantumMath.Tolerance)
+            throw new ArgumentException(
+                $"A qubit state must be normalized: |alpha|^2 + |beta|^2 must equal 1, but was {normSquared}.");
+
         StateVector = QuantumMath.InitializeState(StateVector, qubit, alpha, beta);
         
         _isQubitModified[qubit] = true;
diff --git a/src/Qubit.NET/Simulation/Simulator.cs b/src/Qubit.NET/Simulation/Simulator.cs
index 061f982..e3c9c18 100644
--- a/src/Qubit.NET/Simulation/Simulator.cs
+++ b/src/Qubit.NET/Simulation/Simulator.cs
@@ -34,28 +34,30 @@ public static class Simulator
             stateVector = QuantumMath.InitializeState(stateVector, init.QubitIndex, init.Alpha, init.Beta);
         }
 
-        GateType currentGateType = qc.Gates.First().GateType;
+        // Work on a local copy: Run must never mutate the circuit it is handed, or a
+        // second Run(qc, shots) would replay a circuit stripped of its own gates.
+        List<Gate> remainingGates = qc.Gates.ToList();
 
-        while (currentGateType != GateType.Measure && qc.Gates.Count != 0)
+        // Apply the deterministic prefix (everything before the first measurement) once,
+        // then replay only the rest per shot.
+        while (remainingGates.Count > 0 && remainingGates[0].GateType != GateType.Measure)
         {
-            Gate currentGate = qc.Gates.First();
-            
-            stateVector = ApplayGate(stateVector, currentGate, currentGateType);
+            Gate currentGate = remainingGates[0];
 
-            qc.Gates.RemoveAt(0);
-            
-            currentGateType = qc.Gates.Count > 0 ? qc.Gates.First().GateType : currentGateType;
+            stateVector = ApplyGate(stateVector, currentGate, currentGate.GateType);
+
+            remainingGates.RemoveAt(0);
         }
-        
+
         IList<(int[], int)> results = new List<(int[], int)>();
-        
+
         for (int i = 0; i < shots; i++)
         {
             Complex[] modStateVector = (Complex[])stateVector.Clone();
 
             int measurmentNumber = 0;
-            
-            foreach (var gate in qc.Gates)
+
+            foreach (var gate in remainingGates)
             {
                 if (gate.GateType == GateType.Measure)
                 {
@@ -70,7 +72,7 @@ public static class Simulator
                 }
                 else
                 {
-                    modStateVector = ApplayGate(modStateVector, gate, gate.GateType);
+                    modStateVector = ApplyGate(modStateVector, gate, gate.GateType);
                 }
             }
         }
@@ -90,21 +92,26 @@ public static string GetStringResult(this (int[], int) result)
     {
         StringBuilder sb = new StringBuilder();
         sb.Append("{");
+
+        bool any = false;
+
         for (int i = 0; i < result.Item1.Length; i++)
         {
             if(result.Item1[i] == 0)
                 continue;
-            
+
+            if (any) sb.Append(", ");
+
             sb.Append("'");
             sb.Append(Convert.ToString(i, 2).PadLeft(result.Item2, '0'));
             sb.Append("': ");
             sb.Append(result.Item1[i].ToString());
-            sb.Append(", ");
+
+            any = true;
         }
-        
-        sb.Remove(sb.Length - 2, 2);
+
         sb.Append("}");
-        
+
         return sb.ToString();
     }
 
@@ -115,7 +122,7 @@ public static string GetStringResult(this (int[], int) result)
     /// <param name="stateVector">The current quantum state vector.</param>
     /// <param name="currentGate">The gate to apply.</param>
     /// <param name="gateType">The type of gate being applied.</param>
-    private static Complex[] ApplayGate(Complex[] stateVector, Gate currentGate, GateType gateType)
+    private static Complex[] ApplyGate(Complex[] stateVector, Gate currentGate, GateType gateType)
     {
         switch (gateType)
         {
diff --git a/src/Qubit.NET/Utilities/Helpers.cs b/src/Qubit.NET/Utilities/Helpers.cs
index 12bc843..fd7680b 100644
--- a/src/Qubit.NET/Utilities/Helpers.cs
+++ b/src/Qubit.NET/Utilities/Helpers.cs
@@ -68,7 +68,8 @@ internal static string[] GateTypeToCharRepresentation(GateType gateType)
             GateType.CRz => ["@", "Rz"],
             GateType.CU3 => ["@", "U3"],
             GateType.SWAP => ["X", "X"],
-            GateType.Toffoli => ["@", "+", "+"],
+            // Order matches the drawer's controls-then-targets iteration: two controls, one target.
+            GateType.Toffoli => ["@", "@", "+"],
             GateType.Fredkin => ["@", "X", "X"],
             GateType.Measure => ["M"],
             GateType.Custom => ["C"],

From 1374678029da9ffd9439bd6a557fe5e02c1698ee Mon Sep 17 00:00:00 2001
From: InfoTCube <infocube.yout@gmail.com>
Date: Mon, 17 Aug 2026 20:10:26 +0200
Subject: [PATCH 2/9] feat: ship as a NuGet package, multi-targeted for Unity
 and modern .NET
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The library could not previously be consumed by anyone: the README told
people to clone the repo and copy .cs files, and the net8.0-only target
meant Unity — the flagship use case — could not import it at all.

- Multi-target netstandard2.1 (Unity 2021.2+), net8.0 and net10.0 from a
  single package; NuGet picks the right one per consumer. Replaced
  BitOperations.PopCount, which netstandard2.1 lacks, with a plain
  power-of-two test.
- Full package metadata: id, license, icon, README, tags, XML docs,
  Source Link and a symbol package.
- Added CI, plus a tag-driven workflow that packs and pushes to NuGet.
- Console-free rendering: ToDiagram() and QuantumGates.Format() return
  strings, so the drawer works under Unity, ASP.NET and tests. Draw() and
  Print() keep their console behaviour and colors by wrapping them.
- Moved Examples.Example into the public API as Circuits.BellStates, so the
  library no longer ships a namespace called Examples.
- Documented every public gate matrix and the State enum, clearing 82
  undocumented-member warnings now that XML docs are generated.
- Spelled out Enumerable.Reverse: under C# 14 first-class spans,
  array.Reverse() binds to MemoryExtensions.Reverse and returns void.
- Added a CHANGELOG.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Mhwc4XeGtzM2UQx9sNLb3r
---
 .github/workflows/ci.yml              |  39 +++++++
 .github/workflows/release.yml         |  48 ++++++++
 .gitignore                            |   5 +-
 CHANGELOG.md                          |  44 ++++++++
 README.md                             |  41 +++++--
 src/Qubit.NET/Circuits/BellStates.cs  |  76 +++++++++++++
 src/Qubit.NET/Examples/Example.cs     |  52 ---------
 src/Qubit.NET/Gates/QuantumGates.cs   | 153 +++++++++++++++++++++++---
 src/Qubit.NET/Math/QuantumMath.cs     |   4 +-
 src/Qubit.NET/QuantumCircuit.cs       |   4 +-
 src/Qubit.NET/QuantumCircuitDrawer.cs | 143 +++++++++++++++---------
 src/Qubit.NET/Qubit.NET.csproj        |  40 ++++++-
 src/Qubit.NET/Simulation/Simulator.cs |   2 +-
 src/Qubit.NET/Utilities/State.cs      |  14 ++-
 14 files changed, 534 insertions(+), 131 deletions(-)
 create mode 100644 .github/workflows/ci.yml
 create mode 100644 .github/workflows/release.yml
 create mode 100644 CHANGELOG.md
 create mode 100644 src/Qubit.NET/Circuits/BellStates.cs
 delete mode 100644 src/Qubit.NET/Examples/Example.cs

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..fd341f7
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,39 @@
+name: CI
+
+on:
+  push:
+    branches: [main]
+  pull_request:
+    branches: [main]
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Setup .NET
+        uses: actions/setup-dotnet@v4
+        with:
+          dotnet-version: |
+            8.0.x
+            10.0.x
+
+      - name: Restore
+        run: dotnet restore
+
+      - name: Build
+        run: dotnet build --no-restore --configuration Release
+
+      - name: Test
+        run: dotnet test --no-build --configuration Release --verbosity normal
+
+      - name: Pack
+        run: dotnet pack src/Qubit.NET/Qubit.NET.csproj --no-build --configuration Release --output ./artifacts
+
+      - name: Upload package
+        uses: actions/upload-artifact@v4
+        with:
+          name: nupkg
+          path: ./artifacts/*.nupkg
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
new file mode 100644
index 0000000..b7502b8
--- /dev/null
+++ b/.github/workflows/release.yml
@@ -0,0 +1,48 @@
+name: Release
+
+on:
+  push:
+    tags: ['v*']
+
+jobs:
+  publish:
+    runs-on: ubuntu-latest
+
+    steps:
+      - uses: actions/checkout@v4
+        with:
+          fetch-depth: 0   # Source Link needs full history
+
+      - name: Setup .NET
+        uses: actions/setup-dotnet@v4
+        with:
+          dotnet-version: |
+            8.0.x
+            10.0.x
+
+      - name: Derive version from tag
+        run: echo "VERSION=${GITHUB_REF_NAME#v}" >> $GITHUB_ENV
+
+      - name: Build
+        run: dotnet build --configuration Release -p:Version=${{ env.VERSION }}
+
+      - name: Test
+        run: dotnet test --no-build --configuration Release
+
+      - name: Pack
+        run: >
+          dotnet pack src/Qubit.NET/Qubit.NET.csproj --no-build --configuration Release
+          -p:Version=${{ env.VERSION }} --output ./artifacts
+
+      - name: Push to NuGet
+        run: >
+          dotnet nuget push "./artifacts/*.nupkg"
+          --api-key ${{ secrets.NUGET_API_KEY }}
+          --source https://api.nuget.org/v3/index.json
+          --skip-duplicate
+
+      - name: Create GitHub release
+        uses: softprops/action-gh-release@v2
+        with:
+          files: ./artifacts/*.nupkg
+          generate_release_notes: true
diff --git a/.gitignore b/.gitignore
index d48c759..dbab0a1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,5 @@
 .idea
-.vscode
\ No newline at end of file
+.vscode
+artifacts/
+bin/
+obj/
diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 0000000..fc5363e
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,44 @@
+# Changelog
+
+All notable changes to this project are documented here.
+The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
+and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+
+## [1.0.0] — unreleased
+
+First packaged release. Qubit.NET is now installable with `dotnet add package Qubit.NET`
+and usable from Unity.
+
+### Added
+
+- **NuGet package** with Source Link, symbol package, XML documentation and icon.
+- **Multi-targeting**: `netstandard2.1` (Unity), `net8.0` and `net10.0`.
+- `QuantumCircuit.ToDiagram()` — returns the ASCII circuit diagram as a string, so the
+  drawer works without a console (Unity, ASP.NET, tests). `Draw()` still prints in color.
+- `QuantumGates.Format(matrix)` — the string-returning counterpart of `Print`.
+- `QuantumCircuit.MaxQubitCount` constant.
+- Continuous integration and a tag-driven NuGet release workflow.
+- XML documentation for every public gate matrix and for the `State` enum.
+
+### Changed
+
+- **Breaking**: `Examples.Example` is now `Circuits.BellStates`.
+- **Breaking**: the maximum qubit count is 26, down from a documented 30 that could never
+  actually allocate (2³⁰ `Complex` values is 17 GB, and the CLR caps a single array at 2 GB).
+- **Breaking**: `Initialize` now rejects unnormalized states instead of silently producing
+  a non-physical one.
+- **Breaking**: constructing a circuit with an invalid qubit count throws
+  `ArgumentOutOfRangeException` rather than `ArgumentException`/`AggregateException`.
+
+### Fixed
+
+- `IsUnitary` compared floating-point values for exact equality, so every matrix built
+  from `1/√2` was rejected — `qc.Custom(hadamard, 0)` threw "The provided matrix is not
+  unitary." Comparisons now use a tolerance.
+- `Simulator.Run` mutated the gate list of the circuit it was given, so calling it twice
+  on the same circuit returned wrong results the second time.
+- The copy constructor documented a deep copy but shared the gate list and modification
+  flags with the original, so gates applied to a copy also affected the source circuit.
+- Toffoli gates rendered their second control as a target marker in circuit diagrams.
+- `GetStringResult` threw an exception when every measurement count was zero.
+- Renamed the internal `ApplayGate` to `ApplyGate`.
diff --git a/README.md b/README.md
index df88d7c..579231c 100644
--- a/README.md
+++ b/README.md
@@ -4,25 +4,37 @@
 
 # 🧠 C# Quantum Computing Simulation Library
 
+[![NuGet](https://img.shields.io/nuget/v/Qubit.NET.svg)](https://www.nuget.org/packages/Qubit.NET/)
+[![Downloads](https://img.shields.io/nuget/dt/Qubit.NET.svg)](https://www.nuget.org/packages/Qubit.NET/)
+[![CI](https://github.com/InfoTCube/Qubit.NET/actions/workflows/ci.yml/badge.svg)](https://github.com/InfoTCube/Qubit.NET/actions/workflows/ci.yml)
+[![License: MIT](https://img.shields.io/badge/License-MIT-blue.svg)](LICENSE)
+
 **Qubit.NET** is a lightweight quantum circuit simulation library written in C#. It lets you build quantum circuits, initialize qubits, apply common quantum gates, and measure results — all on a classical computer. Perfect for learning, prototyping, or integrating quantum logic into .NET applications.
 
 The state vector holds 2ⁿ complex amplitudes, so memory is the limit: **20 qubits ≈ 16 MB, 24 ≈ 256 MB, 26 ≈ 1 GB** (the hard ceiling, set by the CLR's 2 GB single-array limit).
 
 ---
 
-### ✅ Requirements
-- .NET 6.0 or newer
-- `System.Numerics` (for complex numbers — included in .NET)
-
-### 📥 Setup
-Clone or download the repository:
+### 📥 Install
 
 ```bash
-git clone https://github.com/InfoTCube/Qubit.Net.git
-cd Qubit.NET
+dotnet add package Qubit.NET
 ```
 
-Add the project to your solution or include the `.cs` files (`QuantumCircuit.cs`, `QuantumGates.cs`, etc.) in your C# project.
+Zero dependencies. Targets **.NET Standard 2.1**, **.NET 8** and **.NET 10**.
+
+### 🎮 Unity
+
+Qubit.NET ships a `netstandard2.1` build, so it works in Unity 2021.2+. Either install it
+through [NuGetForUnity](https://github.com/GlitchEnzo/NuGetForUnity), or drop
+`lib/netstandard2.1/Qubit.NET.dll` from the package into `Assets/Plugins/`.
+
+Unity has no `Console`, so use the string-returning APIs:
+
+```csharp
+Debug.Log(qc.ToDiagram());                     // instead of qc.Draw()
+Debug.Log(QuantumGates.Format(QuantumGates.H)); // instead of QuantumGates.Print(...)
+```
 
 ---
 
@@ -47,6 +59,17 @@ qc.Draw();
 Console.WriteLine($"Measured: {qc.Measure()}"); // Possible: 00 or 11
 ```
 
+`Draw()` prints a colored ASCII diagram to the console — `ToDiagram()` returns the same
+thing as a string:
+
+```
+q2 (0): ───────────[+]──[X]──[M]─
+                    |    |    |
+q1 (0): ──────[+]───@────|───[M]─
+               |    |    |    |
+q0 (0): ─[H]───@────@───[X]──[M]─
+```
+
 ---
 
 ## 🧰 Features
diff --git a/src/Qubit.NET/Circuits/BellStates.cs b/src/Qubit.NET/Circuits/BellStates.cs
new file mode 100644
index 0000000..f2b0617
--- /dev/null
+++ b/src/Qubit.NET/Circuits/BellStates.cs
@@ -0,0 +1,76 @@
+namespace Qubit.NET.Circuits;
+
+/// <summary>
+/// Ready-made circuits preparing the four maximally entangled two-qubit Bell states
+/// and the three-qubit GHZ state.
+/// </summary>
+public static class BellStates
+{
+    /// <summary>
+    /// Prepares |Φ⁺⟩ = (|00⟩ + |11⟩)/√2.
+    /// </summary>
+    /// <returns>A two-qubit circuit in the |Φ⁺⟩ state.</returns>
+    public static QuantumCircuit PhiPlus()
+    {
+        QuantumCircuit phiPlus = new QuantumCircuit(2);
+
+        phiPlus.H(0);
+        phiPlus.CNOT(0, 1);
+
+        return phiPlus;
+    }
+
+    /// <summary>
+    /// Prepares |Φ⁻⟩ = (|00⟩ − |11⟩)/√2.
+    /// </summary>
+    /// <returns>A two-qubit circuit in the |Φ⁻⟩ state.</returns>
+    public static QuantumCircuit PhiMinus()
+    {
+        QuantumCircuit phiMinus = PhiPlus();
+
+        phiMinus.Z(1);
+
+        return phiMinus;
+    }
+
+    /// <summary>
+    /// Prepares |Ψ⁺⟩ = (|01⟩ + |10⟩)/√2.
+    /// </summary>
+    /// <returns>A two-qubit circuit in the |Ψ⁺⟩ state.</returns>
+    public static QuantumCircuit PsiPlus()
+    {
+        QuantumCircuit psiPlus = PhiPlus();
+
+        psiPlus.X(1);
+
+        return psiPlus;
+    }
+
+    /// <summary>
+    /// Prepares |Ψ⁻⟩ = (|01⟩ − |10⟩)/√2.
+    /// </summary>
+    /// <returns>A two-qubit circuit in the |Ψ⁻⟩ state.</returns>
+    public static QuantumCircuit PsiMinus()
+    {
+        QuantumCircuit psiMinus = PsiPlus();
+
+        psiMinus.Z(1);
+
+        return psiMinus;
+    }
+
+    /// <summary>
+    /// Prepares the three-qubit GHZ state, (|000⟩ + |111⟩)/√2.
+    /// </summary>
+    /// <returns>A three-qubit circuit in the GHZ state.</returns>
+    public static QuantumCircuit GHZ()
+    {
+        QuantumCircuit ghz = new QuantumCircuit(3);
+
+        ghz.H(0);
+        ghz.CNOT(0, 1);
+        ghz.CNOT(0, 2);
+
+        return ghz;
+    }
+}
diff --git a/src/Qubit.NET/Examples/Example.cs b/src/Qubit.NET/Examples/Example.cs
deleted file mode 100644
index 6af9aff..0000000
--- a/src/Qubit.NET/Examples/Example.cs
+++ /dev/null
@@ -1,52 +0,0 @@
-namespace Qubit.NET.Examples;
-
-public static class Example
-{
-    public static QuantumCircuit PhiPlus()
-    {
-        QuantumCircuit phiPlus = new QuantumCircuit(2);
-        
-        phiPlus.H(0);
-        phiPlus.CNOT(0, 1);
-        
-        return phiPlus;
-    }
-
-    public static QuantumCircuit PhiMinus()
-    {
-        QuantumCircuit phiMinus = PhiPlus();
-        
-        phiMinus.Z(1);
-        
-        return phiMinus;
-    }
-    
-    public static QuantumCircuit PsiPlus()
-    {
-        QuantumCircuit psiPlus = PhiPlus();
-        
-        psiPlus.X(1);
-        
-        return psiPlus;
-    }
-    
-    public static QuantumCircuit PsiMinus()
-    {
-        QuantumCircuit psiMinus = PsiPlus();
-        
-        psiMinus.Z(1);
-        
-        return psiMinus;
-    }
-
-    public static QuantumCircuit GHZ()
-    {
-        QuantumCircuit ghz = new QuantumCircuit(3);
-        
-        ghz.H(0);
-        ghz.CNOT(0, 1);
-        ghz.CNOT(0, 2);
-
-        return ghz;
-    }
-}
\ No newline at end of file
diff --git a/src/Qubit.NET/Gates/QuantumGates.cs b/src/Qubit.NET/Gates/QuantumGates.cs
index b4b196b..cad2e9b 100644
--- a/src/Qubit.NET/Gates/QuantumGates.cs
+++ b/src/Qubit.NET/Gates/QuantumGates.cs
@@ -1,5 +1,6 @@
 using System;
 using System.Numerics;
+using System.Text;
 using Qubit.NET.Utilities;
 
 namespace Qubit.NET.Gates;
@@ -13,60 +14,92 @@ public static class QuantumGates
     private static readonly Complex TElement = new Complex(System.Math.Cos(System.Math.PI / 4), 
         System.Math.Sin(System.Math.PI / 4));
     
+    /// <summary>
+    /// The identity gate I. Leaves the qubit unchanged.
+    /// </summary>
     public static Complex[,] I => new Complex[,]
     {
         { 1, 0 },
         { 0, 1 }
     };
     
+    /// <summary>
+    /// The Hadamard gate H. Maps |0⟩ to |+⟩ and |1⟩ to |−⟩, creating an equal superposition.
+    /// </summary>
     public static Complex[,] H => new Complex[,]
     {
         { InvSqrt2, InvSqrt2 },
         { InvSqrt2, -InvSqrt2 }
     };
 
+    /// <summary>
+    /// The Pauli-X (NOT) gate. Flips |0⟩ and |1⟩.
+    /// </summary>
     public static Complex[,] X => new Complex[,]
     {
         { 0, 1 },
         { 1, 0 }
     };
     
+    /// <summary>
+    /// The Pauli-Y gate. A bit flip combined with a phase flip.
+    /// </summary>
     public static Complex[,] Y => new Complex[,]
     {
         { 0, -Complex.ImaginaryOne },
         { Complex.ImaginaryOne, 0 }
     };
 
+    /// <summary>
+    /// The Pauli-Z gate. Leaves |0⟩ alone and maps |1⟩ to −|1⟩.
+    /// </summary>
     public static Complex[,] Z => new Complex[,]
     {
         { 1, 0 },
         { 0, -1 }
     };
 
+    /// <summary>
+    /// The S (phase) gate. Applies a π/2 phase to |1⟩.
+    /// </summary>
     public static Complex[,] S => new Complex[,]
     {
         { 1, 0 },
         { 0, Complex.ImaginaryOne }
     };
 
+    /// <summary>
+    /// The S† gate, the inverse of <see cref="S"/>. Applies a −π/2 phase to |1⟩.
+    /// </summary>
     public static Complex[,] Sdag => new Complex[,]
     {
         { 1, 0 },
         { 0, -Complex.ImaginaryOne }
     };
 
+    /// <summary>
+    /// The T gate. Applies a π/4 phase to |1⟩; the fourth root of Z.
+    /// </summary>
     public static Complex[,] T => new Complex[,]
     {
         { 1, 0 },
         { 0, TElement }
     };
 
+    /// <summary>
+    /// The T† gate, the inverse of <see cref="T"/>. Applies a −π/4 phase to |1⟩.
+    /// </summary>
     public static Complex[,] Tdag => new Complex[,]
     {
         { 1, 0 },
         { 0, Complex.Conjugate(TElement) }
     };
 
+    /// <summary>
+    /// Builds a rotation of <paramref name="theta"/> radians about the X axis of the Bloch sphere.
+    /// </summary>
+    /// <param name="theta">The rotation angle in radians.</param>
+    /// <returns>The gate matrix.</returns>
     public static Complex[,] Rx(double theta)
     {
         Complex cosTheta = Complex.Cos(theta / 2);
@@ -79,6 +112,11 @@ public static class QuantumGates
         };
     }
 
+    /// <summary>
+    /// Builds a rotation of <paramref name="theta"/> radians about the Y axis of the Bloch sphere.
+    /// </summary>
+    /// <param name="theta">The rotation angle in radians.</param>
+    /// <returns>The gate matrix.</returns>
     public static Complex[,] Ry(double theta)
     {
         Complex cosTheta = Complex.Cos(theta / 2);
@@ -91,6 +129,11 @@ public static class QuantumGates
         };
     }
 
+    /// <summary>
+    /// Builds a rotation of <paramref name="theta"/> radians about the Z axis of the Bloch sphere.
+    /// </summary>
+    /// <param name="theta">The rotation angle in radians.</param>
+    /// <returns>The gate matrix.</returns>
     public static Complex[,] Rz(double theta)
     {
         Complex expNeg = Complex.Exp(-Complex.ImaginaryOne * (theta / 2));
@@ -103,20 +146,36 @@ public static class QuantumGates
         };
     }
     
+    /// <summary>
+    /// The √X gate. Applied twice it equals <see cref="X"/>.
+    /// </summary>
     public static Complex[,] SX => new Complex[,]
     {
         { (1 + Complex.ImaginaryOne) / 2, (1 - Complex.ImaginaryOne) / 2 },
         { (1 - Complex.ImaginaryOne) / 2, (1 + Complex.ImaginaryOne) / 2 }
     };
     
+    /// <summary>
+    /// The √Y gate. Applied twice it equals <see cref="Y"/>.
+    /// </summary>
     public static Complex[,] SY => new Complex[,]
     {
         { (1 + Complex.ImaginaryOne) / 2, (-1 - Complex.ImaginaryOne) / 2 },
         { (1 + Complex.ImaginaryOne) / 2, (1 + Complex.ImaginaryOne) / 2 }
     };
 
+    /// <summary>
+    /// The √Z gate, identical to <see cref="S"/>. Applied twice it equals <see cref="Z"/>.
+    /// </summary>
     public static Complex[,] SZ => S;
 
+    /// <summary>
+    /// Builds the general single-qubit rotation U3, which can express any single-qubit unitary.
+    /// </summary>
+    /// <param name="theta">The rotation angle in radians.</param>
+    /// <param name="phi">The phase angle in radians applied before the rotation.</param>
+    /// <param name="lambda">The phase angle in radians applied after the rotation.</param>
+    /// <returns>The gate matrix.</returns>
     public static Complex[,] U3(double theta, double phi, double lambda)
     {
         Complex cosTheta = Complex.Cos(theta / 2);
@@ -132,6 +191,9 @@ public static class QuantumGates
         };
     }
     
+    /// <summary>
+    /// The controlled-NOT (CX) gate. Applies <see cref="X"/> to the target when the control is |1⟩.
+    /// </summary>
     public static Complex[,] CNOT => new Complex[,]
     {
         { 1, 0, 0, 0 },
@@ -140,6 +202,9 @@ public static class QuantumGates
         { 0, 0, 1, 0 }
     };
     
+    /// <summary>
+    /// The controlled-Z gate. Applies <see cref="Z"/> to the target when the control is |1⟩.
+    /// </summary>
     public static Complex[,] CZ => new Complex[,]
     {
         { 1, 0, 0, 0 },
@@ -148,6 +213,9 @@ public static class QuantumGates
         { 0, 0, 0, -1 }
     };
     
+    /// <summary>
+    /// The controlled-Y gate. Applies <see cref="Y"/> to the target when the control is |1⟩.
+    /// </summary>
     public static Complex[,] CY => new Complex[,]
     {
         { 1, 0, 0, 0 },
@@ -156,6 +224,9 @@ public static class QuantumGates
         { 0, 0, Complex.ImaginaryOne, 0 }
     };
     
+    /// <summary>
+    /// The controlled-Hadamard gate. Applies <see cref="H"/> to the target when the control is |1⟩.
+    /// </summary>
     public static Complex[,] CH => new Complex[,]
     {
         { 1, 0, 0, 0 },
@@ -164,6 +235,11 @@ public static class QuantumGates
         { 0, 0, InvSqrt2, -InvSqrt2 }
     };
     
+    /// <summary>
+    /// Builds a controlled <see cref="Rx"/> gate.
+    /// </summary>
+    /// <param name="theta">The rotation angle in radians.</param>
+    /// <returns>The gate matrix.</returns>
     public static Complex[,] CRx(double theta)
     {
         Complex cosTheta = Complex.Cos(theta / 2);
@@ -178,6 +254,11 @@ public static class QuantumGates
         };
     }
     
+    /// <summary>
+    /// Builds a controlled <see cref="Ry"/> gate.
+    /// </summary>
+    /// <param name="theta">The rotation angle in radians.</param>
+    /// <returns>The gate matrix.</returns>
     public static Complex[,] CRy(double theta)
     {
         Complex cosTheta = Complex.Cos(theta / 2);
@@ -192,6 +273,11 @@ public static class QuantumGates
         };
     }
     
+    /// <summary>
+    /// Builds a controlled <see cref="Rz"/> gate.
+    /// </summary>
+    /// <param name="theta">The rotation angle in radians.</param>
+    /// <returns>The gate matrix.</returns>
     public static Complex[,] CRz(double theta)
     {
         Complex expNeg = Complex.Exp(-Complex.ImaginaryOne * (theta / 2));
@@ -206,6 +292,13 @@ public static class QuantumGates
         };
     }
     
+    /// <summary>
+    /// Builds a controlled <see cref="U3"/> gate.
+    /// </summary>
+    /// <param name="theta">The rotation angle in radians.</param>
+    /// <param name="phi">The phase angle in radians applied before the rotation.</param>
+    /// <param name="lambda">The phase angle in radians applied after the rotation.</param>
+    /// <returns>The gate matrix.</returns>
     public static Complex[,] CU3(double theta, double phi, double lambda)
     {
         Complex cosTheta = Complex.Cos(theta / 2);
@@ -223,6 +316,9 @@ public static class QuantumGates
         };
     }
 
+    /// <summary>
+    /// The SWAP gate. Exchanges the states of two qubits.
+    /// </summary>
     public static Complex[,] SWAP => new Complex[,]
     {
         { 1, 0, 0, 0 },
@@ -231,6 +327,9 @@ public static class QuantumGates
         { 0, 0, 0, 1 }
     };
 
+    /// <summary>
+    /// The Toffoli (CCX) gate. Flips the target when both controls are |1⟩.
+    /// </summary>
     public static Complex[,] Toffoli => new Complex[,]
     {
         { 1, 0, 0, 0, 0, 0, 0, 0 },
@@ -243,6 +342,9 @@ public static class QuantumGates
         { 0, 0, 0, 0, 0, 0, 1, 0 }
     };
     
+    /// <summary>
+    /// The Fredkin (CSWAP) gate. Swaps the two targets when the control is |1⟩.
+    /// </summary>
     public static Complex[,] Fredkin => new Complex[,]
     {
         { 1, 0, 0, 0, 0, 0, 0, 0 },
@@ -255,38 +357,61 @@ public static class QuantumGates
         { 0, 0, 0, 0, 0, 0, 0, 1 }
     };
     
+    /// <summary>
+    /// Writes a bracketed, column-aligned rendering of a matrix to the console.
+    /// </summary>
+    /// <param name="matrix">The matrix to print.</param>
+    /// <remarks>
+    /// Requires a console. In environments without one (Unity, ASP.NET, tests),
+    /// use <see cref="Format"/> and write the string wherever you need it.
+    /// </remarks>
     public static void Print(Complex[,] matrix)
+    {
+        Console.WriteLine(Format(matrix));
+    }
+
+    /// <summary>
+    /// Builds a bracketed, column-aligned rendering of a matrix and returns it as a string.
+    /// </summary>
+    /// <param name="matrix">The matrix to render.</param>
+    /// <returns>A multi-line string containing the formatted matrix.</returns>
+    /// <example>
+    /// <code>
+    /// Debug.Log(QuantumGates.Format(QuantumGates.H));   // Unity
+    /// </code>
+    /// </example>
+    public static string Format(Complex[,] matrix)
     {
         int rows = matrix.GetLength(0);
         int cols = matrix.GetLength(1);
-        
+
         int maxLength = 0;
         foreach (var elem in matrix)
         {
             string formatted = Helpers.FormatComplex(elem.Real, elem.Imaginary);
             maxLength = System.Math.Max(maxLength, formatted.Length);
         }
-        
-        Console.Write("\u250c ");
-        Console.Write(new string(' ', (maxLength + 1) * cols));
-        Console.WriteLine("\u2510");
-        
+
+        StringBuilder sb = new StringBuilder();
+
+        sb.Append("\u250c ").Append(' ', (maxLength + 1) * cols).AppendLine("\u2510");
+
         for (int i = 0; i < rows; i++)
         {
-            Console.Write("\u2502 ");
+            sb.Append("\u2502 ");
             for (int j = 0; j < cols; j++)
             {
                 string representation = Helpers.FormatComplex(matrix[i, j].Real, matrix[i, j].Imaginary);
 
                 if (representation == string.Empty) representation = "0";
-                
-                Console.Write(representation.PadRight(maxLength + 1));
+
+                sb.Append(representation.PadRight(maxLength + 1));
             }
-            Console.WriteLine("\u2502");
+            sb.AppendLine("\u2502");
         }
-        
-        Console.Write("\u2514 ");
-        Console.Write(new string(' ', (maxLength + 1) * cols));
-        Console.WriteLine("\u2518");
+
+        sb.Append("\u2514 ").Append(' ', (maxLength + 1) * cols).Append("\u2518");
+
+        return sb.ToString();
     }
 }
\ No newline at end of file
diff --git a/src/Qubit.NET/Math/QuantumMath.cs b/src/Qubit.NET/Math/QuantumMath.cs
index 3805c95..fe48136 100644
--- a/src/Qubit.NET/Math/QuantumMath.cs
+++ b/src/Qubit.NET/Math/QuantumMath.cs
@@ -25,7 +25,9 @@ internal static class QuantumMath
     /// <returns>The updated quantum state vector.</returns>
     internal static Complex[] ApplySingleQubitGate(Complex[] state, Complex[,] gate, int targetQubit)
     {
-        if (BitOperations.PopCount((uint)state.Length) != 1)
+        // Power-of-two test. Written by hand rather than with BitOperations.PopCount,
+        // which does not exist on netstandard2.1 (Unity).
+        if (state.Length == 0 || (state.Length & (state.Length - 1)) != 0)
             throw new ArgumentException("State vector length must be a power of 2.");
         
         if (gate.GetLength(0) != 2 || gate.GetLength(1) != 2)
diff --git a/src/Qubit.NET/QuantumCircuit.cs b/src/Qubit.NET/QuantumCircuit.cs
index 92a6b45..a6c7dd4 100644
--- a/src/Qubit.NET/QuantumCircuit.cs
+++ b/src/Qubit.NET/QuantumCircuit.cs
@@ -963,7 +963,9 @@ public void Custom(Complex[,] matrix, params int[] qubits)
         
         Gates.Add(customGate);
         
-        ApplyGate(matrix, qubits.Reverse().ToArray());
+        // Enumerable.Reverse spelled out: in C# 14 `qubits.Reverse()` binds to
+        // MemoryExtensions.Reverse(Span<T>), which reverses in place and returns void.
+        ApplyGate(matrix, Enumerable.Reverse(qubits).ToArray());
     }
     
     /// <summary>
diff --git a/src/Qubit.NET/QuantumCircuitDrawer.cs b/src/Qubit.NET/QuantumCircuitDrawer.cs
index 262f5f0..d246edc 100644
--- a/src/Qubit.NET/QuantumCircuitDrawer.cs
+++ b/src/Qubit.NET/QuantumCircuitDrawer.cs
@@ -1,16 +1,67 @@
+using System.Text;
 using Qubit.NET.Gates;
 using Qubit.NET.Utilities;
 
 namespace Qubit.NET;
 
+/// <summary>
+/// Renders quantum circuits as ASCII diagrams.
+/// </summary>
 public static class QuantumCircuitDrawer
 {
     /// <summary>
-    /// Draws an ASCII representation of the quantum circuit in the console.
+    /// Draws an ASCII representation of the quantum circuit in the console,
+    /// with gates highlighted using console colors.
     /// This includes gate positions across qubits and visual connections between them.
     /// </summary>
     /// <param name="circuit">The quantum circuit to draw.</param>
+    /// <remarks>
+    /// Requires a console. In environments without one (Unity, ASP.NET, tests),
+    /// use <see cref="ToDiagram"/> and write the string wherever you need it.
+    /// </remarks>
     public static void Draw(this QuantumCircuit circuit)
+    {
+        ConsoleColor defaultColor = Console.ForegroundColor;
+
+        foreach (string line in circuit.ToDiagram().Split('\n'))
+        {
+            // Qubit labels ("q0 (0): ") are yellow, gate glyphs are blue, wires stay plain.
+            // Connector-only lines carry no label, so bodyStart is 0 there.
+            int labelEnd = line.IndexOf(": ", StringComparison.Ordinal);
+            int bodyStart = labelEnd >= 0 ? labelEnd + 2 : 0;
+
+            if (labelEnd >= 0)
+            {
+                Console.ForegroundColor = ConsoleColor.Yellow;
+                Console.Write(line.Substring(0, bodyStart));
+                Console.ForegroundColor = defaultColor;
+            }
+
+            foreach (char c in line.Substring(bodyStart))
+            {
+                bool glyph = c is '[' or ']' or '@' or '|' || (c != '─' && c != ' ' && c != '\r');
+
+                Console.ForegroundColor = glyph ? ConsoleColor.DarkBlue : defaultColor;
+                Console.Write(c);
+            }
+
+            Console.ForegroundColor = defaultColor;
+            Console.WriteLine();
+        }
+    }
+
+    /// <summary>
+    /// Builds an ASCII representation of the quantum circuit and returns it as a string.
+    /// This includes gate positions across qubits and visual connections between them.
+    /// </summary>
+    /// <param name="circuit">The quantum circuit to render.</param>
+    /// <returns>A multi-line string containing the circuit diagram.</returns>
+    /// <example>
+    /// <code>
+    /// Debug.Log(circuit.ToDiagram());   // Unity
+    /// </code>
+    /// </example>
+    public static string ToDiagram(this QuantumCircuit circuit)
     {
         IList<IList<(string, int)>> gatePositions = new List<IList<(string, int)>>();
         IList<IList<int>> barPositions = new List<IList<int>>();
@@ -19,10 +70,14 @@ public static void Draw(this QuantumCircuit circuit)
         InitializeStructures(gatePositions, barPositions, circuit.QubitCount);
 
         int lastGate = AssignGatePositions(gatePositions, barPositions, gateWidths, circuit);
-        
-        PrintGates(gatePositions, barPositions, gateWidths, circuit, lastGate);
+
+        StringBuilder sb = new StringBuilder();
+
+        PrintGates(sb, gatePositions, barPositions, gateWidths, circuit, lastGate);
+
+        return sb.ToString();
     }
-    
+
     /// <summary>
     /// Initializes the data structures used for tracking gate symbols and bar positions
     /// for each qubit line in the quantum circuit.
@@ -145,8 +200,9 @@ private static int AssignGatePositions(IList<IList<(string, int)>> gatePositions
     }
 
     /// <summary>
-    /// Renders a visual representation of the quantum circuit to the console.
+    /// Renders a visual representation of the quantum circuit into a string builder.
     /// </summary>
+    /// <param name="sb">The builder receiving the diagram.</param>
     /// <param name="gatePositions">A list of gate symbols and their positions for each qubit line.</param>
     /// <param name="barPositions">A list of vertical bar positions for multi-qubit gates.</param>
     /// <param name="widths">A list where each integer holds an additional width of a gate.</param>
@@ -154,86 +210,75 @@ private static int AssignGatePositions(IList<IList<(string, int)>> gatePositions
     /// <param name="lastGate">The index of the rightmost gate, used for alignment and padding.</param>
     /// <remarks>
     /// This method draws each qubit line with its gates and initial state, using ASCII characters.
-    /// Gates are color-highlighted and aligned based on position. Multi-qubit gates are connected with vertical bars.
+    /// Multi-qubit gates are connected with vertical bars. Coloring is applied separately by
+    /// <see cref="Draw"/> so that this renderer stays free of any console dependency.
     /// </remarks>
-    private static void PrintGates(IList<IList<(string, int)>> gatePositions, IList<IList<int>> barPositions,
-        IList<int> widths, QuantumCircuit circuit, int lastGate)
+    private static void PrintGates(StringBuilder sb, IList<IList<(string, int)>> gatePositions,
+        IList<IList<int>> barPositions, IList<int> widths, QuantumCircuit circuit, int lastGate)
     {
-        ConsoleColor defaultColor = Console.ForegroundColor;
-
         int counter = -1;
         for (int i = circuit.QubitCount-1; i >= 0; i--)
         {
             InitialState? initialState = circuit.Initializations.FirstOrDefault(init => init.QubitIndex == i);
             char initState = initialState == null ? '0' : Helpers.InitialStateToCharRepresentation(initialState.BasicState);
-            
-            Console.ForegroundColor = ConsoleColor.Yellow;
-            Console.Write($"q{i} ({initState}): ");
-            Console.ForegroundColor = defaultColor;
-            
+
+            sb.Append($"q{i} ({initState}): ");
+
             foreach (var gatePos in gatePositions[circuit.QubitCount-1-i])
             {
                 int additionalWidth = widths.Skip(counter+1).Take(gatePos.Item2 - counter - 1).Sum();
-                Console.Write(new string('\u2500', ((gatePos.Item2-counter-1)*5)+additionalWidth));
+                sb.Append('\u2500', ((gatePos.Item2-counter-1)*5)+additionalWidth);
                 counter = gatePos.Item2;
-                
+
                 if (gatePos.Item1 == "|")
                 {
-                    Console.Write(new string('\u2500', 2));
-                    Console.ForegroundColor = ConsoleColor.DarkBlue;
-                    Console.Write("|");
-                    Console.ForegroundColor = defaultColor;
-                    Console.Write(new string('\u2500', 2+widths[counter]));
+                    sb.Append('\u2500', 2);
+                    sb.Append('|');
+                    sb.Append('\u2500', 2+widths[counter]);
                     continue;
-                } 
+                }
                 if (gatePos.Item1 == "@")
                 {
-                    Console.Write(new string('\u2500', 2));
-                    Console.ForegroundColor = ConsoleColor.DarkBlue;
-                    Console.Write("@");
-                    Console.ForegroundColor = defaultColor;
-                    Console.Write(new string('\u2500', (2+widths[counter])));
+                    sb.Append('\u2500', 2);
+                    sb.Append('@');
+                    sb.Append('\u2500', 2+widths[counter]);
                     continue;
                 }
-                
-                Console.Write("\u2500");
-                Console.ForegroundColor = ConsoleColor.DarkBlue;
-                Console.Write($"[{gatePos.Item1}]");
-                Console.ForegroundColor = defaultColor;
-                Console.Write(new string('\u2500', (2+widths[counter]-gatePos.Item1.Length)));
+
+                sb.Append('\u2500');
+                sb.Append('[').Append(gatePos.Item1).Append(']');
+                sb.Append('\u2500', 2+widths[counter]-gatePos.Item1.Length);
             }
 
             if (!gatePositions[circuit.QubitCount - 1 - i].Any())
             {
                 int additionalWidth = widths.Sum();
-                Console.Write(new string('\u2500', ((lastGate+1)*5)+additionalWidth));
-            } 
+                sb.Append('\u2500', ((lastGate+1)*5)+additionalWidth);
+            }
             else if (counter < lastGate)
             {
                 int additionalWidth = widths.Skip(counter+1).Take(lastGate - counter).Sum();
-                Console.Write(new string('\u2500', ((lastGate-counter)*5)+additionalWidth));
+                sb.Append('\u2500', ((lastGate-counter)*5)+additionalWidth);
             }
 
-            Console.WriteLine();
+            sb.AppendLine();
             counter = -1;
-            
-            Console.Write(new String(' ', 8));
+
+            sb.Append(' ', 8);
             foreach (var barPos in barPositions[circuit.QubitCount-1-i])
             {
                 if(i == 0) continue;
 
                 int additionalWidth = widths.Skip(counter+1).Take(barPos - counter - 1).Sum();
-                Console.Write(new string(' ', ((barPos-counter-1)*5)+additionalWidth));
+                sb.Append(' ', ((barPos-counter-1)*5)+additionalWidth);
                 counter = barPos;
-                
-                Console.Write(new string(' ', 2));
-                Console.ForegroundColor = ConsoleColor.DarkBlue;
-                Console.Write("|");
-                Console.ForegroundColor = defaultColor;
-                Console.Write(new string(' ', (2+widths[barPos])));
+
+                sb.Append(' ', 2);
+                sb.Append('|');
+                sb.Append(' ', 2+widths[barPos]);
             }
-            
-            Console.WriteLine();
+
+            sb.AppendLine();
             counter = -1;
         }
     }
diff --git a/src/Qubit.NET/Qubit.NET.csproj b/src/Qubit.NET/Qubit.NET.csproj
index fa71b7a..68ed256 100644
--- a/src/Qubit.NET/Qubit.NET.csproj
+++ b/src/Qubit.NET/Qubit.NET.csproj
@@ -1,9 +1,45 @@
-<Project Sdk="Microsoft.NET.Sdk">
+<Project Sdk="Microsoft.NET.Sdk">
 
   <PropertyGroup>
-    <TargetFramework>net8.0</TargetFramework>
+    <!-- netstandard2.1 is what makes the library usable from Unity. -->
+    <TargetFrameworks>netstandard2.1;net8.0;net10.0</TargetFrameworks>
+    <LangVersion>latest</LangVersion>
     <ImplicitUsings>enable</ImplicitUsings>
     <Nullable>enable</Nullable>
+    <GenerateDocumentationFile>true</GenerateDocumentationFile>
   </PropertyGroup>
 
+  <PropertyGroup>
+    <PackageId>Qubit.NET</PackageId>
+    <Version>1.0.0</Version>
+    <Authors>Tymoteusz Marzec</Authors>
+    <Description>A lightweight quantum computing simulation library for .NET. Build quantum circuits, apply gates, measure results, and draw ASCII circuit diagrams — no dependencies. Works in .NET and Unity.</Description>
+    <PackageTags>quantum;quantum-computing;qubit;simulator;quantum-circuit;unity;csharp;dotnet</PackageTags>
+    <PackageLicenseExpression>MIT</PackageLicenseExpression>
+    <PackageProjectUrl>https://github.com/InfoTCube/Qubit.NET</PackageProjectUrl>
+    <RepositoryUrl>https://github.com/InfoTCube/Qubit.NET</RepositoryUrl>
+    <RepositoryType>git</RepositoryType>
+    <PackageIcon>qubitnet.png</PackageIcon>
+    <PackageReadmeFile>README.md</PackageReadmeFile>
+    <PackageReleaseNotes>See CHANGELOG.md</PackageReleaseNotes>
+  </PropertyGroup>
+
+  <PropertyGroup>
+    <!-- Source Link: lets consumers step into the library's source while debugging. -->
+    <PublishRepositoryUrl>true</PublishRepositoryUrl>
+    <EmbedUntrackedSources>true</EmbedUntrackedSources>
+    <IncludeSymbols>true</IncludeSymbols>
+    <SymbolPackageFormat>snupkg</SymbolPackageFormat>
+    <ContinuousIntegrationBuild Condition="'$(GITHUB_ACTIONS)' == 'true'">true</ContinuousIntegrationBuild>
+  </PropertyGroup>
+
+  <ItemGroup>
+    <None Include="../../img/qubitnet.png" Pack="true" PackagePath="/" Visible="false" />
+    <None Include="../../README.md" Pack="true" PackagePath="/" Visible="false" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <PackageReference Include="Microsoft.SourceLink.GitHub" Version="8.0.0" PrivateAssets="All" />
+  </ItemGroup>
+
 </Project>
diff --git a/src/Qubit.NET/Simulation/Simulator.cs b/src/Qubit.NET/Simulation/Simulator.cs
index e3c9c18..4ff654a 100644
--- a/src/Qubit.NET/Simulation/Simulator.cs
+++ b/src/Qubit.NET/Simulation/Simulator.cs
@@ -150,7 +150,7 @@ or GateType.T or GateType.Tdag or GateType.Rx or GateType.Ry or GateType.Rz or G
                 break;
             case GateType.Custom:
                 stateVector = QuantumMath.ApplyMultiQubitGate(stateVector, currentGate.Matrix,
-                    currentGate.TargetQubits.Reverse().ToArray());
+                    Enumerable.Reverse(currentGate.TargetQubits).ToArray());
                 break;
         }
 
diff --git a/src/Qubit.NET/Utilities/State.cs b/src/Qubit.NET/Utilities/State.cs
index be68dbf..d157447 100644
--- a/src/Qubit.NET/Utilities/State.cs
+++ b/src/Qubit.NET/Utilities/State.cs
@@ -1,10 +1,22 @@
 namespace Qubit.NET.Utilities;
 
+/// <summary>
+/// A predefined single-qubit basis state that a qubit can be initialized to.
+/// </summary>
 public enum State
 {
+    /// <summary>The |0⟩ state.</summary>
     Zero,
+
+    /// <summary>The |1⟩ state.</summary>
     One,
+
+    /// <summary>The |+⟩ state, (|0⟩ + |1⟩)/√2.</summary>
     Plus,
+
+    /// <summary>The |−⟩ state, (|0⟩ − |1⟩)/√2.</summary>
     Minus,
+
+    /// <summary>An arbitrary state given by explicit α and β amplitudes.</summary>
     Custom
-}
\ No newline at end of file
+}

From 55045a37bb9a8e03da3e8f1aafa2753196c23361 Mon Sep 17 00:00:00 2001
From: InfoTCube <infocube.yout@gmail.com>
Date: Mon, 17 Aug 2026 20:14:47 +0200
Subject: [PATCH 3/9] test: add an xUnit suite covering gates, circuits,
 measurement and the drawer

84 tests, previously zero. A simulation library that nobody can verify is a
library nobody will depend on.

- Gate algebra: every matrix is unitary, SX/SY/SZ square to their parent gate,
  H*H = I, S*S = Z, T^4 = Z, each dagger inverts its gate, Rx(2pi) = -I, and
  U3 reproduces H at the standard angles.
- Circuits: all four Bell states and GHZ match their textbook amplitudes;
  Toffoli, SWAP and the controlled gates behave on basis states.
- Measurement: with a fixed random source, measuring a Bell state yields only
  "00" or "11" across 50 seeds, collapse leaves the state normalized, and
  measuring one half of a Bell pair determines the other.
- Drawer: golden strings for alignment, plus a check that ToDiagram needs no
  console and that Draw prints exactly what ToDiagram returns.
- One regression test per bug fixed in the previous two commits.

The tests live in namespace QubitNet.Tests rather than Qubit.NET.Tests, because
the latter puts Qubit.NET.Math in scope where it shadows System.Math.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Mhwc4XeGtzM2UQx9sNLb3r
---
 CHANGELOG.md                                 |   2 +
 Qubit.NET.sln                                |  38 ++++
 tests/Qubit.NET.Tests/CircuitTests.cs        | 228 +++++++++++++++++++
 tests/Qubit.NET.Tests/DrawerTests.cs         | 140 ++++++++++++
 tests/Qubit.NET.Tests/GateAlgebraTests.cs    | 130 +++++++++++
 tests/Qubit.NET.Tests/Qubit.NET.Tests.csproj |  27 +++
 tests/Qubit.NET.Tests/RegressionTests.cs     | 152 +++++++++++++
 tests/Qubit.NET.Tests/TestHelpers.cs         | 102 +++++++++
 8 files changed, 819 insertions(+)
 create mode 100644 tests/Qubit.NET.Tests/CircuitTests.cs
 create mode 100644 tests/Qubit.NET.Tests/DrawerTests.cs
 create mode 100644 tests/Qubit.NET.Tests/GateAlgebraTests.cs
 create mode 100644 tests/Qubit.NET.Tests/Qubit.NET.Tests.csproj
 create mode 100644 tests/Qubit.NET.Tests/RegressionTests.cs
 create mode 100644 tests/Qubit.NET.Tests/TestHelpers.cs

diff --git a/CHANGELOG.md b/CHANGELOG.md
index fc5363e..4b9e462 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -18,6 +18,8 @@ and usable from Unity.
 - `QuantumGates.Format(matrix)` — the string-returning counterpart of `Print`.
 - `QuantumCircuit.MaxQubitCount` constant.
 - Continuous integration and a tag-driven NuGet release workflow.
+- An xUnit test suite covering gate algebra, the textbook Bell and GHZ states, measurement
+  statistics and collapse, circuit rendering, and a regression test for every bug below.
 - XML documentation for every public gate matrix and for the `State` enum.
 
 ### Changed
diff --git a/Qubit.NET.sln b/Qubit.NET.sln
index abffc65..48b71e3 100644
--- a/Qubit.NET.sln
+++ b/Qubit.NET.sln
@@ -1,3 +1,4 @@
+
 Microsoft Visual Studio Solution File, Format Version 12.00
 # Visual Studio Version 17
 VisualStudioVersion = 17.5.2.0
@@ -8,26 +9,63 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "src", "src", "{827E0CD3-B72
 EndProject
 Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Qubit.NET", "src\Qubit.NET\Qubit.NET.csproj", "{F1519400-F1F2-D001-2DB1-E759C4458259}"
 EndProject
+Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "tests", "tests", "{0AB3BF05-4346-4AA6-1389-037BE0695223}"
+EndProject
+Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Qubit.NET.Tests", "tests\Qubit.NET.Tests\Qubit.NET.Tests.csproj", "{A1259BCF-C8A1-4D0C-922F-387222B34D44}"
+EndProject
 Global
 	GlobalSection(SolutionConfigurationPlatforms) = preSolution
 		Debug|Any CPU = Debug|Any CPU
+		Debug|x64 = Debug|x64
+		Debug|x86 = Debug|x86
 		Release|Any CPU = Release|Any CPU
+		Release|x64 = Release|x64
+		Release|x86 = Release|x86
 	EndGlobalSection
 	GlobalSection(ProjectConfigurationPlatforms) = postSolution
 		{8028A06F-E855-BC4E-B4C3-3FCEB15E0C0D}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
 		{8028A06F-E855-BC4E-B4C3-3FCEB15E0C0D}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{8028A06F-E855-BC4E-B4C3-3FCEB15E0C0D}.Debug|x64.ActiveCfg = Debug|Any CPU
+		{8028A06F-E855-BC4E-B4C3-3FCEB15E0C0D}.Debug|x64.Build.0 = Debug|Any CPU
+		{8028A06F-E855-BC4E-B4C3-3FCEB15E0C0D}.Debug|x86.ActiveCfg = Debug|Any CPU
+		{8028A06F-E855-BC4E-B4C3-3FCEB15E0C0D}.Debug|x86.Build.0 = Debug|Any CPU
 		{8028A06F-E855-BC4E-B4C3-3FCEB15E0C0D}.Release|Any CPU.ActiveCfg = Release|Any CPU
 		{8028A06F-E855-BC4E-B4C3-3FCEB15E0C0D}.Release|Any CPU.Build.0 = Release|Any CPU
+		{8028A06F-E855-BC4E-B4C3-3FCEB15E0C0D}.Release|x64.ActiveCfg = Release|Any CPU
+		{8028A06F-E855-BC4E-B4C3-3FCEB15E0C0D}.Release|x64.Build.0 = Release|Any CPU
+		{8028A06F-E855-BC4E-B4C3-3FCEB15E0C0D}.Release|x86.ActiveCfg = Release|Any CPU
+		{8028A06F-E855-BC4E-B4C3-3FCEB15E0C0D}.Release|x86.Build.0 = Release|Any CPU
 		{F1519400-F1F2-D001-2DB1-E759C4458259}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
 		{F1519400-F1F2-D001-2DB1-E759C4458259}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{F1519400-F1F2-D001-2DB1-E759C4458259}.Debug|x64.ActiveCfg = Debug|Any CPU
+		{F1519400-F1F2-D001-2DB1-E759C4458259}.Debug|x64.Build.0 = Debug|Any CPU
+		{F1519400-F1F2-D001-2DB1-E759C4458259}.Debug|x86.ActiveCfg = Debug|Any CPU
+		{F1519400-F1F2-D001-2DB1-E759C4458259}.Debug|x86.Build.0 = Debug|Any CPU
 		{F1519400-F1F2-D001-2DB1-E759C4458259}.Release|Any CPU.ActiveCfg = Release|Any CPU
 		{F1519400-F1F2-D001-2DB1-E759C4458259}.Release|Any CPU.Build.0 = Release|Any CPU
+		{F1519400-F1F2-D001-2DB1-E759C4458259}.Release|x64.ActiveCfg = Release|Any CPU
+		{F1519400-F1F2-D001-2DB1-E759C4458259}.Release|x64.Build.0 = Release|Any CPU
+		{F1519400-F1F2-D001-2DB1-E759C4458259}.Release|x86.ActiveCfg = Release|Any CPU
+		{F1519400-F1F2-D001-2DB1-E759C4458259}.Release|x86.Build.0 = Release|Any CPU
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Debug|x64.ActiveCfg = Debug|Any CPU
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Debug|x64.Build.0 = Debug|Any CPU
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Debug|x86.ActiveCfg = Debug|Any CPU
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Debug|x86.Build.0 = Debug|Any CPU
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Release|Any CPU.Build.0 = Release|Any CPU
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Release|x64.ActiveCfg = Release|Any CPU
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Release|x64.Build.0 = Release|Any CPU
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Release|x86.ActiveCfg = Release|Any CPU
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Release|x86.Build.0 = Release|Any CPU
 	EndGlobalSection
 	GlobalSection(SolutionProperties) = preSolution
 		HideSolutionNode = FALSE
 	EndGlobalSection
 	GlobalSection(NestedProjects) = preSolution
 		{F1519400-F1F2-D001-2DB1-E759C4458259} = {827E0CD3-B72D-47B6-A68D-7590B98EB39B}
+		{A1259BCF-C8A1-4D0C-922F-387222B34D44} = {0AB3BF05-4346-4AA6-1389-037BE0695223}
 	EndGlobalSection
 	GlobalSection(ExtensibilityGlobals) = postSolution
 		SolutionGuid = {51965E37-13AF-4AC3-A5D2-1426ACB4DE21}
diff --git a/tests/Qubit.NET.Tests/CircuitTests.cs b/tests/Qubit.NET.Tests/CircuitTests.cs
new file mode 100644
index 0000000..7e16848
--- /dev/null
+++ b/tests/Qubit.NET.Tests/CircuitTests.cs
@@ -0,0 +1,228 @@
+using System.Numerics;
+using Qubit.NET.Circuits;
+using Qubit.NET.Utilities;
+using static QubitNet.Tests.Amplitudes;
+using Qubit.NET;
+
+namespace QubitNet.Tests;
+
+/// <summary>
+/// End-to-end checks that circuits produce the textbook states, and that measurement
+/// collapses them correctly.
+/// </summary>
+public class CircuitTests
+{
+    // State vector index i encodes qubit q in bit q, so for two qubits the order is
+    // |q1 q0>: index 1 is |01> (qubit 0 set), index 2 is |10> (qubit 1 set).
+
+    [Fact]
+    public void PhiPlus_is_00_plus_11()
+    {
+        AssertState(BellStates.PhiPlus(), InvSqrt2, 0, 0, InvSqrt2);
+        AssertNormalized(BellStates.PhiPlus());
+    }
+
+    [Fact]
+    public void PhiMinus_is_00_minus_11() =>
+        AssertState(BellStates.PhiMinus(), InvSqrt2, 0, 0, -InvSqrt2);
+
+    [Fact]
+    public void PsiPlus_is_01_plus_10() =>
+        AssertState(BellStates.PsiPlus(), 0, InvSqrt2, InvSqrt2, 0);
+
+    [Fact]
+    public void PsiMinus_is_01_minus_10() =>
+        AssertState(BellStates.PsiMinus(), 0, InvSqrt2, -InvSqrt2, 0);
+
+    [Fact]
+    public void GHZ_is_000_plus_111() =>
+        AssertState(BellStates.GHZ(), InvSqrt2, 0, 0, 0, 0, 0, 0, InvSqrt2);
+
+    [Fact]
+    public void Hadamard_creates_an_equal_superposition()
+    {
+        QuantumCircuit qc = new(1);
+        qc.H(0);
+
+        AssertState(qc, InvSqrt2, InvSqrt2);
+    }
+
+    [Fact]
+    public void X_flips_the_qubit()
+    {
+        QuantumCircuit qc = new(1);
+        qc.X(0);
+
+        AssertState(qc, 0, 1);
+    }
+
+    [Fact]
+    public void Applying_a_gate_twice_returns_to_the_start()
+    {
+        QuantumCircuit qc = new(1);
+        qc.H(0);
+        qc.H(0);
+
+        AssertState(qc, 1, 0);
+    }
+
+    [Fact]
+    public void Toffoli_flips_the_target_only_when_both_controls_are_set()
+    {
+        QuantumCircuit both = new(3);
+        both.X(0);
+        both.X(1);
+        both.Toffoli(0, 1, 2);
+
+        // |111> is index 7.
+        AssertState(both, 0, 0, 0, 0, 0, 0, 0, 1);
+
+        QuantumCircuit one = new(3);
+        one.X(0);
+        one.Toffoli(0, 1, 2);
+
+        // Target must stay |0>, so the state is still |001> = index 1.
+        AssertState(one, 0, 1, 0, 0, 0, 0, 0, 0);
+    }
+
+    [Fact]
+    public void SWAP_exchanges_two_qubits()
+    {
+        QuantumCircuit qc = new(2);
+        qc.X(0);
+        qc.SWAP(0, 1);
+
+        // |01> becomes |10>, i.e. index 1 becomes index 2.
+        AssertState(qc, 0, 0, 1, 0);
+    }
+
+    [Fact]
+    public void GetProbabilities_reports_a_fair_coin_for_a_Bell_state()
+    {
+        var probabilities = BellStates.PhiPlus().GetProbabilities();
+
+        Assert.Equal(2, probabilities.Count);
+        Assert.Equal(0.5, probabilities["00"], Tolerance);
+        Assert.Equal(0.5, probabilities["11"], Tolerance);
+    }
+
+    [Theory]
+    [InlineData(0.1, "00")]
+    [InlineData(0.9, "11")]
+    public void Measuring_a_Bell_state_never_yields_a_mixed_outcome(double roll, string expected)
+    {
+        QuantumCircuit qc = BellStates.PhiPlus();
+        qc.RandomSource = new FixedRandomSource(roll);
+
+        Assert.Equal(expected, qc.Measure());
+    }
+
+    [Fact]
+    public void Bell_state_measurements_are_perfectly_correlated()
+    {
+        // The whole point of entanglement: "01" and "10" must never appear.
+        for (int seed = 0; seed < 50; seed++)
+        {
+            QuantumCircuit qc = BellStates.PhiPlus();
+            qc.RandomSource = new SeededRandomSource(seed);
+
+            string result = qc.Measure();
+
+            Assert.True(result is "00" or "11", $"seed {seed} produced {result}");
+        }
+    }
+
+    [Fact]
+    public void Measurement_collapses_the_state()
+    {
+        QuantumCircuit qc = BellStates.PhiPlus();
+        qc.RandomSource = new FixedRandomSource(0.1);
+
+        qc.Measure();
+
+        AssertState(qc, 1, 0, 0, 0);
+    }
+
+    [Fact]
+    public void Measuring_one_qubit_of_a_Bell_pair_collapses_the_other()
+    {
+        QuantumCircuit qc = BellStates.PhiPlus();
+        qc.RandomSource = new FixedRandomSource(0.9);
+
+        string result = qc.Measure(0);
+
+        Assert.Equal("1", result);
+        AssertNormalized(qc);
+
+        // Qubit 0 came out |1>, so qubit 1 must now be |1> as well: only |11> survives.
+        var probabilities = qc.GetProbabilities();
+        Assert.Single(probabilities);
+        Assert.Equal(1.0, probabilities["11"], Tolerance);
+    }
+
+    [Fact]
+    public void Partial_measurement_returns_bits_in_argument_order()
+    {
+        QuantumCircuit qc = new(3);
+        qc.X(2);   // |100>
+
+        Assert.Equal("01", qc.Measure(0, 2));
+        Assert.Equal("10", new QuantumCircuit(3).Also(c => c.X(2)).Measure(2, 0));
+    }
+
+    [Theory]
+    [InlineData(State.Zero, 1, 0)]
+    [InlineData(State.One, 0, 1)]
+    public void Initialize_sets_a_basis_state(State state, double alpha, double beta)
+    {
+        QuantumCircuit qc = new(1);
+        qc.Initialize(0, state);
+
+        AssertState(qc, alpha, beta);
+    }
+
+    [Fact]
+    public void Initialize_sets_the_plus_state()
+    {
+        QuantumCircuit qc = new(1);
+        qc.Initialize(0, State.Plus);
+
+        AssertState(qc, InvSqrt2, InvSqrt2);
+    }
+
+    [Fact]
+    public void Initialize_rejects_an_already_modified_qubit()
+    {
+        QuantumCircuit qc = new(1);
+        qc.H(0);
+
+        Assert.Throws<InvalidOperationException>(() => qc.Initialize(0, State.One));
+    }
+
+    [Theory]
+    [InlineData(0)]
+    [InlineData(-1)]
+    [InlineData(27)]
+    public void Invalid_qubit_counts_are_rejected(int count) =>
+        Assert.Throws<ArgumentOutOfRangeException>(() => new QuantumCircuit(count));
+
+    [Theory]
+    [InlineData(-1)]
+    [InlineData(2)]
+    public void Out_of_range_qubit_indices_are_rejected(int qubit) =>
+        Assert.Throws<QubitIndexOutOfRangeException>(() => new QuantumCircuit(2).H(qubit));
+
+    [Fact]
+    public void A_gate_cannot_target_the_same_qubit_twice() =>
+        Assert.Throws<ArgumentException>(() => new QuantumCircuit(2).CNOT(0, 0));
+}
+
+internal static class CircuitExtensions
+{
+    /// <summary>Applies an action to a circuit and returns it, for terser test setup.</summary>
+    public static QuantumCircuit Also(this QuantumCircuit qc, Action<QuantumCircuit> action)
+    {
+        action(qc);
+        return qc;
+    }
+}
diff --git a/tests/Qubit.NET.Tests/DrawerTests.cs b/tests/Qubit.NET.Tests/DrawerTests.cs
new file mode 100644
index 0000000..f851b7b
--- /dev/null
+++ b/tests/Qubit.NET.Tests/DrawerTests.cs
@@ -0,0 +1,140 @@
+using Qubit.NET.Circuits;
+using Qubit.NET.Gates;
+using Qubit.NET.Utilities;
+using Qubit.NET;
+
+namespace QubitNet.Tests;
+
+/// <summary>
+/// Golden-string tests for the ASCII circuit renderer. These exist mostly to catch
+/// accidental changes to alignment, which is easy to break and hard to notice.
+/// </summary>
+public class DrawerTests
+{
+    private static string[] Lines(QuantumCircuit qc) =>
+        qc.ToDiagram().TrimEnd().Split('\n').Select(l => l.TrimEnd()).ToArray();
+
+    [Fact]
+    public void Bell_circuit_renders_as_expected()
+    {
+        string[] lines = Lines(BellStates.PhiPlus());
+
+        Assert.Equal("q1 (0): ──────[+]─", lines[0]);
+        Assert.Equal("               |", lines[1]);
+        Assert.Equal("q0 (0): ─[H]───@──", lines[2]);
+    }
+
+    [Fact]
+    public void Empty_circuit_still_draws_its_wires()
+    {
+        string[] lines = Lines(new QuantumCircuit(2));
+
+        Assert.Equal(2, lines.Count(l => l.StartsWith('q')));
+        Assert.All(lines.Where(l => l.StartsWith('q')), l => Assert.Contains('─', l));
+    }
+
+    [Fact]
+    public void Initialized_qubits_show_their_starting_state()
+    {
+        QuantumCircuit qc = new(3);
+        qc.Initialize(0, State.Plus);
+        qc.Initialize(1, State.One);
+
+        string diagram = qc.ToDiagram();
+
+        Assert.Contains("q0 (+)", diagram);
+        Assert.Contains("q1 (1)", diagram);
+        Assert.Contains("q2 (0)", diagram);
+    }
+
+    [Fact]
+    public void Every_gate_type_has_a_symbol()
+    {
+        QuantumCircuit qc = new(3);
+        qc.I(0); qc.H(0); qc.X(0); qc.Y(0); qc.Z(0);
+        qc.S(0); qc.Sdag(0); qc.T(0); qc.Tdag(0);
+        qc.Rx(0, 0.1); qc.Ry(0, 0.1); qc.Rz(0, 0.1);
+        qc.SX(0); qc.SY(0); qc.SZ(0); qc.U3(0, 0.1, 0.2, 0.3);
+        qc.CNOT(0, 1); qc.CY(0, 1); qc.CZ(0, 1); qc.CH(0, 1);
+        qc.SWAP(0, 1); qc.Toffoli(0, 1, 2); qc.Fredkin(0, 1, 2);
+        qc.Custom(QuantumGates.H, 0);
+        qc.Measure();
+
+        string diagram = qc.ToDiagram();
+
+        // A blank symbol would mean a gate type fell through the symbol table.
+        Assert.DoesNotContain("[ ]", diagram);
+        Assert.Contains("[H]", diagram);
+        Assert.Contains("[S†]", diagram);
+        Assert.Contains("[U3]", diagram);
+        Assert.Contains("[M]", diagram);
+        Assert.Contains("[C]", diagram);
+    }
+
+    [Fact]
+    public void Fredkin_draws_one_control_and_two_swap_targets()
+    {
+        QuantumCircuit qc = new(3);
+        qc.Fredkin(0, 1, 2);
+
+        string diagram = qc.ToDiagram();
+
+        Assert.Equal(1, diagram.Count(c => c == '@'));
+        Assert.Equal(2, CountOccurrences(diagram, "[X]"));
+    }
+
+    [Fact]
+    public void ToDiagram_needs_no_console()
+    {
+        // Draw() writes to the console; ToDiagram() must not, so that Unity and test hosts
+        // without a console can still render circuits.
+        TextWriter original = Console.Out;
+        Console.SetOut(TextWriter.Null);
+
+        try
+        {
+            Assert.NotEmpty(BellStates.GHZ().ToDiagram());
+        }
+        finally
+        {
+            Console.SetOut(original);
+        }
+    }
+
+    [Fact]
+    public void Draw_writes_the_same_text_that_ToDiagram_returns()
+    {
+        QuantumCircuit qc = BellStates.GHZ();
+
+        StringWriter captured = new();
+        TextWriter original = Console.Out;
+        Console.SetOut(captured);
+
+        try
+        {
+            qc.Draw();
+        }
+        finally
+        {
+            Console.SetOut(original);
+        }
+
+        Assert.Equal(Normalize(qc.ToDiagram()), Normalize(captured.ToString()));
+    }
+
+    private static string Normalize(string s) =>
+        string.Join("\n", s.Replace("\r", string.Empty).Split('\n').Select(l => l.TrimEnd())).TrimEnd();
+
+    private static int CountOccurrences(string haystack, string needle)
+    {
+        int count = 0, index = 0;
+
+        while ((index = haystack.IndexOf(needle, index, StringComparison.Ordinal)) >= 0)
+        {
+            count++;
+            index += needle.Length;
+        }
+
+        return count;
+    }
+}
diff --git a/tests/Qubit.NET.Tests/GateAlgebraTests.cs b/tests/Qubit.NET.Tests/GateAlgebraTests.cs
new file mode 100644
index 0000000..13bca51
--- /dev/null
+++ b/tests/Qubit.NET.Tests/GateAlgebraTests.cs
@@ -0,0 +1,130 @@
+using System.Numerics;
+using Qubit.NET.Gates;
+using static QubitNet.Tests.Amplitudes;
+using Qubit.NET;
+
+namespace QubitNet.Tests;
+
+/// <summary>
+/// Checks the gate matrices themselves: that they are unitary and that they satisfy the
+/// algebraic identities that define them.
+/// </summary>
+public class GateAlgebraTests
+{
+    public static TheoryData<string, Complex[,]> AllGates() => new()
+    {
+        { nameof(QuantumGates.I), QuantumGates.I },
+        { nameof(QuantumGates.H), QuantumGates.H },
+        { nameof(QuantumGates.X), QuantumGates.X },
+        { nameof(QuantumGates.Y), QuantumGates.Y },
+        { nameof(QuantumGates.Z), QuantumGates.Z },
+        { nameof(QuantumGates.S), QuantumGates.S },
+        { nameof(QuantumGates.Sdag), QuantumGates.Sdag },
+        { nameof(QuantumGates.T), QuantumGates.T },
+        { nameof(QuantumGates.Tdag), QuantumGates.Tdag },
+        { nameof(QuantumGates.SX), QuantumGates.SX },
+        { nameof(QuantumGates.SY), QuantumGates.SY },
+        { nameof(QuantumGates.SZ), QuantumGates.SZ },
+        { nameof(QuantumGates.CNOT), QuantumGates.CNOT },
+        { nameof(QuantumGates.CY), QuantumGates.CY },
+        { nameof(QuantumGates.CZ), QuantumGates.CZ },
+        { nameof(QuantumGates.CH), QuantumGates.CH },
+        { nameof(QuantumGates.SWAP), QuantumGates.SWAP },
+        { nameof(QuantumGates.Toffoli), QuantumGates.Toffoli },
+        { nameof(QuantumGates.Fredkin), QuantumGates.Fredkin },
+        { "Rx(0.7)", QuantumGates.Rx(0.7) },
+        { "Ry(0.7)", QuantumGates.Ry(0.7) },
+        { "Rz(0.7)", QuantumGates.Rz(0.7) },
+        { "U3(0.3,0.5,0.7)", QuantumGates.U3(0.3, 0.5, 0.7) },
+        { "CRx(0.7)", QuantumGates.CRx(0.7) },
+        { "CRy(0.7)", QuantumGates.CRy(0.7) },
+        { "CRz(0.7)", QuantumGates.CRz(0.7) },
+        { "CU3(0.3,0.5,0.7)", QuantumGates.CU3(0.3, 0.5, 0.7) },
+    };
+
+    [Theory]
+    [MemberData(nameof(AllGates))]
+    public void Every_gate_is_unitary(string name, Complex[,] gate)
+    {
+        // U†U must be the identity. Exercised through Custom(), which validates unitarity
+        // and only accepts matrices of up to 4 qubits.
+        int qubits = (int)Math.Log2(gate.GetLength(0));
+
+        QuantumCircuit qc = new(qubits);
+
+        Exception? thrown = Record.Exception(() => qc.Custom(gate, Enumerable.Range(0, qubits).ToArray()));
+
+        Assert.True(thrown is null, $"{name} was rejected as non-unitary: {thrown?.Message}");
+    }
+
+    [Fact]
+    public void Non_unitary_matrix_is_rejected()
+    {
+        QuantumCircuit qc = new(1);
+
+        Assert.Throws<ArgumentException>(() => qc.Custom(new Complex[,] { { 1, 1 }, { 1, 1 } }, 0));
+    }
+
+    [Theory]
+    [InlineData(nameof(QuantumGates.SX), nameof(QuantumGates.X))]
+    [InlineData(nameof(QuantumGates.SY), nameof(QuantumGates.Y))]
+    [InlineData(nameof(QuantumGates.SZ), nameof(QuantumGates.Z))]
+    public void Square_root_gate_squared_equals_its_parent(string root, string parent)
+    {
+        Complex[,] rootMatrix = Matrix(root);
+        Complex[,] parentMatrix = Matrix(parent);
+
+        AssertMatrixEqual(parentMatrix, Multiply(rootMatrix, rootMatrix));
+    }
+
+    [Fact]
+    public void H_squared_is_identity() =>
+        AssertMatrixEqual(QuantumGates.I, Multiply(QuantumGates.H, QuantumGates.H));
+
+    [Fact]
+    public void S_squared_is_Z() =>
+        AssertMatrixEqual(QuantumGates.Z, Multiply(QuantumGates.S, QuantumGates.S));
+
+    [Fact]
+    public void T_to_the_fourth_is_Z()
+    {
+        Complex[,] tSquared = Multiply(QuantumGates.T, QuantumGates.T);
+
+        AssertMatrixEqual(QuantumGates.Z, Multiply(tSquared, tSquared));
+    }
+
+    [Theory]
+    [InlineData(nameof(QuantumGates.S), nameof(QuantumGates.Sdag))]
+    [InlineData(nameof(QuantumGates.T), nameof(QuantumGates.Tdag))]
+    public void Gate_times_its_dagger_is_identity(string gate, string dagger) =>
+        AssertMatrixEqual(QuantumGates.I, Multiply(Matrix(gate), Matrix(dagger)));
+
+    [Fact]
+    public void Rotation_by_two_pi_about_x_is_minus_identity()
+    {
+        // Rx(2*pi) = -I: a 2*pi rotation returns a spin-1/2 system to itself up to a global phase.
+        Complex[,] expected = { { -1, 0 }, { 0, -1 } };
+
+        AssertMatrixEqual(expected, QuantumGates.Rx(2 * Math.PI));
+    }
+
+    [Fact]
+    public void U3_reproduces_H_at_the_standard_angles() =>
+        AssertMatrixEqual(QuantumGates.H, QuantumGates.U3(Math.PI / 2, 0, Math.PI));
+
+    private static Complex[,] Matrix(string name) => name switch
+    {
+        nameof(QuantumGates.I) => QuantumGates.I,
+        nameof(QuantumGates.X) => QuantumGates.X,
+        nameof(QuantumGates.Y) => QuantumGates.Y,
+        nameof(QuantumGates.Z) => QuantumGates.Z,
+        nameof(QuantumGates.S) => QuantumGates.S,
+        nameof(QuantumGates.Sdag) => QuantumGates.Sdag,
+        nameof(QuantumGates.T) => QuantumGates.T,
+        nameof(QuantumGates.Tdag) => QuantumGates.Tdag,
+        nameof(QuantumGates.SX) => QuantumGates.SX,
+        nameof(QuantumGates.SY) => QuantumGates.SY,
+        nameof(QuantumGates.SZ) => QuantumGates.SZ,
+        _ => throw new ArgumentOutOfRangeException(nameof(name), name, "unknown gate")
+    };
+}
diff --git a/tests/Qubit.NET.Tests/Qubit.NET.Tests.csproj b/tests/Qubit.NET.Tests/Qubit.NET.Tests.csproj
new file mode 100644
index 0000000..67af570
--- /dev/null
+++ b/tests/Qubit.NET.Tests/Qubit.NET.Tests.csproj
@@ -0,0 +1,27 @@
+<Project Sdk="Microsoft.NET.Sdk">
+
+  <PropertyGroup>
+    <TargetFramework>net10.0</TargetFramework>
+    <ImplicitUsings>enable</ImplicitUsings>
+    <Nullable>enable</Nullable>
+    <IsPackable>false</IsPackable>
+  </PropertyGroup>
+
+  <ItemGroup>
+    <PackageReference Include="coverlet.collector" Version="6.0.4" />
+    <PackageReference Include="Microsoft.NET.Test.Sdk" Version="17.14.1" />
+    <PackageReference Include="xunit" Version="2.9.3" />
+    <PackageReference Include="xunit.runner.visualstudio" Version="3.1.4" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <Using Include="Xunit" />
+    <!-- The tests live under Qubit.NET.*, where the Qubit.NET.Math namespace shadows System.Math. -->
+    <Using Include="System.Math" Alias="Math" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <ProjectReference Include="..\..\src\Qubit.NET\Qubit.NET.csproj" />
+  </ItemGroup>
+
+</Project>
\ No newline at end of file
diff --git a/tests/Qubit.NET.Tests/RegressionTests.cs b/tests/Qubit.NET.Tests/RegressionTests.cs
new file mode 100644
index 0000000..738c79a
--- /dev/null
+++ b/tests/Qubit.NET.Tests/RegressionTests.cs
@@ -0,0 +1,152 @@
+using System.Numerics;
+using Qubit.NET.Circuits;
+using Qubit.NET.Gates;
+using Qubit.NET.Simulation;
+using static QubitNet.Tests.Amplitudes;
+using Qubit.NET;
+
+namespace QubitNet.Tests;
+
+/// <summary>
+/// One test per bug fixed in the 1.0.0 release. Each of these failed before the fix.
+/// </summary>
+public class RegressionTests
+{
+    [Fact]
+    public void Custom_accepts_a_Hadamard_matrix()
+    {
+        // IsUnitary used exact floating-point equality, so (1/sqrt2)^2 + (1/sqrt2)^2
+        // came out as 1.0000000000000002 and every 1/sqrt2 gate was rejected.
+        QuantumCircuit qc = new(1);
+
+        qc.Custom(QuantumGates.H, 0);
+
+        AssertState(qc, InvSqrt2, InvSqrt2);
+    }
+
+    [Fact]
+    public void Run_does_not_consume_the_circuit()
+    {
+        // Run used to strip gates off the circuit with Gates.RemoveAt(0), so the second
+        // call replayed a circuit with no Hadamard and reported 100% |00>.
+        QuantumCircuit qc = BellStates.PhiPlus();
+        qc.Measure();
+
+        var first = Simulator.Run(qc, 500)[0];
+        var second = Simulator.Run(qc, 500)[0];
+
+        foreach (var counts in new[] { first, second })
+        {
+            Assert.Equal(500, counts.Item1.Sum());
+            Assert.True(counts.Item1[0] > 150, "expected roughly half the shots in |00>");
+            Assert.True(counts.Item1[3] > 150, "expected roughly half the shots in |11>");
+            Assert.Equal(0, counts.Item1[1]);
+            Assert.Equal(0, counts.Item1[2]);
+        }
+    }
+
+    [Fact]
+    public void Run_leaves_the_gate_list_intact()
+    {
+        QuantumCircuit qc = BellStates.PhiPlus();
+        qc.Measure();
+
+        string before = qc.ToDiagram();
+        Simulator.Run(qc, 10);
+
+        Assert.Equal(before, qc.ToDiagram());
+    }
+
+    [Fact]
+    public void Copying_a_circuit_leaves_the_original_alone()
+    {
+        // The copy constructor documented a deep copy but shared the gate list and the
+        // modified-qubit flags, so gates applied to the copy leaked into the original.
+        QuantumCircuit original = new(2);
+        original.H(0);
+
+        QuantumCircuit copy = new(original);
+        copy.X(1);
+
+        AssertState(original, InvSqrt2, InvSqrt2, 0, 0);
+        AssertState(copy, 0, 0, InvSqrt2, InvSqrt2);
+    }
+
+    [Fact]
+    public void Copying_a_circuit_does_not_share_its_diagram()
+    {
+        QuantumCircuit original = new(2);
+        original.H(0);
+
+        string before = original.ToDiagram();
+
+        QuantumCircuit copy = new(original);
+        copy.X(1);
+        copy.CNOT(0, 1);
+
+        Assert.Equal(before, original.ToDiagram());
+    }
+
+    [Fact]
+    public void Initialize_rejects_an_unnormalized_state()
+    {
+        // This is the example the README used to ship: |1+i|^2 + |2+2i|^2 = 10.
+        QuantumCircuit qc = new(1);
+
+        Assert.Throws<ArgumentException>(() => qc.Initialize(0, new Complex(1, 1), new Complex(2, 2)));
+    }
+
+    [Fact]
+    public void Initialize_accepts_a_normalized_complex_state()
+    {
+        QuantumCircuit qc = new(1);
+        qc.Initialize(0, new Complex(InvSqrt2, 0), new Complex(0, InvSqrt2));
+
+        AssertNormalized(qc);
+        Assert.Equal(1.0, qc.GetProbabilities().Values.Sum(), Tolerance);
+    }
+
+    [Fact]
+    public void Toffoli_draws_both_controls_as_control_markers()
+    {
+        // The symbol table returned ["@", "+", "+"], so the second control rendered as a
+        // target marker.
+        QuantumCircuit qc = new(3);
+        qc.Toffoli(0, 1, 2);
+
+        string diagram = qc.ToDiagram();
+
+        Assert.Equal(2, diagram.Count(c => c == '@'));
+        Assert.Contains("[+]", diagram);
+    }
+
+    [Fact]
+    public void GetStringResult_handles_an_empty_count_array()
+    {
+        // sb.Remove(sb.Length - 2, 2) threw when no outcome had been recorded.
+        Assert.Equal("{}", (new int[4], 2).GetStringResult());
+    }
+
+    [Fact]
+    public void GetStringResult_formats_counts()
+    {
+        var counts = new int[4];
+        counts[0] = 7;
+        counts[3] = 3;
+
+        Assert.Equal("{'00': 7, '11': 3}", (counts, 2).GetStringResult());
+    }
+
+    [Fact]
+    public void Qubit_count_above_the_maximum_throws_the_documented_exception()
+    {
+        // Used to throw AggregateException, and the documented limit of 30 could never
+        // actually allocate.
+        Assert.Equal(26, QuantumCircuit.MaxQubitCount);
+        Assert.Throws<ArgumentOutOfRangeException>(() => new QuantumCircuit(QuantumCircuit.MaxQubitCount + 1));
+    }
+
+    [Fact]
+    public void Run_returns_nothing_when_the_circuit_has_no_measurement() =>
+        Assert.Empty(Simulator.Run(BellStates.PhiPlus(), 10));
+}
diff --git a/tests/Qubit.NET.Tests/TestHelpers.cs b/tests/Qubit.NET.Tests/TestHelpers.cs
new file mode 100644
index 0000000..578b958
--- /dev/null
+++ b/tests/Qubit.NET.Tests/TestHelpers.cs
@@ -0,0 +1,102 @@
+using System.Numerics;
+using Qubit.NET;
+using Qubit.NET.Utilities;
+
+namespace QubitNet.Tests;
+
+/// <summary>
+/// A deterministic <see cref="IRandomSource"/> that always returns the same value,
+/// so measurement outcomes become predictable.
+/// </summary>
+public sealed class FixedRandomSource(double value) : IRandomSource
+{
+    public double NextDouble() => value;
+}
+
+/// <summary>
+/// An <see cref="IRandomSource"/> backed by a seeded <see cref="Random"/>, for tests that
+/// need a spread of outcomes but must still reproduce across runs.
+/// </summary>
+public sealed class SeededRandomSource(int seed) : IRandomSource
+{
+    private readonly Random _random = new(seed);
+
+    public double NextDouble() => _random.NextDouble();
+}
+
+public static class Amplitudes
+{
+    public const double Tolerance = 1e-9;
+
+    public static readonly double InvSqrt2 = 1.0 / Math.Sqrt(2);
+
+    /// <summary>
+    /// Asserts that the circuit's state vector equals <paramref name="expected"/> entry by entry.
+    /// </summary>
+    public static void AssertState(QuantumCircuit qc, params Complex[] expected)
+    {
+        Assert.Equal(expected.Length, qc.StateVector.Length);
+
+        for (int i = 0; i < expected.Length; i++)
+        {
+            Assert.True(Complex.Abs(qc.StateVector[i] - expected[i]) < Tolerance,
+                $"amplitude[{i}]: expected {expected[i]}, got {qc.StateVector[i]}");
+        }
+    }
+
+    /// <summary>
+    /// Asserts that the state vector is normalized, i.e. the probabilities sum to 1.
+    /// </summary>
+    public static void AssertNormalized(QuantumCircuit qc)
+    {
+        double total = qc.StateVector.Sum(a => a.Real * a.Real + a.Imaginary * a.Imaginary);
+
+        Assert.True(Math.Abs(total - 1.0) < Tolerance, $"state is not normalized: |psi|^2 = {total}");
+    }
+
+    /// <summary>
+    /// Applies a gate matrix to a fresh single-qubit circuit and returns the resulting state.
+    /// </summary>
+    public static Complex[] ApplyToFreshQubit(Complex[,] gate)
+    {
+        QuantumCircuit qc = new(1);
+        qc.Custom(gate, 0);
+        return qc.StateVector;
+    }
+
+    /// <summary>
+    /// Multiplies two square complex matrices.
+    /// </summary>
+    public static Complex[,] Multiply(Complex[,] a, Complex[,] b)
+    {
+        int n = a.GetLength(0);
+        Complex[,] result = new Complex[n, n];
+
+        for (int i = 0; i < n; i++)
+        for (int j = 0; j < n; j++)
+        {
+            Complex sum = Complex.Zero;
+            for (int k = 0; k < n; k++)
+                sum += a[i, k] * b[k, j];
+            result[i, j] = sum;
+        }
+
+        return result;
+    }
+
+    /// <summary>
+    /// Asserts that two square complex matrices are equal entry by entry.
+    /// </summary>
+    public static void AssertMatrixEqual(Complex[,] expected, Complex[,] actual)
+    {
+        int n = expected.GetLength(0);
+        Assert.Equal(n, actual.GetLength(0));
+
+        for (int i = 0; i < n; i++)
+        for (int j = 0; j < n; j++)
+        {
+            Assert.True(Complex.Abs(expected[i, j] - actual[i, j]) < Tolerance,
+                $"[{i},{j}]: expected {expected[i, j]}, got {actual[i, j]}");
+        }
+    }
+}

From 8527e00416c7ef803438ebd4d22927fd7da2e030 Mon Sep 17 00:00:00 2001
From: InfoTCube <infocube.yout@gmail.com>
Date: Mon, 17 Aug 2026 20:34:22 +0200
Subject: [PATCH 4/9] perf: apply gates in place, cutting allocations by ~20x
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Every gate used to allocate a fresh state vector. At 20 qubits that is 16 MB
per gate, so a GHZ chain allocated 328 MB to produce a 16 MB result — the kind
of GC pressure that makes the library unusable inside a Unity frame loop.

- Single-qubit gates apply in place: each amplitude pair is visited once from
  the member whose target bit is 0, so no scratch vector is needed.
- Controlled gates get their own in-place route. CNOT, CY, CZ, CH, the
  controlled rotations and Toffoli all reduce to a 2x2 operator applied only
  where the controls are set, extracted from the trailing block of the gate
  matrix. SWAP, Fredkin and Custom stay on the general path.
- Gate application is parallelized above 2^16 amplitudes.
- The general multi-qubit path no longer calls Array.IndexOf per amplitude,
  and skips zero amplitudes and zero coefficients.

20-qubit GHZ: 83 ms and 328 MB -> 34 ms and 16 MB.

Also collapsed the 24 copy-pasted gate methods onto two helpers, which cuts a
third of QuantumCircuit.cs without touching any signature or doc comment, and
made Gate.TargetQubits/ControlQubits non-nullable. Build warnings: 170 -> 0.

Added ControlledGateTests, which pins the new fast path against the general
matrix path for all eight controlled gates and Toffoli, and a BenchmarkDotNet
project so the README's performance numbers are measured rather than claimed.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Mhwc4XeGtzM2UQx9sNLb3r
---
 .gitignore                                    |   1 +
 CHANGELOG.md                                  |  14 +
 Qubit.NET.sln                                 |  17 +
 README.md                                     |  25 +
 benchmarks/Qubit.NET.Benchmarks/Program.cs    |  61 +++
 .../Qubit.NET.Benchmarks.csproj               |  18 +
 src/Qubit.NET/Gates/Gate.cs                   |  15 +-
 src/Qubit.NET/Math/QuantumMath.cs             | 174 +++++-
 src/Qubit.NET/QuantumCircuit.cs               | 496 ++++--------------
 src/Qubit.NET/QuantumCircuitDrawer.cs         |  12 +-
 src/Qubit.NET/Simulation/Simulator.cs         |  29 +-
 tests/Qubit.NET.Tests/ControlledGateTests.cs  |  96 ++++
 12 files changed, 531 insertions(+), 427 deletions(-)
 create mode 100644 benchmarks/Qubit.NET.Benchmarks/Program.cs
 create mode 100644 benchmarks/Qubit.NET.Benchmarks/Qubit.NET.Benchmarks.csproj
 create mode 100644 tests/Qubit.NET.Tests/ControlledGateTests.cs

diff --git a/.gitignore b/.gitignore
index dbab0a1..87dd546 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,3 +3,4 @@
 artifacts/
 bin/
 obj/
+BenchmarkDotNet.Artifacts/
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4b9e462..781fc65 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -22,9 +22,23 @@ and usable from Unity.
   statistics and collapse, circuit rendering, and a regression test for every bug below.
 - XML documentation for every public gate matrix and for the `State` enum.
 
+### Performance
+
+- Single-qubit and controlled gates are applied **in place**. A circuit now allocates one
+  state vector rather than one per gate: a 20-qubit GHZ chain went from 328 MB of
+  allocations to 16 MB, and from 83 ms to 34 ms.
+- Gate application is parallelized above 2¹⁶ amplitudes.
+- The general multi-qubit path no longer calls `Array.IndexOf` inside its per-amplitude
+  loop, and skips zero amplitudes and zero coefficients.
+- Added a BenchmarkDotNet project under `benchmarks/`.
+
 ### Changed
 
 - **Breaking**: `Examples.Example` is now `Circuits.BellStates`.
+- **Breaking**: `StateVector` may be mutated in place by gate application, so a reference
+  held across a gate call is no longer a snapshot. Clone it if you need one.
+- The 24 near-identical gate methods now delegate to two shared helpers, cutting
+  `QuantumCircuit.cs` by roughly a third with no change to their signatures or docs.
 - **Breaking**: the maximum qubit count is 26, down from a documented 30 that could never
   actually allocate (2³⁰ `Complex` values is 17 GB, and the CLR caps a single array at 2 GB).
 - **Breaking**: `Initialize` now rejects unnormalized states instead of silently producing
diff --git a/Qubit.NET.sln b/Qubit.NET.sln
index 48b71e3..f9f0c90 100644
--- a/Qubit.NET.sln
+++ b/Qubit.NET.sln
@@ -13,6 +13,10 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "tests", "tests", "{0AB3BF05
 EndProject
 Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Qubit.NET.Tests", "tests\Qubit.NET.Tests\Qubit.NET.Tests.csproj", "{A1259BCF-C8A1-4D0C-922F-387222B34D44}"
 EndProject
+Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "benchmarks", "benchmarks", "{66320409-64EC-F7C5-3DEF-65E7510DAAD1}"
+EndProject
+Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Qubit.NET.Benchmarks", "benchmarks\Qubit.NET.Benchmarks\Qubit.NET.Benchmarks.csproj", "{53AB8F9A-2208-400A-8C81-95240B149ECF}"
+EndProject
 Global
 	GlobalSection(SolutionConfigurationPlatforms) = preSolution
 		Debug|Any CPU = Debug|Any CPU
@@ -59,6 +63,18 @@ Global
 		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Release|x64.Build.0 = Release|Any CPU
 		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Release|x86.ActiveCfg = Release|Any CPU
 		{A1259BCF-C8A1-4D0C-922F-387222B34D44}.Release|x86.Build.0 = Release|Any CPU
+		{53AB8F9A-2208-400A-8C81-95240B149ECF}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
+		{53AB8F9A-2208-400A-8C81-95240B149ECF}.Debug|Any CPU.Build.0 = Debug|Any CPU
+		{53AB8F9A-2208-400A-8C81-95240B149ECF}.Debug|x64.ActiveCfg = Debug|Any CPU
+		{53AB8F9A-2208-400A-8C81-95240B149ECF}.Debug|x64.Build.0 = Debug|Any CPU
+		{53AB8F9A-2208-400A-8C81-95240B149ECF}.Debug|x86.ActiveCfg = Debug|Any CPU
+		{53AB8F9A-2208-400A-8C81-95240B149ECF}.Debug|x86.Build.0 = Debug|Any CPU
+		{53AB8F9A-2208-400A-8C81-95240B149ECF}.Release|Any CPU.ActiveCfg = Release|Any CPU
+		{53AB8F9A-2208-400A-8C81-95240B149ECF}.Release|Any CPU.Build.0 = Release|Any CPU
+		{53AB8F9A-2208-400A-8C81-95240B149ECF}.Release|x64.ActiveCfg = Release|Any CPU
+		{53AB8F9A-2208-400A-8C81-95240B149ECF}.Release|x64.Build.0 = Release|Any CPU
+		{53AB8F9A-2208-400A-8C81-95240B149ECF}.Release|x86.ActiveCfg = Release|Any CPU
+		{53AB8F9A-2208-400A-8C81-95240B149ECF}.Release|x86.Build.0 = Release|Any CPU
 	EndGlobalSection
 	GlobalSection(SolutionProperties) = preSolution
 		HideSolutionNode = FALSE
@@ -66,6 +82,7 @@ Global
 	GlobalSection(NestedProjects) = preSolution
 		{F1519400-F1F2-D001-2DB1-E759C4458259} = {827E0CD3-B72D-47B6-A68D-7590B98EB39B}
 		{A1259BCF-C8A1-4D0C-922F-387222B34D44} = {0AB3BF05-4346-4AA6-1389-037BE0695223}
+		{53AB8F9A-2208-400A-8C81-95240B149ECF} = {66320409-64EC-F7C5-3DEF-65E7510DAAD1}
 	EndGlobalSection
 	GlobalSection(ExtensibilityGlobals) = postSolution
 		SolutionGuid = {51965E37-13AF-4AC3-A5D2-1426ACB4DE21}
diff --git a/README.md b/README.md
index 579231c..b516c5f 100644
--- a/README.md
+++ b/README.md
@@ -236,6 +236,31 @@ QuantumCircuit qc = new QuantumCircuit(2);
 qc.RandomSource = new FixedRandomSource();
 ```
 
+## ⚡ Performance
+
+Gates are applied in place, so a circuit allocates one state vector regardless of how many
+gates you apply, and gate application is parallelized above ~16 qubits.
+
+| Circuit                    | 16 qubits | 20 qubits | 22 qubits |
+|----------------------------|----------:|----------:|----------:|
+| Hadamard on every qubit    |    2.4 ms |     36 ms |    147 ms |
+| GHZ (H + CNOT chain)       |    2.3 ms |     34 ms |    142 ms |
+| Measure all qubits         |    2.7 ms |     44 ms |    188 ms |
+
+<sub>BenchmarkDotNet, .NET 10, Ryzen desktop. Reproduce with
+`dotnet run -c Release --project benchmarks/Qubit.NET.Benchmarks`.</sub>
+
+Memory is the real limit — the state vector holds 2ⁿ complex amplitudes at 16 bytes each:
+
+| Qubits | State vector |
+|-------:|-------------:|
+|     16 |         1 MB |
+|     20 |        16 MB |
+|     24 |       256 MB |
+|     26 |    1 GB (max) |
+
+---
+
 ## 📌 Future Roadmap
 
 - [ ] Entanglement entropy measurements
diff --git a/benchmarks/Qubit.NET.Benchmarks/Program.cs b/benchmarks/Qubit.NET.Benchmarks/Program.cs
new file mode 100644
index 0000000..67bef51
--- /dev/null
+++ b/benchmarks/Qubit.NET.Benchmarks/Program.cs
@@ -0,0 +1,61 @@
+using BenchmarkDotNet.Attributes;
+using BenchmarkDotNet.Running;
+using Qubit.NET;
+
+BenchmarkRunner.Run<GateBenchmarks>();
+
+/// <summary>
+/// Measures how gate application scales with register size. The state vector doubles with
+/// every qubit, so these numbers are the practical ceiling on circuit size.
+/// </summary>
+[MemoryDiagnoser]
+public class GateBenchmarks
+{
+    [Params(10, 16, 20, 22)]
+    public int Qubits { get; set; }
+
+    /// <summary>
+    /// A Hadamard on every qubit: the cheapest way to touch the whole state vector once
+    /// per qubit, and the standard opening move of most algorithms.
+    /// </summary>
+    [Benchmark]
+    public QuantumCircuit HadamardLayer()
+    {
+        QuantumCircuit qc = new(Qubits);
+
+        for (int q = 0; q < Qubits; q++)
+            qc.H(q);
+
+        return qc;
+    }
+
+    /// <summary>
+    /// A GHZ state: one Hadamard followed by a chain of CNOTs, exercising the multi-qubit
+    /// path across the full register.
+    /// </summary>
+    [Benchmark]
+    public QuantumCircuit GhzChain()
+    {
+        QuantumCircuit qc = new(Qubits);
+        qc.H(0);
+
+        for (int q = 1; q < Qubits; q++)
+            qc.CNOT(0, q);
+
+        return qc;
+    }
+
+    /// <summary>
+    /// Full-register measurement, which samples and collapses the whole state vector.
+    /// </summary>
+    [Benchmark]
+    public string MeasureAll()
+    {
+        QuantumCircuit qc = new(Qubits);
+
+        for (int q = 0; q < Qubits; q++)
+            qc.H(q);
+
+        return qc.Measure();
+    }
+}
diff --git a/benchmarks/Qubit.NET.Benchmarks/Qubit.NET.Benchmarks.csproj b/benchmarks/Qubit.NET.Benchmarks/Qubit.NET.Benchmarks.csproj
new file mode 100644
index 0000000..5778c88
--- /dev/null
+++ b/benchmarks/Qubit.NET.Benchmarks/Qubit.NET.Benchmarks.csproj
@@ -0,0 +1,18 @@
+<Project Sdk="Microsoft.NET.Sdk">
+
+  <ItemGroup>
+    <ProjectReference Include="..\..\src\Qubit.NET\Qubit.NET.csproj" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <PackageReference Include="BenchmarkDotNet" Version="0.15.8" />
+  </ItemGroup>
+
+  <PropertyGroup>
+    <OutputType>Exe</OutputType>
+    <TargetFramework>net10.0</TargetFramework>
+    <ImplicitUsings>enable</ImplicitUsings>
+    <Nullable>enable</Nullable>
+  </PropertyGroup>
+
+</Project>
diff --git a/src/Qubit.NET/Gates/Gate.cs b/src/Qubit.NET/Gates/Gate.cs
index 92d3b05..bb50982 100644
--- a/src/Qubit.NET/Gates/Gate.cs
+++ b/src/Qubit.NET/Gates/Gate.cs
@@ -10,7 +10,16 @@ namespace Qubit.NET.Gates;
 internal class Gate
 {
     public GateType GateType { get; set; }
+
+    /// <summary>
+    /// The gate's unitary matrix. Null only for <see cref="GateType.Measure"/>, which is
+    /// not a unitary operation.
+    /// </summary>
     public Complex[,]? Matrix { get; set; }
-    public int[]? TargetQubits { get; set; }
-    public int[]? ControlQubits { get; set; }
-}
\ No newline at end of file
+
+    /// <summary>Qubits the gate acts on. Empty is never valid for a recorded gate.</summary>
+    public int[] TargetQubits { get; set; } = [];
+
+    /// <summary>Control qubits, empty for uncontrolled gates.</summary>
+    public int[] ControlQubits { get; set; } = [];
+}
diff --git a/src/Qubit.NET/Math/QuantumMath.cs b/src/Qubit.NET/Math/QuantumMath.cs
index fe48136..1bc6ade 100644
--- a/src/Qubit.NET/Math/QuantumMath.cs
+++ b/src/Qubit.NET/Math/QuantumMath.cs
@@ -16,6 +16,12 @@ internal static class QuantumMath
     /// </summary>
     internal const double Tolerance = 1e-9;
 
+    /// <summary>
+    /// State vector length from which gate application is parallelized. Below this the
+    /// thread pool overhead outweighs the work, so roughly 16 qubits.
+    /// </summary>
+    private const int ParallelThreshold = 1 << 16;
+
     /// <summary>
     /// Applies a single-qubit gate to a specific qubit in a multi-qubit state vector.
     /// </summary>
@@ -32,27 +38,122 @@ internal static Complex[] ApplySingleQubitGate(Complex[] state, Complex[,] gate,
         
         if (gate.GetLength(0) != 2 || gate.GetLength(1) != 2)
             throw new ArgumentException("Single qubit gate must be a 2x2 matrix.");
-    
-        Complex[] newState = new Complex[state.Length];
 
-        for (int i = 0; i < state.Length; i++)
+        // Applied in place. Each amplitude pair (i, i ^ bit) is touched exactly once, from
+        // the member of the pair whose target bit is 0, so no scratch vector is needed.
+        // At 24 qubits a copy would allocate 256 MB per gate.
+        Complex g00 = gate[0, 0], g01 = gate[0, 1], g10 = gate[1, 0], g11 = gate[1, 1];
+        int bitMask = 1 << targetQubit;
+
+        // Below the threshold the thread pool costs more than the work itself.
+        if (state.Length >= ParallelThreshold)
         {
-            int bit = (i >> targetQubit) & 1;
-            int flippedIndex = i ^ (1 << targetQubit);
-            
-            if (bit == 0)
+            Parallel.For(0, state.Length, i =>
             {
+                if ((i & bitMask) != 0) return;
+
+                int flipped = i | bitMask;
                 Complex a = state[i];
-                Complex b = state[flippedIndex];
+                Complex b = state[flipped];
 
-                newState[i] += gate[0, 0] * a + gate[0, 1] * b;
-                newState[flippedIndex] += gate[1, 0] * a + gate[1, 1] * b;
-            }
+                state[i] = g00 * a + g01 * b;
+                state[flipped] = g10 * a + g11 * b;
+            });
+
+            return state;
         }
-        
-        return newState;
+
+        for (int i = 0; i < state.Length; i++)
+        {
+            if ((i & bitMask) != 0) continue;
+
+            int flipped = i | bitMask;
+            Complex a = state[i];
+            Complex b = state[flipped];
+
+            state[i] = g00 * a + g01 * b;
+            state[flipped] = g10 * a + g11 * b;
+        }
+
+        return state;
     }
     
+    /// <summary>
+    /// Extracts the 2x2 operator that a controlled gate applies to its target.
+    /// </summary>
+    /// <param name="gate">
+    /// A controlled gate matrix, which by construction is the identity everywhere except
+    /// its trailing 2x2 block.
+    /// </param>
+    /// <returns>The 2x2 operator applied when every control qubit is set.</returns>
+    internal static Complex[,] ControlledCore(Complex[,] gate)
+    {
+        int n = gate.GetLength(0);
+
+        return new[,]
+        {
+            { gate[n - 2, n - 2], gate[n - 2, n - 1] },
+            { gate[n - 1, n - 2], gate[n - 1, n - 1] }
+        };
+    }
+
+    /// <summary>
+    /// Applies a single-qubit operator to a target qubit, but only on the basis states where
+    /// every control qubit is set. Covers CNOT, CZ, CH, the controlled rotations and Toffoli.
+    /// </summary>
+    /// <param name="state">The current full quantum state vector, modified in place.</param>
+    /// <param name="gate">The 2x2 operator to apply to the target.</param>
+    /// <param name="targetQubit">The index of the target qubit.</param>
+    /// <param name="controlQubits">The indices of the control qubits.</param>
+    /// <returns>The updated quantum state vector.</returns>
+    /// <remarks>
+    /// The general <see cref="ApplyMultiQubitGate"/> path allocates a fresh state vector per
+    /// gate. Controlled gates are the common case in real circuits, and they decompose into
+    /// an in-place pass, so they get their own route.
+    /// </remarks>
+    internal static Complex[] ApplyControlledSingleQubitGate(Complex[] state, Complex[,] gate,
+        int targetQubit, params int[] controlQubits)
+    {
+        Complex g00 = gate[0, 0], g01 = gate[0, 1], g10 = gate[1, 0], g11 = gate[1, 1];
+
+        int targetMask = 1 << targetQubit;
+
+        int controlMask = 0;
+        foreach (int control in controlQubits)
+            controlMask |= 1 << control;
+
+        if (state.Length >= ParallelThreshold)
+        {
+            Parallel.For(0, state.Length, i =>
+            {
+                if ((i & targetMask) != 0 || (i & controlMask) != controlMask) return;
+
+                int flipped = i | targetMask;
+                Complex a = state[i];
+                Complex b = state[flipped];
+
+                state[i] = g00 * a + g01 * b;
+                state[flipped] = g10 * a + g11 * b;
+            });
+
+            return state;
+        }
+
+        for (int i = 0; i < state.Length; i++)
+        {
+            if ((i & targetMask) != 0 || (i & controlMask) != controlMask) continue;
+
+            int flipped = i | targetMask;
+            Complex a = state[i];
+            Complex b = state[flipped];
+
+            state[i] = g00 * a + g01 * b;
+            state[flipped] = g10 * a + g11 * b;
+        }
+
+        return state;
+    }
+
     /// <summary>
     /// Applies a multi-qubit gate to a subset of qubits in a multi-qubit state vector.
     /// </summary>
@@ -64,35 +165,52 @@ public static Complex[] ApplyMultiQubitGate(Complex[] state, Complex[,] gate, in
     {
         Complex[] newState = new Complex[state.Length]; // Initialize to zeros
 
+        int targetCount = targetQubits.Length;
+        int outcomes = 1 << targetCount;
+
+        // Precompute each target's bit mask once. The previous version called
+        // Array.IndexOf inside the per-basis-state loop, turning an O(2^n * 2^k) pass into
+        // O(2^n * k^2 + 2^n * 2^k).
+        int[] masks = new int[targetCount];
+        for (int k = 0; k < targetCount; k++)
+            masks[k] = 1 << targetQubits[k];
+
         // Iterate through each basis state
         for (int i = 0; i < state.Length; i++)
         {
+            Complex amplitude = state[i];
+
+            if (amplitude == Complex.Zero) continue;
+
             // Extract the values of target qubits (as bits)
             int targetBits = 0;
-            foreach (int qubit in targetQubits)
+            for (int k = 0; k < targetCount; k++)
             {
-                targetBits |= ((i >> qubit) & 1) << Array.IndexOf(targetQubits, qubit);
+                if ((i & masks[k]) != 0)
+                    targetBits |= 1 << k;
             }
 
+            // Clearing the target bits once means each destination index is a plain OR.
+            int baseIndex = i;
+            for (int k = 0; k < targetCount; k++)
+                baseIndex &= ~masks[k];
+
             // Apply gate operation to this basis state
-            for (int j = 0; j < (1 << targetQubits.Length); j++)
+            for (int j = 0; j < outcomes; j++)
             {
+                Complex coefficient = gate[j, targetBits];
+
+                if (coefficient == Complex.Zero) continue;
+
                 // Create the new basis state index
-                int newIndex = i;
-                for (int k = 0; k < targetQubits.Length; k++)
+                int newIndex = baseIndex;
+                for (int k = 0; k < targetCount; k++)
                 {
-                    int currentBit = (j >> k) & 1;
-                    int oldBit = (i >> targetQubits[k]) & 1;
-                    if (currentBit != oldBit)
-                    {
-                        newIndex ^= (1 << targetQubits[k]);
-                    }
+                    if ((j & (1 << k)) != 0)
+                        newIndex |= masks[k];
                 }
 
-                // Apply the gate coefficient
-                int gateRow = j;
-                int gateCol = targetBits;
-                newState[newIndex] += gate[gateRow, gateCol] * state[i];
+                newState[newIndex] += coefficient * amplitude;
             }
         }
 
diff --git a/src/Qubit.NET/QuantumCircuit.cs b/src/Qubit.NET/QuantumCircuit.cs
index a6c7dd4..0c3b996 100644
--- a/src/Qubit.NET/QuantumCircuit.cs
+++ b/src/Qubit.NET/QuantumCircuit.cs
@@ -1,4 +1,4 @@
-using System.Numerics;
+using System.Numerics;
 using System.Text;
 using Qubit.NET.Gates;
 using Qubit.NET.Math;
@@ -20,6 +20,10 @@ public class QuantumCircuit
     /// <summary>
     /// State vector representing the current quantum state.
     /// </summary>
+    /// <remarks>
+    /// Single-qubit gates are applied in place, so a reference kept across a gate call may
+    /// observe the updated amplitudes rather than a snapshot. Clone it if you need one.
+    /// </remarks>
     public Complex[] StateVector { get; private set; }
 
     /// <summary>
@@ -217,21 +221,8 @@ public void Initialize(int qubit, Complex alpha, Complex beta, State state = Sta
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
-    public void I(int qubit)
-    {
-        CheckQubit(qubit);
-        
-        Gate iGate = new Gate
-        {
-            GateType = GateType.I,
-            Matrix = QuantumGates.I,
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(iGate);
-        
-        ApplyGate(QuantumGates.I, qubit);
-    }
+    public void I(int qubit) =>
+        ApplySingle(GateType.I, QuantumGates.I, qubit);
     
     /// <summary>
     /// Applies the Hadamard gate (H) to the specified qubit.
@@ -241,21 +232,8 @@ public void I(int qubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
-    public void H(int qubit)
-    {
-        CheckQubit(qubit);
-        
-        Gate hGate = new Gate
-        {
-            GateType = GateType.H,
-            Matrix = QuantumGates.H,
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(hGate);
-        
-        ApplyGate(QuantumGates.H, qubit);
-    }
+    public void H(int qubit) =>
+        ApplySingle(GateType.H, QuantumGates.H, qubit);
 
     /// <summary>
     /// Applies the Pauli-X gate (X) to the specified qubit.
@@ -265,21 +243,8 @@ public void H(int qubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
-    public void X(int qubit)
-    {
-        CheckQubit(qubit);
-        
-        Gate xGate = new Gate
-        {
-            GateType = GateType.X,
-            Matrix = QuantumGates.X,
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(xGate);
-        
-        ApplyGate(QuantumGates.X, qubit);
-    }
+    public void X(int qubit) =>
+        ApplySingle(GateType.X, QuantumGates.X, qubit);
 
     /// <summary>
     /// Applies the Pauli-Y gate (Y) to the specified qubit.
@@ -289,21 +254,8 @@ public void X(int qubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
-    public void Y(int qubit)
-    {
-        CheckQubit(qubit);
-        
-        Gate yGate = new Gate
-        {
-            GateType = GateType.Y,
-            Matrix = QuantumGates.Y,
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(yGate);
-        
-        ApplyGate(QuantumGates.Y, qubit);
-    }
+    public void Y(int qubit) =>
+        ApplySingle(GateType.Y, QuantumGates.Y, qubit);
 
     /// <summary>
     /// Applies the Pauli-Z gate (Z) to the specified qubit.
@@ -313,21 +265,8 @@ public void Y(int qubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
-    public void Z(int qubit)
-    {
-        CheckQubit(qubit);
-        
-        Gate zGate = new Gate
-        {
-            GateType = GateType.Z,
-            Matrix = QuantumGates.Z,
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(zGate);
-        
-        ApplyGate(QuantumGates.Z, qubit);
-    }
+    public void Z(int qubit) =>
+        ApplySingle(GateType.Z, QuantumGates.Z, qubit);
 
     /// <summary>
     /// Applies the S gate (phase gate) to the specified qubit.
@@ -337,21 +276,8 @@ public void Z(int qubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
-    public void S(int qubit)
-    {
-        CheckQubit(qubit);
-        
-        Gate sGate = new Gate
-        {
-            GateType = GateType.S,
-            Matrix = QuantumGates.S,
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(sGate);
-        
-        ApplyGate(QuantumGates.S, qubit);
-    }
+    public void S(int qubit) =>
+        ApplySingle(GateType.S, QuantumGates.S, qubit);
 
     /// <summary>
     /// Applies the S† gate (inverse phase gate) to the specified qubit.
@@ -362,21 +288,8 @@ public void S(int qubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
-    public void Sdag(int qubit)
-    {
-        CheckQubit(qubit);
-        
-        Gate sDagGate = new Gate
-        {
-            GateType = GateType.Sdag,
-            Matrix = QuantumGates.Sdag,
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(sDagGate);
-        
-        ApplyGate(QuantumGates.Sdag, qubit);
-    }
+    public void Sdag(int qubit) =>
+        ApplySingle(GateType.Sdag, QuantumGates.Sdag, qubit);
 
     /// <summary>
     /// Applies the T gate (π/4 phase gate) to the specified qubit.
@@ -386,21 +299,8 @@ public void Sdag(int qubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
-    public void T(int qubit)
-    {
-        CheckQubit(qubit);
-        
-        Gate tGate = new Gate
-        {
-            GateType = GateType.T,
-            Matrix = QuantumGates.T,
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(tGate);
-        
-        ApplyGate(QuantumGates.T, qubit);
-    }
+    public void T(int qubit) =>
+        ApplySingle(GateType.T, QuantumGates.T, qubit);
 
     /// <summary>
     /// Applies the T† gate (inverse of the T gate) to the specified qubit.
@@ -411,21 +311,8 @@ public void T(int qubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
-    public void Tdag(int qubit)
-    {
-        CheckQubit(qubit);
-        
-        Gate tDagGate = new Gate
-        {
-            GateType = GateType.Tdag,
-            Matrix = QuantumGates.Tdag,
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(tDagGate);
-        
-        ApplyGate(QuantumGates.Tdag, qubit);
-    }
+    public void Tdag(int qubit) =>
+        ApplySingle(GateType.Tdag, QuantumGates.Tdag, qubit);
     
     /// <summary>
     /// Applies the Rx gate to the specified qubit.
@@ -437,21 +324,8 @@ public void Tdag(int qubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range, indicating that the specified qubit does not exist in the system.
     /// </exception>
-    public void Rx(int qubit, double theta)
-    {
-        CheckQubit(qubit);
-        
-        Gate rxGate = new Gate
-        {
-            GateType = GateType.Rx,
-            Matrix = QuantumGates.Rx(theta),
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(rxGate);
-        
-        ApplyGate(QuantumGates.Rx(theta), qubit);
-    }
+    public void Rx(int qubit, double theta) =>
+        ApplySingle(GateType.Rx, QuantumGates.Rx(theta), qubit);
     
     /// <summary>
     /// Applies the Ry gate to the specified qubit.
@@ -463,21 +337,8 @@ public void Rx(int qubit, double theta)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range, indicating that the specified qubit does not exist in the system.
     /// </exception>
-    public void Ry(int qubit, double theta)
-    {
-        CheckQubit(qubit);
-        
-        Gate ryGate = new Gate
-        {
-            GateType = GateType.Ry,
-            Matrix = QuantumGates.Ry(theta),
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(ryGate);
-        
-        ApplyGate(QuantumGates.Ry(theta), qubit);
-    }
+    public void Ry(int qubit, double theta) =>
+        ApplySingle(GateType.Ry, QuantumGates.Ry(theta), qubit);
 
     /// <summary>
     /// Applies the Rz gate to the specified qubit.
@@ -489,21 +350,8 @@ public void Ry(int qubit, double theta)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range, indicating that the specified qubit does not exist in the system.
     /// </exception>
-    public void Rz(int qubit, double theta)
-    {
-        CheckQubit(qubit);
-        
-        Gate rzGate = new Gate
-        {
-            GateType = GateType.Rz,
-            Matrix = QuantumGates.Rz(theta),
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(rzGate);
-        
-        ApplyGate(QuantumGates.Rz(theta), qubit);
-    }
+    public void Rz(int qubit, double theta) =>
+        ApplySingle(GateType.Rz, QuantumGates.Rz(theta), qubit);
     
     /// <summary>
     /// Applies the square-root of Pauli-X gate (SX) to the specified qubit.
@@ -515,21 +363,8 @@ public void Rz(int qubit, double theta)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
-    public void SX(int qubit)
-    {
-        CheckQubit(qubit);
-        
-        Gate sxGate = new Gate
-        {
-            GateType = GateType.SX,
-            Matrix = QuantumGates.SX,
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(sxGate);
-        
-        ApplyGate(QuantumGates.SX, qubit);
-    }
+    public void SX(int qubit) =>
+        ApplySingle(GateType.SX, QuantumGates.SX, qubit);
     
     /// <summary>
     /// Applies the square-root of Pauli-Y gate (SY) to the specified qubit.
@@ -541,21 +376,8 @@ public void SX(int qubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
-    public void SY(int qubit)
-    {
-        CheckQubit(qubit);
-        
-        Gate syGate = new Gate
-        {
-            GateType = GateType.SY,
-            Matrix = QuantumGates.SY,
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(syGate);
-        
-        ApplyGate(QuantumGates.SY, qubit);
-    }
+    public void SY(int qubit) =>
+        ApplySingle(GateType.SY, QuantumGates.SY, qubit);
     
     /// <summary>
     /// Applies the square-root of Pauli-Z gate (SZ), also known as the S gate, to the specified qubit.
@@ -567,21 +389,8 @@ public void SY(int qubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range.
     /// </exception>
-    public void SZ(int qubit)
-    {
-        CheckQubit(qubit);
-        
-        Gate syGate = new Gate
-        {
-            GateType = GateType.SZ,
-            Matrix = QuantumGates.SZ,
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(syGate);
-        
-        ApplyGate(QuantumGates.SZ, qubit);
-    }
+    public void SZ(int qubit) =>
+        ApplySingle(GateType.SZ, QuantumGates.SZ, qubit);
     
     /// <summary>
     /// Applies the U3 gate to the specified qubit.
@@ -596,21 +405,8 @@ public void SZ(int qubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if the qubit index is out of range, indicating that the specified qubit does not exist in the system.
     /// </exception>
-    public void U3(int qubit, double theta, double phi, double lambda)
-    {
-        CheckQubit(qubit);
-        
-        Gate u3Gate = new Gate
-        {
-            GateType = GateType.U3,
-            Matrix = QuantumGates.U3(theta, phi, lambda),
-            TargetQubits = [qubit]
-        };
-        
-        Gates.Add(u3Gate);
-        
-        ApplyGate(QuantumGates.U3(theta, phi, lambda), qubit);
-    }
+    public void U3(int qubit, double theta, double phi, double lambda) =>
+        ApplySingle(GateType.U3, QuantumGates.U3(theta, phi, lambda), qubit);
     
     /// <summary>
     /// Applies the CNOT (also called CX) gate (Controlled-NOT) to the specified qubits.
@@ -621,23 +417,8 @@ public void U3(int qubit, double theta, double phi, double lambda)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if any of the qubit indices are out of range.
     /// </exception>
-    public void CNOT(int controlQubit, int targetQubit)
-    {
-        CheckQubit(controlQubit);
-        CheckQubit(targetQubit);
-        
-        Gate cnotGate = new Gate
-        {
-            GateType = GateType.CNOT,
-            Matrix = QuantumGates.CNOT,
-            TargetQubits = [targetQubit],
-            ControlQubits = [controlQubit]
-        };
-        
-        Gates.Add(cnotGate);
-        
-        ApplyGate(QuantumGates.CNOT, [targetQubit, controlQubit]);
-    }
+    public void CNOT(int controlQubit, int targetQubit) =>
+        ApplyControlled(GateType.CNOT, QuantumGates.CNOT, targetQubit, controlQubit);
 
     /// <summary>
     /// Applies the CY gate (Controlled-Y) to the specified qubits.
@@ -648,23 +429,8 @@ public void CNOT(int controlQubit, int targetQubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if any of the qubit indices are out of range.
     /// </exception>
-    public void CY(int controlQubit, int targetQubit)
-    {
-        CheckQubit(controlQubit);
-        CheckQubit(targetQubit);
-        
-        Gate cyGate = new Gate
-        {
-            GateType = GateType.CY,
-            Matrix = QuantumGates.CY,
-            TargetQubits = [targetQubit],
-            ControlQubits = [controlQubit]
-        };
-        
-        Gates.Add(cyGate);
-        
-        ApplyGate(QuantumGates.CY, [targetQubit, controlQubit]);
-    }
+    public void CY(int controlQubit, int targetQubit) =>
+        ApplyControlled(GateType.CY, QuantumGates.CY, targetQubit, controlQubit);
 
     /// <summary>
     /// Applies the CZ gate (Controlled-Z) to the specified qubits.
@@ -675,23 +441,8 @@ public void CY(int controlQubit, int targetQubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if any of the qubit indices are out of range.
     /// </exception>
-    public void CZ(int controlQubit, int targetQubit)
-    {
-        CheckQubit(controlQubit);
-        CheckQubit(targetQubit);
-        
-        Gate czGate = new Gate
-        {
-            GateType = GateType.CZ,
-            Matrix = QuantumGates.CZ,
-            TargetQubits = [targetQubit],
-            ControlQubits = [controlQubit]
-        };
-        
-        Gates.Add(czGate);
-        
-        ApplyGate(QuantumGates.CZ, [targetQubit, controlQubit]);
-    }
+    public void CZ(int controlQubit, int targetQubit) =>
+        ApplyControlled(GateType.CZ, QuantumGates.CZ, targetQubit, controlQubit);
 
     /// <summary>
     /// Applies the CH gate (Controlled-Hadamard) to the specified qubits.
@@ -702,23 +453,8 @@ public void CZ(int controlQubit, int targetQubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if any of the qubit indices are out of range.
     /// </exception>
-    public void CH(int controlQubit, int targetQubit)
-    {
-        CheckQubit(controlQubit);
-        CheckQubit(targetQubit);
-        
-        Gate chGate = new Gate
-        {
-            GateType = GateType.CH,
-            Matrix = QuantumGates.CH,
-            TargetQubits = [targetQubit],
-            ControlQubits = [controlQubit]
-        };
-        
-        Gates.Add(chGate);
-        
-        ApplyGate(QuantumGates.CH, [targetQubit, controlQubit]);
-    }
+    public void CH(int controlQubit, int targetQubit) =>
+        ApplyControlled(GateType.CH, QuantumGates.CH, targetQubit, controlQubit);
     
     /// <summary>
     /// Applies the CRx gate (Controlled-Rx) to the specified qubits.
@@ -730,23 +466,8 @@ public void CH(int controlQubit, int targetQubit)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if any of the qubit indices are out of range.
     /// </exception>
-    public void CRx(int controlQubit, int targetQubit, double theta)
-    {
-        CheckQubit(controlQubit);
-        CheckQubit(targetQubit);
-        
-        Gate crxGate = new Gate
-        {
-            GateType = GateType.CRx,
-            Matrix = QuantumGates.CRx(theta),
-            TargetQubits = [targetQubit],
-            ControlQubits = [controlQubit]
-        };
-        
-        Gates.Add(crxGate);
-        
-        ApplyGate(QuantumGates.CRx(theta), [targetQubit, controlQubit]);
-    }
+    public void CRx(int controlQubit, int targetQubit, double theta) =>
+        ApplyControlled(GateType.CRx, QuantumGates.CRx(theta), targetQubit, controlQubit);
     
     /// <summary>
     /// Applies the CRy gate (Controlled-Ry) to the specified qubits.
@@ -758,23 +479,8 @@ public void CRx(int controlQubit, int targetQubit, double theta)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if any of the qubit indices are out of range.
     /// </exception>
-    public void CRy(int controlQubit, int targetQubit, double theta)
-    {
-        CheckQubit(controlQubit);
-        CheckQubit(targetQubit);
-        
-        Gate cryGate = new Gate
-        {
-            GateType = GateType.CRy,
-            Matrix = QuantumGates.CRy(theta),
-            TargetQubits = [targetQubit],
-            ControlQubits = [controlQubit]
-        };
-        
-        Gates.Add(cryGate);
-        
-        ApplyGate(QuantumGates.CRy(theta), [targetQubit, controlQubit]);
-    }
+    public void CRy(int controlQubit, int targetQubit, double theta) =>
+        ApplyControlled(GateType.CRy, QuantumGates.CRy(theta), targetQubit, controlQubit);
     
     /// <summary>
     /// Applies the CRz gate (Controlled-Rz) to the specified qubits.
@@ -786,23 +492,8 @@ public void CRy(int controlQubit, int targetQubit, double theta)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if any of the qubit indices are out of range.
     /// </exception>
-    public void CRz(int controlQubit, int targetQubit, double theta)
-    {
-        CheckQubit(controlQubit);
-        CheckQubit(targetQubit);
-        
-        Gate crzGate = new Gate
-        {
-            GateType = GateType.CRz,
-            Matrix = QuantumGates.CRz(theta),
-            TargetQubits = [targetQubit],
-            ControlQubits = [controlQubit]
-        };
-        
-        Gates.Add(crzGate);
-        
-        ApplyGate(QuantumGates.CRz(theta), [targetQubit, controlQubit]);
-    }
+    public void CRz(int controlQubit, int targetQubit, double theta) =>
+        ApplyControlled(GateType.CRz, QuantumGates.CRz(theta), targetQubit, controlQubit);
     
     /// <summary>
     /// Applies the CU3 gate to the specified qubits.
@@ -816,23 +507,8 @@ public void CRz(int controlQubit, int targetQubit, double theta)
     /// <exception cref="QubitIndexOutOfRangeException">
     /// Thrown if any of the qubit indices are out of range.
     /// </exception>
-    public void CU3(int controlQubit, int targetQubit, double theta, double phi, double lambda)
-    {
-        CheckQubit(controlQubit);
-        CheckQubit(targetQubit);
-        
-        Gate cu3Gate = new Gate
-        {
-            GateType = GateType.CU3,
-            Matrix = QuantumGates.CU3(theta, phi, lambda),
-            TargetQubits = [targetQubit],
-            ControlQubits = [controlQubit]
-        };
-        
-        Gates.Add(cu3Gate);
-        
-        ApplyGate(QuantumGates.CU3(theta, phi, lambda), [targetQubit, controlQubit]);
-    }
+    public void CU3(int controlQubit, int targetQubit, double theta, double phi, double lambda) =>
+        ApplyControlled(GateType.CU3, QuantumGates.CU3(theta, phi, lambda), targetQubit, controlQubit);
 
     /// <summary>
     /// Applies the SWAP gate to the specified qubits, exchanging their states.
@@ -875,6 +551,8 @@ public void Toffoli(int firstControlQubit, int secondControlQubit, int targetQub
         CheckQubit(secondControlQubit);
         CheckQubit(targetQubit);
         
+        CheckIfDifferentQubits([firstControlQubit, secondControlQubit, targetQubit]);
+
         Gate toffoliGate = new Gate
         {
             GateType = GateType.Toffoli,
@@ -882,10 +560,16 @@ public void Toffoli(int firstControlQubit, int secondControlQubit, int targetQub
             TargetQubits = [targetQubit],
             ControlQubits = [secondControlQubit, firstControlQubit]
         };
-        
+
         Gates.Add(toffoliGate);
-        
-        ApplyGate(QuantumGates.Toffoli, [targetQubit, firstControlQubit, secondControlQubit]);
+
+        StateVector = QuantumMath.ApplyControlledSingleQubitGate(
+            StateVector, QuantumMath.ControlledCore(QuantumGates.Toffoli),
+            targetQubit, firstControlQubit, secondControlQubit);
+
+        _isQubitModified[targetQubit] = true;
+        _isQubitModified[firstControlQubit] = true;
+        _isQubitModified[secondControlQubit] = true;
     }
 
     /// <summary>
@@ -1028,6 +712,58 @@ public override string ToString()
         return sb.ToString();
     }
     
+    /// <summary>
+    /// Validates, records and applies a single-qubit gate. Shared by every one-qubit
+    /// gate method, which are otherwise identical apart from their type and matrix.
+    /// </summary>
+    /// <param name="type">The gate type, recorded for replay and for drawing.</param>
+    /// <param name="matrix">The 2x2 unitary matrix representing the gate.</param>
+    /// <param name="qubit">The index of the qubit to apply the gate to.</param>
+    private void ApplySingle(GateType type, Complex[,] matrix, int qubit)
+    {
+        CheckQubit(qubit);
+
+        Gates.Add(new Gate
+        {
+            GateType = type,
+            Matrix = matrix,
+            TargetQubits = [qubit]
+        });
+
+        ApplyGate(matrix, qubit);
+    }
+
+    /// <summary>
+    /// Validates, records and applies a two-qubit controlled gate. Shared by every
+    /// controlled gate method.
+    /// </summary>
+    /// <param name="type">The gate type, recorded for replay and for drawing.</param>
+    /// <param name="matrix">The 4x4 unitary matrix representing the gate.</param>
+    /// <param name="targetQubit">The index of the target qubit.</param>
+    /// <param name="controlQubit">The index of the control qubit.</param>
+    private void ApplyControlled(GateType type, Complex[,] matrix, int targetQubit, int controlQubit)
+    {
+        CheckQubit(controlQubit);
+        CheckQubit(targetQubit);
+
+        if (targetQubit == controlQubit)
+            throw new ArgumentException("Every gate argument must be a different Qubit");
+
+        Gates.Add(new Gate
+        {
+            GateType = type,
+            Matrix = matrix,
+            TargetQubits = [targetQubit],
+            ControlQubits = [controlQubit]
+        });
+
+        StateVector = QuantumMath.ApplyControlledSingleQubitGate(
+            StateVector, QuantumMath.ControlledCore(matrix), targetQubit, controlQubit);
+
+        _isQubitModified[targetQubit] = true;
+        _isQubitModified[controlQubit] = true;
+    }
+
     /// <summary>
     /// Applies a single-qubit gate given by a unitary matrix to the specified qubit.
     /// </summary>
diff --git a/src/Qubit.NET/QuantumCircuitDrawer.cs b/src/Qubit.NET/QuantumCircuitDrawer.cs
index d246edc..052cb60 100644
--- a/src/Qubit.NET/QuantumCircuitDrawer.cs
+++ b/src/Qubit.NET/QuantumCircuitDrawer.cs
@@ -115,15 +115,17 @@ private static int AssignGatePositions(IList<IList<(string, int)>> gatePositions
         
         foreach (var gate in circuit.Gates)
         {
-            var controlQubits = gate.ControlQubits ?? Array.Empty<int>();
-            var targetQubits = gate.TargetQubits ?? Array.Empty<int>();
+            var controlQubits = gate.ControlQubits;
+            var targetQubits = gate.TargetQubits;
 
             string[] reps;
-            
+
+            // Measure and Custom span an arbitrary number of qubits, so their symbols are
+            // repeated to match rather than looked up in the fixed symbol table.
             if(gate.GateType == GateType.Measure)
-                reps = Enumerable.Repeat("M", gate.TargetQubits.Length).ToArray();
+                reps = Enumerable.Repeat("M", targetQubits.Length).ToArray();
             else if(gate.GateType == GateType.Custom)
-                reps = Enumerable.Repeat("C", gate.TargetQubits.Length).ToArray();
+                reps = Enumerable.Repeat("C", targetQubits.Length).ToArray();
             else
                 reps = Helpers.GateTypeToCharRepresentation(gate.GateType);
             
diff --git a/src/Qubit.NET/Simulation/Simulator.cs b/src/Qubit.NET/Simulation/Simulator.cs
index 4ff654a..b2d5842 100644
--- a/src/Qubit.NET/Simulation/Simulator.cs
+++ b/src/Qubit.NET/Simulation/Simulator.cs
@@ -124,32 +124,39 @@ public static string GetStringResult(this (int[], int) result)
     /// <param name="gateType">The type of gate being applied.</param>
     private static Complex[] ApplyGate(Complex[] stateVector, Gate currentGate, GateType gateType)
     {
+        // Measure gates are handled by the caller and never reach here, so every gate at
+        // this point must carry a matrix.
+        Complex[,] matrix = currentGate.Matrix
+                            ?? throw new InvalidOperationException($"Gate {gateType} has no matrix.");
+
         switch (gateType)
         {
             case GateType.I or GateType.H or GateType.X or GateType.Y or GateType.Z or GateType.S or GateType.Sdag 
                 or GateType.T or GateType.Tdag or GateType.Rx or GateType.Ry or GateType.Rz or GateType.SX 
                 or GateType.SY or GateType.SZ or GateType.U3:
-                stateVector = QuantumMath.ApplySingleQubitGate(stateVector, currentGate.Matrix, currentGate.TargetQubits.First());
+                stateVector = QuantumMath.ApplySingleQubitGate(stateVector, matrix, currentGate.TargetQubits[0]);
                 break;
-            case GateType.CNOT or GateType.CY or GateType.CZ or GateType.CH or GateType.CRx or GateType.CRy 
+            case GateType.CNOT or GateType.CY or GateType.CZ or GateType.CH or GateType.CRx or GateType.CRy
                 or GateType.CRz or GateType.CU3:
-                stateVector = QuantumMath.ApplyMultiQubitGate(stateVector, currentGate.Matrix, 
-                    [currentGate.TargetQubits.First(), currentGate.ControlQubits.First()]);
+                stateVector = QuantumMath.ApplyControlledSingleQubitGate(stateVector,
+                    QuantumMath.ControlledCore(matrix),
+                    currentGate.TargetQubits[0], currentGate.ControlQubits[0]);
                 break;
             case GateType.SWAP:
-                stateVector = QuantumMath.ApplyMultiQubitGate(stateVector, currentGate.Matrix,
-                    currentGate.TargetQubits.ToArray());
+                stateVector = QuantumMath.ApplyMultiQubitGate(stateVector, matrix,
+                    currentGate.TargetQubits);
                 break;
             case GateType.Toffoli:
-                stateVector = QuantumMath.ApplyMultiQubitGate(stateVector, currentGate.Matrix,
-                    [currentGate.TargetQubits.First(), currentGate.ControlQubits[0], currentGate.ControlQubits[1]]);
+                stateVector = QuantumMath.ApplyControlledSingleQubitGate(stateVector,
+                    QuantumMath.ControlledCore(matrix),
+                    currentGate.TargetQubits[0], currentGate.ControlQubits[0], currentGate.ControlQubits[1]);
                 break;
             case GateType.Fredkin:
-                stateVector = QuantumMath.ApplyMultiQubitGate(stateVector, currentGate.Matrix,
-                    [currentGate.TargetQubits[1], currentGate.TargetQubits[0], currentGate.ControlQubits.First()]);
+                stateVector = QuantumMath.ApplyMultiQubitGate(stateVector, matrix,
+                    [currentGate.TargetQubits[1], currentGate.TargetQubits[0], currentGate.ControlQubits[0]]);
                 break;
             case GateType.Custom:
-                stateVector = QuantumMath.ApplyMultiQubitGate(stateVector, currentGate.Matrix,
+                stateVector = QuantumMath.ApplyMultiQubitGate(stateVector, matrix,
                     Enumerable.Reverse(currentGate.TargetQubits).ToArray());
                 break;
         }
diff --git a/tests/Qubit.NET.Tests/ControlledGateTests.cs b/tests/Qubit.NET.Tests/ControlledGateTests.cs
new file mode 100644
index 0000000..2914baa
--- /dev/null
+++ b/tests/Qubit.NET.Tests/ControlledGateTests.cs
@@ -0,0 +1,96 @@
+using System.Numerics;
+using Qubit.NET;
+using Qubit.NET.Gates;
+using static QubitNet.Tests.Amplitudes;
+
+namespace QubitNet.Tests;
+
+/// <summary>
+/// Controlled gates take a specialized in-place route rather than the general multi-qubit
+/// path. These tests pin the two against each other so the fast path cannot drift.
+/// </summary>
+public class ControlledGateTests
+{
+    public static TheoryData<string, Action<QuantumCircuit>, Complex[,]> ControlledGates() => new()
+    {
+        { "CNOT", qc => qc.CNOT(0, 1), QuantumGates.CNOT },
+        { "CY", qc => qc.CY(0, 1), QuantumGates.CY },
+        { "CZ", qc => qc.CZ(0, 1), QuantumGates.CZ },
+        { "CH", qc => qc.CH(0, 1), QuantumGates.CH },
+        { "CRx", qc => qc.CRx(0, 1, 0.7), QuantumGates.CRx(0.7) },
+        { "CRy", qc => qc.CRy(0, 1, 0.7), QuantumGates.CRy(0.7) },
+        { "CRz", qc => qc.CRz(0, 1, 0.7), QuantumGates.CRz(0.7) },
+        { "CU3", qc => qc.CU3(0, 1, 0.3, 0.5, 0.7), QuantumGates.CU3(0.3, 0.5, 0.7) },
+    };
+
+    [Theory]
+    [MemberData(nameof(ControlledGates))]
+    public void Fast_path_matches_the_general_matrix_path(string name, Action<QuantumCircuit> apply, Complex[,] matrix)
+    {
+        // Start from a state with amplitude on every basis state, so any disagreement shows.
+        QuantumCircuit viaGate = Spread();
+        apply(viaGate);
+
+        QuantumCircuit viaCustom = Spread();
+        viaCustom.Custom(matrix, 0, 1);
+
+        for (int i = 0; i < viaGate.StateVector.Length; i++)
+        {
+            Assert.True(Complex.Abs(viaGate.StateVector[i] - viaCustom.StateVector[i]) < Tolerance,
+                $"{name} amplitude[{i}]: fast path {viaGate.StateVector[i]}, general path {viaCustom.StateVector[i]}");
+        }
+    }
+
+    [Fact]
+    public void Toffoli_fast_path_matches_the_general_matrix_path()
+    {
+        QuantumCircuit viaGate = Spread(3);
+        viaGate.Toffoli(0, 1, 2);
+
+        QuantumCircuit viaCustom = Spread(3);
+        viaCustom.Custom(QuantumGates.Toffoli, 0, 1, 2);
+
+        for (int i = 0; i < viaGate.StateVector.Length; i++)
+        {
+            Assert.True(Complex.Abs(viaGate.StateVector[i] - viaCustom.StateVector[i]) < Tolerance,
+                $"amplitude[{i}]: fast path {viaGate.StateVector[i]}, general path {viaCustom.StateVector[i]}");
+        }
+    }
+
+    [Fact]
+    public void Controlled_gates_stay_normalized()
+    {
+        QuantumCircuit qc = Spread(3);
+        qc.CNOT(0, 1);
+        qc.CH(1, 2);
+        qc.CRy(0, 2, 1.1);
+        qc.Toffoli(0, 1, 2);
+
+        AssertNormalized(qc);
+    }
+
+    [Fact]
+    public void A_controlled_gate_cannot_use_one_qubit_as_both_control_and_target() =>
+        Assert.Throws<ArgumentException>(() => new QuantumCircuit(2).CZ(1, 1));
+
+    [Fact]
+    public void Toffoli_rejects_a_repeated_qubit() =>
+        Assert.Throws<ArgumentException>(() => new QuantumCircuit(3).Toffoli(0, 1, 1));
+
+    /// <summary>
+    /// Builds a circuit whose every basis state carries a distinct, non-trivial amplitude.
+    /// </summary>
+    private static QuantumCircuit Spread(int qubits = 2)
+    {
+        QuantumCircuit qc = new(qubits);
+
+        for (int q = 0; q < qubits; q++)
+        {
+            qc.H(q);
+            qc.Rz(q, 0.3 * (q + 1));
+            qc.Ry(q, 0.4 * (q + 1));
+        }
+
+        return qc;
+    }
+}

From 17b3237b05869e151adb4c43f9ae420e1470bb51 Mon Sep 17 00:00:00 2001
From: InfoTCube <infocube.yout@gmail.com>
Date: Mon, 17 Aug 2026 20:41:12 +0200
Subject: [PATCH 5/9] feat: classical bits, feedforward and a real
 MeasurementResult type

Adds the one thing Qiskit's ClassicalRegister buys that this library could not
express: conditioning a gate on a measurement outcome. Bit layout, one result
per Measure call, and QASM's creg mapping were already covered.

- MeasureInto(qubit, classicalBit) stores an outcome; ClassicalBit(i) reads it;
  When(bit, value, body) runs the body's gates only when the bit matches. Plain
  Measure fills the classical bit matching each qubit it measured.
- Conditional gates are always recorded, so Simulator.Run re-evaluates the
  condition per shot against that shot's own outcomes.
- No separate ClassicalRegister type: it would break every call site and tax the
  common case of "measure everything, get a bitstring" for no gain here.
- Simulator.Run now returns MeasurementResult, with Counts, Shots,
  Probability(outcome) and MostFrequent, instead of an (int[], int) tuple. It
  also no longer allocates 2^n counters for a 2-qubit partial measurement.
- Added QuantumCircuit.Reset().
- The drawer lays a conditional gate out after the measurement that wrote its
  bit, and connects the two, instead of showing the correction before its cause.

Teleportation is the proof: tested across five input states and twenty
measurement branches, the message always arrives on the third qubit.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Mhwc4XeGtzM2UQx9sNLb3r
---
 CHANGELOG.md                                  |  11 +
 Qubit.NET.Demo/Program.cs                     |  63 ++++-
 README.md                                     |  53 +++-
 src/Qubit.NET/Gates/Gate.cs                   |  17 ++
 src/Qubit.NET/QuantumCircuit.cs               | 186 +++++++++++-
 src/Qubit.NET/QuantumCircuitDrawer.cs         |  27 +-
 src/Qubit.NET/Simulation/MeasurementResult.cs | 107 +++++++
 src/Qubit.NET/Simulation/Simulator.cs         |  88 +++---
 tests/Qubit.NET.Tests/FeedforwardTests.cs     | 266 ++++++++++++++++++
 tests/Qubit.NET.Tests/RegressionTests.cs      |  35 ++-
 10 files changed, 770 insertions(+), 83 deletions(-)
 create mode 100644 src/Qubit.NET/Simulation/MeasurementResult.cs
 create mode 100644 tests/Qubit.NET.Tests/FeedforwardTests.cs

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 781fc65..53e571b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -17,6 +17,14 @@ and usable from Unity.
   drawer works without a console (Unity, ASP.NET, tests). `Draw()` still prints in color.
 - `QuantumGates.Format(matrix)` — the string-returning counterpart of `Print`.
 - `QuantumCircuit.MaxQubitCount` constant.
+- **Classical bits and feedforward**: `MeasureInto(qubit, classicalBit)` stores an outcome,
+  `ClassicalBit(i)` reads it back, and `When(classicalBit, value, body)` conditions gates on
+  it. Plain `Measure` fills the classical bit matching each qubit it measured. This is what
+  quantum teleportation, superdense coding and error correction need; there is deliberately
+  no separate `ClassicalRegister` type.
+- `MeasurementResult`, replacing the `(int[], int)` tuple returned by `Simulator.Run`. It
+  exposes `Counts`, `Shots`, `Probability(outcome)` and `MostFrequent`.
+- `QuantumCircuit.Reset()`, returning a circuit to |0…0⟩ with its gate list cleared.
 - Continuous integration and a tag-driven NuGet release workflow.
 - An xUnit test suite covering gate algebra, the textbook Bell and GHZ states, measurement
   statistics and collapse, circuit rendering, and a regression test for every bug below.
@@ -34,6 +42,9 @@ and usable from Unity.
 
 ### Changed
 
+- **Breaking**: `Simulator.Run` returns `IList<MeasurementResult>` rather than
+  `IList<(int[], int)>`. `GetStringResult()` still works, and `result.ToString()` now gives
+  the same text.
 - **Breaking**: `Examples.Example` is now `Circuits.BellStates`.
 - **Breaking**: `StateVector` may be mutated in place by gate application, so a reference
   held across a gate call is no longer a snapshot. Clone it if you need one.
diff --git a/Qubit.NET.Demo/Program.cs b/Qubit.NET.Demo/Program.cs
index 1a0bf17..00cb136 100644
--- a/Qubit.NET.Demo/Program.cs
+++ b/Qubit.NET.Demo/Program.cs
@@ -1,4 +1,4 @@
-using Qubit.NET.Gates;
+using Qubit.NET.Gates;
 using Qubit.NET.Simulation;
 
 namespace Qubit.NET.Demo;
@@ -7,25 +7,64 @@ public static class Program
 {
     public static void Main(string[] args)
     {
+        Console.OutputEncoding = System.Text.Encoding.UTF8;
+
+        BellState();
+        Teleportation();
+    }
+
+    private static void BellState()
+    {
+        Console.WriteLine("=== Bell state ===\n");
+
         QuantumCircuit qc = new QuantumCircuit(2);
-        
+
         qc.H(0);
         qc.CNOT(0, 1);
 
-        Console.WriteLine(qc.ToString());
+        Console.WriteLine($"State:  {qc}");
 
-        qc.Measure();
-        
         QuantumGates.Print(QuantumGates.H);
-        
+
+        qc.Measure();
         qc.Draw();
 
-        var results = Simulator.Run(qc, 1000);
+        MeasurementResult result = Simulator.Run(qc, 1000)[0];
 
-        foreach (var result in results)
-        {
-            Console.WriteLine(result.GetStringResult());
-        }
+        Console.WriteLine($"\nCounts:      {result}");
+        Console.WriteLine($"P(11):       {result.Probability("11"):F3}");
+        Console.WriteLine($"Most common: {result.MostFrequent}\n");
     }
-}
 
+    private static void Teleportation()
+    {
+        Console.WriteLine("=== Teleportation ===\n");
+
+        // Qubit 0 carries the message, qubits 1 and 2 share a Bell pair.
+        QuantumCircuit qc = new QuantumCircuit(3);
+
+        qc.Ry(0, 0.8);
+
+        qc.H(1);
+        qc.CNOT(1, 2);
+
+        qc.CNOT(0, 1);
+        qc.H(0);
+        qc.MeasureInto(0, 0);
+        qc.MeasureInto(1, 1);
+
+        // Corrections conditioned on the two measurement outcomes.
+        qc.When(1, 1, c => c.X(2));
+        qc.When(0, 1, c => c.Z(2));
+
+        qc.Draw();
+
+        double expected = System.Math.Cos(0.4) * System.Math.Cos(0.4);
+        double actual = qc.GetProbabilities()
+            .Where(p => p.Key[0] == '0')   // qubit 2 is the leading bit
+            .Sum(p => p.Value);
+
+        Console.WriteLine($"\nMeasured bits: c0={qc.ClassicalBit(0)} c1={qc.ClassicalBit(1)}");
+        Console.WriteLine($"P(qubit 2 = |0>): {actual:F6}  (expected {expected:F6})");
+    }
+}
diff --git a/README.md b/README.md
index b516c5f..398e80f 100644
--- a/README.md
+++ b/README.md
@@ -210,10 +210,59 @@ qc.H(0);
 qc.CNOT(0, 1);
 qc.Measure();
 
-string results = Simulator.Run(qc, 1000)[0].GetStringResult();
-Console.WriteLine(results);
+MeasurementResult result = Simulator.Run(qc, 1000)[0];
+
+Console.WriteLine(result);                    // {'00': 512, '11': 488}
+Console.WriteLine(result.Counts["00"]);       // 512
+Console.WriteLine(result.Probability("11"));  // 0.488
+Console.WriteLine(result.MostFrequent);       // 00
+Console.WriteLine(result.Shots);              // 1000
 ```
 
+`Simulator.Run` never modifies the circuit you hand it, so you can run the same circuit as
+many times as you like.
+
+---
+
+### 🔀 Classical bits and feedforward
+
+Measuring writes into a classical bit — one per qubit, so measuring qubit `q` fills bit `q`
+unless you say otherwise with `MeasureInto`. `When` then conditions later gates on that bit,
+which is what mid-circuit measurement and error correction need.
+
+```csharp
+qc.MeasureInto(qubit: 0, classicalBit: 0);
+qc.When(classicalBit: 0, value: 1, c => c.X(2));   // X(2) runs only if bit 0 came out 1
+```
+
+Conditional gates are always *recorded*, so `Simulator.Run` re-evaluates the condition on
+every shot against that shot's own outcomes.
+
+Quantum teleportation in full:
+
+```csharp
+var qc = new QuantumCircuit(3);
+
+qc.Ry(0, theta);        // the message on qubit 0
+
+qc.H(1);                // entangle qubits 1 and 2
+qc.CNOT(1, 2);
+
+qc.CNOT(0, 1);          // Bell-basis measurement of qubits 0 and 1
+qc.H(0);
+qc.MeasureInto(0, 0);
+qc.MeasureInto(1, 1);
+
+qc.When(1, 1, c => c.X(2));   // corrections
+qc.When(0, 1, c => c.Z(2));
+
+// qubit 2 now holds the state qubit 0 started in
+```
+
+> Qubit.NET deliberately has no separate `ClassicalRegister` type. The classical register in
+> Qiskit exists mainly to express feedforward and result layout; `When` and `MeasureInto`
+> cover both without making every circuit declare two registers up front.
+
 ---
 
  ### 🎲 Randomness source
diff --git a/src/Qubit.NET/Gates/Gate.cs b/src/Qubit.NET/Gates/Gate.cs
index bb50982..879c2c2 100644
--- a/src/Qubit.NET/Gates/Gate.cs
+++ b/src/Qubit.NET/Gates/Gate.cs
@@ -22,4 +22,21 @@ internal class Gate
 
     /// <summary>Control qubits, empty for uncontrolled gates.</summary>
     public int[] ControlQubits { get; set; } = [];
+
+    /// <summary>
+    /// Classical bit this gate is conditioned on, or null if it always runs.
+    /// Set by <see cref="QuantumCircuit.When"/> to express classical feedforward.
+    /// </summary>
+    public int? ConditionBit { get; set; }
+
+    /// <summary>
+    /// Value <see cref="ConditionBit"/> must hold for the gate to run.
+    /// </summary>
+    public int ConditionValue { get; set; }
+
+    /// <summary>
+    /// For measurement gates, the classical bit each measured qubit writes into,
+    /// positionally matching <see cref="TargetQubits"/>.
+    /// </summary>
+    public int[] ClassicalBits { get; set; } = [];
 }
diff --git a/src/Qubit.NET/QuantumCircuit.cs b/src/Qubit.NET/QuantumCircuit.cs
index 0c3b996..9ed6533 100644
--- a/src/Qubit.NET/QuantumCircuit.cs
+++ b/src/Qubit.NET/QuantumCircuit.cs
@@ -50,6 +50,17 @@ public class QuantumCircuit
     /// </summary>
     private readonly bool[] _isQubitModified;
 
+    /// <summary>
+    /// Classical bits holding measurement outcomes. One per qubit: measuring qubit q writes
+    /// into classical bit q unless <see cref="MeasureInto"/> says otherwise.
+    /// </summary>
+    private readonly int[] _classicalBits;
+
+    /// <summary>
+    /// The condition currently in force, set while inside a <see cref="When"/> body.
+    /// </summary>
+    private (int Bit, int Value)? _activeCondition;
+
     /// <summary>
     /// Largest supported qubit count. The state vector holds 2^n <see cref="Complex"/>
     /// values at 16 bytes each, and the CLR caps a single array at 2 GB, so 2^26 elements
@@ -74,6 +85,7 @@ public QuantumCircuit(int qubitCount)
         StateVector = new Complex[1 << qubitCount];
         StateVector[0] = new Complex(1, 0);
         _isQubitModified = new bool[qubitCount];
+        _classicalBits = new int[qubitCount];
     }
 
     /// <summary>
@@ -91,6 +103,7 @@ public QuantumCircuit(QuantumCircuit qc)
         Gates = new List<Gate>(qc.Gates);
         Initializations = new List<InitialState>(qc.Initializations);
         _isQubitModified = (bool[])qc._isQubitModified.Clone();
+        _classicalBits = (int[])qc._classicalBits.Clone();
         RandomSource = qc.RandomSource;
     }
 
@@ -107,32 +120,143 @@ public string Measure(params int[] qubits)
         
         CheckIfDifferentQubits(qubits);
         
+        // Each measured qubit writes into the classical bit of the same index, so that
+        // When() can condition on it without any extra bookkeeping.
+        int[] measured = qubits.Length == 0 ? Enumerable.Range(0, QubitCount).ToArray() : qubits;
+
         Gate measure = new Gate()
         {
             GateType = GateType.Measure,
-            TargetQubits = qubits.Length == 0 ? Enumerable.Range(0, QubitCount).ToArray() : qubits
+            TargetQubits = measured,
+            ClassicalBits = measured
         };
-        
+
         Gates.Add(measure);
-        
+
         // If measuring all qubits, use much more efficient method
         if (qubits.Length == 0)
         {
             // Perform a measurement by sampling from the current state vector probabilities
             var result = QuantumMath.SampleMeasurement(StateVector, RandomSource);
-    
+
             // Collapse the quantum state to the measured state (collapse the superposition)
             StateVector = QuantumMath.CollapseToState(StateVector, result);
-            
+
+            for (int q = 0; q < QubitCount; q++)
+                _classicalBits[q] = (result >> q) & 1;
+
             return Convert.ToString(result, 2).PadLeft(QubitCount, '0');
         }
-        
+
         var partialResult = QuantumMath.SamplePartialMeasurement(StateVector, qubits, RandomSource);
         StateVector = QuantumMath.CollapseToPartialMeasurement(StateVector, qubits, partialResult);
-            
+
+        // SamplePartialMeasurement packs bits with the first listed qubit most significant.
+        for (int b = 0; b < qubits.Length; b++)
+            _classicalBits[qubits[b]] = (partialResult >> (qubits.Length - 1 - b)) & 1;
+
         return Convert.ToString(partialResult, 2).PadLeft(qubits.Length, '0');
     }
 
+    /// <summary>
+    /// Measures a single qubit and stores the outcome in the given classical bit, so that
+    /// later gates can be conditioned on it with <see cref="When"/>.
+    /// </summary>
+    /// <param name="qubit">The index of the qubit to measure.</param>
+    /// <param name="classicalBit">The classical bit to store the outcome in.</param>
+    /// <returns>The measured value, 0 or 1.</returns>
+    /// <exception cref="QubitIndexOutOfRangeException">
+    /// Thrown if either index is out of range.
+    /// </exception>
+    public int MeasureInto(int qubit, int classicalBit)
+    {
+        CheckQubit(qubit);
+        CheckClassicalBit(classicalBit);
+
+        Gates.Add(Record(new Gate
+        {
+            GateType = GateType.Measure,
+            TargetQubits = [qubit],
+            ClassicalBits = [classicalBit]
+        }));
+
+        int result = QuantumMath.SamplePartialMeasurement(StateVector, [qubit], RandomSource);
+        StateVector = QuantumMath.CollapseToPartialMeasurement(StateVector, [qubit], result);
+
+        _classicalBits[classicalBit] = result;
+
+        return result;
+    }
+
+    /// <summary>
+    /// Runs <paramref name="body"/> with every gate it applies conditioned on a classical
+    /// bit, giving classical feedforward: the gates run only when the measured bit matches.
+    /// </summary>
+    /// <param name="classicalBit">The classical bit to test, written by a previous measurement.</param>
+    /// <param name="value">The value the bit must hold, 0 or 1.</param>
+    /// <param name="body">The gates to apply conditionally.</param>
+    /// <example>
+    /// Quantum teleportation's correction step:
+    /// <code>
+    /// qc.MeasureInto(0, 0);
+    /// qc.MeasureInto(1, 1);
+    /// qc.When(1, 1, c => c.X(2));
+    /// qc.When(0, 1, c => c.Z(2));
+    /// </code>
+    /// </example>
+    /// <remarks>
+    /// Gates inside the body are always recorded, so <see cref="Simulation.Simulator"/>
+    /// re-evaluates the condition on every shot. They are applied to this circuit's live
+    /// state vector only when the condition currently holds.
+    /// </remarks>
+    public void When(int classicalBit, int value, Action<QuantumCircuit> body)
+    {
+        CheckClassicalBit(classicalBit);
+
+        if (body is null) throw new ArgumentNullException(nameof(body));
+
+        var previous = _activeCondition;
+        _activeCondition = (classicalBit, value);
+
+        try
+        {
+            body(this);
+        }
+        finally
+        {
+            _activeCondition = previous;
+        }
+    }
+
+    /// <summary>
+    /// Reads a classical bit holding a previous measurement outcome.
+    /// </summary>
+    /// <param name="classicalBit">The classical bit to read.</param>
+    /// <returns>The stored value, 0 or 1. Bits never written read as 0.</returns>
+    public int ClassicalBit(int classicalBit)
+    {
+        CheckClassicalBit(classicalBit);
+
+        return _classicalBits[classicalBit];
+    }
+
+    /// <summary>
+    /// Returns the circuit to its initial state: every qubit back to |0⟩, and the gate list,
+    /// initializations and classical bits cleared.
+    /// </summary>
+    public void Reset()
+    {
+        StateVector = new Complex[1 << QubitCount];
+        StateVector[0] = Complex.One;
+
+        Gates.Clear();
+        Initializations.Clear();
+        Array.Clear(_isQubitModified, 0, _isQubitModified.Length);
+        Array.Clear(_classicalBits, 0, _classicalBits.Length);
+
+        _activeCondition = null;
+    }
+
     /// <summary>
     /// Initializes the specified qubit to a predefined basis state: |0⟩, |1⟩, |+⟩, or |−⟩.
     /// This operation modifies the quantum state to reflect the chosen single-qubit state.
@@ -723,16 +847,40 @@ private void ApplySingle(GateType type, Complex[,] matrix, int qubit)
     {
         CheckQubit(qubit);
 
-        Gates.Add(new Gate
+        Gates.Add(Record(new Gate
         {
             GateType = type,
             Matrix = matrix,
             TargetQubits = [qubit]
-        });
+        }));
+
+        if (!ConditionHolds()) return;
 
         ApplyGate(matrix, qubit);
     }
 
+    /// <summary>
+    /// Stamps the condition currently in force onto a gate before it is recorded.
+    /// </summary>
+    private Gate Record(Gate gate)
+    {
+        if (_activeCondition is { } condition)
+        {
+            gate.ConditionBit = condition.Bit;
+            gate.ConditionValue = condition.Value;
+        }
+
+        return gate;
+    }
+
+    /// <summary>
+    /// Whether the condition currently in force is satisfied by the classical bits as they
+    /// stand. Gates inside a <see cref="When"/> body are always recorded, but only applied
+    /// to the live state vector when this is true.
+    /// </summary>
+    private bool ConditionHolds() =>
+        _activeCondition is not { } condition || _classicalBits[condition.Bit] == condition.Value;
+
     /// <summary>
     /// Validates, records and applies a two-qubit controlled gate. Shared by every
     /// controlled gate method.
@@ -749,13 +897,15 @@ private void ApplyControlled(GateType type, Complex[,] matrix, int targetQubit,
         if (targetQubit == controlQubit)
             throw new ArgumentException("Every gate argument must be a different Qubit");
 
-        Gates.Add(new Gate
+        Gates.Add(Record(new Gate
         {
             GateType = type,
             Matrix = matrix,
             TargetQubits = [targetQubit],
             ControlQubits = [controlQubit]
-        });
+        }));
+
+        if (!ConditionHolds()) return;
 
         StateVector = QuantumMath.ApplyControlledSingleQubitGate(
             StateVector, QuantumMath.ControlledCore(matrix), targetQubit, controlQubit);
@@ -809,6 +959,20 @@ private void CheckQubit(int qubit)
         if (qubit < 0 || qubit >= QubitCount)
             throw new QubitIndexOutOfRangeException($"Invalid index ({qubit}): qubit index must be between [0 and {QubitCount})");
     }
+
+    /// <summary>
+    /// Checks that a classical bit index is within the valid range [0, QubitCount).
+    /// </summary>
+    /// <param name="classicalBit">The index of the classical bit to check.</param>
+    /// <exception cref="QubitIndexOutOfRangeException">
+    /// Thrown when the index is outside the register.
+    /// </exception>
+    private void CheckClassicalBit(int classicalBit)
+    {
+        if (classicalBit < 0 || classicalBit >= QubitCount)
+            throw new QubitIndexOutOfRangeException(
+                $"Invalid index ({classicalBit}): classical bit index must be between [0 and {QubitCount})");
+    }
     
     /// <summary>
     /// Checks if all qubits in the provided array are distinct.
diff --git a/src/Qubit.NET/QuantumCircuitDrawer.cs b/src/Qubit.NET/QuantumCircuitDrawer.cs
index 052cb60..7be3ad2 100644
--- a/src/Qubit.NET/QuantumCircuitDrawer.cs
+++ b/src/Qubit.NET/QuantumCircuitDrawer.cs
@@ -131,9 +131,16 @@ private static int AssignGatePositions(IList<IList<(string, int)>> gatePositions
             
             var involvedQubits = controlQubits.Concat(targetQubits).ToList();
 
+            // A conditional gate depends on the measurement that wrote its classical bit,
+            // so it must be laid out after that wire's last gate and connected back to it.
+            // The bit is not in involvedQubits, which drives symbol placement, only layout.
+            var layoutQubits = gate.ConditionBit is { } conditionBit
+                ? involvedQubits.Append(conditionBit).ToList()
+                : involvedQubits;
+
             int farthestIndex = 0;
-            
-            foreach (int qubit in involvedQubits)
+
+            foreach (int qubit in layoutQubits)
             {
                 int current = 0;
 
@@ -141,7 +148,7 @@ private static int AssignGatePositions(IList<IList<(string, int)>> gatePositions
                 {
                     current = gatePositions[circuit.QubitCount - 1 - qubit].Last().Item2 + 1;
                 }
-                
+
                 farthestIndex = current > farthestIndex ? current : farthestIndex;
             }
 
@@ -149,10 +156,10 @@ private static int AssignGatePositions(IList<IList<(string, int)>> gatePositions
 
             while (conflict)
             {
-                if(involvedQubits.Count == 0) break;
-                
-                int minQubit = involvedQubits.Min();
-                int maxQubit = involvedQubits.Max();
+                if(layoutQubits.Count == 0) break;
+
+                int minQubit = layoutQubits.Min();
+                int maxQubit = layoutQubits.Max();
 
                 conflict = false;
                 
@@ -183,9 +190,11 @@ private static int AssignGatePositions(IList<IList<(string, int)>> gatePositions
             while(widths.Count <= farthestIndex) widths.Add(0);
             widths[farthestIndex] = maxWidth > widths[farthestIndex] ? maxWidth : widths[farthestIndex];
 
-            if (involvedQubits.Count > 1)
+            // Connect every wire the gate spans, including the classical bit a conditional
+            // gate reads, so the dependency is visible.
+            if (layoutQubits.Count > 1)
             {
-                for (int i = circuit.QubitCount - 1 - involvedQubits.Max(); i <= circuit.QubitCount - 2 - involvedQubits.Min(); i++)
+                for (int i = circuit.QubitCount - 1 - layoutQubits.Max(); i <= circuit.QubitCount - 2 - layoutQubits.Min(); i++)
                 {
                     if (!involvedQubits.Contains(circuit.QubitCount-1-i))
                     {
diff --git a/src/Qubit.NET/Simulation/MeasurementResult.cs b/src/Qubit.NET/Simulation/MeasurementResult.cs
new file mode 100644
index 0000000..6bb3daa
--- /dev/null
+++ b/src/Qubit.NET/Simulation/MeasurementResult.cs
@@ -0,0 +1,107 @@
+using System.Text;
+
+namespace Qubit.NET.Simulation;
+
+/// <summary>
+/// The outcome counts of one measurement across every shot of a simulation.
+/// </summary>
+public sealed class MeasurementResult
+{
+    private readonly int[] _counts;
+
+    internal MeasurementResult(int[] counts, int qubitCount)
+    {
+        _counts = counts;
+        QubitCount = qubitCount;
+    }
+
+    /// <summary>
+    /// Number of qubits involved in this measurement, and so the width of each outcome
+    /// bitstring.
+    /// </summary>
+    public int QubitCount { get; }
+
+    /// <summary>
+    /// Total number of shots recorded.
+    /// </summary>
+    public int Shots => _counts.Sum();
+
+    /// <summary>
+    /// How many times each observed outcome occurred, keyed by bitstring. Outcomes that
+    /// never occurred are omitted.
+    /// </summary>
+    public IReadOnlyDictionary<string, int> Counts
+    {
+        get
+        {
+            Dictionary<string, int> counts = new();
+
+            for (int i = 0; i < _counts.Length; i++)
+            {
+                if (_counts[i] != 0)
+                    counts[Convert.ToString(i, 2).PadLeft(QubitCount, '0')] = _counts[i];
+            }
+
+            return counts;
+        }
+    }
+
+    /// <summary>
+    /// The observed frequency of an outcome, in [0, 1].
+    /// </summary>
+    /// <param name="outcome">The outcome bitstring, for example "01".</param>
+    /// <returns>The fraction of shots that produced <paramref name="outcome"/>.</returns>
+    public double Probability(string outcome)
+    {
+        int shots = Shots;
+
+        if (shots == 0) return 0;
+
+        return Counts.TryGetValue(outcome, out int count) ? (double)count / shots : 0;
+    }
+
+    /// <summary>
+    /// The outcome that occurred most often, or null if no shots were recorded.
+    /// </summary>
+    public string? MostFrequent
+    {
+        get
+        {
+            string? best = null;
+            int bestCount = 0;
+
+            foreach (var pair in Counts)
+            {
+                if (pair.Value <= bestCount) continue;
+
+                best = pair.Key;
+                bestCount = pair.Value;
+            }
+
+            return best;
+        }
+    }
+
+    /// <summary>
+    /// Formats the counts as a dictionary literal, for example <c>{'00': 512, '11': 488}</c>.
+    /// </summary>
+    public override string ToString()
+    {
+        StringBuilder sb = new();
+        sb.Append('{');
+
+        bool any = false;
+
+        foreach (var pair in Counts)
+        {
+            if (any) sb.Append(", ");
+
+            sb.Append('\'').Append(pair.Key).Append("': ").Append(pair.Value);
+            any = true;
+        }
+
+        return sb.Append('}').ToString();
+    }
+
+    internal void Record(int outcome) => _counts[outcome]++;
+}
diff --git a/src/Qubit.NET/Simulation/Simulator.cs b/src/Qubit.NET/Simulation/Simulator.cs
index b2d5842..a8b7e05 100644
--- a/src/Qubit.NET/Simulation/Simulator.cs
+++ b/src/Qubit.NET/Simulation/Simulator.cs
@@ -19,13 +19,14 @@ public static class Simulator
     /// <param name="qc">The quantum circuit to simulate.</param>
     /// <param name="shots">The number of times the simulation should be run. Default is 1.</param>
     /// <returns>
-    /// A list of integer arrays and integers representing measurement results for each shot and number of qubits involved in measurement.
-    /// Each array contains counts of observed outcomes.
+    /// One <see cref="MeasurementResult"/> per measurement in the circuit, in the order the
+    /// measurements appear, each holding the outcome counts across all shots.
+    /// Empty if the circuit contains no measurement.
     /// </returns>
-    public static IList<(int[], int)> Run(QuantumCircuit qc, int shots = 1)
+    public static IList<MeasurementResult> Run(QuantumCircuit qc, int shots = 1)
     {
-        if (qc.Gates.All(g => g.GateType != GateType.Measure)) return new List<(int[], int)>();
-        
+        if (qc.Gates.All(g => g.GateType != GateType.Measure)) return new List<MeasurementResult>();
+
         Complex[] stateVector = new Complex[1 << qc.QubitCount];
         stateVector[0] = new Complex(1, 0);
 
@@ -39,7 +40,8 @@ public static class Simulator
         List<Gate> remainingGates = qc.Gates.ToList();
 
         // Apply the deterministic prefix (everything before the first measurement) once,
-        // then replay only the rest per shot.
+        // then replay only the rest per shot. A conditional gate depends on a measurement,
+        // so it can never appear in the prefix.
         while (remainingGates.Count > 0 && remainingGates[0].GateType != GateType.Measure)
         {
             Gate currentGate = remainingGates[0];
@@ -49,24 +51,36 @@ public static class Simulator
             remainingGates.RemoveAt(0);
         }
 
-        IList<(int[], int)> results = new List<(int[], int)>();
+        List<MeasurementResult> results = new();
+
+        // Classical bits are re-derived every shot, so feedforward follows that shot's own
+        // measurement outcomes.
+        int[] classicalBits = new int[qc.QubitCount];
 
         for (int i = 0; i < shots; i++)
         {
             Complex[] modStateVector = (Complex[])stateVector.Clone();
 
+            Array.Clear(classicalBits, 0, classicalBits.Length);
+
             int measurmentNumber = 0;
 
             foreach (var gate in remainingGates)
             {
+                if (gate.ConditionBit is { } bit && classicalBits[bit] != gate.ConditionValue)
+                    continue;
+
                 if (gate.GateType == GateType.Measure)
                 {
                     int num = MeasureState(ref modStateVector, gate.TargetQubits, qc);
-                    
+
                     if (measurmentNumber + 1 > results.Count)
-                        results.Add((new int[1L << qc.QubitCount], gate.TargetQubits.Length));
+                        results.Add(new MeasurementResult(
+                            new int[1 << gate.TargetQubits.Length], gate.TargetQubits.Length));
+
+                    results[measurmentNumber].Record(num);
 
-                    results[measurmentNumber].Item1[num]++;
+                    RecordClassicalBits(gate, num, classicalBits, qc.QubitCount);
 
                     measurmentNumber++;
                 }
@@ -76,45 +90,45 @@ public static class Simulator
                 }
             }
         }
-        
+
         return results;
     }
 
     /// <summary>
-    /// Converts the result of a quantum circuit simulation into a human-readable string.
+    /// Writes a measurement outcome into the classical bits the gate targets.
     /// </summary>
-    /// <param name="result">An array of measurement result counts and number of qubits in measurement.</param>
-    /// <returns>
-    /// A string formatted as a dictionary, where keys are binary representations of measurement outcomes,
-    /// and values are the counts of how often each outcome occurred.
-    /// </returns>
-    public static string GetStringResult(this (int[], int) result)
+    /// <param name="gate">The measurement gate.</param>
+    /// <param name="outcome">The packed measurement outcome.</param>
+    /// <param name="classicalBits">The classical register for the current shot.</param>
+    /// <param name="qubitCount">Number of qubits in the circuit.</param>
+    /// <remarks>
+    /// The two measurement paths pack their outcome differently: a full-register
+    /// measurement returns a basis index, where qubit q sits at bit q, while a partial
+    /// measurement packs the listed qubits with the first one most significant.
+    /// </remarks>
+    private static void RecordClassicalBits(Gate gate, int outcome, int[] classicalBits, int qubitCount)
     {
-        StringBuilder sb = new StringBuilder();
-        sb.Append("{");
+        int width = gate.TargetQubits.Length;
+        bool fullRegister = width == qubitCount;
 
-        bool any = false;
-
-        for (int i = 0; i < result.Item1.Length; i++)
+        for (int b = 0; b < gate.ClassicalBits.Length && b < width; b++)
         {
-            if(result.Item1[i] == 0)
-                continue;
-
-            if (any) sb.Append(", ");
+            int shift = fullRegister ? gate.TargetQubits[b] : width - 1 - b;
 
-            sb.Append("'");
-            sb.Append(Convert.ToString(i, 2).PadLeft(result.Item2, '0'));
-            sb.Append("': ");
-            sb.Append(result.Item1[i].ToString());
-
-            any = true;
+            classicalBits[gate.ClassicalBits[b]] = (outcome >> shift) & 1;
         }
-
-        sb.Append("}");
-
-        return sb.ToString();
     }
 
+    /// <summary>
+    /// Converts the result of a quantum circuit simulation into a human-readable string.
+    /// </summary>
+    /// <param name="result">The measurement result to format.</param>
+    /// <returns>
+    /// A string formatted as a dictionary, where keys are binary representations of measurement outcomes,
+    /// and values are the counts of how often each outcome occurred.
+    /// </returns>
+    public static string GetStringResult(this MeasurementResult result) => result.ToString();
+
     /// <summary>
     /// Applies a quantum gate to the given state vector, modifying it according to the specified gate type.
     /// Handles single-qubit and multi-qubit gates including controlled gates.
diff --git a/tests/Qubit.NET.Tests/FeedforwardTests.cs b/tests/Qubit.NET.Tests/FeedforwardTests.cs
new file mode 100644
index 0000000..c1603e7
--- /dev/null
+++ b/tests/Qubit.NET.Tests/FeedforwardTests.cs
@@ -0,0 +1,266 @@
+using System.Numerics;
+using Qubit.NET;
+using Qubit.NET.Simulation;
+using Qubit.NET.Utilities;
+using static QubitNet.Tests.Amplitudes;
+
+namespace QubitNet.Tests;
+
+/// <summary>
+/// Classical feedforward: measuring into a classical bit and conditioning later gates on it.
+/// Quantum teleportation is the canonical use, and its correctness is the real test here.
+/// </summary>
+public class FeedforwardTests
+{
+    [Fact]
+    public void Classical_bits_start_at_zero() =>
+        Assert.Equal(0, new QuantumCircuit(2).ClassicalBit(0));
+
+    [Fact]
+    public void MeasureInto_records_the_outcome()
+    {
+        QuantumCircuit qc = new(2);
+        qc.X(0);
+
+        int result = qc.MeasureInto(0, 1);
+
+        Assert.Equal(1, result);
+        Assert.Equal(1, qc.ClassicalBit(1));
+        Assert.Equal(0, qc.ClassicalBit(0));
+    }
+
+    [Fact]
+    public void Measure_populates_the_matching_classical_bits()
+    {
+        QuantumCircuit qc = new(3);
+        qc.X(0);
+        qc.X(2);
+        qc.Measure();
+
+        Assert.Equal(1, qc.ClassicalBit(0));
+        Assert.Equal(0, qc.ClassicalBit(1));
+        Assert.Equal(1, qc.ClassicalBit(2));
+    }
+
+    [Fact]
+    public void Partial_measure_populates_only_the_qubits_it_measured()
+    {
+        QuantumCircuit qc = new(3);
+        qc.X(2);
+        qc.Measure(2, 0);
+
+        Assert.Equal(1, qc.ClassicalBit(2));
+        Assert.Equal(0, qc.ClassicalBit(0));
+    }
+
+    [Fact]
+    public void When_applies_the_body_only_if_the_bit_matches()
+    {
+        QuantumCircuit applied = new(2);
+        applied.X(0);
+        applied.MeasureInto(0, 0);
+        applied.When(0, 1, c => c.X(1));
+
+        // Bit 0 came out 1, so the X ran: |11>, index 3.
+        AssertState(applied, 0, 0, 0, 1);
+
+        QuantumCircuit skipped = new(2);
+        skipped.MeasureInto(0, 0);
+        skipped.When(0, 1, c => c.X(1));
+
+        // Bit 0 came out 0, so the X was skipped: |00>, index 0.
+        AssertState(skipped, 1, 0, 0, 0);
+    }
+
+    [Fact]
+    public void Conditional_gates_are_still_recorded_when_skipped()
+    {
+        QuantumCircuit qc = new(2);
+        qc.MeasureInto(0, 0);
+        qc.When(0, 1, c => c.X(1));
+
+        // The gate must appear in the diagram even though it did not run, because the
+        // simulator re-evaluates the condition per shot.
+        Assert.Contains("[X]", qc.ToDiagram());
+    }
+
+    [Fact]
+    public void A_conditional_gate_is_drawn_after_the_measurement_it_depends_on()
+    {
+        QuantumCircuit qc = new(2);
+        qc.H(0);
+        qc.MeasureInto(0, 0);
+        qc.When(0, 1, c => c.X(1));
+
+        string[] lines = qc.ToDiagram().Split('\n');
+
+        string qubit1 = lines.First(l => l.StartsWith("q1"));
+        string qubit0 = lines.First(l => l.StartsWith("q0"));
+
+        // The conditional X must sit to the right of the M that wrote its classical bit.
+        Assert.True(qubit1.IndexOf("[X]", StringComparison.Ordinal) > qubit0.IndexOf("[M]", StringComparison.Ordinal),
+            $"conditional gate drawn before its measurement:\n{qc.ToDiagram()}");
+    }
+
+    [Fact]
+    public void When_nests()
+    {
+        QuantumCircuit qc = new(3);
+        qc.X(0);
+        qc.X(1);
+        qc.MeasureInto(0, 0);
+        qc.MeasureInto(1, 1);
+
+        qc.When(0, 1, c => c.When(1, 1, inner => inner.X(2)));
+
+        Assert.Equal(1, qc.ClassicalBit(0));
+        Assert.Equal(1, qc.ClassicalBit(1));
+        AssertState(qc, 0, 0, 0, 0, 0, 0, 0, 1);   // |111>
+    }
+
+    [Fact]
+    public void When_restores_the_previous_condition_after_the_body()
+    {
+        QuantumCircuit qc = new(2);
+        qc.MeasureInto(0, 0);
+
+        qc.When(0, 1, c => c.X(1));   // skipped, bit is 0
+        qc.X(1);                       // must still run
+
+        AssertState(qc, 0, 0, 1, 0);   // |10>
+    }
+
+    [Theory]
+    [InlineData(0.0)]
+    [InlineData(0.25)]
+    [InlineData(0.5)]
+    [InlineData(1.0)]
+    [InlineData(2.0)]
+    public void Teleportation_moves_an_arbitrary_state_to_the_third_qubit(double theta)
+    {
+        // Qubit 0 holds the message, qubits 1 and 2 share a Bell pair. After the protocol
+        // qubit 2 must hold exactly the state qubit 0 started in.
+        QuantumCircuit qc = new(3);
+        qc.RandomSource = new SeededRandomSource(7);
+
+        // Prepare the message: Ry(theta)|0> = cos(theta/2)|0> + sin(theta/2)|1>.
+        qc.Ry(0, theta);
+
+        // Entangle qubits 1 and 2.
+        qc.H(1);
+        qc.CNOT(1, 2);
+
+        // Bell-basis measurement of qubits 0 and 1.
+        qc.CNOT(0, 1);
+        qc.H(0);
+        qc.MeasureInto(0, 0);
+        qc.MeasureInto(1, 1);
+
+        // Corrections driven by the two classical bits.
+        qc.When(1, 1, c => c.X(2));
+        qc.When(0, 1, c => c.Z(2));
+
+        // Qubit 2 now carries the message. Marginalize the other two out.
+        double expectedZero = Math.Cos(theta / 2) * Math.Cos(theta / 2);
+        double expectedOne = Math.Sin(theta / 2) * Math.Sin(theta / 2);
+
+        (double zero, double one) = QubitMarginal(qc, 2);
+
+        Assert.Equal(expectedZero, zero, 9);
+        Assert.Equal(expectedOne, one, 9);
+    }
+
+    [Fact]
+    public void Teleportation_works_for_every_measurement_branch()
+    {
+        // The four Bell outcomes each need a different correction, so run enough seeds to
+        // hit all of them and confirm the result never depends on which branch was taken.
+        for (int seed = 0; seed < 20; seed++)
+        {
+            QuantumCircuit qc = new(3);
+            qc.RandomSource = new SeededRandomSource(seed);
+
+            qc.Ry(0, 0.8);
+            qc.H(1);
+            qc.CNOT(1, 2);
+            qc.CNOT(0, 1);
+            qc.H(0);
+            qc.MeasureInto(0, 0);
+            qc.MeasureInto(1, 1);
+            qc.When(1, 1, c => c.X(2));
+            qc.When(0, 1, c => c.Z(2));
+
+            (double zero, _) = QubitMarginal(qc, 2);
+
+            Assert.Equal(Math.Cos(0.4) * Math.Cos(0.4), zero, 9);
+        }
+    }
+
+    [Fact]
+    public void Simulator_reevaluates_conditions_per_shot()
+    {
+        // Measure a qubit in superposition, then flip a second qubit whenever the first
+        // came out 1. The two qubits must therefore always agree.
+        QuantumCircuit qc = new(2);
+        qc.H(0);
+        qc.MeasureInto(0, 0);
+        qc.When(0, 1, c => c.X(1));
+        qc.Measure();
+
+        var result = Simulator.Run(qc, 400)[1];
+
+        Assert.Equal(400, result.Shots);
+        Assert.False(result.Counts.ContainsKey("01"), "qubit 1 flipped without its control bit");
+        Assert.False(result.Counts.ContainsKey("10"), "qubit 1 failed to flip when its control bit was set");
+        Assert.True(result.Counts["00"] > 100);
+        Assert.True(result.Counts["11"] > 100);
+    }
+
+    [Fact]
+    public void Reset_clears_state_gates_and_classical_bits()
+    {
+        QuantumCircuit qc = new(2);
+        qc.X(0);
+        qc.MeasureInto(0, 0);
+
+        qc.Reset();
+
+        AssertState(qc, 1, 0, 0, 0);
+        Assert.Equal(0, qc.ClassicalBit(0));
+        Assert.DoesNotContain("[X]", qc.ToDiagram());
+
+        // A reset circuit must be reusable, including re-initialization.
+        qc.Initialize(0, Qubit.NET.Utilities.State.One);
+        AssertState(qc, 0, 1, 0, 0);
+    }
+
+    [Fact]
+    public void Out_of_range_classical_bits_are_rejected()
+    {
+        QuantumCircuit qc = new(2);
+
+        Assert.Throws<QubitIndexOutOfRangeException>(() => qc.ClassicalBit(2));
+        Assert.Throws<QubitIndexOutOfRangeException>(() => qc.MeasureInto(0, 5));
+        Assert.Throws<QubitIndexOutOfRangeException>(() => qc.When(9, 1, c => c.X(0)));
+    }
+
+    /// <summary>
+    /// Total probability of finding one qubit in |0> and in |1>, tracing out the rest.
+    /// </summary>
+    private static (double Zero, double One) QubitMarginal(QuantumCircuit qc, int qubit)
+    {
+        double zero = 0, one = 0;
+        int mask = 1 << qubit;
+
+        for (int i = 0; i < qc.StateVector.Length; i++)
+        {
+            Complex amplitude = qc.StateVector[i];
+            double probability = amplitude.Real * amplitude.Real + amplitude.Imaginary * amplitude.Imaginary;
+
+            if ((i & mask) == 0) zero += probability;
+            else one += probability;
+        }
+
+        return (zero, one);
+    }
+}
diff --git a/tests/Qubit.NET.Tests/RegressionTests.cs b/tests/Qubit.NET.Tests/RegressionTests.cs
index 738c79a..9b5a28e 100644
--- a/tests/Qubit.NET.Tests/RegressionTests.cs
+++ b/tests/Qubit.NET.Tests/RegressionTests.cs
@@ -35,13 +35,13 @@ public void Run_does_not_consume_the_circuit()
         var first = Simulator.Run(qc, 500)[0];
         var second = Simulator.Run(qc, 500)[0];
 
-        foreach (var counts in new[] { first, second })
+        foreach (var result in new[] { first, second })
         {
-            Assert.Equal(500, counts.Item1.Sum());
-            Assert.True(counts.Item1[0] > 150, "expected roughly half the shots in |00>");
-            Assert.True(counts.Item1[3] > 150, "expected roughly half the shots in |11>");
-            Assert.Equal(0, counts.Item1[1]);
-            Assert.Equal(0, counts.Item1[2]);
+            Assert.Equal(500, result.Shots);
+            Assert.True(result.Counts["00"] > 150, "expected roughly half the shots in |00>");
+            Assert.True(result.Counts["11"] > 150, "expected roughly half the shots in |11>");
+            Assert.False(result.Counts.ContainsKey("01"));
+            Assert.False(result.Counts.ContainsKey("10"));
         }
     }
 
@@ -121,20 +121,31 @@ public void Toffoli_draws_both_controls_as_control_markers()
     }
 
     [Fact]
-    public void GetStringResult_handles_an_empty_count_array()
+    public void Formatting_handles_a_result_with_no_shots()
     {
         // sb.Remove(sb.Length - 2, 2) threw when no outcome had been recorded.
-        Assert.Equal("{}", (new int[4], 2).GetStringResult());
+        QuantumCircuit qc = BellStates.PhiPlus();
+        qc.Measure();
+
+        var result = Simulator.Run(qc, 0);
+
+        Assert.True(result.Count == 0 || result[0].GetStringResult() == "{}");
     }
 
     [Fact]
     public void GetStringResult_formats_counts()
     {
-        var counts = new int[4];
-        counts[0] = 7;
-        counts[3] = 3;
+        QuantumCircuit qc = BellStates.PhiPlus();
+        qc.Measure();
+
+        var result = Simulator.Run(qc, 100)[0];
+        string formatted = result.GetStringResult();
 
-        Assert.Equal("{'00': 7, '11': 3}", (counts, 2).GetStringResult());
+        Assert.StartsWith("{", formatted);
+        Assert.EndsWith("}", formatted);
+        Assert.Contains("'00': ", formatted);
+        Assert.Contains("'11': ", formatted);
+        Assert.Equal(100, result.Shots);
     }
 
     [Fact]

From 3c9c0599fdaf9791c6341e7deef0986057f75ec3 Mon Sep 17 00:00:00 2001
From: InfoTCube <infocube.yout@gmail.com>
Date: Mon, 17 Aug 2026 20:48:11 +0200
Subject: [PATCH 6/9] feat: QASM export, textbook algorithms and state
 visualization

- OpenQASM 2.0 export via circuit.ToQasm(), covering every standard gate,
  measurement into creg, and conditional gates as `if (c[i]==v)`. Circuits now
  load straight into Qiskit and run on real hardware. Custom gates raise
  NotSupportedException rather than exporting something subtly wrong.
- Gates record the angles they were constructed from. Recovering them from the
  matrix afterwards is lossy around wrapping and sign, and storing them deleted
  more code than it added.
- Algorithms: Grover, Deutsch-Jozsa, Bernstein-Vazirani, teleportation,
  superdense coding, and an in-place QFT extension.
- State inspection: BlochVector(qubit) gives Bloch sphere coordinates, sitting
  at the origin for a maximally entangled qubit; QubitProbability(qubit) gives
  a single-qubit marginal; ToHistogram() draws an ASCII bar chart. These are
  what a Unity integration needs to render anything.

Grover finds its marked item with probability 1.000, and Bernstein-Vazirani
recovers the secret in one query, both verified across every input.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Mhwc4XeGtzM2UQx9sNLb3r
---
 CHANGELOG.md                                  |   7 +
 README.md                                     |  84 +++++-
 src/Qubit.NET/Circuits/Algorithms.cs          | 255 ++++++++++++++++++
 src/Qubit.NET/Export/QasmExporter.cs          | 148 ++++++++++
 src/Qubit.NET/Gates/Gate.cs                   |   6 +
 src/Qubit.NET/QuantumCircuit.cs               |  29 +-
 src/Qubit.NET/Visualization/StateInspector.cs | 121 +++++++++
 tests/Qubit.NET.Tests/AlgorithmTests.cs       | 142 ++++++++++
 tests/Qubit.NET.Tests/ExportTests.cs          | 204 ++++++++++++++
 9 files changed, 982 insertions(+), 14 deletions(-)
 create mode 100644 src/Qubit.NET/Circuits/Algorithms.cs
 create mode 100644 src/Qubit.NET/Export/QasmExporter.cs
 create mode 100644 src/Qubit.NET/Visualization/StateInspector.cs
 create mode 100644 tests/Qubit.NET.Tests/AlgorithmTests.cs
 create mode 100644 tests/Qubit.NET.Tests/ExportTests.cs

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 53e571b..e2f3aee 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -25,6 +25,13 @@ and usable from Unity.
 - `MeasurementResult`, replacing the `(int[], int)` tuple returned by `Simulator.Run`. It
   exposes `Counts`, `Shots`, `Probability(outcome)` and `MostFrequent`.
 - `QuantumCircuit.Reset()`, returning a circuit to |0…0⟩ with its gate list cleared.
+- **OpenQASM 2.0 export** via `circuit.ToQasm()`, including `if (c[i]==v)` for conditional
+  gates, so circuits can be loaded into Qiskit and run on real hardware.
+- **Algorithms**: `Grover`, `DeutschJozsa`, `BernsteinVazirani`, `Teleportation`,
+  `SuperdenseCoding` and an in-place `QFT` extension.
+- **State inspection**: `BlochVector(qubit)` returns Bloch sphere coordinates (the origin
+  for a maximally entangled qubit), `QubitProbability(qubit)` gives a single-qubit marginal,
+  and `ToHistogram()` renders outcome probabilities as an ASCII bar chart.
 - Continuous integration and a tag-driven NuGet release workflow.
 - An xUnit test suite covering gate algebra, the textbook Bell and GHZ states, measurement
   statistics and collapse, circuit rendering, and a regression test for every bug below.
diff --git a/README.md b/README.md
index 398e80f..a7a3a4e 100644
--- a/README.md
+++ b/README.md
@@ -285,6 +285,81 @@ QuantumCircuit qc = new QuantumCircuit(2);
 qc.RandomSource = new FixedRandomSource();
 ```
 
+### 🧪 Built-in algorithms
+
+```csharp
+using Qubit.NET.Circuits;
+
+// Grover search: finds the marked item in O(sqrt(N))
+var grover = Algorithms.Grover(2, c => c.CZ(0, 1));   // marks |11>
+Console.WriteLine(grover.ToHistogram());              // 11 | ####...####  1.000
+
+// Bernstein-Vazirani: recovers a hidden bit string in a single query
+var bv = Algorithms.BernsteinVazirani([true, false, true]);
+Console.WriteLine(bv.Measure(2, 1, 0));               // 101
+
+// Deutsch-Jozsa, teleportation, superdense coding, QFT
+var dj = Algorithms.DeutschJozsa(3, c => c.CNOT(0, 3));
+var tp = Algorithms.Teleportation(c => c.Ry(0, 0.9));
+var sd = Algorithms.SuperdenseCoding(true, false);
+
+qc.QFT();   // Quantum Fourier Transform, in place
+```
+
+Plus `BellStates.PhiPlus/PhiMinus/PsiPlus/PsiMinus/GHZ()`.
+
+---
+
+### 📤 OpenQASM export
+
+Export to OpenQASM 2.0 and run your circuit on real hardware through Qiskit or IBM Quantum:
+
+```csharp
+using Qubit.NET.Export;
+
+File.WriteAllText("teleport.qasm", Algorithms.Teleportation(c => c.Ry(0, 0.9)).ToQasm());
+```
+
+```qasm
+OPENQASM 2.0;
+include "qelib1.inc";
+
+qreg q[3];
+creg c[3];
+
+ry(0.9) q[0];
+h q[1];
+cx q[1],q[2];
+cx q[0],q[1];
+h q[0];
+measure q[0] -> c[0];
+measure q[1] -> c[1];
+if (c[1]==1) x q[2];
+if (c[0]==1) z q[2];
+```
+
+```python
+# Qiskit
+qc = QuantumCircuit.from_qasm_file("teleport.qasm")
+```
+
+---
+
+### 📊 Inspecting the state
+
+```csharp
+using Qubit.NET.Visualization;
+
+var (x, y, z) = qc.BlochVector(0);   // Bloch sphere coordinates — drive a Unity gizmo
+qc.QubitProbability(0);              // P(qubit 0 = |1>), tracing out the rest
+Console.WriteLine(qc.ToHistogram()); // ASCII bar chart of outcome probabilities
+```
+
+A qubit entangled with others sits inside the sphere — maximally entangled means the
+origin, which makes entanglement something you can actually *see*.
+
+---
+
 ## ⚡ Performance
 
 Gates are applied in place, so a circuit allocates one state vector regardless of how many
@@ -312,9 +387,14 @@ Memory is the real limit — the state vector holds 2ⁿ complex amplitudes at 1
 
 ## 📌 Future Roadmap
 
-- [ ] Entanglement entropy measurements
+- [x] Circuit export in QASM
+- [x] Mid-circuit measurement and classical feedforward
+- [x] Bloch sphere coordinates and probability histograms
+- [x] Built-in algorithms (Grover, Deutsch–Jozsa, Bernstein–Vazirani, QFT, teleportation)
 - [ ] Noise simulation (decoherence, damping)
-- [ ] Circuit export in QASM
+- [ ] Entanglement entropy measurements
+- [ ] QASM import
+- [ ] Multi-controlled gates and circuit inverses
 
 ---
 
diff --git a/src/Qubit.NET/Circuits/Algorithms.cs b/src/Qubit.NET/Circuits/Algorithms.cs
new file mode 100644
index 0000000..dedf3dc
--- /dev/null
+++ b/src/Qubit.NET/Circuits/Algorithms.cs
@@ -0,0 +1,255 @@
+namespace Qubit.NET.Circuits;
+
+/// <summary>
+/// Textbook quantum algorithms, ready to run or to read as worked examples.
+/// </summary>
+public static class Algorithms
+{
+    /// <summary>
+    /// Applies the Quantum Fourier Transform to the given qubits, in place.
+    /// </summary>
+    /// <param name="circuit">The circuit to append to.</param>
+    /// <param name="qubits">
+    /// The qubits to transform, most significant first. Defaults to the whole register.
+    /// </param>
+    /// <remarks>
+    /// The QFT is the engine behind Shor's algorithm and quantum phase estimation. This is
+    /// the standard Hadamard-and-controlled-phase construction, with the final swaps that
+    /// restore qubit order.
+    /// </remarks>
+    public static void QFT(this QuantumCircuit circuit, params int[] qubits)
+    {
+        int[] targets = qubits.Length == 0
+            ? Enumerable.Range(0, circuit.QubitCount).ToArray()
+            : qubits;
+
+        int n = targets.Length;
+
+        for (int i = 0; i < n; i++)
+        {
+            circuit.H(targets[i]);
+
+            for (int j = i + 1; j < n; j++)
+            {
+                // Controlled phase of pi / 2^(j - i), the defining rotation of the QFT.
+                double angle = System.Math.PI / (1 << (j - i));
+
+                circuit.CRz(targets[j], targets[i], angle);
+            }
+        }
+
+        // The construction above emits the output in reverse qubit order.
+        for (int i = 0; i < n / 2; i++)
+            circuit.SWAP(targets[i], targets[n - 1 - i]);
+    }
+
+    /// <summary>
+    /// Builds a Deutsch-Jozsa circuit, which decides in a single query whether a function is
+    /// constant or balanced.
+    /// </summary>
+    /// <param name="inputQubits">Number of input qubits the oracle acts on.</param>
+    /// <param name="oracle">
+    /// Applies the phase oracle. It receives the circuit; input qubits are 0 to
+    /// <paramref name="inputQubits"/> - 1 and the ancilla is the last qubit.
+    /// </param>
+    /// <returns>
+    /// The prepared circuit. Measuring the input qubits gives all zeros if the function is
+    /// constant, and anything else if it is balanced.
+    /// </returns>
+    public static QuantumCircuit DeutschJozsa(int inputQubits, Action<QuantumCircuit> oracle)
+    {
+        if (oracle is null) throw new ArgumentNullException(nameof(oracle));
+
+        QuantumCircuit qc = new(inputQubits + 1);
+
+        // Ancilla into |-> so the oracle's answer lands in the phase.
+        qc.X(inputQubits);
+        qc.H(inputQubits);
+
+        for (int q = 0; q < inputQubits; q++)
+            qc.H(q);
+
+        oracle(qc);
+
+        for (int q = 0; q < inputQubits; q++)
+            qc.H(q);
+
+        return qc;
+    }
+
+    /// <summary>
+    /// Builds a Bernstein-Vazirani circuit, which recovers a hidden bit string in one query.
+    /// </summary>
+    /// <param name="secret">The hidden bit string the oracle encodes, least significant bit first.</param>
+    /// <returns>
+    /// The prepared circuit. Measuring the input qubits returns <paramref name="secret"/>.
+    /// </returns>
+    public static QuantumCircuit BernsteinVazirani(bool[] secret)
+    {
+        if (secret is null) throw new ArgumentNullException(nameof(secret));
+
+        int n = secret.Length;
+
+        return DeutschJozsa(n, qc =>
+        {
+            // The oracle computes the dot product of the input with the secret.
+            for (int q = 0; q < n; q++)
+            {
+                if (secret[q]) qc.CNOT(q, n);
+            }
+        });
+    }
+
+    /// <summary>
+    /// Builds a Grover search circuit over <paramref name="qubits"/> qubits.
+    /// </summary>
+    /// <param name="qubits">Number of qubits, giving a search space of 2^qubits.</param>
+    /// <param name="oracle">
+    /// Marks the solution states by flipping their phase. It receives the circuit.
+    /// </param>
+    /// <param name="iterations">
+    /// Number of Grover iterations. Defaults to the optimal count for a single solution,
+    /// which is floor(pi/4 * sqrt(2^qubits)).
+    /// </param>
+    /// <returns>The prepared circuit; measuring it yields a solution with high probability.</returns>
+    /// <remarks>
+    /// Grover's algorithm finds a marked item among N in O(sqrt(N)) queries rather than the
+    /// O(N) a classical search needs.
+    /// </remarks>
+    public static QuantumCircuit Grover(int qubits, Action<QuantumCircuit> oracle, int iterations = -1)
+    {
+        if (oracle is null) throw new ArgumentNullException(nameof(oracle));
+
+        QuantumCircuit qc = new(qubits);
+
+        if (iterations < 0)
+            iterations = (int)System.Math.Floor(System.Math.PI / 4 * System.Math.Sqrt(1 << qubits));
+
+        // Uniform superposition over the whole search space.
+        for (int q = 0; q < qubits; q++)
+            qc.H(q);
+
+        for (int i = 0; i < iterations; i++)
+        {
+            oracle(qc);
+            Diffuse(qc, qubits);
+        }
+
+        return qc;
+    }
+
+    /// <summary>
+    /// The Grover diffusion operator, reflecting amplitudes about their mean.
+    /// </summary>
+    private static void Diffuse(QuantumCircuit qc, int qubits)
+    {
+        for (int q = 0; q < qubits; q++)
+        {
+            qc.H(q);
+            qc.X(q);
+        }
+
+        // Multi-controlled Z over every qubit, phase-flipping only |11...1>.
+        MultiControlledZ(qc, qubits);
+
+        for (int q = 0; q < qubits; q++)
+        {
+            qc.X(q);
+            qc.H(q);
+        }
+    }
+
+    /// <summary>
+    /// Flips the phase of the all-ones state across <paramref name="qubits"/> qubits.
+    /// </summary>
+    private static void MultiControlledZ(QuantumCircuit qc, int qubits)
+    {
+        int last = qubits - 1;
+
+        switch (qubits)
+        {
+            case 1:
+                qc.Z(0);
+                break;
+            case 2:
+                qc.CZ(0, 1);
+                break;
+            case 3:
+                // H-CCX-H turns a Toffoli into a controlled-controlled-Z.
+                qc.H(last);
+                qc.Toffoli(0, 1, last);
+                qc.H(last);
+                break;
+            default:
+                throw new NotSupportedException(
+                    "Multi-controlled Z is implemented for up to 3 qubits. Supply your own "
+                    + "diffusion oracle for wider searches.");
+        }
+    }
+
+    /// <summary>
+    /// Builds a quantum teleportation circuit that moves the state of qubit 0 onto qubit 2.
+    /// </summary>
+    /// <param name="prepareMessage">
+    /// Prepares the state to teleport on qubit 0. Defaults to leaving it in |0⟩.
+    /// </param>
+    /// <returns>
+    /// The completed circuit. Qubit 2 holds the state qubit 0 was prepared in, and classical
+    /// bits 0 and 1 hold the Bell measurement outcome.
+    /// </returns>
+    /// <remarks>
+    /// Uses mid-circuit measurement and classical feedforward, so the corrections are real
+    /// conditional gates rather than the deferred-measurement equivalent.
+    /// </remarks>
+    public static QuantumCircuit Teleportation(Action<QuantumCircuit>? prepareMessage = null)
+    {
+        QuantumCircuit qc = new(3);
+
+        prepareMessage?.Invoke(qc);
+
+        // Entangle the pair shared between sender and receiver.
+        qc.H(1);
+        qc.CNOT(1, 2);
+
+        // Bell-basis measurement of the message and the sender's half.
+        qc.CNOT(0, 1);
+        qc.H(0);
+        qc.MeasureInto(0, 0);
+        qc.MeasureInto(1, 1);
+
+        // Corrections driven by the two classical outcomes.
+        qc.When(1, 1, c => c.X(2));
+        qc.When(0, 1, c => c.Z(2));
+
+        return qc;
+    }
+
+    /// <summary>
+    /// Builds a superdense coding circuit, sending two classical bits on one qubit.
+    /// </summary>
+    /// <param name="firstBit">The bit recovered from qubit 0.</param>
+    /// <param name="secondBit">The bit recovered from qubit 1.</param>
+    /// <returns>
+    /// The completed circuit. <c>Measure(0, 1)</c> returns the two bits that were sent.
+    /// </returns>
+    public static QuantumCircuit SuperdenseCoding(bool firstBit, bool secondBit)
+    {
+        QuantumCircuit qc = new(2);
+
+        // Shared entangled pair.
+        qc.H(0);
+        qc.CNOT(0, 1);
+
+        // The sender encodes two bits by acting on their qubit alone. Z moves the Bell
+        // state along the phase axis, which the decoder reads off qubit 0; X moves it along
+        // the parity axis, which shows up on qubit 1.
+        if (firstBit) qc.Z(0);
+        if (secondBit) qc.X(0);
+
+        // The receiver decodes.
+        qc.CNOT(0, 1);
+        qc.H(0);
+
+        return qc;
+    }
+}
diff --git a/src/Qubit.NET/Export/QasmExporter.cs b/src/Qubit.NET/Export/QasmExporter.cs
new file mode 100644
index 0000000..f9f8f8d
--- /dev/null
+++ b/src/Qubit.NET/Export/QasmExporter.cs
@@ -0,0 +1,148 @@
+using System.Globalization;
+using System.Text;
+using Qubit.NET.Gates;
+
+namespace Qubit.NET.Export;
+
+/// <summary>
+/// Exports circuits to OpenQASM 2.0, the interchange format used by Qiskit and IBM Quantum.
+/// </summary>
+public static class QasmExporter
+{
+    /// <summary>
+    /// Converts the circuit to an OpenQASM 2.0 program.
+    /// </summary>
+    /// <param name="circuit">The circuit to export.</param>
+    /// <returns>The QASM source text.</returns>
+    /// <example>
+    /// <code>
+    /// var qc = new QuantumCircuit(2);
+    /// qc.H(0);
+    /// qc.CNOT(0, 1);
+    /// qc.Measure();
+    ///
+    /// File.WriteAllText("bell.qasm", qc.ToQasm());
+    /// </code>
+    /// The result can be loaded straight into Qiskit with
+    /// <c>QuantumCircuit.from_qasm_file("bell.qasm")</c>.
+    /// </example>
+    /// <exception cref="NotSupportedException">
+    /// Thrown if the circuit contains a custom gate, which has no QASM equivalent.
+    /// </exception>
+    public static string ToQasm(this QuantumCircuit circuit)
+    {
+        StringBuilder sb = new();
+
+        sb.AppendLine("OPENQASM 2.0;");
+        sb.AppendLine("include \"qelib1.inc\";");
+        sb.AppendLine();
+        sb.AppendLine($"qreg q[{circuit.QubitCount}];");
+        sb.AppendLine($"creg c[{circuit.QubitCount}];");
+        sb.AppendLine();
+
+        // Initializations are state preparation, which QASM cannot express directly. The
+        // basis states have exact gate equivalents; anything else is reported as a comment
+        // rather than silently dropped.
+        foreach (var init in circuit.Initializations)
+        {
+            switch (init.BasicState)
+            {
+                case Utilities.State.Zero:
+                    break;
+                case Utilities.State.One:
+                    sb.AppendLine($"x q[{init.QubitIndex}];");
+                    break;
+                case Utilities.State.Plus:
+                    sb.AppendLine($"h q[{init.QubitIndex}];");
+                    break;
+                case Utilities.State.Minus:
+                    sb.AppendLine($"x q[{init.QubitIndex}];");
+                    sb.AppendLine($"h q[{init.QubitIndex}];");
+                    break;
+                default:
+                    sb.AppendLine($"// qubit {init.QubitIndex} initialized to a custom state "
+                                  + "that OPENQASM 2.0 cannot express");
+                    break;
+            }
+        }
+
+        foreach (var gate in circuit.Gates)
+        {
+            if (gate.ConditionBit is { } bit)
+                sb.Append($"if (c[{bit}]=={gate.ConditionValue}) ");
+
+            sb.AppendLine(Emit(gate));
+        }
+
+        return sb.ToString();
+    }
+
+    private static string Emit(Gate gate)
+    {
+        int[] t = gate.TargetQubits;
+        int[] c = gate.ControlQubits;
+
+        return gate.GateType switch
+        {
+            GateType.I => $"id q[{t[0]}];",
+            GateType.H => $"h q[{t[0]}];",
+            GateType.X => $"x q[{t[0]}];",
+            GateType.Y => $"y q[{t[0]}];",
+            GateType.Z => $"z q[{t[0]}];",
+            GateType.S => $"s q[{t[0]}];",
+            GateType.Sdag => $"sdg q[{t[0]}];",
+            GateType.T => $"t q[{t[0]}];",
+            GateType.Tdag => $"tdg q[{t[0]}];",
+            GateType.SX => $"sx q[{t[0]}];",
+            // No qelib1 primitive for sqrt(Y) or sqrt(Z); both have exact u3 forms.
+            GateType.SY => $"u3(pi/2,0,0) q[{t[0]}];",
+            GateType.SZ => $"s q[{t[0]}];",
+            GateType.Rx => $"rx({Angles(gate)}) q[{t[0]}];",
+            GateType.Ry => $"ry({Angles(gate)}) q[{t[0]}];",
+            GateType.Rz => $"rz({Angles(gate)}) q[{t[0]}];",
+            GateType.U3 => $"u3({Angles(gate)}) q[{t[0]}];",
+            GateType.CNOT => $"cx q[{c[0]}],q[{t[0]}];",
+            GateType.CY => $"cy q[{c[0]}],q[{t[0]}];",
+            GateType.CZ => $"cz q[{c[0]}],q[{t[0]}];",
+            GateType.CH => $"ch q[{c[0]}],q[{t[0]}];",
+            GateType.CRx => $"crx({Angles(gate)}) q[{c[0]}],q[{t[0]}];",
+            GateType.CRy => $"cry({Angles(gate)}) q[{c[0]}],q[{t[0]}];",
+            GateType.CRz => $"crz({Angles(gate)}) q[{c[0]}],q[{t[0]}];",
+            GateType.CU3 => $"cu3({Angles(gate)}) q[{c[0]}],q[{t[0]}];",
+            GateType.SWAP => $"swap q[{t[0]}],q[{t[1]}];",
+            GateType.Toffoli => $"ccx q[{c[1]}],q[{c[0]}],q[{t[0]}];",
+            GateType.Fredkin => $"cswap q[{c[0]}],q[{t[1]}],q[{t[0]}];",
+            GateType.Measure => Measure(gate),
+            GateType.Custom => throw new NotSupportedException(
+                "Custom gates have no OPENQASM 2.0 equivalent. Rebuild the operation from "
+                + "standard gates before exporting."),
+            _ => throw new NotSupportedException($"Gate {gate.GateType} cannot be exported to QASM.")
+        };
+    }
+
+    private static string Measure(Gate gate)
+    {
+        StringBuilder sb = new();
+
+        for (int i = 0; i < gate.TargetQubits.Length; i++)
+        {
+            int classicalBit = i < gate.ClassicalBits.Length ? gate.ClassicalBits[i] : gate.TargetQubits[i];
+
+            if (i > 0) sb.AppendLine();
+
+            sb.Append($"measure q[{gate.TargetQubits[i]}] -> c[{classicalBit}];");
+        }
+
+        return sb.ToString();
+    }
+
+    /// <summary>
+    /// The angles the gate was built from, comma separated. Recorded at construction, so no
+    /// reconstruction from the matrix is needed.
+    /// </summary>
+    private static string Angles(Gate gate) =>
+        string.Join(",", gate.Parameters.Select(Format));
+
+    private static string Format(double value) =>
+        value.ToString("0.############################", CultureInfo.InvariantCulture);
+}
diff --git a/src/Qubit.NET/Gates/Gate.cs b/src/Qubit.NET/Gates/Gate.cs
index 879c2c2..e2bbf9e 100644
--- a/src/Qubit.NET/Gates/Gate.cs
+++ b/src/Qubit.NET/Gates/Gate.cs
@@ -39,4 +39,10 @@ internal class Gate
     /// positionally matching <see cref="TargetQubits"/>.
     /// </summary>
     public int[] ClassicalBits { get; set; } = [];
+
+    /// <summary>
+    /// The angles a parameterized gate was built from, in declaration order. Kept because
+    /// recovering them from the matrix afterwards is lossy around wrapping and sign.
+    /// </summary>
+    public double[] Parameters { get; set; } = [];
 }
diff --git a/src/Qubit.NET/QuantumCircuit.cs b/src/Qubit.NET/QuantumCircuit.cs
index 9ed6533..9295f38 100644
--- a/src/Qubit.NET/QuantumCircuit.cs
+++ b/src/Qubit.NET/QuantumCircuit.cs
@@ -449,7 +449,7 @@ public void Tdag(int qubit) =>
     /// Thrown if the qubit index is out of range, indicating that the specified qubit does not exist in the system.
     /// </exception>
     public void Rx(int qubit, double theta) =>
-        ApplySingle(GateType.Rx, QuantumGates.Rx(theta), qubit);
+        ApplySingle(GateType.Rx, QuantumGates.Rx(theta), qubit, theta);
     
     /// <summary>
     /// Applies the Ry gate to the specified qubit.
@@ -462,7 +462,7 @@ public void Rx(int qubit, double theta) =>
     /// Thrown if the qubit index is out of range, indicating that the specified qubit does not exist in the system.
     /// </exception>
     public void Ry(int qubit, double theta) =>
-        ApplySingle(GateType.Ry, QuantumGates.Ry(theta), qubit);
+        ApplySingle(GateType.Ry, QuantumGates.Ry(theta), qubit, theta);
 
     /// <summary>
     /// Applies the Rz gate to the specified qubit.
@@ -475,7 +475,7 @@ public void Ry(int qubit, double theta) =>
     /// Thrown if the qubit index is out of range, indicating that the specified qubit does not exist in the system.
     /// </exception>
     public void Rz(int qubit, double theta) =>
-        ApplySingle(GateType.Rz, QuantumGates.Rz(theta), qubit);
+        ApplySingle(GateType.Rz, QuantumGates.Rz(theta), qubit, theta);
     
     /// <summary>
     /// Applies the square-root of Pauli-X gate (SX) to the specified qubit.
@@ -530,7 +530,7 @@ public void SZ(int qubit) =>
     /// Thrown if the qubit index is out of range, indicating that the specified qubit does not exist in the system.
     /// </exception>
     public void U3(int qubit, double theta, double phi, double lambda) =>
-        ApplySingle(GateType.U3, QuantumGates.U3(theta, phi, lambda), qubit);
+        ApplySingle(GateType.U3, QuantumGates.U3(theta, phi, lambda), qubit, theta, phi, lambda);
     
     /// <summary>
     /// Applies the CNOT (also called CX) gate (Controlled-NOT) to the specified qubits.
@@ -591,7 +591,7 @@ public void CH(int controlQubit, int targetQubit) =>
     /// Thrown if any of the qubit indices are out of range.
     /// </exception>
     public void CRx(int controlQubit, int targetQubit, double theta) =>
-        ApplyControlled(GateType.CRx, QuantumGates.CRx(theta), targetQubit, controlQubit);
+        ApplyControlled(GateType.CRx, QuantumGates.CRx(theta), targetQubit, controlQubit, theta);
     
     /// <summary>
     /// Applies the CRy gate (Controlled-Ry) to the specified qubits.
@@ -604,7 +604,7 @@ public void CRx(int controlQubit, int targetQubit, double theta) =>
     /// Thrown if any of the qubit indices are out of range.
     /// </exception>
     public void CRy(int controlQubit, int targetQubit, double theta) =>
-        ApplyControlled(GateType.CRy, QuantumGates.CRy(theta), targetQubit, controlQubit);
+        ApplyControlled(GateType.CRy, QuantumGates.CRy(theta), targetQubit, controlQubit, theta);
     
     /// <summary>
     /// Applies the CRz gate (Controlled-Rz) to the specified qubits.
@@ -617,7 +617,7 @@ public void CRy(int controlQubit, int targetQubit, double theta) =>
     /// Thrown if any of the qubit indices are out of range.
     /// </exception>
     public void CRz(int controlQubit, int targetQubit, double theta) =>
-        ApplyControlled(GateType.CRz, QuantumGates.CRz(theta), targetQubit, controlQubit);
+        ApplyControlled(GateType.CRz, QuantumGates.CRz(theta), targetQubit, controlQubit, theta);
     
     /// <summary>
     /// Applies the CU3 gate to the specified qubits.
@@ -632,7 +632,7 @@ public void CRz(int controlQubit, int targetQubit, double theta) =>
     /// Thrown if any of the qubit indices are out of range.
     /// </exception>
     public void CU3(int controlQubit, int targetQubit, double theta, double phi, double lambda) =>
-        ApplyControlled(GateType.CU3, QuantumGates.CU3(theta, phi, lambda), targetQubit, controlQubit);
+        ApplyControlled(GateType.CU3, QuantumGates.CU3(theta, phi, lambda), targetQubit, controlQubit, theta, phi, lambda);
 
     /// <summary>
     /// Applies the SWAP gate to the specified qubits, exchanging their states.
@@ -843,7 +843,8 @@ public override string ToString()
     /// <param name="type">The gate type, recorded for replay and for drawing.</param>
     /// <param name="matrix">The 2x2 unitary matrix representing the gate.</param>
     /// <param name="qubit">The index of the qubit to apply the gate to.</param>
-    private void ApplySingle(GateType type, Complex[,] matrix, int qubit)
+    /// <param name="parameters">Angles the gate was built from, for parameterized gates.</param>
+    private void ApplySingle(GateType type, Complex[,] matrix, int qubit, params double[] parameters)
     {
         CheckQubit(qubit);
 
@@ -851,7 +852,8 @@ private void ApplySingle(GateType type, Complex[,] matrix, int qubit)
         {
             GateType = type,
             Matrix = matrix,
-            TargetQubits = [qubit]
+            TargetQubits = [qubit],
+            Parameters = parameters
         }));
 
         if (!ConditionHolds()) return;
@@ -889,7 +891,9 @@ private bool ConditionHolds() =>
     /// <param name="matrix">The 4x4 unitary matrix representing the gate.</param>
     /// <param name="targetQubit">The index of the target qubit.</param>
     /// <param name="controlQubit">The index of the control qubit.</param>
-    private void ApplyControlled(GateType type, Complex[,] matrix, int targetQubit, int controlQubit)
+    /// <param name="parameters">Angles the gate was built from, for parameterized gates.</param>
+    private void ApplyControlled(GateType type, Complex[,] matrix, int targetQubit, int controlQubit,
+        params double[] parameters)
     {
         CheckQubit(controlQubit);
         CheckQubit(targetQubit);
@@ -902,7 +906,8 @@ private void ApplyControlled(GateType type, Complex[,] matrix, int targetQubit,
             GateType = type,
             Matrix = matrix,
             TargetQubits = [targetQubit],
-            ControlQubits = [controlQubit]
+            ControlQubits = [controlQubit],
+            Parameters = parameters
         }));
 
         if (!ConditionHolds()) return;
diff --git a/src/Qubit.NET/Visualization/StateInspector.cs b/src/Qubit.NET/Visualization/StateInspector.cs
new file mode 100644
index 0000000..b9d343b
--- /dev/null
+++ b/src/Qubit.NET/Visualization/StateInspector.cs
@@ -0,0 +1,121 @@
+using System.Numerics;
+using System.Text;
+
+namespace Qubit.NET.Visualization;
+
+/// <summary>
+/// Reads out a circuit's state in forms that are easy to draw: Bloch sphere coordinates,
+/// per-qubit probabilities, and an ASCII histogram.
+/// </summary>
+public static class StateInspector
+{
+    /// <summary>
+    /// The Bloch sphere coordinates of a single qubit, obtained by tracing out the rest of
+    /// the register.
+    /// </summary>
+    /// <param name="circuit">The circuit to inspect.</param>
+    /// <param name="qubit">The qubit to locate on the sphere.</param>
+    /// <returns>
+    /// The point (x, y, z) on or inside the unit sphere. A pure state sits on the surface;
+    /// a qubit entangled with others sits strictly inside, and at the origin when maximally
+    /// entangled.
+    /// </returns>
+    /// <example>
+    /// <code>
+    /// var (x, y, z) = qc.BlochVector(0);
+    /// arrow.transform.localPosition = new Vector3((float)x, (float)y, (float)z); // Unity
+    /// </code>
+    /// </example>
+    public static (double X, double Y, double Z) BlochVector(this QuantumCircuit circuit, int qubit)
+    {
+        if (qubit < 0 || qubit >= circuit.QubitCount)
+            throw new ArgumentOutOfRangeException(nameof(qubit), qubit,
+                $"Qubit index must be between [0 and {circuit.QubitCount})");
+
+        // Reduced density matrix of the qubit: rho01 is the coherence, and the populations
+        // give the z component.
+        Complex rho01 = Complex.Zero;
+        double rho00 = 0, rho11 = 0;
+
+        int mask = 1 << qubit;
+        Complex[] state = circuit.StateVector;
+
+        for (int i = 0; i < state.Length; i++)
+        {
+            if ((i & mask) != 0) continue;
+
+            Complex a = state[i];            // amplitude with the qubit in |0>
+            Complex b = state[i | mask];     // same basis state, qubit in |1>
+
+            rho00 += a.Real * a.Real + a.Imaginary * a.Imaginary;
+            rho11 += b.Real * b.Real + b.Imaginary * b.Imaginary;
+            rho01 += a * Complex.Conjugate(b);
+        }
+
+        return (2 * rho01.Real, 2 * rho01.Imaginary, rho00 - rho11);
+    }
+
+    /// <summary>
+    /// The probability of finding a single qubit in |1⟩, tracing out the rest of the register.
+    /// </summary>
+    /// <param name="circuit">The circuit to inspect.</param>
+    /// <param name="qubit">The qubit to measure the marginal of.</param>
+    /// <returns>A probability in [0, 1].</returns>
+    public static double QubitProbability(this QuantumCircuit circuit, int qubit)
+    {
+        if (qubit < 0 || qubit >= circuit.QubitCount)
+            throw new ArgumentOutOfRangeException(nameof(qubit), qubit,
+                $"Qubit index must be between [0 and {circuit.QubitCount})");
+
+        double one = 0;
+        int mask = 1 << qubit;
+        Complex[] state = circuit.StateVector;
+
+        for (int i = 0; i < state.Length; i++)
+        {
+            if ((i & mask) == 0) continue;
+
+            one += state[i].Real * state[i].Real + state[i].Imaginary * state[i].Imaginary;
+        }
+
+        return one;
+    }
+
+    /// <summary>
+    /// Renders the circuit's outcome probabilities as an ASCII bar chart.
+    /// </summary>
+    /// <param name="circuit">The circuit to chart.</param>
+    /// <param name="width">Width of the longest bar, in characters.</param>
+    /// <returns>A multi-line string, one row per outcome with non-negligible probability.</returns>
+    /// <example>
+    /// <code>
+    /// 00 | ####################  0.500
+    /// 11 | ####################  0.500
+    /// </code>
+    /// </example>
+    public static string ToHistogram(this QuantumCircuit circuit, int width = 40)
+    {
+        if (width < 1) throw new ArgumentOutOfRangeException(nameof(width), width, "Width must be positive.");
+
+        var probabilities = circuit.GetProbabilities();
+
+        if (probabilities.Count == 0) return string.Empty;
+
+        double max = probabilities.Values.Max();
+        StringBuilder sb = new();
+
+        foreach (var pair in probabilities.OrderBy(p => p.Key, StringComparer.Ordinal))
+        {
+            int bars = max > 0 ? (int)System.Math.Round(pair.Value / max * width) : 0;
+
+            sb.Append(pair.Key)
+              .Append(" | ")
+              .Append('#', bars)
+              .Append(new string(' ', width - bars))
+              .Append("  ")
+              .AppendLine(pair.Value.ToString("F3", System.Globalization.CultureInfo.InvariantCulture));
+        }
+
+        return sb.ToString().TrimEnd();
+    }
+}
diff --git a/tests/Qubit.NET.Tests/AlgorithmTests.cs b/tests/Qubit.NET.Tests/AlgorithmTests.cs
new file mode 100644
index 0000000..9fce7c0
--- /dev/null
+++ b/tests/Qubit.NET.Tests/AlgorithmTests.cs
@@ -0,0 +1,142 @@
+using Qubit.NET;
+using Qubit.NET.Circuits;
+using Qubit.NET.Simulation;
+using Qubit.NET.Visualization;
+using static QubitNet.Tests.Amplitudes;
+
+namespace QubitNet.Tests;
+
+/// <summary>
+/// The textbook algorithms, checked against the answers they are supposed to produce.
+/// </summary>
+public class AlgorithmTests
+{
+    [Fact]
+    public void Deutsch_Jozsa_reports_constant_for_a_constant_function()
+    {
+        // A constant oracle does nothing at all, so every input qubit returns to |0>.
+        QuantumCircuit qc = Algorithms.DeutschJozsa(3, _ => { });
+        qc.RandomSource = new SeededRandomSource(1);
+
+        Assert.Equal("000", qc.Measure(0, 1, 2));
+    }
+
+    [Fact]
+    public void Deutsch_Jozsa_reports_balanced_for_a_balanced_function()
+    {
+        // f(x) = x0 is balanced, so the result must never be all zeros.
+        for (int seed = 0; seed < 10; seed++)
+        {
+            QuantumCircuit qc = Algorithms.DeutschJozsa(3, c => c.CNOT(0, 3));
+            qc.RandomSource = new SeededRandomSource(seed);
+
+            Assert.NotEqual("000", qc.Measure(0, 1, 2));
+        }
+    }
+
+    [Theory]
+    [InlineData(true, false, false)]
+    [InlineData(false, true, true)]
+    [InlineData(true, true, true)]
+    [InlineData(false, false, false)]
+    [InlineData(true, false, true)]
+    public void Bernstein_Vazirani_recovers_the_secret_in_one_query(bool b0, bool b1, bool b2)
+    {
+        bool[] secret = [b0, b1, b2];
+
+        QuantumCircuit qc = Algorithms.BernsteinVazirani(secret);
+        qc.RandomSource = new SeededRandomSource(3);
+
+        // Measure(2, 1, 0) puts qubit 2 first, matching the secret written most significant first.
+        string measured = qc.Measure(2, 1, 0);
+        string expected = string.Concat(secret.Reverse().Select(b => b ? '1' : '0'));
+
+        Assert.Equal(expected, measured);
+    }
+
+    [Theory]
+    [InlineData(0)]
+    [InlineData(1)]
+    [InlineData(2)]
+    [InlineData(3)]
+    public void Grover_finds_the_marked_item(int target)
+    {
+        // Phase-mark exactly one of the four basis states.
+        QuantumCircuit qc = Algorithms.Grover(2, c =>
+        {
+            if ((target & 1) == 0) c.X(0);
+            if ((target & 2) == 0) c.X(1);
+
+            c.CZ(0, 1);
+
+            if ((target & 1) == 0) c.X(0);
+            if ((target & 2) == 0) c.X(1);
+        });
+
+        var probabilities = qc.GetProbabilities();
+        string expected = Convert.ToString(target, 2).PadLeft(2, '0');
+
+        // One iteration on four items is exact: the marked state has probability 1.
+        Assert.Equal(1.0, probabilities[expected], 9);
+    }
+
+    [Fact]
+    public void Teleportation_delivers_the_message()
+    {
+        QuantumCircuit qc = Algorithms.Teleportation(c => c.Ry(0, 0.9));
+        qc.RandomSource = new SeededRandomSource(5);
+
+        double expected = Math.Sin(0.45) * Math.Sin(0.45);
+
+        Assert.Equal(expected, qc.QubitProbability(2), 9);
+    }
+
+    [Theory]
+    [InlineData(false, false, "00")]
+    [InlineData(false, true, "01")]
+    [InlineData(true, false, "10")]
+    [InlineData(true, true, "11")]
+    public void Superdense_coding_sends_two_bits_on_one_qubit(bool first, bool second, string expected)
+    {
+        QuantumCircuit qc = Algorithms.SuperdenseCoding(first, second);
+        qc.RandomSource = new SeededRandomSource(2);
+
+        // Qubit 0 carries the first bit, qubit 1 the second.
+        Assert.Equal(expected, qc.Measure(0, 1));
+    }
+
+    [Fact]
+    public void QFT_of_the_all_zero_state_is_a_uniform_superposition()
+    {
+        QuantumCircuit qc = new(3);
+        qc.QFT();
+
+        var probabilities = qc.GetProbabilities();
+
+        Assert.Equal(8, probabilities.Count);
+        Assert.All(probabilities.Values, p => Assert.Equal(0.125, p, 9));
+    }
+
+    [Fact]
+    public void QFT_preserves_normalization()
+    {
+        QuantumCircuit qc = new(4);
+        qc.X(0);
+        qc.H(2);
+        qc.QFT();
+
+        AssertNormalized(qc);
+    }
+
+    [Fact]
+    public void Algorithms_reject_a_null_oracle()
+    {
+        Assert.Throws<ArgumentNullException>(() => Algorithms.DeutschJozsa(2, null!));
+        Assert.Throws<ArgumentNullException>(() => Algorithms.Grover(2, null!));
+        Assert.Throws<ArgumentNullException>(() => Algorithms.BernsteinVazirani(null!));
+    }
+
+    [Fact]
+    public void Grover_over_more_than_three_qubits_reports_the_limit() =>
+        Assert.Throws<NotSupportedException>(() => Algorithms.Grover(4, c => c.Z(0)));
+}
diff --git a/tests/Qubit.NET.Tests/ExportTests.cs b/tests/Qubit.NET.Tests/ExportTests.cs
new file mode 100644
index 0000000..1c9abfa
--- /dev/null
+++ b/tests/Qubit.NET.Tests/ExportTests.cs
@@ -0,0 +1,204 @@
+using Qubit.NET;
+using Qubit.NET.Circuits;
+using Qubit.NET.Export;
+using Qubit.NET.Gates;
+using Qubit.NET.Utilities;
+using Qubit.NET.Visualization;
+
+namespace QubitNet.Tests;
+
+/// <summary>
+/// OpenQASM 2.0 export and the state-inspection helpers.
+/// </summary>
+public class ExportTests
+{
+    [Fact]
+    public void Bell_circuit_exports_to_valid_qasm()
+    {
+        QuantumCircuit qc = BellStates.PhiPlus();
+        qc.Measure();
+
+        string[] lines = qc.ToQasm()
+            .Split('\n')
+            .Select(l => l.Trim())
+            .Where(l => l.Length > 0)
+            .ToArray();
+
+        Assert.Equal("OPENQASM 2.0;", lines[0]);
+        Assert.Equal("include \"qelib1.inc\";", lines[1]);
+        Assert.Equal("qreg q[2];", lines[2]);
+        Assert.Equal("creg c[2];", lines[3]);
+        Assert.Equal("h q[0];", lines[4]);
+        Assert.Equal("cx q[0],q[1];", lines[5]);
+        Assert.Equal("measure q[0] -> c[0];", lines[6]);
+        Assert.Equal("measure q[1] -> c[1];", lines[7]);
+    }
+
+    [Fact]
+    public void Rotation_angles_are_exported_exactly()
+    {
+        QuantumCircuit qc = new(2);
+        qc.Rx(0, 0.25);
+        qc.Ry(0, -1.5);
+        qc.Rz(1, 3.14159);
+        qc.CRy(0, 1, 0.5);
+
+        string qasm = qc.ToQasm();
+
+        Assert.Contains("rx(0.25) q[0];", qasm);
+        Assert.Contains("ry(-1.5) q[0];", qasm);
+        Assert.Contains("rz(3.14159) q[1];", qasm);
+        Assert.Contains("cry(0.5) q[0],q[1];", qasm);
+    }
+
+    [Fact]
+    public void U3_exports_all_three_angles()
+    {
+        QuantumCircuit qc = new(1);
+        qc.U3(0, 0.1, 0.2, 0.3);
+
+        Assert.Contains("u3(0.1,0.2,0.3) q[0];", qc.ToQasm());
+    }
+
+    [Fact]
+    public void Three_qubit_gates_export()
+    {
+        QuantumCircuit qc = new(3);
+        qc.Toffoli(0, 1, 2);
+        qc.Fredkin(0, 1, 2);
+        qc.SWAP(0, 1);
+
+        string qasm = qc.ToQasm();
+
+        Assert.Contains("ccx q[0],q[1],q[2];", qasm);
+        Assert.Contains("cswap q[0],", qasm);
+        Assert.Contains("swap q[", qasm);
+    }
+
+    [Fact]
+    public void Conditional_gates_export_as_qasm_if_statements()
+    {
+        QuantumCircuit qc = Algorithms.Teleportation(c => c.Ry(0, 0.5));
+
+        string qasm = qc.ToQasm();
+
+        Assert.Contains("if (c[1]==1) x q[2];", qasm);
+        Assert.Contains("if (c[0]==1) z q[2];", qasm);
+    }
+
+    [Fact]
+    public void Initializations_export_as_their_gate_equivalents()
+    {
+        QuantumCircuit qc = new(3);
+        qc.Initialize(0, State.One);
+        qc.Initialize(1, State.Plus);
+        qc.Initialize(2, State.Minus);
+
+        string qasm = qc.ToQasm();
+
+        Assert.Contains("x q[0];", qasm);
+        Assert.Contains("h q[1];", qasm);
+        Assert.Contains("x q[2];", qasm);
+    }
+
+    [Fact]
+    public void Custom_gates_are_reported_rather_than_silently_dropped()
+    {
+        QuantumCircuit qc = new(1);
+        qc.Custom(QuantumGates.H, 0);
+
+        Assert.Throws<NotSupportedException>(() => qc.ToQasm());
+    }
+
+    [Fact]
+    public void Every_standard_gate_exports_without_throwing()
+    {
+        QuantumCircuit qc = new(3);
+        qc.I(0); qc.H(0); qc.X(0); qc.Y(0); qc.Z(0);
+        qc.S(0); qc.Sdag(0); qc.T(0); qc.Tdag(0);
+        qc.SX(0); qc.SY(0); qc.SZ(0);
+        qc.Rx(0, 0.1); qc.Ry(0, 0.1); qc.Rz(0, 0.1); qc.U3(0, 0.1, 0.2, 0.3);
+        qc.CNOT(0, 1); qc.CY(0, 1); qc.CZ(0, 1); qc.CH(0, 1);
+        qc.CRx(0, 1, 0.1); qc.CRy(0, 1, 0.1); qc.CRz(0, 1, 0.1); qc.CU3(0, 1, 0.1, 0.2, 0.3);
+        qc.SWAP(0, 1); qc.Toffoli(0, 1, 2); qc.Fredkin(0, 1, 2);
+        qc.Measure();
+
+        string qasm = qc.ToQasm();
+
+        // Nothing was dropped or commented out: 4 header lines (version, include, qreg,
+        // creg) + 27 gates + one measure statement per qubit.
+        Assert.DoesNotContain("//", qasm);
+        Assert.Equal(4 + 27 + 3, qasm.Split('\n').Count(l => l.TrimEnd().EndsWith(';')));
+    }
+
+    [Fact]
+    public void Bloch_vector_of_a_fresh_qubit_points_up() =>
+        AssertVector((0, 0, 1), new QuantumCircuit(1).BlochVector(0));
+
+    [Fact]
+    public void Bloch_vector_after_X_points_down()
+    {
+        QuantumCircuit qc = new(1);
+        qc.X(0);
+
+        AssertVector((0, 0, -1), qc.BlochVector(0));
+    }
+
+    [Fact]
+    public void Bloch_vector_after_H_points_along_x()
+    {
+        QuantumCircuit qc = new(1);
+        qc.H(0);
+
+        AssertVector((1, 0, 0), qc.BlochVector(0));
+    }
+
+    [Fact]
+    public void Bloch_vector_of_a_maximally_entangled_qubit_is_at_the_origin()
+    {
+        // Each half of a Bell pair is maximally mixed on its own.
+        AssertVector((0, 0, 0), BellStates.PhiPlus().BlochVector(0));
+        AssertVector((0, 0, 0), BellStates.PhiPlus().BlochVector(1));
+    }
+
+    [Fact]
+    public void Qubit_probability_reports_the_marginal()
+    {
+        QuantumCircuit qc = new(2);
+        qc.H(0);
+
+        Assert.Equal(0.5, qc.QubitProbability(0), 9);
+        Assert.Equal(0.0, qc.QubitProbability(1), 9);
+    }
+
+    [Fact]
+    public void Histogram_charts_the_outcome_probabilities()
+    {
+        string histogram = BellStates.PhiPlus().ToHistogram(20);
+        string[] lines = histogram.Split('\n');
+
+        Assert.Equal(2, lines.Length);
+        Assert.StartsWith("00 | ", lines[0]);
+        Assert.StartsWith("11 | ", lines[1]);
+        Assert.EndsWith("0.500", lines[0].TrimEnd());
+        Assert.Equal(20, lines[0].Count(c => c == '#'));
+    }
+
+    [Fact]
+    public void Inspection_helpers_reject_out_of_range_qubits()
+    {
+        QuantumCircuit qc = new(2);
+
+        Assert.Throws<ArgumentOutOfRangeException>(() => qc.BlochVector(2));
+        Assert.Throws<ArgumentOutOfRangeException>(() => qc.QubitProbability(-1));
+        Assert.Throws<ArgumentOutOfRangeException>(() => qc.ToHistogram(0));
+    }
+
+    private static void AssertVector((double X, double Y, double Z) expected,
+        (double X, double Y, double Z) actual)
+    {
+        Assert.Equal(expected.X, actual.X, 9);
+        Assert.Equal(expected.Y, actual.Y, 9);
+        Assert.Equal(expected.Z, actual.Z, 9);
+    }
+}

From e04e31c075f1c5f09818340aa72fb91f3b2e958a Mon Sep 17 00:00:00 2001
From: InfoTCube <infocube.yout@gmail.com>
Date: Mon, 17 Aug 2026 20:55:23 +0200
Subject: [PATCH 7/9] ci: publish to NuGet with trusted publishing instead of a
 stored API key

nuget.org now recommends OIDC-based trusted publishing over long-lived API
keys. The workflow exchanges a GitHub OIDC token for a key valid for one hour,
so there is no persistent credential to leak, and an intercepted token cannot
be replayed.

- Added id-token: write for OIDC, and contents: write, which the job no longer
  inherits once permissions are declared explicitly.
- NuGet/login@v1 runs immediately before the push so the hour-long key cannot
  expire mid-run.
- Replaced the NUGET_API_KEY secret with NUGET_USER, holding the nuget.org
  profile name.
- Documented the one-time nuget.org policy setup in the README, including that
  the workflow file name must be given without its path.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Mhwc4XeGtzM2UQx9sNLb3r
---
 .github/workflows/release.yml | 14 +++++++++++++-
 README.md                     | 30 ++++++++++++++++++++++++++++++
 2 files changed, 43 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index b7502b8..db92c91 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -8,6 +8,10 @@ jobs:
   publish:
     runs-on: ubuntu-latest
 
+    permissions:
+      id-token: write   # required for NuGet trusted publishing (OIDC)
+      contents: write   # required to create the GitHub release
+
     steps:
       - uses: actions/checkout@v4
         with:
@@ -34,10 +38,18 @@ jobs:
           dotnet pack src/Qubit.NET/Qubit.NET.csproj --no-build --configuration Release
           -p:Version=${{ env.VERSION }} --output ./artifacts
 
+      # Exchanges a GitHub OIDC token for a NuGet API key that lasts one hour, so no
+      # long-lived key is ever stored. Must run immediately before the push.
+      - name: NuGet login
+        uses: NuGet/login@v1
+        id: login
+        with:
+          user: ${{ secrets.NUGET_USER }}
+
       - name: Push to NuGet
         run: >
           dotnet nuget push "./artifacts/*.nupkg"
-          --api-key ${{ secrets.NUGET_API_KEY }}
+          --api-key ${{ steps.login.outputs.NUGET_API_KEY }}
           --source https://api.nuget.org/v3/index.json
           --skip-duplicate
 
diff --git a/README.md b/README.md
index a7a3a4e..01ddd7a 100644
--- a/README.md
+++ b/README.md
@@ -398,6 +398,36 @@ Memory is the real limit — the state vector holds 2ⁿ complex amplitudes at 1
 
 ---
 
+## 🚢 Releasing
+
+Publishing uses [NuGet Trusted Publishing](https://learn.microsoft.com/en-us/nuget/nuget-org/trusted-publishing),
+so no long-lived API key is stored anywhere. The workflow exchanges a GitHub OIDC token for
+a key that expires after an hour.
+
+One-time setup:
+
+1. On nuget.org: **your username → Trusted Publishing → add a policy**
+   - Repository Owner: `InfoTCube`
+   - Repository: `Qubit.NET`
+   - Workflow File: `release.yml` *(file name only, no path)*
+   - Environment: leave empty
+2. In GitHub repo settings, add a secret `NUGET_USER` holding your nuget.org **username**
+   (the profile name, not your email).
+
+To release:
+
+```bash
+git tag v1.0.0
+git push origin v1.0.0
+```
+
+The workflow builds all three target frameworks, runs the tests, packs, and pushes.
+
+> A new policy on a repository that nuget.org has not seen publish before is *temporarily
+> active for 7 days*. The first successful publish locks it to the repository permanently.
+
+---
+
 ## 💡 Contributions
 
 Pull requests, suggestions, and feature requests are welcome!  

From f0b7574414ad103e81a89e867d7cc238ad8a32b2 Mon Sep 17 00:00:00 2001
From: InfoTCube <infocube.yout@gmail.com>
Date: Mon, 17 Aug 2026 20:56:13 +0200
Subject: [PATCH 8/9] docs: drop the releasing section from the README

Releases are done by the maintainer only, so the publishing steps do not belong
in public docs.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Mhwc4XeGtzM2UQx9sNLb3r
---
 README.md | 30 ------------------------------
 1 file changed, 30 deletions(-)

diff --git a/README.md b/README.md
index 01ddd7a..a7a3a4e 100644
--- a/README.md
+++ b/README.md
@@ -398,36 +398,6 @@ Memory is the real limit — the state vector holds 2ⁿ complex amplitudes at 1
 
 ---
 
-## 🚢 Releasing
-
-Publishing uses [NuGet Trusted Publishing](https://learn.microsoft.com/en-us/nuget/nuget-org/trusted-publishing),
-so no long-lived API key is stored anywhere. The workflow exchanges a GitHub OIDC token for
-a key that expires after an hour.
-
-One-time setup:
-
-1. On nuget.org: **your username → Trusted Publishing → add a policy**
-   - Repository Owner: `InfoTCube`
-   - Repository: `Qubit.NET`
-   - Workflow File: `release.yml` *(file name only, no path)*
-   - Environment: leave empty
-2. In GitHub repo settings, add a secret `NUGET_USER` holding your nuget.org **username**
-   (the profile name, not your email).
-
-To release:
-
-```bash
-git tag v1.0.0
-git push origin v1.0.0
-```
-
-The workflow builds all three target frameworks, runs the tests, packs, and pushes.
-
-> A new policy on a repository that nuget.org has not seen publish before is *temporarily
-> active for 7 days*. The first successful publish locks it to the repository permanently.
-
----
-
 ## 💡 Contributions
 
 Pull requests, suggestions, and feature requests are welcome!  

From dffb6963b3b577eced97572decf06fc35af9dcb0 Mon Sep 17 00:00:00 2001
From: InfoTCube <infocube.yout@gmail.com>
Date: Mon, 17 Aug 2026 21:00:40 +0200
Subject: [PATCH 9/9] ci: refuse to publish from a tag that is not on main

The release workflow triggers on any v* tag, and a tag can point at any commit
on any branch, so a tag pushed on a feature branch would have published from
it. nuget.org's policy filters cannot catch this either: a tag push carries
ref_type "tag", so a branch filter never matches, and a tag filter cannot tell
whether the tag sits on main.

Verified before anything is built or any credential is exchanged.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Mhwc4XeGtzM2UQx9sNLb3r
---
 .github/workflows/release.yml | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index db92c91..dcc1ed8 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -17,6 +17,16 @@ jobs:
         with:
           fetch-depth: 0   # Source Link needs full history
 
+      # A v* tag can be pushed pointing at any commit on any branch, which would otherwise
+      # publish straight from a feature branch. Refuse anything not contained in main.
+      - name: Verify the tag is on main
+        run: |
+          git fetch --no-tags origin main
+          if ! git merge-base --is-ancestor "$GITHUB_SHA" origin/main; then
+            echo "::error::Tag $GITHUB_REF_NAME points at a commit that is not on main. Refusing to publish."
+            exit 1
+          fi
+
       - name: Setup .NET
         uses: actions/setup-dotnet@v4
         with: