1 // -*- C++ -*-
2 //==============================================================================================
3 //
4 // This file is part of LiDIA --- a library for computational number theory
5 //
6 // Copyright (c) 1994--2001 the LiDIA Group. All rights reserved.
7 //
8 // See http://www.informatik.tu-darmstadt.de/TI/LiDIA/
9 //
10 //----------------------------------------------------------------------------------------------
11 //
12 // $Id$
13 //
14 // Author : Victor Shoup (VS), Thomas Pfahler (TPf)
15 // Changes : See CVS log
16 //
17 //==============================================================================================
18
19
20 #ifndef LIDIA_FP_POLYNOMIAL_H_GUARD_
21 #define LIDIA_FP_POLYNOMIAL_H_GUARD_
22
23
24 #ifndef LIDIA_LIDIA_H_GUARD_
25 # include "LiDIA/LiDIA.h"
26 #endif
27 # ifndef LIDIA_COMPARATOR_H_GUARD_
28 # include "LiDIA/comparator.h"
29 #endif
30 #ifndef LIDIA_BASE_VECTOR_H_GUARD_
31 # include "LiDIA/base_vector.h"
32 #endif
33 #ifndef LIDIA_BIGINT_H_GUARD_
34 # include "LiDIA/bigint.h"
35 #endif
36 #ifndef LIDIA_BIGINT_POLYNOMIAL_H_GUARD_
37 # include "LiDIA/bigint_polynomial.h"
38 #endif
39 #ifndef LIDIA_RESIDUE_CLASS_LIST_H_GUARD_
40 # include "LiDIA/base/residue_class_list.h"
41 #endif
42
43
44
45 #ifdef LIDIA_NAMESPACE
46 namespace LiDIA {
47 # define IN_NAMESPACE_LIDIA
48 #endif
49
50
51 // forward declarations required later on
52
53 class Fp_polynomial;
54 class Fp_pol_ref;
55 class base_fft_rep;
56 class fft_rep;
57 class modular_fft_rep;
58 class Fp_poly_modulus;
59 class Fp_poly_multiplier;
60 class poly_matrix;
61 class poly_mod_rep;
62 class eq_frac_info;
63 class poly_argument;
64
65
66
67 #ifndef HEADBANGER
68
69
70 inline void
Remainder(bigint & x,const bigint & a,const bigint & p)71 Remainder (bigint & x, const bigint & a, const bigint & p)
72 {
73 remainder(x, a, p);
74 if (x.is_lt_zero()) {
75 add(x, x, p);
76 }
77 }
78
79
80
81 inline void
AddMod(bigint & c,const bigint & a,const bigint & b,const bigint & p)82 AddMod (bigint & c, const bigint & a, const bigint & b, const bigint & p)
83 {
84 add(c, a, b);
85 if (c >= p) {
86 subtract(c, c, p);
87 }
88 }
89
90
91
92 inline void
SubMod(bigint & c,const bigint & a,const bigint & b,const bigint & p)93 SubMod (bigint & c, const bigint & a, const bigint & b, const bigint & p)
94 {
95 subtract(c, a, b);
96 if (c.is_lt_zero()) {
97 add(c, c, p);
98 }
99 }
100
101
102
103 inline void
MulMod(bigint & c,const bigint & a,const bigint & b,const bigint & p)104 MulMod (bigint & c, const bigint & a, const bigint & b, const bigint & p)
105 {
106 multiply(c, a, b);
107 Remainder(c, c, p);
108 }
109
110
111
112 inline void
NegateMod(bigint & x,const bigint & a,const bigint & p)113 NegateMod (bigint & x, const bigint & a, const bigint & p)
114 {
115 if (a.is_zero()) {
116 x.assign_zero();
117 }
118 else {
119 subtract(x, p, a);
120 }
121 }
122
123
124
125 inline void
InvMod(bigint & x,const bigint & a,const bigint & p)126 InvMod(bigint & x, const bigint & a, const bigint & p)
127 {
128 // x = a^{-1} mod n, 0 <= x < p
129 // error is raised if inverse not defined
130 if (xgcd_left(x, a, p).is_one()) {
131 Remainder(x, x, p);
132 return;
133 }
134 else {
135 lidia_error_handler("Fp_polynomial", "InvMod(...)::inverse does not exist");
136 }
137 }
138
139 #endif // HEADBANGER
140
141
142
143
144 //************************************************************
145 // class Fp_pol_ref
146 // tries to differentiate between the use of
147 // Fp_polynomial::operator[] as an l-value and
148 // as an r-value
149 //************************************************************
150
151 class Fp_pol_ref
152 {
153 friend class Fp_polynomial; // MM
154
155
156 private:
157
158 class Fp_polynomial &p;
159 lidia_size_t ix;
160
161 Fp_pol_ref(); // not needed
Fp_pol_ref(const Fp_pol_ref & a)162 Fp_pol_ref (const Fp_pol_ref & a) : p(a.p), ix(a.ix)
163 { }
164
165
166 public:
167
Fp_pol_ref(Fp_polynomial & f,lidia_size_t i)168 Fp_pol_ref(Fp_polynomial &f, lidia_size_t i) : p(f), ix(i)
169 { }
170
171 //l-values
172 Fp_pol_ref & operator = (const bigint &a);
173 void assign_zero ();
174 void assign_one ();
175 void assign (const bigint &a);
176
177 //r-value
178 operator bigint ();
179 };
180
181
182
183 #ifndef HEADBANGER
184
185 //************************************************************
186 // class crossover_class
187 // class for computing crossover-points with respect to
188 // bitlengths of used moduli
189 //************************************************************
190
191 #define CROV_NUM_VALUES 10
192
193 class crossover_class
194 {
195 int x[CROV_NUM_VALUES];
196 int y_fftmul[CROV_NUM_VALUES];
197 int y_fftdiv[CROV_NUM_VALUES];
198 int y_fftrem[CROV_NUM_VALUES];
199 int y_inv[CROV_NUM_VALUES];
200 int y_gcd[CROV_NUM_VALUES];
201 int y_xgcd[CROV_NUM_VALUES];
202 double y2_fftmul[CROV_NUM_VALUES];
203 double y2_fftdiv[CROV_NUM_VALUES];
204 double y2_fftrem[CROV_NUM_VALUES];
205 double y2_inv[CROV_NUM_VALUES];
206 double y2_gcd[CROV_NUM_VALUES];
207 double y2_xgcd[CROV_NUM_VALUES];
208 int halfgcd;
209 int log2_newton;
210
211 crossover_class (const crossover_class &); //disable
212
213
214 public :
215
216 crossover_class ();
217 void init(const int x_val[CROV_NUM_VALUES],
218 const int fftmul_val[CROV_NUM_VALUES],
219 const int fftdiv_val[CROV_NUM_VALUES],
220 const int fftrem_val[CROV_NUM_VALUES],
221 const int inv_val[CROV_NUM_VALUES],
222 const int gcd_val[CROV_NUM_VALUES],
223 const int xgcd_val[CROV_NUM_VALUES]);
224
225 int fftmul_crossover (const bigint &modulus) const;
226 int fftdiv_crossover (const bigint &modulus) const;
227 int fftrem_crossover (const bigint &modulus) const;
228 int inv_crossover (const bigint &modulus) const;
229 int halfgcd_crossover (const bigint &modulus) const;
230 int gcd_crossover (const bigint &modulus) const;
231 int xgcd_crossover (const bigint &modulus) const;
232 int log2_newton_crossover (const bigint &modulus) const;
233 };
234
235 #endif // HEADBANGER
236
237
238
239 //************************************************************
240 //
241 // Fp_polynomial
242 //
243 // The class Fp_polynomial implements polynomial arithmetic modulo p.
244 // Polynomials are represented as vectors of bigint with values in the
245 // range 0..p-1.
246 //
247 // If f is a Fp_polynomial, then f.coeff is a bigint*.
248 // The zero polynomial is represented as a zero length vector.
249 // Otherwise, f.coeff[0] is the constant-term, and f.coeff[f.c_length-1]
250 // is the leading coefficient, which is always non-zero.
251 // Use the member function remove_leading_zeros() to strip leading zeros.
252 //
253 //**************************************************************
254
255
256
257 class Fp_polynomial
258 {
259 #ifndef HEADBANGER
260 private:
261 static bigint ZERO; //coeff. value for zero polynomial
262 static residue_class_list< bigint > *L; //list of moduli
263
264 bigint *coeff; //coefficient vector
265 lidia_size_t c_length;
266 lidia_size_t c_alloc;
267 residue_class< bigint > *Mp; //pointer to modulus
268 #endif
269
set_degree(lidia_size_t n)270 void set_degree (lidia_size_t n)
271 {
272 debug_handler_l("Fp_polynomial", "set_degree (lidia_size_t)", 2);
273 set_max_degree(n);
274 c_length = comparator< lidia_size_t >::max(n+1, 0);
275 }
276
277 void read_x_term (std::istream &datei, bigint &coeff, lidia_size_t &exp);
278
279
280 public:
281
282 static crossover_class crossovers; // crossover-points for fft-arithm.
283 // depending on architecture
284
285
286 //***************************************************************
287 //
288 // Constructors, Destructors, and Basic Functions
289 //
290 //***************************************************************
291
Fp_polynomial()292 Fp_polynomial () :
293 coeff(0),
294 c_length(0),
295 c_alloc(0),
296 Mp(0)
297 {
298 debug_handler_l("Fp_polynomial", "Fp_polynomial()", 1);
299 }
300
301 Fp_polynomial (const Fp_polynomial& a);
302 Fp_polynomial (const polynomial< bigint > & f, const bigint & p);
303
304 ~Fp_polynomial();
305
306 void kill();
307
308 void set_modulus (const Fp_polynomial &f);
309 //sets modulus to f.modulus(), assigns zero polynomial
310 void set_modulus (const bigint &p);
311 //sets modulus to p, assigns zero polynomial
312
313 const bigint & modulus() const;
314
degree()315 lidia_size_t degree () const
316 // note that the zero polynomial has degree -1.
317 {
318 debug_handler_l("Fp_polynomial", "degree(void)", 2);
319 return c_length - 1;
320 }
321
322 #ifndef HEADBANGER
323 void comp_modulus (const Fp_polynomial& x, const char* fctname) const;
324 void set_max_degree (lidia_size_t n);
325 #endif // HEADBANGER
326
327 //***************************************************************
328 //
329 // routines for manipulating coefficients
330 //
331 //***************************************************************
332
333 #ifndef HEADBANGER
334 void get_coefficient (bigint& a, lidia_size_t i) const;
335 void set_coefficient (const bigint& a, lidia_size_t i);
336 void set_coefficient (lidia_size_t i);
337 #endif // HEADBANGER
338
339 const bigint& lead_coeff () const;
340 const bigint& const_term () const;
341
342 void remove_leading_zeros ();
343 void make_monic ();
344
345 friend class Fp_pol_ref;
346 Fp_pol_ref operator[] (lidia_size_t i)
347 {
348 return Fp_pol_ref(*this, i);
349 }
350 // operator[] with context sensitivity
351 // tries to differentiate between use as l-value and r-value
352
353 const bigint & operator[] (lidia_size_t i) const;
354 //always r-value
355
356
357 //***************************************************************
358 //
359 // assignments
360 //
361 //****************************************************************/
362
363 Fp_polynomial & operator = (const Fp_polynomial& a);
364 Fp_polynomial & operator = (const bigint& c);
365
366 #ifndef HEADBANGER
367 void assign (const bigint & a);
368 //void assign (const bigmod& a);
369 //void assign (const base_vector < bigmod > & a);
370 void assign (const base_vector< bigint > & a, const bigint &p);
371 void assign (const Fp_polynomial& a);
372 #endif // HEADBANGER
373
assign_zero()374 void assign_zero ()
375 {
376 debug_handler_l ("Fp_polynomial", "assign_zero (void)", 2);
377 if (modulus().is_zero()) {
378 lidia_error_handler("Fp_polynomial", "assign_zero(void)::modulus = 0");
379 }
380
381 c_length = 0; //don't delete the coeff. vector !
382 }
383
assign_one()384 void assign_one ()
385 {
386 debug_handler_l ("Fp_polynomial", "assign_one (void)", 2);
387 if (modulus().is_zero()) {
388 lidia_error_handler("Fp_polynomial", "assign_one(void)::modulus = 0");
389 }
390
391 set_degree(0);
392 coeff[0].assign_one();
393 }
394
395 void assign_x ();
396
397 void assign_zero (const bigint &p); //sets modulus to p
398 void assign_one (const bigint &p); // "
399 void assign_x (const bigint &p); // "
400
401 void randomize (lidia_size_t n); //random polynomial of degree n
402
403
404 //***************************************************************
405 //
406 // comparisons
407 //
408 //***************************************************************
409
410 friend bool operator == (const Fp_polynomial& a, const Fp_polynomial& b);
411 friend bool operator != (const Fp_polynomial& a, const Fp_polynomial& b);
412
413 bool is_zero () const;
414 bool is_one () const;
415 bool is_x () const;
416
is_monic()417 bool is_monic () const
418 {
419 return (lead_coeff().is_one());
420 }
421
422 bool is_binomial() const;
423
424
425 //***************************************************************
426 //
427 // operators
428 //
429 //***************************************************************
430
431 Fp_polynomial & operator += (const Fp_polynomial &a);
432 Fp_polynomial & operator += (const bigint &a);
433 Fp_polynomial & operator -= (const Fp_polynomial &a);
434 Fp_polynomial & operator -= (const bigint &a);
435 Fp_polynomial & operator *= (const Fp_polynomial &a);
436 Fp_polynomial & operator *= (const bigint &a);
437 Fp_polynomial & operator /= (const Fp_polynomial &a);
438 Fp_polynomial & operator /= (const bigint &a);
439 Fp_polynomial & operator %= (const Fp_polynomial &a);
440
441 bigint operator () (const bigint &value) const;
442
443
444 //***************************************************************
445 //
446 // procedural versions for arithmetics
447 //
448 //***************************************************************
449
450
451 #ifndef HEADBANGER
452 void negate ();
453 friend void negate (Fp_polynomial& x, const Fp_polynomial& a);
454 friend void add (Fp_polynomial& x, const Fp_polynomial& a,
455 const Fp_polynomial& b);
456 friend void add (Fp_polynomial & x, const Fp_polynomial& a,
457 const bigint& b);
458 friend void add (Fp_polynomial& x, const bigint& a,
459 const Fp_polynomial& b);
460
461 friend void subtract (Fp_polynomial& x, const Fp_polynomial& a,
462 const Fp_polynomial& b);
463 friend void subtract (Fp_polynomial & x, const Fp_polynomial& a,
464 const bigint& b);
465 friend void subtract (Fp_polynomial& x, const bigint& a,
466 const Fp_polynomial& b);
467
468 friend void multiply (Fp_polynomial& x, const Fp_polynomial& a,
469 const Fp_polynomial& b);
470 friend void multiply (Fp_polynomial & x, const Fp_polynomial& a,
471 const bigint& b);
472 friend void multiply (Fp_polynomial & x, const bigint& a,
473 const Fp_polynomial& b);
474 friend void multiply_by_scalar (Fp_polynomial &c,
475 const Fp_polynomial &a,
476 const bigint &b);
477
478 // These always use "classical" arithmetic
479 friend void plain_mul (Fp_polynomial& x, const Fp_polynomial& a,
480 const Fp_polynomial& b);
481 friend void plain_sqr (Fp_polynomial& x, const Fp_polynomial& a);
482
483 // These always use FFT arithmetic
484 friend void fft_mul (Fp_polynomial& x, const Fp_polynomial& a,
485 const Fp_polynomial& b);
486 friend void fft_sqr (Fp_polynomial& x, const Fp_polynomial& a);
487 #endif // HEADBANGER
488
489 friend void square(Fp_polynomial& x, const Fp_polynomial& a);
490
491
492
493 #ifndef HEADBANGER
494 friend void div_rem (Fp_polynomial& q, Fp_polynomial& r,
495 const Fp_polynomial& a, const Fp_polynomial& b);
496 // q = a/b, r = a%b
497
498 friend void divide (Fp_polynomial& q, const Fp_polynomial& a,
499 const Fp_polynomial& b);
500 friend void divide (Fp_polynomial& q, const bigint& a,
501 const Fp_polynomial& b);
502 friend void divide (Fp_polynomial& q, const Fp_polynomial& a,
503 const bigint& b);
504 // q = a/b
505
506 friend void remainder (Fp_polynomial& r, const Fp_polynomial& a,
507 const Fp_polynomial& b);
508 // r = a%b
509 #endif // HEADBANGER
510
511 friend void invert (Fp_polynomial& c, const Fp_polynomial& a,
512 lidia_size_t m);
513 // computes c = a^{-1} % x^m
514 // constant term must be non-zero
515
516
517 #ifndef HEADBANGER
518 // These always use "classical" arithmetic
519 friend void plain_div_rem (Fp_polynomial& q, Fp_polynomial& r,
520 const Fp_polynomial& a,
521 const Fp_polynomial& b);
522 friend void plain_div (Fp_polynomial& q, const Fp_polynomial& a,
523 const Fp_polynomial& b);
524 friend void plain_rem (Fp_polynomial& r, const Fp_polynomial& a,
525 const Fp_polynomial& b);
526
527 // These always use FFT arithmetic
528 friend void fft_div_rem (Fp_polynomial& q, Fp_polynomial& r,
529 const Fp_polynomial& a,
530 const Fp_polynomial& b);
531 friend void fft_div (Fp_polynomial& q, const Fp_polynomial& a,
532 const Fp_polynomial& b);
533 friend void fft_rem (Fp_polynomial& r, const Fp_polynomial& a,
534 const Fp_polynomial& b);
535
536 friend void plain_inv (Fp_polynomial& x, const Fp_polynomial& a,
537 lidia_size_t m);
538 // always uses "classical" algorithm
539 // ALIAS RESTRICTION: input may not alias output
540
541 friend void newton_inv (Fp_polynomial& x, const Fp_polynomial& a,
542 lidia_size_t m);
543 // uses a Newton Iteration with the FFT.
544 // ALIAS RESTRICTION: input may not alias output
545 #endif // HEADBANGER
546
547
548
549 //***************************************************************
550 //
551 // Miscellaneaous
552 //
553 //****************************************************************/
554
555 friend void swap (Fp_polynomial& x, Fp_polynomial& y);
556 // swap x & y (only pointers are swapped)
557
558 friend void trunc (Fp_polynomial& c, const Fp_polynomial& a,
559 lidia_size_t m);
560 // c = a % x^m
561
562 friend void shift_right (Fp_polynomial& c, const Fp_polynomial& a,
563 lidia_size_t n);
564 // c = a/x^n
565
566 friend void shift_left (Fp_polynomial& c, const Fp_polynomial& a,
567 lidia_size_t n);
568 // c = a*x^n
569
570 friend void derivative (Fp_polynomial& c, const Fp_polynomial& a);
571 // c = derivative of a
572 friend Fp_polynomial derivative (const Fp_polynomial& a);
573
574 friend void add_multiple (Fp_polynomial &f,
575 const Fp_polynomial &g,
576 const bigint &s,
577 lidia_size_t n,
578 const Fp_polynomial &h);
579 // f = g + s*x^n*h, n >= 0
580
581 friend void cyclic_reduce (Fp_polynomial& x, const Fp_polynomial& a,
582 lidia_size_t m);
583 // computes c = a mod x^m-1
584
585 friend void copy_reverse (Fp_polynomial& x, const Fp_polynomial& a,
586 lidia_size_t lo, lidia_size_t hi);
587 // x[0..hi-lo+1] = reverse(a[lo..hi]), with zero fill
588 // input may not alias output
589
590 void build_from_roots (const base_vector< bigint > & a);
591 // computes the polynomial (x-a[0]) ... (x-a[n-1]),
592 // where n = a.length()
593 // implementation: see file "fftbuild.c"
594
595 friend void multiply_by_x_mod (Fp_polynomial& c,
596 const Fp_polynomial& a,
597 const Fp_polynomial& f);
598 // c = (a * x) mod f
599
600
601 //********************************************************************
602 //
603 // I/O
604 //
605 // two I/O formats:
606 //
607 // 1) [a_n ... a_1 a_0] mod p
608 // represents the polynomial a_n*x^n + ... + a_1*x + a_0 mod p.
609 //
610 // 2) a_n*x^n + ... + a_1*x + a_0 mod p
611 // is self-explanatory.
612 // The default is 2).
613 //
614 // On output, all coefficients will be integers between 0 and p-1,
615 // and a_n not zero (the zero polynomial is [ ]).
616 // On input, the coefficients are arbitrary integers which are
617 // then reduced modulo p, and leading zeros stripped. The function
618 // Fp_polynomial::read(std::istream&) can handle both formats; the decision
619 // is made according to the first character being '[' or not.
620 //
621 //*********************************************************************
622
623 void read (std::istream & s);
624 void write (std::ostream & s) const; //writes format 1)
625 void pretty_read (std::istream & s); //reads format 2)
626 void pretty_print (std::ostream & s = std::cout) const; //writes format 2)
627
628
629 //***************************************************************
630 //
631 // special functions that have to be friends
632 //
633 //***************************************************************
634
635 friend class base_fft_rep;
636 friend class fft_rep;
637 friend class modular_fft_rep;
638 friend class Fp_poly_modulus;
639 friend class Fp_poly_multiplier;
640 friend class poly_matrix;
641 friend class poly_mod_rep;
642 friend class eq_frac_info;
643 friend class poly_argument;
644
645
646 friend void remainder (Fp_polynomial &, const Fp_polynomial &, const Fp_poly_modulus &);
647 //set_degree
648
649
650 }; //end class Fp_polynomial
651
652 // declaration of Fp_polynomial's friends
653
654 bool operator == (const Fp_polynomial& a, const Fp_polynomial& b);
655 bool operator != (const Fp_polynomial& a, const Fp_polynomial& b);
656
657 void negate (Fp_polynomial& x, const Fp_polynomial& a);
658 void add (Fp_polynomial& x, const Fp_polynomial& a,
659 const Fp_polynomial& b);
660 void add (Fp_polynomial & x, const Fp_polynomial& a,
661 const bigint& b);
662 void add (Fp_polynomial& x, const bigint& a,
663 const Fp_polynomial& b);
664
665 void subtract (Fp_polynomial& x, const Fp_polynomial& a,
666 const Fp_polynomial& b);
667 void subtract (Fp_polynomial & x, const Fp_polynomial& a,
668 const bigint& b);
669 void subtract (Fp_polynomial& x, const bigint& a,
670 const Fp_polynomial& b);
671
672 void multiply (Fp_polynomial& x, const Fp_polynomial& a,
673 const Fp_polynomial& b);
674 void multiply (Fp_polynomial & x, const Fp_polynomial& a,
675 const bigint& b);
676 void multiply (Fp_polynomial & x, const bigint& a,
677 const Fp_polynomial& b);
678 void multiply_by_scalar (Fp_polynomial &c,
679 const Fp_polynomial &a,
680 const bigint &b);
681
682 // These always use "classical" arithmetic
683 void plain_mul (Fp_polynomial& x, const Fp_polynomial& a,
684 const Fp_polynomial& b);
685 void plain_sqr (Fp_polynomial& x, const Fp_polynomial& a);
686
687 // These always use FFT arithmetic
688 void fft_mul (Fp_polynomial& x, const Fp_polynomial& a,
689 const Fp_polynomial& b);
690 void fft_sqr (Fp_polynomial& x, const Fp_polynomial& a);
691
692 void square(Fp_polynomial& x, const Fp_polynomial& a);
693
694
695
696 #ifndef HEADBANGER
697 void div_rem (Fp_polynomial& q, Fp_polynomial& r,
698 const Fp_polynomial& a, const Fp_polynomial& b);
699 // q = a/b, r = a%b
700
701 void divide (Fp_polynomial& q, const Fp_polynomial& a,
702 const Fp_polynomial& b);
703 void divide (Fp_polynomial& q, const bigint& a,
704 const Fp_polynomial& b);
705 void divide (Fp_polynomial& q, const Fp_polynomial& a,
706 const bigint& b);
707 // q = a/b
708
709 void remainder (Fp_polynomial& r, const Fp_polynomial& a,
710 const Fp_polynomial& b);
711 // r = a%b
712 #endif // HEADBANGER
713
714 void invert (Fp_polynomial& c, const Fp_polynomial& a,
715 lidia_size_t m);
716 // computes c = a^{-1} % x^m
717 // constant term must be non-zero
718
719
720 #ifndef HEADBANGER
721 // These always use "classical" arithmetic
722 void plain_div_rem (Fp_polynomial& q, Fp_polynomial& r,
723 const Fp_polynomial& a,
724 const Fp_polynomial& b);
725 void plain_div (Fp_polynomial& q, const Fp_polynomial& a,
726 const Fp_polynomial& b);
727 void plain_rem (Fp_polynomial& r, const Fp_polynomial& a,
728 const Fp_polynomial& b);
729
730 // These always use FFT arithmetic
731 void fft_div_rem (Fp_polynomial& q, Fp_polynomial& r,
732 const Fp_polynomial& a,
733 const Fp_polynomial& b);
734 void fft_div (Fp_polynomial& q, const Fp_polynomial& a,
735 const Fp_polynomial& b);
736 void fft_rem (Fp_polynomial& r, const Fp_polynomial& a,
737 const Fp_polynomial& b);
738
739 void plain_inv (Fp_polynomial& x, const Fp_polynomial& a,
740 lidia_size_t m);
741 // always uses "classical" algorithm
742 // ALIAS RESTRICTION: input may not alias output
743
744 void newton_inv (Fp_polynomial& x, const Fp_polynomial& a,
745 lidia_size_t m);
746 // uses a Newton Iteration with the FFT.
747 // ALIAS RESTRICTION: input may not alias output
748 #endif // HEADBANGER
749
750
751
752 //***************************************************************
753 //
754 // Miscellaneaous
755 //
756 //****************************************************************/
757
758 void swap (Fp_polynomial& x, Fp_polynomial& y);
759 // swap x & y (only pointers are swapped)
760
761 void trunc (Fp_polynomial& c, const Fp_polynomial& a,
762 lidia_size_t m);
763 // c = a % x^m
764
765 void shift_right (Fp_polynomial& c, const Fp_polynomial& a,
766 lidia_size_t n);
767 // c = a/x^n
768
769 void shift_left (Fp_polynomial& c, const Fp_polynomial& a,
770 lidia_size_t n);
771 // c = a*x^n
772
773 void derivative (Fp_polynomial& c, const Fp_polynomial& a);
774 // c = derivative of a
775 Fp_polynomial derivative (const Fp_polynomial& a);
776
777 void add_multiple (Fp_polynomial &f,
778 const Fp_polynomial &g,
779 const bigint &s,
780 lidia_size_t n,
781 const Fp_polynomial &h);
782 // f = g + s*x^n*h, n >= 0
783
784 void cyclic_reduce (Fp_polynomial& x, const Fp_polynomial& a,
785 lidia_size_t m);
786 // computes c = a mod x^m-1
787
788 void copy_reverse (Fp_polynomial& x, const Fp_polynomial& a,
789 lidia_size_t lo, lidia_size_t hi);
790 // x[0..hi-lo+1] = reverse(a[lo..hi]), with zero fill
791 // input may not alias output
792
793 void multiply_by_x_mod (Fp_polynomial& c,
794 const Fp_polynomial& a,
795 const Fp_polynomial& f);
796 // c = (a * x) mod f
797
798
799 void remainder (Fp_polynomial &, const Fp_polynomial &,
800 const Fp_poly_modulus &);
801
802 //***************************************************************
803 //
804 // operators
805 //
806 //***************************************************************
807
808 Fp_polynomial operator - (const Fp_polynomial &a);
809 Fp_polynomial operator + (const Fp_polynomial &a, const Fp_polynomial &b);
810 Fp_polynomial operator + (const bigint &a, const Fp_polynomial &b);
811 Fp_polynomial operator + (const Fp_polynomial &a, const bigint &b);
812 Fp_polynomial operator - (const Fp_polynomial &a, const Fp_polynomial &b);
813 Fp_polynomial operator - (const bigint &a, const Fp_polynomial &b);
814 Fp_polynomial operator - (const Fp_polynomial &a, const bigint &b);
815 Fp_polynomial operator * (const Fp_polynomial &a, const Fp_polynomial &b);
816 Fp_polynomial operator * (const bigint &a, const Fp_polynomial &b);
817 Fp_polynomial operator * (const Fp_polynomial &a, const bigint &b);
818 Fp_polynomial operator / (const Fp_polynomial &a, const Fp_polynomial &b);
819 Fp_polynomial operator / (const bigint &a, const Fp_polynomial &b);
820 Fp_polynomial operator / (const Fp_polynomial &a, const bigint &b);
821 Fp_polynomial operator % (const Fp_polynomial &a, const Fp_polynomial &b);
822
823
824 // for class factorization < Fp_polynomial > :
825 // used in class single_factor < Fp_polynomial > :
826 bool operator < (const Fp_polynomial &a, const Fp_polynomial &b);
827 bool operator <= (const Fp_polynomial &a, const Fp_polynomial &b);
828 bool operator == (const Fp_polynomial& a, const Fp_polynomial& b);
829 bool operator != (const Fp_polynomial& a, const Fp_polynomial& b);
830
831 void negate (Fp_polynomial& x, const Fp_polynomial& a);
832 void add (Fp_polynomial& x, const Fp_polynomial& a, const Fp_polynomial& b);
833 void add (Fp_polynomial & x, const Fp_polynomial& a, const bigint& b);
834 inline
add(Fp_polynomial & x,const bigint & a,const Fp_polynomial & b)835 void add (Fp_polynomial& x, const bigint& a, const Fp_polynomial& b)
836 {
837 add(x, b, a);
838 }
839
840 void subtract (Fp_polynomial& x, const Fp_polynomial& a,
841 const Fp_polynomial& b);
842 void subtract (Fp_polynomial & x, const Fp_polynomial& a, const bigint& b);
843 inline
subtract(Fp_polynomial & x,const bigint & a,const Fp_polynomial & b)844 void subtract (Fp_polynomial& x, const bigint& a, const Fp_polynomial& b)
845 {
846 LiDIA::negate(x, b);
847 add(x, x, a);
848 }
849
850
851 void multiply (Fp_polynomial& x, const Fp_polynomial& a,
852 const Fp_polynomial& b);
853 inline
multiply(Fp_polynomial & x,const Fp_polynomial & a,const bigint & b)854 void multiply (Fp_polynomial & x, const Fp_polynomial& a, const bigint& b)
855 {
856 multiply_by_scalar(x, a, b);
857 }
858 inline
multiply(Fp_polynomial & x,const bigint & a,const Fp_polynomial & b)859 void multiply (Fp_polynomial & x, const bigint& a, const Fp_polynomial& b)
860 {
861 multiply_by_scalar(x, b, a);
862 }
863 void multiply_by_scalar (Fp_polynomial &c, const Fp_polynomial &a,
864 const bigint &b);
865
866 // These always use "classical" arithmetic
867 void plain_mul (Fp_polynomial& x, const Fp_polynomial& a,
868 const Fp_polynomial& b);
869 void plain_sqr (Fp_polynomial& x, const Fp_polynomial& a);
870
871 // These always use FFT arithmetic
872 void fft_mul (Fp_polynomial& x, const Fp_polynomial& a,
873 const Fp_polynomial& b);
874 void fft_sqr (Fp_polynomial& x, const Fp_polynomial& a);
875
876 void square(Fp_polynomial& x, const Fp_polynomial& a);
877
878 // q = a/b, r = a%b
879 void div_rem (Fp_polynomial& q, Fp_polynomial& r, const Fp_polynomial& a,
880 const Fp_polynomial& b);
881
882 // q = a/b
883 void divide (Fp_polynomial& q, const Fp_polynomial& a, const Fp_polynomial& b);
884 void divide (Fp_polynomial& q, const bigint& a, const Fp_polynomial& b);
885 void divide (Fp_polynomial& q, const Fp_polynomial& a, const bigint& b);
886
887 // r = a%b
888 void remainder (Fp_polynomial& r, const Fp_polynomial& a,
889 const Fp_polynomial& b);
890
891
892 // computes c = a^{-1} % x^m
893 // constant term must be non-zero
894 void invert (Fp_polynomial& c, const Fp_polynomial& a, lidia_size_t m);
895
896 // These always use "classical" arithmetic
897 void plain_div_rem (Fp_polynomial& q, Fp_polynomial& r, const Fp_polynomial& a,
898 const Fp_polynomial& b);
899 void plain_div (Fp_polynomial& q, const Fp_polynomial& a,
900 const Fp_polynomial& b);
901 void plain_rem (Fp_polynomial& r, const Fp_polynomial& a,
902 const Fp_polynomial& b);
903
904 // These always use FFT arithmetic
905 void fft_div_rem (Fp_polynomial& q, Fp_polynomial& r, const Fp_polynomial& a,
906 const Fp_polynomial& b);
907 void fft_div (Fp_polynomial& q, const Fp_polynomial& a,
908 const Fp_polynomial& b);
909 void fft_rem (Fp_polynomial& r, const Fp_polynomial& a,
910 const Fp_polynomial& b);
911
912 // always uses "classical" algorithm
913 // ALIAS RESTRICTION: input may not alias output
914 void plain_inv (Fp_polynomial& x, const Fp_polynomial& a, lidia_size_t m);
915
916 // uses a Newton Iteration with the FFT.
917 // ALIAS RESTRICTION: input may not alias output
918 void newton_inv (Fp_polynomial& x, const Fp_polynomial& a, lidia_size_t m);
919
920
921 //***************************************************************
922 //
923 // Miscellaneaous
924 //
925 //****************************************************************/
926
927 // swap x & y (only pointers are swapped)
928 void swap (Fp_polynomial& x, Fp_polynomial& y);
929
930 // c = a % x^m
931 void trunc (Fp_polynomial& c, const Fp_polynomial& a, lidia_size_t m);
932
933 // c = a/x^n
934 void shift_right (Fp_polynomial& c, const Fp_polynomial& a, lidia_size_t n);
935
936 // c = a*x^n
937 void shift_left (Fp_polynomial& c, const Fp_polynomial& a, lidia_size_t n);
938
939 // c = derivative of a
940 void derivative (Fp_polynomial& c, const Fp_polynomial& a);
941 inline
derivative(const Fp_polynomial & a)942 Fp_polynomial derivative (const Fp_polynomial& a)
943 {
944 Fp_polynomial x;
945
946 derivative(x, a);
947 return x;
948 }
949
950 // f = g + s*x^n*h, n >= 0
951 void add_multiple (Fp_polynomial &f, const Fp_polynomial &g,
952 const bigint &s, lidia_size_t n, const Fp_polynomial &h);
953
954 // computes c = a mod x^m-1
955 void cyclic_reduce (Fp_polynomial& x, const Fp_polynomial& a, lidia_size_t m);
956
957 // x[0..hi-lo+1] = reverse(a[lo..hi]), with zero fill
958 // input may not alias output
959 void copy_reverse (Fp_polynomial& x, const Fp_polynomial& a, lidia_size_t lo,
960 lidia_size_t hi);
961
962 // c = (a * x) mod f
963 void multiply_by_x_mod (Fp_polynomial& c, const Fp_polynomial& a,
964 const Fp_polynomial& f);
965
966 //set_degree
967 void remainder (Fp_polynomial &, const Fp_polynomial &,
968 const Fp_poly_modulus &);
969
970
971 //***************************************************************
972 //
973 // gcd's
974 //
975 //***************************************************************
976
977 // x = gcd(a, b), x is always monic (or zero if a == b == 0).
978 void gcd (Fp_polynomial& x, const Fp_polynomial& a, const Fp_polynomial& b);
979
980 inline Fp_polynomial
gcd(const Fp_polynomial & a,const Fp_polynomial & b)981 gcd(const Fp_polynomial& a, const Fp_polynomial& b)
982 {
983 Fp_polynomial x;
984
985 gcd(x, a, b);
986 return x;
987 }
988
989
990
991 void
992 xgcd (Fp_polynomial& d, Fp_polynomial& s, Fp_polynomial& t,
993 const Fp_polynomial& a, const Fp_polynomial& b);
994 // d = gcd(a, b), a s + b t = d
995
996 void
997 xgcd_left (Fp_polynomial& d, Fp_polynomial& s, const Fp_polynomial& a,
998 const Fp_polynomial& b);
999 // d = gcd(a, b), a s + b ?? = d
1000
1001 inline void
xgcd_right(Fp_polynomial & d,Fp_polynomial & t,const Fp_polynomial & a,const Fp_polynomial & b)1002 xgcd_right (Fp_polynomial& d, Fp_polynomial& t, const Fp_polynomial& a,
1003 const Fp_polynomial& b)
1004 {
1005 xgcd_left(d, t, b, a);
1006 }
1007
1008
1009
1010 // d = gcd(a, b), a ?? + b t = d
1011
1012 void
1013 resultant (bigint &, const Fp_polynomial &, const Fp_polynomial &);
1014
1015
1016 inline bigint
resultant(const Fp_polynomial & a,const Fp_polynomial & b)1017 resultant (const Fp_polynomial& a, const Fp_polynomial& b)
1018 {
1019 bigint x;
1020 resultant (x, a, b);
1021 return x;
1022 }
1023
1024
1025
1026 #ifndef HEADBANGER
1027 // These always use "classical" arithmetic
1028 void
1029 plain_xgcd (Fp_polynomial& d, Fp_polynomial& s, Fp_polynomial& t,
1030 const Fp_polynomial& a, const Fp_polynomial& b);
1031 void
1032 plain_xgcd_left (Fp_polynomial& d, Fp_polynomial& s,
1033 const Fp_polynomial& a, const Fp_polynomial& b);
1034 void
1035 plain_gcd (Fp_polynomial& x, const Fp_polynomial& a, const Fp_polynomial& b);
1036 #endif // HEADBANGER
1037
1038
1039 //***************************************************************
1040 //
1041 // Modular Arithmetic without pre-conditioning
1042 //
1043 //***************************************************************
1044
1045 #ifndef HEADBANGER
1046 // arithmetic mod f.
1047 // all inputs and outputs are polynomials of degree less than deg(f).
1048 // ASSUMPTION: f is assumed monic, and deg(f) > 0.
1049 // NOTE: if you want to do many computations with a fixed f,
1050 // use the Fp_poly_modulus data structure and associated routines below.
1051
1052 void
1053 multiply_mod (Fp_polynomial& c, const Fp_polynomial& a, const Fp_polynomial& b, const Fp_polynomial& f);
1054 // c = (a * b) % f
1055
1056 void
1057 square_mod (Fp_polynomial& c, const Fp_polynomial& a, const Fp_polynomial& f);
1058 // c = a^2 % f
1059
1060 void
1061 multiply_by_x_mod (Fp_polynomial& c, const Fp_polynomial& a, const Fp_polynomial& f);
1062 // c = (a * x) mod f
1063
1064 void
1065 invert_mod (Fp_polynomial& c, const Fp_polynomial& a, const Fp_polynomial& f);
1066 // c = a^{-1} % f, error if a is not invertible
1067
1068 bool
1069 invert_mod_status(Fp_polynomial& c, const Fp_polynomial& a, const Fp_polynomial& f);
1070 // if (a, f) = 1, returns 1 and sets c = a^{-1} % f
1071 // otherwise, returns 0 and sets c = (a, f)
1072
1073 void
1074 power_mod (Fp_polynomial& c, const Fp_polynomial& a, const bigint& e, const Fp_polynomial& f);
1075 // c = a^e % f
1076 // WARNING: obsolete. Use power_mod with Fp_poly_modulus (see below).
1077
1078 void
1079 power_x_mod (Fp_polynomial& c, const bigint& e, const Fp_polynomial& f);
1080 //c = x^e mod f
1081 // WARNING: obsolete. Use power_mod with Fp_poly_modulus (see below).
1082
1083 void
1084 power_x_plus_a_mod (Fp_polynomial& c, const bigint& a, const bigint& e, const Fp_polynomial& f);
1085 // c = (x + a)^e mod f
1086 // WARNING: obsolete. Use power_mod with Fp_poly_modulus (see below).
1087
1088 #endif // HEADBANGER
1089
1090
1091
1092 //***************************************************************
1093 //
1094 // I/O
1095 //
1096 //***************************************************************
1097
1098 std::istream& operator >> (std::istream& s, Fp_polynomial& x);
1099 std::ostream& operator << (std::ostream& s, const Fp_polynomial& a);
1100
1101
1102
1103 //***************************************************************
1104 //
1105 // Miscellaneaous
1106 //
1107 //***************************************************************
1108
1109 void randomize (Fp_polynomial& x, const bigint& p, lidia_size_t n);
1110 // generate a random polynomial of degree = n with coefficients in Z/pZ
1111
1112 void power (Fp_polynomial &x, const Fp_polynomial &a, lidia_size_t e);
1113 // x = a^e, e >= 0
1114
1115
1116
1117 //*********************************************************************
1118 //
1119 // miscellaneous
1120 //
1121 //*********************************************************************
1122
1123
1124 // ************ min_poly, checked_min_poly, irred_poly ************
1125
1126 void
1127 prob_min_poly (Fp_polynomial& h, const Fp_polynomial& g, lidia_size_t m,
1128 const Fp_poly_modulus& F);
1129 // computes the monic minimal polynomial of (g mod f).
1130 // m = a bound on the degree of the minimal polynomial.
1131 // The algorithm is probabilistic, always returns a divisor of
1132 // the minimal polynomial, and returns a proper divisor with
1133 // probability at most m/p.
1134
1135
1136 void min_poly (Fp_polynomial& h, const Fp_polynomial& g, lidia_size_t m,
1137 const Fp_poly_modulus& F);
1138 // same as above, but guarantees that result is correct
1139
1140
1141 void irred_poly (Fp_polynomial& h, const Fp_polynomial& g, lidia_size_t m,
1142 const Fp_poly_modulus& F);
1143 // same as above, but assumes that f is irreducible
1144 // (or at least that the minimal poly of g is itself irreducible).
1145 // The algorithm is deterministic (and hence is always correct).
1146
1147
1148
1149 // ************ Multiplication of several polynomials ************
1150 //void multiply(Fp_polynomial& x, base_vector < Fp_polynomial > & a);
1151 // x = product of all a[i]'s, contents of a[i] are destroyed
1152
1153
1154
1155 // ************ Modular Composition ************
1156 // algorithms for computing g(h) mod f
1157
1158 void compose (Fp_polynomial& x, const Fp_polynomial& g,
1159 const Fp_polynomial& h, const Fp_poly_modulus& F);
1160 // x = g(h) mod f
1161
1162
1163 void compose2 (Fp_polynomial& x1, Fp_polynomial& x2,
1164 const Fp_polynomial& g1, const Fp_polynomial& g2,
1165 const Fp_polynomial& h, const Fp_poly_modulus& F);
1166 // xi = gi(h) mod f (i = 1, 2)
1167 // ALIAS RESTRICTION: xi may not alias gj, for i != j
1168
1169
1170 void compose3 (Fp_polynomial& x1, Fp_polynomial& x2, Fp_polynomial& x3,
1171 const Fp_polynomial& g1, const Fp_polynomial& g2, const Fp_polynomial& g3,
1172 const Fp_polynomial& h, const Fp_poly_modulus& F);
1173 // xi = gi(h) mod f (i = 1..3)
1174 // ALIAS RESTRICTION: xi may not alias gj, for i != j
1175
1176
1177
1178 // ************ update_map ************
1179 void update_map (base_vector< bigint > & x, const base_vector< bigint > & a,
1180 const Fp_poly_multiplier& B, const Fp_poly_modulus& F);
1181 // computes (a, b), (a, (b*X)%f), ..., (a, (b*X^{n-1})%f),
1182 // where (,) denotes the vector inner product.
1183 // This is really a "transposed" MulMod by B.
1184
1185 void update_map (base_vector< bigint > & xx, const base_vector< bigint > & a,
1186 const Fp_polynomial& b, const Fp_polynomial& f);
1187 // same as above, but uses only classical arithmetic
1188
1189
1190
1191 // ************ inner_product ************
1192 void inner_product (bigint& x, const base_vector< bigint > & a,
1193 const Fp_polynomial &b, lidia_size_t offset = 0);
1194
1195
1196
1197 //*********************************************************************
1198 //
1199 // factorization.cc
1200 //
1201 //*********************************************************************
1202
1203 void rec_find_roots (base_vector< bigint > & x, const Fp_polynomial& f);
1204
1205 base_vector< bigint > find_roots (const Fp_polynomial & f, int flag = 0);
1206 // returns the list of roots of f (without multiplicities)
1207 // if (flag != 0), f must be monic and the product of deg(f) distinct
1208 // linear factors
1209 // otherwise, no assumptions on f are made
1210
1211 bigint find_root (const Fp_polynomial & ff);
1212 // finds a single root of ff.
1213 // assumes that ff is monic and splits into distinct linear factors
1214
1215 Fp_polynomial find_factor( const Fp_polynomial & ff, lidia_size_t d );
1216 // Finds an irreducible factor of ff of degree d
1217 // Conditions: ff splits into distinct irreducible factors of degree d.
1218
1219 Fp_polynomial find_factor_degree_d( const Fp_polynomial & ff, lidia_size_t d );
1220 // Finds a factor of ff of degree in [d, 2d]
1221 // Conditions: ff splits into distinct irreducible factors of degree 1.
1222
1223 bool prob_irred_test (const Fp_polynomial & f, lidia_size_t iter = 1);
1224 // performs a fast, probabilistic irreduciblity test
1225 // the test can err only if f is reducible, and the
1226 // error probability is bounded by p^{-iter}.
1227 // Works for any p.
1228
1229
1230 bool det_irred_test (const Fp_polynomial & f);
1231 // performs a recursive deterministic irreducibility test
1232
1233
1234 bool iter_irred_test (const Fp_polynomial & f);
1235 // performs an iterative deterministic irreducibility test,
1236 // based on DDF
1237
1238
1239
1240 void build_irred (Fp_polynomial& f, const bigint & p, lidia_size_t n);
1241 // Build a monic irreducible poly of degree n.
1242
1243
1244 void build_random_irred (Fp_polynomial & f, const Fp_polynomial & g);
1245 // g is a monic irreducible polynomial.
1246 // constructs a random monic irreducible polynomial f of the same degree.
1247
1248 lidia_size_t compute_degree (const Fp_polynomial & h, const Fp_poly_modulus & F,
1249 lidia_size_t d = -1);
1250 // f = F.f is assumed to be an "equal degree" polynomial
1251 // h = x^p mod f
1252 // the common degree of the irreducible factors of f is computed
1253 // d is multiple of common degree, if -1, choose degree F.f
1254 // This routine is useful in counting points on elliptic curves
1255
1256
1257 lidia_size_t prob_compute_degree (const Fp_polynomial & h, const Fp_poly_modulus & F);
1258 // same as above, but uses a slightly faster probabilistic algorithm
1259 // the return value may be 0 or may be too big, but for large p
1260 // (relative to n), this happens with very low probability.
1261
1262
1263
1264 void trace_map (Fp_polynomial & w, const Fp_polynomial & a, lidia_size_t d,
1265 const Fp_poly_modulus & F, const Fp_polynomial & b);
1266 // w = a+a^q+...+^{q^{d-1}} mod f;
1267 // it is assumed that d >= 0, and b = x^q mod f, q a power of p
1268 // Space allocation can be controlled via ComposeBound
1269
1270
1271
1272 void power_compose (Fp_polynomial& w, const Fp_polynomial& b, lidia_size_t d,
1273 const Fp_poly_modulus& F);
1274 // w = x^{q^d} mod f;
1275 // it is assumed that d >= 0, and b = x^q mod f, q a power of p
1276 // Space allocation can be controlled via ComposeBound
1277
1278
1279
1280 #ifdef LIDIA_NAMESPACE
1281 } // end of namespace LiDIA
1282 # undef IN_NAMESPACE_LIDIA
1283 #endif
1284
1285
1286
1287 #define LIDIA_CLASS_FP_POLYNOMIAL
1288
1289 #include "LiDIA/specialization/Fp_polynomial.special"
1290
1291
1292
1293 #endif // LIDIA_FP_POLYNOMIAL_H_GUARD_
1294