How to compute $\sum_{\gamma \in \text{SL}(2,\mathbb Z), \|\gamma\|<10}\frac{1}{\|\gamma \|^4}$ in Sagemath

Question

Sorry for this question, but how to compute

$$\sum_{\gamma \in \text{SL}(2,\mathbb Z), \|\gamma\|<1000}\frac{1}{\|\gamma \|^4}$$

in Sagemath? Here the norm is just Frobebius norm (square root of the sum of squares of entries)

I tried something like

    from sage.all import *
        G = SL(2,ZZ)
        s = 0
        for r in G:
            if r.norm() < 10:
                s += 1/r.norm()^4
        s

but there is an error and I couldn't figure out how to use SL(2,ZZ) properly.

---

To save the computational effort, I changed the range from 1000 to 10.

Answer 1

I don't know Sagemath syntax, but G is an infinite group. So your for loop is asking to loop over infinitely many things. That's an abstract reason why that approach won't work, syntax aside.

You know that all four entries of a matrix that contributes to the sum will be in $\left\{-1000,\cdots,1000\right\}$ . So it's not efficient, but you could sum over $(a,b,c,d)\in\left\{-1000,\cdots,1000\right\}^4$ and only add the $\frac{1}{\sqrt{a^2+b^2+c^2+d^2}}$ contribution when $ ad-bc=1$.

You could improve efficiency by summing like this:

$$\sum_{a=-1000}^{1000}\sum_{b=-\sqrt{1000^2-a^2}}^{\sqrt{1000^2-a^2}}\sum_{c=-\sqrt{1000^2-a^2-b^2}}^{\sqrt{1000^2-a^2-b^2}}\sum_{d=-\sqrt{1000^2-a^2-b^2-c^2}}^{\sqrt{1000^2-a^2-b^2-c^2}}$$

And more by not looking over all those $d$-values. Instead, check if $\frac{1+bc}{a}$ is an integer. If so, that is the only value of $d$ to use. And in that case do include the contribution to the overall sum. (Of course you have to handle the $a=0$ tuples separately.)

Answer 2

Here's a quick and dirty method (there's probably something better, but I don't have the time to look for it right now):

We know that an element $M \in \mathsf{SL}(2,\mathbb{Z})$ of (frobenius) norm $< 1000$ looks like a matrix $\begin{pmatrix} a & b \\ c & d \end{pmatrix}$ with integer entries, $ad-bc = 1$, and $\sqrt{a^2 + b^2 + c^2 + d^2} < 1000$. Now for the crude part: This, in particular, tells us that each $|a|,|b|,|c|,|d| < 1000$, so we can just check all of these.

In code, we get:

N = 1000
total = 0
count = 0
for a in range(-N,N):
  for b in range(-N,N):
    for c in range(-N,N):
      for d in range(-N,N):
        count += 1
        if count % 1000 == 0: print(count, "/", (2*N)^4)
    norm = sqrt(a^2 + b^2 + c^2 + d^2)
    if norm < 1000 and a*d - b*c == 1:
      total += (1/norm^4).n()

though be warned, this will take a very long time! It will print out how far along in the process it is as it goes, but it will have to try $2001^4 \approx 10^{13}$ many matrices. You can be more clever with the endpoints in order to speed this up, but I'll leave that optimization to you.

---

I hope this helps ^_^

Answer 3

Below, the matrix $A$ stays for $A=A(a,b,c,d)$, slight abuse of notation, it is seen as a function in the variables $a,b,c,d\in\Bbb Z$. Its norm will be denoted by $|A|$, so $|A|^2=a^2+b^2+c^2+d^2$. I will set $$ N=1000 $$ for an easy typing. Let $\Gamma=\operatorname{SL}(2,\Bbb Z)$. Let $S=S(N)$ be the sum to be computed: $$ S =S(N)=\sum _{A\in\Gamma\\ |A|<N} \frac 1{|A^4|} \ , $$ First of all, we compute the contribution of all those $A$ that involve at least one zero entry. If $a=0$, then $bc=-1$, there are two cases for $b$, $b=\pm1$, and $c$ is adjusted accordingly as $c=-1/b=\-p1$, and the value of $d$ obeys $-N<d<N$. The corresponding sum is $$ S_{a=0} =2\sum_{-N<d<N} \frac 1{(0+1+1+d^2)^2} =4\sum_{0<k<N} \frac 1{(0+1+1+k^2)^2} +\frac 2{(0+1+1+0)^2} \ , $$ and it can be easily computed numerically. In a similar manner, introduce $S_{d=0}$, then $S_{a=0}=S_{d=0}$. The term with $a=d=0$ is considered in both sums, so the contribution to $S$ from the matrices $A$ with $ad=0$ is $$ S_{ad=0} = 8\sum_{0<k<N} \frac 1{(0+1+1+k^2)^2} +\frac 2{(0+1+1+0)^2} \ . $$ A similar story concerns the case $bc=0$, and the contribution of both cases is $$ S_{abcd=0} = 16\sum_{0<k<N} \frac 1{(0+1+1+k^2)^2} + 1 \ . $$ This "special case" is cleared. Consider now the contribution from all other cases.

