twistedlseries-maarten.sagews

from psage.eulerprod import LSeriesEllipticCurve, LSeriesTwist

def rep(alpha, primes):
    vals = []
    for p in primes:
        vals += [p.ideallog(alpha)[0]]
    return vals

def _iter__OK_mod_p(OK,P):
    M = P.pari_hnf().sage().diagonal()
    from itertools import product
    reps = product(*map(range,M))
    B = OK.basis()
    return [sum(i*j for i,j in zip(r,B)) for r in reps]


def totally_positive_generator(p):
    K = p.ring()
    units = K.unit_group().gens()
    x = p.gens_reduced()
    assert len(x)==1, "p not pricipally generated!"
    x = x[0]
    M = Matrix(GF(2),[[phi(g)<0 for phi in K.real_embeddings()] for g in units])
    v = vector(GF(2),[phi(x)<0 for phi in K.real_embeddings()])
    x *= prod(g for g,e in zip(units,M.solve_left(v)) if e)
    assert x.is_totally_positive(), "Totally positive generator does not exist!"
    return x

class NFChar():
    def __init__(self, factored_conductor, rootlist, order, zeta = None, prec = None):
        """
        Currently we allow only tamely ramified characters.
        """
        assert all(p.is_prime() for p in factored_conductor), "factored_conductor contains an element that is not prime"
        assert len(set(factored_conductor)) == len(factored_conductor), "factored_conductor contains duplicate elements"
        self._factored_conductor = factored_conductor
        self._primes_conductor = [p.smallest_integer() for p in factored_conductor]
        self._conductor = prod(factored_conductor)
        self._m = totally_positive_generator(self._conductor)
        self._rootlist = rootlist
        assert gcd(rootlist+[d])==1, "the order of the character is not %s"%(order)
        self._order = order
        self._field = factored_conductor[0].number_field()
        self._ring = self._field.maximal_order()
        self._degree = self._field.absolute_degree()
        if prec == None:
            self._CC = CC
        else:
            self._CC = ComplexField(prec)
        if zeta is None:
            self._zeta = self._CC.zeta(d)
        else:
            #assert zeta**d==1
            self._zeta = zeta
            
        assert self._field.class_number()==1
        
        for g in self._field.unit_group().gens():
            assert self._call_nf_elt(self._field(g)) == 1, " character is not 1 on units"

    def __call__(self, p):
        x = p.gens_reduced()
        assert len(x)==1, "Do not know how to evaluate on a non principal ideal!"
        return self._call_nf_elt(p.gens_reduced()[0])

    @cached_method
    def _iter_rep_mod_m(self):
        return _iter__OK_mod_p(self._ring,self._conductor)
    

    
    def gauss_sum(self):
        assert self._ring.discriminant()==5, "only implemented for sqrt5"
        tau = 0
        sqrt5 = (self._field*5).factor()[0][0].gens_reduced()[0]
        print sqrt5
        for alpha in self._iter_rep_mod_m():
            #print tau
            a = alpha/self._m#*sqrt5
            a = a.trace()
            den = a.denominator()
            num = a.numerator()
            tau += self._call_nf_elt(alpha)*self._CC.zeta(den)^(num%den)
        return tau

    def _call_nf_elt(self,alpha):
        fact = self._factored_conductor
        rootlist = self._rootlist

        for p,q in zip(self._primes_conductor,fact):
            if p.divides(alpha.norm()) and alpha in q:
                return 0
        #return m.ideallog(p.gens_reduced()[0])
        reps = rep(alpha, fact)
        #print "#rootlist %s, #reps %s"%(len(rootlist),len(reps))
        return self._zeta**(sum([rootlist[i]*reps[i] for i in range(len(rootlist))]) % self._order)


    
    def is_primitive(self):
        return all(i%self._order != 0 for i in self._rootlist)

    def conductor(self):
        return prod(self._factored_conductor).norm()

    def order(self):
        return self._order

