Fixed dynamic phase gradient gate decohering qubits due to precision …

…loss

src/gates/MultiplyAccumulateGates.js

@@ -154,4 +154,4 @@ MultiplyAccumulateGates.all = [
     ...MultiplyAccumulateGates.SquareSubtractInputFamily.all,
 ];
 
-export {MultiplyAccumulateGates, BIG_MUL_MOD_SHADER_CODE}
+export {MultiplyAccumulateGates, BIG_MUL_MOD_SHADER_CODE, MUL_STEP}

src/gates/PhaseGradientGates.js

@@ -12,15 +12,32 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
+import {Config} from "src/Config.js";
 import {Gate} from "src/circuit/Gate.js"
 import {GatePainting} from "src/draw/GatePainting.js"
 import {ketArgs, ketShaderPhase} from "src/circuit/KetShaderUtil.js"
+import {MUL_STEP} from "src/gates/MultiplyAccumulateGates.js"
 import {WglArg} from "src/webgl/WglArg.js"
 
 const PHASE_GRADIENT_SHADER = ketShaderPhase(
-    'uniform float factor;',
     `
-        float angle = out_id * factor;
+        uniform float factor;
+
+        /// Scales an angle by an integer factor.
+        /// Performs the multiplication gradually, to avoid losing precision.
+        float angle_mul(float base_angle, float whole_factor) {
+            float result = 0.0;
+            for (int k = 0; k < ${Math.ceil(Config.MAX_WIRE_COUNT/MUL_STEP)}; k++) {
+                result += base_angle * mod(whole_factor, ${1<<MUL_STEP}.0);
+                result = mod(result, 6.283185307179586476925286766559);
+                whole_factor = floor(whole_factor / ${1<<MUL_STEP}.0);
+                base_angle = mod(base_angle * ${1<<MUL_STEP}.0, 6.283185307179586476925286766559);
+            }
+            return result;
+        }
+    `,
+    `
+        float angle = angle_mul(factor, out_id);
         return vec2(cos(angle), sin(angle));
     `);
 

test/circuit/CircuitStats.test.js

@@ -16,6 +16,7 @@ import {Suite, assertThat, assertTrue} from "test/TestUtil.js"
 import {CircuitStats} from "src/circuit/CircuitStats.js"
 
 import {CircuitDefinition} from "src/circuit/CircuitDefinition.js"
+import {Complex} from "src/math/Complex.js"
 import {GateColumn} from "src/circuit/GateColumn.js"
 import {Gates} from "src/gates/AllGates.js"
 import {Matrix} from "src/math/Matrix.js"
@@ -41,7 +42,7 @@ const circuit = (diagram, ...extras) => CircuitDefinition.fromTextDiagram(new Ma
     ['-', undefined],
     ['+', undefined],
     ['|', undefined],
-    ['/', undefined]
+    ['/', null]
 ]), diagram);
 
 suite.testUsingWebGL("empty", () => {
@@ -117,7 +118,7 @@ suite.testUsingWebGL('incoherent-amplitude-display', () => {
 suite.testUsingWebGL('coherent-amplitude-display', () => {
     let c = circuit(`-H-•-a/--
                      ---X-//--
-                     -H-------`, ['a', Gates.Displays.AmplitudeDisplayFamily.ofSize(2)]);
+                     -H-------`, ['a', Gates.Displays.AmplitudeDisplayFamily]);
     let stats = CircuitStats.fromCircuitAtTime(c, 0);
     assertThat(stats.qubitDensityMatrix(Infinity, 0)).isApproximatelyEqualTo(Matrix.square(0.5, 0, 0, 0.5));
     assertThat(stats.qubitDensityMatrix(Infinity, 1)).isApproximatelyEqualTo(Matrix.square(0.5, 0, 0, 0.5));
@@ -141,7 +142,7 @@ suite.testUsingWebGL('conditional-bloch-display', () => {
 
 suite.testUsingWebGL('probability-display', () => {
     let c = circuit(`-H-•-%-
-                     ---X-/-`, ['%', Gates.Displays.ProbabilityDisplayFamily.ofSize(2)]);
+                     ---X-/-`, ['%', Gates.Displays.ProbabilityDisplayFamily]);
     let stats = CircuitStats.fromCircuitAtTime(c, 0);
     assertThat(stats.qubitDensityMatrix(Infinity, 0)).isApproximatelyEqualTo(Matrix.square(0.5, 0, 0, 0.5));
     assertThat(stats.customStatsForSlot(5, 0)).isApproximatelyEqualTo(
@@ -151,15 +152,15 @@ suite.testUsingWebGL('probability-display', () => {
 suite.testUsingWebGL('controlled-multi-probability-display', () => {
     let c = circuit(`---◦-
                      -H-%-
-                     ---/-`, ['%', Gates.Displays.ProbabilityDisplayFamily.ofSize(2)]);
+                     ---/-`, ['%', Gates.Displays.ProbabilityDisplayFamily]);
     let stats = CircuitStats.fromCircuitAtTime(c, 0);
     assertThat(stats.customStatsForSlot(3, 1)).isApproximatelyEqualTo(
         Matrix.col(0.5, 0.5, 0, 0));
 });
 
 suite.testUsingWebGL('density-display', () => {
     let c = circuit(`-d/-
-                     -//-`, ['d', Gates.Displays.DensityMatrixDisplayFamily.ofSize(2)]);
+                     -//-`, ['d', Gates.Displays.DensityMatrixDisplayFamily]);
     let stats = CircuitStats.fromCircuitAtTime(c, 0);
     assertThat(stats.customStatsForSlot(1, 0)).isApproximatelyEqualTo(
         Matrix.square(
@@ -172,7 +173,7 @@ suite.testUsingWebGL('density-display', () => {
 suite.testUsingWebGL('shifted-density-display', () => {
     let c = circuit(`----
                      -d/-
-                     -//-`, ['d', Gates.Displays.DensityMatrixDisplayFamily.ofSize(2)]);
+                     -//-`, ['d', Gates.Displays.DensityMatrixDisplayFamily]);
     let stats = CircuitStats.fromCircuitAtTime(c, 0);
     assertThat(stats.customStatsForSlot(1, 1)).isApproximatelyEqualTo(
         Matrix.square(
@@ -262,3 +263,31 @@ suite.testUsingWebGL('survival-rates-controlled-postselection', () => {
     assertThat(stats.survivalRate(13)).isApproximatelyEqualTo(0.5);
     assertThat(stats.survivalRate(14)).isApproximatelyEqualTo(0.5);
 });
+
+suite.testUsingWebGL('dynamic-phase-gradient-keeps-qubits-coherent', () => {
+    let stats = CircuitStats.fromCircuitAtTime(
+        circuit(`-H-P-
+                 -H-/-
+                 -H-/-
+                 -H-/-
+                 -H-/-
+                 -H-/-
+                 -H-/-
+                 -H-/-
+                 -H-/-
+                 -H-/-
+                 -H-/-
+                 -H-/-
+                 -H-/-
+                 -H-/-
+                 -H-/-
+                 -H-/-`, ['P', Gates.PhaseGradientGates.DynamicPhaseGradientFamily]),
+        0.9);
+
+    // Check coherence of each qubit.
+    for (let i = 0; i < 16; i++) {
+        let [x, y, z] = stats.qubitDensityMatrix(Infinity, i).qubitDensityMatrixToBlochVector();
+        let r = x*x + y*y + z*z;
+        assertThat(r).withInfo({i, x, y, z}).isApproximatelyEqualTo(1);
+    }
+});