Ticket #3542: 3542-dharvey-bernoulli-numbers-multimodular.patch

sage/rings/arith.py

@@ -145 +145 @@
 
 algebraic_dependency = algdep
 
-def bernoulli(n, algorithm='pari'):
+def bernoulli(n, algorithm='default', num_threads=1):
     r"""
     Return the n-th Bernoulli number, as a rational number.
 
     INPUT:
         n -- an integer
         algorithm:
-            'pari' -- (default) use the PARI C library, which is
-                      by *far* the fastest.
+            'default' -- (default) use 'pari' for n <= 30000, and
+                         'bernmm' for n > 30000 (this is just a heuristic,
+                         and not guaranteed to be optimal on all hardware)
+            'pari' -- use the PARI C library
             'gap'  -- use GAP
             'gp'   -- use PARI/GP interpreter
             'magma' -- use MAGMA (optional)
             'python' -- use pure Python implementation
+            'bernmm' -- use bernmm package (a multimodular algorithm)
+        num_threads -- positive integer, number of threads to use
+                       (only used for bernmm algorithm)
 
     EXAMPLES:
         sage: bernoulli(12)
@@ -176 +181 @@
         -691/2730
         sage: bernoulli(12, algorithm='python')
         -691/2730
+        sage: bernoulli(12, algorithm='bernmm')
+        -691/2730
+        sage: bernoulli(12, algorithm='bernmm', num_threads=4)
+        -691/2730
 
     \note{If $n>50000$ then algorithm = 'gp' is used instead of
     algorithm = 'pari', since the C-library interface to PARI
@@ -185 +194 @@
     """
     from sage.rings.all import Integer, Rational
     n = Integer(n)
+
+    if algorithm == 'default':
+        algorithm = 'pari' if n <= 30000 else 'bernmm'
+
     if n > 50000 and algorithm == 'pari':
         algorithm = 'gp'
+
     if algorithm == 'pari':
         x = pari(n).bernfrac()         # Use the PARI C library
         return Rational(x)
@@ -205 +219 @@
     elif algorithm == 'python':
         import sage.rings.bernoulli
         return sage.rings.bernoulli.bernoulli_python(n)
+    elif algorithm == 'bernmm':
+        import sage.rings.bernmm
+        return sage.rings.bernmm.bernmm_bern_rat(n, num_threads)
     else:
         raise ValueError, "invalid choice of algorithm"
 

(a) /dev/null vs. (b) b/sage/rings/bernmm.pyx

@@ +1 @@
+r"""
+Cython wrapper for bernmm library
+
+AUTHOR: 
+    - David Harvey (2008-06): initial version
+"""
+
+#*****************************************************************************
+#       Copyright (C) 2008 William Stein <wstein@gmail.com>
+#                     2008 David Harvey <dmharvey@math.harvard.edu>
+#
+#  Distributed under the terms of the GNU General Public License (GPL)
+#                  http://www.gnu.org/licenses/
+#*****************************************************************************
+
+include "../ext/cdefs.pxi"
+include "../ext/interrupt.pxi"
+
+
+cdef extern from "bernmm/bern_rat.h":
+    void bern_rat "bernmm::bern_rat" (mpq_t res, long k, int num_threads)
+
+cdef extern from "bernmm/bern_modp.h":
+    long bern_modp "bernmm::bern_modp" (long p, long k)
+
+
+
+from sage.rings.rational cimport Rational
+
+
+def bernmm_bern_rat(long k, int num_threads = 1):
+    r"""
+    Computes k-th Bernoulli number using a multimodular algorithm.
+    (Wrapper for bernmm library.)
+
+    INPUT:
+        k -- non-negative integer
+        num_threads -- integer >= 1, number of threads to use
+
+    COMPLEXITY:
+        Pretty much quadratic in $k$. See the paper ``A multimodular algorithm
+        for computing Bernoulli numbers'', David Harvey, 2008, for more details.
+
+    EXAMPLES:
+        sage: from sage.rings.bernmm import bernmm_bern_rat
+
+        sage: bernmm_bern_rat(0)
+        1
+        sage: bernmm_bern_rat(1)
+        -1/2
+        sage: bernmm_bern_rat(2)
+        1/6
+        sage: bernmm_bern_rat(3)
+        0
+        sage: bernmm_bern_rat(100)
+        -94598037819122125295227433069493721872702841533066936133385696204311395415197247711/33330
+        sage: bernmm_bern_rat(100, 3)
+        -94598037819122125295227433069493721872702841533066936133385696204311395415197247711/33330
+        
+    TESTS:
+        sage: lst1 = [ bernoulli(2*k, algorithm='bernmm', num_threads=2) for k in [2932, 2957, 3443, 3962, 3973] ]
+        sage: lst2 = [ bernoulli(2*k, algorithm='pari') for k in [2932, 2957, 3443, 3962, 3973] ]
+        sage: lst1 == lst2
+        True
+        sage: [ Zmod(101)(t) for t in lst1 ]
+        [77, 72, 89, 98, 86]
+        sage: [ Zmod(101)(t) for t in lst2 ]
+        [77, 72, 89, 98, 86]
+    """
+    cdef Rational x
+
+    if k < 0:
+        raise ValueError, "k must be non-negative"
+
+    x = Rational()
+    _sig_on
+    bern_rat(x.value, k, num_threads)
+    _sig_off
+
+    return x
+
+
+def bernmm_bern_modp(long p, long k):
+    r"""
+    Computes $B_k \mod p$, where $B_k$ is the k-th Bernoulli number.
+
+    If $B_k$ is not $p$-integral, returns -1.
+
+    INPUT:
+        p -- a prime
+        k -- non-negative integer
+
+    COMPLEXITY:
+        Pretty much linear in $p$.
+
+    EXAMPLES:
+        sage: from sage.rings.bernmm import bernmm_bern_modp
+
+        sage: bernoulli(0) % 5, bernmm_bern_modp(5, 0)
+        (1, 1)
+        sage: bernoulli(1) % 5, bernmm_bern_modp(5, 1)
+        (2, 2)
+        sage: bernoulli(2) % 5, bernmm_bern_modp(5, 2)
+        (1, 1)
+        sage: bernoulli(3) % 5, bernmm_bern_modp(5, 3)
+        (0, 0)
+        sage: bernoulli(4), bernmm_bern_modp(5, 4)
+        (-1/30, -1)
+        sage: bernoulli(18) % 5, bernmm_bern_modp(5, 18)
+        (4, 4)
+        sage: bernoulli(19) % 5, bernmm_bern_modp(5, 19)
+        (0, 0)
+
+        sage: p = 10000019; k = 1000
+        sage: bernoulli(k) % p
+        1972762
+        sage: bernmm_bern_modp(p, k)
+        1972762
+
+    """
+    cdef long x
+
+    if k < 0:
+        raise ValueError, "k must be non-negative"
+
+    _sig_on
+    x = bern_modp(p, k)
+    _sig_off
+
+    return x
+    
+
+# ============ end of file

(a) /dev/null vs. (b) b/sage/rings/bernmm/README.txt

@@ +1 @@
+This directory contains the bernmm library, by David Harvey.
+
+bernmm is maintained as a separate project, see my website dharvey.net.
+
+If you make any changes/improvements, please tell me about it!
+I might want to merge those changes into the upstream version.
+
+2008/06/19
+-------------------------------------

(a) /dev/null vs. (b) b/sage/rings/bernmm/bern_modp.cpp

