check

Author: jlapeyre

pygridsynth/gridsynth.py

@@ -246,7 +246,6 @@ def gridsynth(
                 epsilon=epsilon,
                 dps=dps,
                 loop_controller=loop_controller,
-                phase=phase,
                 verbose=verbose,
                 measure_time=measure_time,
                 show_graph=show_graph)
@@ -256,7 +255,6 @@ def gridsynth(
                 epsilon=epsilon,
                 dps=dps,
                 loop_controller=loop_controller,
-                phase=phase,
                 verbose=verbose,
                 measure_time=measure_time,
                 show_graph=show_graph)
@@ -268,7 +266,6 @@ def _gridsynth_exact(
     epsilon,
     dps=None,
     loop_controller=None,
-    phase=False,
     verbose=False,
     measure_time=False,
     show_graph=False,
@@ -350,7 +347,6 @@ def _gridsynth_upto_phase(
     epsilon,
     dps=None,
     loop_controller=None,
-    phase=False,
     verbose=False,
     measure_time=False,
     show_graph=False,
@@ -393,39 +389,35 @@ def _gridsynth_upto_phase(
         if measure_time:
             start = time.time()
 
-        sol = solve_TDGP(
+        sol0 = _candidates(
             epsilon_region,
             unit_disk,
-            *transformed,
-            k,
-            phase=phase,
+            transformed,
+            phase=False,
             verbose=verbose,
             show_graph=show_graph,
         )
 
-        sol1 = solve_TDGP(
+        sol1 = _candidates(
             epsilon_region1,
             unit_disk1,
-            *transformed1,
-            k,
-            phase=phase,
+            transformed1,
+            phase=True,
             verbose=verbose,
             show_graph=show_graph,
         )
 
-        if measure_time:
-            time_of_solve_TDGP += time.time() - start
+        candidates = all_candidates(sol0, sol1)
 
         start = time.time()
 
-#        for (z, z1) in zip(sol, sol1):
-        for z in sol1:
-
-            z = z * inv_delta
-
+        for (tcount, z, phase) in candidates():
+            if (z * z.conj).residue == 0 and not phase:
+                continue
             xi = 1 - DRootTwo.fromDOmega(z.conj * z)
             w = diophantine_dyadic(xi, loop_controller=loop_controller)
             if w != NO_SOLUTION:
+#                print(f"{tcount} {phase}")
                 z = z.reduce_denomexp()
                 w = w.reduce_denomexp()
 
@@ -434,10 +426,16 @@ def _gridsynth_upto_phase(
                 elif z.k < w.k:
                     z = z.renew_denomexp(w.k)
 
-                if (z + w).reduce_denomexp().k < z.k:
-                    u_approx = DOmegaUnitary(z, w, 7)
+                if phase:
+                    if (z + w).reduce_denomexp().k < z.k:
+                        u_approx = DOmegaUnitary(z, w, 7)
+                    else:
+                        u_approx = DOmegaUnitary(z, w.mul_by_omega_inv(), 7)
                 else:
-                    u_approx = DOmegaUnitary(z, w.mul_by_omega_inv(), 7)
+                    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
@@ -448,12 +446,71 @@ def _gridsynth_upto_phase(
                     )
                 if verbose:
                     print(f"{z=}, {w=}")
+                    print(f"tcount={tcount} phase={phase}")
                     print("------------------")
+
                 return u_approx
         if measure_time:
             time_of_diophantine_dyadic += time.time() - start
-        k += 1
 
+def _candidates(
+    epsilon_region,
+    disk,
+    transformed,
+    phase=False,
+    verbose=False,
+    show_graph=False,
+#    measure_time=False,
+):
+    if phase:
+        inv_delta = DOmega(ZOmega(0, -1, 1, 0), 1)
+    def icandidates():
+        k = 0
+        while True:
+            sol = solve_TDGP(
+                epsilon_region,
+                disk,
+                *transformed,
+                k,
+                phase=phase,
+                verbose=verbose,
+                show_graph=show_graph,
+            )
+            # if measure_time:
+            #     time_of_solve_TDGP += time.time() - start
+
+            for z in sol:
+                klde = z.k
+                if klde > 0:
+                    tcount = 2*klde - 2
+                else:
+                    tcount = 0
+                if phase:
+                    tcount += 1
+                    z = z * inv_delta
+                yield (tcount, z, phase)
+            k += 1
+
+    return icandidates
+
+def all_candidates(icandidates0, icandidates1):
+    c0 = icandidates0()
+    c1 = icandidates1()
+    def itall_candidates():
+        tup0 = next(c0)
+        tup1 = next(c1)
+
+        while True:
+            tcount0 = tup0[0]
+            tcount1 = tup1[0]
+
+            if tcount0 < tcount1:
+                yield tup0
+                tup0 = next(c0)
+            else:
+                yield tup1
+                tup1 = next(c1)
+    return itall_candidates
 
 def gridsynth_circuit(
     theta,