Q-HG with b=1.sagews

######
###  The functions here work only for hypergeometric data {alpha, beta={1....1}, z=1}, which is defined over Q.


############# Setup
P = Primes()
P = P.unrank
UCF = UniversalCyclotomicField()
E = UCF.gen

def Sdsets(d):
    Sd=divisors(d)
    PS=[]
    NS=[]
    for n in Sd:
        if moebius(d/n)!=0:
            if moebius(d/n) == 1:
                PS.append(n)
            else:
                NS.append(n)
    NN=prod(m^(-m) for m in PS)*prod(m^(m) for m in NS)
    return PS, NS, NN


############


# rising factorial (a)_n
def ps(a,n):
    return prod(a+k for k in srange(n))

# Truncated hypergeometric function F(a1,\cdots,an;1,b2,\cdots,bn;t)_n
def Hyp(A,B,t,n):
    return sum(prod(ps(a,k) for a in A)/prod(ps(b,k) for b in B)*t^k for k in srange(n))



#  Apply the multiplication formula


def SS(d,i):
    PS, NS, NN=Sdsets(d)
    aa=prod(factorial(m*i) for m in PS)/prod(factorial(m*i) for m in NS)
    aa=aa/factorial(i)^(euler_phi(d))*NN^i
    return aa



def trunFF(L,n):
    hp=0
    for i in srange(n):
        aa=1
        for d in L:
            if d==2:
                aa=aa*SS(d,i)^2
            else:
                aa=aa*SS(d,i)
        hp=hp+aa
    return hp



# FF= GF(p^r)FF.multiplicative_generator()

###############
#### Finite Field

def C(p,i,x): #Modify multiplicative characters such that chi(0)=0 for each one including the trivial character
    D=DirichletGroup(p)
    chi=D[i]
    if i==0:
        if x==0:
            return 0
        else:
            return 1
    else:
        return chi(x)

# define Gauss sum
def Gsum(p,i):
    w=E(p)
    return sum(w^t*C(p,i,t) for t in srange(p))


def S(p,d,i):
    D=DirichletGroup(p)
    chi=D[i]
    PS, NS, NN=Sdsets(d)
    aa=prod(Gsum(p,(m*i)%(p-1)) for m in PS)/prod(Gsum(p,(m*i)%(p-1)) for m in NS)
    aa=aa*Gsum(p,(p-1-i)%(p-1))^euler_phi(d)*C(p,i,NN)
    return aa

#def Hp(L,t,p):
#      return sum(prod(S(p,d,i) for d in L)*C(p,i,t) for i in srange(p-1))/(1-p)

def Hp(L,p):
    hp=0
    for i in srange(p-1):
        aa=1
        for d in L:
            if d==2:
                aa=aa*S(p,d,i)^2
            else:
                aa=aa*S(p,d,i)
        hp=hp+aa
    return hp/(1-p)

##############
#using the supercongruence 4F3(r1,1-r1,r2,1-r2;1,1,1,1;1)\equiv ap(f) mod p^3 and |ap(f)|le 2p^{3/2} to determine the ap values
def ap(L,p): 
    x=trunFF(L,p)
    a=x%p^3
    if abs(a)< 2*p^(3/2):
        return a
    else:
        return a-p^3

##############
#Below identify the corresponding newform
#Serrebdd is the upperbound for level, given by Serre, will search all its divisors

def Serrebdd(L):
    N=1
    LL={d for d in L}
    for a in LL:
        dd=factor(a)
        for i in srange(len(dd)):
            d=dd[i][0]
            if d==2:
                N=N*2^8
            else:
                if d==3:
                    N=N*3^5
                else:
                    N=N*d^2
    return N

def CuspformIdentifier(LevelUpper,weight,P,Ap): #LevelUpper is the upperbound for level, will search all its divisors, P is a list of primes and Ap is the list of corresponding ap value
    Candidate=[]
    S=divisors(LevelUpper)
    for  N in S:
        M=Newforms(N,weight,names='a')
        for j in srange(len(M)):
            f=M[j]
            if [f[p] for p in P]==Ap:
                    Candidate.append([N,j])
    return Candidate

#Compute a database of MF, QQ coefficients, weight=4, level divides 2^8*3^5
R=Newforms(8,4)[0].base_ring() ## ser R =Rational Field
MF=[]
Serrebdd=2^8*3^5
S=divisors(Serrebdd)
for  N in S:
        chi=DirichletGroup(N)[0]
        M=Newforms(chi,4,names='a')
        NQ=[f for f in M if f.base_ring()==R]
        if len(NQ)>0:
          MF.append([N,NQ])
