Elementary Matrix Decomposition and The Computation of Darmon Points with Higher Conductor

We extend the algorithm of Darmon-Green and Darmon-Pollack for computing p-adic Darmon points on elliptic curves to the case of composite conductor. We also extend the algorithm of Darmon-Logan for computing ATR Darmon points to treat curves of nontrivial conductor. Both cases involve an algorithmic decomposition into elementary matrices in congruence subgroups {\Gamma}(N) for ideals N in certain rings of S-integers. We use these extensions to provide additional evidence in support of the conjectures on the rationality of Darmon points.


Introduction
Let E be an elliptic curve over Q of conductor pM , with p a prime not dividing M . Let K be a real quadratic field in which p is inert and all the primes dividing M are split, and denote by H p = P 1 (K p ) \ P 1 (Q p ) the K p -points of the p-adic upper half plane.
A construction of Darmon [Dar01] associates to every τ ∈ H p a local point P τ ∈ E(K p ), which is defined as a certain period of the modular form f corresponding to E under the Modularity Theorem. The points P τ are conjectured to be rational over ring class fields of K, and to behave in many aspects as Heegner points.
An algorithm for the effective calculation of p-adic Darmon points was given in [DG02] and was improved in [DP06]. Some P τ 's were computed in concrete examples and checked to be padically close to global points, providing extensive numerical evidence in support of the conjectures. However, due to the restrictions imposed by the algorithm, only elliptic curves of prime conductor p (that is, with M = 1) could be treated.
In the articles [DG02] and [DP06] it is crucial to assume that M = 1 when applying the "continued fraction trick" (see [DG02,p. 42]) in order to transform certain semi-indefinite integrals into double integrals. The present article provides a different procedure for performing this step when M > 1 provided that E satisfies a mild condition (see Assumption 3.2 for its precise statement).
In other words, we extend the algorithm of [DG02] and [DP06] to a much larger class of curves. As an application, we compute p-adic Darmon points on curves of composite conductor and we check that they are close to global points, which provides new experimental evidence in support of the validity of Darmon's construction in its stated generality.
The method for transforming semi-indefinite integrals into double integrals is based on an algorithmic decomposition into elementary matrices in congruence subgroups Γ 1 (N) for ideals N in certain rings of S-integers. This can be seen as an effective version of the congruence subgroup problem. In particular we improve on a theorem of Cooke-Weinberger (see Corollary 2.4).
In addition, essentially the same method of elementary matrix decompositions can also be applied to ATR points, a different instance of Darmon  elliptic curves over totally real fields. Although the construction of ATR points is radically different (for instance, they are defined by means of complex periods), their explicit computation has some formal similarities with respect to the p-adic setting. In particular the methods used until the present have all used a "continued fraction trick" which in this case only applies to curves of trivial conductor (cf. [Gär11,p. 108] for a discussion of this issue). The present article provides also a method for computing ATR Darmon points in curves of non-trivial conductor.
The rest of the article is organized as follows. In Section 2 we introduce an algorithm for computing elementary matrix decompositions in certain congruence subgroups. We state it in a level of generality so that it can be applied both to the p-adic and the ATR setting. In Section 3 we recall the definition of p-adic Darmon points, we discuss the algorithm for computing them in curves of composite level and we include some tables of numerical computations performed using it. Finally, in Section 4 we briefly recall ATR points and we make explicit the method for computing them in curves of non-trivial conductor, as well as a detailed example of a numerical verification of Darmon's conjecture in this case.
Acknowledgments. It is a pleasure to thank John Voight for helpful comments on elementary matrix decompositions, Robert Pollack for clarifying to us some details on his implementation of the overconvergent modular symbols algorithm, and Henri Darmon for many valuable observations. Guitart wants to thank the Max Planck Institute for Mathematics for their hospitality and financial support during his stay at the Institute, where part of the present work has been carried out.