R.<x> = QQ[]
K.<phi>=NumberField(x^2 -x-1)
E = EllipticCurve(K, [1, phi + 1, phi, phi, 0])
m = K*13
fact = [prime for (prime, i) in m.factor()]
d = 3

foundlist = len(fact)*[1]
for t in Integers(d)^len(m.factor()):
    try:
        if 0 in list(t):
            continue
        rootlist = list(t.lift())
        chi = NFChar(fact,rootlist,d,prec=1000)
        foundlist = rootlist

        print rootlist
        x = chi.gauss_sum();
        print
        print x.abs(),x
    except:
        pass

[1] -2*phi + 1 13.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 11.3407561726636152699368162869855117983731253191903113812827727652124965276757522368493039894057371978570715000611418213033343214441382821742420971004022136131129716550720474800350052564067483123669351089814797540077428600188606608171977226854509076353420967694562397343116132495426591104271695619234 + 6.35509633539823757082843897576929635164123227418046373191723445222979581773290744537857331917542382091652147604151185486267441580032280296429497772415175454074114827575353571701097120106343964546518212091823462667259056273414807931764641831386718012747116711234367738464348211389670063607384427700001*I [2] -2*phi + 1 13.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 11.3407561726636152699368162869855117983731253191903113812827727652124965276757522368493039894057371978570715000611418213033343214441382821742420971004022136131129716550720474800350052564067483123669351089814797540077428600188606608171977226854509076353420967694562397343116132495426591104271695619234 - 6.35509633539823757082843897576929635164123227418046373191723445222979581773290744537857331917542382091652147604151185486267441580032280296429497772415175454074114827575353571701097120106343964546518212091823462667259056273414807931764641831386718012747116711234367738464348211389670063607384427700001*I

x

11.3407561726636152699368162869855117983731253191903113812827727652124965276757522368493039894057371978570715000611418213033343214441382821742420971004022136131129716550720474800350052564067483123669351089814797540077428600188606608171977226854509076353420967694562397343116132495426591104271695619234 - 6.35509633539823757082843897576929635164123227418046373191723445222979581773290744537857331917542382091652147604151185486267441580032280296429497772415175454074114827575353571701097120106343964546518212091823462667259056273414807931764641831386718012747116711234367738464348211389670063607384427700001*I

CC1 = ComplexField(1000)
I = CC1.0
omega = prod([E.period_lattice(phi).basis(prec=1000)[0] for phi in K.embeddings(RR)])

#magma LValues for cubic twist of conductor 13 prec e-74

l1,l2 =[ 1.5457956238192788683557199097437729233665818071610358735968195553567568888501269808125460910379551385036304017986579868697753000481659891333945500875045540441514902685477548634003391025711156928736075 -
0.86622796175514932234718383126453984183577137642394022731205644771173084410073936907261536775497232496651092083801009064197224548531397645425498120806474732224747486052910601468664988353184945233634762*I,
1.5457956238192788683557199097437729233665818071610358735968195553567568888501269808125460910379551385036304017986579868697753000481659891333945500875045540441514902685477548634003391025711156928736077 +
0.86622796175514932234718383126453984183577137642394022731205644771173084410073936907261536775497232496651092083801009064197224548531397645425498120806474732224747486052910601468664988353184945233634753*I ]

13.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

#mamgma LValues for cubic twist of conductor 47 prec e74
l1,l2 = [
0.24505801487087417016068854563803448172521943510270968636805221995783706098457669244070126558322132754114006237893747704077409789994033154666761763371421995661540622180222286491609840322732421598650475 +
0.42445293274705574412757500887055189552420282612221945972766975611605748380421821057087580584770905794678215958395810249645747226385011392829719555699234058808465069312539939557361013672581656837204405*I,
0.24505801487087417016068854563803448172521943510270968636805221995783706098457669244070126558322132754114006237893747704077409789994033154666761763371421995661540622180222286491609840322732421598650477 -
0.42445293274705574412757500887055189552420282612221945972766975611605748380421821057087580584770905794678215958395810249645747226385011392829719555699234058808465069312539939557361013672581656837204386*I ]

