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

- Adds a discussion section in rotation_synthesis/README.md about how to
select the digits of precision
- updates to_quirk and related methods to have a `allow_global_phase`
parameter that affects which quirk gates are used (e.g. `T` vs
`Rz(pi/4)`)

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 $10\log_{10}{1/\epsilon}$ (obtained experimentally). 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,19 +132,28 @@ 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") -> 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 False, then each gate (except H) gets replaced by SU2 version:
+                    - S -> Rz(pi/2)
+                    - T -> Rz(pi/4)
+                    - X -> Rx(-pi)
+                    - Y -> Ry(-pi)
+                    - Z -> Rz(-pi)
+                Then a prefix composed of SXSX that corrects the phase difference caused by H gates
+                is added.
         Returns:
             A quirk link.
         Raises:
             ValueError: If twirl=True
         """
         if self.twirl:
             raise ValueError("to_quirk is not supported when twirl=True")
-        gates = ctr.to_quirk(self.to_matrix(), fmt)
+        gates = ctr.to_quirk(self.to_matrix(), fmt, allow_global_phase)
         cols = '[' + ','.join(f'[{g}]' for g in gates) + ']'
         return "https://algassert.com/quirk#circuit={\"cols\":%s}" % cols
 
@@ -277,30 +286,43 @@ def to_cirq(self, fmt: str = "xz", qs: Optional[Sequence[cirq.Qid]] = None) -> c
             ),
         )
 
-    def to_quirk(self, fmt: str = "xz") -> 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 False then each gate (except H) gets replaced by SU2 version:
+                    - S -> Rz(pi/2)
+                    - T -> Rz(pi/4)
+                    - X -> Rx(-pi)
+                    - Y -> Ry(-pi)
+                    - Z -> Rz(-pi)
+                Then a prefix composed of SXSX that corrects the phase difference caused by H gates
+                is added.
         Returns:
             A quirk link.
         """
         correction = self.correction
         if not isinstance(correction, UnitaryChannel):
             raise ValueError(f"to_quirk is not supported for correction of type {type(correction)}")
-        rot = ctr.to_quirk(self.rotation.to_matrix(), fmt)
-        cor = ctr.to_quirk(correction.to_matrix(), fmt)
+        rot = ctr.to_quirk(self.rotation.to_matrix(), fmt, allow_global_phase)
+        cor = ctr.to_quirk(correction.to_matrix(), fmt, allow_global_phase)
+        if allow_global_phase:
+            xgate = '"X"'
+        else:
+            xgate = '{"id":"Rzft","arg":"-pi"}'
         first_row = []
         second_row = []
         # CNOT
         first_row.append("\"•\"")
-        second_row.append("\"X\"")
+        second_row.append(xgate)
         # rotation
         first_row.extend(rot)
         second_row.extend("1" for _ in rot)
         # CNOT
         first_row.append("\"•\"")
-        second_row.append("\"X\"")
+        second_row.append(xgate)
         # measure
         first_row.append("1")
         second_row.append("\"Measure\"")

qualtran/rotation_synthesis/matrix/clifford_t_repr.py

@@ -52,26 +52,32 @@ def _xyz_sequence(matrix: su2_ct.SU2CliffordT) -> tuple[str, ...]:
 ]
 
 
-def _xz_sequence(matrix: su2_ct.SU2CliffordT, use_hs: bool = True) -> Optional[tuple[str, ...]]:
+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]
     if use_hs:
         cliffords.append(su2_ct.HSqrt2)
         cliffords.append(su2_ct.HSqrt2 @ su2_ct.SSqrt2)
     candidates = []
+    pref = prv.removesuffix('*')
     for name, t in _T_list:
+        if name.startswith(pref):
+            continue
         for c in cliffords:
             new = c.adjoint() @ matrix
             new = t.adjoint() @ new
             new = new.scale_down()
             if new is None or not new.is_valid():
                 continue
-            seq = _xz_sequence(new, False)
+            seq = _xz_sequence(new, False, name)
             if seq is None:
                 continue
             gates = cast(tuple[str, ...], c.gates)
             candidates.append(seq + (name,) + gates)
+            break
     if not candidates:
         return None
     return min(candidates, key=len)
@@ -113,19 +119,34 @@ def to_sequence(matrix: su2_ct.SU2CliffordT, fmt: str = 'xyz') -> tuple[str, ...
 }
 
 
-def _to_quirk_name(name: str) -> str:
-    if name == "I":
-        return "1"
-    if name in ("X", "Y", "Z", "H"):
-        return "\"" + name + "\""
-    if name == "S":
-        return "\"Z^½\""
-    if name == "S*":
-        return "\"Z^-½\""
-    if name.startswith("T"):
-        if name.endswith("*"):
-            return "\"" + name[1].upper() + "^-¼" + "\""
-        return "\"" + name[1].upper() + "^¼" + "\""
+def _to_quirk_name(name: str, allow_global_phase: bool = False) -> str:
+    if allow_global_phase:
+        if name == "I":
+            return "1"
+        if name in ("X", "Y", "Z", "H"):
+            return "\"" + name + "\""
+        if name == "S":
+            return "\"Z^½\""
+        if name == "S*":
+            return "\"Z^-½\""
+        if name.startswith("T"):
+            if name.endswith("*"):
+                return "\"" + name[1].upper() + "^-¼" + "\""
+            return "\"" + name[1].upper() + "^¼" + "\""
+    else:
+        if name == "I":
+            return "1"
+        if name == "H":
+            return "\"H\""
+        if name in ("X", "Y", "Z"):
+            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)
     raise ValueError(f"{name=} is not supported")
 
 
@@ -145,13 +166,29 @@ 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) -> 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:
         matrix: The matrix to represent.
         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:
+                - S -> Rz(pi/2)
+                - T -> Rz(pi/4)
+                - X -> Rx(-pi)
+                - Y -> Ry(-pi)
+                - Z -> Rz(-pi)
+            Then a prefix composed of SXSX that corrects the phase difference caused by H gates
+            is added.
     Returns:
         A tuple quirk symbols.
     """
-    return tuple(_to_quirk_name(name) for name in to_sequence(matrix, fmt))
+    sequence = to_sequence(matrix, fmt)
+    phase_correction = tuple[str, ...]()
+    if not allow_global_phase:
+        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

@@ -43,9 +43,18 @@ def test_to_xyz_seq(g):
 
 
 @pytest.mark.parametrize("g", _make_random_su(50, 5, random_cliffords=True, seed=0))
-def test_to_xz_seq(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
+    for i in range(len(seq)):
+        if seq[i].startswith('T'):
+            first_t = i
+            break
+    if first_t is not None:
+        ts = 'Tx', 'Tx*', 'Tz', 'Tz*'
+        assert all(s in ts for s in seq[first_t:-2]), f'{seq=}'
     got = su2_ct.SU2CliffordT.from_sequence(seq)
     assert got == g
 

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)