Preliminaries: elementary matrix decomposition in Γ 1
Let F be a number field with ring of integers O, and let S be a finite set of places of F containing the archimedean ones. Let O S denote the subring of F consisting of those elements whose valuation is non-negative at all the places outside S.
For an ideal N of O S , we denote by Γ 1 (N) the subgroup of SL 2 (O S ) defined as and by E 1,N the subgroup of Γ 1 (N) generated by the elementary matrices of the form (2.2) 1 0 y 1 with y ∈ N and If the group of units O × S is infinite, then Γ 1 (N) = E 1,N ; see for instance [Vas72] for a proof of this result, and also for its relation with the Congruence Subgroup Problem. Therefore, when O × S is infinite every matrix in Γ 1 (N) can be factored into a product of elementary matrices of type (2.2). However, the proof given in [Vas72] is not explicit, so it does not give rise to an algorithm for systematically performing such decomposition. In this section we see that the results and techniques introduced in [CW75] can be adapted to provide -assuming GRH-such an algorithm, in the particular case where F has at least a real archimedean place. Thus from now on we will assume that O × S is infinite and that F has at least a real place. The following lemma is stated and proved in [BMS67, Lemma 2.2 (b)] using the notation of Mennicke symbols. For the convenience of the reader we restate it in terms of explicit matrix formulas.
Lemma 2.1. Let γ = a b c d be an element in Γ 1 (N). Suppose that c = u + ta for some unit u ∈ O × S and some t ∈ O S , and let T be the matrix Then T γ = ( 1 x 0 1 ) for some x ∈ O S . In particular, Proof. A direct computation shows that the first column of T γ is ( 1 0 ) . Since T γ belongs to Γ 1 (N), we see that Observe that, since a − 1 and c belong to N , identity (2.3) already expresses γ as a product of elementary matrices of type (2.2). The next step is to show that, assuming GRH, one can reduce to the case where a is congruent to a unit mod c by multiplying γ by an elementary matrix. Before stating the result of [CW75] that grants this, we recall some terminology related to ray class groups. Let I m be the multiplicative group of ideals of O relatively prime to m, and let There is a natural map i : F 1,m → I m given by x → (x). The quotient is the ray class group modulo m. The following result is [CW75, Theorem 2.13].  (1) Iterate over the elements λ ∈ O to find λ such that a ′ = a + λc generates a prime ideal and (4) Use Lemma 2.1 to find an expression of γ ′ as a product of elementary matrices Proof. We can, and do, choose e ∈ O × S such that ea and ec belong to O. By applying Theorem 2.2 with m = (ec), we see that there exists a prime ideal q in the same class of is integral and principal, we see that q is also principal. Therefore, q = (q) for some q ≡ ea (mod (ec)), so that q = ea + λec for some λ ∈ O. But q and This justifies that step (1) of the algorithm can be accomplished. Now the class of c ′ in (O S /a ′ ) × can be represented by some unit u ∈ O × S , and this justifies step (3).
In the particular case N = O S , the result in [CW75, Th. 2.14] asserts that any matrix γ ∈ SL 2 (O S ) can be expressed as a product of at most 7 elementary matrices. Theorem 2.3 above and expression (2.3) give the following generalization to arbitrary ideals N, which also slightly improves the number of elementary matrices needed to 5.
Corollary 2.4. Under the assumption of GRH, every matrix in Γ 1 (N) is a product of at most 5 elementary matrices of type (2.2).