@@ +1 @@
+/*
+   bern_modp.cpp:  computing isolated Bernoulli numbers modulo p
+   
+   Copyright (C) 2008, David Harvey
+   
+   This file is part of the bernmm package (version 1.0).
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <limits.h>
+#include <signal.h>
+#include <gmp.h>
+#include <NTL/ZZ.h>
+#include "bern_modp_util.h"
+#include "bern_modp.h"
+
+
+NTL_CLIENT;
+
+
+using namespace std;
+
+
+namespace bernmm {
+
+
+/******************************************************************************
+
+   Computing the main sum (general case)
+
+******************************************************************************/
+
+/*
+   Returns (1 - g^k) B_k / 2k mod p.
+   
+   PRECONDITIONS:
+      5 <= p < NTL_SP_BOUND, p prime
+      2 <= k <= p-3, k even
+      pinv = 1 / ((double) p)
+      g = a multiplicative generator of GF(p), in [0, p)
+*/
+long bernsum_powg(long p, double pinv, long k, long g)
+{
+   long half_gm1 = (g + ((g & 1) ? 0 : p) - 1) / 2;    // (g-1)/2 mod p
+   long g_to_jm1 = 1;
+   long g_to_km1 = PowerMod(g, k-1, p, pinv);
+   long g_to_km1_to_j = g_to_km1;
+   long sum = 0;
+   double g_pinv = ((double) g) / ((double) p);
+   mulmod_precon_t g_to_km1_pinv = PrepMulModPrecon(g_to_km1, p, pinv);
+   
+   for (long j = 1; j <= (p-1)/2; j++)
+   {
+      // at this point,
+      //    g_to_jm1 holds g^(j-1) mod p
+      //    g_to_km1_to_j holds (g^(k-1))^j mod p
+   
+      // update g_to_jm1 and compute q = (g*(g^(j-1) mod p) - (g^j mod p)) / p
+      long q;
+      g_to_jm1 = MulDivRem(q, g_to_jm1, g, p, g_pinv);
+      
+      // compute h = -h_g(g^j) = q - (g-1)/2
+      long h = SubMod(q, half_gm1, p);
+      
+      // add h_g(g^j) * (g^(k-1))^j to running total
+      sum = SubMod(sum, MulMod(h, g_to_km1_to_j, p, pinv), p);
+      
+      // update g_to_km1_to_j
+      g_to_km1_to_j = MulModPrecon(g_to_km1_to_j, g_to_km1, p, g_to_km1_pinv);
+   }
+
+   return sum;
+}
+
+
+
+/******************************************************************************
+
+   Computing the main sum (c = 1/2 case)
+
+******************************************************************************/
+
+
+/*
+   The Expander class stores precomputed information for a fixed integer p,
+   that subsequently permits fast computation of the binary expansion of s/p
+   for any 0 < s < p.
+   
+   The constructor takes p and max_words as input. Must have 1 <= max_words <=
+   MAX_INV. It computes an approximation to 1/p.
+   
+   The function expand(word_t* res, long s, long n) computes n words of s/p.
+   Must have 0 < s < p and 1 <= n <= max_words. The output is written to res.
+   The first word of output is junk. The next n words are the digits of s/p,
+   from least to most significant. The buffer must be at least n+2 words long
+   (even though the first and last words are never used for output).
+
+   A "word" is a word_t, and contains WORD_BITS bits. On most systems this
+   will be an mp_limb_t.
+*/
+
+#define MAX_INV 256
+
+#if (GMP_NAIL_BITS == 0) && (GMP_LIMB_BITS >= LONG_BIT)
+// fast mpn-based version
+
+typedef mp_limb_t word_t;
+#define WORD_BITS GMP_LIMB_BITS
+
+class Expander
+{
+private:
+   // Approximation to 1/p. We store (max_words + 1) limbs.
+   mp_limb_t pinv[MAX_INV + 2];
+   mp_limb_t p;
+   int max_words;
+
+public:
+   Expander(long p, int max_words)
+   {
+      assert(max_words >= 1);
+      assert(max_words <= MAX_INV);
+      
+      this->max_words = max_words;
+      this->p = p;
+      mp_limb_t one = 1;
+      mpn_divrem_1(pinv, max_words + 1, &one, 1, p);
+   }
+
+   void expand(word_t* res, long s, int n)
+   {
+      assert(s > 0 && s < p);
+      assert(n >= 1);
+      assert(n <= max_words);
+   
+      if (s == 1)
+      {
+         // already have 1/p; just copy it
+         for (int i = 1; i <= n; i++)
+            res[i] = pinv[max_words - n + i];
+      }
+      else
+      {
+         mpn_mul_1(res, pinv + max_words - n, n + 1, (mp_limb_t) s);
+
+         // If the first output limb is really close to 0xFFFF..., then there's
+         // a possibility of overflow, so fall back on doing division directly.
+         // This should happen extremely rarely --- essentially never on a
+         // 64-bit system, and very occasionally on a 32-bit system.
+         if (res[0] > -((mp_limb_t) s))
+         {
+            mp_limb_t ss = s;
+            mpn_divrem_1(res, n + 1, &ss, 1, p);
+         }
+      }
+   }
+};
+
+
+#else
+// slow mpz-based version, since GMP is using nails, or mp_limb_t is
+// absurdly narrow
+
+typedef unsigned long word_t;
+#define WORD_BITS LONG_BIT
+
+class Expander
+{
+private:
+   mp_limb_t p;
+   mpz_t temp;
+
+public:
+   Expander(long p, int max_words)
+   {
+      this->p = p;
+      mpz_init(temp);
+   }
+   
+   ~Expander()
+   {
+      mpz_clear(temp);
+   }
+
+   void expand(word_t* res, long s, int n)
+   {
+      assert(s > 0 && s < p);
+      assert(n >= 1);
+      
+      mpz_set_ui(temp, s);
+      mpz_mul_2exp(temp, temp, WORD_BITS * n);
+      mpz_fdiv_q_ui(temp, temp, p);
+      mpz_export(res + 1, NULL, -1, sizeof(word_t), 0, 0, temp);
+   }
+};
+
+#endif
+
+
+
+/*
+   Returns (2^(-k) - 1) 2 B_k / k  mod p.
+   
+   (Note: this is useless if 2^k = 1 mod p.)
+      
+   PRECONDITIONS:
+      5 <= p < NTL_SP_BOUND, p prime
+      2 <= k <= p-3, k even
+      pinv = 1 / ((double) p)
+      g = a multiplicative generator of GF(p), in [0, p)
+      n = multiplicative order of 2 in GF(p)
+*/
+
+#define TABLE_LG_SIZE 8
+#define TABLE_SIZE (((word_t) 1) << TABLE_LG_SIZE)
+#define TABLE_MASK (TABLE_SIZE - 1)
+#define NUM_TABLES (WORD_BITS / TABLE_LG_SIZE)
+
+#if WORD_BITS % TABLE_LG_SIZE != 0
+#error Number of bits in a long must be divisible by TABLE_LG_SIZE
+#endif
+
+long bernsum_pow2(long p, double pinv, long k, long g, long n)
+{
+   // In the main summation loop we accumulate data into the _tables_ array;
+   // tables[y][z] contributes to the final answer with a weight of
+   // 
+   // sum(-(-1)^z[t] * (2^(k-1))^(WORD_BITS - 1 - y * TABLE_LG_SIZE - t) :
+   //                                                  0 <= t < TABLE_LG_SIZE),
+   //
+   // where z[t] denotes the t-th binary digit of z (LSB is t = 0).
+   // The memory footprint for _tables_ is 4KB on a 32-bit machine, or 16KB
+   // on a 64-bit machine, so should fit easily into L1 cache.
+   long tables[NUM_TABLES][TABLE_SIZE];
+   memset(tables, 0, sizeof(long) * NUM_TABLES * TABLE_SIZE);
+   
+   long m = (p-1) / n;
+
+   // take advantage of symmetry (n' and m' from the paper)
+   if (n & 1)
+      m >>= 1;
+   else
+      n >>= 1;
+
+   // g^(k-1)
+   long g_to_km1 = PowerMod(g, k-1, p, pinv);
+   // 2^(k-1)
+   long two_to_km1 = PowerMod(2, k-1, p, pinv);
+   // B^(k-1), where B = 2^WORD_BITS
+   long B_to_km1 = PowerMod(two_to_km1, WORD_BITS, p, pinv);
+   // B^(MAX_INV)
+   long s_jump = PowerMod(2, MAX_INV * WORD_BITS, p, pinv);
+
+   // help speed up modmuls
+   mulmod_precon_t g_pinv = PrepMulModPrecon(g, p, pinv);
+   mulmod_precon_t g_to_km1_pinv = PrepMulModPrecon(g_to_km1, p, pinv);
+   mulmod_precon_t two_to_km1_pinv = PrepMulModPrecon(two_to_km1, p, pinv);
+   mulmod_precon_t B_to_km1_pinv = PrepMulModPrecon(B_to_km1, p, pinv);
+   mulmod_precon_t s_jump_pinv = PrepMulModPrecon(s_jump, p, pinv);
+
+   long g_to_km1_to_i = 1;
+   long g_to_i = 1;
+   long sum = 0;
+
+   // Precompute some of the binary expansion of 1/p; at most MAX_INV words,
+   // or possibly less if n is sufficiently small
+   Expander expander(p, (n >= MAX_INV * WORD_BITS)
+                                       ? MAX_INV : ((n - 1) / WORD_BITS + 1));
+
+   // =========== phase 1: main summation loop
+   
+   // loop over outer sum
+   for (long i = 0; i < m; i++)
+   {
+      // s keeps track of g^i*2^j mod p
+      long s = g_to_i;
+      // x keeps track of (g^i*2^j)^(k-1) mod p
+      long x = g_to_km1_to_i;
+      
+      // loop over inner sum; break it up into chunks of length at most
+      // MAX_INV * WORD_BITS. If n is large, this allows us to do most of
+      // the work with mpn_mul_1 instead of mpn_divrem_1, and also improves
+      // memory locality.
+      for (long nn = n; nn > 0; nn -= MAX_INV * WORD_BITS)
+      {
+         word_t s_over_p[MAX_INV + 2];
+         long bits, words;
+
+         if (nn >= MAX_INV * WORD_BITS)
+         {
+            // do one chunk of length exactly MAX_INV * WORD_BITS
+            bits = MAX_INV * WORD_BITS;
+            words = MAX_INV;
+         }
+         else
+         {
+            // last chunk of length less than MAX_INV * WORD_BITS
+            bits = nn;
+            words = (nn - 1) / WORD_BITS + 1;
+         }
+         
+         // compute some bits of the binary expansion of s/p
+         expander.expand(s_over_p, s, words);
+         word_t* next = s_over_p + words;
+
+         // loop over whole words
+         for (; bits >= WORD_BITS; bits -= WORD_BITS, next--)
+         {
+            word_t y = *next;
+
+#if NUM_TABLES != 8 && NUM_TABLES != 4
+            // generic version
+            for (long h = 0; h < NUM_TABLES; h++)
+            {
+               long& target = tables[h][y & TABLE_MASK];
+               target = SubMod(target, x, p);
+               y >>= TABLE_LG_SIZE;
+            }
+#else
+            // unrolled versions for 32-bit/64-bit machines
+            long& target0 = tables[0][y & TABLE_MASK];
+            target0 = SubMod(target0, x, p);
+
+            long& target1 = tables[1][(y >> TABLE_LG_SIZE) & TABLE_MASK];
+            target1 = SubMod(target1, x, p);
+
+            long& target2 = tables[2][(y >> (2*TABLE_LG_SIZE)) & TABLE_MASK];
+            target2 = SubMod(target2, x, p);
+
+            long& target3 = tables[3][(y >> (3*TABLE_LG_SIZE)) & TABLE_MASK];
+            target3 = SubMod(target3, x, p);
+#if NUM_TABLES == 8
+            long& target4 = tables[4][(y >> (4*TABLE_LG_SIZE)) & TABLE_MASK];
+            target4 = SubMod(target4, x, p);
+
+            long& target5 = tables[5][(y >> (5*TABLE_LG_SIZE)) & TABLE_MASK];
+            target5 = SubMod(target5, x, p);
+
+            long& target6 = tables[6][(y >> (6*TABLE_LG_SIZE)) & TABLE_MASK];
+            target6 = SubMod(target6, x, p);
+
+            long& target7 = tables[7][(y >> (7*TABLE_LG_SIZE)) & TABLE_MASK];
+            target7 = SubMod(target7, x, p);
+#endif
+#endif
+
+            x = MulModPrecon(x, B_to_km1, p, B_to_km1_pinv);
+         }
+         
+         // loop over remaining bits in the last word
+         word_t y = *next;
+         for (; bits > 0; bits--)
+         {
+            if (y & (((word_t) 1) << (WORD_BITS - 1)))
+               sum = SubMod(sum, x, p);
+            else
+               sum = AddMod(sum, x, p);
+
+            x = MulModPrecon(x, two_to_km1, p, two_to_km1_pinv);
+            y <<= 1;
+         }
+         
+         // update s
+         s = MulModPrecon(s, s_jump, p, s_jump_pinv);
+      }
+
+      // update g^i and (g^(k-1))^i
+      g_to_i = MulModPrecon(g_to_i, g, p, g_pinv);
+      g_to_km1_to_i = MulModPrecon(g_to_km1_to_i, g_to_km1, p, g_to_km1_pinv);
+   }
+
+   // =========== phase 2: consolidate table data
+
+   // compute weights[z] = sum((-1)^z[t] * (2^(k-1))^(TABLE_LG_SIZE - 1 - t) :
+   //                                                  0 <= t < TABLE_LG_SIZE).
+
+   long weights[TABLE_SIZE];
+   weights[0] = 0;
+   for (long h = 0, x = 1; h < TABLE_LG_SIZE;
+        h++, x = MulModPrecon(x, two_to_km1, p, two_to_km1_pinv))
+   {
+      for (long i = (1L << h) - 1; i >= 0; i--)
+      {
+         weights[2*i+1] = SubMod(weights[i], x, p);
+         weights[2*i]   = AddMod(weights[i], x, p);
+      }
+   }
+
+   // combine table data with weights
+
+   long x_jump = PowerMod(two_to_km1, TABLE_LG_SIZE, p, pinv);
+
+   for (long h = NUM_TABLES - 1, x = 1; h >= 0; h--)
+   {
+      mulmod_precon_t x_pinv = PrepMulModPrecon(x, p, pinv);
+
+      for (long i = 0; i < TABLE_SIZE; i++)
+      {
+         long y = MulMod(tables[h][i], weights[i], p, pinv);
+         y = MulModPrecon(y, x, p, x_pinv);
+         sum = SubMod(sum, y, p);
+      }
+
+      x = MulModPrecon(x_jump, x, p, x_pinv);
+   }
+
+   return sum;
+}
+
+
+/******************************************************************************
+
+   Computing the main sum (c = 1/2 case, with REDC arithmetic)
+   
+   Throughout this section F denotes 2^(LONG_BIT / 2).
+
+******************************************************************************/
+
+
+/*
+   Returns x/F mod n. Output is in [0, 2n), i.e. *not* reduced completely
+   into [0, n).
+
+   PRECONDITIONS:
+      3 <= n < F, n odd
+      0 <= x < nF    (if n < F/2)
+      0 <= x < nF/2  (if n > F/2)
+      ninv2 = -1/n mod F
+*/
+#define LOW_MASK ((1L << (LONG_BIT / 2)) - 1)
+static inline long RedcFast(long x, long n, long ninv2)
+{
+   unsigned long y = (x * ninv2) & LOW_MASK;
+   unsigned long z = x + (n * y);
+   return z >> (LONG_BIT / 2);
+}
+
+
+/*
+   Same as RedcFast(), but reduces output into [0, n).
+*/
+static inline long Redc(long x, long n, long ninv2)
+{
+   long y = RedcFast(x, n, ninv2);
+   if (y >= n)
+      y -= n;
+   return y;
+}
+
+
+/*
+   Computes -1/n mod F, in [0, F).
+   
+   PRECONDITIONS:
+      3 <= n < F, n odd
+*/
+long PrepRedc(long n)
+{
+   long ninv2 = -n;   // already correct mod 8
+   
+   // newton's method for 2-adic inversion
+   for (long bits = 3; bits < LONG_BIT/2; bits *= 2)
+      ninv2 = 2*ninv2 + n * ninv2 * ninv2;
+      
+   return ninv2 & LOW_MASK;
+}
+
+
+/*
+   Same as bernsum_pow2(), but uses REDC arithmetic, and various delayed
+   reduction strategies.
+
+   PRECONDITIONS:
+      Same as bernsum_pow2(), and in addition:
+      p < 2^(LONG_BIT/2 - 1)
+      
+   (See bernsum_pow2() for code comments; we only add comments here where
+   something is different from bernsum_pow2())
+*/
+long bernsum_pow2_redc(long p, double pinv, long k, long g, long n)
+{
+   long pinv2 = PrepRedc(p);
+   long F = (1L << (LONG_BIT/2)) % p;
+   
+   long tables[NUM_TABLES][TABLE_SIZE];
+   memset(tables, 0, sizeof(long) * NUM_TABLES * TABLE_SIZE);
+
+   long m = (p-1) / n;
+
+   if (n & 1)
+      m >>= 1;
+   else
+      n >>= 1;
+   
+   long g_to_km1 = PowerMod(g, k-1, p, pinv);
+   long two_to_km1 = PowerMod(2, k-1, p, pinv);
+   long B_to_km1 = PowerMod(two_to_km1, WORD_BITS, p, pinv);
+   long s_jump = PowerMod(2, MAX_INV * WORD_BITS, p, pinv);
+   
+   long g_redc = MulMod(g, F, p, pinv);
+   long g_to_km1_redc = MulMod(g_to_km1, F, p, pinv);
+   long two_to_km1_redc = MulMod(two_to_km1, F, p, pinv);
+   long B_to_km1_redc = MulMod(B_to_km1, F, p, pinv);
+   long s_jump_redc = MulMod(s_jump, F, p, pinv);
+
+   long g_to_km1_to_i = 1;    // always in [0, 2p)
+   long g_to_i = 1;           // always in [0, 2p)
+   long sum = 0;
+
+   Expander expander(p, (n >= MAX_INV * WORD_BITS)
+                                       ? MAX_INV : ((n - 1) / WORD_BITS + 1));
+
+   // =========== phase 1: main summation loop
+   
+   for (long i = 0; i < m; i++)
+   {
+      long s = g_to_i;           // always in [0, p)
+      if (s >= p)
+         s -= p;
+      
+      long x = g_to_km1_to_i;    // always in [0, 2p)
+      
+      for (long nn = n; nn > 0; nn -= MAX_INV * WORD_BITS)
+      {
+         word_t s_over_p[MAX_INV + 2];
+         long bits, words;
+         
+         if (nn >= MAX_INV * WORD_BITS)
+         {
+            bits = MAX_INV * WORD_BITS;
+            words = MAX_INV;
+         }
+         else
+         {
+            bits = nn;
+            words = (nn - 1) / WORD_BITS + 1;
+         }
+
+         expander.expand(s_over_p, s, words);
+         word_t* next = s_over_p + words;
+         
+         for (; bits >= WORD_BITS; bits -= WORD_BITS, next--)
+         {
+            word_t y = *next;
+            
+            // note: we add the values into tables *without* reduction mod p
+
+#if NUM_TABLES != 8 && NUM_TABLES != 4
+            // generic version
+            for (long h = 0; h < NUM_TABLES; h++)
+            {
+               tables[h][y & TABLE_MASK] += x;
+               y >>= TABLE_LG_SIZE;
+            }
+#else
+            // unrolled versions for 32-bit/64-bit machines
+            tables[0][ y                       & TABLE_MASK] += x;
+            tables[1][(y >>    TABLE_LG_SIZE ) & TABLE_MASK] += x;
+            tables[2][(y >> (2*TABLE_LG_SIZE)) & TABLE_MASK] += x;
+            tables[3][(y >> (3*TABLE_LG_SIZE)) & TABLE_MASK] += x;
+#if NUM_TABLES == 8
+            tables[4][(y >> (4*TABLE_LG_SIZE)) & TABLE_MASK] += x;
+            tables[5][(y >> (5*TABLE_LG_SIZE)) & TABLE_MASK] += x;
+            tables[6][(y >> (6*TABLE_LG_SIZE)) & TABLE_MASK] += x;
+            tables[7][(y >> (7*TABLE_LG_SIZE)) & TABLE_MASK] += x;
+#endif            
+#endif
+
+            x = RedcFast(x * B_to_km1_redc, p, pinv2);
+         }
+         
+         // bring x into [0, p) for next loop
+         if (x >= p)
+            x -= p;
+
+         word_t y = *next;
+         for (; bits > 0; bits--)
+         {
+            if (y & (((word_t) 1) << (WORD_BITS - 1)))
+               sum = SubMod(sum, x, p);
+            else
+               sum = AddMod(sum, x, p);
+
+            x = Redc(x * two_to_km1_redc, p, pinv2);
+            y <<= 1;
+         }
+         
+         s = Redc(s * s_jump_redc, p, pinv2);
+      }
+
+      g_to_i = RedcFast(g_to_i * g_redc, p, pinv2);
+      g_to_km1_to_i = RedcFast(g_to_km1_to_i * g_to_km1_redc, p, pinv2);
+   }
+   
+   // At this point, each table entry is at most p^2 (since x was always
+   // in [0, 2p), and the inner loop was called at most (p/2) / WORD_BITS
+   // times, and 2p * p/2 / WORD_BITS * TABLE_LG_SIZE <= p^2).
+   
+   // =========== phase 2: consolidate table data
+   
+   long weights[TABLE_SIZE];
+   weights[0] = 0;
+   // we store the weights multiplied by a factor of 2^(3*LONG_BIT/2) to
+   // compensate for the three rounds of REDC reduction in the loop below
+   for (long h = 0, x = PowerMod(2, 3*LONG_BIT/2, p, pinv);
+        h < TABLE_LG_SIZE; h++, x = Redc(x * two_to_km1_redc, p, pinv2))
+   {
+      for (long i = (1L << h) - 1; i >= 0; i--)
+      {
+         weights[2*i+1] = SubMod(weights[i], x, p);
+         weights[2*i]   = AddMod(weights[i], x, p);
+      }
+   }
+
+   long x_jump = PowerMod(two_to_km1, TABLE_LG_SIZE, p, pinv);
+   long x_jump_redc = MulMod(x_jump, F, p, pinv);
+
+   for (long h = NUM_TABLES - 1, x = 1; h >= 0; h--)
+   {
+      for (long i = 0; i < TABLE_SIZE; i++)
+      {
+         long y;
+         y = RedcFast(tables[h][i], p, pinv2);
+         y = RedcFast(y * weights[i], p, pinv2);
+         y = RedcFast(y * x, p, pinv2);
+         sum += y;
+      }
+
+      x = Redc(x * x_jump_redc, p, pinv2);
+   }
+
+   return sum % p;
+}
+
+
+
+/******************************************************************************
+
+   Wrappers for bernsum_*
+
+******************************************************************************/
+
+
+/*
+   Returns B_k/k mod p, in the range [0, p).
+   
+   PRECONDITIONS:
+      5 <= p < NTL_SP_BOUND, p prime
+      2 <= k <= p-3, k even
+      pinv = 1 / ((double) p)
+   
+   Algorithm: uses bernsum_powg() to compute the main sum.
+*/
+long _bern_modp_powg(long p, double pinv, long k)
+{
+   Factorisation F(p-1);
+   long g = primitive_root(p, pinv, F);
+
+   // compute main sum
+   long x = bernsum_powg(p, pinv, k, g);
+
+   // divide by (1 - g^k) and multiply by 2
+   long g_to_k = PowerMod(g, k, p, pinv);
+   long t = InvMod(p + 1 - g_to_k, p);
+   x = MulMod(x, t, p, pinv);
+   x = AddMod(x, x, p);
+   
+   return x;
+}
+
+
+/*
+   Returns B_k/k mod p, in the range [0, p).
+   
+   PRECONDITIONS:
+      5 <= p < NTL_SP_BOUND, p prime
+      2 <= k <= p-3, k even
+      pinv = 1 / ((double) p)
+      2^k != 1 mod p
+
+   Algorithm: uses bernsum_pow2() (or bernsum_pow2_redc() if p is small
+   enough) to compute the main sum.
+*/
+long _bern_modp_pow2(long p, double pinv, long k)
+{
+   Factorisation F(p-1);
+   long g = primitive_root(p, pinv, F);
+   long n = order(2, p, pinv, F);
+
+   // compute main sum
+   long x;
+   if (p < (1L << (LONG_BIT/2 - 1)))
+      x = bernsum_pow2_redc(p, pinv, k, g, n);
+   else
+      x = bernsum_pow2(p, pinv, k, g, n);
+
+   // divide by 2*(2^(-k) - 1)
+   long t = PowerMod(2, -k, p, pinv) - 1;
+   t = AddMod(t, t, p);
+   t = InvMod(t, p);
+   x = MulMod(x, t, p, pinv);
+   
+   return x;
+}
+
+
+
+/*
+   Returns B_k/k mod p, in the range [0, p).
+   
+   PRECONDITIONS:
+      5 <= p < NTL_SP_BOUND, p prime
+      2 <= k <= p-3, k even
+      pinv = 1 / ((double) p)
+*/
+long _bern_modp(long p, double pinv, long k)
+{
+   if (PowerMod(2, k, p, pinv) != 1)
+      // 2^k != 1 mod p, so we use the faster version
+      return _bern_modp_pow2(p, pinv, k);
+   else
+      // forced to use slower version
+      return _bern_modp_powg(p, pinv, k);
+}
+
+
+
+/******************************************************************************
+
+   Main bern_modp() routine
+
+******************************************************************************/
+
+long bern_modp(long p, long k)
+{
+   assert(k >= 0);
+   assert(2 <= p && p < NTL_SP_BOUND);
+
+   // B_0 = 1
+   if (k == 0)
+      return 1;
+
+   // B_1 = -1/2 mod p
+   if (k == 1)
+   {
+      if (p == 2)
+         return -1;
+      return (p-1)/2;
+   }
+
+   // B_k = 0 for odd k >= 3
+   if (k & 1)
+      return 0;
+
+   // denominator of B_k is always divisible by 6 for k >= 2
+   if (p <= 3)
+      return -1;
+      
+   // use Kummer's congruence (k = m mod p-1  =>  B_k/k = B_m/m mod p)
+   long m = k % (p-1);
+   if (m == 0)
+      return -1;
+
+   double pinv = 1 / ((double) p);
+   long x = _bern_modp(p, pinv, m);    // = B_m/m mod p
+   return MulMod(x, k, p, pinv);
+}
+
+
+};    // end namespace
+
+
+
+// end of file ================================================================

