diff --git a/MANIFEST.in b/MANIFEST.in index a2857ec..517630b 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,3 +1,4 @@ include darmonpoints/*.pyx -include darmonpoints/Fdoms/* -include darmonpoints/KleinianGroups-1.0/* +recursive-include darmonpoints/Fdoms * +recursive-include darmonpoints/KleinianGroups-1.0 * + diff --git a/darmonpoints/arithgroup_generic.py b/darmonpoints/arithgroup_generic.py index 2e2f198..1d7efbe 100644 --- a/darmonpoints/arithgroup_generic.py +++ b/darmonpoints/arithgroup_generic.py @@ -704,8 +704,8 @@ def cusp_reduction_table(self): ## Define new function on the fly to pick which of Q/more general field we work in ## lift_to_matrix takes parameters c,d, then lifts (c:d) to a 2X2 matrix over the NF representing it - lift_to_matrix = ( - lambda c, d: lift_to_sl2z(c, d, P.N()) + lift_to_matrix = lambda c, d: ( + lift_to_sl2z(c, d, P.N()) if K.degree() == 1 else lift_to_sl2_Ok(P.N(), c, d) ) diff --git a/darmonpoints/cohomology_arithmetic.py b/darmonpoints/cohomology_arithmetic.py index d37a625..b725165 100644 --- a/darmonpoints/cohomology_arithmetic.py +++ b/darmonpoints/cohomology_arithmetic.py @@ -527,6 +527,7 @@ class ArithCoh(CohomologyGroup, UniqueRepresentation): Initialised by inputting an arithmetic group G, and the coefficient module. """ + Element = ArithCohElement def __init__(self, G, V=None, **kwargs): diff --git a/darmonpoints/divisors.py b/darmonpoints/divisors.py index fd03727..a622541 100644 --- a/darmonpoints/divisors.py +++ b/darmonpoints/divisors.py @@ -112,6 +112,9 @@ def restrict(self, condition): ] ) + def intersects(self, right): + return any(hP in right._data for hP in self._data) + def __iter__(self): return iter(((self._ptdict[hP], n) for hP, n in self._data.items())) @@ -203,17 +206,19 @@ def _acted_upon_(self, g, on_left): def left_act_by_matrix(self, g): a, b, c, d = g.list() gp = self.parent() - K = self.parent().base_ring() newdict = defaultdict(ZZ) new_ptdata = {} for P, n in self: if P == Infinity: try: - new_pt = K(a) / K(c) + new_pt = a / c except (PrecisionError, ZeroDivisionError): new_pt = Infinity else: - new_pt = (K(a) * P + K(b)) / (K(c) * P + K(d)) + try: + new_pt = (a * P + b) / (c * P + d) + except (PrecisionError, ZeroDivisionError): + new_pt = Infinity hnew_pt = _hash(new_pt) newdict[hnew_pt] += n new_ptdata[hnew_pt] = new_pt @@ -254,7 +259,8 @@ def rational_function(self, as_map=False, z=None): return lambda z: prod(((1 - z / P) ** n for P, n in self), z.parent()(1)) else: if z is None: - z = PolynomialRing(self.parent()._field, names="z").gen() + K = self.parent()._field + z = K["z"].gen() return prod(((1 - z / P) ** n for P, n in self), z.parent()(1)) def as_list_of_differences(self): diff --git a/darmonpoints/meromorphic.py b/darmonpoints/meromorphic.py index 3684acc..e212cfb 100644 --- a/darmonpoints/meromorphic.py +++ b/darmonpoints/meromorphic.py @@ -39,9 +39,13 @@ from .divisors import Divisors + def evalpoly(poly, x): # Initialize result - result = poly[-1] + try: + result = poly[-1] + except IndexError: + return x.parent()(0) # Evaluate value of polynomial # using Horner's method for a in reversed(poly): @@ -49,13 +53,13 @@ def evalpoly(poly, x): result += a return result + class MeromorphicFunctionsElement(ModuleElement): def __init__(self, parent, data, parameter, check=True): - # verbose(f'{check = }') ModuleElement.__init__(self, parent) + prec = parent._prec if check: K = parent.base_ring() - prec = parent._prec Ps = parent._Ps if isinstance(data.parent(), Divisors): self._value = divisor_to_pseries(parameter, Ps, data, prec) @@ -73,57 +77,84 @@ def __init__(self, parent, data, parameter, check=True): else: self._value = data self._moments = self._value + self._value = [o.add_bigoh(prec) for o in self._value] + while len(self._value) < prec: + self._value.append(0) + r = self._value[0] + if r != 1: + rinv = ~r + self._value = [rinv * o for o in self._value] + self._value = vector(self._value) self._parameter = Matrix(parent.base_ring(), 2, 2, parameter) def power_series(self): - return self.parent().power_series_ring()(self._value) + return self.parent().power_series_ring()(list(self._value)) def __call__(self, D): return self.evaluate(D) def evaluate(self, D): # meromorphic functions K = self.parent().base_ring() - phi = eval_pseries_map(self._parameter) + if type(D) in (int, Integer): + D = K(D) a, b, c, d = self._parameter.list() + phi = lambda Q: a / c if Q == Infinity else (a * Q + b) / (c * Q + d) try: pt = phi(self.normalize_point) pole = -d / c - valinf = evalpoly(self._value,pt) * (self.normalize_point - pole) + valinf = evalpoly(self._value, pt) * (self.normalize_point - pole) except AttributeError: - pt = a / c pole = None - valinf = evalpoly(self._value,pt) + if c == 0: + valinf = 1 + else: + valinf = evalpoly(self._value, a / c) assert pole is None + + def ev(P): + fac = (P - pole) if pole is not None else 1 + phiP = phi(P) + return fac * evalpoly(self._value, phiP) / valinf + # except PrecisionError: + # return fac + if isinstance(D.parent(), Divisors): - return prod( - (evalpoly(self._value,phi(P)) / valinf * ((P - pole) if pole is not None else 1)) ** n - for P, n in D - ) + return prod(ev(P) ** n for P, n in D) else: - return evalpoly(self._value,phi(D)) / valinf * ((D - pole) if pole is not None else 1) + return ev(D) def eval_derivative(self, D): + K = self.parent().base_ring() + if type(D) in (int, Integer): + D = K(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() + phi = lambda Q: a / c if Q == Infinity else (a * Q + b) / (c * Q + d) 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) + if c == 0: + valinf = 1 + else: + valinf = evalpoly(self._value, a / c) assert pole is None - chainrule = (a*d-b*c) / (c*D+d)**2 - return evalpoly(valder,phi(D)) * chainrule / valinf * ((D - pole) if pole is not None else 1) + chainrule = (a * d - b * c) / (c * D + d) ** 2 + return ( + evalpoly(valder, phi(D)) + * chainrule + / valinf + * ((D - pole) if pole is not None else 1) + ) def _cmp_(self, right): - a = self.parent().power_series_ring()(self._value) - b = self.parent().power_series_ring()(right._value) + a = self.power_series() + b = right.power_series() return (a > b) - (a < b) def __nonzero__(self): @@ -131,40 +162,35 @@ def __nonzero__(self): def valuation(self, p=None): K = self.parent().base_ring() - a = (self._value[0]-1).valuation(p) + a = (self._value[0] - 1).valuation(p) b = min([o.valuation(p) for o in self._value[1:]]) - return min(a,b) + return min(a, b) def _add_(self, right): # multiplicative! 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] + ans = (self.power_series() * right.power_series()).list()[:prec] return self.__class__(self.parent(), ans, self._parameter, check=False) def _sub_(self, right): # multiplicative! 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] + ans = (self.power_series() / right.power_series()).list()[:prec] return self.__class__(self.parent(), ans, self._parameter, check=False) def _neg_(self): # multiplicative! - ps = self.parent().power_series_ring() prec = self.parent()._prec - ans = (ps(self._value)**-1).list()[:prec+1] + ans = (~self.power_series()).list()[:prec] return self.__class__(self.parent(), ans, self._parameter, check=False) def scale_by(self, k): # multiplicative! - ps = self.parent().power_series_ring() - ans = (ps(self._value)**k).list()[:prec+1] + ans = (self.power_series() ** k).list()[:prec] return self.__class__(self.parent(), ans, self._parameter, check=False) def _repr_(self): - ps = self.parent().power_series_ring() - return repr(ps(self._value)) + return repr(self.power_series()) def _acted_upon_(self, g, on_left): assert not on_left @@ -173,27 +199,23 @@ def _acted_upon_(self, g, on_left): else: return self.left_act_by_matrix(g) + def fast_act(self, key): + zz_ps_vec, param = key + return self.__class__( + self.parent(), self._value * zz_ps_vec, param, check=False + ) + def left_act_by_matrix(self, g, param=None): # meromorphic functions - Ps = self.parent().power_series_ring() - K = self.parent().base_ring() - prec = self.parent()._prec - p = self.parent()._p + t = cputime() + parent = self.parent() + K = parent.base_ring() + prec = parent._prec + p = parent._p # Below we code the action which is compatible with the usual action # 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 - 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).add_bigoh(prec) for o in ans[:prec+1]] # normalize - return MeromorphicFunctions(K, p=p, prec=prec)(ans, param, check=False) - - -def eval_pseries_map(parameter): - a, b, c, d = parameter.list() - return lambda Q: (Q.parent()(a) * Q + Q.parent()(b)) / ( - Q.parent()(c) * Q + Q.parent()(d) - ) + zz_ps_vec = parent.get_action_data(g, self._parameter, param) + ans = self._value * zz_ps_vec + return self.__class__(self.parent(), ans, param, check=False) def divisor_to_pseries(parameter, Ps, data, prec): @@ -204,9 +226,9 @@ def divisor_to_pseries(parameter, Ps, data, prec): for Q, n in data: K = Q.parent() if n > 0: - num *= (1 - K(((c*Q+d)/(a*Q+b)))*t)**ZZ(n) + num *= (1 - K(((c * Q + d) / (a * Q + b))) * t) ** ZZ(n) else: - den *= (1 - K(((c*Q+d)/(a*Q+b)))*t)**ZZ(-n) + den *= (1 - K(((c * Q + d) / (a * Q + b))) * t) ** ZZ(-n) ans = Ps(num).add_bigoh(prec) * ~(Ps(den).add_bigoh(prec)) return ans @@ -216,6 +238,7 @@ class MeromorphicFunctions(Parent, UniqueRepresentation): TESTS: """ + Element = MeromorphicFunctionsElement @staticmethod @@ -250,13 +273,25 @@ def get_action_data(self, g, oldparam, param, prec=None): assert param is not None tg = param.change_ring(K) a, b, c, d = (oldparam * (tg * g).adjugate()).list() - zz = (self._Ps([b, a]) / self._Ps([d, c])).map_coefficients(lambda x : x.add_bigoh(prec)).truncate(prec) + zz = ( + (self._Ps([b, a]) / self._Ps([d, c])) + .map_coefficients(lambda x: x.add_bigoh(prec)) + .truncate(prec) + ) ans = [zz.parent()(1), zz] - for _ in range(prec - 1): - zz_ps = (zz * ans[-1]).map_coefficients(lambda x : x.add_bigoh(prec)).truncate(prec + 1) + while len(ans) < prec: # DEBUG - was prec + 1 + zz_ps = ( + (zz * ans[-1]) + .map_coefficients(lambda x: x.add_bigoh(prec)) + .truncate(prec) # DEBUG - was prec + 1 + ) ans.append(zz_ps) set_verbose(verb_level) - return [o.list() for o in ans] + m = Matrix(K, prec, prec, 0) + for i, zz_ps in enumerate(ans): + for j, aij in enumerate(zz_ps): + m[i, j] = aij # i, j entry contains j-th term of zz^i + return m def base_ring(self): return self._base_ring diff --git a/darmonpoints/ocmodule.py b/darmonpoints/ocmodule.py index 056630d..4097304 100644 --- a/darmonpoints/ocmodule.py +++ b/darmonpoints/ocmodule.py @@ -690,6 +690,7 @@ class AddMeromorphicFunctions(Parent, CachedRepresentation): sage: A == B True """ + Element = AddMeromorphicFunctionsElement def __init__(self, K, twisting_matrix=None): diff --git a/darmonpoints/padicperiods.py b/darmonpoints/padicperiods.py index 9cd9759..e2f4ae8 100644 --- a/darmonpoints/padicperiods.py +++ b/darmonpoints/padicperiods.py @@ -4,12 +4,14 @@ from itertools import chain, groupby, islice, product, starmap, tee import pyximport +from sage.arith.srange import srange from sage.matrix.constructor import Matrix, block_diagonal_matrix, block_matrix, matrix from sage.misc.banner import version as sage_version -from sage.misc.persist import db +from sage.misc.persist import db, load, save from sage.modules.fg_pid.fgp_module import FGP_Module, FGP_Module_class from sage.rings.integer_ring import ZZ from sage.rings.padics.precision_error import PrecisionError +from sage.rings.power_series_ring import PowerSeriesRing from .cohomology_arithmetic import * from .homology import get_homology_kernel, inverse_shapiro, lattice_homology_cycle @@ -54,22 +56,22 @@ def Thetas(p1, p2, p3, q1, q2, q3, prec=None): c = 0 # Define the different conditions and term forms resdict = {} - resdict["1p1m"] = resdict["1p2m"] = 0 - resdict["2p2m"] = resdict["1m2p"] = 0 - resdict["3p3m"] = resdict["2m3p"] = 0 - resdict["2p3m"] = resdict["3m1p"] = resdict["3p1m"] = 0 - res = 0 + resdict["1p1m"] = resdict["1p2m"] = QQ(0) + resdict["2p2m"] = resdict["1m2p"] = QQ(0) + resdict["3p3m"] = resdict["2m3p"] = QQ(0) + resdict["2p3m"] = resdict["3m1p"] = resdict["3p1m"] = QQ(0) + res = QQ(0) assert imax > 0 jdict = {} p1dict = {} p3dict = {} p2dict = precompute_powers(p2, q2, imax) - for i in range(-imax, imax + 1): + for i in srange(-imax, imax + 1): newjmax = RR(2 * prec - 0.25 - i**2 + RR(i).abs()) if newjmax >= 0: newjmax = (newjmax.sqrt() + 0.5).ceiling() jdict[i] = newjmax - jrange = range(-newjmax, newjmax + 1) + jrange = srange(-newjmax, newjmax + 1) jlist = list( set([j**2 - j for j in jrange] + [j for j in jrange]).difference( set(p1dict.keys()) @@ -83,17 +85,17 @@ def Thetas(p1, p2, p3, q1, q2, q3, prec=None): ).difference(set(p3dict.keys())) ) for j in jlist: - tmp = q1 ** QQ(j / 2).floor() + tmp = q1 ** ZZ(QQ(j / 2).floor()) if j % 2 == 1: tmp *= p1 p1dict[j] = tmp for k in klist: - tmp = q3 ** QQ(k / 2).floor() + tmp = q3 ** ZZ(QQ(k / 2).floor()) if k % 2 == 1: tmp *= p3 p3dict[k] = tmp for i, jmax in jdict.items(): - for j in range(-jmax, jmax + 1): + for j in srange(-jmax, jmax + 1): P = p1dict[j**2 - j] * p2dict[i**2 - i] * p3dict[(i - j) ** 2 - (i - j)] p11 = p1dict[j] p22 = p2dict[i] @@ -214,13 +216,13 @@ def Theta(p1, p2, p3, version, prec=None): resdict["3p1m"] += p1lp3l_inv_term else: res += newterm * term(i, j) - return res + return resdict if version is None else res def lambdavec(p1, p2, p3, prec=None, theta=None, prec_pseries=None): if theta is None: assert prec is not None - th = Thetas(p1, p2, p3, prec=prec) + th = Thetas(p1, p2, p3, p1**2, p2**2, p3**3, prec=prec) else: th = theta @@ -285,8 +287,14 @@ def xvec_padic(p1, p2, p3, q1, q2, q3, prec=None): def igusa_clebsch_absolute_from_half_periods( - p1, p2, p3, q1, q2, q3, prec=None, padic=True + p1, p2, p3, q1=None, q2=None, q3=None, prec=None, padic=True ): + if q1 is None: + q1 = p1**2 + if q2 is None: + q2 = p2**2 + if q3 is None: + q3 = p3**2 I2, I4, I6, I10 = igusa_clebsch_from_half_periods( p1, p2, p3, q1, q2, q3, prec=prec, padic=padic ) @@ -1545,12 +1553,25 @@ def jacobian_matrix(fvec): def compute_lvec_and_Mlist(prec): - R = PowerSeriesRing(QQ, names="p", num_gens=3, default_prec=prec) - p1, p2, p3 = R.gens() - theta = Theta(p1, p2, p3, version=None, prec=prec) - lvec = lambdavec(p1, p2, p3, theta=theta, prec_pseries=prec) - Mlist = compute_twisted_jacobian_data(lvec) - lvec = [o.polynomial() for o in lvec.list()] + R0 = PolynomialRing(QQ, 3, names="p") + R = R0.fraction_field() + PS = PowerSeriesRing(QQ, names="p", num_gens=3) + p1, p2, p3 = R0.gens() + p1, p2, p3 = R(p1), R(p2), R(p3) + # theta = Theta(p1, p2, p3, version=None, prec=prec) + lvec = matrix( + PS, + 3, + 1, + [ + PS(o.numerator()) / PS(o.denominator()) + for o in lambdavec( + p1, p2, p3, prec=prec, theta=None, prec_pseries=prec + ).list() + ], + ) + p1, p2, p3 = PS.gens() + Mlist = compute_twisted_jacobian_data(lvec, p1, p2, p3, prec) save((lvec, Mlist), "lvec_and_Mlist_%s.sobj" % prec) return (lvec, Mlist) @@ -1558,29 +1579,29 @@ def compute_lvec_and_Mlist(prec): def load_lvec_and_Mlist(prec): try: lvec, Mlist = load("lvec_and_Mlist_%s.sobj" % prec) - except OSError: + except FileNotFoundError: lvec, Mlist = compute_lvec_and_Mlist(prec) return (lvec, Mlist) def evaluate_twisted_jacobian_matrix(p1, p2, p3, Mlist): - retvec = [Mlist[i](p1, p2, p3) for i in range(6)] - h1 = Mlist[6](p1, p2, p3) - h2 = Mlist[7](p1, p2, p3) - h3 = Mlist[8](p1, p2, p3) + mlist = [o.polynomial()(p1, p2, p3) for o in Mlist] + retvec = mlist[:6] + h1 = mlist[6] + h2 = mlist[7] + h3 = mlist[8] ll = ((1 - p3) / (1 + p3)) ** 2 h1 = ll * h1 h2 = ll * h2 - h3 = -4 * (1 - p3) / (1 + p3) ** 3 * Mlist[9](p1, p2, p3) + ll * h3 + h3 = -4 * (1 - p3) / (1 + p3) ** 3 * mlist[9] + ll * h3 retvec.extend([h1, h2, h3]) return Matrix(3, 3, retvec) -def compute_twisted_jacobian_data(fvec): +def compute_twisted_jacobian_data(fvec, x, y, z, prec): f1, f2, f3 = fvec.list() - x, y, z = f1.variables() Mlist = [ - o.polynomial() + o.truncate(prec) for o in [ f1.derivative(x), f1.derivative(y), @@ -1591,19 +1612,20 @@ def compute_twisted_jacobian_data(fvec): f3.derivative(x), f3.derivative(y), f3.derivative(z), + f3, ] - ] + [f3.polynomial()] + ] return Mlist def find_initial_approx(L1, L2, L3, lvec): # Evaluates a matrix M with entries in Z[[x,y,z]] at points x0,y0,z0 def ev(x0, y0, z0): - return [ + return ( lvec[0](x0, y0, z0), lvec[1](x0, y0, z0), ((1 - z0) / (1 + z0)) ** 2 * lvec[2](x0, y0, z0), - ] + ) K = L1.parent() n_tries = 0 @@ -1615,7 +1637,7 @@ def ev(x0, y0, z0): FP = ev(*P) if min([(u - v).valuation() for u, v in zip(FP, [L1, L2, L3])]) > 1 / 2: good_P = P - print(P, [(u - v).valuation() for u, v in zip(FP, [L1, L2, L3])]) + # print(P, [(u - v).valuation() for u, v in zip(FP, [L1, L2, L3])]) return good_P @@ -1631,6 +1653,11 @@ def HalfPeriodsInTermsOfLambdas( else: lvec, Mlist = lvec_and_Mlist + try: + lvec = [o.polynomial() for o in lvec.list()] + except AttributeError: + pass + # Evaluates a matrix M with entries in Z[[x,y,z]] at points x0,y0,z0 def ev(x0, y0, z0): return [ @@ -1643,14 +1670,15 @@ def ev(x0, y0, z0): HP0 = [K(0), K(0), K(0)] Pn = Matrix(K, 3, 1, HP0) # 0th approximation n_iters = 0 - while n_iters < 20: + while n_iters < max_iters: n_iters += 1 Jinv = evaluate_twisted_jacobian_matrix(*Pn.list(), Mlist=Mlist).inverse() FPn = matrix(3, 1, ev(*Pn.list())) Pnn = Pn - Jinv * (FPn - L0) print( - "(%s)" % n_iters, - [(u - v).valuation() for u, v in zip(Pn.list(), Pnn.list())], + "(%s)" + % n_iters + # [(u - v).valuation() for u, v in zip(Pn.list(), Pnn.list())], ) if all([u == v for u, v in zip(Pn.list(), Pnn.list())]): return Pn diff --git a/darmonpoints/rationalfunctions.py b/darmonpoints/rationalfunctions.py index b120842..7b062e9 100644 --- a/darmonpoints/rationalfunctions.py +++ b/darmonpoints/rationalfunctions.py @@ -155,6 +155,7 @@ def left_act_by_matrix(self, g): # rational functions class RationalFunctions(Parent, CachedRepresentation): r""" """ + Element = RationalFunctionsElement def __init__(self, K): diff --git a/darmonpoints/representations.py b/darmonpoints/representations.py index fac91f0..a0df41b 100644 --- a/darmonpoints/representations.py +++ b/darmonpoints/representations.py @@ -193,6 +193,7 @@ class CoIndModule(Parent): sage: from darmonpoints.sarithgroup import * sage: G = BigArithGroup(5,6,1,outfile='/tmp/darmonpoints.tmp') # optional - magma """ + Element = CoIndElement def __init__(self, G, V): diff --git a/darmonpoints/schottky.py b/darmonpoints/schottky.py index 3d2056b..e6d3ca5 100644 --- a/darmonpoints/schottky.py +++ b/darmonpoints/schottky.py @@ -2,11 +2,15 @@ from itertools import chain from sage.categories.groups import Groups +from sage.combinat.rooted_tree import LabelledRootedTree as LTR from sage.functions.generalized import sgn +from sage.graphs.graph import Graph from sage.groups.matrix_gps.linear import GL +from sage.matrix.constructor import Matrix from sage.matrix.matrix_space import MatrixSpace from sage.misc.cachefunc import cached_method from sage.misc.latex import latex +from sage.misc.lazy_attribute import lazy_attribute from sage.misc.verbose import verbose from sage.modules.module import Module from sage.rings.all import ZZ, IntegerRing @@ -18,26 +22,15 @@ from sage.structure.richcmp import richcmp from sage.structure.sage_object import SageObject from sage.structure.unique_representation import UniqueRepresentation -from sage.misc.lazy_attribute import lazy_attribute + from .divisors import Divisors, DivisorsElement from .meromorphic import * from .theta import * +from .util import muted infinity = Infinity -def find_midpoint(K, v0, v1): - t = K(v0.a - v1.a).valuation() - if t < v0.r_valuation and t < v1.r_valuation: - distance = (v0.r_valuation - t) + (v1.r_valuation - t) - else: - distance = abs(v0.r_valuation - v1.r_valuation) - v, comp = (v1, False) if v1.r_valuation > v0.r_valuation else (v0, True) - return Balls(K)( - v.a, v.r_valuation - QQ(distance) / 2, is_open=not comp, is_complement=comp - ) - - def reduce_word(w): r = [] for i in w: @@ -89,16 +82,15 @@ def uniq(lst): ans.append(o) return ans + +@muted @cached_function def find_eigenvector_matrix(g): t = g.trace() n = g.determinant() - delta = (t*t - 4*n).sqrt() - vaps = sorted([(t+delta)/2, (t-delta)/2], key=lambda o: o.valuation()) - verb_level = get_verbose() - set_verbose(0) + delta = (t * t - 4 * n).sqrt() + vaps = sorted([(t + delta) / 2, (t - delta) / 2], key=lambda o: o.valuation()) veps = sum(((g - l).right_kernel().basis() for l in vaps), []) - set_verbose(verb_level) return g.matrix_space()(veps).transpose() @@ -106,7 +98,7 @@ def find_eigenvector_matrix(g): def find_parameter(g, ball, pi=None, check=True): if pi is None: pi = g.parent().base_ring().uniformizer() - g = g.apply_map(lambda o:o.lift_to_precision()) + g = g.apply_map(lambda o: o.lift_to_precision()) ans = copy(find_eigenvector_matrix(g).adjugate()) r = ball.radius # column 0 means that the boundary of B_g (the ball attached to g) @@ -122,8 +114,25 @@ def find_parameter(g, ball, pi=None, check=True): if check: if act(ans, z0) != 1: raise RuntimeError - if ans * ball.complement() != ball.parent()(0, 0).closure(): - 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 + + +@cached_function +def find_parameter_new(g, ball, pi=None, check=True): + K = g.parent().base_ring() + if pi is None: + pi = K.uniformizer() + g = g.apply_map(lambda o: o.lift_to_precision()) + r = ball.radius + z0 = ball.center.lift_to_precision() + if ball.is_complement: + ans = Matrix(K, 2, 2, [0, 1, 1, -z0]) + else: + ans = Matrix(K, 2, 2, [1, -z0, 0, 1]) return ans @@ -132,6 +141,7 @@ def __init__(self, K, generators): self.K = K self.pi = K.uniformizer() self._generators = tuple([o.change_ring(K) for o in generators]) + self._inverse_generators = tuple(v.adjugate() for v in self._generators) self._G = Groups().free(len(generators)) def base_ring(self): @@ -143,9 +153,13 @@ def base_field(self): def uniformizer(self): return self.pi + @muted def element_to_matrix(self, gamma): + verb_level = get_verbose() + set_verbose(0) ans = self._G(gamma)(self._generators) ans.set_immutable() + set_verbose(verb_level) return ans def group(self): @@ -158,32 +172,32 @@ def generators(self): return self._generators def inverse_generators(self): - return tuple(~v for v in self._generators) + 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( - self.generators(), self.inverse_generators(), length + self._generators, self._inverse_generators, length ) def all_elements_up_to_length(self, N): return chain( *[ - enumerate_group_elements(self.generators(), self.inverse_generators(), i) + enumerate_group_elements(self._generators, self._inverse_generators, i) for i in range(N + 1) ] ) def theta_naive(self, m, a=ZZ(0), b=ZZ(1), z=None, gamma=None, **kwargs): + K = self.K if z is Infinity: - return 1 + return K(1) if gamma is not None: if b != ZZ(1): raise ValueError("If gamma is specified, then b should not") b = act(gamma, a) - K = z.parent() max_length = kwargs.pop("max_length", -1) num = K(1) den = K(1) @@ -191,17 +205,22 @@ def theta_naive(self, m, a=ZZ(0), b=ZZ(1), z=None, gamma=None, **kwargs): for i in NN: new_num = K(1) new_den = K(1) - if i > 1 and max_length == -1: + if i > 2 and max_length == -1: can_stop = True for _, g in self.enumerate_group_elements(i): ga = act(g, a) - tmp1 = K(1) if ga is Infinity else z - K(ga) + tmp1 = K(1) if ga is Infinity else z - ga gb = act(g, b) - tmp2 = K(1) if gb is Infinity else z - K(gb) + tmp2 = K(1) if gb is Infinity else z - gb new_num *= tmp1 new_den *= tmp2 - val = (tmp1 / tmp2 - K(1)).valuation() - if val < m - tmp1.valuation(): + try: + val = (tmp1 / tmp2 - 1).valuation() + tmp1val = tmp1.valuation() + except TypeError: + val = -Infinity + tmp1val = 0 + if val < m - tmp1val: can_stop = False num *= new_num den *= new_den @@ -277,7 +296,7 @@ def action_table(self): self._action_table = [] VT = self.NJtree.vertex_table for g in self._generators: - for h in [g, ~g]: + for h in [g, g.adjugate()]: new_list = [] for v in self.NJtree.vertices(): hv = h * v @@ -291,21 +310,30 @@ def action_table(self): self._action_table.append(new_list) return self._action_table + @lazy_attribute + def base(self): + for g in self._generators: + try: + b0, b1 = find_eigenvector_matrix(g).column(0) + return b0 / b1 + except (ZeroDivisionError, NotImplementedError): + pass + raise RuntimeError("Base not found") + def construct_tree(self, level): - g0 = self._generators[0] - b0, b1 = find_eigenvector_matrix(g0).column(0) - base = b0 / b1 T = NeighborJoiningTree( - self.K, [act(g, base) for wd, g in self.all_elements_up_to_length(level)] + self.K, + [act(g, self.base) for wd, g in self.all_elements_up_to_length(level)], ) if level == 1: pp = [choose_leaf(tt) for tt in T.tree[1][:2]] - self.root_vertex = T.BT([T.root, *pp]) + self.root_vertex = T.BT.vertex_from_three_points([T.root, *pp]) self.NJtree = T self._action_table = None return self.NJtree @cached_method + @muted def fundamental_domain(self, level=1, check=True): while True: level += 1 @@ -316,10 +344,11 @@ def fundamental_domain(self, level=1, check=True): except ValueError as e: verbose(str(e)) + @muted def _fundamental_domain(self, i0=None, check=True): tree = self.NJtree if i0 is None: - i0 = tree.vertex_index(self.root_vertex) + i0 = tree.vertices().index(self.root_vertex) adj = tree.adjacency_list() edge_classes = self.edge_classes() open_edges = [(None, i0)] @@ -351,16 +380,35 @@ def _fundamental_domain(self, i0=None, check=True): raise ValueError("Not a correct fundamental domain") good_generators = [] balls = {} - verts = self.NJtree.vertices() - for i, (e0, e1, word) in enumerate(pairing): + T = self.NJtree + verts = T.vertices() + mid = lambda x, y: T.BT.find_midpoint(x, y) + for i, ((v00, v01), (v10, v11), word) in enumerate(pairing): w = self._G.one() for l in word: w = w * (self._G.gen(l // 2) ** ((-1) ** l)) good_generators.append(w) - B0 = find_midpoint(self.K, verts[e0[0]], verts[e0[1]]) - B1 = find_midpoint(self.K, verts[e1[0]], verts[e1[1]]) + g = self.element_to_matrix(w) + B0 = mid(verts[v00], verts[v01]) + B1 = mid(verts[v10], verts[v11]) + if ( + verts[v00].radius == verts[v01].radius + or verts[v10].radius == verts[v11].radius + ): + assert ( + g * B0.complement() != B1.closure() + or g.adjugate() * B1.complement() != B0.closure() + ) + B0 = mid(verts[v00], B0) + B1 = (g * B0.closure()).complement() balls[i + 1] = B1 balls[-(i + 1)] = B0 + if check: + try: + newgens = [self.element_to_matrix(o) for o in good_generators] + test_fundamental_domain(newgens, balls) + except RuntimeError: + raise ValueError("Not a correct fundamental domain (test)") return good_generators, balls, pairing def certify(self, pairing, i0): @@ -426,10 +474,7 @@ def edge_classes(self, action_table=None): action_table = self.action_table() adj = self.NJtree.adjacency_list() n_vertices = len(adj) - edges = set() - for i in range(n_vertices): - for j in adj[i]: - edges.add((i, j)) + edges = set(((i, j) for i in range(n_vertices) for j in adj[i])) r = dict() visited = set() for i in edges: @@ -463,48 +508,31 @@ def __init__( transformation=None, parameters=None, keep_generators=True, - check=False, + check=True, ): if balls is not None: - self.balls = balls + self._balls = balls if transformation is not None: self._transformation = transformation if parameters is not None: self.parameters = parameters self._keep_generators = keep_generators - # if balls is None: - # G = PreSchottkyGroup(K, generators) - # gp = G.group() - # good_gens, good_balls, _ = G.fundamental_domain() - # if all(len(o.Tietze()) == 1 for o in good_gens): - # self.balls = {} - # for j0, gg in enumerate(good_gens): - # i = gg.Tietze()[0] - # 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 parameters is not None: - # self.parameters = parameters - if check: - test_fundamental_domain(generators, self.balls) super().__init__(K, generators) + if check: + self.test_fundamental_domain() + + def test_fundamental_domain(self): + balls = self.balls() + gens = self._generators + test_fundamental_domain(gens, balls) - @lazy_attribute def balls(self): + try: + return self._balls + except AttributeError: + pass K = self.base_ring() - generators = self.generators() + generators = self._generators G = PreSchottkyGroup(K, generators) gp = G.group() good_gens, good_balls, _ = G.fundamental_domain() @@ -522,13 +550,17 @@ def balls(self): else: # print(good_gens) self._generators = [G.element_to_matrix(g) for g in good_gens] + self._inverse_generators = tuple( + [o.adjugate() for o in self._generators] + ) balls = good_balls self._transformation = good_gens + self._balls = balls return balls @cached_method def gens_extended(self): - gens = [(i + 1, o) for i, o in enumerate(self.generators())] + gens = [(i + 1, o) for i, o in enumerate(self._generators)] gens.extend([(-i, o.determinant() ** -1 * o.adjugate()) for i, o in gens]) for i, o in gens: o.set_immutable() @@ -536,49 +568,58 @@ def gens_extended(self): @lazy_attribute def parameters(self): - balls = self.balls + balls = self.balls() ans = [find_parameter(o, balls[i], self.pi) for i, o in self.gens_extended()] for t in ans: t.set_immutable() return ans - @cached_method - def a_point(self, in_interior=True): + def a_point(self): K = self.base_ring() - x = K.random_element() n_iters = 0 while n_iters < 100: n_iters += 1 + x = K.random_element() try: x0 = self.to_fundamental_domain(x)[0] - if in_interior: - if any(x0 in B.closure() for B in self.balls.values()): - raise ValueError - return x0 except ValueError: - x = K.random_element() + continue + if self.in_fundamental_domain(x0, strict=True): + return x0 raise RuntimeError("Reached maximum iterations (100)") def find_point(self, gamma, eps=1, idx=None): - balls = self.balls + balls = self.balls() + gens = self.gens_extended() if idx is not None: B = balls[-idx] if idx > 0: - assert self.generators()[idx - 1] == gamma + if self._generators[idx - 1] != gamma: + B1 = next(balls[-i] for i, g in gens if g == gamma) + # print(f'!! {B = } but {B1 = }') + # print(balls) + B = B1 else: - assert self.generators()[-idx - 1] ** -1 == gamma + if self._generators[-idx - 1] ** -1 != gamma: + B1 = next(balls[-i] for i, g in gens if g == gamma) + # print(f'!! {B = } but {B1 = }') + # print(balls) + B = B1 else: - gens = self.gens_extended() B = next(balls[-i] for i, g in gens if g == gamma) ans = B.center.lift_to_precision() + eps * self.pi ** ZZ(B.radius) - test = lambda x : self.in_fundamental_domain(x,strict=False) + test = lambda x: self.in_fundamental_domain(x, strict=False) try: if not test(ans): - raise RuntimeError(f'{ans = } should be in fundamental domain, but it is not.') - ga = act(gamma,ans) + raise RuntimeError( + f"{ans = } should be in fundamental domain, but it is not." + ) + ga = act(gamma, ans) if not test(ga): - raise RuntimeError(f'gamma * ans = {ga} should be in fundamental domain, but it is not.') + raise RuntimeError( + f"gamma * ans = {ga} should be in fundamental domain, but it is not." + ) except RuntimeError: new_idx = -idx if idx is not None else None ginv = gamma**-1 @@ -592,8 +633,8 @@ def to_fundamental_domain(self, x): and a word w = [i1,...,ik] (in Tietze notation) such that gi1 * ... * gik * z = x """ - gens = self.generators() - inv_gens = self.inverse_generators() + gens = self._generators + inv_gens = self._inverse_generators gens = [None] + list(gens) + list(reversed(inv_gens)) word = [] i = self.in_which_ball(x) @@ -606,45 +647,52 @@ def to_fundamental_domain(self, x): ) x = x1 i = self.in_which_ball(x) - return x, word + return x, tuple(word) 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 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 = }") + found = False + while not found: + z0 = self.a_point() + z1 = act(gamma, z0) + z2, word = self.to_fundamental_domain(z1) + if z0 == z2: + found = True + break + else: # DEBUG: why is this needed? + gens = [(i + 1, o) for i, o in enumerate(self._generators)] + gens.extend( + [(-i, o.determinant() ** -1 * o.adjugate()) for i, o in gens] + ) + for i, g in gens: + gz2 = act(g, z2) + if z0 == gz2: + word.append(-i) + found = True + break + # if not found: + # raise RuntimeError(f"Something went wrong! {z0 = }, {z2 = } {word = }") return tuple(word) + @muted def matrix_to_element(self, g): return self._G(self.word_problem(g)) def in_which_ball(self, x, closure=False): if closure: - return next((i for i, B in self.balls.items() if x in B.closure()), None) + return next((i for i, B in self.balls().items() if x in B.closure()), None) else: - return next((i for i, B in self.balls.items() if x in B), None) + return next((i for i, B in self.balls().items() if x in B), None) def find_equivalent_divisor(self, D): # Compute an equivalent divisor to D with proper support # We replace (a) with (a0) - (g z) + (z) for a0 and z in # the closure of the fundamental domain. p = self.pi - gens = self.generators() - balls = self.balls + gens = self._generators + balls = self.balls() wd = [0 for _ in range(self.genus())] Div = D.parent() ans = Div([]) @@ -653,10 +701,10 @@ def find_equivalent_divisor(self, D): for i in new_wd: wd[abs(i) - 1] += n * sgn(i) ans += n * Div([(1, P0)]) - for g, e in zip(gens, wd): + for i, (g, e) in enumerate(zip(gens, wd)): if e == 0: continue - zz = self.find_point(g) + zz = self.find_point(g, idx=i + 1) ans -= e * Div([(1, zz), (-1, act(g, zz))]) return ans @@ -684,82 +732,83 @@ def theta(self, prec, a, b=None, **kwargs): True """ + K = self.base_ring() 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, base_ring=K, **kwargs) + recursive = kwargs.get("recursive", True) + if recursive: + D0 = self.find_equivalent_divisor(D0) + ans = ThetaOC(self, a=D0, b=None, prec=prec, base_ring=K, **kwargs) z = kwargs.pop("z", None) improve = kwargs.pop("improve", True) if improve or z is not None: ans = ans.improve(prec) return ans if z is None else ans(z) - @cached_method def u_function(self, gamma, prec, a=None, **kwargs): r""" Compute u_gamma """ 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" + gens = self._generators + z = kwargs.get("z", None) + if len(wd) > 1 or wd[0] < 0: + if z is None: + kwargs.pop("z", None) + return lambda z: prod( + self._u_function(gens[abs(i) - 1], prec, a=a)(z) ** ZZ(sgn(i)) + for i in wd ) return prod( - self.u_function(gens[abs(i) - 1], prec, a=a, **kwargs) ** sgn(i) + self._u_function(gens[abs(i) - 1], prec, a=a)(z) ** ZZ(sgn(i)) for i in wd ) + if z is None: + return lambda z: self._u_function(gens[wd[0] - 1], prec, a=a)(z) + else: + return self._u_function(gens[wd[0] - 1], prec, a=a)(z) + @cached_method + def _u_function(self, gamma, prec, a): + r""" + Compute u_gamma + """ + wd = self.word_problem(gamma) + gens = self._generators + assert len(wd) == 1 and wd[0] > 0 if a is None: - a = self.find_point(gamma) # Should pass idx=wd[0] DEBUG ?? + a = self.find_point(gamma, idx=wd[0]) 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) + ans = ThetaOC(self, a=D, b=None, prec=prec, base_ring=K) + ans = ans.improve(prec) + return ans @cached_method def period_matrix(self, prec, **kwargs): - g = len(self.generators()) + 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) + g1 = self._generators[i] + T = self.u_function(g1, prec, a=z1, **kwargs) for j in range(i, g): - g2 = self.generators()[j] + 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) + num = T(z2) + den = T(g2_z2) M[i, j] = num / den M[j, i] = num / den return M @@ -792,40 +841,46 @@ def period(self, i, j, prec, **kwargs): sage: (q11g/q11-1).valuation() > prec True """ - g = len(self.generators()) + g = len(self._generators) if i in ZZ: assert j in ZZ - g1 = self.generators()[i] - g2 = self.generators()[j] + g1 = self._generators[i] + g2 = self._generators[j] else: - a = word_to_abelian(i, g) - b = word_to_abelian(j, g) + ilist = i + jlist = j m = self.period_matrix(prec, **kwargs) - return a * m * b + ans = m.parent().base_ring()(1) + for i in ilist: + i0, isgn = abs(i) - 1, sgn(i) + for j in jlist: + j0, jsgn = abs(j) - 1, sgn(j) + ans *= m[i0, j0] ** ZZ(isgn * jsgn) + return ans z1 = kwargs.pop("z1", None) if z1 is None: z1 = self.a_point() z2 = kwargs.pop("z2", None) if z2 is None: - z2 = self.find_point(g2, eps=1 + self.pi, idx=j) + z2 = self.find_point(g2, eps=1 + self.pi, idx=j + 1) g2_z2 = act(g2, z2) - T = self.u_function(g1, prec, a=z1, **kwargs).improve(prec) - num = T(z2, recursive=False) - den = T(g2_z2, recursive=False) + T = self.u_function(g1, prec, a=z1, **kwargs) + num = T(z2) + den = T(g2_z2) verbose(f"{num = }") verbose(f"{den = }") return num / den def period_naive(self, i, j, prec, **kwargs): - g1 = self.generators()[i] - g2 = self.generators()[j] - z1 = self.find_point(g1, idx=i) + g1 = self._generators[i] + g2 = self._generators[j] + z1 = self.find_point(g1, idx=i + 1) if i == j: z1 = act(g1.adjugate(), z1) - z2 = self.find_point(g2, eps=self.pi + 1, idx=j) + z2 = self.find_point(g2, eps=self.pi + 1, idx=j + 1) else: - z2 = self.find_point(g2, idx=j) + z2 = self.find_point(g2, idx=j + 1) 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 = }") @@ -835,12 +890,20 @@ def period_naive(self, i, j, prec, **kwargs): class Ball(Element): def __init__(self, parent, center, radius, is_open=True, is_complement=False): - self.radius = radius - self.center = parent.K(center).add_bigoh(self.radius + 1) + self.radius = QQ(radius) + self.center = parent.K(center).add_bigoh(ZZ(self.radius.ceil()) + 1) self.is_complement = is_complement self.is_open = is_open Element.__init__(self, parent) + def __hash__(self): + return hash(str(self.get_key())) + + @cached_method + def get_key(self): + r = self.radius + return r if self.center.is_zero(r) else (self.center.add_bigoh(r), r) + def _repr_(self): comp = "P^1 - " if self.is_complement else "" opn = "" if self.is_open else "+" @@ -883,16 +946,7 @@ def interior(self): return C(self.parent(), self.center, self.radius, is_open, self.is_complement) def __eq__(self, other): - if not ( - self.radius == other.radius - and self.is_open == other.is_open - and self.is_complement == other.is_complement - ): - return False - if self.is_complement: - return self.complement() == other.complement() - else: - return self.center in other and other.center in self + return self.parent().are_equal(self, other) def __ne__(self, other): return not self == other @@ -908,13 +962,22 @@ def _acted_upon_(self, g, on_left): a, b, c, d = g.list() x = self.center x = x.lift_to_precision() - try: - if (x + d / c).valuation() > self.radius: - inf_in_image = True - else: + if self.is_open: + try: + if (x + d / c).valuation() > self.radius: + inf_in_image = True + else: + inf_in_image = False + except ZeroDivisionError: + inf_in_image = False + else: + try: + if (x + d / c).valuation() >= self.radius: + inf_in_image = True + else: + inf_in_image = False + except ZeroDivisionError: inf_in_image = False - except ZeroDivisionError: - inf_in_image = False if inf_in_image: gpx = (g.determinant() / c**2).valuation() @@ -935,121 +998,103 @@ def _acted_upon_(self, g, on_left): class Balls(UniqueRepresentation, Parent): Element = Ball - def __init__(self, K): + def __init__(self, K, oriented=True): self.K = K + self._oriented = oriented Parent.__init__(self) + def base_ring(self): + return self.K + def _repr_(self): return "Set of balls for P^1(K), with K = " + str(self.K) def _element_constructor_(self, center, radius, is_open=True, is_complement=False): + if not self._oriented: + is_open = False + is_complement = False return self.element_class(self, center, radius, is_open, is_complement) - -class BTVertex(ModuleElement): - def __init__(self, parent, a, r_valuation): - KK = parent.base_ring() - try: - ln = len(a) - ps = [KK(i) for i in a if i is not Infinity] - if len(ps) == 2: - self.r_valuation = (ps[0] - ps[1]).valuation() - self.a = ps[0].add_bigoh(self.r_valuation) + def are_equal(self, left, right): + if self._oriented: + if not ( + left.radius == right.radius + and left.is_open == right.is_open + and left.is_complement == right.is_complement + ): + return False + if left.is_complement: + return left.complement() == right.complement() else: - d12 = (ps[0] - ps[1]).valuation() - d23 = (ps[1] - ps[2]).valuation() - d13 = (ps[0] - ps[2]).valuation() - maxval = max([d12, d23, d13]) - self.r_valuation = maxval - self.a = ( - ps[0].add_bigoh(maxval) if d12 == maxval or d13 == maxval else ps[1] - ) - except TypeError: - try: - self.r_valuation = ZZ(r_valuation) - self.a = KK(a).add_bigoh(r_valuation) - except TypeError: - self.r_valuation = Infinity - self.a = KK(a) - ModuleElement.__init__(self, parent) - - def contained_in(self, other): - return (self.a - other.a).valuation() <= other.r_valuation - - def __hash__(self): - return hash(self.get_key()) - - @cached_method - def is_zero(self): - return self.a.is_zero(self.r_valuation) - - @cached_method - def get_key(self): - return ( - self.r_valuation - if self.is_zero() - else (self.a.add_bigoh(self.r_valuation), self.r_valuation) - ) - - def _repr_(self): - return f"Ball({self.a},|π|^{self.r_valuation})" - - def __eq__(self, other): - return self.get_key() == other.get_key() - - def _add_(self, other): - C = self.__class__ - return C( - self.parent(), self.a + other.a, min(self.r_valuation, other.r_valuation) - ) - - def _lmul_(self, b): - C = self.__class__ - KK = self.parent().base_ring() - assert b != 0 - return C(self.parent(), self.a * b, self.r_valuation + KK(b).valuation()) - - def inv(self): - C = self.__class__ - if self.is_zero(): - return C(self.parent(), 0, -self.r_valuation) + return left.center in right and right.center in left else: - return C( - self.parent(), 1 / self.a, self.r_valuation - 2 * self.a.valuation() + return ( + left.radius == right.radius + and (left.center - right.center).valuation() >= left.radius ) - def _acted_upon_(self, g, on_left): - if hasattr(g, "matrix"): - g = g.matrix() - if g[1, 0] == 0: - return 1 / g[1, 1] * (g[0, 0] * self + g[0, 1]) + def find_midpoint(self, v0, v1): + K = self.base_ring() + t = K(v0.center - v1.center).valuation() + if t < v0.radius and t < v1.radius: + distance = (v0.radius - t) + (v1.radius - t) else: - t = (g[1, 0] * self + g[1, 1]).inv() - r = (g[0, 1] - g[0, 0] * g[1, 1] / g[1, 0]) * t + g[0, 0] / g[1, 0] - return r - - -class BTTree(UniqueRepresentation, Module): - Element = BTVertex - - def __init__(self, K): - self.prec = K.precision_cap() - Module.__init__(self, base=K) + distance = abs(v0.radius - v1.radius) + v, comp = (v1, False) if v1.radius > v0.radius else (v0, True) + BK = Balls(K) + return BK( + v.center, v.radius - QQ(distance) / 2, is_open=not comp, is_complement=comp + ) - def _element_constructor_(self, a, rval=None): - if rval is None: - rval = Infinity - return self.element_class(self, a, rval) + def vertex_from_three_points(self, a): + KK = self.base_ring() + ln = len(a) + ps = [KK(i) for i in a if i is not Infinity] + if len(ps) == 2: + val = (ps[0] - ps[1]).valuation() + return self.element_class(self, ps[0].add_bigoh(val), val) + else: + assert len(ps) == 3 + d12 = (ps[0] - ps[1]).valuation() + d23 = (ps[1] - ps[2]).valuation() + d13 = (ps[0] - ps[2]).valuation() + maxval = max([d12, d23, d13]) + center = ( + ps[0].add_bigoh(maxval) if d12 == maxval or d13 == maxval else ps[1] + ) + return self.element_class(self, center, maxval) + + +def cross_ratio(points): + [p1, p2, p3, p4] = points + if p1 is Infinity: + return (p3 - p2) / (p4 - p2) + elif p2 is Infinity: + return (p3 - p1) / (p4 - p1) + elif p3 is Infinity: + return (p4 - p2) / (p4 - p1) + elif p4 is Infinity: + return (p3 - p1) / (p3 - p2) + else: + return (p3 - p1) * (p4 - p2) / ((p4 - p1) * (p3 - p2)) - def _an_element_(self): - return self.element_class(self, 1, 2) - def _coerce_map_from_(self, S): - if self.base().has_coerce_map_from(S): - return True +def four_point_configuration(K, pts): + p1, p2, p3, p4 = pts + v1 = abs(K(cross_ratio([p1, p2, p3, p4])).valuation()) + v2 = abs(K(cross_ratio([p1, p3, p2, p4])).valuation()) + v3 = abs(K(cross_ratio([p1, p4, p2, p3])).valuation()) + if v2 > 0 and v3 > 0: + return p3, v2 + elif v1 > 0 and v3 > 0: + return p2, v1 + elif v1 > 0 and v2 > 0: + return p1, v1 + else: + return None, 0 -def four_point_configuration(K, pts): +def four_point_configuration_works(K, pts): V = [[K(pi - pj).valuation() for pj in pts] for pi in pts] v2 = abs(V[0][1] + V[2][3] - V[1][2] - V[0][3]) v3 = abs(V[0][1] + V[2][3] - V[1][3] - V[0][2]) @@ -1076,7 +1121,7 @@ def choose_leaf(tree): class NeighborJoiningTree(SageObject): def __init__(self, K, leaves): self.K = K - self.BT = BTTree(self.K) + self.BT = Balls(self.K, oriented=False) leaves = uniq(leaves) self.root = leaves[0] self.tree = [Infinity, [[None, leaves[1]], [None, leaves[2]]]] @@ -1086,6 +1131,28 @@ def __init__(self, K, leaves): self.add_leaf(leaf) self.update_vertex_table() + def to_graph(self): + tree = self.tree + G = Graph(weighted=True) + G.add_vertex(0) + V = [(0, tt) for tt in tree[1]] + W = [] + newlabel = 0 + Gdict = {} + while len(V) > 0: + for label, T in V: + newlabel += 1 + G.add_vertex(newlabel) + if leaf(T): + G.add_edge(label, newlabel, label=str(T[1].add_bigoh(3))) + Gdict[newlabel] = T[1] + else: + G.add_edge(label, newlabel, label=str(T[0])) + W.extend([(newlabel, tt) for tt in T[1]]) + V = W + W = [] + return G, Gdict + def to_string(self, subtree): s = "--(" + str(subtree[0]) + ")--" if not leaf(subtree): @@ -1143,17 +1210,15 @@ def add_leaf(self, new_leaf, subtree=None): def vertices(self, subtree=None): if subtree is None: subtree = self.tree - if not leaf(subtree): - p1 = choose_leaf(subtree[1][0]) - p2 = choose_leaf(subtree[1][1]) - return sum( - (self.vertices(tt) for tt in subtree[1]), [self.BT([self.root, p1, p2])] - ) - else: + if leaf(subtree): return [] - - def vertex_index(self, v0, subtree=None): - return next((i for i, v in enumerate(self.vertices(subtree)) if v == v0)) + children = subtree[1] + p0 = choose_leaf(children[0]) + p1 = choose_leaf(children[1]) + V = [self.BT.vertex_from_three_points([self.root, p0, p1])] + for tt in children: + V.extend(self.vertices(tt)) + return V def update_vertex_table(self): self.vertex_table = {hash_vertex(v): i for i, v in enumerate(self.vertices())} @@ -1163,14 +1228,8 @@ def adjacency_list(self): return self._adjacency_list r = tuple([] for _ in self.vertices()) for v1, v2 in self.adjacency(): - try: - i1 = self.vertex_table[hash_vertex(v1)] - except KeyError: - i1 = self.vertices().index(v1) - try: - i2 = self.vertex_table[hash_vertex(v2)] - except KeyError: - i2 = self.vertices().index(v2) + i1 = self.vertex_table[hash_vertex(v1)] + i2 = self.vertex_table[hash_vertex(v2)] r[i1].append(i2) r[i2].append(i1) self._adjacency_list = r @@ -1183,7 +1242,7 @@ def adjacency(self, subtree=None, root=None): if not leaf(subtree): p1 = choose_leaf(subtree[1][0]) p2 = choose_leaf(subtree[1][1]) - v = self.BT([self.root, p1, p2]) + v = self.BT.vertex_from_three_points([self.root, p1, p2]) result = [[root, v]] if root is not None and not leaf(subtree) else [] return sum((self.adjacency(tt, v) for tt in subtree[1]), result) else: @@ -1192,7 +1251,7 @@ def adjacency(self, subtree=None, root=None): def test_fundamental_domain(gens, balls): all_gens = [(i + 1, g) for i, g in enumerate(gens)] + [ - (-(i + 1), g**-1) for i, g in enumerate(gens) + (-(i + 1), g.adjugate()) for i, g in enumerate(gens) ] fails = [] for i, g in all_gens: @@ -1200,14 +1259,21 @@ def test_fundamental_domain(gens, balls): if i == j: continue if (balls[i].closure()).intersects(balls[j].closure()): - fails.append((i, j)) - verbose(f"Test *failed* for balls {i = } and {j = }") - else: - verbose(f"Test passed for balls {i = } and {j = }") - if g * balls[-i].complement() != balls[i].closure(): + fails.append( + ( + i, + j, + ) + ) + # verbose(f"Test *failed* for balls {i = } and {j = }") + # else: + # verbose(f"Test passed for balls {i = } and {j = }") + B0 = g * balls[-i].complement() + B1 = balls[i].closure() + if B0 != B1: fails.append(i) - verbose(f"Test *failed* for {i = }") - else: - verbose(f"Test passed for {i = }") + # verbose(f"Test *failed* for {i = }. {B0 = }, {B1 = }") + # else: + # verbose(f"Test passed for {i = }") if len(fails) > 0: raise RuntimeError(f"Some tests failed. ({fails})") diff --git a/darmonpoints/theta.py b/darmonpoints/theta.py index 7b34518..65e3273 100644 --- a/darmonpoints/theta.py +++ b/darmonpoints/theta.py @@ -25,6 +25,38 @@ infinity = Infinity +def eval_rat(D, z): + if z == Infinity: + return 1 + ans = 1 + fails = 0 + for P, n in D: + if P == Infinity: + continue + zP = z - P + if zP == 0: + fails += n + else: + ans *= zP**n + # assert fails == 0 # DEBUG ! + return ans + + +def eval_rat_derivative(D, z): + ans = 0 + fails = 0 + for P, n in D: + if P == Infinity: + continue + zP = z - P + if zP == 0: + fails += n + else: + ans += n * zP**-1 + # assert fails == 0 # DEBUG ! + return ans * eval_rat(D, z) + + def act(mtrx, z): """ Compute the action of a 2x2 matrix on an element. @@ -36,16 +68,10 @@ def act(mtrx, z): if z is Infinity: return Infinity if c == 0 else a / c try: - return (a * z + b) / (c * z + d) #a / c + (b - a * d / c) / (c * z + d) + return (a * z + b) / (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: - return (a * z + b) / (c * z + d) - except PrecisionError: - return Infinity + def element_to_matrix(wd, generators): ans = wd(generators) @@ -73,7 +99,7 @@ def __init__(self, G, a=0, b=None, prec=None, **kwargs): self.G = G self.p = G.pi generators = G.generators() - balls = G.balls + balls = G.balls() if prec is None: try: self.prec = K.precision_cap() @@ -89,45 +115,64 @@ def __init__(self, G, a=0, b=None, prec=None, **kwargs): "Must specify a degree-0 divisor, or parameters a and b" ) else: - D = self.Div([(1, a), (-1, b)]) + D = self.Div([(1, K(a)), (-1, K(b))]) self.a = a self.b = b self.D = D self.m = 1 - self.z = K['z'].gen() + self.z = K["z"].gen() params = G.parameters gens_ext = G.gens_extended() - # self.val will contain the 0 and 1 terms D1 = sum((g * D for (i, g), tau in zip(gens_ext, params)), self.Div([])) - self.val = self.z.parent()(1) - self.val *= prod((self.z - P) ** n for P, n in (D + D1)) + self.val = D + D1 + # self.val = self.z.parent()(1) + # self.val *= prod((self.z - P) ** n for P, n in (D + D1) if P != Infinity) + self.MM = MeromorphicFunctions(self.K, self.p, self.prec) self.Fnlist = [{}] - D1dict = {i : g * D for (i, g), tau in zip(gens_ext,params)} - for (i, g), tau in zip(gens_ext,params): - gD = sum( - g * val for j, val in D1dict.items() if j != -i - ) - self.Fnlist[0][i] = MeromorphicFunctions(self.K, self.p, self.prec)(gD, tau) + D1dict = {i: g * D for (i, g), tau in zip(gens_ext, params)} + for (i, g), tau in zip(gens_ext, params): + gD = sum(g * val for j, val in D1dict.items() if j != -i) + self.Fnlist[0][i] = self.MM(gD, tau) def improve(self, m): gens_ext = self.G.gens_extended() params = self.G.parameters + action_data = {} + for (i, gi), tau in zip(gens_ext, params): + for j, Fj in self.Fnlist[-1].items(): + if i != -j: + action_data[i, j] = ( + self.MM.get_action_data(gi, Fj._parameter, tau), + tau, + ) for it in range(m): if self.m >= self.prec: + self.Fnlist = [ + { + ky: sum((F[ky] for F in self.Fnlist[1:]), self.Fnlist[0][ky]) + for ky in self.Fnlist[0] + } + ] return self tmp = {} for (i, gi), tau in zip(gens_ext, params): for j, Fj in self.Fnlist[-1].items(): if i != -j: - vl = Fj.left_act_by_matrix(gi, tau) + vl = Fj.fast_act(action_data[i, j]) try: tmp[i] += vl except KeyError: tmp[i] = vl self.Fnlist.append(tmp) self.m += 1 + self.Fnlist = [ + { + ky: sum((F[ky] for F in self.Fnlist[1:]), self.Fnlist[0][ky]) + for ky in self.Fnlist[0] + } + ] return self def improve_one(self): @@ -177,7 +222,7 @@ 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) + recursive = kwargs.get("recursive", True) ans = 1 if recursive: z0, wd = G.to_fundamental_domain(z) @@ -189,17 +234,16 @@ def evaluate(self, z, **kwargs): ) else: z0 = z - ans *= self.val(z0) + ans *= eval_rat(self.val, z0) ans *= prod(F(z0) for FF in self.Fnlist for ky, F in FF.items()) return ans - def eval_derivative(self, z, recursive=False): + def eval_derivative(self, z, recursive=False, return_value=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) + vz = eval_rat(self.val, z) Fnz = {} for FF in self.Fnlist: for ky, F in FF.items(): @@ -210,9 +254,15 @@ def eval_derivative(self, z, recursive=False): Fnz[ky] = Fz # val' * Fn(z) Fnzall = prod((o for o in Fnz.values())) - valder = self.val.derivative()(z) + valder = eval_rat_derivative(self.val, z) - v0 = self.val(z) * Fnzall + v0 = vz * Fnzall # need to add the terms of val * Fn' - tmp = sum(f.eval_derivative(z) / f(z) for FF in self.Fnlist for f in FF.values()) - return valder * Fnzall + tmp * v0 + tmp = sum( + f.eval_derivative(z) / f(z) for FF in self.Fnlist for f in FF.values() + ) + ans = valder * Fnzall + tmp * v0 + if return_value: + return ans, v0 + else: + return ans diff --git a/darmonpoints/util.py b/darmonpoints/util.py index c1623cd..4fdc67b 100644 --- a/darmonpoints/util.py +++ b/darmonpoints/util.py @@ -1,4 +1,5 @@ import configparser +import functools import sys import types from functools import reduce @@ -39,6 +40,18 @@ from sage.structure.factorization import Factorization +def muted(func): + @functools.wraps(func) + def ff(*args, **kwargs): + verb_level = get_verbose() + set_verbose(0) + ans = func(*args, **kwargs) + set_verbose(verb_level) + return ans + + return ff + + def is_smooth(x, B): for p in B: x /= p ** (valuation(x, p)) @@ -1958,6 +1971,7 @@ def recognize_J( CEloc, _ = get_C_and_C2(Eloc, qEpows, Cp, precp) EH = E.change_ring(HCF) + success = False for twopow in twopowlist: success = False power = QQ(twopow * known_multiple) ** -1 @@ -2281,18 +2295,22 @@ def print_table_latex(self, header_string=None): s += "}" + frame_str + "\n" # first row s += " & ".join( - LatexExpr(x) - if isinstance(x, (str, LatexExpr)) - else "$" + latex(x).strip() + "$" + ( + LatexExpr(x) + if isinstance(x, (str, LatexExpr)) + else "$" + latex(x).strip() + "$" + ) for x in rows[0] ) s += " \\\\" + frame_str + head_row_str + "\n" # other rows for row in rows[1:]: s += " & ".join( - LatexExpr(x) - if isinstance(x, (str, LatexExpr)) - else "$" + latex(x).strip() + "$" + ( + LatexExpr(x) + if isinstance(x, (str, LatexExpr)) + else "$" + latex(x).strip() + "$" + ) for x in row ) s += " \\\\" + frame_str + "\n" diff --git a/pyproject.toml b/pyproject.toml index 8b21c4a..b1a66eb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,5 @@ [build-system] -requires = ["setuptools>=61.2","Cython", "wheel"] +requires = ["setuptools","Cython", "wheel"] build-backend = "setuptools.build_meta" [project]