Cythonize hypergeometric trace formula

kedlaya commented 4 years ago

The goal of this ticket is to improve upon this comparison between Magma:

> H := HypergeometricData([3],[4]);
> time EulerFactor(H, 1/2, 1000003);                                           
1000003*$.1^2 - 788*$.1 + 1
Time: 2.840
> H := HypergeometricData([10],[5]);  
> time EulerFactor(H, 1/2, 1009);
1018081*$.1^4 + 8072*$.1^3 + 590*$.1^2 + 8*$.1 + 1
Time: 4.890
> H := HypergeometricData([[10,6],[5,4]]);
> time EulerFactor(H, 1/2, 61);
11694146092834141*$.1^6 - 14271143697997*$.1^5 - 13901224364*$.1^4 + 
    104106138*$.1^3 - 61244*$.1^2 - 277*$.1 + 1
Time: 4.670

and Sage:

sage: from sage.modular.hypergeometric_motive import HypergeometricData as HG
sage: H = HG(cyclotomic=[[3],[4]])
sage: time H.euler_factor(2, 1000003)
Wall time: 1min 5s
1000003*T^2 - 788*T + 1
sage: H = HG(cyclotomic=[[10],[5]])
sage: time H.euler_factor(2, 1009)
Wall time: 49 s
1018081*T^4 + 8072*T^3 + 590*T^2 + 8*T + 1
sage: H = HG(cyclotomic=[[10,6],[5,4]])
sage: time H.euler_factor(2, 61)
Wall time: 14.8 s
11694146092834141*T^6 - 14271143697997*T^5 - 13901224364*T^4 + 104106138*T^3 - 61244*T^2 - 277*T + 1

by moving two bottleneck steps into Cython.

Computing Gauss sums via the Gross-Koblitz formula, in the gauss_table function in rings/padics/padic_generic_element.pyx.
- We refactor Sage's implementation of the p-adic Gamma function to expose some internals as C functions.
- When possible, we do the p-adic computations using C longs.
- We short-circuit the computation when the working precision equals 1. This benefits the case of a degree-2 L-function.
- We enable global caching of the Gauss sums, as in Magma's HypergeometricMotiveSaveLimit function.
Computing products of Gauss sums needed for the hypergeometric trace formula, in the hgm_coeffs function in modular/hypergeometric_misc.pyx.
- We again use C longs when possible.

CC: @roed314

Component: modular forms

Keywords: hypergeometric motives

Author: Kiran Kedlaya

Branch/Commit: a7bfd25

Reviewer: David Roe, Frédéric Chapoton

Issue created by migration from https://trac.sagemath.org/ticket/28902

kedlaya commented 4 years ago

Branch: u/kedlaya/cythonize_hgm

kedlaya commented 4 years ago

comment:2

One additional change: I refined the precision bound a bit. Particularly in degree 2, it makes a substantial difference to work with one p-adic digit rather than two (the computation of p-adic Gamma scales noticeably linearly in the number of digits).

Timings I recorded on this branch (same machine as before):

sage: from sage.modular.hypergeometric_motive import HypergeometricData as HG
sage: H = HG(cyclotomic=[[3],[4]])
sage: time H.euler_factor(2, 1000003)
Wall time: 13 s
1000003*T^2 - 788*T + 1
sage: H = HG(cyclotomic=[[10],[5]])
sage: time H.euler_factor(2, 1009)
Wall time: 14.6 s
1018081*T^4 + 8072*T^3 + 590*T^2 + 8*T + 1
sage: H = HG(cyclotomic=[[10,6],[5,4]])
sage: time H.euler_factor(2, 61)
Wall time: 4.87 s
11694146092834141*T^6 - 14271143697997*T^5 - 13901224364*T^4 + 104106138*T^3 - 61244*T^2 - 277*T + 1

New commits:

d23d2e8 Cythonize hypergeometric trace formula

kedlaya commented 4 years ago

Commit: d23d2e8

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

0948ef0 Move evaluation of Dwork expansion to a cdef function

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Changed commit from d23d2e8 to 0948ef0

kedlaya commented 4 years ago

comment:4

In this commit, I moved the evaluation of the Mahler series for p-adic Gamma into a helper cdef function; this eliminates some Python call overhead. As a bonus, I can undo the change in call syntax to dwork_expansion and thereby restore backward compatibility.

I also made some additional optimizations in the case where the working precision is 1, as in the first of my examples.

Finally, I moved the precomputation of the Mahler series into padics/misc.py. This seems like a step backward since that is not a Cython file, but conceptually it seemed to make sense, and it doesn't have much effect on the timings.

With this commit, the timings are down to 9.15s, 11s, 3.81s. Some cleanup and doctests still needed.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

f3064b8 Add early abort in main loop of hgm_coeffs

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Changed commit from 0948ef0 to f3064b8

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

6051397 Add doctests

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Changed commit from f3064b8 to 6051397

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Changed commit from 6051397 to 04d14f2

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

04d14f2 Add cdef local variable to dwork_mahler_coeffs