(a) /dev/null vs. (b) b/sage/rings/bernmm/bern_modp.h

@@ +1 @@
+/*
+   bern_modp.cpp:  computing isolated Bernoulli numbers modulo p
+   
+   Copyright (C) 2008, David Harvey
+   
+   This file is part of the bernmm package (version 1.0).
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef BERNMM_BERN_MODP_H
+#define BERNMM_BERN_MODP_H
+
+
+namespace bernmm {
+
+
+/*
+   Returns B_k mod p, in [0, p), or -1 if B_k is not p-integral.
+
+   PRECONDITIONS:
+      2 <= p < NTL_SP_BOUND, p prime
+      k >= 0
+*/
+long bern_modp(long p, long k);
+
+
+/*
+   Exported for testing.
+*/
+long _bern_modp_powg(long p, double pinv, long k);
+long _bern_modp_pow2(long p, double pinv, long k);
+
+
+};
+
+
+#endif
+
+// end of file ================================================================

(a) /dev/null vs. (b) b/sage/rings/bernmm/bern_modp_util.cpp

@@ +1 @@
+/*
+   bern_modp_util.cpp:  number-theoretic utility functions
+   
+   Copyright (C) 2008, David Harvey
+   
+   This file is part of the bernmm package (version 1.0).
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <NTL/ZZ.h>
+#include "bern_modp_util.h"
+
+
+NTL_CLIENT;
+
+
+namespace bernmm {
+
+
+long PowerMod(long a, long ee, long n, double ninv)
+{
+   long x, y;
+
+   unsigned long e;
+
+   if (ee < 0)
+      e = - ((unsigned long) ee);
+   else
+      e = ee;
+
+   x = 1;
+   y = a;
+   while (e) {
+      if (e & 1) x = MulMod(x, y, n, ninv);
+      y = MulMod(y, y, n, ninv);
+      e = e >> 1;
+   }
+
+   if (ee < 0) x = InvMod(x, n);
+
+   return x;
+}
+
+
+
+void Factorisation::helper(long k, long m)
+{
+   if (m == 1)
+      return;
+
+   for (long i = k + 1; i * i <= m; i++)
+   {
+      if (m % i == 0)
+      {
+         // found a factor
+         factors.push_back(i);
+         // remove that factor entirely
+         for (m /= i; m % i == 0; m /= i);
+         // recurse
+         helper(i, m);
+         return;
+      }
+   }
+
+   // no more factors
+   factors.push_back(m);
+}
+
+
+Factorisation::Factorisation(long n)
+{
+   this->n = n;
+   helper(1, n);
+}
+
+
+PrimeTable::PrimeTable(long bound)
+{
+   long size = (bound - 1) / LONG_BIT + 1;   // = ceil(bound / LONG_BIT)
+   data.resize(size);
+   
+   for (long i = 2; i * i < bound; i++)
+      if (is_prime(i))
+         for (long j = 2*i; j < bound; j += i)
+            set(j);
+}
+
+
+long order(long x, long p, double pinv, const Factorisation& F)
+{
+   // in the loop below, m is always some multiple of the order of x
+   long m = p - 1;
+
+   // try to remove factors from m until we can't remove any more
+   for (int i = 0; i < F.factors.size(); i++)
+   {
+      long q = F.factors[i];
+
+      while (m % q == 0)
+      {
+         long mm = m / q;
+         if (PowerMod(x, mm, p, pinv) != 1)
+            break;
+         m = mm;
+      }
+   }
+
+   return m;
+}
+
+
+
+long primitive_root(long p, double pinv, const Factorisation& F)
+{
+   if (p == 2)
+      return 1;
+
+   long g = 2;
+   for (; g < p; g++)
+      if (order(g, p, pinv, F) == p - 1)
+         return g;
+         
+   // no generator exists!?
+   abort();
+}
+
+
+
+};    // end namespace
+
+
+// end of file ================================================================

(a) /dev/null vs. (b) b/sage/rings/bernmm/bern_modp_util.h

@@ +1 @@
+/*
+   bern_modp_util.h:  number-theoretic utility functions
+   
+   Copyright (C) 2008, David Harvey
+   
+   This file is part of the bernmm package (version 1.0).
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef BERNMM_BERN_MODP_UTIL_H
+#define BERNMM_BERN_MODP_UTIL_H
+
+
+#include <vector>
+#include <cassert>
+
+
+namespace bernmm {
+
+
+/*
+   Same as NTL's PowerMod, but also accepts an _ninv_ parameter, which is the
+   same as the ninv parameter for NTL's MulMod routines, i.e. should have
+   ninv = 1 / ((double) n).
+   
+   (Implementation is adapted from ZZ.c in NTL 5.4.1.)
+*/
+long PowerMod(long a, long ee, long n, double ninv);
+
+
+/*
+   Represents the factorisation of an integer n into distinct prime factors.
+   
+   (Very naive implementation!)
+*/
+class Factorisation
+{
+protected:
+   /*
+      Finds distinct prime factors of m in the range k < p <= m.
+      Assumes that m does not have any prime factors p <= k.
+      Appends factors found to _factors_.
+   */
+   void helper(long k, long m);
+
+public:
+   // the integer
+   long n;
+
+   // the distinct factors (in increasing order)
+   std::vector<long> factors;
+
+   // initialises with given integer
+   Factorisation(long n);
+};
+
+
+
+class PrimeTable
+{
+private:
+   std::vector<long> data;   // bit-vector; 0 means prime, 1 means composite
+   
+   // read bit from index i
+   inline bool get(long i) const
+   {
+      return (data[i / LONG_BIT] >> (i % LONG_BIT)) & 1;
+   }
+
+   // set bit at index i
+   inline void set(long i)
+   {
+      data[i / LONG_BIT] |= (1L << (i % LONG_BIT));
+   }
+
+
+public:
+   // initialise with primes up to given bound
+   PrimeTable(long bound);
+   
+   // test whether n is prime by table lookup
+   inline bool is_prime(long n) const
+   {
+      return !get(n);
+   }
+   
+   // returns smallest prime p that is larger than n
+   long next_prime(long n) const
+   {
+      for (n++; !is_prime(n); n++);
+      return n;
+   }
+};
+
+
+
+/*
+   Returns 1 if n is prime.
+*/
+int is_prime(long n);
+
+
+/*
+   Returns smallest prime larger than p.
+*/
+long next_prime(long p);
+
+
+/*
+   Computes order of x mod p, given the factorisation F of p-1.
+*/
+long order(long x, long p, double pinv, const Factorisation& F);
+
+
+/*
+   Finds the smallest primitive root mod p, given the factorisation F of p-1.
+*/
+long primitive_root(long p, double pinv, const Factorisation& F);
+
+
+};    // end namespace
+
+
+#endif
+
+// end of file ================================================================