omega = prod([E.period_lattice(phi).basis(prec=1000)[0] for phi in K.embeddings(RR)])

c1 = l1*47*CC1(sqrt(5))/omega
c2 = l2*47.conjugate() *CC1(sqrt(5))/omega

c1,c2

(0.9999999995653503062529576591124754896370504363434188055463447835913000245126430422276394574958995879729683943441266440810124106473200871203105837568300491479678574046284886618745443521709166986551736 + 1.732050807586195067941157851961213937336650671437163840656006336307774130379298105089823768284308618777381632838540174758579615723491693438187234643639901697333509474185442000836983390248527316165499*I, 0.9999999995653503062529576591124754896370504363434188055463447835913000245126430422276394574958995879729683943441266440810124106473200871203105837568300491479678574046284886618745443521709166986551737 - 1.732050807586195067941157851961213937336650671437163840656006336307774130379298105089823768284308618777381632838540174758579615723491693438187234643639901697333509474185442000836983390248527316165498*I)

x.conjugate()

11.3407561726636152699368162869855117983731253191903113812827727652124965276757522368493039894057371978570715000611418213033343214441382821742420971004022136131129716550720474800350052564067483123669351089814797540077428600188606608171977226854509076353420967694562397343116132495426591104271695619234 - 6.35509633539823757082843897576929635164123227418046373191723445222979581773290744537857331917542382091652147604151185486267441580032280296429497772415175454074114827575353571701097120106343964546518212091823462667259056273414807931764641831386718012747116711234367738464348211389670063607384427700001*I

1/c1;1/c2;

0.4999999999999999999999999999999999999999999999999999999999999999999999999934956231565666574343088924407419896235212233272735000147214157716184640643318561169074097207246090575533634257256677689681113 - 9.789279086055920531037253731441653548566460090198335763674624714507240225895177402992985655384451524612420671959350614999606032086253312249764159225105263460065631461269354850272811495630508557477052e-75*I 0.2610221927089224350365631712415115144995162664082569178958147777508118762544302932184219095265161745286366069433103338206566987549141885123237098148801496534879660519721023329528268952448575047226335 - 0.4264591597250857028365957388617469633260920057454412677493119493117616729542769917171234683549142379885358787137148975693919254844183663218964443829974242399178036667969166857387601436367253872148488*I

c1, c2;

(0.9999999995653503062529576591124754896370504363434188055463447835913000245126430422276394574958995879729683943441266440810124106473200871203105837568300491479678574046284886618745443521709166986551736 + 1.732050807586195067941157851961213937336650671437163840656006336307774130379298105089823768284308618777381632838540174758579615723491693438187234643639901697333509474185442000836983390248527316165499*I, 0.9999999995653503062529576591124754896370504363434188055463447835913000245126430422276394574958995879729683943441266440810124106473200871203105837568300491479678574046284886618745443521709166986551737 - 1.732050807586195067941157851961213937336650671437163840656006336307774130379298105089823768284308618777381632838540174758579615723491693438187234643639901697333509474185442000836983390248527316165498*I)

c1*c2-4;c1+c2-2

-8.093088569877350958144983871293747909312945882759218380628346650543790566425014413834834802396732744152051493872087143862134447449917591500152645793415664931919724196162966488372064625571333503992959e-10 + 9.144941236411822604101211875669864313850200804989490725343365823468150650983299713086503530538041448730974292603275057833515443317510437438730834983953523435630821104602005456475703659443420968491774e-199*I -8.692993874940846817750490207258991273131623889073104328173999509747139155447210850082008240540632113117467118379751787053598257593788324863399017040642851907430226762509112956581666026896526704936875e-10 + 7.707879042118536194885307152350314207388026392776856468503694051208869834400209758172910118596349221073249760908474691602534445081901654412644560915046541181460263502450261741886664512959454816300209e-199*I