Computation of p-adic Darmon points on curves with composite conductor
In this section we explain the algorithm for computing p-adic Darmon points in curves of composite conductor and some of the related computational issues. After briefly reviewing the definition of p-adic Darmon points in §3.1, in §3.2 we describe the algorithm to transform the semi-indefinite integrals appearing in the definition of these points into double multiplicative integrals. In §3.3 we explain an efficient way of computing these double integrals, and finally in §3.4 we comment on the calculations that have been carried out in support of Darmon's conjecture.
3.1. Review of p-adic Darmon points. Our presentation of the necessary background and the definition of p-adic Stark-Heegner points follows closely [DP06,§1], to which we refer the reader for more details and, in fact, for an excellent account of this material in the prime level case.
Let E be an elliptic curve over Q of conductor N = pM , with p a prime not dividing M . Let f (z) = n≥1 a n e 2πinz be the weight two newform on Γ 0 (N ) whose L-series coincides with that of E. The coefficient a p is 1 (resp. −1) if E has split (resp. non-split) multiplicative reduction at p.
Let R be the order in M 2 (Z[1/p]) consisting of matrices that are upper triangular modulo M , and let Γ = R × 1 denote its group of units of determinant 1. Let K be a real quadratic field in which p is inert and all primes dividing M are split. Set H p = P 1 (K p ) \ P 1 (Q p ), in which Γ acts by Möbius transformations. The p-adic Darmon point construction yields a map given in terms of certain p-adic integrals, and whose definition ultimately relies on the Z-modular symbol attached to E.
For w ∞ ∈ {±1} we denote by I f : P 1 (Q) × P 1 (Q) → Z the Z-valued modular symbol attached to E and w ∞ . That is to say where ω f = 2πif (z)dz and Ω + , Ω − ∈ R >0 are the unique periods with the property that the map I f defined in this way takes values in Z and in no proper ideal of Z. To simplify the exposition we assume for the rest of the article that ω ∞ = 1, but the construction works very similarly for ω ∞ = −1.
A Z-valued measure on P 1 (Q p ) is a finitely additive function µ from the set of compact open subsets of P 1 (Q p ) to Z. In [Dar01] Darmon attaches to each pair r, s ∈ P 1 (Q) a Z-valued measure µ f {r → s} on P 1 (Q p ) with total measure 0 by defining where the limit is taken over increasingly finer coverings U of P 1 (Q p ) by compact open subsets U α , and x α is any point in U α . If µ is Z-valued the multiplicative integral can be defined by replacing the Riemann sum by a Riemann product: For τ 1 , τ 2 ∈ H p and r, s ∈ P 1 (Q p ) Darmon defines a K p -valued double integral as where log denotes a fixed branch of the p-adic logarithm. Since These two integrals are related by the formula Therefore, it gives rise to a group homomorphism where the subscript Γ denotes the subgroup of Γ-coinvariants. Let Z[Γ] denote the group ring of Γ and let I Γ denote the augmentation ideal, defined by the exact sequence Tensoring with I Γ over Z and taking Γ-coinvariants gives Since I Γ is generated over Z by elements of the form γ − 1, choosing base points τ ∈ H p and x ∈ P 1 (Q) one can define integration maps x e Γ ω f := Int × τ,x (e Γ · 1 ⊗ (γ − 1)) ∈ K × p /Λ. This semi-indefinite integral satisfies the following properties: e Γ ω f , for all τ ∈ H p , r, s, t ∈ Γx; (2) × 3.1.4. Darmon points. Let q denote the p-adic period of E and let Φ Tate : K × p /q Z → E(K p ) be Tate's uniformization map.
Given τ ∈ H p , let O τ be the ring of matrices in R that have the vector (τ, 1) as eigenvector, which is isomorphic to a Z[1/p] order of K. Let H τ denote its ring class field which we can, and do, view as a subfield of K p by choosing a prime of H τ above p. The stabilizer Γ τ of τ in Γ is a cyclic group of infinite order isomorphic to O × τ,1 / ±1 . Let γ τ be a generator of Γ τ , and define Conjecture 3.1 (Darmon). The local point P τ = Φ Tate (J τ ) belongs to E(H τ ).

