birds_other.rst

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 and zn_poly (a C library in Sage for
fast arithmetic in :math:`(\ZZ/n\ZZ)[x]`).

We import the hypellfrob function and call it on a polynomial over
:math:`\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
    [ 0 + O(101)  0 + O(101) 93 + O(101) 62 + O(101)]
    [ 0 + O(101)  0 + O(101) 55 + O(101) 19 + O(101)]
    [ 0 + O(101)  0 + O(101) 65 + O(101) 42 + O(101)]
    [ 0 + O(101)  0 + O(101) 89 + O(101) 29 + O(101)]

We do the same calculation but in :math:`\ZZ/101^4\ZZ`,
which gives enough precision to recognize the exact characteristic
polynomial in :math:`\ZZ[x]` of Frobenius as an element of the
endomorphism ring. This computation is still very fast, taking only a
fraction of a second.

.. link

::

    sage: M = hypellfrob(p, 4, f)   # about 0.25 seconds
    sage: M[0,0]
    91844754 + O(101^4)

The characteristic polynomial of Frobenius is :math:`x^4 + 7x^3 +
167x^2 + 707x + 10201`, which determines the :math:`\zeta` function of
the curve :math:`y^2= f(x)`.

.. link

::

    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))

.. 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 :math:`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 :math:`M` of weight :math:`4` modular
symbols for a certain congruence subgroup :math:`\Gamma_H(13)`
of level :math:`13`. Then we compute a basis for this space,
expressed in terms of Manin symbols. Finally, we compute the Hecke
operator :math:`T_2` acting on :math:`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 
    and over Rational Field
    sage: M.basis()
    ([X^2,(0,1)], [X^2,(0,7)], [X^2,(2,5)], [X^2,(2,8)], [X^2,(2,9)], 
     [X^2,(2,10)], [X^2,(2,11)], [X^2,(2,12)], [X^2,(4,0)], [X^2,(4,3)], 
     [X^2,(4,6)], [X^2,(4,8)], [X^2,(4,12)], [X^2,(7,1)])
    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

{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
:math:`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

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 :math:`n>1` with discriminant :math:`d
\leq B`.

We compute the totally real quadratic fields of discriminant
:math:`\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]]

We compute all totally real quintic fields of discriminant
:math:`\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]]


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 :math:`\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
:math:`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] :math:`B_k` for :math:`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
    :math:`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


Polynomial Arithmetic
=====================

FLINT: Univariate Polynomial Arithmetic
---------------------------------------

Sage uses Bill Hart and David Harvey's GPL'd Flint C library for
arithmetic in :math:`\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: ...'

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)^20
    sage: q = p+1
    sage: timeit('p*q')   # random output
    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: ...'