for d in range(1,13):
    algdep(c1,d,known_digits=10).factor()

1 x^2 - 2*x + 4 x^2 - 2*x + 4 x^2 - 2*x + 4 x^2 - 2*x + 4 x^2 - 2*x + 4 x^2 - 2*x + 4 x^2 - 2*x + 4 x^2 - 2*x + 4 x^2 - 2*x + 4 x^2 - 2*x + 4 x^2 - 2*x + 4

for d in range(1,13):
    algdep(c2,d,known_digits=60).factor()

x - 2 x - 2 x - 2 x - 2 x - 2 x - 2 x - 2 x - 2 x - 2 x - 2 x - 2 x - 2

R.<x> = QQ[]
K.<phi>=NumberField(x^2 -x-1)
E = EllipticCurve(K, [1, phi + 1, phi, phi, 0])
m = K*47
fact = [prime for (prime, i) in m.factor()]
d = 3

foundlist = len(fact)*[1]
for t in Integers(d)^len(m.factor()):
    try:
        if 0 in list(t):
            continue
        rootlist = list(t.lift())
        chi = NFChar(fact,rootlist,d)
        foundlist = rootlist

        print rootlist
        x = chi.gauss_sum();
        print
        print x.abs(),x
    except:
        pass

[1] -2*phi + 1 47.0000000000000 47.0000000000000 + 1.30007116183606e-13*I [2] -2*phi + 1 46.9999999999999 46.9999999999999 + 3.00870439673417e-14*I

newchi = NFChar(fact, foundlist, d)
newchi(K.prime_above(2))

1.00000000000000

x = newchi.gauss_sum();
print x.norm(),x

169.000000000000 11.3407561726636 - 6.35509633539824*I

p = K.prime_above(2)
K = p.ring()
units = K.unit_group().gens()
x = p.gens_reduced()
print x
assert len(x)==1, "p not pricipally generated!"
x = x[0]
M = Matrix(GF(2),[[phi(g)<0 for phi in K.real_embeddings()] for g in units])
v = vector(GF(2),[phi(x)<0 for phi in K.real_embeddings()])
print v
x *= prod(g for g,e in zip(units,M.solve_left(v)) if e)
print x
assert x.is_totally_positive(), "Totally positive generator does not exist!"

Error in lines 1-1 Traceback (most recent call last): File "/projects/sage/sage-7.5/local/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 995, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> NameError: name 'K' is not defined

K

Error in lines 1-1 Traceback (most recent call last): File "/projects/sage/sage-7.5/local/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 995, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> NameError: name 'K' is not defined

M.solve_left(v)

(0, 1)

v;x

(1, 0) 2*phi

totally_positive_generator??

   File: /projects/68c8b2b8-03ba-44d4-a0d1-5d771c8cb465/code
Unable to read source code (source code not available)

totally_positive_generator(K.prime_above(2))

2 (0, 0) (0, 0) 2

L = LSeriesEllipticCurve(E)
print "done1"
LEchi =  LSeriesTwist(L, newchi, epsilon='solve', conductor=ZZ(5^2*E.conductor().norm()*(newchi.conductor()^2)))
print E.conductor().norm()
print newchi.conductor()
LEchi
print "done2"

1.00000000000000 done1 31 169 Twist of L-series of Elliptic Curve defined by y^2 + x*y + phi*y = x^3 + (phi+1)*x^2 + phi*x over Number Field in phi with defining polynomial x^2 - x - 1 by <__builtin__.NFChar instance at 0x7fd03418dd40> done2

