Ticket #3542: bernmm.2.patch

sage/rings/arith.py

@@ -145 +145 @@
 
 algebraic_dependency = algdep
 
-def bernoulli(n, algorithm='pari', num_threads=1):
+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)
@@ -192 +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)

sage/rings/bernmm.pyx

@@ -29 +29 @@
 
 
 def bernmm_bern_rat(long k, int num_threads = 1):
-    """
+    r"""
     Computes k-th Bernoulli number using a multimodular algorithm.
     (Wrapper for bernmm library.)
 
@@ -39 +39 @@
 
     COMPLEXITY:
         Pretty much quadratic in $k$. See the paper ``A multimodular algorithm
-        for computing isolated Bernoulli numbers'', David Harvey, 2008, for
-        more details.
+        for computing Bernoulli numbers'', David Harvey, 2008, for more details.
 
     EXAMPLES:
         sage: from sage.rings.bernmm import bernmm_bern_rat
@@ -73 +72 @@
 
 
 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.

sage/rings/bernmm/bern_modp.cpp

@@ -436 +436 @@
 
    PRECONDITIONS:
       3 <= n < F, n odd
-      0 <= x < nF
+      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)
@@ -638 +639 @@
          sum += y;
       }
 
-      x = RedcFast(x * x_jump_redc, p, pinv2);
+      x = Redc(x * x_jump_redc, p, pinv2);
    }
 
    return sum % p;

sage/rings/bernmm/bernmm-test.cpp

@@ -189 +189 @@
       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++)
    {
@@ -210 +210 @@
       success = success & testcase__bern_modp_pow2(p, 10);
    }
    
-   // one just below the REDC barrier
+   // 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);
-   }
-
-   // and one above the REDC barrier
-   {
-      long p = 1L << (LONG_BIT/2 - 1);
-      while (!ProbPrime(p))
-         p++;
       success = success && testcase__bern_modp_pow2(p, 10);
    }