kedlaya commented 4 years ago

comment:8

This batch of commits implements a speedup to hgm_coeffs: it turns out that some of the coefficients are forced to be zero for combinatorial reasons, and it is helpful to order the computation so that this vanishing can be detected for an early abort. The relevance of this will vary with the hypergeometric data.

For the examples I am testing, the latest timings are 5.48s, 7.81s, 2.77s.

I also added doctests to restore 100% coverage in the affected files.

New commits:

04d14f2 Add cdef local variable to dwork_mahler_coeffs

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Changed commit from 04d14f2 to 7c1e45a

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

7c1e45a In HGM, get p_ring from an always defined entity

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

3160926 Optimizations in inner loop of gauss_table

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Changed commit from 7c1e45a to 3160926

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Changed commit from 3160926 to ebad8fb

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

`a1519c3`	`Refactor: dwork_mahler_coeffs to padic_generic_element.pyx`
`ebad8fb`	`Use longs when possible in gauss_table`

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Changed commit from ebad8fb to e5138ac

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

`48dfb02`	`Make _gauss_table a class attribute to simplify cache management`
`56e1b7e`	`Pass integer part of gauss_table as an array of longs`
`e5138ac`	`Use longs when possible in hgm_coeffs`

kedlaya commented 4 years ago

Changed branch from u/kedlaya/cythonize_hgm to u/kedlaya/cythonize_hgm_v2

kedlaya commented 4 years ago

Changed commit from e5138ac to none

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Commit: 9fc8f1b

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

9fc8f1b Summary commit of Cythonization of HGMs

kedlaya commented 4 years ago

comment:15

To clean things up, I squashed the preceding commits and the next one into a single commit on a new branch. Timings in the three examples are now quite comparable with Magma: 3.57s, 4.75s, 2.73s.

To test things more deeply (particularly the cacheing functionality), I reran a computation from the LMFDB: this iterates over all hypergeometric data of degree 2 through 10 and computes various Euler factors. All answers agree with Magma.

At this point I think I'm ready to put this to review.

kedlaya commented 4 years ago

Author: Kiran Kedlaya

kedlaya commented 4 years ago

Description changed:

--- 
+++ 
@@ -32,9 +32,13 @@
 Wall time: 14.8 s
 11694146092834141*T^6 - 14271143697997*T^5 - 13901224364*T^4 + 104106138*T^3 - 61244*T^2 - 277*T + 1

-by moving some bottleneck steps into Cython. +by moving two bottleneck steps into Cython.

-- dwork_expansion, used to compute the p-adic Gamma function. -- Computing Gauss sums via the Gross-Koblitz formula. These should be cached if one is going to compute the trace formula multiple times for the same prime (Magma already provides this option). -- Computing the coefficients of the hypergeometric trace formula. These should be cached if one is going to fix the hypergeometric data and the prime, and loop over parameter values. +- Computing Gauss sums via the Gross-Koblitz formula, in the gauss_table function in rings/padics/padic_generic_element.pyx.

- We refactor Sage's implementation of the p-adic Gamma function to expose some internals as C functions.
- When possible, we do the p-adic computations using C longs.
- We short-circuit the computation when the working precision equals 1. This benefits the case of a degree-2 L-function.
- We enable global caching of the Gauss sums, as in Magma's HypergeometricMotiveSaveLimit function. +- Computing products of Gauss sums needed for the hypergeometric trace formula, in the hgm_coeffs function in modular/hypergeometric_misc.pyx.
- We again use C longs when possible.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Changed commit from 9fc8f1b to 98fc6b7

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

`848c515`	`Summary commit of Cythonization of HGMs`
`65d3cac`	`Fix pyflakes issues in hypergeometric_motive.py`
`98fc6b7`	`Merge remote-tracking branch 'refs/remotes/trac/u/kedlaya/cythonize_hgm_v2' into cythonize_hgm_v2`

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Changed commit from 98fc6b7 to 3b3873d

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

3b3873d Remove non-ASCII character from docstring

kedlaya commented 4 years ago

comment:18

Some of the patchbots found a non-ASCII issue. This should fix that.

roed314 commented 4 years ago

Changed branch from u/kedlaya/cythonize_hgm_v2 to u/roed/cythonize_hgm_v2

roed314 commented 4 years ago

Reviewer: David Roe

roed314 commented 4 years ago

comment:20

I made a few minor changes and one somewhat substantial one (switching from a Python list to a C array of long longs, which gave a 33% speedup in the second example when I tested). If you're happy with the changes and tests pass, I'm happy to give this ticket a positive review.

New commits:

d70c842 Minor tweaks and change to long long* for dwork coeffs when using longs

roed314 commented 4 years ago

Changed commit from 3b3873d to d70c842

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

0cb6788 Remove unnecessary import

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Changed commit from d70c842 to 0cb6788

fchapoton commented 4 years ago

comment:22

looks good to me

side note (no action needed) : pep8 does not require ** to be framed with spaces, but does require that for / and *

fchapoton commented 4 years ago

