twistedlseries.sage

from psage.number_fields.sqrt5.misc import *
from psage.lseries.eulerprod import LSeriesEllipticCurve, LSeriesEllipticCurveSqrt5, LSeriesTwist

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

# Given an element of a number field and a modulus prod primes returns
def rep(alpha, primes):
    return [p.ideallog(alpha)[0] for p in primes]

# Hecke characters of number fields
class NFChar():
    # Initialise a new Hecke character of order order, factored conductor should be a list of primes,
    # rootlist a list of integers of the same length that describe the power of zeta that the gene
    def __init__(self, factored_conductor, rootlist, order, zeta = None, prec = None):
        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"
        assert len(rootlist) == len(factored_conductor), "rootlist and factored_conductor should be the same length"
        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+[order])==1, "the order of the character is not %s"%(order)
        self._field = factored_conductor[0].number_field()
        self._ring = self._field.maximal_order()
        self._degree = self._field.absolute_degree()
        self._order = order
        if prec == None:
            self._CC = CC
        else:
            self._CC = ComplexField(prec)
        if zeta is None:
            self._zeta = self._CC.zeta(order)
        else:
            #assert zeta**d==1,"order of zeta does not appear to be %s"%order
            self._zeta = zeta

        assert self._field.class_number()==1,"Class number of base field is not 1"

        for p in factored_conductor:
            assert ZZ(order).divides(p.norm() - 1), "order %s does not divide order of (O/%s)^*, which is %s"%(order,p,p.norm()-1)

        self._pows = [self._zeta**d for d in range(order)]

        if False:
            self._logs = dict()
            for p in self._factored_conductor:
                self._logs[p] = dict()
                kp = p.residue_field()
                for x in kp:
                    if x == 0:
                        continue
                    self._logs[p][x] = p.ideallog(kp.lift(x))[0]

        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):
        p = self._field.ideal(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 _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])
        if False: # new 
            reps = [self._logs[p][p.residue_field()(alpha)] for p in fact]
        reps = rep(alpha, [prime for prime in fact])

        return self._pows[sum([rootlist[i]*reps[i] for i in range(len(rootlist))]) % self._order]

    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 conjugate(self):
        return NFChar(self._factored_conductor, [-d for d in self._rootlist], self._order, self._zeta, self._CC.prec())

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

    def conductor(self):
        K = self._factored_conductor[0].number_field()
        return ZZ(prod(self._factored_conductor).norm())#*K.discriminant())

    def order(self):
        return self._order

    def __repr__(self):
        return "Hecke character of modulus %s and order %s"%(self._conductor,self._order)

# Search for all characters of the base field of E (which for now must be Qsqrt5) that are 
# of order D and have conductor norm less than bound, for each pair of complex conjugate characters, only 1 is used.
# When a character is found the twist is computed (and optionally the functional equation checked)
# then the twist is evaluated at 1 and the result is printed and added to a list, which is returned.
def value_search(E, D, bound, checkfe=False, lower=0, prec=None):
    LE = LSeriesEllipticCurveSqrt5(E)
    omega = prod([E.period_lattice(phi).basis(prec=1000)[0] for phi in F.embeddings(RR)])
    ideals = E.base_field().ideals_of_bdd_norm(bound)
    chis = []
    #print len(ideals),ideals
    if prec is None:
        one = CC(1)
    else:
        one  = ComplexField(prec)(1)
    # we could also use for m in E.base_field().primes_of_bounded_norm_iter(bound):
    for nm in ideals:
        if nm < lower:
            continue
        if nm % 5 == 0 or nm % 31 == 0:
            continue
        for m in ideals[nm]:
            if nm == 1:
                continue
            fact = [prime for (prime, i) in m.factor()]
            if max([i for (prime, i) in m.factor()]) > 1:
                continue
            for t in Integers(D)^len(fact):
                try:
                    if 0 in list(t):
                        continue
                    if t[0].lift() > D/2: # Only take one character in each conjugacy class
                        continue
                    rootlist = list(t.lift())
                    chi = NFChar(fact,rootlist,D, CC.zeta(D))
                    print m,m.norm(),m.factor(), rootlist
                    chis += [chi]

                except AssertionError: # Character didn't work, try the next
                    pass
                except RuntimeError:
                    print "Error, maybe conductor was wrong?!"

    print len(chis),"characters to be evaluated"
    # "prime" the base L-series by computing the biggest character so as much as possible is cached at the start
    try_chi(LE,chis[-1],omega,checkfe=checkfe,one=one)
    # TODO don't recompute this in the parallel computation!
    out = [(inp[0][1],)+out for inp,out in try_chi([((LE,chi,omega),{"checkfe":checkfe,"one":one}) for chi in chis])]
    # parallel might have reordered chis so we need outchis too
    outchis, values, n_values = zip(*out)
    return (values, n_values, outchis)

@parallel#(p_iter="reference")
def try_chi(LE, chi, omega, one=None,checkfe=False):
    LEchi = LE.twist(chi,epsilon='solve')
    if one == None:
        one = CC(1)
    if checkfe:
        LEchi.check_functional_equation(1.2)
        LEchi.check_functional_equation(1.1)
    val = LEchi(one)
    print "raw value",val
    normalised = val*chi.conjugate().gauss_sum()*CC(sqrt(5))/omega
    print "normalised",normalised
    print "algdeprts",algdep(normalised, 2, known_bits=53).roots(CC)
    return (val, normalised)