LEchi.anlist(100)
L.conductor()
L.epsilon()
L.number_of_coefficients(2)
L.degree()
L.hodge_numbers()
L(3/2)
LEchi.number_of_coefficients(5)
LEchi.degree()
LEchi.hodge_numbers()
LEchi(ComplexField(5)(3/2))

[0, 1, -0.000000000000000, -0.000000000000000, -3.00000000000000, 2.00000000000000, 0.000000000000000, -0.000000000000000, -0.000000000000000, 2.00000000000000, -0.000000000000000, 0.000000000000000, 0.000000000000000, 0, 0.000000000000000, -0.000000000000000, 5.00000000000000, 0, -0.000000000000000, 0.000000000000000, -6.00000000000000, 0.000000000000000, 0.000000000000000, 0, 0.000000000000000, -1.00000000000000, 0.000000000000000, -0.000000000000000, 0.000000000000000, -4.00000000000000, 0.000000000000000, -7.00000000000000, -0.000000000000000, 0.000000000000000, 0.000000000000000, -0.000000000000000, -6.00000000000000, 0, 0.000000000000000, 0.000000000000000, -0.000000000000000, 12.0000000000000, 0.000000000000000, 0, 0.000000000000000, 4.00000000000000, 0.000000000000000, 0, -0.000000000000000, 2.00000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, -8.00000000000000, 0.000000000000000, 4.00000000000000, 0.000000000000000, -0.000000000000000, -3.00000000000000, 0.000000000000000, 0.000000000000000, 0, 0.000000000000000, 0.000000000000000, 0.000000000000000, 8.00000000000000, -0.000000000000000, 0, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 16.0000000000000, 10.0000000000000, -5.00000000000000, -0.000000000000000, 0, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, -4.00000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0.000000000000000, 0, -0.000000000000000, 0.000000000000000, 3.00000000000000] 775 1 534 2 [0, 0, 1, 1] 0.607704196335035 99729 2 [0, 0, 1, 1]

#better
def is_prime(n):
    """returns True if n is a prime false otherwise"""
    if n<2:
        return False
    
    for i in range(2,floor(sqrt(n))+1):
        if n%i == 0:
            return False
    
    return True

def primes_in_range(a,b):
    """returns the list of all primes in the range(a,b)"""
    return [i for i in range(a,b) if is_prime(i)]

print primes_in_range(2,20)

[2, 3, 5, 7, 11, 13, 17, 19]

L = LSeriesEllipticCurve(E)
import cProfile
t = cputime()
cProfile.runctx("L.anlist(100000)",None, locals(), filename="timing_data/profile_test", sort='cumulative')
print cputime(t)

478.148 478.148

LEchi =  LSeriesTwist(L, newchi, epsilon='solve', conductor=ZZ(5^2*E.conductor().norm()*(newchi.conductor()^2)))
t = cputime()
cProfile.runctx("LEchi.anlist(100000)",None, locals(), filename="timing_data/profile_test1", sort='cumulative')
print cputime(t)

26.52

print cputime(t)

L = LSeriesEllipticCurve(E)
L.anlist(100)
print "L done"
LEchi =  LSeriesTwist(L, newchi, epsilon='solve', conductor=ZZ(5^2*E.conductor().norm()*(newchi.conductor()^2)))
LEchi.anlist(100)

%sh
python gprof2dot/gprof2dot.py -f pstats timing_data/profile_test | dot -Tpng -o timing_data/profile_LSeriesEllipticCurve_anlist100000_cachetest_opt.png 2> /dev/null

%sh
python gprof2dot/gprof2dot.py -f pstats timing_data/profile_test1 | dot -Tpng -o timing_data/profile_LSeriesTwist_anlist100000_cachtest_opt.png 2> /dev/null

L

L-series of Elliptic Curve defined by y^2 + x*y + phi*y = x^3 + (phi+1)*x^2 + phi*x over Number Field in phi with defining polynomial x^2 - x - 1

︠96c00861-087e-464f-b834-00be02de9d9f︠