(a) /dev/null vs. (b) b/sage/rings/bernmm/bern_rat.cpp

@@ +1 @@
+/*
+   bern_rat.cpp:  multi-modular algorithm for computing Bernoulli numbers
+   
+   Copyright (C) 2008, David Harvey
+   
+   This file is part of the bernmm package (version 1.0).
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <gmp.h>
+#include <NTL/ZZ.h>
+#include <cmath>
+#include <vector>
+#include <set>
+#include "bern_modp_util.h"
+#include "bern_modp.h"
+#include "bern_rat.h"
+
+#ifdef USE_THREADS
+#include <pthread.h>
+#endif
+
+
+using namespace std;
+using namespace NTL;
+
+
+namespace bernmm {
+
+
+/*
+   Computes the denominator of B_k using Clausen/von Staudt.
+*/
+void bern_den(mpz_t res, long k, const PrimeTable& table)
+{
+   mpz_set_ui(res, 1);
+
+   // loop through factors of k
+   for (long f = 1; f*f <= k; f++)
+   {
+      // if f divides k....
+      if (k % f == 0)
+      {
+         // ... then both f + 1 and k/f + 1 are candidates for primes
+         // dividing the denominator of B_k
+         if (table.is_prime(f + 1))
+            mpz_mul_ui(res, res, f + 1);
+
+         if (f*f != k)
+            if (table.is_prime(k/f + 1))
+               mpz_mul_ui(res, res, k/f + 1);
+      }
+   }
+}
+
+
+// width of interval for each block
+#define BLOCK_SIZE 1000
+
+
+/*
+   Represents that B_k is congruent to _residue_ modulo _modulus_.
+*/
+struct Item
+{
+   mpz_t modulus;
+   mpz_t residue;
+   
+   Item()
+   {
+      mpz_init(modulus);
+      mpz_init(residue);
+   }
+   
+   ~Item()
+   {
+      mpz_clear(residue);
+      mpz_clear(modulus);
+   }
+};
+
+
+/*
+   Items get sorted by modulus.
+*/
+struct Item_cmp
+{
+   bool operator()(const Item* x, const Item* y)
+   {
+      return mpz_cmp(x->modulus, y->modulus) < 0;
+   }
+};
+
+
+/*
+   Returns new Item that combines information from op1 and op2 via CRT.
+*/
+Item* CRT(Item* op1, Item* op2)
+{
+   Item* res = new Item;
+
+   // let n1, n2 be the moduli, and r1, r2 be the residues
+   
+   // res->modulus = t, where t = 0 mod n1, t = 1 mod n2
+   mpz_invert(res->modulus, op1->modulus, op2->modulus);
+   mpz_mul(res->modulus, res->modulus, op1->modulus);
+
+   // res->residue = r2 - r1
+   mpz_sub(res->residue, op2->residue, op1->residue);
+   // res->residue = t * (r2 - r1)
+   mpz_mul(res->residue, res->residue, res->modulus);
+   // res->residue = r1 + t * (r2 - r1)
+   mpz_add(res->residue, res->residue, op1->residue);
+   // res->modulus = n1 * n2
+   mpz_mul(res->modulus, op1->modulus, op2->modulus);
+   // res->residue = r1 mod n1, r2 = mod n2
+   mpz_mod(res->residue, res->residue, res->modulus);
+   
+   return res;
+}
+
+
+struct State
+{
+   long k;
+   long bound;   // only use primes less than this bound
+   const PrimeTable* table;
+   
+   // index of block that should be processed next
+   long next;
+   
+   std::set<Item*, Item_cmp> items;
+#ifdef USE_THREADS
+   pthread_mutex_t lock;
+#endif
+   
+   State(long k, long bound, const PrimeTable& table)
+   {
+      this->k = k;
+      this->bound = bound;
+      this->next = 0;
+      this->table = &table;
+#ifdef USE_THREADS
+      pthread_mutex_init(&lock, NULL);
+#endif
+   }
+   
+   ~State()
+   {
+#ifdef USE_THREADS
+      pthread_mutex_destroy(&lock);
+#endif
+   }
+};
+
+
+void* worker(void* arg)
+{
+   State& state = *((State*) arg);
+   long k = state.k;
+
+#ifdef USE_THREADS
+   pthread_mutex_lock(&state.lock);
+#endif
+
+   while (1)
+   {
+      if (state.next * BLOCK_SIZE < state.bound)
+      {
+         // need to generate more modular data
+         
+         long next = state.next++;
+#ifdef USE_THREADS
+         pthread_mutex_unlock(&state.lock);
+#endif
+
+         Item* item = new Item;
+
+         mpz_set_ui(item->modulus, 1);
+         mpz_set_ui(item->residue, 0);
+         
+         for (long p = max(5, state.table->next_prime(next * BLOCK_SIZE));
+              p < state.bound && p < (next+1) * BLOCK_SIZE;
+              p = state.table->next_prime(p))
+         {
+            if (k % (p-1) == 0)
+               continue;
+
+            // compute B_k mod p
+            long b = bern_modp(p, k);
+            
+            // CRT into running total
+            long x = MulMod(SubMod(b, mpz_fdiv_ui(item->residue, p), p),
+                            InvMod(mpz_fdiv_ui(item->modulus, p), p), p);
+            mpz_addmul_ui(item->residue, item->modulus, x);
+            mpz_mul_ui(item->modulus, item->modulus, p);
+         }
+         
+#ifdef USE_THREADS
+         pthread_mutex_lock(&state.lock);
+#endif
+         state.items.insert(item);
+      }
+      else
+      {
+         // all modular data has been generated
+
+         if (state.items.size() <= 1)
+         {
+            // no more CRTs for this thread to perform
+#ifdef USE_THREADS
+            pthread_mutex_unlock(&state.lock);
+#endif
+            return NULL;
+         }
+         
+         // CRT two smallest items together
+         Item* item1 = *(state.items.begin());
+         state.items.erase(state.items.begin());
+         Item* item2 = *(state.items.begin());
+         state.items.erase(state.items.begin());
+#ifdef USE_THREADS
+         pthread_mutex_unlock(&state.lock);
+#endif
+         
+         Item* item3 = CRT(item1, item2);
+         delete item1;
+         delete item2;
+
+#ifdef USE_THREADS
+         pthread_mutex_lock(&state.lock);
+#endif
+         state.items.insert(item3);
+      }
+   }
+}
+
+
+void bern_rat(mpq_t res, long k, int num_threads)
+{
+   // special cases
+
+   if (k == 0)
+   {
+      // B_0 = 1
+      mpq_set_ui(res, 1, 1);
+      return;
+   }
+   
+   if (k == 1)
+   {
+      // B_1 = -1/2
+      mpq_set_si(res, -1, 2);
+      return;
+   }
+   
+   if (k == 2)
+   {
+      // B_2 = 1/6
+      mpq_set_si(res, 1, 6);
+      return;
+   }
+   
+   if (k & 1)
+   {
+      // B_k = 0 if k is odd
+      mpq_set_ui(res, 0, 1);
+      return;
+   }
+
+   if (num_threads <= 0)
+      num_threads = 1;
+
+   mpz_t num, den;
+   mpz_init(num);
+   mpz_init(den);
+   
+   const double log2 =    0.69314718055994528622676;
+   const double invlog2 = 1.44269504088896340735992;   // = 1/log(2)
+
+   // compute preliminary prime bound and build prime table
+   long bound1 = (long) max(37.0, ceil((k + 0.5) * log(k) * invlog2));
+   PrimeTable table(bound1);
+
+   // compute denominator of B_k
+   bern_den(den, k, table);
+
+   // compute number of bits we need to resolve the numerator
+   long bits = (long) ceil((k + 0.5) * log(k) * invlog2 - 4.094 * k + 2.470
+                                      + log(mpz_get_d(den)) * invlog2);
+
+   // compute tighter prime bound
+   // (note: we can safely get away with double-precision here. It would
+   // only start being insufficient around k = 10^13 or so, which is totally
+   // impractical at present.)
+   double prod = 1.0;
+   long prod_bits = 0;
+   long p;
+   for (p = 5; prod_bits < bits + 1; p = table.next_prime(p))
+   {
+      if (p >= NTL_SP_BOUND)
+         abort();   // !!!!! not sure what else we can do here...
+      if (k % (p-1) != 0)
+         prod *= (double) p;
+      int exp;
+      prod = frexp(prod, &exp);
+      prod_bits += exp;
+   }
+   long bound2 = p;
+
+   State state(k, bound2, table);
+
+#ifdef USE_THREADS
+   vector<pthread_t> threads(num_threads - 1);
+
+   // spawn worker threads to process blocks
+   for (long i = 0; i < num_threads - 1; i++)
+      pthread_create(&threads[i], NULL, worker, &state);
+#endif
+   
+   worker(&state);    // make this thread a worker too
+
+#ifdef USE_THREADS
+   for (long i = 0; i < num_threads - 1; i++)
+      pthread_join(threads[i], NULL);
+#endif
+      
+   // reconstruct B_k as a rational number
+   Item* item = *(state.items.begin());
+   mpz_mul(num, item->residue, den);
+   mpz_mod(num, num, item->modulus);
+   
+   if (k % 4 == 0)
+   {
+      // B_k is negative
+      mpz_sub(num, item->modulus, num);
+      mpz_neg(num, num);
+   }
+
+   delete item;
+   
+   mpz_swap(num, mpq_numref(res));
+   mpz_swap(den, mpq_denref(res));
+   
+   mpz_clear(num);
+   mpz_clear(den);
+}
+
+
+
+};    // end namespace
+
+
+// end of file ================================================================

