Line data Source code
1 : /* Median/middle product.
2 :
3 : Copyright 2003, 2004, 2005, 2006, 2007, 2008 Laurent Fousse, Paul Zimmermann,
4 : Alexander Kruppa, Dave Newman.
5 :
6 : This file is part of the ECM Library.
7 :
8 : The ECM Library is free software; you can redistribute it and/or modify
9 : it under the terms of the GNU Lesser General Public License as published by
10 : the Free Software Foundation; either version 3 of the License, or (at your
11 : option) any later version.
12 :
13 : The ECM Library is distributed in the hope that it will be useful, but
14 : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 : or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 : License for more details.
17 :
18 : You should have received a copy of the GNU Lesser General Public License
19 : along with the ECM Library; see the file COPYING.LIB. If not, see
20 : http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 : 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22 :
23 : /* Reference:
24 : [1] Tellegen's Principle into Practice, by A. Bostan, G. Lecerf and E. Schost,
25 : Proc. of ISSAC'03, Philadelphia, 2003.
26 : */
27 :
28 : #include <stdio.h>
29 : #include "ecm-impl.h"
30 :
31 : #ifndef MAX
32 : #define MAX(a,b) (((a) > (b)) ? (a) : (b))
33 : #endif
34 :
35 : #ifndef MIN
36 : #define MIN(a,b) (((a) < (b)) ? (a) : (b))
37 : #endif
38 :
39 : extern unsigned int Fermat;
40 :
41 : static void list_add_wrapper (listz_t, listz_t, listz_t, unsigned int,
42 : unsigned int);
43 : static void list_sub_wrapper (listz_t, listz_t, listz_t, unsigned int,
44 : unsigned int);
45 : static unsigned int TKarMul (listz_t, unsigned int, listz_t, unsigned int,
46 : listz_t, unsigned int, listz_t);
47 : static void list_sub_safe (listz_t, listz_t, listz_t, unsigned int,
48 : unsigned int, unsigned int);
49 : static void list_add_safe (listz_t, listz_t, listz_t, unsigned int,
50 : unsigned int, unsigned int);
51 : static unsigned int TToomCookMul (listz_t, unsigned int, listz_t, unsigned int,
52 : listz_t, unsigned int, listz_t);
53 : static unsigned int TToomCookMul_space (unsigned int, unsigned int,
54 : unsigned int);
55 :
56 : static void
57 22938804 : list_add_wrapper (listz_t p, listz_t q, listz_t r, unsigned int n,
58 : unsigned int max_r)
59 : {
60 22938804 : list_add (p, q, r, MIN (n, max_r));
61 22938804 : if (n > max_r)
62 431674 : list_set (p + max_r, q + max_r, n - max_r);
63 22938804 : }
64 :
65 : static void
66 22938804 : list_sub_wrapper (listz_t p, listz_t q, listz_t r, unsigned int n,
67 : unsigned int max_r)
68 : {
69 22938804 : list_sub (p, q, r, MIN (n, max_r));
70 22938804 : if (n > max_r)
71 12185 : list_set (p + max_r, q + max_r, n - max_r);
72 22938804 : }
73 :
74 : /* Given a[0..m] and c[0..l], puts in b[0..n] the coefficients
75 : of degree m to n+m of rev(a)*c, i.e.
76 : b[0] = a[0]*c[0] + ... + a[i]*c[i] with i = min(m, l)
77 : ...
78 : b[k] = a[0]*c[k] + ... + a[i]*c[i+k] with i = min(m, l-k)
79 : ...
80 : b[n] = a[0]*c[n] + ... + a[i]*c[i+n] with i = min(m, l-n) [=l-n].
81 : Using auxiliary memory in t.
82 : Implements algorithm TKarMul of [1].
83 : Assumes deg(c) = l <= m+n.
84 : */
85 :
86 : static unsigned int
87 111101393 : TKarMul (listz_t b, unsigned int n,
88 : listz_t a, unsigned int m, listz_t c, unsigned int l, listz_t t)
89 : {
90 : unsigned int k, mu, nu, h;
91 : unsigned int s1;
92 111101393 : unsigned tot_muls = 0;
93 : #ifdef DEBUG
94 : fprintf (ECM_STDOUT, "Enter TKarMul.\nm = %d\nn = %d\nl = %d\n", m, n, l);
95 : fprintf (ECM_STDOUT, "a = ");
96 : print_list (a, m + 1);
97 : fprintf (ECM_STDOUT, "\nc = ");
98 : print_list (c, l + 1);
99 : fprintf (ECM_STDOUT, "\n");
100 : #endif
101 :
102 :
103 111101393 : if (n == 0)
104 : {
105 : #ifdef DEBUG
106 : fprintf (ECM_STDOUT, "Case n = 0.\n");
107 : #endif
108 87486333 : mpz_mul (b[0], a[0], c[0]);
109 88137197 : for (k = 1; (k <= m) && (k <= l); k++)
110 650864 : mpz_addmul (b[0], a[k], c[k]);
111 : #ifdef DEBUG
112 : fprintf (ECM_STDOUT, "Exit TKarMul.\n");
113 : #endif
114 87486333 : return MIN (m, l) + 1;
115 : }
116 :
117 23615060 : if (m == 0)
118 : {
119 : #ifdef DEBUG
120 : fprintf (ECM_STDOUT, "Case m = 0.\n");
121 : #endif
122 717314 : for (k = 0; (k <= l) && (k <= n); k++)
123 500049 : mpz_mul (b[k], a[0], c[k]);
124 217265 : for (k = l + 1; k <= n; k++)
125 0 : mpz_set_ui (b[k], 0);
126 : #ifdef DEBUG
127 : fprintf (ECM_STDOUT, "Exit TKarMul.\n");
128 : #endif
129 217265 : return MIN (n, l) + 1;
130 : }
131 :
132 23397795 : mu = (m / 2) + 1; /* 1 <= mu <= m */
133 23397795 : nu = (n / 2) + 1; /* 1 <= nu <= n */
134 23397795 : h = MAX (mu, nu); /* h >= 1 */
135 :
136 : #ifdef DEBUG
137 : fprintf (ECM_STDOUT, "mu = %d\nnu = %d\nh = %d\n", mu, nu, h);
138 : #endif
139 :
140 23397795 : if (mu > n)
141 : {
142 : #ifdef DEBUG
143 : fprintf (ECM_STDOUT, "Case mu > n.\n");
144 : #endif
145 :
146 137353 : tot_muls += TKarMul (b, n, a, mu - 1, c, l, t);
147 137353 : if (l >= mu)
148 : {
149 : /* we have to check l-mu <= n + (m-mu), i.e. l <= n+m */
150 137353 : tot_muls += TKarMul (t, n, a + mu, m - mu, c + mu, l - mu, t + n + 1);
151 137353 : list_add (b, b, t, n + 1);
152 : }
153 : #ifdef DEBUG
154 : fprintf (ECM_STDOUT, "Exit TKarMul.\n");
155 : #endif
156 137353 : return tot_muls;
157 : }
158 :
159 23260442 : if (nu > m)
160 : {
161 : #ifdef DEBUG
162 : fprintf (ECM_STDOUT, "Case nu > m.\n");
163 : #endif
164 :
165 : /* we have to check MIN(l,m+nu-1) <= nu-1+m: trivial */
166 321638 : tot_muls += TKarMul (b, nu - 1, a, m, c, MIN (l, m + nu - 1), t);
167 :
168 : /* Description broken in reference. Should be a list
169 : * concatenation, not an addition.
170 : * Fixed now.
171 : */
172 :
173 321638 : if (l >= nu)
174 : {
175 : /* we have to check l-nu <= n-nu+m, i.e. l <= n+m: trivial */
176 321638 : tot_muls += TKarMul (b + nu, n - nu, a, m, c + nu, l - nu, t);
177 : }
178 : else
179 0 : list_zero (b + nu, n - nu + 1);
180 : #ifdef DEBUG
181 : fprintf (ECM_STDOUT, "Exit TKarMul.\n");
182 : #endif
183 321638 : return tot_muls;
184 : }
185 :
186 : /* We want nu = mu */
187 :
188 22938804 : mu = nu = h;
189 :
190 : #ifdef DEBUG
191 : fprintf (ECM_STDOUT, "Base Case.\n");
192 : #endif
193 :
194 22938804 : s1 = MIN (l + 1, n + mu);
195 22938804 : if (l + 1 > nu)
196 22938804 : list_sub_wrapper (t, c, c + nu, s1, l - nu + 1);
197 : else
198 0 : list_set (t, c, s1);
199 : #ifdef DEBUG
200 : fprintf (ECM_STDOUT, "DEBUG c - c[nu].\n");
201 : print_list (t, s1);
202 : fprintf (ECM_STDOUT, "We compute (1) - (3)\n");
203 : #endif
204 22938804 : tot_muls += TKarMul (b, nu - 1, a, mu - 1, t, s1 - 1, t + s1);
205 : /* (1) - (3) */
206 : #ifdef DEBUG
207 : print_list (b, nu);
208 : fprintf (ECM_STDOUT, "We compute (2) - (4)\n");
209 : #endif
210 22938804 : if (s1 >= nu + 1) { /* nu - 1 */
211 22938804 : tot_muls += TKarMul (b + nu, n - nu, a + mu, m - mu,
212 22938804 : t + nu, s1 - nu - 1, t + s1);
213 : /* (2) - (4) */
214 : }
215 : else {
216 0 : list_zero (b + nu, n - nu + 1);
217 : }
218 : #ifdef DEBUG
219 : print_list (b + nu, n - nu + 1);
220 : #endif
221 22938804 : list_add_wrapper (t, a, a + mu, mu, m + 1 - mu);
222 : #ifdef DEBUG
223 : fprintf (ECM_STDOUT, "We compute (2) + (3)\n");
224 : #endif
225 22938804 : if (l >= nu) {
226 22938804 : tot_muls += TKarMul (t + mu, nu - 1, t, mu - 1, c + nu, l - nu,
227 22938804 : t + mu + nu);
228 : }
229 : else
230 0 : list_zero (t + mu, nu);
231 : /* (2) + (3) */
232 : #ifdef DEBUG
233 : print_list (t + mu, nu);
234 : #endif
235 22938804 : list_add (b, b, t + mu, nu);
236 22938804 : list_sub (b + nu, t + mu, b + nu, n - nu + 1);
237 22938804 : return tot_muls;
238 : }
239 :
240 : /* Computes the space needed for TKarMul of b[0..n],
241 : * a[0..m] and c[0..l]
242 : */
243 :
244 : static unsigned int
245 1388269 : TKarMul_space (unsigned int n, unsigned int m, unsigned int l)
246 : {
247 : unsigned int mu, nu, h;
248 : unsigned int s1;
249 :
250 : unsigned int r1, r2;
251 :
252 1388269 : if (n == 0)
253 870686 : return 0;
254 :
255 517583 : if (m == 0)
256 60036 : return 0;
257 :
258 457547 : mu = (m / 2) + 1;
259 457547 : nu = (n / 2) + 1;
260 457547 : h = MAX (mu, nu);
261 :
262 457547 : if (mu > n)
263 : {
264 60187 : r1 = TKarMul_space (n, mu - 1, l);
265 60187 : if (l >= mu)
266 : {
267 59048 : r2 = TKarMul_space (n, m - mu, l - mu) + n + 1;
268 59048 : r1 = MAX (r1, r2);
269 : }
270 60187 : return r1;
271 : }
272 :
273 397360 : if (nu > m)
274 : {
275 68659 : r1 = TKarMul_space (nu - 1, m, MIN (l, m + nu - 1));
276 :
277 68659 : if (l >= nu)
278 : {
279 67923 : r2 = TKarMul_space (n - nu, m,l - nu);
280 67923 : r1 = MAX (r1, r2);
281 : }
282 68659 : return r1;
283 : }
284 :
285 328701 : mu = nu = h;
286 :
287 328701 : s1 = MIN (l + 1, n + mu);
288 328701 : r1 = TKarMul_space (nu - 1, mu - 1, s1 - 1) + s1;
289 328701 : if (s1 >= nu + 1) {
290 325972 : r2 = TKarMul_space (n - nu, m - mu, s1 - nu - 1) + s1;
291 325972 : r1 = MAX (r1, r2);
292 : }
293 328701 : if (l >= nu) {
294 325972 : r2 = TKarMul_space (nu - 1, mu - 1, l - nu) + mu + nu;
295 325972 : r1 = MAX (r1, r2);
296 : }
297 328701 : return r1;
298 : }
299 :
300 : /* list_sub with bound checking
301 : */
302 :
303 : static void
304 37891204 : list_sub_safe (listz_t ret, listz_t a, listz_t b,
305 : unsigned int sizea, unsigned int sizeb,
306 : unsigned int needed)
307 : {
308 : unsigned int i;
309 : unsigned int safe;
310 37891204 : safe = MIN(sizea, sizeb);
311 37891204 : safe = MIN(safe, needed);
312 :
313 37891204 : list_sub (ret, a, b, safe);
314 :
315 37891204 : i = safe;
316 48793495 : while (i < needed)
317 : {
318 10902291 : if (i < sizea)
319 : {
320 464424 : if (i < sizeb)
321 0 : mpz_sub (ret[i], a[i], b[i]);
322 : else
323 464424 : mpz_set (ret[i], a[i]);
324 : }
325 : else
326 : {
327 10437867 : if (i < sizeb)
328 10437867 : mpz_neg (ret[i], b[i]);
329 : else
330 0 : mpz_set_ui (ret[i], 0);
331 : }
332 10902291 : i++;
333 : }
334 37891204 : }
335 :
336 : /* list_add with bound checking
337 : */
338 :
339 : static void
340 9472801 : list_add_safe (listz_t ret, listz_t a, listz_t b,
341 : unsigned int sizea, unsigned int sizeb,
342 : unsigned int needed)
343 : {
344 : unsigned int i;
345 : unsigned int safe;
346 9472801 : safe = MIN(sizea, sizeb);
347 9472801 : safe = MIN(safe, needed);
348 :
349 9472801 : list_add (ret, a, b, i = safe);
350 :
351 9472801 : while (i < needed)
352 : {
353 0 : if (i < sizea)
354 : {
355 0 : if (i < sizeb)
356 0 : mpz_add (ret[i], a[i], b[i]);
357 : else
358 0 : mpz_set (ret[i], a[i]);
359 : }
360 : else
361 : {
362 0 : if (i < sizeb)
363 0 : mpz_set (ret[i], b[i]);
364 : else
365 0 : mpz_set_ui (ret[i], 0);
366 : }
367 0 : i++;
368 : }
369 9472801 : }
370 :
371 : static unsigned int
372 51179360 : TToomCookMul (listz_t b, unsigned int n,
373 : listz_t a, unsigned int m, listz_t c, unsigned int l,
374 : listz_t tmp)
375 : {
376 : unsigned int nu, mu, h;
377 : unsigned int i;
378 : unsigned int btmp;
379 51179360 : unsigned int tot_muls = 0;
380 :
381 51179360 : nu = n / 3 + 1;
382 51179360 : mu = m / 3 + 1;
383 :
384 : /* ensures n + 1 > 2 * nu */
385 51179360 : if ((n < 2 * nu) || (m < 2 * mu))
386 : {
387 : #ifdef DEBUG
388 : fprintf (ECM_STDOUT, "Too small operands, calling TKara.\n");
389 : #endif
390 41366999 : return TKarMul (b, n, a, m, c, l, tmp);
391 : }
392 :
393 : /* First strip unnecessary trailing coefficients of c:
394 : */
395 :
396 9812361 : l = MIN(l, n + m);
397 :
398 : /* Now the degenerate cases. We want 2 * nu <= m.
399 : *
400 : */
401 :
402 9812361 : if (m < 2 * nu)
403 : {
404 : #ifdef DEBUG
405 : fprintf (ECM_STDOUT, "Degenerate Case 1.\n");
406 : #endif
407 166221 : tot_muls += TToomCookMul (b, nu - 1, a, m, c, l, tmp);
408 166221 : if (l >= nu)
409 166221 : tot_muls += TToomCookMul (b + nu, nu - 1, a, m,
410 166221 : c + nu, l - nu, tmp);
411 : else
412 0 : list_zero (b + nu, nu);
413 166221 : if (l >= 2 * nu) /* n >= 2 * nu is assured. Hopefully */
414 166221 : tot_muls += TToomCookMul (b + 2 * nu, n - 2 * nu, a, m,
415 166221 : c + 2 * nu, l - 2 * nu, tmp);
416 : else
417 0 : list_zero (b + 2 * nu, n - 2 * nu + 1);
418 166221 : return tot_muls;
419 : }
420 :
421 : /* Second degenerate case. We want 2 * mu <= n.
422 : */
423 :
424 9646140 : if (n < 2 * mu)
425 : {
426 : #ifdef DEBUG
427 : fprintf (ECM_STDOUT, "Degenerate Case 2.\n");
428 : #endif
429 173339 : tot_muls += TToomCookMul (b, n, a, mu - 1, c, l, tmp);
430 173339 : if (l >= mu)
431 : {
432 346678 : tot_muls += TToomCookMul (tmp, n, a + mu, mu - 1,
433 173339 : c + mu, l - mu, tmp + n + 1);
434 173339 : list_add (b, b, tmp, n + 1);
435 : }
436 173339 : if (l >= 2 * mu)
437 : {
438 346678 : tot_muls += TToomCookMul (tmp, n, a + 2 * mu, m - 2 * mu,
439 173339 : c + 2 * mu, l - 2 * mu, tmp + n + 1);
440 173339 : list_add (b, b, tmp, n + 1);
441 : }
442 173339 : return tot_muls;
443 : }
444 :
445 : #ifdef DEBUG
446 : fprintf (ECM_STDOUT, "Base Case.\n");
447 : fprintf (ECM_STDOUT, "a = ");
448 : print_list (a, m + 1);
449 :
450 : fprintf (ECM_STDOUT, "\nc = ");
451 : print_list (c, l + 1);
452 : #endif
453 9472801 : h = MAX(nu, mu);
454 9472801 : nu = mu = h;
455 :
456 9472801 : list_sub_safe (tmp, c + 3 * h, c + h,
457 9472801 : (l + 1 > 3 * h ? l + 1 - 3 * h : 0),
458 9472801 : (l + 1 > h ? l + 1 - h : 0), 2 * h - 1);
459 9472801 : list_sub_safe (tmp + 2 * h - 1, c, c + 2 * h,
460 9472801 : l + 1, (l + 1 > 2 * h ? l + 1 - 2 * h : 0),
461 9472801 : 2 * h - 1);
462 46867054 : for (i = 0; i < 2 * h - 1; i++)
463 37394253 : mpz_mul_2exp (tmp[2 * h - 1 + i], tmp[2 * h - 1 + i], 1);
464 :
465 : #ifdef DEBUG
466 : print_list (tmp, 4 * h - 2);
467 : #endif
468 :
469 : /* --------------------------------
470 : * | 0 .. 2*h-2 | 2*h-1 .. 4*h-3 |
471 : * --------------------------------
472 : * | c3 - c1 | 2(c0 - c2) |
473 : * --------------------------------
474 : */
475 :
476 9472801 : list_add (tmp + 2 * h - 1, tmp + 2 * h - 1, tmp, 2 * h - 1);
477 :
478 18945602 : tot_muls += TToomCookMul (b, h - 1, a, h - 1, tmp + 2 * h - 1,
479 9472801 : 2 * h - 2, tmp + 4 * h - 2);
480 :
481 : /* b[0 .. h - 1] = 2 * m0 */
482 :
483 : #ifdef DEBUG
484 : fprintf (ECM_STDOUT, "2 * m0 = ");
485 : print_list (b, h);
486 : #endif
487 :
488 9472801 : list_add (tmp + 2 * h - 1, a, a + h, h);
489 :
490 9472801 : list_add (tmp + 2 * h - 1, tmp + 2 * h - 1, a + 2 * h,
491 9472801 : MIN(h, m + 1 - 2 * h));
492 :
493 : /* tmp[2*h-1 .. 3*h-2] = a0 + a1 + a2 */
494 :
495 : #ifdef DEBUG
496 : fprintf (ECM_STDOUT, "\na0 + a1 + a2 = ");
497 : print_list (tmp + 2 * h - 1, h);
498 : #endif
499 :
500 9472801 : list_sub_safe (tmp + 3 * h - 1, c + 2 * h, c + 3 * h,
501 9472801 : (l + 1 > 2 * h ? l + 1 - 2 * h : 0),
502 9472801 : (l + 1 > 3 * h ? l + 1 - 3 * h : 0),
503 9472801 : 2 * h - 1);
504 :
505 : /* -------------------------------------------------
506 : * | 0 .. 2*h-2 | 2*h-1 .. 3*h-2 | 3*h-1 .. 5*h-3 |
507 : * -------------------------------------------------
508 : * | c3 - c1 | a0 + a1 + a2 | c2 - c3 |
509 : * -------------------------------------------------
510 : */
511 :
512 9472801 : btmp = (l + 1 > h ? l + 1 - h : 0);
513 9472801 : btmp = MIN(btmp, 2 * h - 1);
514 46867054 : for (i = 0; i < btmp; i++)
515 : {
516 37394253 : mpz_mul_2exp (tmp[5 * h - 2 + i], c[h + i], 1);
517 37394253 : mpz_add (tmp[5 * h - 2 + i], tmp[5 * h - 2 + i], tmp[3 * h - 1 + i]);
518 : }
519 9472801 : while (i < 2 * h - 1)
520 : {
521 0 : mpz_set (tmp[5 * h - 2 + i], tmp[3 * h - 1 + i]);
522 0 : i++;
523 : }
524 :
525 18945602 : tot_muls += TToomCookMul (b + h, h - 1, tmp + 2 * h - 1, h - 1,
526 9472801 : tmp + 5 * h - 2, 2 * h - 2,
527 9472801 : tmp + 7 * h - 3);
528 :
529 : /* b[h .. 2 * h - 1] = 2 * m1 */
530 : #ifdef DEBUG
531 : fprintf (ECM_STDOUT, "\n2 * m1 = ");
532 : print_list (b + h, h);
533 : #endif
534 :
535 : /* ------------------------------------------------------------------
536 : * | 0 .. 2*h-2 | 2*h-1 .. 3*h-2 | 3*h-1 .. 5*h-3 | 5*h-2 .. 7*h-4 |
537 : * ------------------------------------------------------------------
538 : * | c3 - c1 | a0 + a1 + a2 | c2 - c3 | c2 - c3 + 2c1 |
539 : * ------------------------------------------------------------------
540 : */
541 :
542 :
543 32906328 : for (i = 0; i < h; i++)
544 : {
545 23433527 : mpz_add (tmp[2 * h - 1 + i], tmp[2 * h - 1 + i], a[i + h]);
546 23433527 : if (2 * h + i <= m)
547 18454047 : mpz_addmul_ui (tmp[2 * h - 1 + i], a[2 * h + i], 3);
548 : }
549 18945602 : tot_muls += TToomCookMul (tmp + 5 * h - 2, h - 1,
550 9472801 : tmp + 2 * h - 1, h - 1,
551 9472801 : tmp, 2 * h - 2, tmp + 6 * h - 2);
552 :
553 : /* tmp[5*h-2 .. 6*h - 3] = 6 * m2 */
554 :
555 : #ifdef DEBUG
556 : fprintf (ECM_STDOUT, "\n6 * m2 = ");
557 : print_list (tmp + 5 * h - 2, h);
558 : #endif
559 32906328 : for (i = 0; i < h; i++)
560 : {
561 23433527 : mpz_sub (tmp[2 * h - 1 + i], a[i], a[h + i]);
562 23433527 : if (i + 2 * h <= m)
563 18454047 : mpz_add (tmp[2 * h - 1 + i], tmp[2 * h - 1 + i], a[2 * h + i]);
564 : }
565 :
566 46867054 : for (i = 0; i < 2 * h - 1; i++)
567 : {
568 37394253 : mpz_mul_ui (tmp[3 * h - 1 + i], tmp[3 * h - 1 + i], 3);
569 37394253 : mpz_mul_2exp (tmp[i], tmp[i], 1);
570 : }
571 :
572 9472801 : list_add (tmp + 3 * h - 1, tmp + 3 * h - 1, tmp, 2 * h - 1);
573 :
574 18945602 : tot_muls += TToomCookMul (tmp + 6 * h - 2, h - 1,
575 9472801 : tmp + 2 * h - 1, h - 1,
576 9472801 : tmp + 3 * h - 1, 2 * h - 2,
577 9472801 : tmp + 7 * h - 2);
578 :
579 : /* tmp[6h-2 .. 7h - 3] = 6 * mm1 */
580 :
581 : #ifdef DEBUG
582 : fprintf (ECM_STDOUT, "\n6 * mm1 = ");
583 : print_list (tmp + 6 * h - 2, h);
584 : #endif
585 9472801 : list_add_safe (tmp, tmp, c + 2 * h,
586 : 2 * h,
587 9472801 : (l + 1 > 2 * h ? l + 1 - 2 * h : 0),
588 9472801 : 2 * h - 1);
589 :
590 9472801 : list_sub_safe (tmp, c + 4 * h, tmp,
591 9472801 : (l + 1 > 4 * h ? l + 1 - 4 * h : 0),
592 9472801 : 2 * h - 1, 2 * h - 1);
593 :
594 18945602 : tot_muls += TToomCookMul (b + 2 * h, n - 2 * h, a + 2 * h, m - 2 * h,
595 9472801 : tmp, 2 * h - 1, tmp + 7 * h - 2);
596 :
597 : /* b[2 * h .. n] = minf */
598 :
599 : #ifdef DEBUG
600 : fprintf (ECM_STDOUT, "\nminf = ");
601 : print_list (b + 2 * h, n + 1 - 2 * h);
602 : #endif
603 :
604 : /* Layout of b :
605 : * ---------------------------------------
606 : * | 0 ... h-1 | h ... 2*h-1 | 2*h ... n |
607 : * ---------------------------------------
608 : * | 2 * m0 | 2 * m1 | minf |
609 : * ---------------------------------------
610 : *
611 : * Layout of tmp :
612 : * ---------------------------------------------------
613 : * | 0 ... 5*h-1 | 5*h-2 ... 6*h-3 | 6*h-2 ... 7*h-3 |
614 : * ---------------------------------------------------
615 : * | ?????? | 6 * m2 | 6 * mm1 |
616 : * ---------------------------------------------------
617 : */
618 :
619 9472801 : list_add (tmp, tmp + 5 * h - 2, tmp + 6 * h - 2, h);
620 32906328 : for (i = 0; i < h; i++)
621 23433527 : mpz_divby3_1op (tmp[i]);
622 :
623 : /* t1 = 2 (m2 + mm1)
624 : * tmp[0 .. h - 1] = t1
625 : */
626 :
627 9472801 : list_add (b, b, b + h, h);
628 9472801 : list_add (b, b, tmp, h);
629 32906328 : for (i = 0; i < h; i++)
630 23433527 : mpz_tdiv_q_2exp (b[i], b[i], 1);
631 :
632 : /* b_{low} should be correct */
633 :
634 9472801 : list_add (tmp + h, b + h, tmp, h);
635 :
636 : /* t2 = t1 + 2 m1
637 : * tmp[h .. 2h - 1] = t2
638 : */
639 :
640 9472801 : list_add (b + h, tmp, tmp + h, h);
641 9472801 : list_sub (b + h, b + h, tmp + 6 * h - 2, h);
642 32906328 : for (i = 0; i < h; i++)
643 23433527 : mpz_tdiv_q_2exp (b[h + i], b[h + i], 1);
644 :
645 : /* b_{mid} should be correct */
646 :
647 9472801 : list_add (tmp + h, tmp + h, tmp + 5 * h - 2, n + 1 - 2 * h);
648 27912365 : for (i = 0; i < n + 1 - 2 * h; i++)
649 18439564 : mpz_tdiv_q_2exp (tmp[h + i], tmp[h + i], 1);
650 :
651 9472801 : list_add (b + 2 * h, b + 2 * h, tmp + h, n + 1 - 2 * h);
652 : /* b_{high} should be correct */
653 :
654 9472801 : return tot_muls;
655 : }
656 :
657 : /* Returns space needed by TToomCookMul */
658 :
659 : unsigned int
660 304951 : TToomCookMul_space (unsigned int n, unsigned int m, unsigned int l)
661 :
662 : {
663 : unsigned int nu, mu, h;
664 : unsigned int stmp1, stmp2;
665 :
666 304951 : nu = n / 3 + 1;
667 304951 : mu = m / 3 + 1;
668 :
669 304951 : stmp1 = stmp2 = 0;
670 :
671 : /* ensures n + 1 > 2 * nu */
672 304951 : if ((n < 2 * nu) || (m < 2 * mu))
673 151807 : return TKarMul_space (n, m, l);
674 :
675 : /* First strip unnecessary trailing coefficients of c:
676 : */
677 :
678 153144 : l = MIN(l, n + m);
679 :
680 : /* Now the degenerate cases. We want 2 * nu < m.
681 : *
682 : */
683 :
684 153144 : if (m <= 2 * nu)
685 : {
686 82321 : stmp1 = TToomCookMul_space (nu - 1, m, l);
687 82321 : if (l >= 2 * nu)
688 78677 : stmp2 = TToomCookMul_space (n - 2 * nu, m, l - 2 * nu);
689 3644 : else if (l >= nu)
690 1568 : stmp2 = TToomCookMul_space (nu - 1, m, l - nu);
691 82321 : return MAX(stmp1, stmp2);
692 : }
693 :
694 : /* Second degenerate case. We want 2 * mu < n.
695 : */
696 :
697 70823 : if (n <= 2 * mu)
698 : {
699 62904 : stmp1 += TToomCookMul_space (n, mu - 1, l);
700 62904 : if (l >= 2 * mu)
701 60130 : stmp2 = TToomCookMul_space (n, m - 2 * mu, l - 2 * mu) + n + 1;
702 2774 : else if (l >= mu)
703 1131 : stmp2 = TToomCookMul_space (n, mu - 1, l - mu) + n + 1;
704 62904 : return MAX(stmp1, stmp2);
705 : }
706 :
707 7919 : h = MAX(nu, mu);
708 :
709 7919 : stmp1 = TToomCookMul_space (h - 1, h - 1, 2 * h - 2);
710 7919 : stmp2 = stmp1 + 7 * h - 2;
711 7919 : stmp1 = stmp1 + 6 * h - 2;
712 7919 : stmp1 = MAX(stmp1, stmp2);
713 7919 : stmp2 = TToomCookMul_space (n - 2 * h, m - 2 * h, 2 * h - 1) + 7*h-2;
714 7919 : return MAX(stmp1, stmp2);
715 : }
716 :
717 : /* Given a[0..m] and c[0..l], puts in b[0..n] the coefficients
718 : of degree m to n+m of rev(a)*c, i.e.
719 : b[0] = a[0]*c[0] + ... + a[i]*c[i] with i = min(m, l)
720 : ...
721 : b[k] = a[0]*c[k] + ... + a[i]*c[i+k] with i = min(m, l-k)
722 : ...
723 : b[n] = a[0]*c[n] + ... + a[i]*c[i+n] with i = min(m, l-n) [=l-n].
724 : Using auxiliary memory in tmp.
725 :
726 : Assumes n <= l.
727 :
728 : Returns number of multiplications if known, 0 if not known,
729 : and -1 for error.
730 : */
731 : int
732 2895765 : TMulGen (listz_t b, unsigned int n, listz_t a, unsigned int m,
733 : listz_t c, unsigned int l, listz_t tmp, mpz_t modulus)
734 : {
735 : ASSERT (n <= l);
736 :
737 2895765 : if (Fermat)
738 : {
739 : unsigned int i;
740 295194 : for (i = l + 1; i > 1 && (i&1) == 0; i >>= 1);
741 : ASSERT(i == 1);
742 : ASSERT(n + 1 == (l + 1) / 2);
743 : ASSERT(m == l - n || m + 1 == l - n);
744 98424 : return F_mul_trans (b, a, c, m + 1, l + 1, Fermat, tmp);
745 : }
746 :
747 2797341 : if ((double) n * (double) mpz_sizeinbase (modulus, 2) >= KS_TMUL_THRESHOLD)
748 : {
749 666 : if (TMulKS (b, n, a, m, c, l, modulus, 1)) /* Non-zero means error */
750 0 : return -1;
751 666 : return 0; /* We have no mul count so we return 0 */
752 : }
753 :
754 2796675 : return TToomCookMul (b, n, a, m, c, l, tmp);
755 : }
756 :
757 :
758 : unsigned int
759 2421 : TMulGen_space (unsigned int n, unsigned int m, unsigned int l)
760 : {
761 2421 : if (Fermat)
762 39 : return 2 * (l + 1);
763 : else
764 2382 : return TToomCookMul_space (n, m, l);
765 : }
|