From cbd7972d007ef71b0c029ca7e71c0f6e3da6c043 Mon Sep 17 00:00:00 2001 From: Marc Masdeu Date: Mon, 19 Feb 2024 16:55:37 +0100 Subject: [PATCH] Debugged code to work with hyperelliptic Mumford project --- darmonpoints/divisors.py | 2 +- darmonpoints/meromorphic.py | 100 ++++++--- darmonpoints/schottky.py | 402 +++++++++++++----------------------- darmonpoints/theta.py | 258 +++++++++++++++++++++++ 4 files changed, 468 insertions(+), 294 deletions(-) create mode 100644 darmonpoints/theta.py diff --git a/darmonpoints/divisors.py b/darmonpoints/divisors.py index 571bf1f..cfab0db 100644 --- a/darmonpoints/divisors.py +++ b/darmonpoints/divisors.py @@ -210,7 +210,7 @@ def left_act_by_matrix(self, g): if P == Infinity: try: new_pt = K(a) / K(c) - except ZeroDivisionError: + except (PrecisionError,ZeroDivisionError): new_pt = Infinity else: new_pt = (K(a) * P + K(b)) / (K(c) * P + K(d)) diff --git a/darmonpoints/meromorphic.py b/darmonpoints/meromorphic.py index b25f81e..f9f5baf 100644 --- a/darmonpoints/meromorphic.py +++ b/darmonpoints/meromorphic.py @@ -39,6 +39,15 @@ from .divisors import Divisors +def evalpoly(poly, x): + # Initialize result + result = poly[-1] + # Evaluate value of polynomial + # using Horner's method + for a in reversed(poly): + result *= x + result += a + return result class MeromorphicFunctionsElement(ModuleElement): def __init__(self, parent, data, parameter, check=True): @@ -54,18 +63,19 @@ def __init__(self, parent, data, parameter, check=True): for Q, n in data: self._value *= phi(K(Q)) ** n self._value /= self._value[0] + self._value = self._value.list() if len(data.support()) == 1: self.normalize_point = Q + 1 # DEBUG elif data == 0: - self._value = Ps(1) # multiplicative! + self._value = Ps(1).list() # multiplicative! elif data.parent() == parent: - self._value = deepcopy(data._value) + self._value = data._value #deepcopy(data._value) else: val = Ps(data) val /= val[0] - self._value = val.add_bigoh(prec) + self._value = val.add_bigoh(prec).list() else: - self._value = deepcopy(data) + self._value = data # deepcopy(data) self._moments = self._value self._parameter = Matrix(parent.base_ring(), 2, 2, parameter) @@ -78,58 +88,85 @@ def __call__(self, D): def evaluate(self, D): # meromorphic functions K = self.parent().base_ring() phi = eval_pseries_map(self._parameter) - poly = self._value.polynomial() a, b, c, d = self._parameter.list() try: pt = phi(self.normalize_point) pole = -d / c - valinf = poly(pt) * (self.normalize_point - pole) + valinf = evalpoly(self._value,pt) * (self.normalize_point - pole) except AttributeError: pt = a / c pole = None - valinf = poly(pt) + valinf = evalpoly(self._value,pt) if isinstance(D.parent(), Divisors): return prod( - (poly(phi(P)) / valinf * ((P - pole) if pole is not None else 1)) ** n + (evalpoly(self._value,phi(P)) / valinf * ((P - pole) if pole is not None else 1)) ** n for P, n in D ) else: - return poly(phi(D)) / valinf * ((D - pole) if pole is not None else 1) + return evalpoly(self._value,phi(D)) / valinf * ((D - pole) if pole is not None else 1) + + def eval_derivative(self, D): + if isinstance(D.parent(), Divisors): + raise NotImplementedError + K = self.parent().base_ring() + phi = eval_pseries_map(self._parameter) + a, b, c, d = self._parameter.list() + valder = self.power_series().derivative().list() + try: + pt = phi(self.normalize_point) + pole = -d / c + valinf = evalpoly(self._value, pt) * (self.normalize_point - pole) + except AttributeError: + pt = a / c + pole = None + valinf = evalpoly(self._value, pt) + chainrule = (a*d-b*c) / (c*D+d) + return evalpoly(valder,phi(D)) * chainrule / valinf * ((D - pole) if pole is not None else 1) def _cmp_(self, right): - return (self._value > right._value) - (self._value < right._value) + a = self.parent().power_series_ring()(self._value) + b = self.parent().power_series_ring()(right._value) + return (a > b) - (a < b) def __nonzero__(self): - return self._value != 1 + return not (self._value[0] == 1 and all(o == 0 for o in self._value[1:])) def valuation(self, p=None): K = self.parent().base_ring() - return min( - [Infinity] + [o.valuation(p) for o in (self._value - 1).coefficients()] - ) + a = (self._value[0]-1).valuation(p) + b = min([o.valuation(p) for o in self._value[1:]]) + return min(a,b) def _add_(self, right): # multiplicative! - ans = self._value * right._value if self._parameter != right._parameter: raise RuntimeError("Trying to add incompatible series") + ps = self.parent().power_series_ring() + prec = self.parent()._prec + ans = (ps(self._value) * ps(right._value)).list()[:prec+1] return self.__class__(self.parent(), ans, self._parameter, check=False) def _sub_(self, right): # multiplicative! - ans = self._value / right._value if self._parameter != right._parameter: raise RuntimeError("Trying to subtract incompatible series") + ps = self.parent().power_series_ring() + prec = self.parent()._prec + ans = (ps(self._value) / ps(right._value)).list()[:prec+1] return self.__class__(self.parent(), ans, self._parameter, check=False) def _neg_(self): # multiplicative! - ans = self._value**-1 + ps = self.parent().power_series_ring() + prec = self.parent()._prec + ans = (ps(self._value)**-1).list()[:prec+1] return self.__class__(self.parent(), ans, self._parameter, check=False) def scale_by(self, k): # multiplicative! - ans = self._value**k + ps = self.parent().power_series_ring() + ans = (ps(self._value)**k).list()[:prec+1] return self.__class__(self.parent(), ans, self._parameter, check=False) def _repr_(self): - return repr(self._value) + ps = self.parent().power_series_ring() + return repr(ps(self._value)) def _acted_upon_(self, g, on_left): assert not on_left @@ -139,7 +176,7 @@ def _acted_upon_(self, g, on_left): return self.left_act_by_matrix(g) def left_act_by_matrix(self, g, param=None): # meromorphic functions - Ps = self._value.parent() + Ps = self.parent().power_series_ring() K = self.parent().base_ring() prec = self.parent()._prec p = self.parent()._p @@ -147,13 +184,18 @@ def left_act_by_matrix(self, g, param=None): # meromorphic functions # P |-> (aP+b)/(cP+D) # radius = tuple((ky, val) for ky, val in self.parent().radius.items()) zz_ps_vec = self.parent().get_action_data(g, self._parameter, param) - poly = self._value.polynomial() - ans = sum( - a * zz_ps_vec[e] for a, e in zip(poly.coefficients(), poly.exponents()) - ) - ans /= ans[0] # always normalize - ans = ans.add_bigoh(prec) - return MeromorphicFunctions(K, p=p, prec=prec)(ans, param) + poly = Ps(self._value).polynomial() + # print([type(zz_ps_vec[e]) for e in poly.exponents()]) + # ans = sum( + # a * Ps(zz_ps_vec[i]) for i, a in enumerate(self._value[:prec+1]) + # ).list() + # print(f'{prec = }') + # print(len(self._value)) + # print([len(o) for o in zz_ps_vec]) + ans = [sum(a * zz_ps_vec[i][j] for i, a in enumerate(self._value[:prec+1]) if j < len(zz_ps_vec[i])) for j in range(prec+1)] + r = ans[0] + ans = [o / r for o in ans[:prec+1]] # normalize + return MeromorphicFunctions(K, p=p, prec=prec)(ans, param, check=False) def eval_pseries_map(parameter): @@ -211,12 +253,12 @@ def get_action_data(self, g, oldparam, param, prec=None): a, b, c, d = (oldparam * (tg * g).adjugate()).list() zz = (self._Ps([b, a]) / self._Ps([d, c])).truncate(prec) Ps = self._Ps - ans = [Ps(1), zz] + ans = [zz.parent()(1), zz] for _ in range(prec - 1): zz_ps = (zz * ans[-1]).truncate(prec + 1) ans.append(zz_ps) set_verbose(verb_level) - return ans + return [o.list() for o in ans] def base_ring(self): return self._base_ring diff --git a/darmonpoints/schottky.py b/darmonpoints/schottky.py index 70401b4..98ed7ad 100644 --- a/darmonpoints/schottky.py +++ b/darmonpoints/schottky.py @@ -6,6 +6,7 @@ from sage.groups.matrix_gps.linear import GL from sage.matrix.matrix_space import MatrixSpace from sage.misc.cachefunc import cached_method +from sage.misc.latex import latex from sage.misc.verbose import verbose from sage.modules.module import Module from sage.rings.all import ZZ, IntegerRing @@ -17,56 +18,14 @@ from sage.structure.richcmp import richcmp from sage.structure.sage_object import SageObject from sage.structure.unique_representation import UniqueRepresentation -from sage.misc.latex import latex + from .divisors import Divisors, DivisorsElement from .meromorphic import * +from .theta import * infinity = Infinity -@cached_function -def find_eigenvector_matrix(g): - vaps = sorted([o for o, _ in g.charpoly().roots()], key=lambda o: o.valuation()) - verb_level = get_verbose() - set_verbose(0) - veps = sum(((g - l).right_kernel().basis() for l in vaps), []) - set_verbose(verb_level) - return g.matrix_space()(veps).transpose() - - -@cached_function -def find_parameter(g, r, pi=None, ball=None): - if pi is None: - pi = g.parent().base_ring().uniformizer() - if ball is not None: - assert r == ball.radius - ans = copy(find_eigenvector_matrix(g).adjugate()) - # column 0 means that the boundary of B_g (the ball attached to g) - # gets sent to the circle of radius 1. - # Concretely, if B_g = B(a, p^-r), we send a+p^r |-> 1. - # As a consequence, P^1-B_g gets sent to the closed ball B(0,1). - z0 = ball.center.lift_to_precision() + pi ** ZZ(ball.radius) - assert z0 in ball.closure() and z0 not in ball - lam = (ans[0, 0] * z0 + ans[0, 1]) / (ans[1, 0] * z0 + ans[1, 1]) - ans.rescale_row(0, 1 / lam) - ans.set_immutable() - assert act(ans, z0) == 1 - assert ans * ball.complement() == ball.parent()(0, 0).closure() - return ans - - -def act(g, z): - r""" - Compute the action of a 2 by 2 matrix g on an element z - """ - if g[1][0] == 0: - return (g[0][0] * z + g[0][1]) / g[1][1] - else: - return g[0][0] / g[1][0] + (g[0][1] - g[0][0] * g[1][1] / g[1][0]) / ( - g[1][0] * z + g[1][1] - ) - - def find_midpoint(K, v0, v1): t = K(v0.a - v1.a).valuation() if t < v0.r_valuation and t < v1.r_valuation: @@ -130,169 +89,6 @@ def uniq(lst): ans.append(o) return ans - -class ThetaOC(SageObject): - def __init__(self, G, a=0, b=None, prec=None, **kwargs): - K = G.K - self.G = G - self.p = G.pi - generators = G.generators() - balls = G.balls - self.K = K - if prec is None: - try: - self.prec = K.precision_cap() - except AttributeError: - self.prec = None - else: - self.prec = prec - self.Div = Divisors(K) - gens = [(i + 1, o) for i, o in enumerate(generators)] - gens.extend([(-i, o**-1) for i, o in gens]) - self.gens = tuple( - (i, o, find_parameter(o, balls[i].radius, self.p, ball=balls[i])) - for i, o in gens - ) - for i, o, t in self.gens: - o.set_immutable() - t.set_immutable() - if b is None: - D = self.Div(a) - if D.degree() != 0: - raise ValueError( - "Must specify a degree-0 divisor, or parameters a and b" - ) - else: - D = self.Div(K(a)) - self.Div(K(b)) - # if not all(G.in_fundamental_domain(P, closure=True) for P, n in D): - # raise ValueError("Divisor should be supported on fundamental domain.") - self.a = a - self.b = b - self.D = D - self.m = 1 - PP = PolynomialRing(K, names="z") - self.z = PP.gen() - - # self.val will contain the 0 and 1 terms - D1 = sum((g * D for i, g, tau in self.gens), self.Div([])) - self.val = prod(((self.z - P) ** n for P, n in (D + D1)), PP(1)) - Fdom = [(i,) for i, g, tau in self.gens] - self.Fn0 = {} - self.radius = [] - self.parameters = {} - D1dict = {} - for wd in Fdom: - g0 = G.element_to_matrix((wd[-1],)) - g1 = G.element_to_matrix(tuple(wd[:-1])) - assert g1 == 1 - g = g1 * g0 - g.set_immutable() - r = (g1 * balls[wd[-1]]).radius - self.radius.append((g, r)) - tau = find_parameter(g, r, self.p, ball=balls[wd[-1]]) - self.parameters[wd] = tau - D1dict[wd] = g * D - self.radius = tuple(self.radius) - for wd, (g, r) in zip(Fdom, self.radius): - gD = sum( - g * val for ky, val in D1dict.items() if ky[0] != -wd[0] - ) # DEBUG - only works when wd is a 1-tuple! - tau = self.parameters[wd] - self.Fn0[wd] = MeromorphicFunctions(self.K, self.p, self.prec)(gD, tau) - self.Fn = deepcopy(self.Fn0) - - def improve(self, m): - for it in range(m): - if self.m >= self.prec: - return self - newFn = deepcopy(self.Fn0) - for i, gi, _ in self.gens: - tau = self.parameters[(i,)] - for w, Fw in self.Fn.items(): - if i != -w[0]: - newFn[(i,)] += Fw.left_act_by_matrix(gi, tau) - self.Fn = newFn - self.m += 1 - return self - - def improve_one(self): - return self.improve(1) - - def _repr_(self): - a, b = self.a, self.b - if b is None: - try: - lst = a.as_list_of_differences() - if len(lst) == 1: - a, b = lst[0] - except AttributeError: - pass - try: - a = a.lift() - b = b.lift() - except AttributeError: - pass - return f"Θ(z;{a},{b})_{{{self.m}}}" - - def _latex_(self): - a, b = self.a, self.b - if b is None: - try: - lst = a.as_list_of_differences() - if len(lst) == 1: - a, b = lst[0] - except AttributeError: - pass - try: - a = a.lift() - b = b.lift() - except AttributeError: - pass - try: - a = a.rational_reconstruction() - b = b.rational_reconstruction() - except AttributeError: - pass - return f"\\Theta(z;{latex(a)},{latex(b)})_{{{latex(self.m)}}}" - - def __call__(self, z, recursive=True): - return self.evaluate(z, recursive=recursive) - - def evaluate(self, z, recursive=True): - if isinstance(z, DivisorsElement): - return prod(self(P, recursive=recursive) ** n for P, n in z) - G = self.G - ans = self.val(z) - ans *= prod(F(z) for ky, F in self.Fn.items()) - if recursive: - z0, wd = G.to_fundamental_domain(z) - wdab = [[g, 0] for g in G.generators()] - for i in wd: - wdab[abs(i) - 1] += sgn(i) - ans *= prod( - G.u_function(g, self.prec).evaluate(self.D, recursive=False) ** i - for g, i in wdab - ) - return ans - - def eval_derivative(self, z, recursive=True): - if recursive and not G.in_fundamental_domain(z, closure=True): - raise NotImplementedError("Recursivity not implemented for derivative") - if isinstance(z, DivisorsElement): - return prod(self.eval_derivative(P, recursive=recursive) ** n for P, n in z) - v0 = self.val(z) - Fnz = {} - for ky, F in self.Fn.items(): - Fnz[ky] = F(z) - ans = prod((F(z) for F in self.Fn.values()), self.val.derivative()(z)) - for ky0 in self.Fn: - ans *= prod( - (Fnz[ky] for ky, F in Fnz.items() if ky != ky0), - self.Fn[ky0].eval_derivative(z), - ) - return ans - - class SchottkyGroup_abstract(SageObject): def __init__(self, K, generators): self.K = K @@ -328,7 +124,7 @@ def inverse_generators(self): return self._inverse_generators def _repr_(self): - return "Schottky group with %s generators"%len(self.generators()) + return "Schottky group with %s generators" % len(self.generators()) def enumerate_group_elements(self, length): return enumerate_group_elements( @@ -343,7 +139,7 @@ def all_elements_up_to_length(self, N): ] ) - def theta_naive(self, z, m, a=ZZ(0), b=ZZ(1), gamma=None, **kwargs): + def theta_naive(self, m, a=ZZ(0), b=ZZ(1), z=None, gamma=None, **kwargs): if z is Infinity: return 1 if gamma is not None: @@ -376,35 +172,37 @@ def theta_naive(self, z, m, a=ZZ(0), b=ZZ(1), gamma=None, **kwargs): break return num / den - def theta_derivative_naive(self, z, m, a=ZZ(0), b=ZZ(1), **kwargs): - max_length = kwargs.get("max_length", -1) - K = kwargs.get("base_field", z.parent()) - z = K(z) - num = K(1) - den = K(1) + def theta_derivative_naive(self, m, a, b, z=None, max_length=-1, base_field=None, s=None): + if s is not None: + A = self.theta_naive(m,a,b,z=z,max_length=max_length,base_field=base_field) + B = self.theta_naive(m,act(s,a),act(s,b),z=z,max_length=max_length,base_field=base_field) + Ap = self.theta_derivative_naive(m,a,b,z=z,max_length=max_length,base_field=base_field,s=None) + Bp = self.theta_derivative_naive(m,act(s,a),act(s,b),z=z,max_length=max_length,base_field=base_field,s=None) + return Ap*B + A * Bp - second_term = K(0) - old_ans = K(1) + num = 1 + den = 1 + + second_term = 0 + old_ans = 1 val = 0 for i in NN: - # verbose(f'{i = }') + verbose(f'{i = }') for _, g in self.enumerate_group_elements(i): - ga = K(act(g, a)) - gb = K(act(g, b)) + ga = act(g, a) + gb = act(g, b) new_num = z - ga new_den = z - gb num *= new_num den *= new_den - new_second_term = (ga - gb) / K(new_num * new_den) + new_second_term = (ga - gb) / (new_num * new_den) second_term += new_second_term new_ans = (num / den) * second_term - can_stop = i > 0 and val >= (new_ans / old_ans - 1).valuation() - if i == max_length or (max_length == -1 and can_stop): + if i == max_length or (i > 0 and val >= (new_ans / old_ans -1).valuation()): break old_ans = new_ans - val = (new_ans / old_ans - 1).valuation() - return new_ans.add_bigoh(val + new_ans.valuation()) - + val = (new_ans / old_ans -1).valuation() + return new_ans.add_bigoh(val+new_ans.valuation()) class PreSchottkyGroup(SchottkyGroup_abstract): def __init__(self, K, generators): @@ -596,7 +394,7 @@ def to_schottky(self): class SchottkyGroup(SchottkyGroup_abstract): - def __init__(self, K, generators, balls=None, keep_generators=True): + def __init__(self, K, generators, balls=None, transformation=None, keep_generators=True, check=False): if balls is None: G = PreSchottkyGroup(K, generators) gp = G.group() @@ -608,12 +406,19 @@ def __init__(self, K, generators, balls=None, keep_generators=True): idx = sgn(i) * (j0 + 1) self.balls[i] = good_balls[idx] self.balls[-i] = good_balls[-idx] + self._transformation = [w for w in gp.gens()] else: if keep_generators: raise ValueError("Generators are not in good position") else: + # print(good_gens) generators = [G.element_to_matrix(g) for g in good_gens] self.balls = good_balls + self._transformation = good_gens + else: + self.balls = balls + self._transformation = transformation + if check: test_fundamental_domain(generators, self.balls) super().__init__(K, generators) @@ -634,11 +439,18 @@ def a_point(self, in_interior=True): x = K.random_element() raise RuntimeError("Reached maximum iterations (100)") - def find_point(self, gamma, eps=1): - gens = [(i + 1, o) for i, o in enumerate(self.generators())] - gens.extend([(-i, o.determinant() ** -1 * o.adjugate()) for i, o in gens]) + def find_point(self, gamma, eps=1, idx=None): balls = self.balls - B = next(balls[-i] for i, g in gens if g == gamma) + if idx is not None: + B = balls[-idx] + if idx > 0: + assert self.generators()[idx-1] == gamma + else: + assert self.generators()[-idx-1]**-1 == gamma + else: + gens = [(i + 1, o) for i, o in enumerate(self.generators())] + gens.extend([(-i, o.determinant() ** -1 * o.adjugate()) for i, o in gens]) + B = next(balls[-i] for i, g in gens if g == gamma) return B.center.lift_to_precision() + eps * self.pi ** ZZ(B.radius) def to_fundamental_domain(self, x): @@ -663,21 +475,25 @@ def to_fundamental_domain(self, x): i = self.in_which_ball(x) return x, word - def in_fundamental_domain(self, x, closure=False): - y, wd = self.to_fundamental_domain(x) - if len(wd) > 0: - return False - if closure: - return True - else: - return self.in_which_ball(y) is None + def in_fundamental_domain(self, x, strict=False): + return self.in_which_ball(x, closure=strict) is None def word_problem(self, gamma): z0 = self.a_point() z1 = act(gamma, z0) z2, word = self.to_fundamental_domain(z1) - if z0 != z2: - raise RuntimeError(f"Something went wrong! {z0 = }, {z2 = } {word = }") + if not z0.is_equal_to(z2, 10): # DEBUG + gens = [(i + 1, o) for i, o in enumerate(self.generators())] + gens.extend([(-i, o.determinant() ** -1 * o.adjugate()) for i, o in gens]) + found = False + for i, g in gens: + gz2 = act(g, z2) + if z0.is_equal_to(gz2, 10): + word.append(i) + found = True + break + if not found: + raise RuntimeError(f"Something went wrong! {z0 = }, {z2 = } {word = }") return tuple(word) def matrix_to_element(self, g): @@ -703,16 +519,16 @@ def find_equivalent_divisor(self, D): P0, new_wd = self.to_fundamental_domain(P) for i in new_wd: wd[abs(i) - 1] += n * sgn(i) - ans += n * Div(P0) + ans += n * Div([(1,P0)]) for i, (g, e) in enumerate(zip(gens, wd)): if e == 0: continue gamma = gens[i] zz = self.find_point(gamma) - ans -= e * (Div(zz) - Div(act(gamma, zz))) + ans -= e * Div([(1,zz), (-1, act(gamma,zz))]) return ans - def theta(self, prec, a=ZZ(0), b=ZZ(1), **kwargs): + def theta(self, prec, a, b=None, **kwargs): r""" Compute the Theta function @@ -730,39 +546,89 @@ def theta(self, prec, a=ZZ(0), b=ZZ(1), **kwargs): sage: b = 14 sage: z0 = K(8) sage: m = 10 - sage: Tg = G.theta_naive(z0, m, a=a,b=b) + sage: Tg = G.theta_naive(m, z=z0, a=a,b=b) sage: T = G.theta(prec, a, b).improve(m) sage: (T(z0) / Tg - 1).valuation() > m True """ - K = self.base_ring() - DK = Divisors(K) gens = self.generators() - D0 = DK(a) - DK(b) + if b is not None: + try: + K = a.parent() + except AttributeError: + K = self.base_ring() + if K.is_exact(): # DEBUG, set to avoid 0 and 1... + K = self.base_ring() + DK = Divisors(K) + D0 = DK([(1,a),(-1,b)]) + else: + D0 = a + try: + K = a.parent().base_ring() + except AttributeError: + K = self.base_ring() + if K.is_exact(): # DEBUG, set to avoid 0 and 1... + K = self.base_ring() + s = kwargs.pop("s", None) if s is not None: D0 += s * D0 D = self.find_equivalent_divisor(D0) - ans = ThetaOC(self, a=D, b=None, prec=prec, **kwargs) + ans = ThetaOC(self, a=D, b=None, prec=prec, base_ring=K, **kwargs) + z = kwargs.pop('z', None) improve = kwargs.pop("improve", True) - if improve: + if improve or z is not None: ans = ans.improve(prec) - return ans + return ans if z is None else ans(z) @cached_method def u_function(self, gamma, prec, a=None, **kwargs): r""" Compute u_gamma """ - K = self.base_ring() - DK = Divisors(K) + wd = self.word_problem(gamma) + gens = self.generators() + if len(wd) > 1: + if kwargs.get('z', None) is None: + raise NotImplementedError('Need to specify a point at which we will evaluate') + return prod(self.u_function(gens[abs(i)-1], prec, a=a, **kwargs)**sgn(i) for i in wd) + if a is None: - a = self.find_point(gamma) - assert self.in_fundamental_domain(a, closure=True) - assert self.in_fundamental_domain(act(gamma, a), closure=True) - D = DK(K(a)) - DK(K(act(gamma, a))) - return ThetaOC(self, a=D, b=None, prec=prec, **kwargs).improve(prec) + a = self.find_point(gamma) # Should pass idx=wd[0] DEBUG ?? + assert self.in_fundamental_domain(a, strict=False) + assert self.in_fundamental_domain(act(gamma, a), strict=False) + K = a.parent() + DK = Divisors(K) + D = DK([(1,a), (-1, act(gamma, a))]) + ans = ThetaOC(self, a=D, b=None, prec=prec, base_ring=K, **kwargs) + improve = kwargs.pop("improve", True) + if improve: + ans = ans.improve(prec) + z = kwargs.pop('z', None) + return ans if z is None else ans(z) + + @cached_method + def period_matrix(self, prec, **kwargs): + g = len(self.generators()) + M = MatrixSpace(self.base_ring(),g,g)(0) + z1 = kwargs.get("z1", None) + if z1 is None: + z1 = self.a_point() + for i in range(g): + g1 = self.generators()[i] + T = self.u_function(g1, prec, a=z1, **kwargs).improve(prec) + for j in range(i, g): + g2 = self.generators()[j] + z2 = kwargs.get("z2", None) + if z2 is None: + z2 = self.find_point(g2, eps=1 + self.pi) + g2_z2 = act(g2, z2) + num = T(z2, recursive=False) + den = T(g2_z2, recursive=False) + M[i,j] = num / den + M[j,i] = num / den + return M @cached_method def period(self, i, j, prec, **kwargs): @@ -792,8 +658,17 @@ def period(self, i, j, prec, **kwargs): sage: (q11g/q11-1).valuation() > prec True """ - g1 = self.generators()[i] - g2 = self.generators()[j] + g = len(self.generators()) + if i in ZZ: + assert j in ZZ + g1 = self.generators()[i] + g2 = self.generators()[j] + else: + a = word_to_abelian(i, g) + b = word_to_abelian(j, g) + m = self.period_matrix(prec, **kwargs) + return a * m * b + z1 = kwargs.pop("z1", None) if z1 is None: z1 = self.a_point() @@ -801,7 +676,7 @@ def period(self, i, j, prec, **kwargs): if z2 is None: z2 = self.find_point(g2, eps=1 + self.pi) g2_z2 = act(g2, z2) - T = self.u_function(g1, prec, a=z1, **kwargs) + T = self.u_function(g1, prec, a=z1, **kwargs).improve(prec) num = T(z2, recursive=False) den = T(g2_z2, recursive=False) verbose(f"{num = }") @@ -817,8 +692,8 @@ def period_naive(self, i, j, prec, **kwargs): z2 = self.find_point(g2, eps=self.pi + 1) else: z2 = self.find_point(g2) - num = self.theta_naive(z2, prec, a=z1, gamma=g1, **kwargs) - den = self.theta_naive(act(g2, z2), prec, a=z1, gamma=g1, **kwargs) + num = self.theta_naive(prec, z=z2, a=z1, gamma=g1, **kwargs) + den = self.theta_naive(prec, z=act(g2, z2), a=z1, gamma=g1, **kwargs) verbose(f"{num = }") verbose(f"{den = }") return num / den @@ -847,7 +722,6 @@ def __contains__(self, b): return b not in self.complement() if b is Infinity: return False - b = self.parent().K(b) d = b - self.center if self.is_open: in_ball = d.valuation() > self.radius diff --git a/darmonpoints/theta.py b/darmonpoints/theta.py new file mode 100644 index 0000000..a82aa6e --- /dev/null +++ b/darmonpoints/theta.py @@ -0,0 +1,258 @@ +from copy import copy, deepcopy +from itertools import chain + +from sage.categories.groups import Groups +from sage.functions.generalized import sgn +from sage.groups.matrix_gps.linear import GL +from sage.matrix.matrix_space import MatrixSpace +from sage.misc.cachefunc import cached_method +from sage.misc.latex import latex +from sage.misc.verbose import verbose +from sage.modules.module import Module +from sage.rings.all import ZZ, IntegerRing +from sage.rings.infinity import Infinity +from sage.rings.semirings.non_negative_integer_semiring import NN +from sage.sets.set import Set +from sage.structure.element import Element, ModuleElement +from sage.structure.parent import Parent +from sage.structure.richcmp import richcmp +from sage.structure.sage_object import SageObject +from sage.structure.unique_representation import UniqueRepresentation + +from .divisors import Divisors, DivisorsElement +from .meromorphic import * + +infinity = Infinity + +def act(mtrx,z): + """ + Compute the action of a 2x2 matrix on an element. + + If the base field is Qp, we first "normalize" mtrx to avoid + precision loss. + """ + if mtrx == 1: + return z + br = mtrx.base_ring() + if br.is_subring(QQ): + a, b, c, d = mtrx.list() + else: + pp = br.prime() + min_val = min([o.valuation() for o in mtrx.list()]) + new_mtrx = pp**(-min_val) * mtrx + a, b, c, d = new_mtrx.list() + + if z is Infinity: + return Infinity if c == 0 else a / c + try: + ans = (a*z+b)/d if c == 0 else a/c + (b-a*d/c)/(c*z+d) + except PrecisionError: + return Infinity + except TypeError: + emb = z.parent().embedding_from_subfield + a,b,c,d = emb(a), emb(b), emb(c), emb(d) + try: + ans = (a*z+b)/d if c == 0 else a/c + (b-a*d/c)/(c*z+d) + except PrecisionError: + return Infinity + return ans + +@cached_function +def find_eigenvector_matrix(g): + vaps = sorted([o for o, _ in g.charpoly().roots()], key=lambda o: o.valuation()) + verb_level = get_verbose() + set_verbose(0) + veps = sum(((g - l).right_kernel().basis() for l in vaps), []) + set_verbose(verb_level) + return g.matrix_space()(veps).transpose() + +@cached_function +def find_parameter(g, r, pi=None, ball=None, check=True): + if pi is None: + pi = g.parent().base_ring().uniformizer() + if ball is not None: + assert r == ball.radius + ans = copy(find_eigenvector_matrix(g).adjugate()) + # column 0 means that the boundary of B_g (the ball attached to g) + # gets sent to the circle of radius 1. + # Concretely, if B_g = B(a, p^-r), we send a+p^r |-> 1. + # As a consequence, P^1-B_g gets sent to the closed ball B(0,1). + z0 = ball.center.lift_to_precision() + pi ** ZZ(ball.radius) + if check: + assert z0 in ball.closure() and z0 not in ball + lam = (ans[0, 0] * z0 + ans[0, 1]) / (ans[1, 0] * z0 + ans[1, 1]) + ans.rescale_row(0, 1 / lam) + ans.set_immutable() + if check: + if act(ans, z0) != 1: + raise RuntimeError + if ans * ball.complement() != ball.parent()(0, 0).closure(): + # print(ans * ball.complement()) + # print(ball.parent()(0,0).closure()) + raise RuntimeError + return ans + +def element_to_matrix(wd, generators): + ans = wd(generators) + ans.set_immutable() + return ans + +def word_to_abelian(wd, g): + try: + wd = wd.Tietze() + except AttributeError: + pass + wdab = [0 for i in range(g)] + for i in wd: + wdab[abs(i) - 1] += sgn(i) + return (ZZ**g)(wdab) + +class ThetaOC(SageObject): + def __init__(self, G, a=0, b=None, prec=None, **kwargs): + K = kwargs.get('base_ring',None) + if K is None: + K = G.K + self.K = K + self.G = G + self.p = G.pi + generators = G.generators() + balls = G.balls + if prec is None: + try: + self.prec = K.precision_cap() + except AttributeError: + self.prec = None + else: + self.prec = prec + self.Div = Divisors(K) + gens = [(i + 1, o) for i, o in enumerate(generators)] + gens.extend([(-i, o**-1) for i, o in gens]) + self.gens = tuple( + (i, o, find_parameter(o, balls[i].radius, self.p, ball=balls[i])) + for i, o in gens + ) + for i, o, t in self.gens: + o.set_immutable() + t.set_immutable() + if b is None: + D = self.Div(a) + if D.degree() != 0: + raise ValueError( + "Must specify a degree-0 divisor, or parameters a and b" + ) + else: + D = self.Div([(1,a), (-1,b)]) + # if not all(G.in_fundamental_domain(P, strict=False) for P, n in D): + # raise ValueError("Divisor should be supported on fundamental domain.") + self.a = a + self.b = b + self.D = D + self.m = 1 + PP = PolynomialRing(K, names="z") + self.z = PP.gen() + + # self.val will contain the 0 and 1 terms + D1 = sum((g * D for i, g, tau in self.gens), self.Div([])) + self.val = prod(((self.z - P) ** n for P, n in (D + D1)), PP(1)) + self.Fn0 = {} + D1dict = {} + for i, g, tau in self.gens: + D1dict[i] = g * D + for i, g, tau in self.gens: + gD = sum( + g * val for ky, val in D1dict.items() if ky != -i + ) # DEBUG - only works when wd is a 1-tuple! + self.Fn0[i] = MeromorphicFunctions(self.K, self.p, self.prec)(gD, tau) + self.Fn = self.Fn0 + + def improve(self, m): + for it in range(m): + if self.m >= self.prec: + return self + newFn = deepcopy(self.Fn0) + for i, gi, tau in self.gens: + for j, Fj in self.Fn.items(): + if i != -j: + newFn[i] += Fj.left_act_by_matrix(gi, tau) + self.Fn = newFn + self.m += 1 + return self + + def improve_one(self): + return self.improve(1) + + def _repr_(self): + a, b = self.a, self.b + if b is None: + try: + lst = a.as_list_of_differences() + if len(lst) == 1: + a, b = lst[0] + except AttributeError: + pass + try: + a = a.lift() + b = b.lift() + except AttributeError: + pass + return f"Θ(z;{a},{b})_{{{self.m}}}" + + def _latex_(self): + a, b = self.a, self.b + if b is None: + try: + lst = a.as_list_of_differences() + if len(lst) == 1: + a, b = lst[0] + except AttributeError: + pass + try: + a = a.lift() + b = b.lift() + except AttributeError: + pass + try: + a = a.rational_reconstruction() + b = b.rational_reconstruction() + except AttributeError: + pass + return f"\\Theta(z;{latex(a)},{latex(b)})_{{{latex(self.m)}}}" + + def __call__(self, z, **kwargs): + return self.evaluate(z, **kwargs) + + def evaluate(self, z, **kwargs): + if isinstance(z, DivisorsElement): + return prod(self(P, **kwargs) ** n for P, n in z) + G = self.G + recursive = kwargs.get('recursive', False) + ans = 1 + if recursive: + z0, wd = G.to_fundamental_domain(z) + wdab = word_to_abelian(wd, len(G.generators())) + ans *= prod( + G.u_function(g, self.prec).evaluate(self.D, recursive=False) ** i + for g, i in zip(G.generators(), wdab) if i != 0 + ) + else: + z0 = z + ans *= self.val(z0) + ans *= prod(F(z0) for ky, F in self.Fn.items()) + return ans + + def eval_derivative(self, z, recursive=False): + if recursive and not G.in_fundamental_domain(z, strict=False): + raise NotImplementedError("Recursivity not implemented for derivative") + if isinstance(z, DivisorsElement): + return prod(self.eval_derivative(P, recursive=recursive) ** n for P, n in z) + v0 = self.val(z) + Fnz = {} + for ky, F in self.Fn.items(): + Fnz[ky] = F(z) + ans = prod((o for o in Fnz.values()), self.val.derivative()(z)) + for ky0 in self.Fn: + ans += v0 * prod( + (Fnz[ky] for ky, F in Fnz.items() if ky != ky0), + self.Fn[ky0].eval_derivative(z), + ) + return ans