Computation of semi-indefinite integrals.
In order to effectively compute the points J τ one needs to compute the semi-indefinite integrals The method used in [DG02] and [DP06] boils down to using properties (1), (2) and (3) of semiindefinite integrals to express them in terms of double integrals, which can be effectively computed either via Riemann products as in [DG02] or, more efficiently, via overconvergent modular symbols as in [DP06].
The algorithm for expressing semi-indefinite integrals in terms of double integrals of [DG02] and [DP06] is based on the continued fraction algorithm, and it only works under the assumption that M = 1 (i.e., that E has conductor equal to p). In this section we introduce an algorithm that, assuming GRH and under some mild assumptions on E, works for all levels M . In Section 3.3 we will see that the resulting definite double integrals obtained by this method can also be computed using overconvergent modular symbols, by suitably adapting the techniques of [DP06].
From now on we assume that the modular form f satisfies the following condition, which is non-vacuous only when M is prime.
Assumption 3.2. There is a d > 1, d | M such that f has eigenvalue 1 with respect to the Atkin-Lehner operator W d .

Remark 3.3. Under Assumption 3.2, the double multiplicative integral is also invariant under the
LetΓ be the subgroup of PGL 2 (Q) generated by Γ and w d . Then, by replacing Γ byΓ in the argument of Section 3.3 one can extend the definition of the semi-indefinite integrals × τ s r e Γ ω f to all pairs r, s ∈ P 1 (Q) lying in the sameΓ-orbit.
Let Γ 1 be the congruence subgroup of SL 2 (Z[1/p]) defined as Γ 1 = {γ ∈ Γ : γ ≡ ( 1 * 0 1 ) (mod M )} ⊂ Γ. Using the properties of the multiplicative integral it is easy to see that Therefore, replacing P τ by a multiple of it if necessary we can always assume that γ τ belongs to Γ 1 (but see also Remark 3.4).
Observe that Γ 1 is one of the groups treated in Section 2. Indeed, with the notation as in that section, if we let F = Q, S = {∞, p} and N = M · Z[1/p] we have that Γ 1 = Γ 1 (N). In particular, the algorithm described in Theorem 2.3 gives an effective method (under our running assumption of GRH) for computing a decomposition of γ τ of the form where the matrices U i and L i are of the form In particular, L i and U i belong to Γ. Then, for G ∈ Γ we have that (3.5) In view of decomposition (3.3), repeated application of (3.4) and (3.5) transforms the semi-indefinite integral (3.2) into a product of double multiplicative integrals. Observe that Assumption 3.2 is crucial in (3.5) because of the integral × τ 0 ∞ e Γ ω f . Indeed, the cusps 0 and ∞ are never in the same Γ-orbit if M > 1, but w d ∞ = 0 so they are in the sameΓ-orbit and thanks to Remark 3.3 the integral × τ 0 ∞ e Γ ω f is well defined. Remark 3.4. As we already mentioned, one can overcome the fact that γ τ = a b c d does not generally belong to Γ 1 by computing an appropriate power (J τ ) m . In some cases one can apply an alternative procedure instead, which turns out to be more convenient in the actual computations and it allows for the computation of J τ itself. Namely, if a ≡ p n (mod M ) for some integer n, then the matrix g = p −n 0 0 p n belongs to Γ and with gγ τ belonging to Γ 1 .

