1 //==============================================================================================
2 //
3 // This file is part of LiDIA --- a library for computational number theory
4 //
5 // Copyright (c) 1994--2001 the LiDIA Group. All rights reserved.
6 //
7 // See http://www.informatik.tu-darmstadt.de/TI/LiDIA/
8 //
9 //----------------------------------------------------------------------------------------------
10 //
11 // $Id$
12 //
13 // Author : Thomas Pfahler (TPf)
14 // mainly based on V. Shoup's Fp_poly_modulus for Fp_polynomial
15 // Changes : See CVS log
16 //
17 //==============================================================================================
18
19
20 #ifdef HAVE_CONFIG_H
21 # include "config.h"
22 #endif
23 #include "LiDIA/gf_polynomial.h"
24
25
26
27 #ifdef LIDIA_NAMESPACE
28 namespace LiDIA {
29 #endif
30
31
32
33 //******************************************************************
34 //
35 // class gf_poly_modulus
36 //
37 //******************************************************************
38
39 // If you need to do a lot of arithmetic modulo a fixed f,
40 // build gf_poly_modulus F for f. This pre-computes information about f
41 // that speeds up the computation a great deal.
42 // f should be monic, and deg(f) > 0.
43 // For mathematical background, see V. Shoup's paper
44
45
gf_poly_modulus(const gf_poly_modulus & p)46 gf_poly_modulus::gf_poly_modulus(const gf_poly_modulus &p) :
47 F_plain(p.F_plain), F(p.F), F_recip(p.F_recip)
48 { }
49
50
51
gf_poly_modulus(const gf_polynomial & f)52 gf_poly_modulus::gf_poly_modulus(const gf_polynomial &f)
53 {
54 build(f);
55 }
56
57
58
build(const gf_polynomial & f)59 void gf_poly_modulus::build(const gf_polynomial &f)
60 {
61 debug_handler("gf_poly_modulus", "build(gf_polynomial&)");
62
63 if (f.degree() <= 0 || !f.is_monic())
64 lidia_error_handler("gf_poly_modulus", "build(gf_polynomial&)::"
65 "argument must be monic and of degree > 0");
66
67 F_plain.assign(f);
68
69 int fd = f.get_field().degree();
70 if (f.degree() == 1 || (f.get_field().characteristic() == 2 &&
71 (fd > 8 || (fd > 1 && f.degree() < 100*fd)))) {
72 use_Fp = false;
73 return;
74 }
75
76 use_Fp = true;
77
78 gf_polynomial P1, P2;
79 copy_reverse(P1, f, 0, f.degree());
80 invert(P2, P1, f.degree());
81 copy_reverse(P1, P2, 0, f.degree()-2);
82
83 to_Kronecker(F_recip, P1, 0, P1.degree());
84 to_Kronecker(F, f, 0, f.degree());
85 }
86
87
88
89 // x = a % f
90 // deg(a) <= 2(n-1), where n = f.degree()
rem21(gf_polynomial & x,const gf_polynomial & a) const91 void gf_poly_modulus::rem21(gf_polynomial& x, const gf_polynomial& a) const
92 {
93 debug_handler("gf_poly_modulus", "rem21(gf_polynomial&, gf_polynomial&)");
94 lidia_size_t da, n;
95 da = a.degree();
96 n = F_plain.degree();
97
98 if (da > 2*n-2)
99 lidia_error_handler("gf_poly_modulus", "rem21(gf_polynomial&, gf_polynomial&)::bad_args");
100
101 if (da < n) {
102 x.assign(a);
103 return;
104 }
105 if (!use_Fp) {
106 remainder(x, a, F_plain);
107 return;
108 }
109
110 Fp_polynomial r1, r2;
111 gf_polynomial P1;
112 P1.ffield = a.ffield;
113 gf_polynomial::build_frame(a.ffield);
114
115 to_Kronecker(r1, a, n, 2*(n-1));
116 multiply(r1, F_recip, r1);
117 from_Kronecker(P1, r1, n-2, 2*n-4);
118
119 to_Kronecker(r1, P1, 0, P1.degree());
120 multiply(r1, F, r1);
121 from_Kronecker(P1, r1, 0, n-1);
122
123 lidia_size_t i, ds = P1.degree();
124 x.ffield = a.ffield;
125 x.set_degree(n-1);
126 for (i = 0; i < n; i++) {
127 if (i <= ds) subtract(x[i], a[i], P1[i]);
128 else x[i].assign(a[i]);
129 }
130 x.remove_leading_zeros();
131 gf_polynomial::delete_frame();
132 }
133
134
135
136 // x = a % f, no restrictions on deg(a); makes repeated calls to rem21
remainder(gf_polynomial & x,const gf_polynomial & a,const gf_poly_modulus & F)137 void remainder(gf_polynomial& x, const gf_polynomial& a, const gf_poly_modulus& F)
138 {
139 debug_handler("gf_poly_modulus", "remainder(gf_polynomial&, gf_polynomial&, gf_poly_modulus&)");
140
141 lidia_size_t da = a.degree();
142 lidia_size_t n = F.modulus().degree();
143
144 if (da <= 2*n-2) {
145 F.rem21(x, a);
146 return;
147 }
148
149 if (!F.use_Fp) {
150 remainder(x, a, F.F_plain);
151 return;
152 }
153
154 gf_polynomial buf;
155 buf.ffield = (gf_polynomial::common_field(a.ffield, F.F_plain.ffield));
156 gf_polynomial::build_frame(buf.ffield);
157
158 lidia_size_t a_len = da+1;
159
160 while (a_len > 0) {
161 lidia_size_t old_buf_len = buf.degree() + 1;
162 lidia_size_t amt = comparator< lidia_size_t >::min(2*n - 1 - old_buf_len, a_len);
163
164 buf.set_degree(old_buf_len+amt-1);
165
166 lidia_size_t i;
167
168 for (i = old_buf_len+amt-1; i >= amt; i--)
169 buf[i].assign(buf[i-amt]);
170
171 for (i = amt-1; i >= 0; i--)
172 buf[i].assign(a[a_len-amt+i]);
173 buf.remove_leading_zeros();
174
175 F.rem21(buf, buf);
176 a_len -= amt;
177 }
178 gf_polynomial::delete_frame();
179 x.assign(buf);
180 }
181
182
183
184 // x = (a * b) % f
185 // deg(a), deg(b) < deg_f
multiply(gf_polynomial & x,const gf_polynomial & a,const gf_polynomial & b,const gf_poly_modulus & F)186 void multiply(gf_polynomial& x, const gf_polynomial& a, const gf_polynomial& b, const gf_poly_modulus& F)
187 {
188 debug_handler("gf_poly_modulus", "multiply(gf_polynomial&, gf_polynomial&, gf_polynomial&, gf_poly_modulus&)");
189
190 lidia_size_t da, db, n;
191
192 da = a.degree();
193 db = b.degree();
194 n = F.modulus().degree();
195
196 if (da >= n || db >= n)
197 lidia_error_handler("gf_poly_modulus", "multiply(gf_polynomial&, gf_polynomial&, gf_polynomial&, gf_poly_modulus&)::degree of gf_polynomials must be < degree of gf_poly_modulus");
198
199 if (!F.use_Fp) {
200 remainder(x, a*b, F.F_plain);
201 return;
202 }
203
204 x.ffield = (gf_polynomial::common_field(
205 gf_polynomial::common_field(a.ffield, b.ffield),
206 F.F_plain.ffield));
207 gf_polynomial::build_frame(x.ffield);
208
209 Fp_polynomial r1, r2, prod;
210 gf_polynomial P1;
211 P1.ffield = x.ffield;
212
213 #if 0
214 //HERE IS THE ORIGINAL AND READABLE CODE ...
215 to_Kronecker(r1, a, 0, a.degree());
216 to_Kronecker(r2, b, 0, b.degree());
217 multiply(prod, r1, r2);
218 from_Kronecker(P1, prod, n, da + db);
219
220 to_Kronecker(r1, P1, 0, P1.degree());
221 multiply(r2, F.F_recip, r1);
222 from_Kronecker(P1, r2, n-2, 2*n-4);
223
224 to_Kronecker(r1, P1, 0, P1.degree());
225 multiply(r2, r1, F.F);
226 from_Kronecker(x, r2, 0, n-1);
227
228 from_Kronecker(P1, prod, 0, n-1);
229 subtract(x, P1, x);
230 gf_polynomial::delete_frame();
231 #endif
232 //NOW, THE OPTIMIZED VERSION (EQUIVALENT TO THE UPPER ONE)
233
234 lidia_size_t N = a.get_field().degree();
235 //cmp. Kronecker substitution:
236 //one coeff. of gf_polynomial equals (2*N-1) coeffs. of Fp_polynomial
237
238 to_Kronecker(r1, a, 0, a.degree());
239 to_Kronecker(r2, b, 0, b.degree());
240 multiply(prod, r1, r2);
241 from_Kronecker(P1, prod, n, da + db);
242 trunc(prod, prod, (2*N-1)*n);
243 //deletes the coefficients from prod which we have already converted
244
245 to_Kronecker(r1, P1, 0, P1.degree());
246 multiply(r2, F.F_recip, r1);
247 from_Kronecker(P1, r2, n-2, 2*n-4);
248
249 to_Kronecker(r1, P1, 0, P1.degree());
250 multiply(r2, r1, F.F);
251
252 //from_Kronecker(x, r2, 0, n-1);
253 //from_Kronecker(P1, prod, 0, n-1);
254 //subtract(x, P1, x);
255
256 trunc(r2, r2, n*(2*N-1));
257 subtract(prod, prod, r2);
258 from_Kronecker(x, prod, 0, n-1);
259
260 gf_polynomial::delete_frame();
261 }
262
263
264
265 // x = a^2 % f a.degree() < f.degree()
square(gf_polynomial & x,const gf_polynomial & a,const gf_poly_modulus & F)266 void square(gf_polynomial& x, const gf_polynomial& a, const gf_poly_modulus& F)
267 {
268 debug_handler("gf_poly_modulus", "square(gf_polynomial&, gf_polynomial&, gf_poly_modulus&)");
269
270
271 lidia_size_t da, n;
272
273 da = a.degree();
274 n = F.F_plain.degree();
275
276 if (da >= n)
277 debug_handler("gf_poly_modulus", "square(gf_polynomial&, gf_polynomial&, gf_poly_modulus&)::degree of gf_polynomial must be < degree of gf_poly_modulus");
278
279 if (!F.use_Fp) {
280 square(x, a);
281 remainder(x, x, F.F_plain);
282 return;
283 }
284
285 //compare :
286 //multiply(gf_polynomial&, gf_polynomial&, gf_polynomial&, gf_poly_modulus&)
287
288 Fp_polynomial r1, r2, prod;
289 gf_polynomial P1;
290 x.ffield = (gf_polynomial::common_field(a.ffield, F.F_plain.ffield));
291 P1.ffield = x.ffield;
292 gf_polynomial::build_frame(x.ffield);
293
294 lidia_size_t N = a.get_field().degree();
295 //cmp. Kronecker substitution:
296 //one coeff. of gf_polynomial equals (2*N-1) coeffs. of Fp_polynomial
297
298 to_Kronecker(r1, a, 0, a.degree());
299 square(prod, r1);
300 from_Kronecker(P1, prod, n, 2*a.degree());
301 trunc(prod, prod, (2*N-1)*n);
302 //deletes the coefficients from prod which we have already converted
303
304 to_Kronecker(r1, P1, 0, P1.degree());
305 multiply(r2, F.F_recip, r1);
306 from_Kronecker(P1, r2, n-2, 2*n-4);
307
308 to_Kronecker(r1, P1, 0, P1.degree());
309 multiply(r2, r1, F.F);
310
311 //from_Kronecker(x, r2, 0, n-1);
312 //from_Kronecker(P1, prod, 0, n-1);
313 //subtract(x, P1, x);
314
315 trunc(r2, r2, n*(2*N-1));
316 subtract(prod, prod, r2);
317 from_Kronecker(x, prod, 0, n-1);
318
319 gf_polynomial::delete_frame();
320
321 }
322
323
324
325 // x = a^e % f
power(gf_polynomial & h,const gf_polynomial & g,const bigint & e,const gf_poly_modulus & F)326 void power(gf_polynomial& h, const gf_polynomial& g, const bigint& e, const gf_poly_modulus& F)
327 {
328 debug_handler("gf_poly_modulus", "power(gf_polynomial&, gf_polynomial&, bigint&, gf_poly_modulus&)");
329
330 gf_polynomial lg(g);
331
332 if (e.is_negative())
333 lidia_error_handler("gf_poly_modulus", "power(gf_polynomial&, gf_polynomial&, bigint&, gf_poly_modulus&)::exponent must be positive");
334
335 int i, n = e.bit_length();
336
337 //### gf_poly_multiplier G(lg, F);
338
339 h.assign_one(gf_polynomial::common_field(g.ffield, F.F_plain.ffield));
340
341 for (i = n - 1; i >= 0; i--) {
342 square(h, h, F);
343 if (e.bit(i))
344 //### multiply(h, h, G, F); //gf_poly_multiplier
345 multiply(h, h, lg, F);
346 }
347
348 }
349
350
351
352 // x = X^e % f
power_x(gf_polynomial & h,const bigint & e,const gf_poly_modulus & F)353 void power_x(gf_polynomial& h, const bigint& e, const gf_poly_modulus& F)
354 {
355 debug_handler("gf_poly_modulus", "power_x(gf_polynomial&, bigint&, gf_poly_modulus&)");
356
357 if (e.is_negative())
358 lidia_error_handler("gf_poly_modulus", "power_x(gf_polynomial&, bigint&, gf_poly_modulus&)::exponent must be positive");
359
360 int i, n = e.bit_length();
361
362 h.assign_one(F.F_plain.get_field());
363
364 for (i = n - 1; i >= 0; i--) {
365 square(h, h, F);
366 if (e.bit(i))
367 multiply_by_x_mod(h, h, F.modulus());
368 }
369
370 }
371
372
373
374 // x = (X + a)^e % f
power_x_plus_a(gf_polynomial & h,const gf_element & a,const bigint & e,const gf_poly_modulus & F)375 void power_x_plus_a(gf_polynomial& h, const gf_element& a, const bigint& e, const gf_poly_modulus& F)
376 {
377 debug_handler("gf_poly_modulus", "power_x_plus_a(gf_polynomial&, gf_element&, bigint&, gf_poly_modulus&)");
378
379 if (e.is_negative())
380 lidia_error_handler("gf_poly_modulus", "power_x_plus_a(gf_polynomial&, gf_element&, bigint&, gf_poly_modulus&)::exponent must be positive");
381
382
383 gf_polynomial t1, t2;
384 gf_element la(a); //allows input to alias output
385
386 int i, n = e.bit_length();
387
388 h.assign_one(gf_polynomial::common_field(a.get_field(),
389 F.F_plain.get_field()));
390
391 for (i = n - 1; i >= 0; i--) {
392 square(h, h, F);
393 if (e.bit(i)) {
394 multiply_by_x_mod(t1, h, F.modulus());
395 multiply(t2, h, la);
396 add(h, t1, t2);
397 }
398 }
399
400 }
401
402
403
404 #if 0
405 //
406 // The idea of gf_poly_multiplier (the analogon to 'poly_multiplier' for
407 // Fp_polynomial) does not speed up the multiplication significantly.
408 // Time is spent (with or without this class) in three multiplications of
409 // Fp_polynomial's. Conversions (see Kronecker substitution) do hardly cost
410 // anything - a this is the only thing we could save.
411 // Bad luck.
412 //
413
414 //******************************************************************
415 //
416 // class gf_poly_multiplier
417 //
418 //******************************************************************
419
420 // If you need to compute a * b % f for a fixed b, but for many a's
421 // (for example, computing powers of b modulo f), it is
422 // more efficient to first build a gf_poly_multiplier B for b,
423 // and then use the routine below.
424 // For mathematical background, see V. Shoup's paper
425
426 gf_poly_multiplier::gf_poly_multiplier() : mod(0)
427 { }
428
429
430
431 gf_poly_multiplier::gf_poly_multiplier(const gf_poly_multiplier &B)
432 {
433 if (mod != 0)
434 build(B.modulus(), *mod);
435 }
436
437
438
439 gf_poly_multiplier::gf_poly_multiplier(const gf_polynomial& b,
440 const gf_poly_modulus& F)
441 {
442 build(b, F);
443 }
444
445
446
447 gf_poly_multiplier::~gf_poly_multiplier()
448 { }
449
450
451
452 void
453 gf_poly_multiplier::build(const gf_polynomial& b, const gf_poly_modulus& F)
454 {
455 mod = &F;
456 remainder(B_plain, b, F);
457 to_Kronecker(B, b, 0, b.degree()); // store b
458
459 lidia_size_t n = F.modulus().degree();
460 Fp_polynomial g;
461 gf_polynomial h;
462 h.ffield = b.ffield;
463 gf_polynomial::build_frame(b.ffield);
464 multiply(g, B, F.F_recip);
465 from_Kronecker(h, g, n-1, 2*n-3);
466 h.delete_frame();
467 to_Kronecker(B_div_F, h, 0, h.degree()); // store (b*F_recip)/X^n
468 }
469
470
471
472 void multiply(gf_polynomial &x, const gf_polynomial& a,
473 gf_poly_multiplier& B, gf_poly_modulus & F)
474 {
475 if (B.mod != &F || B.mod == 0)
476 lidia_error_handler("gf_poly_multiplier", "multiply(...)::"
477 "wrong Fp_poly_modulus, or poly_multiplier not initialized");
478
479 lidia_size_t n = F.modulus().degree();
480 lidia_size_t da = a.degree();
481 if (da >= n)
482 lidia_error_handler("gf_poly_multiplier", "multiply(...)::"
483 "degree of Fp_polynomial must be < degree of Fp_poly_modulus");
484
485 Fp_polynomial r1, r2, r3, r4, r5, r6;
486 gf_polynomial P1;
487 x.ffield = (a.ffield);
488 P1.ffield = (a.ffield);
489 gf_polynomial::build_frame(a.ffield);
490
491 to_Kronecker(r1, a, 0, da);
492 multiply(r2, r1, B.B_div_F);
493 from_Kronecker(P1, r2, n-1, 2*n-3);
494
495 multiply(r3, r1, B.B);
496 to_Kronecker(r4, P1, 0, P1.degree());
497 multiply(r5, r4, F.F);
498
499 subtract(r6, r3, r5);
500 from_Kronecker(x, r6, 0, n-1);
501 P1.delete_frame();
502 x.delete_frame();
503 }
504 #endif
505
506
507
508 #ifdef LIDIA_NAMESPACE
509 } // end of namespace LiDIA
510 #endif
511