add a discussion section in rotation_synthesis/README.md about digits of precision and add a parameter to control global phase of to_quirk

Author: NoureldinYosri

qualtran/rotation_synthesis/README.md

@@ -116,3 +116,15 @@ expected number of T gates is 18.982
 >>> 'actual diamond distance: %e'%mixed_fallback.diamond_norm_distance_to_rz(theta, config)
 'actual diamond distance: 5.397363e-10'
 ```
+
+## Effect of digits of precision
+The number of digits of precision used (i.e. rs.with_dps(digits_of_precision)) affects the result of the synthesis as follows:
+
+- very low => A math error will be raised by one of the checks|
+- low => A valid solution may be missed, in other words you get either a solution that has more T gates or None|
+- just right => A correct solution that has a number of T gates on par with the current state of the art|
+- high => same solution as above but in more time|
+
+Essentially, the code will either return a valid synthesis or None. If the code returns a result then it may be improved by increasing the number of digits of precision and if the code returns None then we need to increase the number of digits of precisions or `max_n` or both.
+
+As a rule of thumb, the number of digits of precision should be close to $9\log_10{1/\epsilon}$. This works for large $epsilon$ and is an upperbound for very small $\epsilon$, for example for $\epsilon=10^{-50}$ we need 400 digits and for $\epsilon=10^{-100}$ we need 800 digits.

qualtran/rotation_synthesis/channels/channel.py

@@ -132,13 +132,13 @@ def to_cirq(self, fmt: str = "xz", qs: Optional[Sequence[cirq.Qid]] = None) -> c
             q = cirq.q(0)
         return cirq.Circuit(ctr.to_cirq(self.to_matrix(), fmt, q))
 
-    def to_quirk(self, fmt: str = "xz", allow_global_phase: bool = False) -> str:
+    def to_quirk(self, fmt: str = "xz", allow_global_phase: bool = True) -> str:
         """Retruns a quirk link representing the channel operation.
 
         Args:
             fmt: The gates to use (see the documentation of to_sequence).
             allow_global_phase: whether the result can have a global phase or not.
-                If this is set then each gate (except H) gets replaced by SU2 version:
+                If this is False, then each gate (except H) gets replaced by SU2 version:
                     - S -> Rz(pi/2)
                     - T -> Rz(pi/4)
                     - X -> Rx(-pi)
@@ -286,13 +286,13 @@ def to_cirq(self, fmt: str = "xz", qs: Optional[Sequence[cirq.Qid]] = None) -> c
             ),
         )
 
-    def to_quirk(self, fmt: str = "xz", allow_global_phase: bool = False) -> str:
+    def to_quirk(self, fmt: str = "xz", allow_global_phase: bool = True) -> str:
         """Retruns a quirk link representing the channel operation.
 
         Args:
             fmt: The gates to use (see the documentation of to_sequence).
             allow_global_phase: whether the result can have a global phase or not.
-                If this is set then each gate (except H) gets replaced by SU2 version:
+                If this is False then each gate (except H) gets replaced by SU2 version:
                     - S -> Rz(pi/2)
                     - T -> Rz(pi/4)
                     - X -> Rx(-pi)

qualtran/rotation_synthesis/matrix/clifford_t_repr.py

@@ -52,9 +52,9 @@ def _xyz_sequence(matrix: su2_ct.SU2CliffordT) -> tuple[str, ...]:
 ]
 
 
-def _xz_sequence(matrix: su2_ct.SU2CliffordT, use_hs: bool = True, prv: str = 'dummy') -> Optional[tuple[str, ...]]:
-    # if use_hs:
-        # matrix = matrix.reduce()
+def _xz_sequence(
+    matrix: su2_ct.SU2CliffordT, use_hs: bool = True, prv: str = 'dummy'
+) -> Optional[tuple[str, ...]]:
     if matrix.det() == 2:
         return clifford(matrix)
     cliffords = [su2_ct.ISqrt2]
@@ -64,7 +64,8 @@ def _xz_sequence(matrix: su2_ct.SU2CliffordT, use_hs: bool = True, prv: str = 'd
     candidates = []
     pref = prv.removesuffix('*')
     for name, t in _T_list:
-        if name.startswith(pref): continue
+        if name.startswith(pref):
+            continue
         for c in cliffords:
             new = c.adjoint() @ matrix
             new = t.adjoint() @ new
@@ -138,14 +139,14 @@ def _to_quirk_name(name: str, allow_global_phase: bool = False) -> str:
         if name == "H":
             return "\"H\""
         if name in ("X", "Y", "Z"):
-            return '{"id":"R%sft","arg":"-pi"}'%(name.lower())
+            return '{"id":"R%sft","arg":"-pi"}' % (name.lower())
         if name == "S":
             return '{"id":"Rzft","arg":"pi/2"}'
         if name == "S*":
             return '{"id":"Rzft","arg":"-pi/2"}'
         if name.startswith("T"):
             angle = ['pi/4', '-pi/4'][name.endswith('*')]
-            return '{"id":"R%sft","arg":"%s"}'%(name[1].lower(), angle)
+            return '{"id":"R%sft","arg":"%s"}' % (name[1].lower(), angle)
     raise ValueError(f"{name=} is not supported")
 
 
@@ -165,7 +166,9 @@ def to_cirq(
     return tuple(_CIRQ_GATE_MAP[g](q) for g in to_sequence(matrix, fmt))
 
 
-def to_quirk(matrix: su2_ct.SU2CliffordT, fmt: str, allow_global_phase: bool = False) -> tuple[str, ...]:
+def to_quirk(
+    matrix: su2_ct.SU2CliffordT, fmt: str, allow_global_phase: bool = False
+) -> tuple[str, ...]:
     """Retruns a representation of the matrix as a sequence of quirk symbols.
 
     Args:
@@ -184,8 +187,8 @@ def to_quirk(matrix: su2_ct.SU2CliffordT, fmt: str, allow_global_phase: bool = F
         A tuple quirk symbols.
     """
     sequence = to_sequence(matrix, fmt)
-    phase_correction = ()
+    phase_correction = tuple[str, ...]()
     if not allow_global_phase:
-        phase = sum(g=='H' for g in sequence) % 4
+        phase = sum(g == 'H' for g in sequence) % 4
         phase_correction = ('"Z^½"', '"X"', '"Z^½"', '"X"') * phase
     return phase_correction + tuple(_to_quirk_name(name, allow_global_phase) for name in sequence)

qualtran/rotation_synthesis/matrix/clifford_t_repr_test.py

@@ -41,18 +41,10 @@ def test_to_xyz_seq(g):
     got = su2_ct.SU2CliffordT.from_sequence(seq)
     assert got == g
 
-def reduce(g: su2_ct.SU2CliffordT):
-    import qualtran.rotation_synthesis.rings.zsqrt2 as z2
-    import qualtran.rotation_synthesis.rings.zw as zw
-    if not(g.det() > 2 * z2.LAMBDA_KLIUCHNIKOV):
-        return g
-    if not all(v.is_divisible_by(zw.LAMBDA_KLIUCHNIKOV) for v in g.matrix.flat):
-        return g
-    return reduce(su2_ct.SU2CliffordT([[v//zw.LAMBDA_KLIUCHNIKOV for v in r] for r in g.matrix]))
 
 @pytest.mark.parametrize("g", _make_random_su(50, 5, random_cliffords=True, seed=0))
-def test_to_xz_seq(g):
-    g = reduce(g)
+def test_to_xz_seq(g: su2_ct.SU2CliffordT):
+    g = g.rescale()
     seq = ctr.to_sequence(g, 'xz')
     assert not any('Ty' in g for g in seq)
     first_t = None

qualtran/rotation_synthesis/matrix/su2_ct.py

@@ -221,9 +221,12 @@ def rescale(self) -> 'SU2CliffordT':
         while u.det() > 2 * zsqrt2.LAMBDA_KLIUCHNIKOV:
             if not all(a.is_divisible_by(zw.LAMBDA_KLIUCHNIKOV) for a in u.matrix.flat):
                 break
-            u = SU2CliffordT(
+            new_u = SU2CliffordT(
                 [[x // zw.LAMBDA_KLIUCHNIKOV for x in row] for row in u.matrix], u.gates
             )
+            if not new_u.is_valid():
+                break
+            u = new_u
         return u
 
     def num_t_gates(self) -> int:

qualtran/rotation_synthesis/protocols/clifford_t_synthesis.py

@@ -62,8 +62,6 @@ def _solve(
         m = n
         real_bound_fn = protocol.make_real_bound_fn(m, config)
         if not real_bound_fn(p):
-            # if verbose:
-            #     print('skip', p)
             bad += new
             if verbose and bad % 10**5 == 0:
                 print(f"at {n=}", "seen", bad, "out of", total, "points, ratio =", bad / total)