3.3.
Computation of the definite double integrals. As we have seen in Subsection 3.2 the computation of the period J τ is reduced to products of integrals of the form Write µ for the measure µ f {r → s}, and consider a decomposition of P 1 (Q p ) into a disjoint union of L open balls of the form which will be fixed later. This yields in turn a decomposition Fix such an i and let g = g i be the corresponding matrix, written g = a b c d . We are thus reduced to calculating Apply the change of variables x = g · t to get Let log p be the unique homomorphism K × p → K p such that log p (1 − t) = − ∞ n=1 t n /n and log p (p) = 0. It is surjective, with kernel ker log p : where U is the group of roots of unity in K × p . Suppose that we can express the integrand as a power series in t of the form with α n belonging to Z p for all n ≥ 1. Then the expression in (3.8) converges for t ∈ Z p and is constant modulo p vp(α0)+1 . Therefore the expression of (3.6) can be determined modulo this power of p by performing L evaluations. The logarithm log p (J τ ) of the period is evaluated by noting that Interchanging the infinite sum with the integral, finding log p (J τ ) boils down to calculating which is the nth moment of µ at gZ p . This data can be efficiently computed in time polynomial in the number of p-adic digits of precision, thanks to the methods of Darmon and Pollack (see for instance [DP06,display 23]). Finally, one recovers the original period via the formula where ζ is the Teichmuller lift of the unit part modulo p of the multiplicative integral. Note also that v p (J τ ) is the sum of the valuations of the α 0 appearing in the decomposition (3.7). In order to find the power series in (3.8), we calculate: where g is the matrix Note that for any two matrices g, h we have hg = gh −1 . Also, note that for any choice of h, a decomposition gives rise to another decomposition Therefore by choosing an appropriate h ∈ GL 2 (Q p ) ∩ M 2 (Z) we can assume that v p (τ 1 − a) = 0 for all a = 0, 1, . . . , p − 1 and that v p (τ 2 ) ≥ 0. Note that in order to obtain a power series as in (3.8), the matrices g = a b c d that we consider should satisfy The conditions on τ 1 and τ 2 imply that the matrix corresponding to the contribution of P 1 (Q p ) \ Z p satisfies (3.9), and we concentrate on the integral on Z p . Let r ≥ 1 be the largest integer such that τ 2 is congruent to some integer modulo p r . Let t 2 ∈ Z be a representative for the class of τ 2 (mod p r ). Write also t (i) for the representative of t 2 (mod p i ) in the range 0, . . . , p i − 1. We can then use the decomposition given by the matrices g in the set This decomposition consists of exactly p + 1 + r(p − 1) opens.

Numerical computations.
To test our methods we have written a Sage implementation of the above algorithms, modifying an existing implementation written by Robert Pollack which in turn adapted part of the code originally written in Magma by Darmon and Pollack ([DP06]). The code can be found in the second author's webpage. Given an elliptic curve of conductor N = pM and a quadratic field K, the code finds all the optimal embeddings of level N of K into M 2 (Z[ 1 p ]), and computes the Stark-Heegner period corresponding to the fixed point of K acting on H p via each embedding to a prescribed precision. The Tate parametrization yields the coordinates of the Stark-Heegner point on E(K p ) which are then recognized as algebraic coordinates using standard routines.
Apart from gathering numerical evidence in support of Darmon's conjecture, it is also worth remarking that the relative large height of the points thus found would make it impossible to find them using naive point search methods. This is, therefore, the only known method to finding points of infinite order on such curves.
The rest of this subsection contains the evidence that we have collected in support of Conjecture 3.1. Although we do not intend to be exhaustive, we provide examples of curves of small composite conductor which satisfy Assumption 3.2 and for which we have been able to positively test the conjecture. For each of these curves we consider all the real quadratic fields K of discriminant D < 200 allowed by the splitting conditions on p and M ; for each such field, we consider τ ∈ H p such that H τ equals the Hilbert class field of K, and we are able to recognize in all the cases P τ as an algebraic point defined over the Hilbert class field of K. For those fields with nontrivial class group, we give the relative minimal polynomial h D of the X-coordinate of the point.

