1
2 #ifndef NTL_zz_pX__H
3 #define NTL_zz_pX__H
4
5 #include <NTL/vector.h>
6 #include <NTL/lzz_p.h>
7 #include <NTL/vec_lzz_p.h>
8 #include <NTL/Lazy.h>
9 #include <NTL/SmartPtr.h>
10 #include <NTL/mat_lzz_p.h>
11
12 NTL_OPEN_NNS
13
14 // some cross-over points
15
16 #define NTL_zz_pX_MOD_CROSSOVER (zz_pX_mod_crossover[zz_pInfo->PrimeCnt])
17 #define NTL_zz_pX_MUL_CROSSOVER (zz_pX_mul_crossover[zz_pInfo->PrimeCnt])
18 #define NTL_zz_pX_NEWTON_CROSSOVER (zz_pX_newton_crossover[zz_pInfo->PrimeCnt])
19 #define NTL_zz_pX_DIV_CROSSOVER (zz_pX_div_crossover[zz_pInfo->PrimeCnt])
20 #define NTL_zz_pX_HalfGCD_CROSSOVER (zz_pX_halfgcd_crossover[zz_pInfo->PrimeCnt])
21 #define NTL_zz_pX_GCD_CROSSOVER (zz_pX_gcd_crossover[zz_pInfo->PrimeCnt])
22 #define NTL_zz_pX_BERMASS_CROSSOVER (zz_pX_bermass_crossover[zz_pInfo->PrimeCnt])
23 #define NTL_zz_pX_TRACE_CROSSOVER (zz_pX_trace_crossover[zz_pInfo->PrimeCnt])
24
25 extern const long zz_pX_mod_crossover[];
26 extern const long zz_pX_mul_crossover[];
27 extern const long zz_pX_newton_crossover[];
28 extern const long zz_pX_div_crossover[];
29 extern const long zz_pX_halfgcd_crossover[];
30 extern const long zz_pX_gcd_crossover[];
31 extern const long zz_pX_bermass_crossover[];
32 extern const long zz_pX_trace_crossover[];
33
34
35
36 /************************************************************
37
38 zz_pX
39
40 The class zz_pX implements polynomial arithmetic modulo p.
41 Polynomials are represented as vec_zz_p's.
42 If f is a zz_pX, then f.rep is a vec_zz_p.
43 The zero polynomial is represented as a zero length vector.
44 Otherwise. f.rep[0] is the constant-term, and f.rep[f.rep.length()-1]
45 is the leading coefficient, which is always non-zero.
46 The member f.rep is public, so the vector representation is fully
47 accessible.
48 Use the member function normalize() to strip leading zeros.
49
50 **************************************************************/
51
52
53 class zz_pE; // forward declaration
54 class zz_pXModulus;
55 class fftRep;
56 class zz_pXMultiplier;
57
58
59 class zz_pX {
60 public:
61 typedef zz_p coeff_type;
62 typedef zz_pE residue_type;
63 typedef zz_pXModulus modulus_type;
64 typedef zz_pXMultiplier multiplier_type;
65 typedef fftRep fft_type;
66
67
68 vec_zz_p rep;
69
70 typedef vec_zz_p VectorBaseType;
71
72
73
74 /***************************************************************
75
76 Constructors, Destructors, and Assignment
77
78 ****************************************************************/
79
80
zz_pX()81 zz_pX() {}
82 // initial value 0
83
zz_pX(long a)84 explicit zz_pX(long a) { *this = a; }
zz_pX(zz_p a)85 explicit zz_pX(zz_p a) { *this = a; }
86
87
zz_pX(INIT_SIZE_TYPE,long n)88 zz_pX(INIT_SIZE_TYPE, long n) { rep.SetMaxLength(n); }
89
90
91 inline zz_pX(long i, zz_p c);
92 inline zz_pX(long i, long c);
93
94 inline zz_pX(INIT_MONO_TYPE, long i, zz_p c);
95 inline zz_pX(INIT_MONO_TYPE, long i, long c);
96 inline zz_pX(INIT_MONO_TYPE, long i);
97
98
99 inline zz_pX& operator=(long a);
100 inline zz_pX& operator=(zz_p a);
101
102 // default copy constructor and assignment
103 // default destructor
104
105 void normalize();
106 // strip leading zeros
107
SetMaxLength(long n)108 void SetMaxLength(long n)
109 // pre-allocate space for n coefficients.
110 // Value is unchanged
111
112 { rep.SetMaxLength(n); }
113
114
kill()115 void kill()
116 // free space held by this polynomial. Value becomes 0.
117
118 { rep.kill(); }
119
120
121
SetLength(long n)122 void SetLength(long n) { rep.SetLength(n); }
123 zz_p& operator[](long i) { return rep[i]; }
124 const zz_p& operator[](long i) const { return rep[i]; }
125
126
swap(zz_pX & x)127 void swap(zz_pX& x)
128 {
129 rep.swap(x.rep);
130 }
131
132
133 static const zz_pX& zero();
134
zz_pX(zz_pX & x,INIT_TRANS_TYPE)135 zz_pX(zz_pX& x, INIT_TRANS_TYPE) : rep(x.rep, INIT_TRANS) { }
136
137 };
138
139
140 NTL_DECLARE_RELOCATABLE((zz_pX*))
141
142
143
144
145 /********************************************************************
146
147 input and output
148
149 I/O format:
150
151 [a_0 a_1 ... a_n],
152
153 represents the polynomial a_0 + a_1*X + ... + a_n*X^n.
154
155 On output, all coefficients will be integers between 0 and p-1,
156 amd a_n not zero (the zero polynomial is [ ]).
157 On input, the coefficients are arbitrary integers which are
158 then reduced modulo p, and leading zeros stripped.
159
160 *********************************************************************/
161
162
163 NTL_SNS istream& operator>>(NTL_SNS istream& s, zz_pX& x);
164 NTL_SNS ostream& operator<<(NTL_SNS ostream& s, const zz_pX& a);
165
166
167
168
169 /**********************************************************
170
171 Some utility routines
172
173 ***********************************************************/
174
175
deg(const zz_pX & a)176 inline long deg(const zz_pX& a) { return a.rep.length() - 1; }
177 // degree of a polynomial.
178 // note that the zero polynomial has degree -1.
179
180 const zz_p coeff(const zz_pX& a, long i);
181 // zero if i not in range
182
183 void GetCoeff(zz_p& x, const zz_pX& a, long i);
184 // x = a[i], or zero if i not in range
185
186 const zz_p LeadCoeff(const zz_pX& a);
187 // zero if a == 0
188
189 const zz_p ConstTerm(const zz_pX& a);
190 // zero if a == 0
191
192 void SetCoeff(zz_pX& x, long i, zz_p a);
193 // x[i] = a, error is raised if i < 0
194
195 void SetCoeff(zz_pX& x, long i, long a);
196 // x[i] = a, error is raised if i < 0
197
198 void SetCoeff(zz_pX& x, long i);
199 // x[i] = 1, error is raised if i < 0
200
zz_pX(long i,zz_p a)201 inline zz_pX::zz_pX(long i, zz_p a) { SetCoeff(*this, i, a); }
zz_pX(long i,long a)202 inline zz_pX::zz_pX(long i, long a) { SetCoeff(*this, i, a); }
203
zz_pX(INIT_MONO_TYPE,long i,zz_p a)204 inline zz_pX::zz_pX(INIT_MONO_TYPE, long i, zz_p a) { SetCoeff(*this, i, a); }
zz_pX(INIT_MONO_TYPE,long i,long a)205 inline zz_pX::zz_pX(INIT_MONO_TYPE, long i, long a) { SetCoeff(*this, i, a); }
zz_pX(INIT_MONO_TYPE,long i)206 inline zz_pX::zz_pX(INIT_MONO_TYPE, long i) { SetCoeff(*this, i); }
207
208 void SetX(zz_pX& x);
209 // x is set to the monomial X
210
211 long IsX(const zz_pX& a);
212 // test if x = X
213
clear(zz_pX & x)214 inline void clear(zz_pX& x)
215 // x = 0
216
217 { x.rep.SetLength(0); }
218
set(zz_pX & x)219 inline void set(zz_pX& x)
220 // x = 1
221
222 { x.rep.SetLength(1); set(x.rep[0]); }
223
swap(zz_pX & x,zz_pX & y)224 inline void swap(zz_pX& x, zz_pX& y)
225 // swap x & y (only pointers are swapped)
226
227 { x.swap(y); }
228
229 void random(zz_pX& x, long n);
random_zz_pX(long n)230 inline zz_pX random_zz_pX(long n)
231 { zz_pX x; random(x, n); NTL_OPT_RETURN(zz_pX, x); }
232
233 // generate a random polynomial of degree < n
234
235 void trunc(zz_pX& x, const zz_pX& a, long m);
236 // x = a % X^m
237
trunc(const zz_pX & a,long m)238 inline zz_pX trunc(const zz_pX& a, long m)
239 { zz_pX x; trunc(x, a, m); NTL_OPT_RETURN(zz_pX, x); }
240
241 void RightShift(zz_pX& x, const zz_pX& a, long n);
242 // x = a/X^n
243
RightShift(const zz_pX & a,long n)244 inline zz_pX RightShift(const zz_pX& a, long n)
245 { zz_pX x; RightShift(x, a, n); NTL_OPT_RETURN(zz_pX, x); }
246
247 void LeftShift(zz_pX& x, const zz_pX& a, long n);
248 // x = a*X^n
249
LeftShift(const zz_pX & a,long n)250 inline zz_pX LeftShift(const zz_pX& a, long n)
251 { zz_pX x; LeftShift(x, a, n); NTL_OPT_RETURN(zz_pX, x); }
252
253
254 #ifndef NTL_TRANSITION
255
256 inline zz_pX operator>>(const zz_pX& a, long n)
257 { zz_pX x; RightShift(x, a, n); NTL_OPT_RETURN(zz_pX, x); }
258
259 inline zz_pX operator<<(const zz_pX& a, long n)
260 { zz_pX x; LeftShift(x, a, n); NTL_OPT_RETURN(zz_pX, x); }
261
262 inline zz_pX& operator<<=(zz_pX& x, long n)
263 { LeftShift(x, x, n); return x; }
264
265 inline zz_pX& operator>>=(zz_pX& x, long n)
266 { RightShift(x, x, n); return x; }
267
268 #endif
269
270
271
272 void diff(zz_pX& x, const zz_pX& a);
273 // x = derivative of a
274
diff(const zz_pX & a)275 inline zz_pX diff(const zz_pX& a)
276 { zz_pX x; diff(x, a); NTL_OPT_RETURN(zz_pX, x); }
277
278 void MakeMonic(zz_pX& x);
279 // makes x monic
280
281 void reverse(zz_pX& c, const zz_pX& a, long hi);
282
reverse(const zz_pX & a,long hi)283 inline zz_pX reverse(const zz_pX& a, long hi)
284 { zz_pX x; reverse(x, a, hi); NTL_OPT_RETURN(zz_pX, x); }
285
reverse(zz_pX & c,const zz_pX & a)286 inline void reverse(zz_pX& c, const zz_pX& a)
287 { reverse(c, a, deg(a)); }
288
reverse(const zz_pX & a)289 inline zz_pX reverse(const zz_pX& a)
290 { zz_pX x; reverse(x, a); NTL_OPT_RETURN(zz_pX, x); }
291
292
VectorCopy(vec_zz_p & x,const zz_pX & a,long n)293 inline void VectorCopy(vec_zz_p& x, const zz_pX& a, long n)
294 { VectorCopy(x, a.rep, n); }
295
VectorCopy(const zz_pX & a,long n)296 inline vec_zz_p VectorCopy(const zz_pX& a, long n)
297 { return VectorCopy(a.rep, n); }
298
299
300
301
302 /*******************************************************************
303
304 conversion routines
305
306 ********************************************************************/
307
308
309
310 void conv(zz_pX& x, long a);
311
to_zz_pX(long a)312 inline zz_pX to_zz_pX(long a)
313 { zz_pX x; conv(x, a); NTL_OPT_RETURN(zz_pX, x); }
314
315
316 void conv(zz_pX& x, const ZZ& a);
317
to_zz_pX(const ZZ & a)318 inline zz_pX to_zz_pX(const ZZ& a)
319 { zz_pX x; conv(x, a); NTL_OPT_RETURN(zz_pX, x); }
320
321 void conv(zz_pX& x, zz_p a);
322
to_zz_pX(zz_p a)323 inline zz_pX to_zz_pX(zz_p a)
324 { zz_pX x; conv(x, a); NTL_OPT_RETURN(zz_pX, x); }
325
326
327 void conv(zz_pX& x, const vec_zz_p& a);
328
to_zz_pX(const vec_zz_p & a)329 inline zz_pX to_zz_pX(const vec_zz_p& a)
330 { zz_pX x; conv(x, a); NTL_OPT_RETURN(zz_pX, x); }
331
332 inline zz_pX& zz_pX::operator=(zz_p a)
333 { conv(*this, a); return *this; }
334
335 inline zz_pX& zz_pX::operator=(long a)
336 { conv(*this, a); return *this; }
337
338
339 /* additional legacy conversions for v6 conversion regime */
340
conv(zz_pX & x,const zz_pX & a)341 inline void conv(zz_pX& x, const zz_pX& a)
342 { x = a; }
343
conv(vec_zz_p & x,const zz_pX & a)344 inline void conv(vec_zz_p& x, const zz_pX& a)
345 { x = a.rep; }
346
347
348 /* ------------------------------------- */
349
350
351
352 /*************************************************************
353
354 Comparison
355
356 **************************************************************/
357
358 long IsZero(const zz_pX& a);
359
360 long IsOne(const zz_pX& a);
361
362 inline long operator==(const zz_pX& a, const zz_pX& b)
363 {
364 return a.rep == b.rep;
365 }
366
367 inline long operator!=(const zz_pX& a, const zz_pX& b)
368 { return !(a == b); }
369
370 long operator==(const zz_pX& a, long b);
371 long operator==(const zz_pX& a, zz_p b);
372
373 inline long operator==(long a, const zz_pX& b) { return b == a; }
374 inline long operator==(zz_p a, const zz_pX& b) { return b == a; }
375
376 inline long operator!=(const zz_pX& a, long b) { return !(a == b); }
377 inline long operator!=(const zz_pX& a, zz_p b) { return !(a == b); }
378 inline long operator!=(long a, const zz_pX& b) { return !(a == b); }
379 inline long operator!=(zz_p a, const zz_pX& b) { return !(a == b); }
380
381
382
383 /***************************************************************
384
385 Addition
386
387 ****************************************************************/
388
389 void add(zz_pX& x, const zz_pX& a, const zz_pX& b);
390 // x = a + b
391
392 void sub(zz_pX& x, const zz_pX& a, const zz_pX& b);
393 // x = a - b
394
395 void negate(zz_pX& x, const zz_pX& a);
396 // x = -a
397
398 // scalar versions
399
400 void add(zz_pX & x, const zz_pX& a, zz_p b); // x = a + b
add(zz_pX & x,const zz_pX & a,long b)401 inline void add(zz_pX& x, const zz_pX& a, long b) { add(x, a, to_zz_p(b)); }
402
add(zz_pX & x,zz_p a,const zz_pX & b)403 inline void add(zz_pX& x, zz_p a, const zz_pX& b) { add(x, b, a); }
add(zz_pX & x,long a,const zz_pX & b)404 inline void add(zz_pX& x, long a, const zz_pX& b) { add(x, b, a); }
405
406 void sub(zz_pX & x, const zz_pX& a, zz_p b); // x = a - b
sub(zz_pX & x,const zz_pX & a,long b)407 inline void sub(zz_pX& x, const zz_pX& a, long b) { sub(x, a, to_zz_p(b)); }
408
409 void sub(zz_pX& x, zz_p a, const zz_pX& b);
sub(zz_pX & x,long a,const zz_pX & b)410 inline void sub(zz_pX& x, long a, const zz_pX& b) { sub(x, to_zz_p(a), b); }
411
412 inline zz_pX operator+(const zz_pX& a, const zz_pX& b)
413 { zz_pX x; add(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
414
415 inline zz_pX operator+(const zz_pX& a, zz_p b)
416 { zz_pX x; add(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
417
418 inline zz_pX operator+(const zz_pX& a, long b)
419 { zz_pX x; add(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
420
421 inline zz_pX operator+(zz_p a, const zz_pX& b)
422 { zz_pX x; add(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
423
424 inline zz_pX operator+(long a, const zz_pX& b)
425 { zz_pX x; add(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
426
427
428 inline zz_pX operator-(const zz_pX& a, const zz_pX& b)
429 { zz_pX x; sub(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
430
431 inline zz_pX operator-(const zz_pX& a, zz_p b)
432 { zz_pX x; sub(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
433
434 inline zz_pX operator-(const zz_pX& a, long b)
435 { zz_pX x; sub(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
436
437 inline zz_pX operator-(zz_p a, const zz_pX& b)
438 { zz_pX x; sub(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
439
440 inline zz_pX operator-(long a, const zz_pX& b)
441 { zz_pX x; sub(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
442
443
444 inline zz_pX& operator+=(zz_pX& x, const zz_pX& b)
445 { add(x, x, b); return x; }
446
447 inline zz_pX& operator+=(zz_pX& x, zz_p b)
448 { add(x, x, b); return x; }
449
450 inline zz_pX& operator+=(zz_pX& x, long b)
451 { add(x, x, b); return x; }
452
453 inline zz_pX& operator-=(zz_pX& x, const zz_pX& b)
454 { sub(x, x, b); return x; }
455
456 inline zz_pX& operator-=(zz_pX& x, zz_p b)
457 { sub(x, x, b); return x; }
458
459 inline zz_pX& operator-=(zz_pX& x, long b)
460 { sub(x, x, b); return x; }
461
462
463 inline zz_pX operator-(const zz_pX& a)
464 { zz_pX x; negate(x, a); NTL_OPT_RETURN(zz_pX, x); }
465
466 inline zz_pX& operator++(zz_pX& x) { add(x, x, 1); return x; }
467 inline void operator++(zz_pX& x, int) { add(x, x, 1); }
468 inline zz_pX& operator--(zz_pX& x) { sub(x, x, 1); return x; }
469 inline void operator--(zz_pX& x, int) { sub(x, x, 1); }
470
471
472
473 /*****************************************************************
474
475 Multiplication
476
477 ******************************************************************/
478
479
480 void mul(zz_pX& x, const zz_pX& a, const zz_pX& b);
481 // x = a * b
482
483 void sqr(zz_pX& x, const zz_pX& a);
sqr(const zz_pX & a)484 inline zz_pX sqr(const zz_pX& a)
485 { zz_pX x; sqr(x, a); NTL_OPT_RETURN(zz_pX, x); }
486 // x = a^2
487
488 void mul(zz_pX& x, const zz_pX& a, zz_p b);
mul(zz_pX & x,const zz_pX & a,long b)489 inline void mul(zz_pX& x, const zz_pX& a, long b) { mul(x, a, to_zz_p(b)); }
490
mul(zz_pX & x,zz_p a,const zz_pX & b)491 inline void mul(zz_pX& x, zz_p a, const zz_pX& b) { mul(x, b, a); }
mul(zz_pX & x,long a,const zz_pX & b)492 inline void mul(zz_pX& x, long a, const zz_pX& b) { mul(x, b, a); }
493
494
495 inline zz_pX operator*(const zz_pX& a, const zz_pX& b)
496 { zz_pX x; mul(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
497
498 inline zz_pX operator*(const zz_pX& a, zz_p b)
499 { zz_pX x; mul(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
500
501 inline zz_pX operator*(const zz_pX& a, long b)
502 { zz_pX x; mul(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
503
504 inline zz_pX operator*(zz_p a, const zz_pX& b)
505 { zz_pX x; mul(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
506
507 inline zz_pX operator*(long a, const zz_pX& b)
508 { zz_pX x; mul(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
509
510 inline zz_pX& operator*=(zz_pX& x, const zz_pX& b)
511 { mul(x, x, b); return x; }
512
513 inline zz_pX& operator*=(zz_pX& x, zz_p b)
514 { mul(x, x, b); return x; }
515
516 inline zz_pX& operator*=(zz_pX& x, long b)
517 { mul(x, x, b); return x; }
518
519
520 void PlainMul(zz_pX& x, const zz_pX& a, const zz_pX& b);
521 // always uses the "classical" algorithm
522
523 void PlainSqr(zz_pX& x, const zz_pX& a);
524 // always uses the "classical" algorithm
525
526
527 void FFTMul(zz_pX& x, const zz_pX& a, const zz_pX& b);
528 // always uses the FFT
529
530 void FFTSqr(zz_pX& x, const zz_pX& a);
531 // always uses the FFT
532
533 void MulTrunc(zz_pX& x, const zz_pX& a, const zz_pX& b, long n);
534 // x = a * b % X^n
535
MulTrunc(const zz_pX & a,const zz_pX & b,long n)536 inline zz_pX MulTrunc(const zz_pX& a, const zz_pX& b, long n)
537 { zz_pX x; MulTrunc(x, a, b, n); NTL_OPT_RETURN(zz_pX, x); }
538
539 void PlainMulTrunc(zz_pX& x, const zz_pX& a, const zz_pX& b, long n);
540 void FFTMulTrunc(zz_pX& x, const zz_pX& a, const zz_pX& b, long n);
541
542 void SqrTrunc(zz_pX& x, const zz_pX& a, long n);
543 // x = a^2 % X^n
544
SqrTrunc(const zz_pX & a,long n)545 inline zz_pX SqrTrunc(const zz_pX& a, long n)
546 { zz_pX x; SqrTrunc(x, a, n); NTL_OPT_RETURN(zz_pX, x); }
547
548 void PlainSqrTrunc(zz_pX& x, const zz_pX& a, long n);
549 void FFTSqrTrunc(zz_pX& x, const zz_pX& a, long n);
550
551 void power(zz_pX& x, const zz_pX& a, long e);
power(const zz_pX & a,long e)552 inline zz_pX power(const zz_pX& a, long e)
553 { zz_pX x; power(x, a, e); NTL_OPT_RETURN(zz_pX, x); }
554
555
556
557
558
559
560 // The following data structures and routines allow one
561 // to hand-craft various algorithms, using the FFT convolution
562 // algorithms directly.
563 // Look in the file zz_pX.c for examples.
564
565
566
567
568 // FFT representation of polynomials
569
570 class fftRep {
571
572 public:
573 long k; // a 2^k point representation
574 long MaxK; // maximum space allocated
575 long len; // length of truncated FFT
576 long NumPrimes;
577 UniqueArray<long> tbl[4];
578
fftRep()579 fftRep() : k(-1), MaxK(-1), len(0), NumPrimes(0) { }
580
fftRep(const fftRep & R)581 fftRep(const fftRep& R) : k(-1), MaxK(-1), len(0), NumPrimes(0)
582 { *this = R; }
583
fftRep(INIT_SIZE_TYPE,long InitK)584 fftRep(INIT_SIZE_TYPE, long InitK) : k(-1), MaxK(-1), len(0), NumPrimes(0)
585 { SetSize(InitK); }
586
587 fftRep& operator=(const fftRep&);
588 void SetSize(long NewK);
589 void DoSetSize(long NewK, long NewNumPrimes);
590 };
591
592
593
594
595 void TofftRep_trunc(fftRep& y, const zz_pX& x, long k, long len,
596 long lo, long hi);
597
TofftRep_trunc(fftRep & y,const zz_pX & x,long k,long len)598 inline void TofftRep_trunc(fftRep& y, const zz_pX& x, long k, long len)
599 { TofftRep_trunc(y, x, k, len, 0, deg(x)); }
600
601 inline
TofftRep(fftRep & y,const zz_pX & x,long k,long lo,long hi)602 void TofftRep(fftRep& y, const zz_pX& x, long k, long lo, long hi)
603 // computes an n = 2^k point convolution of x[lo..hi].
604 { TofftRep_trunc(y, x, k, 1L << k, lo, hi); }
605
TofftRep(fftRep & y,const zz_pX & x,long k)606 inline void TofftRep(fftRep& y, const zz_pX& x, long k)
607
608 { TofftRep(y, x, k, 0, deg(x)); }
609
610 void RevTofftRep(fftRep& y, const vec_zz_p& x,
611 long k, long lo, long hi, long offset);
612 // computes an n = 2^k point convolution of X^offset*x[lo..hi] mod X^n-1
613 // using "inverted" evaluation points.
614
615
616
617 void FromfftRep(zz_pX& x, fftRep& y, long lo, long hi);
618 // converts from FFT-representation to coefficient representation
619 // only the coefficients lo..hi are computed
620 // NOTE: this version destroys the data in y
621
622 // non-destructive versions of the above
623
624 void NDFromfftRep(zz_pX& x, const fftRep& y, long lo, long hi, fftRep& temp);
625 void NDFromfftRep(zz_pX& x, const fftRep& y, long lo, long hi);
626
627 void RevFromfftRep(vec_zz_p& x, fftRep& y, long lo, long hi);
628
629 // converts from FFT-representation to coefficient representation
630 // using "inverted" evaluation points.
631 // only the coefficients lo..hi are computed
632
633
634
635
636 void FromfftRep(zz_p* x, fftRep& y, long lo, long hi);
637 // convert out coefficients lo..hi of y, store result in x.
638 // no normalization is done.
639
640
641 // direct manipulation of FFT reps
642
643 void mul(fftRep& z, const fftRep& x, const fftRep& y);
644 void sub(fftRep& z, const fftRep& x, const fftRep& y);
645 void add(fftRep& z, const fftRep& x, const fftRep& y);
646
647 void reduce(fftRep& x, const fftRep& a, long k);
648 // reduces a 2^l point FFT-rep to a 2^k point FFT-rep
649
650 void AddExpand(fftRep& x, const fftRep& a);
651 // x = x + (an "expanded" version of a)
652
653
654
655
656
657
658
659 /*************************************************************
660
661 Division
662
663 **************************************************************/
664
665 void DivRem(zz_pX& q, zz_pX& r, const zz_pX& a, const zz_pX& b);
666 // q = a/b, r = a%b
667
668 void div(zz_pX& q, const zz_pX& a, const zz_pX& b);
669 // q = a/b
670
671
672 void div(zz_pX& q, const zz_pX& a, zz_p b);
div(zz_pX & q,const zz_pX & a,long b)673 inline void div(zz_pX& q, const zz_pX& a, long b)
674 { div(q, a, to_zz_p(b)); }
675
676 void rem(zz_pX& r, const zz_pX& a, const zz_pX& b);
677 // r = a%b
678
679 long divide(zz_pX& q, const zz_pX& a, const zz_pX& b);
680 // if b | a, sets q = a/b and returns 1; otherwise returns 0
681
682 long divide(const zz_pX& a, const zz_pX& b);
683 // if b | a, sets q = a/b and returns 1; otherwise returns 0
684
685
686 void InvTrunc(zz_pX& x, const zz_pX& a, long m);
687 // computes x = a^{-1} % X^m
688 // constant term must be non-zero
689
InvTrunc(const zz_pX & a,long m)690 inline zz_pX InvTrunc(const zz_pX& a, long m)
691 { zz_pX x; InvTrunc(x, a, m); NTL_OPT_RETURN(zz_pX, x); }
692
693
694
695
696 // These always use "classical" arithmetic
697 void PlainDivRem(zz_pX& q, zz_pX& r, const zz_pX& a, const zz_pX& b);
698 void PlainDiv(zz_pX& q, const zz_pX& a, const zz_pX& b);
699 void PlainRem(zz_pX& r, const zz_pX& a, const zz_pX& b);
700
701
702 // These always use FFT arithmetic
703 void FFTDivRem(zz_pX& q, zz_pX& r, const zz_pX& a, const zz_pX& b);
704 void FFTDiv(zz_pX& q, const zz_pX& a, const zz_pX& b);
705 void FFTRem(zz_pX& r, const zz_pX& a, const zz_pX& b);
706
707 void PlainInvTrunc(zz_pX& x, const zz_pX& a, long m);
708 // always uses "classical" algorithm
709 // ALIAS RESTRICTION: input may not alias output
710
711 void NewtonInvTrunc(zz_pX& x, const zz_pX& a, long m);
712 // uses a Newton Iteration with the FFT.
713 // ALIAS RESTRICTION: input may not alias output
714
715
716 inline zz_pX operator/(const zz_pX& a, const zz_pX& b)
717 { zz_pX x; div(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
718
719 inline zz_pX operator/(const zz_pX& a, zz_p b)
720 { zz_pX x; div(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
721
722 inline zz_pX operator/(const zz_pX& a, long b)
723 { zz_pX x; div(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
724
725 inline zz_pX& operator/=(zz_pX& x, zz_p b)
726 { div(x, x, b); return x; }
727
728 inline zz_pX& operator/=(zz_pX& x, long b)
729 { div(x, x, b); return x; }
730
731 inline zz_pX& operator/=(zz_pX& x, const zz_pX& b)
732 { div(x, x, b); return x; }
733
734
735 inline zz_pX operator%(const zz_pX& a, const zz_pX& b)
736 { zz_pX x; rem(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
737
738 inline zz_pX& operator%=(zz_pX& x, const zz_pX& b)
739 { rem(x, x, b); return x; }
740
741
742
743
744 /***********************************************************
745
746 GCD's
747
748 ************************************************************/
749
750
751 void GCD(zz_pX& x, const zz_pX& a, const zz_pX& b);
752 // x = GCD(a, b), x is always monic (or zero if a==b==0).
753
GCD(const zz_pX & a,const zz_pX & b)754 inline zz_pX GCD(const zz_pX& a, const zz_pX& b)
755 { zz_pX x; GCD(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
756
757
758 void XGCD(zz_pX& d, zz_pX& s, zz_pX& t, const zz_pX& a, const zz_pX& b);
759 // d = gcd(a,b), a s + b t = d
760
761 void PlainXGCD(zz_pX& d, zz_pX& s, zz_pX& t, const zz_pX& a, const zz_pX& b);
762 // same as above, but uses classical algorithm
763
764
765 void PlainGCD(zz_pX& x, const zz_pX& a, const zz_pX& b);
766 // always uses "cdlassical" arithmetic
767
768
769 class zz_pXMatrix {
770 private:
771
772 zz_pXMatrix(const zz_pXMatrix&); // disable
773 zz_pX elts[2][2];
774
775 public:
776
zz_pXMatrix()777 zz_pXMatrix() { }
778
779 void operator=(const zz_pXMatrix&);
operator()780 zz_pX& operator() (long i, long j) { return elts[i][j]; }
operator()781 const zz_pX& operator() (long i, long j) const { return elts[i][j]; }
782 };
783
784
785 void HalfGCD(zz_pXMatrix& M_out, const zz_pX& U, const zz_pX& V, long d_red);
786 // deg(U) > deg(V), 1 <= d_red <= deg(U)+1.
787 //
788 // This computes a 2 x 2 polynomial matrix M_out such that
789 // M_out * (U, V)^T = (U', V')^T,
790 // where U', V' are consecutive polynomials in the Euclidean remainder
791 // sequence of U, V, and V' is the polynomial of highest degree
792 // satisfying deg(V') <= deg(U) - d_red.
793
794 void XHalfGCD(zz_pXMatrix& M_out, zz_pX& U, zz_pX& V, long d_red);
795
796 // same as above, except that U is replaced by U', and V by V'
797
798
799 /*************************************************************
800
801 Modular Arithmetic without pre-conditioning
802
803 **************************************************************/
804
805 // arithmetic mod f.
806 // all inputs and outputs are polynomials of degree less than deg(f).
807 // ASSUMPTION: f is assumed monic, and deg(f) > 0.
808 // NOTE: if you want to do many computations with a fixed f,
809 // use the zz_pXModulus data structure and associated routines below.
810
811
812
813 void MulMod(zz_pX& x, const zz_pX& a, const zz_pX& b, const zz_pX& f);
814 // x = (a * b) % f
815
MulMod(const zz_pX & a,const zz_pX & b,const zz_pX & f)816 inline zz_pX MulMod(const zz_pX& a, const zz_pX& b, const zz_pX& f)
817 { zz_pX x; MulMod(x, a, b, f); NTL_OPT_RETURN(zz_pX, x); }
818
819 void SqrMod(zz_pX& x, const zz_pX& a, const zz_pX& f);
820 // x = a^2 % f
821
SqrMod(const zz_pX & a,const zz_pX & f)822 inline zz_pX SqrMod(const zz_pX& a, const zz_pX& f)
823 { zz_pX x; SqrMod(x, a, f); NTL_OPT_RETURN(zz_pX, x); }
824
825 void MulByXMod(zz_pX& x, const zz_pX& a, const zz_pX& f);
826 // x = (a * X) mod f
827
MulByXMod(const zz_pX & a,const zz_pX & f)828 inline zz_pX MulByXMod(const zz_pX& a, const zz_pX& f)
829 { zz_pX x; MulByXMod(x, a, f); NTL_OPT_RETURN(zz_pX, x); }
830
831 void InvMod(zz_pX& x, const zz_pX& a, const zz_pX& f);
832 // x = a^{-1} % f, error is a is not invertible
833
InvMod(const zz_pX & a,const zz_pX & f)834 inline zz_pX InvMod(const zz_pX& a, const zz_pX& f)
835 { zz_pX x; InvMod(x, a, f); NTL_OPT_RETURN(zz_pX, x); }
836
837 long InvModStatus(zz_pX& x, const zz_pX& a, const zz_pX& f);
838 // if (a, f) = 1, returns 0 and sets x = a^{-1} % f
839 // otherwise, returns 1 and sets x = (a, f)
840
841
842
843 /******************************************************************
844
845 Modular Arithmetic with Pre-conditioning
846
847 *******************************************************************/
848
849
850 // If you need to do a lot of arithmetic modulo a fixed f,
851 // build zz_pXModulus F for f. This pre-computes information about f
852 // that speeds up the computation a great deal.
853
854
855 class zz_pXModulus {
856 public:
zz_pXModulus()857 zz_pXModulus() : UseFFT(0), n(-1) { }
858
859 zz_pX f; // the modulus
860 long UseFFT;// flag indicating whether FFT should be used.
861 long n; // n = deg(f)
862 long k; // least k s/t 2^k >= n
863 long l; // least l s/t 2^l >= 2n-3
864 fftRep FRep; // 2^k point rep of f
865 // H = rev((rev(f))^{-1} rem X^{n-1})
866 fftRep HRep; // 2^l point rep of H
867
868 OptionalVal< Lazy<vec_zz_p> > tracevec;
869 // extra level of indirection to ensure relocatability
870
871 zz_pXModulus(const zz_pX& ff);
872
873 operator const zz_pX& () const { return f; }
val()874 const zz_pX& val() const { return f; }
875
876 };
877
878
879 NTL_DECLARE_RELOCATABLE((zz_pXModulus*))
880
881
deg(const zz_pXModulus & F)882 inline long deg(const zz_pXModulus& F) { return F.n; }
883
884 void build(zz_pXModulus& F, const zz_pX& f);
885 // deg(f) > 0
886
887
888 void rem21(zz_pX& x, const zz_pX& a, const zz_pXModulus& F);
889 // x = a % f
890 // deg(a) <= 2(n-1), where n = F.n = deg(f)
891
892 void rem(zz_pX& x, const zz_pX& a, const zz_pXModulus& F);
893 // x = a % f, no restrictions on deg(a); makes repeated calls to rem21
894
895 inline zz_pX operator%(const zz_pX& a, const zz_pXModulus& F)
896 { zz_pX x; rem(x, a, F); NTL_OPT_RETURN(zz_pX, x); }
897
898 inline zz_pX& operator%=(zz_pX& x, const zz_pXModulus& F)
899 { rem(x, x, F); return x; }
900
901
902 void DivRem(zz_pX& q, zz_pX& r, const zz_pX& a, const zz_pXModulus& F);
903
904 void div(zz_pX& q, const zz_pX& a, const zz_pXModulus& F);
905
906 inline zz_pX operator/(const zz_pX& a, const zz_pXModulus& F)
907 { zz_pX x; div(x, a, F); NTL_OPT_RETURN(zz_pX, x); }
908
909 inline zz_pX& operator/=(zz_pX& x, const zz_pXModulus& F)
910 { div(x, x, F); return x; }
911
912
913
914 void MulMod(zz_pX& x, const zz_pX& a, const zz_pX& b, const zz_pXModulus& F);
915 // x = (a * b) % f
916 // deg(a), deg(b) < n
917
MulMod(const zz_pX & a,const zz_pX & b,const zz_pXModulus & F)918 inline zz_pX MulMod(const zz_pX& a, const zz_pX& b, const zz_pXModulus& F)
919 { zz_pX x; MulMod(x, a, b, F); NTL_OPT_RETURN(zz_pX, x); }
920
921
922 void SqrMod(zz_pX& x, const zz_pX& a, const zz_pXModulus& F);
923 // x = a^2 % f
924 // deg(a) < n
925
926
SqrMod(const zz_pX & a,const zz_pXModulus & F)927 inline zz_pX SqrMod(const zz_pX& a, const zz_pXModulus& F)
928 { zz_pX x; SqrMod(x, a, F); NTL_OPT_RETURN(zz_pX, x); }
929
930
931 void PowerMod(zz_pX& x, const zz_pX& a, const ZZ& e, const zz_pXModulus& F);
932 // x = a^e % f, e >= 0
933
PowerMod(const zz_pX & a,const ZZ & e,const zz_pXModulus & F)934 inline zz_pX PowerMod(const zz_pX& a, const ZZ& e, const zz_pXModulus& F)
935 { zz_pX x; PowerMod(x, a, e, F); NTL_OPT_RETURN(zz_pX, x); }
936
PowerMod(zz_pX & x,const zz_pX & a,long e,const zz_pXModulus & F)937 inline void PowerMod(zz_pX& x, const zz_pX& a, long e, const zz_pXModulus& F)
938 { PowerMod(x, a, ZZ_expo(e), F); }
939
PowerMod(const zz_pX & a,long e,const zz_pXModulus & F)940 inline zz_pX PowerMod(const zz_pX& a, long e, const zz_pXModulus& F)
941 { zz_pX x; PowerMod(x, a, e, F); NTL_OPT_RETURN(zz_pX, x); }
942
943
944
945 void PowerXMod(zz_pX& x, const ZZ& e, const zz_pXModulus& F);
946 // x = X^e % f, e >= 0
947
PowerXMod(const ZZ & e,const zz_pXModulus & F)948 inline zz_pX PowerXMod(const ZZ& e, const zz_pXModulus& F)
949 { zz_pX x; PowerXMod(x, e, F); NTL_OPT_RETURN(zz_pX, x); }
950
PowerXMod(zz_pX & x,long e,const zz_pXModulus & F)951 inline void PowerXMod(zz_pX& x, long e, const zz_pXModulus& F)
952 { PowerXMod(x, ZZ_expo(e), F); }
953
PowerXMod(long e,const zz_pXModulus & F)954 inline zz_pX PowerXMod(long e, const zz_pXModulus& F)
955 { zz_pX x; PowerXMod(x, e, F); NTL_OPT_RETURN(zz_pX, x); }
956
957 void PowerXPlusAMod(zz_pX& x, zz_p a, const ZZ& e, const zz_pXModulus& F);
958 // x = (X + a)^e % f, e >= 0
959
PowerXPlusAMod(zz_p a,const ZZ & e,const zz_pXModulus & F)960 inline zz_pX PowerXPlusAMod(zz_p a, const ZZ& e, const zz_pXModulus& F)
961 { zz_pX x; PowerXPlusAMod(x, a, e, F); NTL_OPT_RETURN(zz_pX, x); }
962
PowerXPlusAMod(zz_pX & x,zz_p a,long e,const zz_pXModulus & F)963 inline void PowerXPlusAMod(zz_pX& x, zz_p a, long e, const zz_pXModulus& F)
964 { PowerXPlusAMod(x, a, ZZ_expo(e), F); }
965
966
PowerXPlusAMod(zz_p a,long e,const zz_pXModulus & F)967 inline zz_pX PowerXPlusAMod(zz_p a, long e, const zz_pXModulus& F)
968 { zz_pX x; PowerXPlusAMod(x, a, e, F); NTL_OPT_RETURN(zz_pX, x); }
969
970 // If you need to compute a * b % f for a fixed b, but for many a's
971 // (for example, computing powers of b modulo f), it is
972 // much more efficient to first build a zz_pXMultiplier B for b,
973 // and then use the routine below.
974
975 class zz_pXMultiplier {
976 public:
zz_pXMultiplier()977 zz_pXMultiplier() : UseFFT(0) { }
978 zz_pXMultiplier(const zz_pX& b, const zz_pXModulus& F);
979
980
981 zz_pX b;
982 long UseFFT;
983 fftRep B1;
984 fftRep B2;
985
val()986 const zz_pX& val() const { return b; }
987
988 };
989
990 void build(zz_pXMultiplier& B, const zz_pX& b, const zz_pXModulus& F);
991
992
993
994 void MulMod(zz_pX& x, const zz_pX& a, const zz_pXMultiplier& B,
995 const zz_pXModulus& F);
996
997 // x = (a * b) % f
998
MulMod(const zz_pX & a,const zz_pXMultiplier & B,const zz_pXModulus & F)999 inline zz_pX MulMod(const zz_pX& a, const zz_pXMultiplier& B,
1000 const zz_pXModulus& F)
1001 { zz_pX x; MulMod(x, a, B, F); NTL_OPT_RETURN(zz_pX, x); }
1002
1003
1004
1005
1006 /*******************************************************
1007
1008 Evaluation and related problems
1009
1010 ********************************************************/
1011
1012
1013 void BuildFromRoots(zz_pX& x, const vec_zz_p& a);
1014 // computes the polynomial (X-a[0]) ... (X-a[n-1]), where n = a.length()
1015
BuildFromRoots(const vec_zz_p & a)1016 inline zz_pX BuildFromRoots(const vec_zz_p& a)
1017 { zz_pX x; BuildFromRoots(x, a); NTL_OPT_RETURN(zz_pX, x); }
1018
1019
1020
1021 void eval(zz_p& b, const zz_pX& f, zz_p a);
1022 // b = f(a)
1023
1024
eval(const zz_pX & f,zz_p a)1025 inline zz_p eval(const zz_pX& f, zz_p a)
1026 { zz_p x; eval(x, f, a); return x; }
1027
1028
1029 void eval(vec_zz_p& b, const zz_pX& f, const vec_zz_p& a);
1030 // b[i] = f(a[i])
1031
eval(const zz_pX & f,const vec_zz_p & a)1032 inline vec_zz_p eval(const zz_pX& f, const vec_zz_p& a)
1033 { vec_zz_p x; eval(x, f, a); NTL_OPT_RETURN(vec_zz_p, x); }
1034
1035
1036 void interpolate(zz_pX& f, const vec_zz_p& a, const vec_zz_p& b);
1037 // computes f such that f(a[i]) = b[i]
1038
interpolate(const vec_zz_p & a,const vec_zz_p & b)1039 inline zz_pX interpolate(const vec_zz_p& a, const vec_zz_p& b)
1040 { zz_pX x; interpolate(x, a, b); NTL_OPT_RETURN(zz_pX, x); }
1041
1042
1043
1044 /*****************************************************************
1045
1046 vectors of zz_pX's
1047
1048 *****************************************************************/
1049
1050 typedef Vec<zz_pX> vec_zz_pX;
1051
1052
1053 /**********************************************************
1054
1055 Modular Composition and Minimal Polynomials
1056
1057 ***********************************************************/
1058
1059
1060 // algorithms for computing g(h) mod f
1061
1062
1063
1064
1065 void CompMod(zz_pX& x, const zz_pX& g, const zz_pX& h, const zz_pXModulus& F);
1066 // x = g(h) mod f
1067
CompMod(const zz_pX & g,const zz_pX & h,const zz_pXModulus & F)1068 inline zz_pX CompMod(const zz_pX& g, const zz_pX& h,
1069 const zz_pXModulus& F)
1070 { zz_pX x; CompMod(x, g, h, F); NTL_OPT_RETURN(zz_pX, x); }
1071
1072
1073 void Comp2Mod(zz_pX& x1, zz_pX& x2, const zz_pX& g1, const zz_pX& g2,
1074 const zz_pX& h, const zz_pXModulus& F);
1075 // xi = gi(h) mod f (i=1,2)
1076
1077 void Comp3Mod(zz_pX& x1, zz_pX& x2, zz_pX& x3,
1078 const zz_pX& g1, const zz_pX& g2, const zz_pX& g3,
1079 const zz_pX& h, const zz_pXModulus& F);
1080 // xi = gi(h) mod f (i=1..3)
1081
1082
1083
1084 // The routine build (see below) which is implicitly called
1085 // by the various compose and UpdateMap routines builds a table
1086 // of polynomials.
1087 // If zz_pXArgBound > 0, then the table is limited in
1088 // size to approximamtely that many KB.
1089 // If zz_pXArgBound <= 0, then it is ignored, and space is allocated
1090 // so as to maximize speed.
1091 // Initially, zz_pXArgBound = 0.
1092
1093
1094 // If a single h is going to be used with many g's
1095 // then you should build a zz_pXArgument for h,
1096 // and then use the compose routine below.
1097 // build computes and stores h, h^2, ..., h^m mod f.
1098 // After this pre-computation, composing a polynomial of degree
1099 // roughly n with h takes n/m multiplies mod f, plus n^2
1100 // scalar multiplies.
1101 // Thus, increasing m increases the space requirement and the pre-computation
1102 // time, but reduces the composition time.
1103 // If zz_pXArgBound > 0, a table of size less than m may be built.
1104
1105 struct zz_pXArgument {
1106 vec_zz_pX H;
1107 };
1108
1109 extern
1110 NTL_CHEAP_THREAD_LOCAL
1111 long zz_pXArgBound;
1112
1113
1114 void build(zz_pXArgument& H, const zz_pX& h, const zz_pXModulus& F, long m);
1115
1116 // m must be > 0, otherwise an error is raised
1117
1118 void CompMod(zz_pX& x, const zz_pX& g, const zz_pXArgument& H,
1119 const zz_pXModulus& F);
1120
1121 inline zz_pX
CompMod(const zz_pX & g,const zz_pXArgument & H,const zz_pXModulus & F)1122 CompMod(const zz_pX& g, const zz_pXArgument& H, const zz_pXModulus& F)
1123 { zz_pX x; CompMod(x, g, H, F); NTL_OPT_RETURN(zz_pX, x); }
1124
1125
1126 // New alternative CompMod strategy that just reduces to
1127 // matrix multiplication ...
1128
1129
1130 struct zz_pXNewArgument {
1131 Mat<zz_p> mat;
1132 zz_pX poly;
1133 };
1134
1135 void build(zz_pXNewArgument& H, const zz_pX& h, const zz_pXModulus& F, long m);
1136 void CompMod(zz_pX& x, const zz_pX& g, const zz_pXNewArgument& H,
1137 const zz_pXModulus& F);
1138 void reduce(zz_pXNewArgument& H, const zz_pXModulus& F);
1139
1140 void ProjectPowers(vec_zz_p& x, const vec_zz_p& a, long k,
1141 const zz_pXNewArgument& H, const zz_pXModulus& F);
1142
1143
1144
1145
1146
1147
1148
1149 #ifndef NTL_TRANSITION
1150
1151 void UpdateMap(vec_zz_p& x, const vec_zz_p& a,
1152 const zz_pXMultiplier& B, const zz_pXModulus& F);
1153
1154 inline vec_zz_p
UpdateMap(const vec_zz_p & a,const zz_pXMultiplier & B,const zz_pXModulus & F)1155 UpdateMap(const vec_zz_p& a,
1156 const zz_pXMultiplier& B, const zz_pXModulus& F)
1157 { vec_zz_p x; UpdateMap(x, a, B, F);
1158 NTL_OPT_RETURN(vec_zz_p, x); }
1159
1160 #endif
1161
1162
1163 /* computes (a, b), (a, (b*X)%f), ..., (a, (b*X^{n-1})%f),
1164 where ( , ) denotes the vector inner product.
1165
1166 This is really a "transposed" MulMod by B.
1167 */
1168
1169 void PlainUpdateMap(vec_zz_p& x, const vec_zz_p& a,
1170 const zz_pX& b, const zz_pX& f);
1171
1172
1173 // same as above, but uses only classical arithmetic
1174
1175
1176 void ProjectPowers(vec_zz_p& x, const vec_zz_p& a, long k,
1177 const zz_pX& h, const zz_pXModulus& F);
1178
1179 // computes (a, 1), (a, h), ..., (a, h^{k-1} % f)
1180 // this is really a "transposed" compose.
1181
ProjectPowers(const vec_zz_p & a,long k,const zz_pX & h,const zz_pXModulus & F)1182 inline vec_zz_p ProjectPowers(const vec_zz_p& a, long k,
1183 const zz_pX& h, const zz_pXModulus& F)
1184 {
1185 vec_zz_p x;
1186 ProjectPowers(x, a, k, h, F);
1187 NTL_OPT_RETURN(vec_zz_p, x);
1188 }
1189
1190
1191 void ProjectPowers(vec_zz_p& x, const vec_zz_p& a, long k,
1192 const zz_pXArgument& H, const zz_pXModulus& F);
1193
ProjectPowers(const vec_zz_p & a,long k,const zz_pXArgument & H,const zz_pXModulus & F)1194 inline vec_zz_p ProjectPowers(const vec_zz_p& a, long k,
1195 const zz_pXArgument& H, const zz_pXModulus& F)
1196 {
1197 vec_zz_p x;
1198 ProjectPowers(x, a, k, H, F);
1199 NTL_OPT_RETURN(vec_zz_p, x);
1200 }
1201
1202 // same as above, but uses a pre-computed zz_pXArgument
1203
project(zz_p & x,const vec_zz_p & a,const zz_pX & b)1204 inline void project(zz_p& x, const vec_zz_p& a, const zz_pX& b)
1205 { InnerProduct(x, a, b.rep); }
1206
project(const vec_zz_p & a,const zz_pX & b)1207 inline zz_p project(const vec_zz_p& a, const zz_pX& b)
1208 { zz_p x; project(x, a, b); return x; }
1209
1210
1211 void MinPolySeq(zz_pX& h, const vec_zz_p& a, long m);
1212 // computes the minimum polynomial of a linealy generated sequence;
1213 // m is a bound on the degree of the polynomial;
1214 // required: a.length() >= 2*m
1215
MinPolySeq(const vec_zz_p & a,long m)1216 inline zz_pX MinPolySeq(const vec_zz_p& a, long m)
1217 { zz_pX x; MinPolySeq(x, a, m); NTL_OPT_RETURN(zz_pX, x); }
1218
1219 void ProbMinPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F, long m);
1220
ProbMinPolyMod(const zz_pX & g,const zz_pXModulus & F,long m)1221 inline zz_pX ProbMinPolyMod(const zz_pX& g, const zz_pXModulus& F, long m)
1222 { zz_pX x; ProbMinPolyMod(x, g, F, m); NTL_OPT_RETURN(zz_pX, x); }
1223
1224
ProbMinPolyMod(zz_pX & h,const zz_pX & g,const zz_pXModulus & F)1225 inline void ProbMinPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F)
1226 { ProbMinPolyMod(h, g, F, F.n); }
1227
ProbMinPolyMod(const zz_pX & g,const zz_pXModulus & F)1228 inline zz_pX ProbMinPolyMod(const zz_pX& g, const zz_pXModulus& F)
1229 { zz_pX x; ProbMinPolyMod(x, g, F); NTL_OPT_RETURN(zz_pX, x); }
1230
1231
1232 // computes the monic minimal polynomial if (g mod f).
1233 // m = a bound on the degree of the minimal polynomial.
1234 // If this argument is not supplied, it defaults to deg(f).
1235 // The algorithm is probabilistic, always returns a divisor of
1236 // the minimal polynomial, and returns a proper divisor with
1237 // probability at most m/p.
1238
1239 void MinPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F, long m);
1240
MinPolyMod(const zz_pX & g,const zz_pXModulus & F,long m)1241 inline zz_pX MinPolyMod(const zz_pX& g, const zz_pXModulus& F, long m)
1242 { zz_pX x; MinPolyMod(x, g, F, m); NTL_OPT_RETURN(zz_pX, x); }
1243
1244
MinPolyMod(zz_pX & h,const zz_pX & g,const zz_pXModulus & F)1245 inline void MinPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F)
1246 { MinPolyMod(h, g, F, F.n); }
1247
MinPolyMod(const zz_pX & g,const zz_pXModulus & F)1248 inline zz_pX MinPolyMod(const zz_pX& g, const zz_pXModulus& F)
1249 { zz_pX x; MinPolyMod(x, g, F); NTL_OPT_RETURN(zz_pX, x); }
1250
1251
1252 // same as above, but guarantees that result is correct
1253
1254 void IrredPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F, long m);
1255
IrredPolyMod(const zz_pX & g,const zz_pXModulus & F,long m)1256 inline zz_pX IrredPolyMod(const zz_pX& g, const zz_pXModulus& F, long m)
1257 { zz_pX x; IrredPolyMod(x, g, F, m); NTL_OPT_RETURN(zz_pX, x); }
1258
1259
IrredPolyMod(zz_pX & h,const zz_pX & g,const zz_pXModulus & F)1260 inline void IrredPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F)
1261 { IrredPolyMod(h, g, F, F.n); }
1262
IrredPolyMod(const zz_pX & g,const zz_pXModulus & F)1263 inline zz_pX IrredPolyMod(const zz_pX& g, const zz_pXModulus& F)
1264 { zz_pX x; IrredPolyMod(x, g, F); NTL_OPT_RETURN(zz_pX, x); }
1265
1266
1267 // same as above, but assumes that f is irreducible,
1268 // or at least that the minimal poly of g is itself irreducible.
1269 // The algorithm is deterministic (and is always correct).
1270
1271
1272
1273 /*****************************************************************
1274
1275 Traces, norms, resultants
1276
1277 ******************************************************************/
1278
1279 void TraceVec(vec_zz_p& S, const zz_pX& f);
1280
TraceVec(const zz_pX & f)1281 inline vec_zz_p TraceVec(const zz_pX& f)
1282 { vec_zz_p x; TraceVec(x, f); NTL_OPT_RETURN(vec_zz_p, x); }
1283
1284
1285 void FastTraceVec(vec_zz_p& S, const zz_pX& f);
1286 void PlainTraceVec(vec_zz_p& S, const zz_pX& f);
1287
1288 void TraceMod(zz_p& x, const zz_pX& a, const zz_pXModulus& F);
1289
TraceMod(const zz_pX & a,const zz_pXModulus & F)1290 inline zz_p TraceMod(const zz_pX& a, const zz_pXModulus& F)
1291 { zz_p x; TraceMod(x, a, F); return x; }
1292
1293 void TraceMod(zz_p& x, const zz_pX& a, const zz_pX& f);
1294
TraceMod(const zz_pX & a,const zz_pX & f)1295 inline zz_p TraceMod(const zz_pX& a, const zz_pX& f)
1296 { zz_p x; TraceMod(x, a, f); return x; }
1297
1298
1299
1300
1301 void ComputeTraceVec(const zz_pXModulus& F);
1302
1303 void NormMod(zz_p& x, const zz_pX& a, const zz_pX& f);
1304
1305
NormMod(const zz_pX & a,const zz_pX & f)1306 inline zz_p NormMod(const zz_pX& a, const zz_pX& f)
1307 { zz_p x; NormMod(x, a, f); return x; }
1308
1309 void resultant(zz_p& rres, const zz_pX& a, const zz_pX& b);
1310
resultant(const zz_pX & a,const zz_pX & b)1311 inline zz_p resultant(const zz_pX& a, const zz_pX& b)
1312 { zz_p x; resultant(x, a, b); return x; }
1313
1314 void CharPolyMod(zz_pX& g, const zz_pX& a, const zz_pX& f);
1315 // g = char poly of (a mod f)
1316 // only implemented for p >= deg(f)+1
1317
CharPolyMod(const zz_pX & a,const zz_pX & f)1318 inline zz_pX CharPolyMod(const zz_pX& a, const zz_pX& f)
1319 { zz_pX x; CharPolyMod(x, a, f); NTL_OPT_RETURN(zz_pX, x); }
1320
1321
1322
1323
1324
1325 NTL_CLOSE_NNS
1326
1327 #endif
1328