|
36 | 36 |
|
37 | 37 | import math |
38 | 38 | import os |
| 39 | +from fractions import Fraction |
39 | 40 |
|
40 | 41 | import numpy as np |
41 | 42 |
|
|
58 | 59 |
|
59 | 60 | # --- the engine: bake the whole parameter loop into one homotopy -------------------------------- |
60 | 61 | # |
61 | | -# Node arithmetic promotes ints, exposes the exact constants bertini.Pi and bertini.I, and honours |
62 | | -# ** -- so the tree is written almost verbatim: x**degree, 2*bertini.Pi, bertini.I. The one thing |
63 | | -# the library will not do implicitly is turn a Python float/complex into a coefficient (a 16-digit |
64 | | -# literal would silently cap the arbitrary-precision tree), so a float parameter goes through its |
65 | | -# exact decimal string via _const. |
66 | | - |
67 | | -def _const(z): |
68 | | - """An exact constant node from a Python number, via coefficient's exact-string path: a real |
69 | | - number -> its decimal; a complex -> re + im * bertini.I.""" |
70 | | - z = complex(z) |
71 | | - node = bertini.coefficient(repr(z.real)) |
72 | | - return node if z.imag == 0 else node + bertini.coefficient(repr(z.imag)) * bertini.I |
73 | | - |
74 | | -def loom_homotopy(degree, center, radius, phi): |
75 | | - """H(x, t) = x^d - d*x - c(t), c(t) = center + radius*exp(i*(theta + phi)), theta = 2pi(1-t). |
76 | | -
|
77 | | - At t=1 theta=0 (the start configuration); at t=0 theta=2pi (the loop has closed). Tracking a |
78 | | - root of the start configuration from t=1 to t=0 carries it once around the loop. |
| 62 | +# Every constant here is EXACT and nothing reaches for a float: the loop parameters are rationals |
| 63 | +# (fractions.Fraction) and the phase is a rational multiple of pi. Node arithmetic promotes ints, |
| 64 | +# exposes the exact constants bertini.Pi and bertini.I, honours **, and bertini.coefficient accepts |
| 65 | +# exact values -- so the homotopy is written verbatim, with no coercion to defeat. |
| 66 | + |
| 67 | +def _c_of_theta(theta, center, radius): |
| 68 | + """The loop's moving coefficient c = center + radius * (cos theta + i sin theta), as a node. |
| 69 | + ``center`` / ``radius`` are exact (int or fractions.Fraction); ``theta`` is a node.""" |
| 70 | + return (bertini.coefficient(center) |
| 71 | + + bertini.coefficient(radius) * (bertini.cos(theta) + bertini.I * bertini.sin(theta))) |
| 72 | + |
| 73 | +def loom_homotopy(degree, center, radius, phi_over_pi): |
| 74 | + """H(x, t) = x^d - d*x - c(t), c(t) = center + radius*exp(i*(theta + phi)), |
| 75 | + theta = 2pi(1-t), phi = phi_over_pi * pi. |
| 76 | +
|
| 77 | + At t=1 theta=phi (the start configuration); at t=0 the loop has closed. Tracking a root of the |
| 78 | + start configuration from t=1 to t=0 carries it once around the loop. |
79 | 79 | """ |
80 | 80 | x = bertini.Variable('x') |
81 | 81 | t = bertini.Variable('t') |
82 | | - theta = 2 * bertini.Pi * (1 - t) + _const(phi) |
83 | | - c_t = _const(center) + _const(radius) * (bertini.cos(theta) + bertini.I * bertini.sin(theta)) |
| 82 | + theta = 2 * bertini.Pi * (1 - t) + bertini.coefficient(phi_over_pi) * bertini.Pi |
84 | 83 | sys = bertini.System() |
85 | | - sys.add_function(x**degree - degree * x - c_t) |
| 84 | + sys.add_function(x**degree - degree * x - _c_of_theta(theta, center, radius)) |
86 | 85 | sys.add_path_variable(t) |
87 | 86 | sys.add_variable_group(bertini.VariableGroup([x])) |
88 | 87 | return sys |
89 | 88 |
|
90 | | -def start_configuration(degree, center, radius, phi): |
91 | | - """The roots of the start configuration f(x) = x^d - d*x - c0, c0 = c(theta=0).""" |
92 | | - c0 = center + radius * complex(math.cos(phi), math.sin(phi)) |
| 89 | +def start_configuration(degree, center, radius, phi_over_pi): |
| 90 | + """The roots of the start configuration f(x) = x^d - d*x - c(t=1), where theta = phi at t=1.""" |
93 | 91 | x = bertini.Variable('x') |
| 92 | + theta_start = bertini.coefficient(phi_over_pi) * bertini.Pi |
94 | 93 | sys = bertini.System() |
95 | 94 | sys.add_variable_group(bertini.VariableGroup([x])) |
96 | | - sys.add_function(x**degree - degree * x - _const(c0)) |
| 95 | + sys.add_function(x**degree - degree * x - _c_of_theta(theta_start, center, radius)) |
97 | 96 | solver = bertini.nag_algorithm.ZeroDimSolver(sys, mptype='adaptive') |
98 | 97 | solver.solve() |
99 | 98 | return [complex(s[0]) for s in solver.all_solutions()] |
100 | 99 |
|
101 | | -def branch_values(degree): |
102 | | - """The (d-1) branch values of x^d - d*x - c: c = x^d - d*x at each critical point (the |
103 | | - (d-1)-th roots of unity, where d*x^{d-1} - d = 0). All share one modulus.""" |
104 | | - crit = [complex(math.cos(2 * math.pi * k / (degree - 1)), |
105 | | - math.sin(2 * math.pi * k / (degree - 1))) for k in range(degree - 1)] |
106 | | - return [z**degree - degree * z for z in crit] |
107 | | - |
108 | | -def track_loop(degree, center, radius, phi, tol=1e-10): |
| 100 | +def track_loop(degree, center, radius, phi_over_pi, tol=1e-10): |
109 | 101 | """Track every strand once around the loop. Returns a list of per-strand dicts with the loop |
110 | 102 | angle ``theta``, the complex position ``x``, and the tracker diagnostics along the path.""" |
111 | | - H = loom_homotopy(degree, center, radius, phi) |
112 | | - roots = start_configuration(degree, center, radius, phi) |
| 103 | + H = loom_homotopy(degree, center, radius, phi_over_pi) |
| 104 | + roots = start_configuration(degree, center, radius, phi_over_pi) |
113 | 105 |
|
114 | 106 | tracker = bertini.AMPTracker(H) |
115 | 107 | tracker.setup(bertini.tracking.Predictor.RK4, tol, 1e6, |
@@ -224,20 +216,20 @@ def stress(ls): |
224 | 216 | def teaching_frame(out): |
225 | 217 | """x^3 - 3x - c: a loop around ONE branch point (c = +2) -> a single transposition.""" |
226 | 218 | bertini.random.set_random_seed(1) # deterministic start-root ordering -> stable render |
227 | | - # branch points at c = +/-2; put the loop centre right of +2 so it encircles +2, excludes -2, |
228 | | - # and its closest approach (the pinch) lands at theta = pi, mid-loop. |
229 | | - radius, graze = 1.5, 0.02 |
230 | | - center = 2.0 + (radius - graze) |
231 | | - strands, _ = track_loop(3, complex(center, 0.0), radius, phi=0.0) |
| 219 | + # x^d - d*x has branch points at c = +/-(d-1); for d=3 that is +/-2. Centre the loop right of |
| 220 | + # +2 so it encircles +2, excludes -2, and its closest approach (the pinch) is at theta = pi. |
| 221 | + graze, radius = Fraction(1, 50), Fraction(3, 2) # 0.02, 1.5 |
| 222 | + center = (3 - 1) + (radius - graze) # 2 + (radius - graze) |
| 223 | + strands, _ = track_loop(3, center, radius, phi_over_pi=0) |
232 | 224 | render(strands, out, |
233 | 225 | 'The Monodromy Loom — teaching case: $x^3 - 3x - c$', |
234 | 226 | 'loop one branch point → two roots SWAP (a transposition); the third rides straight') |
235 | 227 |
|
236 | 228 | def showpiece_frame(out): |
237 | 229 | """x^5 - 5x - c: a loop around ALL FOUR branch points -> a full 5-cycle, four pinches.""" |
238 | 230 | bertini.random.set_random_seed(1) # deterministic start-root ordering -> stable render |
239 | | - radius = abs(branch_values(5)[0]) + 0.15 # circle |c| = R just outside the branch orbit |
240 | | - strands, _ = track_loop(5, 0.0 + 0.0j, radius, phi=math.pi / 4) # phi keeps pinches off the seam |
| 231 | + radius = (5 - 1) + Fraction(3, 20) # circle |c| = R just outside the branch orbit |c| = 4 |
| 232 | + strands, _ = track_loop(5, 0, radius, phi_over_pi=Fraction(1, 4)) # phi = pi/4 keeps pinches off the seam |
241 | 233 | render(strands, out, |
242 | 234 | 'The Monodromy Loom — showpiece: $x^5 - 5x - c$', |
243 | 235 | 'loop encircles four branch points → a full 5-cycle; four pinches where the tracker sweats') |
|
0 commit comments