Computation of ATR Darmon points on curves of non-trivial conductor
In Section 3 we have seen that the algorithm of Theorem 2.3 can be used in the computation of the semi-indefinite integrals entering the definition of p-adic Darmon points. It is a substitute for the continued fractions trick of [DG02] and [DP06].
There is another type of Darmon points, called ATR, whose definition also relies in certain semiindefinite integrals. Although the framework is different (e.g., they are points on elliptic curves over number fields, and the integrals are complex instead of p-adic), the formal properties satisfied by the semi-indefinite integrals are the same in both settings. In the ATR case, the continued fraction algorithm over number fields had been used for computing ATR points on curves with trivial conductor (cf. [DL03], [GM12]). Using a method analogous to that of Section 3.2, Theorem 2.3 also allows for the computation of ATR points on curves with non-trivial conductor.
To be more precise, let F be a real quadratic number field of narrow class number 1 and let O denote its ring of integers. Let E be an elliptic curve over F of conductor N, and let Γ be the congruence subgroup consisting of matrices in SL 2 (O) that are upper triangular modulo N.
Assuming that E is modular, there is a Hilbert modular form f of parallel weight two and level N whose L-series coincides with that of E. Let ω f denote the corresponding Γ-invariant differential 2-form on H × H (with Γ acting on it via the product of the two embeddings of F into R).
The following is analogous to Assumption 3.2.
where the matrices U i and L i are of the form In particular, they belong to Γ and the same expressions of (3.4) and (3.5) (changing the multiplicative by the additive notation) express τ γτ ∞ ∞ ω + f in terms of usual double integrals of the form τ2 τ1 ∞ 0 ω + f , which in principle can be evaluated by integrating the Fourier expansion of ω + f . It is worth remarking that, although the method described above certainly expresses the semiindefinite integrals in terms of definite ones, the running time to directly compute the resulting double integrals to a useful accuracy often turns out to be too high. The problem is that if the limits of integration are too close to the real axis, then the number of Fourier coefficients needed to sum the series to an accurate precision is too high.
These kind of computational difficulties seem to be inherent to the ATR setting, as they were also to some extent present in the initial work of Darmon and Logan [DL03]. In [GM12] some methods for accelerating the computation of double integrals in the trivial level case were introduced. The authors believe that similar techniques should be applied to the non-trivial level setting in order to perform a systematical calculation similar to the one of Section 3.4.
In spite of this, in some simple examples it is possible to directly compute the integrals provided by Theorem 2.3, and hence to compute approximations to the ATR points. The following is an example of this, which we detail because it provides numerical evidence of the validity of Darmon's conjecture in elliptic curves of non-prime conductor.
The conductor of E is N = ( √ 5 + 6), which has norm 31. This curve was previously considered in [Gre11] and [Gär11] (but note the typo in the displayed equation in both references). Let α = 1− √ 5 and let K = F ( √ α). The embedding ϕ : K ֒→ M 2 (F ) sending ω to the matrix W = 3 − ω −1 8 − 3ω −3 + ω is an optimal embedding of level N. Under the embedding F ֒→ R sending √ 5 to the positive square root of 5, the fixed point of W acting on C × is τ = 0.439291418991 + i · 0.353408129753, and the image of the unit (−3 + 2ω)ω + 4 − 3ω ∈ O × K under ϕ is The ATR point attached to (the maximal order of) K is The determinant of γ τ is ω + 1, which is a unit. However its upper left entry is not congruent to 1 modulo N. If we let u = − √ 5 − 2, which is a fundamental unit of F , then −4 + 3ω ≡ u (mod N).
We use this decomposition to transform (4.2) into a sum of usual double integrals. The resulting integrals have limits not too close to the real axis (the smallest imaginary part is ≃ 0.011). Integrating the Fourier series with coefficients a m with norm of m up to 180, 000 gives J τ to an accuracy of approximately 12 digits: J τ ≃ −4.828954817077 + i · 4.534696532333.
There is a point P of infinite order in E(K) having x-coordinate equal to 18883/2420α−16127/2420 (this was found using naive search algorithms); let J denote its corresponding image in C/Λ E . Then the equality J τ = −2J (mod Λ E ) holds up to the computed accuracy, giving numerical evidence of the equality P τ = −2P and, therefore, of the rationality of P τ .