import manin.ellcurve.ellcurve as ec
import manin.plot.graph as G
import math
import manin.rings.arith as arith

P = arith.eratosthenes(10000)

g = G.Graph("bsdprod",(0,0),(13,6), unit=1.5)

def plot(label, C=G.Color(0,0,1)):
    e = ec.EllipticCurve(label)
    N = e.conductor()
    v = e.anlist(10000)

    g.axes()

    f = 1
    j0 = 0
    points = []
    for i in xrange(len(P)):
        p = P[i]
        if N % p != 0:
            f *= float(p)/float(p+1 - v[p])
        if i%50==0:
            points.append(G.Point(math.log(p), 5*f))
        j = int(math.log(p))
        if j > j0:
            g.text((j,-0.5),"$e^{%s}$"%(j))

    g.text((-0.5,5),1)

    g.color(C)
    g.path(points)
    #g.color(0)
    #for p in points:
    #    g.dot(p,0.02)

    g.color(0)
    g.text(points[len(points)-1] + G.Point(0.2,-0.1), "%s: $%s$"%(label,e))
    g.ps()

def go():
    plot("11A",(0,1,0))
    plot("14A",(0,0.5,0))
    plot("37A",(1,0,0))
    plot("37B",(0,0,1))
    plot("389A",(0,1,1))
    plot("5077A",(1,1,0))
    plot("33A",(1,0,0.5))


def plot2(label, C=G.Color(0,0,1), B=1000):
    global P
    e = ec.EllipticCurve(label)
    N = e.conductor()

    print "counting points"
    v = e.anlist(B)

    print "making prime list"
    if B > 10000:
        P = arith.eratosthenes(B)
    print "done making prime list"

    g.axes()

    f = 1
    j = 0
    points = []
    for i in xrange(160):
        k = i/16.0
        x = math.exp(k)
        while j < B and P[j] < x and x < B:
            p = P[j]
            if N % p != 0:
                f *= float(p)/float(p+1 - v[p])
            j += 1
        if j >= B or x >= B:
            break
        points.append(G.Point(k, 5*f))
        if i%16 == 0:
            g.text((i/16,-0.5),"$e^{%s}$"%(i/16))

    g.text((-0.5,5),1)

    g.color(C)
    g.path(points)
    #g.color(0.2)
    #for p in points:
    #    g.dot(p,0.04)

    g.color(0)
    g.text(points[len(points)-1] + G.Point(0.2,-0.1), "%s: $%s$"%(label,e))
    g.ps()

def go2():
    plot2("11A",(0,1,0))
    plot2("14A",(0,0.5,0))
    plot2("37A",(1,0,0))
    plot2("37B",(0,0,1))
    plot2("389A",(0,1,1))
    plot2("5077A",(1,1,0))
    plot2("681B",(0.3,0.7,0.1))    
    plot2("33A",(1,0,0.5))