The Matrix of Frobenius on Hyperelliptic Curves¶
Sage has a highly optimized implementation of the Harvey-Kedlaya algorithm for computing the matrix of Frobenius associated to a curve over a finite field. This is an implementation by David Harvey, which is GPL’d and depends only on NTL (a C library in Sage for fast arithmetic in \((\ZZ/n\ZZ)[x]\)).
We import the hypellfrob function and call it on a polynomial over \(\ZZ\).
sage: from sage.schemes.hyperelliptic_curves.hypellfrob import hypellfrob
sage: R.<x> = PolynomialRing(ZZ)
sage: f = x^5 + 2*x^2 + x + 1; p = 101
sage: M = hypellfrob(p, 1, f); M
[ O(101) O(101) 93 + O(101) 62 + O(101)]
[ O(101) O(101) 55 + O(101) 19 + O(101)]
[ O(101) O(101) 65 + O(101) 42 + O(101)]
[ O(101) O(101) 89 + O(101) 29 + O(101)]
>>> from sage.all import *
>>> from sage.schemes.hyperelliptic_curves.hypellfrob import hypellfrob
>>> R = PolynomialRing(ZZ, names=('x',)); (x,) = R._first_ngens(1)
>>> f = x**Integer(5) + Integer(2)*x**Integer(2) + x + Integer(1); p = Integer(101)
>>> M = hypellfrob(p, Integer(1), f); M
[ O(101) O(101) 93 + O(101) 62 + O(101)]
[ O(101) O(101) 55 + O(101) 19 + O(101)]
[ O(101) O(101) 65 + O(101) 42 + O(101)]
[ O(101) O(101) 89 + O(101) 29 + O(101)]
from sage.schemes.hyperelliptic_curves.hypellfrob import hypellfrob R.<x> = PolynomialRing(ZZ) f = x^5 + 2*x^2 + x + 1; p = 101 M = hypellfrob(p, 1, f); M
We do the same calculation but in \(\ZZ/101^4\ZZ\), which gives enough precision to recognize the exact characteristic polynomial in \(\ZZ[x]\) of Frobenius as an element of the endomorphism ring. This computation is still very fast, taking only a fraction of a second.
sage: M = hypellfrob(p, 4, f) # about 0.25 seconds
sage: M[0,0]
91844754 + O(101^4)
>>> from sage.all import *
>>> M = hypellfrob(p, Integer(4), f) # about 0.25 seconds
>>> M[Integer(0),Integer(0)]
91844754 + O(101^4)
M = hypellfrob(p, 4, f) # about 0.25 seconds M[0,0]
The characteristic polynomial of Frobenius is \(x^4 + 7x^3 + 167x^2 + 707x + 10201\), which determines the \(\zeta\) function of the curve \(y^2= f(x)\).
sage: M.charpoly()
(1 + O(101^4))*x^4 + (7 + O(101^4))*x^3 + (167 + O(101^4))*x^2 + (707 + O(101^4))*x + 10201 + O(101^4)
>>> from sage.all import *
>>> M.charpoly()
(1 + O(101^4))*x^4 + (7 + O(101^4))*x^3 + (167 + O(101^4))*x^2 + (707 + O(101^4))*x + 10201 + O(101^4)
M.charpoly()
Note
Exercise: Write down zeta function explicitly, count points over some finite fields and see that things match up.
Modular Symbols¶
Modular symbols play a key role in algorithms for computing with modular forms, special values of \(L\)-functions, elliptic curves, and modular abelian varieties. Sage has the most general implementation of modular symbols available, thanks to work of myself, Jordi Quer (of Barcelona) and Craig Citro (a student of Hida). Moreover, computation with modular symbols is by far my most favorite part of computational mathematics. There is still a lot of tuning and optimization work to be done for modular symbols in Sage, in order for it to be across the board the fastest implementation in the world, since my Magma implementation is still better in some important cases.
Note
I will talk much more about modular symbols in my lecture on Friday, which will be about modular forms and related objects.
We create the space \(M\) of weight \(4\) modular symbols for a certain congruence subgroup \(\Gamma_H(13)\) of level \(13\). Then we compute a basis for this space, expressed in terms of Manin symbols. Finally, we compute the Hecke operator \(T_2\) acting on \(M\), find its characteristic polynomial and factor it. We also compute the dimension of the cuspidal subspace.
sage: M = ModularSymbols(GammaH(13,[3]), weight=4)
sage: M
Modular Symbols space of dimension 14 for Congruence Subgroup
Gamma_H(13) with H generated by [3] of weight 4 with sign 0
over Rational Field
sage: M.basis()
([X^2,(0,4)],
[X^2,(0,7)],
[X^2,(4,10)],
[X^2,(4,11)],
[X^2,(4,12)],
[X^2,(7,3)],
[X^2,(7,5)],
[X^2,(7,6)],
[X^2,(7,7)],
[X^2,(7,8)],
[X^2,(7,9)],
[X^2,(7,10)],
[X^2,(7,11)],
[X^2,(7,12)])
sage: factor(charpoly(M.T(2)))
(x - 7) * (x + 7) * (x - 9)^2 * (x + 5)^2
* (x^2 - x - 4)^2 * (x^2 + 9)^2
sage: dimension(M.cuspidal_subspace())
10
>>> from sage.all import *
>>> M = ModularSymbols(GammaH(Integer(13),[Integer(3)]), weight=Integer(4))
>>> M
Modular Symbols space of dimension 14 for Congruence Subgroup
Gamma_H(13) with H generated by [3] of weight 4 with sign 0
over Rational Field
>>> M.basis()
([X^2,(0,4)],
[X^2,(0,7)],
[X^2,(4,10)],
[X^2,(4,11)],
[X^2,(4,12)],
[X^2,(7,3)],
[X^2,(7,5)],
[X^2,(7,6)],
[X^2,(7,7)],
[X^2,(7,8)],
[X^2,(7,9)],
[X^2,(7,10)],
[X^2,(7,11)],
[X^2,(7,12)])
>>> factor(charpoly(M.T(Integer(2))))
(x - 7) * (x + 7) * (x - 9)^2 * (x + 5)^2
* (x^2 - x - 4)^2 * (x^2 + 9)^2
>>> dimension(M.cuspidal_subspace())
10
M = ModularSymbols(GammaH(13,[3]), weight=4) M M.basis() factor(charpoly(M.T(2))) dimension(M.cuspidal_subspace())
{Cremona’s Modular Symbols Library} Sage includes John Cremona’s specialized and insanely fast implementation of modular symbols for weight 2 and trivial character. We illustrate below computing the space of modular symbols of level 20014, which has dimension \(5005\), along with a Hecke operator on this space. The whole computation below takes only a few seconds; a similar computation takes a few minutes using Sage’s generic modular symbols code. Moreover, Cremona has done computations at levels over 200,000 using his library, so the code is known to scale well to large problems. The new code in Sage for modular symbols is much more general, but doesn’t scale nearly so well (yet).
sage: M = CremonaModularSymbols(20014) # few seconds
sage: M
Cremona Modular Symbols space of dimension 5005 for
Gamma_0(20014) of weight 2 with sign 0
sage: t = M.hecke_matrix(3) # few seconds
>>> from sage.all import *
>>> M = CremonaModularSymbols(Integer(20014)) # few seconds
>>> M
Cremona Modular Symbols space of dimension 5005 for
Gamma_0(20014) of weight 2 with sign 0
>>> t = M.hecke_matrix(Integer(3)) # few seconds
M = CremonaModularSymbols(20014) # few seconds M t = M.hecke_matrix(3) # few seconds
Enumerating Totally Real Number Fields¶
As part of his project to enumerate Shimura curves, John Voight has contributed code to Sage for enumerating totally real number fields. The algorithm isn’t extremely complicated, but it involves some “inner loops” that have to be coded to run very quickly. Using Cython, Voight was able to implement exactly the variant of Newton iteration that he needed for his problem.
The function enumerate_totallyreal_fields_prim(n, B, ...)
enumerates without using a database (!) primitive (no proper subfield)
totally real fields of degree \(n>1\) with discriminant \(d
\leq B\).
We compute the totally real quadratic fields of discriminant \(\leq 50\). The calculation below, which is almost instant, is done in real time and is not a table lookup.
sage: enumerate_totallyreal_fields_prim(2,50)
[[5, x^2 - x - 1], [8, x^2 - 2], [12, x^2 - 3], [13, x^2 - x - 3],
[17, x^2 - x - 4], [21, x^2 - x - 5], [24, x^2 - 6], [28, x^2 - 7],
[29, x^2 - x - 7], [33, x^2 - x - 8], [37, x^2 - x - 9],
[40, x^2 - 10], [41, x^2 - x - 10], [44, x^2 - 11]]
>>> from sage.all import *
>>> enumerate_totallyreal_fields_prim(Integer(2),Integer(50))
[[5, x^2 - x - 1], [8, x^2 - 2], [12, x^2 - 3], [13, x^2 - x - 3],
[17, x^2 - x - 4], [21, x^2 - x - 5], [24, x^2 - 6], [28, x^2 - 7],
[29, x^2 - x - 7], [33, x^2 - x - 8], [37, x^2 - x - 9],
[40, x^2 - 10], [41, x^2 - x - 10], [44, x^2 - 11]]
enumerate_totallyreal_fields_prim(2,50)
We compute all totally real quintic fields of discriminant \(\leq 10^5\). Again, this is done in real time - it’s not a table lookup!
sage: enumerate_totallyreal_fields_prim(5,10^5)
[[14641, x^5 - x^4 - 4*x^3 + 3*x^2 + 3*x - 1],
[24217, x^5 - 5*x^3 - x^2 + 3*x + 1],
[36497, x^5 - 2*x^4 - 3*x^3 + 5*x^2 + x - 1],
[38569, x^5 - 5*x^3 + 4*x - 1],
[65657, x^5 - x^4 - 5*x^3 + 2*x^2 + 5*x + 1],
[70601, x^5 - x^4 - 5*x^3 + 2*x^2 + 3*x - 1],
[81509, x^5 - x^4 - 5*x^3 + 3*x^2 + 5*x - 2],
[81589, x^5 - 6*x^3 + 8*x - 1],
[89417, x^5 - 6*x^3 - x^2 + 8*x + 3]]
>>> from sage.all import *
>>> enumerate_totallyreal_fields_prim(Integer(5),Integer(10)**Integer(5))
[[14641, x^5 - x^4 - 4*x^3 + 3*x^2 + 3*x - 1],
[24217, x^5 - 5*x^3 - x^2 + 3*x + 1],
[36497, x^5 - 2*x^4 - 3*x^3 + 5*x^2 + x - 1],
[38569, x^5 - 5*x^3 + 4*x - 1],
[65657, x^5 - x^4 - 5*x^3 + 2*x^2 + 5*x + 1],
[70601, x^5 - x^4 - 5*x^3 + 2*x^2 + 3*x - 1],
[81509, x^5 - x^4 - 5*x^3 + 3*x^2 + 5*x - 2],
[81589, x^5 - 6*x^3 + 8*x - 1],
[89417, x^5 - 6*x^3 - x^2 + 8*x + 3]]
enumerate_totallyreal_fields_prim(5,10^5)
Bernoulli Numbers¶
Mathematica and Pari¶
From the Mathematica website:
“Today We Broke the Bernoulli Record: From the Analytical Engine to Mathematica April 29, 2008 Oleksandr Pavlyk, Kernel Technology A week ago, I took our latest development version of Mathematica, and I typed
BernoulliB[10^7]. And then I waited. Yesterday–5 days, 23 hours, 51 minutes, and 37 seconds later–I got the result!”
Tom Boothby did that same computation in Sage, which uses Pari’s bernfrac command that uses evaluation of \(\zeta\) and factorial to high precision, and it took 2 days, 12 hours.
David Harvey’s bernmm¶
Then David Harvey came up with an entirely new algorithm that parallelizes well. He gives these timings for computing \(B_{10^7}\) on his machine (it takes 59 minutes, 57 seconds on my 16-core 1.8ghz Opteron box):
PARI: 75 h, Mathematica: 142 h
bernmm (1 core) = 11.1 h, bernmm (10 cores) = 1.3 h
“Running on 10 cores for 5.5 days, I [David Harvey] computed [the Bernoulli number] \(B_k\) for \(k = 10^8\), which I believe is a new record. Essentially it’s the multimodular algorithm I suggested earlier on this thread, but I figured out some tricks to optimise the crap out of the computation of \(B_k \text{mod} p\).”
So now Sage is the fastest in the world for large Bernoulli numbers. The timings below are on a 24-core 2.66Ghz Xeon box.
sage: w1 = bernoulli(100000, num_threads=16) # 0.9 seconds wall time
sage: w2 = bernoulli(100000, algorithm='pari') # long time (6s on sage.math, 2011)
sage: w1 == w2 # long time
True
>>> from sage.all import *
>>> w1 = bernoulli(Integer(100000), num_threads=Integer(16)) # 0.9 seconds wall time
>>> w2 = bernoulli(Integer(100000), algorithm='pari') # long time (6s on sage.math, 2011)
>>> w1 == w2 # long time
True
w1 = bernoulli(100000, num_threads=16) # 0.9 seconds wall time w2 = bernoulli(100000, algorithm='pari') # long time (6s on sage.math, 2011) w1 == w2 # long time
Polynomial Arithmetic¶
FLINT: Univariate Polynomial Arithmetic¶
Sage uses Bill Hart and David Harvey’s GPL’d Flint C library for arithmetic in \(\ZZ[x]\). Its main claim to fame is that it is the world’s fastest for polynomial multiplication, e.g., in the benchmark below it is faster than NTL and Magma on some systems (though such benchmarks of course change as software improves). Behind the scenes Flint contains some carefully tuned discrete Fourier transform code.
sage: Rflint = PolynomialRing(ZZ, 'x')
sage: f = Rflint([ZZ.random_element(2^64) for _ in [1..32]])
sage: g = Rflint([ZZ.random_element(2^64) for _ in [1..32]])
sage: timeit('f*g') # random output
625 loops, best of 3: 105 microseconds per loop
sage: Rntl = PolynomialRing(ZZ, 'x', implementation='NTL')
sage: f = Rntl([ZZ.random_element(2^64) for _ in [1..32]])
sage: g = Rntl([ZZ.random_element(2^64) for _ in [1..32]])
sage: timeit('f*g') # random output
625 loops, best of 3: 310 microseconds per loop
sage: ff = magma(f); gg = magma(g) #optional - magma
sage: s = 'time v := [%s * %s : i in [1..10^5]];'%(ff.name(), gg.name()) #optional - magma
sage: magma.eval(s) #optional - magma
'Time: ...'
>>> from sage.all import *
>>> Rflint = PolynomialRing(ZZ, 'x')
>>> f = Rflint([ZZ.random_element(Integer(2)**Integer(64)) for _ in (ellipsis_range(Integer(1),Ellipsis,Integer(32)))])
>>> g = Rflint([ZZ.random_element(Integer(2)**Integer(64)) for _ in (ellipsis_range(Integer(1),Ellipsis,Integer(32)))])
>>> timeit('f*g') # random output
625 loops, best of 3: 105 microseconds per loop
>>> Rntl = PolynomialRing(ZZ, 'x', implementation='NTL')
>>> f = Rntl([ZZ.random_element(Integer(2)**Integer(64)) for _ in (ellipsis_range(Integer(1),Ellipsis,Integer(32)))])
>>> g = Rntl([ZZ.random_element(Integer(2)**Integer(64)) for _ in (ellipsis_range(Integer(1),Ellipsis,Integer(32)))])
>>> timeit('f*g') # random output
625 loops, best of 3: 310 microseconds per loop
>>> ff = magma(f); gg = magma(g) #optional - magma
>>> s = 'time v := [%s * %s : i in [1..10^5]];'%(ff.name(), gg.name()) #optional - magma
>>> magma.eval(s) #optional - magma
'Time: ...'
Rflint = PolynomialRing(ZZ, 'x')
f = Rflint([ZZ.random_element(2^64) for _ in [1..32]])
g = Rflint([ZZ.random_element(2^64) for _ in [1..32]])
timeit('f*g') # random output
Rntl = PolynomialRing(ZZ, 'x', implementation='NTL')
f = Rntl([ZZ.random_element(2^64) for _ in [1..32]])
g = Rntl([ZZ.random_element(2^64) for _ in [1..32]])
timeit('f*g') # random output
ff = magma(f); gg = magma(g) #optional - magma
s = 'time v := [%s * %s : i in [1..10^5]];'%(ff.name(), gg.name()) #optional - magma
magma.eval(s) #optional - magma
Singular: Multivariate Polynomial Arithmetic¶
Multivariate polynomial arithmetic in many cases uses Singular in library mode (due to Martin Albrecht), which is quite fast. For example, below we do the Fateman benchmark over the finite field of order 32003, and compare the timing with Magma.
sage: P.<x,y,z> = GF(32003)[]
sage: p = (x+y+z+1)^10
sage: q = p+1
sage: timeit('p*q') # random output
125 loops, best of 3: 1.53 ms per loop
sage: p = (x+y+z+1)^20
sage: q = p+1
sage: timeit('p*q') # not tested - timeout if SAGE_DEBUG=yes
5 loops, best of 3: 384 ms per loop
sage: pp = magma(p); qq = magma(q) #optional - magma
sage: s = 'time w := %s*%s;'%(pp.name(),qq.name()) #optional - magma
sage: magma.eval(s) #optional - magma
'Time: ...'
>>> from sage.all import *
>>> P = GF(Integer(32003))['x, y, z']; (x, y, z,) = P._first_ngens(3)
>>> p = (x+y+z+Integer(1))**Integer(10)
>>> q = p+Integer(1)
>>> timeit('p*q') # random output
125 loops, best of 3: 1.53 ms per loop
>>> p = (x+y+z+Integer(1))**Integer(20)
>>> q = p+Integer(1)
>>> timeit('p*q') # not tested - timeout if SAGE_DEBUG=yes
5 loops, best of 3: 384 ms per loop
>>> pp = magma(p); qq = magma(q) #optional - magma
>>> s = 'time w := %s*%s;'%(pp.name(),qq.name()) #optional - magma
>>> magma.eval(s) #optional - magma
'Time: ...'
P.<x,y,z> = GF(32003)[]
p = (x+y+z+1)^10
q = p+1
timeit('p*q') # random output
p = (x+y+z+1)^20
q = p+1
timeit('p*q') # not tested - timeout if SAGE_DEBUG=yes
pp = magma(p); qq = magma(q) #optional - magma
s = 'time w := %s*%s;'%(pp.name(),qq.name()) #optional - magma
magma.eval(s) #optional - magma