diff -up src/sage/rings/bernmm/bernmm-test.cpp.orig src/sage/rings/bernmm/bernmm-test.cpp --- src/sage/rings/bernmm/bernmm-test.cpp.orig 2015-10-11 18:17:42.808860675 -0300 +++ src/sage/rings/bernmm/bernmm-test.cpp 2015-10-11 18:18:44.684863044 -0300 @@ -70,7 +70,7 @@ void bern_naive(mpq_t* res, long n) */ int testcase__bern_modp_powg(long p, long k, mpq_t b) { - double pinv = 1 / ((double) p); + mulmod_t pinv = PrepMulMod(p); // compute B_k mod p using _bern_modp_powg() long x = _bern_modp_powg(p, pinv, k); @@ -147,7 +147,7 @@ int test__bern_modp_powg() */ int testcase__bern_modp_pow2(long p, long k) { - double pinv = 1 / ((double) p); + mulmod_t pinv = PrepMulMod(p); if (PowerMod(2, k, p, pinv) == 1) return 1; diff -up src/sage/rings/bernmm/bern_modp.cpp.orig src/sage/rings/bernmm/bern_modp.cpp --- src/sage/rings/bernmm/bern_modp.cpp.orig 2015-10-11 18:17:42.814860675 -0300 +++ src/sage/rings/bernmm/bern_modp.cpp 2015-10-11 18:20:28.077867003 -0300 @@ -43,14 +43,14 @@ namespace bernmm { 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 bernsum_powg(long p, mulmod_t 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); + muldivrem_t g_pinv = PrepMulDivRem(g, p); mulmod_precon_t g_to_km1_pinv = PrepMulModPrecon(g_to_km1, p, pinv); for (long j = 1; j <= (p-1)/2; j++) @@ -224,7 +224,7 @@ public: #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) +long bernsum_pow2(long p, mulmod_t 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 @@ -481,7 +481,7 @@ long PrepRedc(long n) (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 bernsum_pow2_redc(long p, mulmod_t pinv, long k, long g, long n) { long pinv2 = PrepRedc(p); long F = (1L << (ULONG_BITS/2)) % p; @@ -655,7 +655,7 @@ long bernsum_pow2_redc(long p, double pi Algorithm: uses bernsum_powg() to compute the main sum. */ -long _bern_modp_powg(long p, double pinv, long k) +long _bern_modp_powg(long p, mulmod_t pinv, long k) { Factorisation F(p-1); long g = primitive_root(p, pinv, F); @@ -685,7 +685,7 @@ long _bern_modp_powg(long p, double pinv 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) +long _bern_modp_pow2(long p, mulmod_t pinv, long k) { Factorisation F(p-1); long g = primitive_root(p, pinv, F); @@ -765,7 +765,7 @@ long bern_modp(long p, long k) if (m == 0) return -1; - double pinv = 1 / ((double) p); + mulmod_t pinv = PrepMulMod(p); long x = _bern_modp(p, pinv, m); // = B_m/m mod p return MulMod(x, k, p, pinv); } diff -up src/sage/rings/bernmm/bern_modp.h.orig src/sage/rings/bernmm/bern_modp.h --- src/sage/rings/bernmm/bern_modp.h.orig 2015-10-11 18:17:42.820860675 -0300 +++ src/sage/rings/bernmm/bern_modp.h 2015-10-11 18:20:53.453867975 -0300 @@ -12,6 +12,7 @@ #ifndef BERNMM_BERN_MODP_H #define BERNMM_BERN_MODP_H +#include namespace bernmm { @@ -29,8 +30,8 @@ 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); +long _bern_modp_powg(long p, NTL::mulmod_t pinv, long k); +long _bern_modp_pow2(long p, NTL::mulmod_t pinv, long k); }; diff -up src/sage/rings/bernmm/bern_modp_util.cpp.orig src/sage/rings/bernmm/bern_modp_util.cpp --- src/sage/rings/bernmm/bern_modp_util.cpp.orig 2015-10-11 18:17:42.825860675 -0300 +++ src/sage/rings/bernmm/bern_modp_util.cpp 2015-10-11 18:21:24.653869170 -0300 @@ -20,7 +20,7 @@ NTL_CLIENT; namespace bernmm { -long PowerMod(long a, long ee, long n, double ninv) +long PowerMod(long a, long ee, long n, mulmod_t ninv) { long x, y; @@ -89,7 +89,7 @@ PrimeTable::PrimeTable(long bound) } -long order(long x, long p, double pinv, const Factorisation& F) +long order(long x, long p, mulmod_t pinv, const Factorisation& F) { // in the loop below, m is always some multiple of the order of x long m = p - 1; @@ -113,7 +113,7 @@ long order(long x, long p, double pinv, -long primitive_root(long p, double pinv, const Factorisation& F) +long primitive_root(long p, mulmod_t pinv, const Factorisation& F) { if (p == 2) return 1; diff -up src/sage/rings/bernmm/bern_modp_util.h.orig src/sage/rings/bernmm/bern_modp_util.h --- src/sage/rings/bernmm/bern_modp_util.h.orig 2015-10-11 18:17:42.830860676 -0300 +++ src/sage/rings/bernmm/bern_modp_util.h 2015-10-11 18:21:58.044870449 -0300 @@ -17,6 +17,7 @@ #include #include #include +#include #if ULONG_MAX == 4294967295U @@ -39,7 +40,7 @@ namespace bernmm { (Implementation is adapted from ZZ.c in NTL 5.4.1.) */ -long PowerMod(long a, long ee, long n, double ninv); +long PowerMod(long a, long ee, long n, NTL::mulmod_t ninv); /* @@ -123,13 +124,13 @@ 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); +long order(long x, long p, NTL::mulmod_t 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); +long primitive_root(long p, NTL::mulmod_t pinv, const Factorisation& F); }; // end namespace diff -up src/sage/rings/bernmm/bern_modp.cpp.orig src/sage/rings/bernmm/bern_modp.cpp --- src/sage/rings/bernmm/bern_modp.cpp.orig 2015-11-02 23:58:56.669503117 -0200 +++ src/sage/rings/bernmm/bern_modp.cpp 2015-11-02 23:59:15.683503846 -0200 @@ -717,7 +717,7 @@ long _bern_modp_pow2(long p, mulmod_t pi 2 <= k <= p-3, k even pinv = 1 / ((double) p) */ -long _bern_modp(long p, double pinv, long k) +long _bern_modp(long p, mulmod_t pinv, long k) { if (PowerMod(2, k, p, pinv) != 1) // 2^k != 1 mod p, so we use the faster version