r"""

Testing L-series for Siegel modular forms
Usage:
 1) copy the correct file to gram_data.wtXX.sage
 2) init_forms(k,prec)
     - init generators and spaces of Siegel modular forms for weight k
       with precision (number of coefficients) prec
 3) set_data()  
    - sets the correct gram matrix to gram_data
 4) getXD()
    - Sets the l-series FD_FG
"""

#attach 'siegel_modular_form.sage'
attach 'siegel_test_fns.sage'
Weight=22
PREC=1000
## Note that in previous versions the loading had to be done outside this file to get correct precision
load 'gd/gram_data.22.sage'
NMAX=6
do_print=True # False
def get_fns22(gd,nmax):
    print "Compute functions for weight 22!"
    #m = Morp([A,B,C,D])
    #R.<a,b,c,d>=PolynomialRing (ZZ)
    #p = -61/80870400*a^3*c - 1/16174080*b^2*c + 1/2695680*a*b*d + c*d
    # Construct Ups22 if it doesn't exist
    # Interesting form
    try:
        F=Ups22=SiegelModularForm("Ups22.sobj")
    except IOError:
        # Ups22: Rational Field
        F=Ups22 = -61/80870400*A^3*C - 1/16174080*B^2*C + 1/2695680*A*B*D + C*D
        f=open("Ups22.sobj","w")
        F.pickle(f)
        f.close()
    print "F is Hecke=",test_hecke_form(F,2)

    # Maass forms:  Number Field in a 
    K.<a>=NumberField(x^3 - 9092496*x^2 + 21143242137600*x - 8184459256292966400)
    try:
        M=SiegelModularForm("M22.sobj")
    except IOError:
        M=(427/1442874228640579584000*a^2 - 33926969/30059879763345408000*a - 16461509/81542642587200)*A^3*C + (-1259/2545230139321982386176*a^2 + 135397057/53025627902541299712*a - 3315701059/2876824430476416)*B^2*C + (1/2992628029773053952*a^2 - 357371/133599465614868480*a + 65130251/15221293282944)*A*B*D + C*D
        fp=open("M22.sobj","w")
        M.pickle(fp)
        fp.close()
    M1 = deepcopy(M)
    M2 = deepcopy(M)
    M3 = deepcopy(M)
    s1, s2, s3 = K.embeddings(RF)
    c1, c2, c3 = M1.coeffs(), M2.coeffs(),M3.coeffs()
    for f in c1:
        c1[f] = s1(M[f])
        c2[f] = s2(M[f])
        c3[f] = s3(M[f])
    # Test that the Maass forms are Hecke
    print "M1 is Hecke=",test_hecke_form(M1,2)
    print "M2 is Hecke=",test_hecke_form(M2,2)
    print "M3 is Hecke=",test_hecke_form(M3,2)

    # Eisenstein Series E: Rational Field
    try:
        E=SiegelModularForm("Eisen22.sobj")
    except IOError:
        E = 591317983484129542246799/878315225547236344903387054080000*A^4*B + 62324003400640902910331/258224676310887485401595793899520*A*B^3 - 28986074191926471709/593235613963048029235200*A^3*C - 13016791830932420887/1560521893993323184272384*B^2*C - 312089362102485871273/7059503806160271547898880*A*B*D + C*D
        f=open("Eisen22.sobj","w")
        E.pickle(f)
        f.close()
    # Klingen Eisenstein Series: Rational Field
    try:
        KL=SiegelModularForm("KEisen22.sobj")
    except IOError:
        KL=-1388543/4218120437760000*A^4*B + 1388543/4218120437760000*A*B^3 + 1317539/33903360000*A^3*C - 45607/2847882240*B^2*C - 657293/16951680000*A*B*D + C*D
        f=open("KEisen22.sobj","w")
        KL.pickle(f)
        f.close()
    print "Eisen22 is Hecke=",test_hecke_form(E,2)
    print "K is Hecke=",test_hecke_form(KL,2)
        
    ## Define the L-series

    # Interesting Spinor L-function
    X_F = X_(F,gd,nmax)  # Spinor  L-fn of F
    # For testing that the Gram matrices are ok
    R_FM1 = R_( F,M1, gd,nmax) # Should be equal to X_F
    R_FM2 = R_( F,M2, gd,nmax) # Should be equal to X_F
    R_FM3 = R_( F,M3, gd,nmax) # Should be equal to X_F
    # Interesting convolution series
    R_FF = R_(F, F, gd,nmax)    
    # Maass Spinor L-functions
    X_M1 = X_(M1, gd,nmax) # Spinor L-fn of M1
    X_M2 = X_(M2, gd,nmax) # Spinor L-fn of M2
    X_M3 = X_(M3, gd,nmax) # Spinor L-fn of M2
    # Eisenstein series
    X_E=X_(E,gd,nmax)
    X_KL=X_(KL,gd,nmax)    
    # Completed interesting convolution series: for printing
    D_FF = D_(F, F, gd,Weight)

    FL=[F]
    # The completed convolution-type L-series: only for printing
    DL=[D_FF]
    # Incomplete Spinor L-functions and conv. L-series: for testing gram matrices 
    TL=[X_F,R_FM1,R_FM2,R_FM3]
    # Incomplete L-series: for checking linear independence 
    BL=[X_M1,X_M2,X_M3,X_E,X_KL,R_FF]
    return [DL,TL,BL,FL]

try:
    [D_FF]=DL
    [X_F,R_FM1,R_FM2,R_FM3]=TL
    [X_M1,X_M2,X_M3,X_E,X_KL,R_FF]=BL
    [F]=FL
except:
    [DL,TL,BL,FL]=get_fns22(GRAM_DATA,NMAX)
    [D_FF]=DL
    [X_F,R_FM1,R_FM2,R_FM3]=TL
    [X_M1,X_M2,X_M3,X_E,X_KL,R_FF]=BL
    [F]=FL

if(do_print):
    print "######### L-series for weight = %s ###########" % (Weight)
    #  number of digits to print
    utd=50
    print "######### D_FF ###########" 
    print_Lcoeff(D_FF,utd)


## Testing (experimentally) the error in the gram matrix. 
## By comparing the spinor L-function of Fj with the convolution L-series D_{FjMj}
if(do_print):
    print "\n"
    print "Comparison: R_FM1 <---> X_F1"
    print_cmp_rel_er(R_FM1,X_F)
    print "Comparison: R_FM2 <---> X_F1"
    print_cmp_rel_er(R_FM2,X_F)
    print "Comparison: R_FM3 <---> X_F1"
    print_cmp_rel_er(R_FM3,X_F)

me1=get_cmp_rel_er(R_FM1,X_F,NMAX)
me2=get_cmp_rel_er(R_FM2,X_F,NMAX)
me3=get_cmp_rel_er(R_FM3,X_F,NMAX)
maxerr=max(me1,me2,me3)
print "Maxerror(difference)=",mppr(maxerr)
eps1=max(MAXERR.values())
#eps1=mpmath.mpf(1E-49)
n1=L2_norm(F,GRAM_DATA,NMAX)
print "max_norm=",mppr(n1)
print "dps=",mpmath.mp.dps
norm=n1**2
eps2=norm*eps1
print "eps1=",mppr(eps1)
print "eps2=",mppr(eps2)
print "Testing all:::::"
T=check_total_indep(BL,eps2)
print T[0],mppr(T[1])

## The true ep2
eps2=proper_eps2([F],[R_FF],GRAM_DATA,NMAX,eps1)
print "True eps2=",mppr(eps2)
T=check_total_indep(BL,eps2)
print T[0],mppr(T[1])