Made functions less dependent on base field

Author: mmasdeu

darmonpoints/arithgroup.py

@@ -85,7 +85,7 @@ def geodesic_circle(alpha, beta, return_equation=True):
     assert C.imag().abs() < 10**-10, C
     try:
         r2 = (alpha - C).norm()
-    except AttributError:
+    except AttributeError:
         r2 = (alpha - C) ** 2
     x, y = PolynomialRing(C.parent(), 2, names="x,y").gens()
     return (
@@ -284,12 +284,12 @@ def get_embgquats_twisted(self, P=None):
         if P is None:
             CC = ComplexField(800)
             P = CC(90) / CC(100) * CC.gen()
-        embgquats = self.embgquats
         Pconj = P.conjugate()
         m = Matrix([[Pconj, -P], [1, -1]])
         madj = m.adjugate()
         return [None] + [tuple((madj * o * m).list()) for o in self.embgquats[1:]]
 
+
     def get_word_rep(self, delta, P=None):  # fuchsian generic
         while True:
             try:
@@ -390,7 +390,7 @@ def magma_word_problem(self, Gm, x):
     def draw_mat_list(self, x0, g1, mlist, color="blue"):
         phi = self.get_archimedean_embedding(1000)
         P = hyperbolic_arc(x0, act_flt(phi(g1.quaternion_rep), x0), color="red")
-        vlist = G._fundamental_domain
+        vlist = self.fundamental_domain()
         for g in mlist:
             new_polygon = [act_flt(phi(g), v) for v in vlist]
             P += hyperbolic_polygon(new_polygon, color=color)
@@ -401,19 +401,18 @@ def mat_list(self, x1, x2, check_fundom=True):  # generic
         Returns a list S of matrices such that the geodesic (x1, x2) is contained in the union
         of the translates s*D (with s in S) of the standard fundamental domain D.
         """
+        CC = x1.parent()
+        emb = self.get_archimedean_embedding(CC.precision())
         verbose("Calling mat_list with x1 = %s and x2 = %s" % (x1, x2))
-        x1_orig = x1
-        x2_orig = x2
         n = 0
         g = 1
         if check_fundom and not self.is_in_fundom(x1):
             t0, g = self.find_fundom_rep(x1)
-            x1, x2 = t0, self(g**-1) * x2
+            x1, x2 = t0, act_flt(emb(self(g**-1).quaternion_rep), x2)
 
         # Here we can assume that x1 is in the fundamental domain
         ans = [self(g)]
         while not self.is_in_fundom(x2):
-            found = False
             intersection_candidates = []
             for i, (v1, v2) in enumerate(self.fundamental_domain_data()):
                 z = intersect_geodesic_arcs(v1, v2, x1, x2)
@@ -423,29 +422,22 @@ def mat_list(self, x1, x2, check_fundom=True):  # generic
                     z1, z2 = perturb_point_on_arc(x1, x2, z, eps)
                     if not self.is_in_fundom(z1):
                         z1, z2 = z2, z1
-                    try:
-                        assert self.is_in_fundom(
-                            z1
-                        ), "z1 and z2 are both outside of fundom!"
-                        assert not self.is_in_fundom(
-                            z2
-                        ), "z1 and z2 are both inside of fundom!"
+                    if self.is_in_fundom(z1) and not self.is_in_fundom(z2):
                         intersection_candidates.append(((z1, z2), i, (v1, v2)))
-                    except AssertionError:
-                        pass
             assert len(intersection_candidates) > 0, "No intersection candidates!!"
             if len(intersection_candidates) > 1:
                 verbose(intersection_candidates)
             z1, z2 = intersection_candidates[
                 ZZ.random_element(len(intersection_candidates))
             ][0]
+            found = False
             for g in self.gquats[1:]:
-                t0 = self(g) ** -1 * z2
+                t0 = act_flt(emb((self(g) ** -1).quaternion_rep), z2)
                 if self.is_in_fundom(t0):
                     assert not found, "Found more than one!"
                     found = True
                     x1 = t0
-                    x2 = self(g) ** -1 * x2
+                    x2 = act_flt(emb((self(g) ** -1).quaternion_rep), x2)
                     verbose("x1 = %s, x2 = %s" % (x1, x2))
                     ans.append(ans[-1] * self(g))
                     break  # DEBUG
@@ -472,7 +464,7 @@ def _fix_sign(self, x, N):
         return x
 
     def find_fundom_rep(
-        self, z0_H, v0=None, return_alternative=False, max_iters=100
+        self, z0_H, v0 =None, return_alternative=False, max_iters=100
     ):  # generic -- Take z0 to the fundamental domain
         r"""
         Returns t0 and g such that g * t0 = z0_H, and t0 is in the fundamental domain
@@ -541,8 +533,8 @@ def find_fundom_rep(
             oldji, oldgg = ji, gg
             wd.append(gg)
             verbose("New g = %s\t delta = %s\t z0=%s" % (gg, delta, z0))
-            t0 = self(delta) * z0_H
-            t1 = self(wd[-1] ** -1 * delta) * z0_H
+            t0 = act_flt(emb(self(delta).quaternion_rep),z0_H)
+            t1 = act_flt(emb(self(wd[-1] ** -1 * delta).quaternion_rep),z0_H)
         delta_inv = delta**-1
         if return_alternative:
             return (t0, delta_inv), (t1, delta_inv * wd[-1])
@@ -1288,17 +1280,17 @@ def mat_list(
         of the translates s*D (with s in S) of the standard fundamental domain D.
         """
         CC = x1.parent()
-        x1 += CC.random_element(10**-1) * x1.imag() ** 2
-        x2 += CC.random_element(10**-1) * x2.imag() ** 2
+        x1 += CC.random_element(0.1) * x1.imag() ** 2
+        x2 += CC.random_element(0.1) * x2.imag() ** 2
         if debug:
             from sage.repl.rich_output.pretty_print import show
 
             self._debug_plot = hyperbolic_polygon(
                 [
-                    CC(-0.5, 10 ^ 2),
+                    CC(-0.5, 10 ** 2),
                     self.fundamental_domain()[1],
                     self.fundamental_domain()[2],
-                    CC(0.5, 10 ^ 2),
+                    CC(0.5, 10 ** 2),
                 ]
             )
         if v0 is None:

darmonpoints/arithgroup_element.py

@@ -11,29 +11,17 @@
 
 from sage.algebras.quatalg.all import QuaternionAlgebra
 from sage.functions.generalized import sgn
-from sage.functions.trig import arctan
-from sage.groups.group import AlgebraicGroup
-from sage.matrix.all import Matrix, matrix
-from sage.misc.all import cached_method, lazy_attribute, walltime
+from sage.misc.all import cached_method, lazy_attribute
 from sage.misc.misc_c import prod
-from sage.misc.persist import db
-from sage.modules.all import vector
-from sage.modules.free_module import FreeModule_generic
 from sage.rings.all import (
     QQ,
-    RR,
     ZZ,
-    ComplexField,
     NumberField,
     PolynomialRing,
     Qp,
     QuadraticField,
-    RealField,
 )
 from sage.structure.element import MultiplicativeGroupElement
-from sage.structure.parent import Parent
-from sage.structure.richcmp import richcmp
-from sage.structure.sage_object import SageObject, load, save
 
 from .util import *
 
@@ -211,7 +199,8 @@ def __iter__(self):
     def matrix(self, prec=-1):
         return self.embed(prec)
 
-    def acton(self, x):
+    def _act_on_(self, x, on_left):
+        assert on_left == True
         return act_flt(self.matrix(), x)
 
     def conjugate_by(self, w):

darmonpoints/darmonvonk.py

@@ -804,7 +804,6 @@ def __init__(self, G, x0, gamma_tau, tau_p):
 
         self._tau = tau_p
         tau = self._tau
-        tau_bar = tau.trace() - tau
         self._gamma_tau = gamma_tau
 
         action_map = lambda g, x: act_flt(vinf(g.quaternion_rep), x)
@@ -835,14 +834,13 @@ def evaluate(self, gamma):
         Div = Divisors(self._tau.parent())
         if gamma.quaternion_rep == 1 or gamma.quaternion_rep == -1:
             return Div(0)
-        wlist = set([])
         x0 = self._x0
         gamma_of_x0 = gamma * x0
         tau = x0.parent()(self._tau_inf)
         tau_bar = x0.parent()(self._tau_inf_bar)
         ans = Div(0)
+        c1 = HH.get_geodesic(x0, gamma_of_x0)        
         for g in self.all_matrix_candidates(gamma):
-            c1 = HH.get_geodesic(x0, gamma_of_x0)
             gtau = g * tau
             gtau_bar = g * tau_bar
             c2 = HH.get_geodesic(gtau, gtau_bar)

darmonpoints/meromorphic.py

@@ -98,7 +98,7 @@ def __call__(self, D):
     def polynomial_approximation(self, z, m):
         K = self.parent().base_ring()
         a, b, c, d = self._parameter.list()
-        phi = lambda Q: K(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
@@ -121,7 +121,7 @@ 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: K(a) / c if Q == Infinity else K(a * Q + b) / K(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
@@ -152,7 +152,7 @@ def eval_derivative(self, D):
             raise NotImplementedError
         K = self.parent().base_ring()
         a, b, c, d = self._parameter.list()
-        phi = lambda Q: K(a) / c if Q == Infinity else K(a * Q + b) / K(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)
@@ -249,9 +249,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 - ((c * Q + d) / (a * Q + b)) * t) ** ZZ(n)
         else:
-            den *= (1 - K(((c * Q + d) / (a * Q + b))) * t) ** ZZ(-n)
+            den *= (1 - ((c * Q + d) / (a * Q + b)) * t) ** ZZ(-n)
     ans = Ps(num).add_bigoh(prec) * ~(Ps(den).add_bigoh(prec))
     return ans