CoCalc Public Fileswww / talks / bernoulli / genber.sage
Author: William A. Stein
1
2def genber_numerical(n, chi, prec, primes):
3    if not chi.is_primitive():
4        raise ValueError, "chi (=%s) must be primitive"%chi
5    CC = ComplexField(prec)
6    K = genber_K(n, chi, CC)
7    a = CC.pi()^n
8    g = chi.gauss_sum_numerical(prec)
9    L = genber_L(n, chi^(-1), CC, primes)
10    d = genber_d(n, chi, prec)
11    print "d = ", d
12    print "K = ", K
13    print "L = ", L
14    print "a = ", a
15    print "g = ", g
16    num = d*K*L/(a*g)
17    return num
18
19def genber_d(n, chi, prec):
20    f = chi.conductor()
21    if f == 1:
22        raise NotImplementedError
23    F = f.factor()
24    assert len(F) >= 1
25    if len(F) > 1:
26        return 1
27    elif f == 4:
28        return 2
29    P = F[0]
30    if P[0] == 2 and P[1] > 2:
31        return 1
32    elif P[1] == 1:
33        return n*P[0]
34    elif P[1] > 1:
35        R = chi.base_ring()
36        phi = R.complex_embedding(prec)
37        return phi(1 - chi(1+p))
38
39
40def genber_K(n, chi, CC):
41    f = chi.modulus()
42    i = CC.0
43    z = 2*factorial(n)*(f/(2*i))^n
44    if n % 2 == 0:
45        z *= -1
46    return z
47
48def genber_L(n, chi, CC, max_prime):
49    R = chi.base_ring()
50    m = chi.modulus()
51    phi = R.complex_embedding(CC.prec())
52    n = phi(n)
53    z = 1
54    zero = 2^(-CC.prec())
55    z_prev = z
56    for p in prime_range(max_prime):
57        h = z * phi(chi(p)) / (CC(p)^n)
58        z = z - h
59        print p, abs(z - z_prev)
60        z_prev = z
61    return 1/z
62
63def test1(prec, primes):
64    e = DirichletGroup(13).0
65    print "group: ", e.parent()
66    print "chi = ", e
67    B = genber_numerical(389, e, prec, primes)
68    print "B3 = ", B
69    return e, B
70
71
72