Many related to theta calculations

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 *
+

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)
         )

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):

darmonpoints/divisors.py

@@ -260,7 +260,7 @@ def rational_function(self, as_map=False, z=None):
         else:
             if z is None:
                 K = self.parent()._field
-                z = K['z'].gen()
+                z = K["z"].gen()
             return prod(((1 - z / P) ** n for P, n in self), z.parent()(1))
 
     def as_list_of_differences(self):

darmonpoints/homology.py

@@ -55,9 +55,7 @@ def lattice_homology_cycle(
     xi2 = H1({})
     for x, a in xlist:
         xi1 += H1({G(x.quaternion_rep): Div(tau)}).mult_by(a)
-        xi2 += H1({G(wp**-1 * x.quaternion_rep * wp): wpinv_mat * Div(tau)}).mult_by(
-            a
-        )
+        xi2 += H1({G(wp**-1 * x.quaternion_rep * wp): wpinv_mat * Div(tau)}).mult_by(a)
     xi10 = xi1
     xi20 = xi2
     while True:

darmonpoints/meromorphic.py

@@ -39,6 +39,7 @@
 
 from .divisors import Divisors
 
+
 def evalpoly(poly, x):
     # Initialize result
     try:
@@ -97,22 +98,26 @@ def evaluate(self, D):  # meromorphic functions
         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)
+        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)
         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)
+            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(ev(P) ** n for P, n in D)
         else:
@@ -126,16 +131,18 @@ def eval_derivative(self, D):
             raise NotImplementedError
         K = self.parent().base_ring()
         a, b, c, d = self._parameter.list()
-        phi = lambda Q : a / c if Q == Infinity else (a * Q + b) / (c * Q + d)
+        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 (
@@ -163,23 +170,23 @@ def _add_(self, right):  # multiplicative!
         if self._parameter != right._parameter:
             raise RuntimeError("Trying to add incompatible series")
         prec = self.parent()._prec
-        ans = (self.power_series() * right.power_series()).list()[: prec]
+        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")
         prec = self.parent()._prec
-        ans = (self.power_series() / right.power_series()).list()[: prec]
+        ans = (self.power_series() / right.power_series()).list()[:prec]
         return self.__class__(self.parent(), ans, self._parameter, check=False)
 
     def _neg_(self):  # multiplicative!
         prec = self.parent()._prec
-        ans = (~self.power_series()).list()[: prec]
+        ans = (~self.power_series()).list()[:prec]
         return self.__class__(self.parent(), ans, self._parameter, check=False)
 
     def scale_by(self, k):  # multiplicative!
-        ans = (self.power_series() ** k).list()[: prec]
+        ans = (self.power_series() ** k).list()[:prec]
         return self.__class__(self.parent(), ans, self._parameter, check=False)
 
     def _repr_(self):
@@ -194,7 +201,9 @@ def _acted_upon_(self, g, on_left):
 
     def fast_act(self, key):
         zz_ps_vec, param = key
-        return self.__class__(self.parent(), self._value * zz_ps_vec, param, check=False)
+        return self.__class__(
+            self.parent(), self._value * zz_ps_vec, param, check=False
+        )
 
     def left_act_by_matrix(self, g, param=None):  # meromorphic functions
         t = cputime()
@@ -208,6 +217,7 @@ def left_act_by_matrix(self, g, param=None):  # meromorphic functions
         ans = self._value * zz_ps_vec
         return self.__class__(self.parent(), ans, param, check=False)
 
+
 def divisor_to_pseries(parameter, Ps, data, prec):
     t = Ps.gen()
     a, b, c, d = parameter.list()
@@ -228,6 +238,7 @@ class MeromorphicFunctions(Parent, UniqueRepresentation):
     TESTS:
 
     """
+
     Element = MeromorphicFunctionsElement
 
     @staticmethod
@@ -268,18 +279,18 @@ def get_action_data(self, g, oldparam, param, prec=None):
             .truncate(prec)
         )
         ans = [zz.parent()(1), zz]
-        while len(ans) < prec: # DEBUG - was 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
+                .truncate(prec)  # DEBUG - was prec + 1
             )
             ans.append(zz_ps)
         set_verbose(verb_level)
         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
+                m[i, j] = aij  # i, j entry contains j-th term of zz^i
         return m
 
     def base_ring(self):

darmonpoints/ocmodule.py

@@ -690,6 +690,7 @@ class AddMeromorphicFunctions(Parent, CachedRepresentation):
     sage: A == B
     True
     """
+
     Element = AddMeromorphicFunctionsElement
 
     def __init__(self, K, twisting_matrix=None):

darmonpoints/padicperiods.py

@@ -4,14 +4,15 @@
 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 save, load, db
-from sage.arith.srange import srange
+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
 from .integrals import integrate_H1
@@ -1553,11 +1554,21 @@ def jacobian_matrix(fvec):
 def compute_lvec_and_Mlist(prec):
     R0 = PolynomialRing(QQ, 3, names="p")
     R = R0.fraction_field()
-    PS = PowerSeriesRing(QQ,names='p',num_gens=3)
+    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()])
+    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)
@@ -1573,7 +1584,7 @@ def load_lvec_and_Mlist(prec):
 
 
 def evaluate_twisted_jacobian_matrix(p1, p2, p3, Mlist):
-    mlist = [o.polynomial()(p1,p2,p3) for o in Mlist]
+    mlist = [o.polynomial()(p1, p2, p3) for o in Mlist]
     retvec = mlist[:6]
     h1 = mlist[6]
     h2 = mlist[7]
@@ -1586,7 +1597,7 @@ def evaluate_twisted_jacobian_matrix(p1, p2, p3, Mlist):
     return Matrix(3, 3, retvec)
 
 
-def compute_twisted_jacobian_data(fvec, x,y,z, prec):
+def compute_twisted_jacobian_data(fvec, x, y, z, prec):
     f1, f2, f3 = fvec.list()
     Mlist = [
         o.truncate(prec)
@@ -1600,7 +1611,7 @@ def compute_twisted_jacobian_data(fvec, x,y,z, prec):
             f3.derivative(x),
             f3.derivative(y),
             f3.derivative(z),
-            f3
+            f3,
         ]
     ]
     return Mlist
@@ -1609,7 +1620,11 @@ def compute_twisted_jacobian_data(fvec, x,y,z, prec):
 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 lvec[0](x0, y0, z0), lvec[1](x0, y0, z0), ((1 - z0) / (1 + z0)) ** 2 * lvec[2](x0, y0, z0)
+        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
@@ -1641,6 +1656,7 @@ def HalfPeriodsInTermsOfLambdas(
         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 [
@@ -1659,8 +1675,9 @@ def ev(x0, y0, z0):
         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

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):

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):