MF 

#Example, verifying the supercongruence of F(1/2,1/2,1/3,2/3;1,1,1,1;1)= a_p(f_{36.4.a.a}) \mod p^3
L=[2,3]
f=Newforms(36,4)[0] #f=f36.4.a.a
for p in primes(5,30):
    Z=Zp(p,4)
    p,Z(trunFF(L,p)-f[p]) #print p, F(1,2,1/2,1/3,2/3;1,1,1,1;1)- a_p(f_{36.4.a.a}) in Zp
"end"

(5, 4*5^3 + 5^4 + 4*5^5 + O(5^7)) (7, 7^3 + 5*7^4 + 4*7^5 + 4*7^6 + O(7^7)) (11, 10*11^3 + 8*11^4 + 11^5 + 3*11^6 + O(11^7)) (13, 4*13^3 + 12*13^4 + 5*13^5 + 10*13^6 + O(13^7)) (17, 16*17^3 + 5*17^4 + 9*17^5 + 16*17^6 + O(17^7)) (19, 13*19^3 + 5*19^4 + 2*19^5 + 3*19^6 + O(19^7)) (23, 15*23^3 + 8*23^4 + 2*23^5 + 5*23^6 + O(23^7)) (29, 21*29^3 + 26*29^4 + 7*29^5 + 13*29^6 + O(29^7)) 'end'

## M is the output of the traces of deg-2 Galois rep of HGM  from Magma
M=[ -25, 15, 20, -72, 2, -114, -30, 101, -430, 30, 110, 330, -621, 660, -376,
-250, 360, 785, 488, -489, 450, -1105, -1425, -1060, -1485, -862, -690, 1865 ]

L=[3,3]
for i in srange(3,30):
    p=P(i)
    aa=trunFF(L, p)%(p^3)
    bb= M[i-3]
    cc=Mod(bb-aa,p^3)
    if cc!=0:
        print(p,aa,bb,cc)


#Finding the Hecke cuspidal eigenform (1/3,1/3), inputs are LevelUpper=3^5,weight=4, where ap with p\in {5,7,11, 17, 19} are given as above
PP=[7,11,13,17,19]
Ap=[ap(L,p) for p in PP]
CuspformIdentifier(Serrebdd(L),4,PP,Ap)

[[27, 0]]

#the outcome says the only candidate is the first level 27, weight 4 modular in the list Newforms(27,4)
Newforms(27,4,names='a')[0]

q - 3*q^2 + q^4 - 15*q^5 + O(q^6)

L=[2,3]
for i in srange(3,8):
    p=P(i)
    aa=trunFF(L, p)%(p^3)
    bb= Hp(L,p)
    cc=Mod(bb-(aa+kronecker(3,p)*p),p^3)
    if cc!=0:
        print(p,aa,bb,cc)

Ap=[ap(L,p) for p in PP]

print(Ap)
CuspformIdentifier(Serrebdd(L),4,PP,Ap)

[8, -36, -10, -18, -100]

#for (1/8,3/8) it is the second form in Newforms(128,4),i.e. f_{128.4.a.b}

L=[5]

for i in srange(3,8):
    p=P(i)
    aa=trunFF(L, p)%(p^3)
    bb= Hp(L,p)
    cc=Mod(bb-(aa+kronecker(5,p)*p),p^3)
    if cc!=0:
        print(p,aa,bb,cc)

Ap=[ap(L,p) for p in PP]
CuspformIdentifier(Serrebdd(L),4,PP,Ap)

[6, -43, -28, 91, -35] [[25, 1]]

[[864, 0]]

L=[8]
for i in srange(3,8):
    p=P(i)
    aa=trunFF(L, p)%(p^3)
    bb= Hp(L,p)
    cc=Mod(bb-(aa+kronecker(2,p)*p),p^3)
    if cc!=0:
        print(p,aa,bb,cc)

L=[10]
for i in srange(3,8):
    p=P(i)
    aa=trunFF(L, p)%(p^3)
    bb= Hp(L,p)
    cc=Mod(bb-(aa+kronecker(1,p)*p),p^3)
    if cc!=0:
        print(p,aa,bb,cc)