Changed reviewer from David Roe to David Roe, Frédéric Chapoton

kedlaya commented 4 years ago

Changed branch from u/roed/cythonize_hgm_v2 to u/kedlaya/cythonize_hgm_v2

kedlaya commented 4 years ago

Changed commit from 0cb6788 to c1aa950

kedlaya commented 4 years ago

comment:25

If we're going to use an array of long longs, I think it makes more sense to use a Python array (which is almost as efficient) so that we can send it directly from one subroutine to the other without conversion to a Python list.

Besides making that change, I also found an overflow bug that I had introduced in this branch and added a doctest for it.

Timings for my three examples are now 2.0s, 3.4s, 2.9s. In particular we are now beating Magma in all of these cases (though they are far from exhaustive).

I also enabled the use of the long-based code when cacheing, which wasn't happening before. Not that this makes any appreciable difference when cacheing, but I wanted my LMFDB test run to stress-test this modification.

Patchbots seem satisfied with the earlier versions, but maybe it's worth waiting to see some runs from this latest commit.

New commits:

`a14aedb`	`Summary commit of Cythonization of HGMs`
`4090795`	`Fix pyflakes issues in hypergeometric_motive.py`
`a8f085e`	`Remove non-ASCII character from docstring`
`e0de275`	`Merge branch 'u/kedlaya/cythonize_hgm_v2' of git://trac.sagemath.org/sage into cythonize_hgm_v2`
`90c7d44`	`Merge branch 'u/roed/cythonize_hgm_v2' of git://trac.sagemath.org/sage into cythonize_hgm_v2`
`62f9b7c`	`Pass arrays of long longs back and forth when possible`
`bbaa0f7`	`Fix overflow in Gauss sums; use longs when cacheing`
`c1aa950`	`Pass table precision manually to hgm_coeffs`

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

de82d2e Fix doctest in hypergeometric_misc

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Changed commit from c1aa950 to de82d2e

kedlaya commented 4 years ago

comment:28

Fixed a doctest failure caught by patchbot.

fchapoton commented 4 years ago

comment:29

green bot, I am setting back to positive

vbraun commented 4 years ago

comment:30

I'm getting this on the python2 buildbot

**********************************************************************
File "src/sage/modular/hypergeometric_misc.pyx", line 19, in sage.modular.hypergeometric_misc.?
Failed example:
    H.euler_factor(2, 7, cache_p=True)
Exception raised:
    Traceback (most recent call last):
      File "/var/lib/buildbot/slave/sage2_git/build/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 681, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/var/lib/buildbot/slave/sage2_git/build/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 1123, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.modular.hypergeometric_misc.?[4]>", line 1, in <module>
        H.euler_factor(Integer(2), Integer(7), cache_p=True)
      File "sage/misc/cachefunc.pyx", line 1948, in sage.misc.cachefunc.CachedMethodCaller.__call__ (build/cythonized/sage/misc/cachefunc.c:10273)
        w = self._instance_call(*args, **kwds)
      File "sage/misc/cachefunc.pyx", line 1824, in sage.misc.cachefunc.CachedMethodCaller._instance_call (build/cythonized/sage/misc/cachefunc.c:9758)
        return self.f(self._instance, *args, **kwds)
      File "/var/lib/buildbot/slave/sage2_git/build/local/lib/python2.7/site-packages/sage/modular/hypergeometric_motive.py", line 1498, in euler_factor
        for i in range(bound)]
      File "sage/misc/cachefunc.pyx", line 1948, in sage.misc.cachefunc.CachedMethodCaller.__call__ (build/cythonized/sage/misc/cachefunc.c:10273)
        w = self._instance_call(*args, **kwds)
      File "sage/misc/cachefunc.pyx", line 1824, in sage.misc.cachefunc.CachedMethodCaller._instance_call (build/cythonized/sage/misc/cachefunc.c:9758)
        return self.f(self._instance, *args, **kwds)
      File "/var/lib/buildbot/slave/sage2_git/build/local/lib/python2.7/site-packages/sage/modular/hypergeometric_motive.py", line 1234, in padic_H_value
        trcoeffs = hgm_coeffs(p, f, prec, gamma, m, D, gtab, use_longs)
      File "sage/modular/hypergeometric_misc.pyx", line 6, in sage.modular.hypergeometric_misc.hgm_coeffs (build/cythonized/sage/modular/hypergeometric_misc.c:2809)
        cpdef hgm_coeffs(long long p, int f, int prec, gamma, array.array m, int D,
      File "sage/modular/hypergeometric_misc.pyx", line 37, in sage.modular.hypergeometric_misc.hgm_coeffs (build/cythonized/sage/modular/hypergeometric_misc.c:1674)
        cdef array.array digit_count = array.array('Q', [0]) * q1
    ValueError: bad typecode (must be c, b, B, u, h, H, i, I, l, L, f or d)

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 4 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

a7bfd25 Use Py2-compatible arrays

sagemath / sage

Cythonize hypergeometric trace formula #28902