3301.sagews

We will first look at the sequence a_n = ( coefficient of x^(n-1) in (x^3+ax+b)^((n-1)/2) ).  

We could expand the polynomial (x^3+ax+b)^((n-1)/2) and read off the coefficient, but it is faster to compute only the coefficient we want.

Below is the code for the long computation.  We introduce the polynomial ring R, which is the ring of polynomials in x with integer coefficients.  (ZZ means the integers, QQ means the rational numbers, RR means the real numbers.)

Then we compute (x^3+x+1)^5 and call it u; and finally extract the coefficient of x^10.  (I don't know of a simpler command to do this.  I first shift the polynomial by x^(-10), which makes the coefficient of x^10 the constant term and deletes the terms with degree less than 10; then I plug in x=0.)

R.<x>=ZZ['x']
print R

Univariate Polynomial Ring in x over Integer Ring

u=(x^3+x+1)^5
print u
print u.shift(-10)(0)

x^15 + 5*x^13 + 5*x^12 + 10*x^11 + 20*x^10 + 20*x^9 + 30*x^8 + 35*x^7 + 30*x^6 + 31*x^5 + 25*x^4 + 15*x^3 + 10*x^2 + 5*x + 1 20

The faster way to compute a_n is to add up the correct binomial and multinomial coefficients.  Two commands are useful for this; but they work in slightly different ways.  "binomial(m,n)" means "m choose n"=(m!)/(n!(m-n)!).  "multinomial(n1,n2,...,nk)" means "n1+n2+...+nk choose n1,n2,...,nk"=((n1+n2+...+nk)!)/(n1!n2!...nk!).

print binomial(5,2)
print multinomial(3,2)

10 10

To compute just the coefficient of x^(n-1), we add up all products of (n-1)/2 choices of x^3, a*x, and b that yield x^(n-1).  We can choose at most (n-1)/3 of x^3, and then enough of a*x to reach x^(n-1), choosing b for the remaining terms.  But we cannot choose more than (n-1)/2 terms total; this means we must choose at least (n-1)/4 of x^3.  We use the floor and ceiling functions ("floor" and "ceil") to limit the range of i.

def ap(a,b,n):
# This function returns the coefficient of x^(n-1) in (x^3+a*x+b)^((n-1)/2).
    ap=0
    if n%2==0:
        return 0
    else:
        k=(n-1)/2
        for i in srange(ceil((n-1)/4),floor((n-1)/3)+1):
            ap=ap+multinomial(i,n-1-3*i,2*i-k)*a^(n-1-3*i)*b^(2*i-k)
        return ap

ap(1,1,15)

105

To see the Atkin and Swinnerton-Dyer congruences, we look at terms a_(mp^n), where p is a prime (not 2 or 3).  For instance, we can look at a_(5^n).  The numbers a_n get large very quickly.

for i in srange(6):
    print ap(1,1,5^i)

1 2 9339 2848456737948291676796723 197289322971585384883256188658365039688581036394317503893254806591079231499924302356181041749886539784432271684534013810594254386 7864874119195563212913146276683596075400747940102741202222310984308352345821212669506940215381879705948820429945024579667480215019774768795241623580576917952577768376071825851856917143440595221941020710276426690903282148717656779976961397624325097913429985764192574070756697110893695453874376297596548363146533831061449742616781047597060122523310852507687259740104435860571257533105212668307499077053245964859597042864373548707641966422353278392244056805111363720390181345218357631375764254848810100017845820021291471766751947206300708410708810068410343237255164872569830771249122067382323014024944439952152134753201067912915493241460035526066493852

To get more manageable numbers, we can take a_(5^n) modulo 5^(n+1).  The ASD congruences say that 

a_(5^(n+1))-A_5*a_(5^n)+5*a_(5^(n-1)) = 0 modulo 5^(n+1)

For the curve y^2=x^3+x+1, A_5=-3.  (You can see by computing ap(1,1,5) that A_5 = 2 (mod 5).  It could be either 2 or -3, since it must have absolute value less than or equal to 2*sqrt(5).  Either by counting points or by testing both possibilities, you can see that A_5=-3.)

for i in srange(6):
    print ap(1,1,5^i)%5^(i+1)

1 2 89 473 1261 9477

for i in srange(1,6):
    print (ap(1,1,5^(i+1))+3*ap(1,1,5^i)+5*ap(1,1,5^(i-1)))%5^(i+1)

0 0 0 0 0

for i in srange(1,6):
    print (ap(1,1,3*5^(i+1))+3*ap(1,1,3*5^i)+5*ap(1,1,3*5^(i-1)))%5^(i+1)

0 0 0 0 0

Notice how long these calculations take; if we tried to compute ap(1,1,5^7) or ap(1,1,5^8), it could take a very long time.  If a computation is taking too long, you can interrupt it using the "action" toolbar at the top.  You can use the modulo operator "%" to test the congruences; you can also use the function N.digits(p) to list the base p digits.  In the calculation below, the first 3 digits are 0, which indicates that the value is congruent to 0 modulo 5^3, as required by the ASD congruences.

(ap(1,1,125)+3*ap(1,1,25)+5*ap(1,1,5)).digits(5)

[0, 0, 0, 3, 4, 3, 1, 3, 2, 3, 2, 3, 2, 0, 3, 4, 2, 3, 2, 2, 4, 4, 3, 4, 3, 4, 1, 3, 2, 2, 3, 1, 2, 4, 4]

Now we look for supercongruences, choices of a so that the coefficients a_n for the curve y^2=x^3+ax satisfy congruences stronger than those given by the Atkin and Swinnerton-Dyer congruences.

def apa(a,n):
# This function returns the coefficient of x^(n-1) in (x^3+a*x)^((n-1)/2).
    if n%4!=1:
        return 0
    else:
        k=(n-1)/4
        return binomial(2*k,k)*a^k

We will look for choices of a so that 
a_(mp^(n+1))-A_p*a_(mp^n)+p*a_(mp^(n-1))=0 (mod p^(2*n+1))
rather than just =0 (mod p^(n+1)).  

To do this we can check all choices of a modulo p^(2*n+1); but if you do, you will notice that there are a lot of solutions, and you get more and more as n increases.  In fact, if a is a solution, then a+k*p^(n+1) is also a solution.  So we only need to check all choices of a modulo p^(n+1).  For values of a such that we have supercongruences, we print a.digits(p), which will help us see the patterns more easily.

(Note that we don't we rule out a=0 mod p, because then we don't have an elliptic curve over the field of p elements.  Also, notice that if p=3 mod 4, we have A_p=a_p=0 always; so you might expect different behavior for p=1 mod 4 and p=3 mod 4.)

p=5; n=3
for a in srange(p^(n+1)):
    if a%p!=0:
        Np=1
        for x in srange(p):
            Np=Np+((x^3+(a%p)*x)^((p-1)/2)+1)%p
        Ap=p+1-Np
        if (apa(a,p^(n+1))-Ap*apa(a,p^n)+p*apa(a,p^(n-1)))%p^(2*n+1)==0:
            print a.digits(p), Ap

[2, 3] 4 [3, 1, 4] -4 [2, 3, 0, 1] 4 [1, 1, 1, 1] 2 [3, 1, 4, 1] -4 [2, 3, 0, 2] 4 [3, 1, 4, 2] -4 [2, 3, 0, 3] 4 [4, 3, 3, 3] -2 [3, 1, 4, 3] -4 [2, 3, 0, 4] 4 [3, 1, 4, 4] -4

For p=5, we get 12 solutions (no matter what we choose for n).  Notice from the patterns of digits that these all come from 4 solutions: [2,3,0,*], [3,1,4,*], [1,1,1,1], and [4,3,3,3].  

You should read these values as [4,3,3,3]=4+3*5+3*5^2+3*5^3=469.  Also, notice that modulo 5^5, [2,3,0,*]=-[3,1,4,*] and [1,1,1,1]=-[4,3,3,3].  

We could run the code above with larger values of n to see more base-5 digits of these values; but it will take a very long time to run for larger values of n.  A better way to find more digits of these values of a yielding supercongruences is simply to test the 5 possibilities for the next digit.  The code below will do this for a=[1,1,1,1].

a=156
for d in srange(5):
    b=a+d*5^4
    if (apa(b,5^5)-2*apa(b,5^4)+5*apa(b,5^3))%5^9==0:
        print b.digits(5)

[1, 1, 1, 1, 1]

Next, we compute several more digits at once, precomputing the values "binomial((5^n-1)/2,(5^n-1)/4)".  We store these values in the list "B"; to use the nth value, we type "B[n]".  (The values are indexed beginning with 0, so we have B[0] through B[9] below.)

B=[]
for n in srange(10):
    B=B+[binomial((5^n-1)/2,(5^n-1)/4)]
print B[0], B[1], B[2]

1 2 924

a=156
for i in srange(4,9):
    for d in srange(5):
        b=a+d*5^i
        if (B[i+1]*b^((5^(i+1)-1)/4)-2*B[i]*b^((5^i-1)/4)+5*B[i-1]*b^((5^(i-1)-1)/4))%5^(2*i+1)==0:
            print b.digits(5)
            atemp=b
    a=atemp

[1, 1, 1, 1, 1] [1, 1, 1, 1, 1, 1] [1, 1, 1, 1, 1, 1, 1] [1, 1, 1, 1, 1, 1, 1, 1] [1, 1, 1, 1, 1, 1, 1, 1, 1]

It was clear what the pattern was before getting all of these extra digits; but I was not aware of the value of "a" beginning with digits [2,3,0,...] before doing the above calculation.  To understand more about this value of "a", I would find more digits using the process above, hoping to find a repeating pattern.

Since the base-5 digits repeat, a=[1,1,1,1,1,...] can be identified with a rational number (just like the case where the digits in a decimal expansion repeat).  Notice that we have the equation  1+5*a=a.  That is, a=-1/4.  For the elliptic curve y^2=x^3-1/4x, we get supercongruences for all primes p=1 (mod 4).  

You can investigate other values of a that yield supercongruences for the family of curves y^2=x^3+ax; you can also investigate the family of curves y^2=x^3+b, for primes p=1 (mod 6).  For the family y^2=x^3+b, we define a function apb.

def apb(b,n):
# This function returns the coefficient of x^(n-1) in (x^3+b)^((n-1)/2).
    if n%6!=1:
        return 0
    else:
        k=(n-1)/6
        return binomial(3*k,k)*b^k

We also define a function for computing the truncated hypergeometric series 2F1(1/2,1/2,1;-1)_((n-1)/2).

def hyp(n):
# This function returns the truncated hypergeometric series 2F1(1/2,1/2,1;-1)_((n-1)/2).
    if n%2!=1:
        return 0
    else:
        an=0
        u=1
        for i in srange((n-1)/2+1):
            an=an+u
            u=u*((2*i+1)/2)^2*(1/(i+1))^2*(-1)
        return an

for i in srange(4):
    print hyp(5^i)

1 57/64 14902132418929/17592186044416 369784744815754696223693890940151011135950432665713409537446921010974193/441711766194596082395824375185729628956870974218904739530401550323154944

Now we have values for a_n that are fractions, but the denominators are only powers of 2.  So, the congruences still make sense for all primes p not equal to 2.  There are two ways to handle this; one is to use the command "numerator"; the other is to use the command "N.inverse_mod(m)", which returns the multiplicative inverse of N modulo m.  From the examples below, you can see that this truncated hypergeometric sequence satisfies supercongruences.  There are many variations of this truncated hypergeometric series that also satisfy congruences.

for i in srange(1,4):
    print numerator(hyp(5^(i+1))+2*hyp(5^i)+5*hyp(5^(i-1)))%5^(2*i+1)

0 0 0

a5n=[]
for i in srange(5):
    b=hyp(5^i)
    a5n=a5n+[(b*denominator(b)*denominator(b).inverse_mod(5^(2*i+2)))%5^(2*i+2)]
print a5n

[1, 538, 14044, 331722, 5438211]

for i in srange(1,4):
    print (a5n[i+1]+2*a5n[i]+5*a5n[i-1])%5^(2*i+1)

0 0 0