-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathriley.py
190 lines (149 loc) · 7.14 KB
/
riley.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
""" Methods to compute globally and locally with the Riley slice.
In this file are found methods to compute with the moduli space known as the `Riley slice'.
"""
import farey
from mpmath import mp
mp.dps = 100
import math
import scipy.optimize
from numpy.polynomial import Polynomial as P
try:
import mpsolve
mpsolve_avail = True
except ImportError:
mpsolve_avail = False
try:
import sympy
sympy_avail = True
except ImportError:
sympy_avail = False
try:
import matlab.engine
matlab_avail = True
matlab_eng = None
except ImportError:
matlab_avail = False
def poly_solve(poly, solver='mpsolve' if mpsolve_avail else 'scipy', max_iter=100, tol=1e-2, try_int=False):
""" Solve a polynomial numerically.
There are three possible solvers: the mpsolve solver (see https://numpi.dm.unipi.it/software/mpsolve), the solver built
in to scipy followed by an attempt to improve the results using Newton's algorithm, or sympy.
Arguments:
poly --- a scipy polynomial
solver -- one of 'mpsolve', 'scipy', 'sympy'
The following arguments are only used by the scipy solver:
max_iter -- maximum number of iterations for Newton's algorithm (default 100)
tol -- allowable error of each point (default 1e-2)
The following arguments are only used by the mpsolve solver:
try_int -- assume that the coefficients of poly are integral (default False)
"""
q = poly.degree()
if solver == 'mpsolve':
if not mpsolve_avail:
raise RuntimeError('mpsolve selected even though it seems not to be installed')
mpsolve_ctx = mpsolve.Context()
mpsolve_poly = mpsolve.MonomialPoly(mpsolve_ctx, q)
for d in range(0,q+1):
if try_int:
mpsolve_poly.set_coefficient(d, int(poly.coef[d]))
else:
mpsolve_poly.set_coefficient(d, float(mp.re(poly.coef[d])))
return mpsolve_ctx.solve(mpsolve_poly)
elif solver == 'sympy':
if not sympy_avail:
raise RuntimeError('sympy selected even though it seems not to be installed')
ZZ = sympy.Symbol('Z')
sympy_poly = 0
for d in range(0,q+1):
sympy_poly = sympy_poly + poly.coef[d] * ZZ**d
return [sympy.N(root) for root in sympy.solve(sympy_poly,ZZ)]
elif solver == 'scipy':
poly = P(list(map(float,poly.coef))) # Need to cast to float in case we have mpmath but not mpsolve
roots_bad = poly.roots()
return [scipy.optimize.newton(poly, root, poly.deriv(),fprime2=poly.deriv(2),maxiter=max_iter,tol=tol) for root in roots_bad]
elif solver == 'matlab':
if not matlab_avail:
raise RuntimeError('matlab selected even though it seems not to be installed')
global matlab_eng
if matlab_eng == None:
matlab_eng = matlab.engine.start_matlab()
roots = matlab_eng.roots(matlab.double([complex(x) for x in poly.coef],is_complex=True))
try:
return list(roots)
except TypeError:
return [roots]
raise RuntimeError(f'unknown solver {solver}')
def riley_slice(a, b, max_denom, solver='mpsolve' if mpsolve_avail else 'scipy', **kwargs):
""" Return an accurate approximation to the Riley slice.
There are three possible solvers: the mpsolve solver (see https://numpi.dm.unipi.it/software/mpsolve), the solver built
in to scipy followed by an attempt to improve the results using Newton's algorithm, or sympy.
Arguments:
a, b -- the order of the cone points represented by X and Y respectively. Use mp.inf for the parabolic case (or 1, since exp(2*pi*i/1) = exp(0) = 1).
max_denom -- the maximum denominator Farey polynomials to compute
solver -- one of 'mpsolve', 'scipy', 'sympy'
Further keyword arguments are passed directly to poly_solve(), i.e. tol and max_iter for scipy.
"""
alpha = 1 if a == mp.inf else mp.exp(1j*mp.pi/a)
beta = 1 if b == mp.inf else mp.exp(1j*mp.pi/b)
points = []
for q in range(1,max_denom+1):
for p in range(1,q+1):
if math.gcd(p,q) == 1:
poly = farey.polynomial_coefficients_fast(p, q, alpha, beta) + 2
try_int = True if alpha == 1 and beta == 1 else False
try:
points.extend(poly_solve(poly, solver, try_int=try_int, **kwargs))
except Exception as e:
raise RuntimeError(f'p = {p} q={q}') from e
return points
def riley_centre(a,b):
""" Return an approximation to the centre of symmetry of the Riley slice.
Arguments:
a,b -- orders of X and Y respectively
"""
alpha = 1 if a == mp.inf else mp.exp(1j*mp.pi/a)
beta = 1 if b == mp.inf else mp.exp(1j*mp.pi/b)
poly_left = farey.polynomial_coefficients_fast(1, 1, alpha, beta, int if (alpha == 1 and beta == 1) else mp.mpf) + 2
roots_left = poly_left.roots()
poly_right = farey.polynomial_coefficients_fast(0, 1, alpha, beta, int if (alpha == 1 and beta == 1) else mp.mpf) + 2
roots_right = poly_right.roots()
return (roots_right[0] + roots_left[0])/2
def cusp_point(a, b, p, q, solver='mpsolve' if mpsolve_avail else 'scipy', **kwargs):
""" Return an approximation to the p/q-cusp point on the Riley slice boundary.
Let Phi_{p/q} be the p/q-Farey polynomial; the p/q-pleating ray is then the connected
component of Phi_{p/q}^{-1}((-\\infty,-2)) with asymptotic slope pi*p/q, and the p/q-cusp
is the point on this branch which is the inverse image of -2.
Arguments:
a,b -- orders of X and Y respectively
p,q -- integers representing the slope of the desired cusp; p and q must be coprime and nonnegative and p/q must lie in [0,2).
solver -- one of 'mpsolve', 'scipy', 'sympy'
Further keyword arguments are passed directly to poly_solve(), i.e. tol and max_iter for scipy.
"""
if p/q > 1:
return mp.conj(cusp_point(a,b,2*q-p,q))
alpha = 1 if a == mp.inf else mp.exp(1j*mp.pi/a)
beta = 1 if b == mp.inf else mp.exp(1j*mp.pi/b)
poly = farey.polynomial_coefficients_fast(p, q, alpha, beta, int if (alpha == 1 and beta == 1) else mp.mpf) + 2
roots = poly_solve(poly, solver, **kwargs)
if q == 1:
return roots[0]
centre = riley_centre(a,b)
if q > 1:
(r1,s1),(r2,s2) = farey.neighbours(p,q)
def right_phase(angle):
if angle < 0:
return angle + 2*mp.pi
else:
return angle
left_cusp_angle = right_phase(mp.arg(cusp_point(a,b,r1,s1) - centre))
right_cusp_angle = right_phase(mp.arg(cusp_point(a,b,r2,s2) - centre))
right_argument_roots = []
for root in roots:
argument = mp.arg(root - centre)
#print(f'{left_cusp_angle} < {argument} < {right_cusp_angle}?')
if left_cusp_angle < argument and argument < right_cusp_angle:
right_argument_roots.append(root)
if right_argument_roots == []:
raise RuntimeError('failed to find cusp: couldn\'t bound argument')
return max(right_argument_roots, key=mp.fabs)
else:
raise ValueError('q < 1?')