(a) /dev/null vs. (b) b/sage/rings/bernmm/bern_rat.h

@@ +1 @@
+/*
+   bern_rat.h:  multi-modular algorithm for computing Bernoulli numbers
+   
+   Copyright (C) 2008, David Harvey
+   
+   This file is part of the bernmm package (version 1.0).
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef BERNMM_BERN_RAT_H
+#define BERNMM_BERN_RAT_H
+
+#include <gmp.h>
+
+
+namespace bernmm {
+
+
+/*
+   Returns B_k as a rational number, stored at _res_.
+   
+   k must be >= 0.
+   
+   Uses _num_threads_ threads. (If USE_THREADS is not #defined, the code just
+   uses one thread.)
+*/
+void bern_rat(mpq_t res, long k, int num_threads);
+
+
+
+};    // end namespace
+
+
+#endif
+
+// end of file ================================================================

(a) /dev/null vs. (b) b/sage/rings/bernmm/bernmm-test.cpp

@@ +1 @@
+/*
+   bernmm-test.cpp:  test module
+   
+   Copyright (C) 2008, David Harvey
+   
+   This file is part of the bernmm package (version 1.0).
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 2 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <iostream>
+#include <NTL/ZZ.h>
+#include <gmp.h>
+#include "bern_modp_util.h"
+#include "bern_modp.h"
+#include "bern_rat.h"
+
+
+NTL_CLIENT;
+
+
+using namespace bernmm;
+using namespace std;
+
+
+/*
+   Computes B_0, B_1, ..., B_{n-1} using naive algorithm, writes them to res.
+*/
+void bern_naive(mpq_t* res, long n)
+{
+   mpq_t t, u;
+   mpq_init(t);
+   mpq_init(u);
+
+   // compute res[j] = B_j / j! for 0 <= j < n
+   if (n > 0)
+      mpq_set_si(res[0], 1, 1);
+
+   for (long j = 1; j < n; j++)
+   {
+      mpq_set_si(res[j], 0, 1);
+      mpq_set_ui(t, 1, 1);
+      for (long k = 0; k < j; k++)
+      {
+         mpz_mul_ui(mpq_denref(t), mpq_denref(t), k + 2);
+         mpq_mul(u, res[j - 1 - k], t);
+         mpq_sub(res[j], res[j], u);
+      }
+   }
+
+   // multiply through by j! for 0 <= j < n
+   mpq_set_ui(t, 1, 1);
+   for (long j = 2; j < n; j++)
+   {
+      mpz_mul_ui(mpq_numref(t), mpq_numref(t), j);
+      mpq_mul(res[j], res[j], t);
+   }
+
+   mpq_clear(u);
+   mpq_clear(t);
+}
+
+
+/*
+   Tests _bern_modp_powg() for a given p and k by comparing against the
+   rational number B_k (must be supplied in b).
+   
+   Returns 1 on success.
+*/
+int testcase__bern_modp_powg(long p, long k, mpq_t b)
+{
+   double pinv = 1 / ((double) p);
+   
+   // compute B_k mod p using _bern_modp_powg()
+   long x = _bern_modp_powg(p, pinv, k);
+   x = MulMod(x, k, p, pinv);
+   
+   // compute B_k mod p from rational B_k
+   long y = mpz_fdiv_ui(mpq_numref(b), p);
+   long z = mpz_fdiv_ui(mpq_denref(b), p);
+   return y == MulMod(z, x, p, pinv);
+}
+
+
+
+/*
+   Tests _bern_modp_powg() by comparing against naive computation of B_k
+   (as a rational) for a range of small p and k.
+   
+   Returns 1 on success.
+*/
+int test__bern_modp_powg()
+{
+   int success = 1;
+
+   const long MAX = 300;
+   mpq_t bern[MAX];
+   
+   // compute B_k's as rational numbers using naive algorithm
+   for (long i = 0; i < MAX; i++)
+      mpq_init(bern[i]);
+   bern_naive(bern, MAX);
+
+   // try a range of k's
+   for (long k = 2; k < MAX && success; k += 2)
+   {
+      // try a range of small p's
+      for (long p = k + 3; p < 2*MAX && success; p += 2)
+      {
+         if (!ProbPrime(p))
+            continue;
+         success = success && testcase__bern_modp_powg(p, k, bern[k]);
+      }
+      
+      // try a single larger p
+      success = success && testcase__bern_modp_powg(1000003, k, bern[k]);
+   }
+
+   // if we're on a 32-bit machine, try a single example with p right near
+   // NTL's boundary (this is infeasible on a 64-bit machine)
+   if (NTL_SP_NBITS <= 32)
+   {
+      long p = NTL_SP_BOUND - 1;
+      while (!ProbPrime(p))
+         p--;
+
+      long k = (MAX/2)*2 - 2;
+      success = success && testcase__bern_modp_powg(p, k, bern[k]);
+   }
+
+   for (long i = 0; i < MAX; i++)
+      mpq_clear(bern[i]);
+   
+   return success;
+}
+
+
+
+/*
+   Tests _bern_modp_pow2() for a given p and k by comparing against result
+   from _bern_modp_powg().
+
+   Returns 1 on success.
+   
+   If 2^k = 1 mod p, then _bern_modp_pow2() won't work, so it just returns 1.
+*/
+int testcase__bern_modp_pow2(long p, long k)
+{
+   double pinv = 1 / ((double) p);
+
+   if (PowerMod(2, k, p, pinv) == 1)
+      return 1;
+
+   long x = _bern_modp_powg(p, pinv, k);
+   long y = _bern_modp_pow2(p, pinv, k);
+   
+   return x == y;
+}
+
+
+
+/*
+   Tests _bern_modp_pow2() by comparing against _bern_modp_powg() for
+   a range of p and k.
+   
+   Returns 1 on success.
+*/
+int test__bern_modp_pow2()
+{
+   int success = 1;
+
+   // exhaustive comparison over some small p and k
+   for (long p = 5; p < 2000 && success; p += 2)
+   {
+      if (!ProbPrime(p))
+         continue;
+
+      for (long k = 2; k <= p - 3 && success; k += 2)
+         success = success && testcase__bern_modp_pow2(p, k);
+   }
+
+   // a few larger values of p
+   for (long p = 1000000; p < 1030000; p++)
+   {
+      if (!ProbPrime(p))
+         continue;
+      
+      long k = 2 * (rand() % ((p-3)/2)) + 2;
+      success = success && testcase__bern_modp_pow2(p, k);
+   }
+   
+   // if we're on a 32-bit machine, try a single example with p right near
+   // NTL's boundary (this is infeasible on a 64-bit machine)
+   if (NTL_SP_NBITS <= 32)
+   {
+      long p = NTL_SP_BOUND - 1;
+      while (!ProbPrime(p))
+         p--;
+      success = success & testcase__bern_modp_pow2(p, 10);
+   }
+   
+   // try a few just below the REDC barrier
+   if (LONG_BIT == 32)
+   {
+      long boundary = 1L << (LONG_BIT/2 - 1);
+      for (long p = boundary - 1000; p < boundary && success; p++)
+      {
+         if (ProbPrime(p))
+         {
+            for (long trial = 0; trial < 1000 && success; trial++)
+            {
+               long k = 2 * (rand() % ((p-3)/2)) + 2;
+               success = success && testcase__bern_modp_pow2(p, k);
+            }
+         }
+      }
+   }
+   else
+   {
+      // on a 64-bit machine, only try one, since these are huge!
+      long p = 1L << (LONG_BIT/2 - 1);
+      while (!ProbPrime(p))
+         p--;
+      success = success && testcase__bern_modp_pow2(p, 10);
+   }
+
+   return success;
+}
+
+
+/*
+   Tests bern_rat() by comparing against the naive algorithm for several small
+   k, and testing against bern_modp() for a couple of larger k.
+   
+   Returns 1 on success.
+*/
+int test_bern_rat()
+{
+   int success = 1;
+
+   const long MAX = 300;
+   mpq_t bern[MAX];
+   
+   // compute B_k's as rational numbers using naive algorithm
+   for (long i = 0; i < MAX; i++)
+      mpq_init(bern[i]);
+   bern_naive(bern, MAX);
+
+   mpq_t x;
+   mpq_init(x);
+
+   // exhaustive test for small k
+   for (long k = 0; k < MAX && success; k++)
+   {
+      bern_rat(x, k, 4);    // try with 4 threads just for fun
+      success = success && mpq_equal(x, bern[k]);
+   }
+   
+   // try a few larger k
+   for (long i = 0; i < 50 && success; i++)
+   {
+      long k = ((random() % 20000) / 2) * 2;
+      bern_rat(x, k, 4);
+      
+      // compare with modular information
+      long p = 1000003;
+      long num = mpz_fdiv_ui(mpq_numref(x), p);
+      long den = mpz_fdiv_ui(mpq_denref(x), p);
+      success = success && (MulMod(bern_modp(p, k), den, p) == num);
+   }
+
+   mpq_clear(x);
+   for (long i = 0; i < MAX; i++)
+      mpq_clear(bern[i]);
+   
+   return success;
+}
+
+
+void report(int success)
+{
+   if (success)
+      cout << "ok" << endl;
+   else
+   {
+      cout << "failed!" << endl;
+      abort();
+   }
+}
+
+
+int main(int argc, char* argv[])
+{
+   if (argc == 1)
+   {
+      cout << "bernmm test module" << endl;
+      cout << endl;
+      cout << "   bernmm-test --test" << endl;
+      cout << "        runs test suite" << endl;
+      cout << "   bernmm-test --rational <k> <threads>" << endl;
+      cout << "        computes B_k with <threads> threads" << endl;
+      cout << "   bernmm-test --modular <p> <k>" << endl;
+      cout << "        computes B_k mod p" << endl;
+      return 0;
+   }
+   
+   if (!strcmp(argv[1], "--test"))
+   {
+      cout << "testing _bern_modp_powg()... " << flush;
+      report(test__bern_modp_powg());
+
+      cout << "testing _bern_modp_pow2()... " << flush;
+      report(test__bern_modp_pow2());
+
+      cout << "testing bern_rat()... " << flush;
+      report(test_bern_rat());
+   }
+   else if (!strcmp(argv[1], "--rational"))
+   {
+      if (argc <= 3)
+      {
+         cout << "not enough arguments" << endl;
+         return 0;
+      }
+      long k = atol(argv[2]);
+      long threads = atol(argv[3]);
+      mpq_t r;
+      mpq_init(r);
+      bern_rat(r, k, threads);
+      gmp_printf("%Zd/%Zd\n", mpq_numref(r), mpq_denref(r));
+      mpq_clear(r);
+   }
+   else if (!strcmp(argv[1], "--modular"))
+   {
+      if (argc <= 3)
+      {
+         cout << "not enough arguments" << endl;
+         return 0;
+      }
+      long p = atol(argv[2]);
+      long k = atol(argv[3]);
+      cout << bern_modp(p, k) << endl;
+   }
+   else
+   {
+      cout << "unknown command" << endl;
+   }
+
+   return 0;
+}
+
+
+// end of file ================================================================

