1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
|
--- ./src/sage/rings/bernmm/bernmm-test.cpp.orig 2015-02-16 17:15:12.000000000 -0700
+++ ./src/sage/rings/bernmm/bernmm-test.cpp 2015-05-07 21:39:58.565251320 -0600
@@ -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);
+ wide_double pinv = wide_double(1) / wide_double(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);
+ wide_double pinv = wide_double(1) / wide_double(p);
if (PowerMod(2, k, p, pinv) == 1)
return 1;
--- ./src/sage/rings/bernmm/bern_modp.cpp.orig 2015-02-16 17:15:12.000000000 -0700
+++ ./src/sage/rings/bernmm/bern_modp.cpp 2015-05-07 20:17:37.680381004 -0600
@@ -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, wide_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);
+ wide_double g_pinv = wide_double(g) / wide_double(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, wide_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
@@ -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, wide_double 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, wide_double 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, wide_double pinv, long k)
{
Factorisation F(p-1);
long g = primitive_root(p, pinv, F);
@@ -717,7 +717,7 @@ long _bern_modp_pow2(long p, double pinv
2 <= k <= p-3, k even
pinv = 1 / ((double) p)
*/
-long _bern_modp(long p, double pinv, long k)
+long _bern_modp(long p, wide_double pinv, long k)
{
if (PowerMod(2, k, p, pinv) != 1)
// 2^k != 1 mod p, so we use the faster version
@@ -765,7 +765,7 @@ long bern_modp(long p, long k)
if (m == 0)
return -1;
- double pinv = 1 / ((double) p);
+ wide_double pinv = wide_double(1) / wide_double (p);
long x = _bern_modp(p, pinv, m); // = B_m/m mod p
return MulMod(x, k, p, pinv);
}
--- ./src/sage/rings/bernmm/bern_modp.h.orig 2015-02-16 17:15:12.000000000 -0700
+++ ./src/sage/rings/bernmm/bern_modp.h 2015-05-09 08:06:39.732529882 -0600
@@ -12,6 +12,7 @@
#ifndef BERNMM_BERN_MODP_H
#define BERNMM_BERN_MODP_H
+#include <NTL/ZZ.h>
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::wide_double pinv, long k);
+long _bern_modp_pow2(long p, NTL::wide_double pinv, long k);
};
--- ./src/sage/rings/bernmm/bern_modp_util.cpp.orig 2015-02-16 17:15:12.000000000 -0700
+++ ./src/sage/rings/bernmm/bern_modp_util.cpp 2015-05-07 21:38:06.662182003 -0600
@@ -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, wide_double 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, wide_double 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, wide_double pinv, const Factorisation& F)
{
if (p == 2)
return 1;
--- ./src/sage/rings/bernmm/bern_modp_util.h.orig 2015-02-16 17:15:12.000000000 -0700
+++ ./src/sage/rings/bernmm/bern_modp_util.h 2015-05-09 08:58:22.618458475 -0600
@@ -17,6 +17,7 @@
#include <vector>
#include <cassert>
#include <climits>
+#include <NTL/ZZ.h>
#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::wide_double 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::wide_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);
+long primitive_root(long p, NTL::wide_double pinv, const Factorisation& F);
}; // end namespace
|
