allow theta and epsilon to be float

pygridsynth/cli.py

@@ -1,7 +1,5 @@
 import argparse
 
-import mpmath
-
 from .diophantine import set_random_seed
 from .gridsynth import gridsynth_gates
 from .loop_controller import LoopController
@@ -40,25 +38,16 @@ def main():
     # Set random seed for deterministic results
     set_random_seed(args.seed)
 
-    # Use the same heuristic as Haskell gridsynth for setting dps
-    if args.dps is None:
-        epsilon1 = mpmath.mpmathify(args.epsilon)
-        dps_of_result = -mpmath.log10(epsilon1)
-        args.dps = 15 + 2.5 * dps_of_result
-    mpmath.mp.dps = args.dps
-    epsilon = mpmath.mpmathify(args.epsilon)
-    mpmath.mp.pretty = True
-    theta = mpmath.mpmathify(args.theta)
-
     loop_controller = LoopController(
         dloop=args.dloop,
         floop=args.floop,
         dtimeout=args.dtimeout,
         ftimeout=args.ftimeout,
     )
     gates = gridsynth_gates(
-        theta=theta,
-        epsilon=epsilon,
+        theta=args.theta,
+        epsilon=args.epsilon,
+        dps=args.dps,
         loop_controller=loop_controller,
         verbose=args.verbose,
         measure_time=args.time,

pygridsynth/gridsynth.py

@@ -1,4 +1,5 @@
 import time
+import warnings
 
 import mpmath
 
@@ -84,11 +85,15 @@ def generate_complex_unitary(sol):
 
 
 def error(theta, gates):
-    u_target = Rz(theta)
-    u_approx = DOmegaUnitary.from_gates(gates).to_complex_matrix
-    return 2 * sqrt(
-        mpmath.re(1 - mpmath.re(mpmath.conj(u_target[0, 0]) * u_approx[0, 0]) ** 2)
-    )
+    tcount = gates.count("T")
+    epsilon = 2 ** (-tcount // 3)
+    dps = _dps_for_epsilon(epsilon)
+    with mpmath.workdps(dps):
+        u_target = Rz(theta)
+        u_approx = DOmegaUnitary.from_gates(gates).to_complex_matrix
+        return 2 * sqrt(
+            mpmath.re(1 - mpmath.re(mpmath.conj(u_target[0, 0]) * u_approx[0, 0]) ** 2)
+        )
 
 
 def check(theta, gates):
@@ -102,126 +107,236 @@ def check(theta, gates):
     print(f"{e=}")
 
 
+def get_synthesized_unitary(gates):
+    tcount = gates.count("T")
+    epsilon = 2 ** (-tcount // 3)
+    dps = _dps_for_epsilon(epsilon)
+    with mpmath.workdps(dps):
+        return DOmegaUnitary.from_gates(gates).to_complex_matrix
+
+
 def gridsynth(
     theta,
     epsilon,
+    dps=None,
     loop_controller=None,
     verbose=False,
     measure_time=False,
     show_graph=False,
 ):
-    if loop_controller is None:
-        loop_controller = LoopController()
+    if isinstance(theta, float):
+        warnings.warn(
+            (
+                f"pygridsynth is synthesizing the angle {theta}. "
+                "Please verify that this is the intended value. "
+                "Using float may introduce precision errors; "
+                "consider using mpmath.mpf for exact precision."
+            ),
+            UserWarning,
+            stacklevel=2,
+        )
 
-    epsilon_region = EpsilonRegion(theta, epsilon)
-    unit_disk = UnitDisk()
-    k = 0
+    if isinstance(epsilon, float):
+        warnings.warn(
+            (
+                f"pygridsynth is using epsilon={epsilon} as the tolerance. "
+                "Please verify that this is the intended value. "
+                "Using float may introduce precision errors; "
+                "consider using mpmath.mpf for exact precision."
+            ),
+            UserWarning,
+            stacklevel=2,
+        )
+
+    if dps is None:
+        dps = _dps_for_epsilon(epsilon)
+    with mpmath.workdps(dps):
+        theta = mpmath.mpmathify(theta)
+        epsilon = mpmath.mpmathify(epsilon)
+        if loop_controller is None:
+            loop_controller = LoopController()
+
+        epsilon_region = EpsilonRegion(theta, epsilon)
+        unit_disk = UnitDisk()
+        k = 0
 
-    start = time.time() if measure_time else 0.0
-    transformed = to_upright_set_pair(
-        epsilon_region, unit_disk, verbose=verbose, show_graph=show_graph
-    )
-    if measure_time:
-        print(f"to_upright_set_pair: {time.time() - start} s")
-    if verbose:
-        print("------------------")
-
-    time_of_solve_TDGP = 0
-    time_of_diophantine_dyadic = 0
-    while True:
         if measure_time:
             start = time.time()
-        sol = solve_TDGP(
-            epsilon_region,
-            unit_disk,
-            *transformed,
-            k,
-            verbose=verbose,
-            show_graph=show_graph,
+        transformed = to_upright_set_pair(
+            epsilon_region, unit_disk, verbose=verbose, show_graph=show_graph
         )
         if measure_time:
-            time_of_solve_TDGP += time.time() - start
-            start = time.time()
-        for z in sol:
-            if (z * z.conj).residue == 0:
-                continue
-            xi = 1 - DRootTwo.fromDOmega(z.conj * z)
-            w = diophantine_dyadic(xi, loop_controller=loop_controller)
-            if w != NO_SOLUTION:
-                z = z.reduce_denomexp()
-                w = w.reduce_denomexp()
-                if z.k > w.k:
-                    w = w.renew_denomexp(z.k)
-                elif z.k < w.k:
-                    z = z.renew_denomexp(w.k)
-                if (z + w).reduce_denomexp().k < z.k:
-                    u_approx = DOmegaUnitary(z, w, 0)
-                else:
-                    u_approx = DOmegaUnitary(z, w.mul_by_omega(), 0)
-                if measure_time:
-                    time_of_diophantine_dyadic += time.time() - start
-                    print(f"time of solve_TDGP: {time_of_solve_TDGP * 1000} ms")
-                    print(
-                        "time of diophantine_dyadic: "
-                        f"{time_of_diophantine_dyadic * 1000} ms"
-                    )
-                if verbose:
-                    print(f"{z=}, {w=}")
-                    print("------------------")
-                return u_approx
-        if measure_time:
-            time_of_diophantine_dyadic += time.time() - start
-        k += 1
+            print(f"to_upright_set_pair: {time.time() - start} s")
+        if verbose:
+            print("------------------")
+
+        time_of_solve_TDGP = 0
+        time_of_diophantine_dyadic = 0
+        while True:
+            if measure_time:
+                start = time.time()
+            sol = solve_TDGP(
+                epsilon_region,
+                unit_disk,
+                *transformed,
+                k,
+                verbose=verbose,
+                show_graph=show_graph,
+            )
+            if measure_time:
+                time_of_solve_TDGP += time.time() - start
+                start = time.time()
+
+            for z in sol:
+                if (z * z.conj).residue == 0:
+                    continue
+                xi = 1 - DRootTwo.fromDOmega(z.conj * z)
+                w = diophantine_dyadic(xi, loop_controller=loop_controller)
+                if w != NO_SOLUTION:
+                    z = z.reduce_denomexp()
+                    w = w.reduce_denomexp()
+                    if z.k > w.k:
+                        w = w.renew_denomexp(z.k)
+                    elif z.k < w.k:
+                        z = z.renew_denomexp(w.k)
+                    if (z + w).reduce_denomexp().k < z.k:
+                        u_approx = DOmegaUnitary(z, w, 0)
+                    else:
+                        u_approx = DOmegaUnitary(z, w.mul_by_omega(), 0)
+                    if measure_time:
+                        time_of_diophantine_dyadic += time.time() - start
+                        print(f"time of solve_TDGP: {time_of_solve_TDGP * 1000} ms")
+                        print(
+                            "time of diophantine_dyadic: "
+                            f"{time_of_diophantine_dyadic * 1000} ms"
+                        )
+                    if verbose:
+                        print(f"{z=}, {w=}")
+                        print("------------------")
+                    return u_approx
+            if measure_time:
+                time_of_diophantine_dyadic += time.time() - start
+            k += 1
 
 
 def gridsynth_circuit(
     theta,
     epsilon,
     wires=[0],
     decompose_phase_gate=True,
+    dps=None,
     loop_controller=None,
     verbose=False,
     measure_time=False,
     show_graph=False,
 ):
-    start_total = time.time() if measure_time else 0.0
-    u_approx = gridsynth(
-        theta=theta,
-        epsilon=epsilon,
-        loop_controller=loop_controller,
-        verbose=verbose,
-        measure_time=measure_time,
-        show_graph=show_graph,
-    )
+    if isinstance(theta, float):
+        warnings.warn(
+            (
+                f"pygridsynth is synthesizing the angle {theta}. "
+                "Please verify that this is the intended value. "
+                "Using float may introduce precision errors; "
+                "consider using mpmath.mpf for exact precision."
+            ),
+            UserWarning,
+            stacklevel=2,
+        )
 
-    start = time.time() if measure_time else 0.0
-    circuit = decompose_domega_unitary(
-        u_approx, wires=wires, decompose_phase_gate=decompose_phase_gate
-    )
-    if measure_time:
-        print(f"time of decompose_domega_unitary: {(time.time() - start) * 1000} ms")
-        print(f"total time: {(time.time() - start_total) * 1000} ms")
+    if isinstance(epsilon, float):
+        warnings.warn(
+            (
+                f"pygridsynth is using epsilon={epsilon} as the tolerance. "
+                "Please verify that this is the intended value. "
+                "Using float may introduce precision errors; "
+                "consider using mpmath.mpf for exact precision."
+            ),
+            UserWarning,
+            stacklevel=2,
+        )
 
-    return circuit
+    if dps is None:
+        dps = _dps_for_epsilon(epsilon)
+    with mpmath.workdps(dps):
+        theta = mpmath.mpmathify(theta)
+        epsilon = mpmath.mpmathify(epsilon)
+        start_total = time.time() if measure_time else 0.0
+        u_approx = gridsynth(
+            theta=theta,
+            epsilon=epsilon,
+            dps=dps,
+            loop_controller=loop_controller,
+            verbose=verbose,
+            measure_time=measure_time,
+            show_graph=show_graph,
+        )
+
+        start = time.time() if measure_time else 0.0
+        circuit = decompose_domega_unitary(
+            u_approx, wires=wires, decompose_phase_gate=decompose_phase_gate
+        )
+        if measure_time:
+            print(
+                f"time of decompose_domega_unitary: {(time.time() - start) * 1000} ms"
+            )
+            print(f"total time: {(time.time() - start_total) * 1000} ms")
+
+        return circuit
 
 
 def gridsynth_gates(
     theta,
     epsilon,
     decompose_phase_gate=True,
+    dps=None,
     loop_controller=None,
     verbose=False,
     measure_time=False,
     show_graph=False,
 ):
-    circuit = gridsynth_circuit(
-        theta=theta,
-        epsilon=epsilon,
-        wires=[0],
-        decompose_phase_gate=decompose_phase_gate,
-        loop_controller=loop_controller,
-        verbose=verbose,
-        measure_time=measure_time,
-        show_graph=show_graph,
-    )
-    return circuit.to_simple_str()
+    if isinstance(theta, float):
+        warnings.warn(
+            (
+                f"pygridsynth is synthesizing the angle {theta}. "
+                "Please verify that this is the intended value. "
+                "Using float may introduce precision errors; "
+                "consider using mpmath.mpf for exact precision."
+            ),
+            UserWarning,
+            stacklevel=2,
+        )
+
+    if isinstance(epsilon, float):
+        warnings.warn(
+            (
+                f"pygridsynth is using epsilon={epsilon} as the tolerance. "
+                "Please verify that this is the intended value. "
+                "Using float may introduce precision errors; "
+                "consider using mpmath.mpf for exact precision."
+            ),
+            UserWarning,
+            stacklevel=2,
+        )
+
+    if dps is None:
+        dps = _dps_for_epsilon(epsilon)
+    with mpmath.workdps(dps):
+        theta = mpmath.mpmathify(theta)
+        epsilon = mpmath.mpmathify(epsilon)
+        circuit = gridsynth_circuit(
+            theta=theta,
+            epsilon=epsilon,
+            wires=[0],
+            decompose_phase_gate=decompose_phase_gate,
+            dps=dps,
+            loop_controller=loop_controller,
+            verbose=verbose,
+            measure_time=measure_time,
+            show_graph=show_graph,
+        )
+        return circuit.to_simple_str()
+
+
+def _dps_for_epsilon(epsilon) -> int:
+    e = mpmath.mpmathify(epsilon)
+    k = -mpmath.log10(e)
+    return int(15 + 2.5 * int(mpmath.ceil(k)))  # used in newsynth

pygridsynth/synthesis_of_cliffordT.py

@@ -1,8 +1,7 @@
 from .normal_form import NormalForm
+from .quantum_gate import W_PHASE, HGate, QuantumCircuit, SGate, SXGate, TGate, WGate
 from .ring import OMEGA_POWER
 from .unitary import DOmegaUnitary
-from .quantum_gate import W_PHASE, QuantumCircuit, HGate, TGate, SGate, SXGate, WGate
-
 
 BIT_SHIFT = [0, 0, 1, 0, 2, 0, 1, 3, 3, 3, 0, 2, 2, 1, 0, 0]
 BIT_COUNT = [0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4]