sage/rings/bernoulli_mod_p.pyx

@@ -6 +6 @@
     - William Stein (2006-07-28): some touch up.  
     - David Harvey (2006-08-06): new, faster algorithm, also using faster NTL interface
     - David Harvey (2007-08-31): algorithm for a single bernoulli number mod p
+    - David Harvey (2008-06): added interface to bernmm, removed old code
 """
 
 #*****************************************************************************
@@ -29 +30 @@
 from sage.libs.ntl.ntl_ZZ_pX cimport ntl_ZZ_pX
 import sage.libs.pari.gen
 from sage.rings.integer_mod_ring import Integers
+from sage.rings.bernmm import bernmm_bern_modp
 
 
 
@@ -208 +210 @@
 
 
 
-def bernoulli_mod_p_single(int p, int k):
+def bernoulli_mod_p_single(long p, long k):
     r"""
     Returns the bernoulli number $B_k$ mod $p$.
 
+    If $B_k$ is not $p$-integral, an ArithmeticError is raised.
+
     INPUT:
         p -- integer, a prime
-        k -- even integer in the range $0 \leq k \leq p-3$
+        k -- non-negative integer
     
     OUTPUT:
         The $k$-th bernoulli number mod $p$.
-
-    ALGORITHM:
-        Uses the identity
-          $$ (1-g^k) B_k/k = 2\sum_{r=1}^{(p-1)/2} g^{r(k-1)} ( [g^r/p] - g [g^(r-1)/p] + (g-1)/2 ), $$
-        where $g$ is a primitive root mod $p$, and where square brackets
-        denote the fractional part. This identity can be derived from
-        Theorem 2.3, chapter 2 of Lang's book "Cyclotomic fields".
-    
-    PERFORMANCE:
-        Linear in $p$. In particular the running time doesn't depend on k.
-        
-        It's much faster than computing *all* bernoulli numbers by using
-        bernoulli_mod_p(). For p = 1000003, the latter takes about 3s on my
-        laptop, whereas this function takes only 0.06s.
-        
-        It may or may not be faster than computing literally bernoulli(k) % p,
-        depending on how big k and p are relative to each other. For example on
-        my laptop, computing bernoulli(2000) % p only takes 0.01s. But
-        computing bernoulli(100000) % p takes 40s, whereas this function still
-        takes only 0.06s.
 
     EXAMPLES:
         sage: bernoulli_mod_p_single(1009, 48)
