set_verbose(ZZ(2))



# Compute B_{m,\chi} for some character chi with conductor f to within eps

def bern_complex(m, chi, eps):
    """
    Return Bernoulli B_{m,chi}.
    """
    #prec = 16*m;
    chi = chi.bar()
    prec = ZZ(500)
    C = ComplexField(prec)
    R = RealField(prec)
    eps = R(eps)
    i = C.gen(0)
    if not chi.is_primitive():
        raise ValueError, "chi (=%s) must be primitive"%chi
    f = chi.conductor()
    mfact = factorial(m)
    K = (-ZZ(1))**(m-ZZ(1)) * ZZ(2) * mfact * C(f/(ZZ(2)*i))**C(m)
    P = R(pi)
    PM = P**m
    GS = chi.gauss_sum_numerical(prec)
    delta = abs((PM*GS)/K)*eps    # accuracy that we have to compute L(chi,m) to.
    N = int( (delta**(-ZZ(1)/(m-ZZ(1)))).ceil() )
    verbose("N = %s"%N)
    phi = chi.base_ring().complex_embedding(prec)
    L = sum([phi(chi(n))/(n**m) for n in range(ZZ(1),N)])
    verbose("L = %s"%L)
    verbose('Lbar = %s'%sum([phi(chi(n))/(n**m) for n in range(ZZ(1),N)]))
    return K/(PM*GS)*L

def bern_charpoly(m, chi, eps):
    G = chi.galois_orbit()
    r = [bern_complex(m, e, eps) for e in G]
    C = parent(r[0])
    C = ComplexField(ZZ(53))
    R, x = C['x'].objgen()
    f = prod([x - a for a in r])
    return f
    
def test_bern_charpoly(n=ZZ(11), m=ZZ(5), eps=RR('0.001')):
    chi = DirichletGroup(n).gen(0)
    f0 = bern_charpoly(m,chi,eps)
    f1 = chi.bernoulli(m).charpoly()
    return f0, f1, chi


### Test It Out
def test_bern_complex(n=ZZ(5), m=ZZ(5), eps=RR('0.001')):
    chi = DirichletGroup(n).gen(0)
    B = bern_complex(m,chi,eps)
    B0 = chi.bernoulli(m)
    phi = parent(B0).complex_embedding()
    B1 = phi(B0)
    print "B = ", B
    print "B1 = ",B1
    err = abs(CC(B) - B1)
    print err<eps