# Program for computing the adjoint ( = symmetric square) L-value L(ad f,1) for # f a newform, using Dokchitser's algorithm. load("NewformLocalData.sage") # Mandatory input parameters # f: a newform # Optional input parameters # phi: an embedding of the number field Kf that f is defined over into CC, or # alternatively an index to pick from Kf.embeddings(CCC). Default is 0. # globalmin: data of a global minimal twist of f from NewformLocalData; # if not already calculated pass 0 and calculate it here. # CompError (default 20): How many decimal places of accuracy we want # manual_case (default 0): specify a case we can't automatically determine. # manual_case = 1 : supercuspidal of conductor p^e such that the adjoint # also has conductor p^e & Euler factor (1+1/p)^{-1} def AdjointLValue(f, phi=0, globalmin=0, CompError=20,verbose=false, manual_case=0): BitsOfAccuracy = ceil(CompError*log(10)/log(2)) CCC = ComplexField(BitsOfAccuracy) # Quit if f is not a newform; establish phi correctly if type(f) is not sage.modular.modform.element.Newform: print "Error: f not a newform" return Kf = f.base_ring() if phi == 0: phi = Kf.embeddings(CC)[0] # Compute global min and establish its embedding to CCC if globalmin == 0: globalmin = FindMinimalTwistGlobal(f,phi) Nmin = globalmin[2] fmin = globalmin[0] Kfmin = fmin.base_ring() phiminCC = globalmin[3] CCCEmbs = Kfmin.embeddings(CCC) emb_found = false; for i in range (0,len(CCCEmbs)): if abs(CCCEmbs[i](Kfmin.gen()) - phiminCC(Kfmin.gen())) < 1e-10: emb_found = true phimin = CCCEmbs[i] break if emb_found == false: print "Error: Could not match embeddings!" # Set up basic parameters for our modular form / L-function k = f.weight() Nmin = fmin.level() chimin = fmin.character() Nchimin = chimin.conductor() gammafactors = [1,k-1,k] # Set up some variables to be worked with below. # EulerFactorsParameters[p][0] should be the number of (1 - a p^{-s})\1 # terms and then the following entries determine the a's EulerFactorParameters = {} AdjN = 1 # Step 1: Factor the level and go prime-by-prime determining the Euler # factor there and the conductor of the adjoint representation; do this # using the NewformLocalData function Nminfact = list(factor(Nmin)) for n in range (0,len(Nminfact)): p = Nminfact[n][0] e = valuation(Nmin,p) if e == 0: alpha = CCC((phimin(fmin[p]) + sqrt(phimin(fmin[p]^2 - 4*chimin(p)*p^(k-1))))/2) beta = CCC((phimin(fmin[p]) - sqrt(phimin(fmin[p]^2 - 4*chimin(p)*p^(k-1))))/2) EulerFactorParameters[p] = [3,alpha,beta] elif e == 1 and Nchimin % p != 0: # p exactly divides the level, trivial character at p => special AdjN = AdjN*p^2 EulerFactorParameters[p] = [1,1/p] elif e >= 1 and (Nmin/Nchimin) % p != 0: # half-ramified principal series AdjN = AdjN*p^(2*e) EulerFactorParameters[p] = [1,1] elif e > 1 and (Nmin/Nchimin) % p == 0 and e % 2 == 1: # must be supercuspidal, not stable under unramified twist, and # adjoint conductor is e+1. AdjN = AdjN*p^(e+1) EulerFactorParameters[p] = [0] elif manual_case == 1: AdjN = AdjN*p^e EulerFactorParameters[p] = [1,-1] else: print "Error: Euler factor not implemented at bad prime",p return # Step 2: Initiate the L-function we want L = Dokchitser(conductor=AdjN, gammaV=gammafactors, weight=1, eps=1, prec=BitsOfAccuracy) Cutoff = L.num_coeffs()+1; if verbose: print "Bad primes considered; conductor:",AdjN," Coefficients needed:",Cutoff # Step 3: Factor Hecke polynomials at good primes we need. P = list(primes(1,Cutoff+1)) for n in range (0,len(P)): p = P[n] if not EulerFactorParameters.has_key(P[n]): alpha = CCC((phimin(fmin[p]) + sqrt(phimin(fmin[p]^2 - 4*chimin(p)*p^(k-1))))/2) beta = CCC((phimin(fmin[p]) - sqrt(phimin(fmin[p]^2 - 4*chimin(p)*p^(k-1))))/2) EulerFactorParameters[p] = [3,alpha,beta] if verbose: print "Good primes considered..." # Step 4: Compute the inverse of the Euler factor at each prime for as # many coefficients as needed. PrimeCoeffs = {} for n in range (0,len(P)): p = P[n] a = EulerFactorParameters[p][0] M = floor(log(Cutoff)/log(p)) PrimeCoeffs[p] = [1] if a == 0: for i in range (1,M+1): PrimeCoeffs[p].append(0) elif a == 1: alpha = EulerFactorParameters[p][1] for i in range (1,M+1): PrimeCoeffs[p].append(alpha^i) elif a == 3: # If a=3 this means we're in the unramified case and are passed # alpha, beta when the Euler factors are ultimately alpha/beta, # beta/alpha, and 1. alpha = EulerFactorParameters[p][1] beta = EulerFactorParameters[p][2] for i in range (1,M+1): coeff = 0 # The p^i-th coefficient is a sum of powers of (alpha/beta)^j for # -i <= j <= i, and can work out that the multiplicity of each power # given by the formula below... for j in range (-i,i+1): coeff += (1+floor((i-abs(j))/2))*(alpha/beta)^j PrimeCoeffs[p].append(coeff) else: print "Error: Euler factor expansion went wrong at prime",p return if verbose: print "Euler factors expanded at all primes..." # Step 5: Put this together to the full list of coefficients. Coeffs = [] for n in range (1,Cutoff+1): # Factor our integer n F = factor(n) coeff = 1 for i in range (0,len(F)): # for each factor p^m retrieve that coefficient. p = F[i][0] m = F[i][1] ppart = PrimeCoeffs[p][m] coeff *= ppart Coeffs.append(coeff) if verbose: print "Coefficients of L-function expanded out." # Step 6: Initiate the L-function, check the functional equation, and return # the value. L.init_coeffs(Coeffs) print " Functional Equation check (should be small):",L.check_functional_equation() Lval = L(1) return Lval