@@ -256 +240 @@
         ValueError: p (=100) must be a prime
         
         sage: bernoulli_mod_p_single(19, 5)
-        Traceback (most recent call last):
-        ...
-        ValueError: k (=5) must be even
+        0
         
         sage: bernoulli_mod_p_single(19, 18)
         Traceback (most recent call last):
         ...
-        ValueError: k (=18) must be non-negative, and at most p-3
+        ArithmeticError: B_k is not integral at p
         
         sage: bernoulli_mod_p_single(19, -4)
         Traceback (most recent call last):
         ...
-        ValueError: k (=-4) must be non-negative, and at most p-3
+        ValueError: k must be non-negative
     
     Check results against bernoulli_mod_p:
     
@@ -299 +281 @@
          
     AUTHOR:
         -- David Harvey (2007-08-31)
+        -- David Harvey (2008-06): rewrote to use bernmm library
 
     """
     if p <= 2:
@@ -309 +292 @@
 
     R = Integers(p)
 
-    if k == 0:
-        return R(1)
-        
-    if k & 1:
-        raise ValueError, "k (=%s) must be even" % k
-        
-    if k < 0 or k > p-3:
-        raise ValueError, "k (=%s) must be non-negative, and at most p-3" % k
-    
-    g = R.multiplicative_generator()
-    cdef llong g_lift = g.lift()
-    cdef llong g_to_km1 = (g**(k-1)).lift()
-    cdef llong g_to_km1_pow = g_to_km1
-    cdef llong c = ((g-1)/2).lift()
-    cdef llong g_pow = 1
-    cdef llong g_pow_new, quot
-    cdef llong sum = 0
-    cdef int r
+    cdef long x = bernmm_bern_modp(p, k)
+    if x == -1:
+        raise ArithmeticError, "B_k is not integral at p"
+    return x
 
-    for r from 0 <= r < (p-1)/2:
-        g_pow_new = g_pow * g_lift
-        quot = g_pow_new / p
-        sum = sum + g_to_km1_pow * (p + c - quot)
-        sum = sum % p
-        g_pow = g_pow_new % p
-        g_to_km1_pow = (g_to_km1_pow * g_to_km1) % p
-    
-    return R(sum) * 2 * k / (1 - g**k)
-    
 
 # ============ end of file

setup.py

@@ -812 +812 @@
               libraries=['ntl','stdc++'],
               language = 'c++',
               include_dirs=debian_include_dirs + ['sage/libs/ntl/']), \
+
+    Extension('sage.rings.bernmm',
+              sources = ['sage/rings/bernmm.pyx',
+                         'sage/rings/bernmm/bern_modp.cpp',
+                         'sage/rings/bernmm/bern_modp_util.cpp',
+                         'sage/rings/bernmm/bern_rat.cpp'],
+              libraries = ['gmp', 'ntl', 'stdc++', 'pthread'],
+              language = 'c++',
+              define_macros=[('USE_THREADS', '1')]), \
 
     Extension('sage.schemes.hyperelliptic_curves.hypellfrob',
                  sources = ['sage/schemes/hyperelliptic_curves/hypellfrob.pyx',