Can we have two (or more) entries of same maximal modulus. If two such entries are in the same line or column, then this common modulus value divides the determinant, so it is $\pm 1$. So all entries are $\pm 1$ by maximality. Taking the determinant modulo two we can replace all $\pm 1$ by $1$ modulo two, so the determinant is zero modulo two. Contradiction. So the two values with maximal modulus are on the same diagonal. Let $k$ be this modulus, so the product on this diagonal is $k^2$ in modulus. Then on the other diagonal we have in modulus at most $(k-1)^2$, and the difference can not be $\pm 1$.

So there is exactly one maximal modulus entry, call it key. If this key is in the main diagonal $(a,d)$, it's ok, else we exchange columns, and multiply the first column with $-1$. So take the contribution below times $2$. If the key is $A$ it's ok, else exchange $a$ and $d$. So take the contribution below times $4$. The key is now $a$. If the key $a$ is negative, replace $A$ by $-A$. So the key $a$ becomes positive. And we take the contribution below times $8$. In the diagonal $(b,c)$ insist that $b>0$ by possibly exchanging signs simultaneously. So contribution below times $16$.

Since i like positive numbers (when writing code), i'd like to have also $d>0$. Well, if $d<0$ we change signs in the row $(c,d)$, thus changing the determinant from $+1$ to $-1$. But now we have a clear situation, $c$ is also positive, so $a,b,c,d>0$, $ad-bc =\pm1$, and $a>b,c,d$. I case of $b=c$ we consider the term once, else $b\ne c$, and we insist that $b>c$ after exchanging, and considering this term times two.

Now there is a simple path to implement. We let $a$ run in the interval $(1,N)$. Then $d$ runs in the interval $(0,a)$. We build $ad\pm 1$ and factorize it in positive integers as $bc$ with $c\le b<a$. And of course, we restrict to the case $a^2+b^2+c^2+d^2 < N^2$. In case of $b=c$ we consider the term $1/(a^2 +b^2 +c^2+d^2)$ times one, else times two.

---

Code:

R = RealField(300)
def S0(N):
    return R(1) + 16*sum([ 1/(R(2) + k^2)^2 for k in [1..N-1]])
def S1(N):
    S = R(0)
    for a in [1..N-1]:
        for d in [1..a-1]:
            if a^2 + d^2 > N^2:    break    # the d loop
            for bc in (ad - 1, ad + 1):
                for c in bc.divisors():
                    b = ZZ(bc / c)
                    if b < c:    break    # the c loop
                    norm2 = a^2 + b^2 + c^2 + d^2
                    if b < a and norm2 < N^2:
                        S += R(1) / norm2^2 * (1 if b==c else 2)
    return 16*S

And we get:

sage: S0(1000) + S1(1000)
4.18456554053875394382856565430709947660176138858177166092612407947397830461918770929226114

I hope this is ok, no time for structural checks.

---

Note: Alternative implementation in pari/gp:

{S0(N) = 1. + 16.*sum(k=1, N-1, 1./(2. + k^2)^2)}
{S1(N) = local(S,a,b,c,d,norm2);
 S = 0.;
 for(a=2, N-1, 
   for(d=1, a-1,
     if(a^2 + d^2 < N^2,
       for(pow=0, 1, bc = a*d + (-1)^pow;
         divs = divisors(bc);
         for(k=1, length(divs),
           c = divs[k]; b = bc / c; norm2 = a^2 + b^2 + c^2 + d^2;
           if(c < a && b < a && norm2 < N^2,
             S += 1 / norm2^2; 
           ))))));
 16*S;}

This gives:

? \p 100
   realprecision = 115 significant digits (100 digits displayed)
? #
   timer = 1 (on)
? S0(1000) + S1(1000)
cpu time = 9,337 ms, real time = 9,347 ms.
%24 = 4.184565540538753943828565654307099476601761388581771660926124079473978304619187709292261024171064063