[2cfffe] | 1 | |
---|
| 2 | #include <NTL/lzz_p.h> |
---|
| 3 | |
---|
| 4 | #include <NTL/new.h> |
---|
| 5 | |
---|
| 6 | NTL_START_IMPL |
---|
| 7 | |
---|
| 8 | zz_pInfoT::zz_pInfoT(long NewP, long maxroot) |
---|
| 9 | { |
---|
| 10 | ref_count = 1; |
---|
| 11 | |
---|
[09da99] | 12 | if (maxroot < 0) Error("zz_pContext: maxroot may not be negative"); |
---|
| 13 | |
---|
[2cfffe] | 14 | if (NewP <= 1) Error("zz_pContext: p must be > 1"); |
---|
| 15 | if (NumBits(NewP) > NTL_SP_NBITS) Error("zz_pContext: modulus too big"); |
---|
| 16 | |
---|
| 17 | ZZ P, B, M, M1, MinusM; |
---|
| 18 | long n, i; |
---|
| 19 | long q, t; |
---|
| 20 | |
---|
| 21 | p = NewP; |
---|
| 22 | |
---|
| 23 | pinv = 1/double(p); |
---|
| 24 | |
---|
| 25 | index = -1; |
---|
| 26 | |
---|
| 27 | conv(P, p); |
---|
| 28 | |
---|
| 29 | sqr(B, P); |
---|
| 30 | LeftShift(B, B, maxroot+NTL_FFTFudge); |
---|
| 31 | |
---|
| 32 | set(M); |
---|
| 33 | n = 0; |
---|
| 34 | while (M <= B) { |
---|
| 35 | UseFFTPrime(n); |
---|
| 36 | q = FFTPrime[n]; |
---|
| 37 | n++; |
---|
| 38 | mul(M, M, q); |
---|
| 39 | } |
---|
| 40 | |
---|
| 41 | if (n > 4) Error("zz_pInit: too many primes"); |
---|
| 42 | |
---|
| 43 | NumPrimes = n; |
---|
| 44 | PrimeCnt = n; |
---|
| 45 | MaxRoot = CalcMaxRoot(q); |
---|
| 46 | |
---|
| 47 | if (maxroot < MaxRoot) |
---|
| 48 | MaxRoot = maxroot; |
---|
| 49 | |
---|
| 50 | negate(MinusM, M); |
---|
| 51 | MinusMModP = rem(MinusM, p); |
---|
| 52 | |
---|
[09da99] | 53 | if (!(CoeffModP = (long *) NTL_MALLOC(n, sizeof(long), 0))) |
---|
[2cfffe] | 54 | Error("out of space"); |
---|
| 55 | |
---|
[09da99] | 56 | if (!(x = (double *) NTL_MALLOC(n, sizeof(double), 0))) |
---|
[2cfffe] | 57 | Error("out of space"); |
---|
| 58 | |
---|
[09da99] | 59 | if (!(u = (long *) NTL_MALLOC(n, sizeof(long), 0))) |
---|
[2cfffe] | 60 | Error("out of space"); |
---|
| 61 | |
---|
| 62 | for (i = 0; i < n; i++) { |
---|
| 63 | q = FFTPrime[i]; |
---|
| 64 | |
---|
| 65 | div(M1, M, q); |
---|
| 66 | t = rem(M1, q); |
---|
| 67 | t = InvMod(t, q); |
---|
| 68 | mul(M1, M1, t); |
---|
| 69 | CoeffModP[i] = rem(M1, p); |
---|
| 70 | x[i] = ((double) t)/((double) q); |
---|
| 71 | u[i] = t; |
---|
| 72 | } |
---|
| 73 | } |
---|
| 74 | |
---|
| 75 | zz_pInfoT::zz_pInfoT(long Index) |
---|
| 76 | { |
---|
| 77 | ref_count = 1; |
---|
| 78 | |
---|
| 79 | index = Index; |
---|
| 80 | |
---|
| 81 | if (index < 0) |
---|
| 82 | Error("bad FFT prime index"); |
---|
| 83 | |
---|
| 84 | // allows non-consecutive indices...I'm not sure why |
---|
| 85 | while (NumFFTPrimes < index) |
---|
| 86 | UseFFTPrime(NumFFTPrimes); |
---|
| 87 | |
---|
| 88 | UseFFTPrime(index); |
---|
| 89 | |
---|
| 90 | p = FFTPrime[index]; |
---|
| 91 | pinv = FFTPrimeInv[index]; |
---|
| 92 | |
---|
| 93 | NumPrimes = 1; |
---|
| 94 | PrimeCnt = 0; |
---|
| 95 | |
---|
| 96 | MaxRoot = CalcMaxRoot(p); |
---|
| 97 | } |
---|
| 98 | |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | |
---|
| 102 | zz_pInfoT::~zz_pInfoT() |
---|
| 103 | { |
---|
| 104 | if (index < 0) { |
---|
| 105 | free(CoeffModP); |
---|
| 106 | free(x); |
---|
| 107 | free(u); |
---|
| 108 | } |
---|
| 109 | } |
---|
| 110 | |
---|
| 111 | |
---|
| 112 | zz_pInfoT *zz_pInfo = 0; |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | typedef zz_pInfoT *zz_pInfoPtr; |
---|
| 116 | |
---|
[287cc8] | 117 | static |
---|
[2cfffe] | 118 | void CopyPointer(zz_pInfoPtr& dst, zz_pInfoPtr src) |
---|
| 119 | { |
---|
| 120 | if (src == dst) return; |
---|
| 121 | |
---|
| 122 | if (dst) { |
---|
| 123 | dst->ref_count--; |
---|
| 124 | |
---|
[287cc8] | 125 | if (dst->ref_count < 0) |
---|
[2cfffe] | 126 | Error("internal error: negative zz_pContext ref_count"); |
---|
| 127 | |
---|
| 128 | if (dst->ref_count == 0) delete dst; |
---|
| 129 | } |
---|
| 130 | |
---|
| 131 | if (src) { |
---|
[287cc8] | 132 | if (src->ref_count == NTL_MAX_LONG) |
---|
[2cfffe] | 133 | Error("internal error: zz_pContext ref_count overflow"); |
---|
[09da99] | 134 | |
---|
| 135 | src->ref_count++; |
---|
[2cfffe] | 136 | } |
---|
| 137 | |
---|
| 138 | dst = src; |
---|
| 139 | } |
---|
[287cc8] | 140 | |
---|
[2cfffe] | 141 | |
---|
| 142 | void zz_p::init(long p, long maxroot) |
---|
| 143 | { |
---|
| 144 | zz_pContext c(p, maxroot); |
---|
| 145 | c.restore(); |
---|
| 146 | |
---|
| 147 | } |
---|
| 148 | |
---|
| 149 | void zz_p::FFTInit(long index) |
---|
| 150 | { |
---|
| 151 | zz_pContext c(INIT_FFT, index); |
---|
| 152 | c.restore(); |
---|
| 153 | } |
---|
| 154 | |
---|
| 155 | zz_pContext::zz_pContext(long p, long maxroot) |
---|
| 156 | { |
---|
| 157 | ptr = NTL_NEW_OP zz_pInfoT(p, maxroot); |
---|
| 158 | } |
---|
| 159 | |
---|
| 160 | zz_pContext::zz_pContext(INIT_FFT_TYPE, long index) |
---|
| 161 | { |
---|
| 162 | ptr = NTL_NEW_OP zz_pInfoT(index); |
---|
| 163 | } |
---|
| 164 | |
---|
| 165 | zz_pContext::zz_pContext(const zz_pContext& a) |
---|
| 166 | { |
---|
| 167 | ptr = 0; |
---|
| 168 | CopyPointer(ptr, a.ptr); |
---|
| 169 | } |
---|
| 170 | |
---|
| 171 | |
---|
| 172 | zz_pContext& zz_pContext::operator=(const zz_pContext& a) |
---|
| 173 | { |
---|
| 174 | CopyPointer(ptr, a.ptr); |
---|
| 175 | return *this; |
---|
| 176 | } |
---|
| 177 | |
---|
| 178 | |
---|
| 179 | zz_pContext::~zz_pContext() |
---|
| 180 | { |
---|
| 181 | CopyPointer(ptr, 0); |
---|
| 182 | } |
---|
| 183 | |
---|
| 184 | void zz_pContext::save() |
---|
| 185 | { |
---|
| 186 | CopyPointer(ptr, zz_pInfo); |
---|
| 187 | } |
---|
| 188 | |
---|
| 189 | void zz_pContext::restore() const |
---|
| 190 | { |
---|
| 191 | CopyPointer(zz_pInfo, ptr); |
---|
| 192 | } |
---|
| 193 | |
---|
| 194 | |
---|
| 195 | |
---|
| 196 | zz_pBak::~zz_pBak() |
---|
| 197 | { |
---|
| 198 | if (MustRestore) |
---|
| 199 | CopyPointer(zz_pInfo, ptr); |
---|
| 200 | |
---|
| 201 | CopyPointer(ptr, 0); |
---|
| 202 | } |
---|
| 203 | |
---|
| 204 | void zz_pBak::save() |
---|
| 205 | { |
---|
| 206 | MustRestore = 1; |
---|
| 207 | CopyPointer(ptr, zz_pInfo); |
---|
| 208 | } |
---|
| 209 | |
---|
| 210 | |
---|
| 211 | void zz_pBak::restore() |
---|
| 212 | { |
---|
| 213 | MustRestore = 0; |
---|
| 214 | CopyPointer(zz_pInfo, ptr); |
---|
| 215 | } |
---|
| 216 | |
---|
| 217 | |
---|
| 218 | |
---|
| 219 | |
---|
[287cc8] | 220 | static inline |
---|
[09da99] | 221 | long reduce(long a, long p) |
---|
[2cfffe] | 222 | { |
---|
| 223 | if (a >= 0 && a < p) |
---|
| 224 | return a; |
---|
| 225 | else { |
---|
| 226 | a = a % p; |
---|
| 227 | if (a < 0) a += p; |
---|
| 228 | return a; |
---|
| 229 | } |
---|
| 230 | } |
---|
| 231 | |
---|
| 232 | |
---|
| 233 | zz_p to_zz_p(long a) |
---|
| 234 | { |
---|
| 235 | return zz_p(reduce(a, zz_p::modulus()), INIT_LOOP_HOLE); |
---|
| 236 | } |
---|
| 237 | |
---|
| 238 | void conv(zz_p& x, long a) |
---|
| 239 | { |
---|
| 240 | x._zz_p__rep = reduce(a, zz_p::modulus()); |
---|
| 241 | } |
---|
| 242 | |
---|
| 243 | zz_p to_zz_p(const ZZ& a) |
---|
| 244 | { |
---|
| 245 | return zz_p(rem(a, zz_p::modulus()), INIT_LOOP_HOLE); |
---|
| 246 | } |
---|
| 247 | |
---|
| 248 | void conv(zz_p& x, const ZZ& a) |
---|
| 249 | { |
---|
| 250 | x._zz_p__rep = rem(a, zz_p::modulus()); |
---|
| 251 | } |
---|
| 252 | |
---|
| 253 | |
---|
| 254 | NTL_END_IMPL |
---|