# HG changeset patch # User William Stein and Craig Citro # Date 1313520357 25200 # Node ID 70d3779139a7742d7d2b520e63b0523d04112d53 # Parent f5872366083fe8704a6f7a4954ff5230ca263316 trac 11375: speed up computation of level one eisenstein series diff --git a/sage/modular/modform/eis_series.py b/sage/modular/modform/eis_series.py --- a/sage/modular/modform/eis_series.py +++ b/sage/modular/modform/eis_series.py @@ -23,7 +23,7 @@ is_FiniteField, ZZ, QQ, Integer, divisors, LCM, is_squarefree) from sage.rings.power_series_ring import PowerSeriesRing -from eis_series_cython import eisenstein_series_poly +from eis_series_cython import eisenstein_series_poly, Ek_ZZ def eisenstein_series_qexp(k, prec = 10, K=QQ, var='q') : r""" @@ -91,16 +91,22 @@ R = PowerSeriesRing(K, var) - if K is QQ : - return a0fac*R(eisenstein_series_poly(k, prec).list(), prec=prec, check=False) + if K == QQ: + ls = Ek_ZZ(k, prec) + # The following is *dramatically* faster than doing the more natural + # "R(ls)" would be: + E = ZZ[var](ls, prec=prec, check=False).change_ring(QQ) + if len(ls)>0: + E._unsafe_mutate(0, a0) + return R(E, prec) + # The following is an older slower alternative to the above three lines: + #return a0fac*R(eisenstein_series_poly(k, prec).list(), prec=prec, check=False) else: # this is a temporary fix due to a change in the # polynomial constructor over finite fields; this # is a notable speed regression, to be fixed soon. return a0fac*R(eisenstein_series_poly(k, prec).list(), prec=prec, check=True) -###################################################################### - def __common_minimal_basering(chi, psi): """ Find the smallest basering over which chi and psi are valued, and diff --git a/sage/modular/modform/eis_series_cython.pyx b/sage/modular/modform/eis_series_cython.pyx --- a/sage/modular/modform/eis_series_cython.pyx +++ b/sage/modular/modform/eis_series_cython.pyx @@ -6,10 +6,136 @@ from sage.rings.rational_field import QQ from sage.rings.power_series_ring import PowerSeriesRing +from sage.rings.integer cimport Integer from sage.rings.arith import primes, bernoulli -from sage.rings.integer cimport Integer +from sage.rings.fast_arith cimport prime_range from sage.libs.flint.fmpz_poly cimport Fmpz_poly +cpdef Ek_ZZ(int k, int prec=10): + """ + Return list of prec integer coefficients of the weight k + Eisenstein series of level 1, normalized so the coefficient of q + is 1, except that the 0th coefficient is set to 1 instead of its + actual value. + + INPUT: + + - `k` -- int + - ``prec`` -- int + + OUTPUT: + + - list of Sage Integers. + + EXAMPLES:: + + sage: from sage.modular.modform.eis_series_cython import Ek_ZZ + sage: Ek_ZZ(4,10) + [1, 1, 9, 28, 73, 126, 252, 344, 585, 757] + sage: [sigma(n,3) for n in [1..9]] + [1, 9, 28, 73, 126, 252, 344, 585, 757] + sage: Ek_ZZ(10,10^3) == [1] + [sigma(n,9) for n in range(1,10^3)] + True + """ + cdef Integer pp + cdef mpz_t q, current_sum, q_plus_one + + cdef unsigned long p + cdef Py_ssize_t i, ind + cdef bint continue_flag + + cdef list power_sum_ls + + cdef unsigned long max_power_sum, temp_index + cdef unsigned long remainder, prev_index + cdef unsigned long additional_p_powers + + mpz_init(q) + mpz_init(current_sum) + mpz_init(q_plus_one) + + # allocate the list for the result + cdef list val = [] + for i from 0 <= i < prec: + pp = (PY_NEW(Integer)) + mpz_set_si(pp.value, 1) + val.append(pp) + + # no need to reallocate this list every time -- just reuse the + # Integers in it + power_sum_ls = [] + max_power_sum = prec + while max_power_sum: + max_power_sum >>= 1 + pp = (PY_NEW(Integer)) + mpz_set_si(pp.value, 1) + power_sum_ls.append(pp) + + for pp in prime_range(prec): + p = mpz_get_ui((pp).value) + mpz_ui_pow_ui(q, p, k - 1) + mpz_add_ui(q_plus_one, q, 1) + mpz_set(current_sum, q_plus_one) + + # NOTE: I (wstein) did benchmarks, and the use of + # PyList_GET_ITEM in the code below is worth it since it leads + # to a significant speedup by not doing bounds checking. + + # set power_sum_ls[1] = q+1 + mpz_set(((PyList_GET_ITEM(power_sum_ls, 1))).value, + current_sum) + max_power_sum = 1 + + ind = 0 + while True: + continue_flag = 0 + # do the first p-1 + for i from 0 < i < p: + ind += p + if (ind >= prec): + continue_flag = 1 + break + mpz_mul(((PyList_GET_ITEM(val, ind))).value, + ((PyList_GET_ITEM(val, ind))).value, + q_plus_one) + ind += p + if (ind >= prec or continue_flag): + break + + # now do the pth one, which is harder. + + # compute the valuation of n at p + additional_p_powers = 0 + temp_index = ind / p + remainder = 0 + while not remainder: + additional_p_powers += 1 + prev_index = temp_index + temp_index = temp_index / p + remainder = prev_index - p*temp_index + + # if we need a new sum, it has to be the next uncomputed one. + if (additional_p_powers > max_power_sum): + mpz_mul(current_sum, current_sum, q) + mpz_add_ui(current_sum, current_sum, 1) + max_power_sum += 1 + + mpz_set(((PyList_GET_ITEM(power_sum_ls, + max_power_sum))).value, + current_sum) + + # finally, set the coefficient + mpz_mul(((PyList_GET_ITEM(val, ind))).value, + ((PyList_GET_ITEM(val, ind))).value, + ((PyList_GET_ITEM(power_sum_ls, + additional_p_powers))).value) + + mpz_clear(q) + mpz_clear(current_sum) + mpz_clear(q_plus_one) + + return val + cpdef eisenstein_series_poly(int k, int prec = 10) : r""" @@ -21,9 +147,9 @@ :func:`~sage.modular.modform.eis_series.eisenstein_series_qexp` and :func:`~sage.modular.modform.vm_basis.victor_miller_basis`; see the docstring of the former for further details. - - EXAMPLES : - + + EXAMPLES:: + sage: from sage.modular.modform.eis_series_cython import eisenstein_series_poly sage: eisenstein_series_poly(12, prec=5) 5 691 65520 134250480 11606736960 274945048560 @@ -91,7 +217,7 @@ fmpz_poly_set_coeff_mpz(res.poly, prec-1, val[prec-1]) for i from 1 <= i < prec - 1 : fmpz_poly_set_coeff_mpz(res.poly, i, val[i]) - + fmpz_poly_scalar_mul_mpz(res.poly, res.poly, ((a0.denominator())).value) fmpz_poly_set_coeff_mpz(res.poly, 0, ((a0.numerator())).value) diff --git a/sage/rings/fast_arith.pxd b/sage/rings/fast_arith.pxd --- a/sage/rings/fast_arith.pxd +++ b/sage/rings/fast_arith.pxd @@ -1,5 +1,7 @@ include "../ext/cdefs.pxi" +cpdef prime_range(start, stop=*, algorithm=*, bint py_ints=*) + cdef class arith_int: cdef int abs_int(self, int x) except -1 cdef int sign_int(self, int n) except -2