← Back to branch summary

~ubuntu-branches/ubuntu/quantal/gclcvs/quantal

« back to all changes in this revision

Viewing changes to gmp3/mpn/generic/mul_n.c

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2004-06-24 15:13:46 UTC
  • Revision ID: james.westby@ubuntu.com-20040624151346-xh0xaaktyyp7aorc
Tags: 2.7.0-26
C_GC_OFFSET is 2 on m68k-linux

Show diffs side-by-side

added added

removed removed

Lines of Context:
gmp3/mpn/generic/mul_n.c
 
1
/* mpn_mul_n and helper function -- Multiply/square natural numbers.
 
2
 
 
3
   THE HELPER FUNCTIONS IN THIS FILE (meaning everything except mpn_mul_n)
 
4
   ARE INTERNAL FUNCTIONS WITH MUTABLE INTERFACES.  IT IS ONLY SAFE TO REACH
 
5
   THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST GUARANTEED
 
6
   THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
 
7
 
 
8
 
 
9
Copyright 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000, 2001 Free Software
 
10
Foundation, Inc.
 
11
 
 
12
This file is part of the GNU MP Library.
 
13
 
 
14
The GNU MP Library is free software; you can redistribute it and/or modify
 
15
it under the terms of the GNU Lesser General Public License as published by
 
16
the Free Software Foundation; either version 2.1 of the License, or (at your
 
17
option) any later version.
 
18
 
 
19
The GNU MP Library is distributed in the hope that it will be useful, but
 
20
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 
21
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 
22
License for more details.
 
23
 
 
24
You should have received a copy of the GNU Lesser General Public License
 
25
along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
 
26
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 
27
MA 02111-1307, USA. */
 
28
 
 
29
#include "gmp.h"
 
30
#include "gmp-impl.h"
 
31
#include "longlong.h"
 
32
 
 
33
 
 
34
#ifndef USE_MORE_MPN
 
35
#if !defined (__alpha) && !defined (__mips)
 
36
/* For all other machines, we want to call mpn functions for the compund
 
37
   operations instead of open-coding them.  */
 
38
#define USE_MORE_MPN 1
 
39
#endif
 
40
#endif
 
41
 
 
42
/*== Function declarations =================================================*/
 
43
 
 
44
static void evaluate3 _PROTO ((mp_ptr, mp_ptr, mp_ptr,
 
45
                               mp_ptr, mp_ptr, mp_ptr,
 
46
                               mp_srcptr, mp_srcptr, mp_srcptr,
 
47
                               mp_size_t, mp_size_t));
 
48
static void interpolate3 _PROTO ((mp_srcptr,
 
49
                                  mp_ptr, mp_ptr, mp_ptr,
 
50
                                  mp_srcptr,
 
51
                                  mp_ptr, mp_ptr, mp_ptr,
 
52
                                  mp_size_t, mp_size_t));
 
53
static mp_limb_t add2Times _PROTO ((mp_ptr, mp_srcptr, mp_srcptr, mp_size_t));
 
54
 
 
55
 
 
56
/*-- mpn_kara_mul_n ---------------------------------------------------------------*/
 
57
 
 
58
/* Multiplies using 3 half-sized mults and so on recursively.
 
59
 * p[0..2*n-1] := product of a[0..n-1] and b[0..n-1].
 
60
 * No overlap of p[...] with a[...] or b[...].
 
61
 * ws is workspace.
 
62
 */
 
63
 
 
64
void
 
65
mpn_kara_mul_n (mp_ptr p, mp_srcptr a, mp_srcptr b, mp_size_t n, mp_ptr ws)
 
66
{
 
67
  mp_limb_t w, w0, w1;
 
68
  mp_size_t n2;
 
69
  mp_srcptr x, y;
 
70
  mp_size_t i;
 
71
  int sign;
 
72
 
 
73
  n2 = n >> 1;
 
74
  ASSERT (n2 > 0);
 
75
 
 
76
  if ((n & 1) != 0)
 
77
    {
 
78
      /* Odd length. */
 
79
      mp_size_t n1, n3, nm1;
 
80
 
 
81
      n3 = n - n2;
 
82
 
 
83
      sign = 0;
 
84
      w = a[n2];
 
85
      if (w != 0)
 
86
        w -= mpn_sub_n (p, a, a + n3, n2);
 
87
      else
 
88
        {
 
89
          i = n2;
 
90
          do
 
91
            {
 
92
              --i;
 
93
              w0 = a[i];
 
94
              w1 = a[n3+i];
 
95
            }
 
96
          while (w0 == w1 && i != 0);
 
97
          if (w0 < w1)
 
98
            {
 
99
              x = a + n3;
 
100
              y = a;
 
101
              sign = ~0;
 
102
            }
 
103
          else
 
104
            {
 
105
              x = a;
 
106
              y = a + n3;
 
107
            }
 
108
          mpn_sub_n (p, x, y, n2);
 
109
        }
 
110
      p[n2] = w;
 
111
 
 
112
      w = b[n2];
 
113
      if (w != 0)
 
114
        w -= mpn_sub_n (p + n3, b, b + n3, n2);
 
115
      else
 
116
        {
 
117
          i = n2;
 
118
          do 
 
119
            {
 
120
              --i;
 
121
              w0 = b[i]; 
 
122
              w1 = b[n3+i];
 
123
            }
 
124
          while (w0 == w1 && i != 0);
 
125
          if (w0 < w1)
 
126
            {
 
127
              x = b + n3;
 
128
              y = b;
 
129
              sign = ~sign;
 
130
            }
 
131
          else
 
132
            {
 
133
              x = b;
 
134
              y = b + n3;
 
135
            }
 
136
          mpn_sub_n (p + n3, x, y, n2);
 
137
        }
 
138
      p[n] = w;
 
139
 
 
140
      n1 = n + 1;
 
141
      if (n2 < KARATSUBA_MUL_THRESHOLD)
 
142
        {
 
143
          if (n3 < KARATSUBA_MUL_THRESHOLD)
 
144
            {
 
145
              mpn_mul_basecase (ws, p, n3, p + n3, n3);
 
146
              mpn_mul_basecase (p, a, n3, b, n3);
 
147
            }
 
148
          else
 
149
            {
 
150
              mpn_kara_mul_n (ws, p, p + n3, n3, ws + n1);
 
151
              mpn_kara_mul_n (p, a, b, n3, ws + n1);
 
152
            }
 
153
          mpn_mul_basecase (p + n1, a + n3, n2, b + n3, n2);
 
154
        }
 
155
      else
 
156
        {
 
157
          mpn_kara_mul_n (ws, p, p + n3, n3, ws + n1);
 
158
          mpn_kara_mul_n (p, a, b, n3, ws + n1);
 
159
          mpn_kara_mul_n (p + n1, a + n3, b + n3, n2, ws + n1);
 
160
        }
 
161
 
 
162
      if (sign)
 
163
        mpn_add_n (ws, p, ws, n1);
 
164
      else
 
165
        mpn_sub_n (ws, p, ws, n1);
 
166
 
 
167
      nm1 = n - 1;
 
168
      if (mpn_add_n (ws, p + n1, ws, nm1))
 
169
        {
 
170
          mp_limb_t x = ws[nm1] + 1;
 
171
          ws[nm1] = x;
 
172
          if (x == 0)
 
173
            ++ws[n];
 
174
        }
 
175
      if (mpn_add_n (p + n3, p + n3, ws, n1))
 
176
        {
 
177
          mp_limb_t x;
 
178
          i = n1 + n3;
 
179
          do {
 
180
            x = p[i] + 1;
 
181
            p[i] = x;
 
182
            ++i;
 
183
          } while (x == 0);
 
184
        }
 
185
    }
 
186
  else
 
187
    {
 
188
      /* Even length. */
 
189
      i = n2;
 
190
      do
 
191
        {
 
192
          --i;
 
193
          w0 = a[i];
 
194
          w1 = a[n2+i];
 
195
        }
 
196
      while (w0 == w1 && i != 0);
 
197
      sign = 0;
 
198
      if (w0 < w1)
 
199
        {
 
200
          x = a + n2;
 
201
          y = a;
 
202
          sign = ~0;
 
203
        }
 
204
      else
 
205
        {
 
206
          x = a;
 
207
          y = a + n2;
 
208
        }
 
209
      mpn_sub_n (p, x, y, n2);
 
210
 
 
211
      i = n2;
 
212
      do 
 
213
        {
 
214
          --i;
 
215
          w0 = b[i];
 
216
          w1 = b[n2+i];
 
217
        }
 
218
      while (w0 == w1 && i != 0);
 
219
      if (w0 < w1)
 
220
        {
 
221
          x = b + n2;
 
222
          y = b;
 
223
          sign = ~sign;
 
224
        }
 
225
      else
 
226
        {
 
227
          x = b;
 
228
          y = b + n2;
 
229
        }
 
230
      mpn_sub_n (p + n2, x, y, n2);
 
231
 
 
232
      /* Pointwise products. */
 
233
      if (n2 < KARATSUBA_MUL_THRESHOLD)
 
234
        {
 
235
          mpn_mul_basecase (ws, p, n2, p + n2, n2);
 
236
          mpn_mul_basecase (p, a, n2, b, n2);
 
237
          mpn_mul_basecase (p + n, a + n2, n2, b + n2, n2);
 
238
        }
 
239
      else
 
240
        {
 
241
          mpn_kara_mul_n (ws, p, p + n2, n2, ws + n);
 
242
          mpn_kara_mul_n (p, a, b, n2, ws + n);
 
243
          mpn_kara_mul_n (p + n, a + n2, b + n2, n2, ws + n);
 
244
        }
 
245
 
 
246
      /* Interpolate. */
 
247
      if (sign)
 
248
        w = mpn_add_n (ws, p, ws, n);
 
249
      else
 
250
        w = -mpn_sub_n (ws, p, ws, n);
 
251
      w += mpn_add_n (ws, p + n, ws, n);
 
252
      w += mpn_add_n (p + n2, p + n2, ws, n);
 
253
      MPN_INCR_U (p + n2 + n, 2*n - (n2 + n), w);
 
254
    }
 
255
}
 
256
 
 
257
void
 
258
mpn_kara_sqr_n (mp_ptr p, mp_srcptr a, mp_size_t n, mp_ptr ws)
 
259
{
 
260
  mp_limb_t w, w0, w1;
 
261
  mp_size_t n2;
 
262
  mp_srcptr x, y;
 
263
  mp_size_t i;
 
264
 
 
265
  n2 = n >> 1;
 
266
  ASSERT (n2 > 0);
 
267
 
 
268
  if ((n & 1) != 0)
 
269
    {
 
270
      /* Odd length. */
 
271
      mp_size_t n1, n3, nm1;
 
272
 
 
273
      n3 = n - n2;
 
274
 
 
275
      w = a[n2];
 
276
      if (w != 0)
 
277
        w -= mpn_sub_n (p, a, a + n3, n2);
 
278
      else
 
279
        {
 
280
          i = n2;
 
281
          do
 
282
            {
 
283
              --i;
 
284
              w0 = a[i];
 
285
              w1 = a[n3+i];
 
286
            }
 
287
          while (w0 == w1 && i != 0);
 
288
          if (w0 < w1)
 
289
            {
 
290
              x = a + n3;
 
291
              y = a;
 
292
            }
 
293
          else
 
294
            {
 
295
              x = a;
 
296
              y = a + n3;
 
297
            }
 
298
          mpn_sub_n (p, x, y, n2);
 
299
        }
 
300
      p[n2] = w;
 
301
 
 
302
      n1 = n + 1;
 
303
 
 
304
      /* n2 is always either n3 or n3-1 so maybe the two sets of tests here
 
305
         could be combined.  But that's not important, since the tests will
 
306
         take a miniscule amount of time compared to the function calls.  */
 
307
      if (BELOW_THRESHOLD (n3, BASECASE_SQR_THRESHOLD))
 
308
        {
 
309
          mpn_mul_basecase (ws, p, n3, p, n3);
 
310
          mpn_mul_basecase (p,  a, n3, a, n3);
 
311
        }
 
312
      else if (BELOW_THRESHOLD (n3, KARATSUBA_SQR_THRESHOLD))
 
313
        {
 
314
          mpn_sqr_basecase (ws, p, n3);
 
315
          mpn_sqr_basecase (p,  a, n3);
 
316
        }
 
317
      else
 
318
        {
 
319
          mpn_kara_sqr_n   (ws, p, n3, ws + n1);         /* (x-y)^2 */
 
320
          mpn_kara_sqr_n   (p,  a, n3, ws + n1);         /* x^2     */
 
321
        }
 
322
      if (BELOW_THRESHOLD (n2, BASECASE_SQR_THRESHOLD))
 
323
        mpn_mul_basecase (p + n1, a + n3, n2, a + n3, n2);
 
324
      else if (BELOW_THRESHOLD (n2, KARATSUBA_SQR_THRESHOLD))
 
325
        mpn_sqr_basecase (p + n1, a + n3, n2);
 
326
      else
 
327
        mpn_kara_sqr_n   (p + n1, a + n3, n2, ws + n1);  /* y^2     */
 
328
 
 
329
 
 
330
      /* Since x^2+y^2-(x-y)^2 = 2xy >= 0 there's no need to track the
 
331
         borrow from mpn_sub_n.  If it occurs then it'll be cancelled by a
 
332
         carry from ws[n].  Further, since 2xy fits in n1 limbs there won't
 
333
         be any carry out of ws[n] other than cancelling that borrow. */
 
334
 
 
335
      mpn_sub_n (ws, p, ws, n1);             /* x^2-(x-y)^2 */
 
336
 
 
337
      nm1 = n - 1;
 
338
      if (mpn_add_n (ws, p + n1, ws, nm1))   /* x^2+y^2-(x-y)^2 = 2xy */
 
339
        {
 
340
          mp_limb_t x = ws[nm1] + 1;
 
341
          ws[nm1] = x;
 
342
          if (x == 0)
 
343
            ++ws[n];
 
344
        }
 
345
      if (mpn_add_n (p + n3, p + n3, ws, n1))
 
346
        {
 
347
          mp_limb_t x;
 
348
          i = n1 + n3;
 
349
          do {
 
350
            x = p[i] + 1;
 
351
            p[i] = x;
 
352
            ++i;
 
353
          } while (x == 0);
 
354
        }
 
355
    }
 
356
  else
 
357
    {
 
358
      /* Even length. */
 
359
      i = n2;
 
360
      do 
 
361
        {
 
362
          --i;
 
363
          w0 = a[i];
 
364
          w1 = a[n2+i];
 
365
        }
 
366
      while (w0 == w1 && i != 0);
 
367
      if (w0 < w1)
 
368
        {
 
369
          x = a + n2;
 
370
          y = a;
 
371
        }
 
372
      else
 
373
        {
 
374
          x = a;
 
375
          y = a + n2;
 
376
        }
 
377
      mpn_sub_n (p, x, y, n2);
 
378
 
 
379
      /* Pointwise products. */
 
380
      if (BELOW_THRESHOLD (n2, BASECASE_SQR_THRESHOLD))
 
381
        {
 
382
          mpn_mul_basecase (ws,    p,      n2, p,      n2);
 
383
          mpn_mul_basecase (p,     a,      n2, a,      n2);
 
384
          mpn_mul_basecase (p + n, a + n2, n2, a + n2, n2);
 
385
        }
 
386
      else if (BELOW_THRESHOLD (n2, KARATSUBA_SQR_THRESHOLD))
 
387
        {
 
388
          mpn_sqr_basecase (ws,    p,      n2);
 
389
          mpn_sqr_basecase (p,     a,      n2);
 
390
          mpn_sqr_basecase (p + n, a + n2, n2);
 
391
        }
 
392
      else
 
393
        {
 
394
          mpn_kara_sqr_n (ws,    p,      n2, ws + n);
 
395
          mpn_kara_sqr_n (p,     a,      n2, ws + n);
 
396
          mpn_kara_sqr_n (p + n, a + n2, n2, ws + n);
 
397
        }
 
398
 
 
399
      /* Interpolate. */
 
400
      w = -mpn_sub_n (ws, p, ws, n);
 
401
      w += mpn_add_n (ws, p + n, ws, n);
 
402
      w += mpn_add_n (p + n2, p + n2, ws, n);
 
403
      MPN_INCR_U (p + n2 + n, 2*n - (n2 + n), w);
 
404
    }
 
405
}
 
406
 
 
407
/*-- add2Times -------------------------------------------------------------*/
 
408
 
 
409
/* z[] = x[] + 2 * y[]
 
410
   Note that z and x might point to the same vectors.
 
411
   FIXME: gcc won't inline this because it uses alloca. */
 
412
#if USE_MORE_MPN
 
413
static inline mp_limb_t
 
414
add2Times (mp_ptr z, mp_srcptr x, mp_srcptr y, mp_size_t n)
 
415
{
 
416
  mp_ptr t;
 
417
  mp_limb_t c;
 
418
  TMP_DECL (marker);
 
419
  TMP_MARK (marker);
 
420
  t = (mp_ptr) TMP_ALLOC (n * BYTES_PER_MP_LIMB);
 
421
  c = mpn_lshift (t, y, n, 1);
 
422
  c += mpn_add_n (z, x, t, n);
 
423
  TMP_FREE (marker);
 
424
  return c;
 
425
}
 
426
#else
 
427
 
 
428
static mp_limb_t
 
429
add2Times (mp_ptr z, mp_srcptr x, mp_srcptr y, mp_size_t n)
 
430
{
 
431
  mp_limb_t c, v, w;
 
432
 
 
433
  ASSERT (n > 0);
 
434
  v = *x; w = *y;
 
435
  c = w >> (BITS_PER_MP_LIMB - 1);
 
436
  w <<= 1;
 
437
  v += w;
 
438
  c += v < w;
 
439
  *z = v;
 
440
  ++x; ++y; ++z;
 
441
  while (--n)
 
442
    {
 
443
      v = *x;
 
444
      w = *y;
 
445
      v += c;
 
446
      c = v < c;
 
447
      c += w >> (BITS_PER_MP_LIMB - 1);
 
448
      w <<= 1;
 
449
      v += w;
 
450
      c += v < w;
 
451
      *z = v;
 
452
      ++x; ++y; ++z;
 
453
    }
 
454
 
 
455
  return c;
 
456
}
 
457
#endif
 
458
 
 
459
/*-- evaluate3 -------------------------------------------------------------*/
 
460
 
 
461
/* Evaluates:
 
462
 *   ph := 4*A+2*B+C
 
463
 *   p1 := A+B+C
 
464
 *   p2 := A+2*B+4*C
 
465
 * where:
 
466
 *   ph[], p1[], p2[], A[] and B[] all have length len,
 
467
 *   C[] has length len2 with len-len2 = 0, 1 or 2.
 
468
 * Returns top words (overflow) at pth, pt1 and pt2 respectively.
 
469
 */
 
470
#if USE_MORE_MPN
 
471
static void
 
472
evaluate3 (mp_ptr ph, mp_ptr p1, mp_ptr p2, mp_ptr pth, mp_ptr pt1, mp_ptr pt2,
 
473
           mp_srcptr A, mp_srcptr B, mp_srcptr C, mp_size_t len,mp_size_t len2)
 
474
{
 
475
  mp_limb_t c, d, e;
 
476
  
 
477
  ASSERT (len - len2 <= 2);
 
478
 
 
479
  e = mpn_lshift (p1, B, len, 1);
 
480
 
 
481
  c = mpn_lshift (ph, A, len, 2);
 
482
  c += e + mpn_add_n (ph, ph, p1, len);
 
483
  d = mpn_add_n (ph, ph, C, len2);
 
484
  if (len2 == len) c += d; else c += mpn_add_1 (ph + len2, ph + len2, len-len2, d);
 
485
  ASSERT (c < 7);
 
486
  *pth = c;
 
487
 
 
488
  c = mpn_lshift (p2, C, len2, 2);
 
489
#if 1
 
490
  if (len2 != len) { p2[len-1] = 0; p2[len2] = c; c = 0; }
 
491
  c += e + mpn_add_n (p2, p2, p1, len);
 
492
#else
 
493
  d = mpn_add_n (p2, p2, p1, len2);
 
494
  c += d;
 
495
  if (len2 != len) c = mpn_add_1 (p2+len2, p1+len2, len-len2, c);
 
496
  c += e;
 
497
#endif
 
498
  c += mpn_add_n (p2, p2, A, len);
 
499
  ASSERT (c < 7);
 
500
  *pt2 = c;
 
501
 
 
502
  c = mpn_add_n (p1, A, B, len);
 
503
  d = mpn_add_n (p1, p1, C, len2);
 
504
  if (len2 == len) c += d;
 
505
  else c += mpn_add_1 (p1+len2, p1+len2, len-len2, d);
 
506
  ASSERT (c < 3);
 
507
  *pt1 = c;
 
508
 
 
509
}
 
510
 
 
511
#else
 
512
 
 
513
static void
 
514
evaluate3 (mp_ptr ph, mp_ptr p1, mp_ptr p2, mp_ptr pth, mp_ptr pt1, mp_ptr pt2,
 
515
           mp_srcptr A, mp_srcptr B, mp_srcptr C, mp_size_t l, mp_size_t ls)
 
516
{
 
517
  mp_limb_t a,b,c, i, t, th,t1,t2, vh,v1,v2;
 
518
 
 
519
  ASSERT (l - ls <= 2);
 
520
 
 
521
  th = t1 = t2 = 0;
 
522
  for (i = 0; i < l; ++i)
 
523
    {
 
524
      a = *A;
 
525
      b = *B;
 
526
      c = i < ls ? *C : 0;
 
527
 
 
528
      /* TO DO: choose one of the following alternatives. */
 
529
#if 0
 
530
      t = a << 2;
 
531
      vh = th + t;
 
532
      th = vh < t;
 
533
      th += a >> (BITS_PER_MP_LIMB - 2);
 
534
      t = b << 1;
 
535
      vh += t;
 
536
      th += vh < t;
 
537
      th += b >> (BITS_PER_MP_LIMB - 1);
 
538
      vh += c;
 
539
      th += vh < c;
 
540
#else
 
541
      vh = th + c;
 
542
      th = vh < c;
 
543
      t = b << 1;
 
544
      vh += t;
 
545
      th += vh < t;
 
546
      th += b >> (BITS_PER_MP_LIMB - 1);
 
547
      t = a << 2;
 
548
      vh += t;
 
549
      th += vh < t;
 
550
      th += a >> (BITS_PER_MP_LIMB - 2);
 
551
#endif
 
552
 
 
553
      v1 = t1 + a;
 
554
      t1 = v1 < a;
 
555
      v1 += b;
 
556
      t1 += v1 < b;
 
557
      v1 += c;
 
558
      t1 += v1 < c;
 
559
 
 
560
      v2 = t2 + a;
 
561
      t2 = v2 < a;
 
562
      t = b << 1;
 
563
      v2 += t;
 
564
      t2 += v2 < t;
 
565
      t2 += b >> (BITS_PER_MP_LIMB - 1);
 
566
      t = c << 2;
 
567
      v2 += t;
 
568
      t2 += v2 < t;
 
569
      t2 += c >> (BITS_PER_MP_LIMB - 2);
 
570
 
 
571
      *ph = vh;
 
572
      *p1 = v1;
 
573
      *p2 = v2;
 
574
 
 
575
      ++A; ++B; ++C;
 
576
      ++ph; ++p1; ++p2;
 
577
    }
 
578
 
 
579
  ASSERT (th < 7);
 
580
  ASSERT (t1 < 3);
 
581
  ASSERT (t2 < 7);
 
582
 
 
583
  *pth = th;
 
584
  *pt1 = t1;
 
585
  *pt2 = t2;
 
586
}
 
587
#endif
 
588
 
 
589
 
 
590
/*-- interpolate3 ----------------------------------------------------------*/
 
591
 
 
592
/* Interpolates B, C, D (in-place) from:
 
593
 *   16*A+8*B+4*C+2*D+E
 
594
 *   A+B+C+D+E
 
595
 *   A+2*B+4*C+8*D+16*E
 
596
 * where:
 
597
 *   A[], B[], C[] and D[] all have length l,
 
598
 *   E[] has length ls with l-ls = 0, 2 or 4.
 
599
 *
 
600
 * Reads top words (from earlier overflow) from ptb, ptc and ptd,
 
601
 * and returns new top words there.
 
602
 */
 
603
 
 
604
#if USE_MORE_MPN
 
605
static void
 
606
interpolate3 (mp_srcptr A, mp_ptr B, mp_ptr C, mp_ptr D, mp_srcptr E,
 
607
              mp_ptr ptb, mp_ptr ptc, mp_ptr ptd, mp_size_t len,mp_size_t len2)
 
608
{
 
609
  mp_ptr ws;
 
610
  mp_limb_t t, tb,tc,td;
 
611
  TMP_DECL (marker);
 
612
  TMP_MARK (marker);
 
613
 
 
614
  ASSERT (len - len2 == 0 || len - len2 == 2 || len - len2 == 4);
 
615
 
 
616
  /* Let x1, x2, x3 be the values to interpolate.  We have:
 
617
   *         b = 16*a + 8*x1 + 4*x2 + 2*x3 +    e
 
618
   *         c =    a +   x1 +   x2 +   x3 +    e
 
619
   *         d =    a + 2*x1 + 4*x2 + 8*x3 + 16*e
 
620
   */
 
621
 
 
622
  ws = (mp_ptr) TMP_ALLOC (len * BYTES_PER_MP_LIMB);
 
623
 
 
624
  tb = *ptb; tc = *ptc; td = *ptd;
 
625
 
 
626
 
 
627
  /* b := b - 16*a -    e
 
628
   * c := c -    a -    e
 
629
   * d := d -    a - 16*e
 
630
   */
 
631
 
 
632
  t = mpn_lshift (ws, A, len, 4);
 
633
  tb -= t + mpn_sub_n (B, B, ws, len);
 
634
  t = mpn_sub_n (B, B, E, len2);
 
635
  if (len2 == len) tb -= t;
 
636
  else tb -= mpn_sub_1 (B+len2, B+len2, len-len2, t);
 
637
 
 
638
  tc -= mpn_sub_n (C, C, A, len);
 
639
  t = mpn_sub_n (C, C, E, len2);
 
640
  if (len2 == len) tc -= t;
 
641
  else tc -= mpn_sub_1 (C+len2, C+len2, len-len2, t);
 
642
 
 
643
  t = mpn_lshift (ws, E, len2, 4);
 
644
  t += mpn_add_n (ws, ws, A, len2);
 
645
#if 1
 
646
  if (len2 != len) t = mpn_add_1 (ws+len2, A+len2, len-len2, t);
 
647
  td -= t + mpn_sub_n (D, D, ws, len);
 
648
#else
 
649
  t += mpn_sub_n (D, D, ws, len2);
 
650
  if (len2 != len) {
 
651
    t = mpn_sub_1 (D+len2, D+len2, len-len2, t);
 
652
    t += mpn_sub_n (D+len2, D+len2, A+len2, len-len2);
 
653
  } /* end if/else */
 
654
  td -= t;
 
655
#endif
 
656
 
 
657
 
 
658
  /* b, d := b + d, b - d */
 
659
 
 
660
#ifdef HAVE_MPN_ADD_SUB_N
 
661
  /* #error TO DO ... */
 
662
#else
 
663
  t = tb + td + mpn_add_n (ws, B, D, len);  
 
664
  td = tb - td - mpn_sub_n (D, B, D, len);
 
665
  tb = t;
 
666
  MPN_COPY (B, ws, len);
 
667
#endif
 
668
  
 
669
  /* b := b-8*c */
 
670
  t = 8 * tc + mpn_lshift (ws, C, len, 3);
 
671
  tb -= t + mpn_sub_n (B, B, ws, len);
 
672
 
 
673
  /* c := 2*c - b */
 
674
  tc = 2 * tc + mpn_lshift (C, C, len, 1);
 
675
  tc -= tb + mpn_sub_n (C, C, B, len);
 
676
 
 
677
  /* d := d/3 */
 
678
  td = (td - mpn_divexact_by3 (D, D, len)) * MODLIMB_INVERSE_3;
 
679
 
 
680
  /* b, d := b + d, b - d */
 
681
#ifdef HAVE_MPN_ADD_SUB_N
 
682
  /* #error TO DO ... */
 
683
#else
 
684
  t = tb + td + mpn_add_n (ws, B, D, len);  
 
685
  td = tb - td - mpn_sub_n (D, B, D, len);
 
686
  tb = t;
 
687
  MPN_COPY (B, ws, len);
 
688
#endif
 
689
 
 
690
      /* Now:
 
691
       *         b = 4*x1
 
692
       *         c = 2*x2
 
693
       *         d = 4*x3
 
694
       */
 
695
 
 
696
  ASSERT(!(*B & 3));
 
697
  mpn_rshift (B, B, len, 2);
 
698
  B[len-1] |= tb<<(BITS_PER_MP_LIMB-2);
 
699
  ASSERT((mp_limb_signed_t)tb >= 0);
 
700
  tb >>= 2;
 
701
 
 
702
  ASSERT(!(*C & 1));
 
703
  mpn_rshift (C, C, len, 1);
 
704
  C[len-1] |= tc<<(BITS_PER_MP_LIMB-1);
 
705
  ASSERT((mp_limb_signed_t)tc >= 0);
 
706
  tc >>= 1;
 
707
 
 
708
  ASSERT(!(*D & 3));
 
709
  mpn_rshift (D, D, len, 2);
 
710
  D[len-1] |= td<<(BITS_PER_MP_LIMB-2);
 
711
  ASSERT((mp_limb_signed_t)td >= 0);
 
712
  td >>= 2;
 
713
 
 
714
#if WANT_ASSERT
 
715
  ASSERT (tb < 2);
 
716
  if (len == len2)
 
717
    {
 
718
      ASSERT (tc < 3);
 
719
      ASSERT (td < 2);
 
720
    }
 
721
  else
 
722
    {
 
723
      ASSERT (tc < 2);
 
724
      ASSERT (!td);
 
725
    }
 
726
#endif
 
727
 
 
728
  *ptb = tb;
 
729
  *ptc = tc;
 
730
  *ptd = td;
 
731
 
 
732
  TMP_FREE (marker);
 
733
}
 
734
 
 
735
#else
 
736
 
 
737
static void
 
738
interpolate3 (mp_srcptr A, mp_ptr B, mp_ptr C, mp_ptr D, mp_srcptr E,
 
739
              mp_ptr ptb, mp_ptr ptc, mp_ptr ptd, mp_size_t l, mp_size_t ls)
 
740
{
 
741
  mp_limb_t a,b,c,d,e,t, i, sb,sc,sd, ob,oc,od;
 
742
  const mp_limb_t maskOffHalf = (~(mp_limb_t) 0) << (BITS_PER_MP_LIMB >> 1);
 
743
 
 
744
#if WANT_ASSERT
 
745
  t = l - ls;
 
746
  ASSERT (t == 0 || t == 2 || t == 4);
 
747
#endif
 
748
 
 
749
  sb = sc = sd = 0;
 
750
  for (i = 0; i < l; ++i)
 
751
    {
 
752
      mp_limb_t tb, tc, td, tt;
 
753
 
 
754
      a = *A;
 
755
      b = *B;
 
756
      c = *C;
 
757
      d = *D;
 
758
      e = i < ls ? *E : 0;
 
759
 
 
760
      /* Let x1, x2, x3 be the values to interpolate.  We have:
 
761
       *         b = 16*a + 8*x1 + 4*x2 + 2*x3 +    e
 
762
       *         c =    a +   x1 +   x2 +   x3 +    e
 
763
       *         d =    a + 2*x1 + 4*x2 + 8*x3 + 16*e
 
764
       */
 
765
 
 
766
      /* b := b - 16*a -    e
 
767
       * c := c -    a -    e
 
768
       * d := d -    a - 16*e
 
769
       */
 
770
      t = a << 4;
 
771
      tb = -(a >> (BITS_PER_MP_LIMB - 4)) - (b < t);
 
772
      b -= t;
 
773
      tb -= b < e;
 
774
      b -= e;
 
775
      tc = -(c < a);
 
776
      c -= a;
 
777
      tc -= c < e;
 
778
      c -= e;
 
779
      td = -(d < a);
 
780
      d -= a;
 
781
      t = e << 4;
 
782
      td = td - (e >> (BITS_PER_MP_LIMB - 4)) - (d < t);
 
783
      d -= t;
 
784
 
 
785
      /* b, d := b + d, b - d */
 
786
      t = b + d;
 
787
      tt = tb + td + (t < b);
 
788
      td = tb - td - (b < d);
 
789
      d = b - d;
 
790
      b = t;
 
791
      tb = tt;
 
792
 
 
793
      /* b := b-8*c */
 
794
      t = c << 3;
 
795
      tb = tb - (tc << 3) - (c >> (BITS_PER_MP_LIMB - 3)) - (b < t);
 
796
      b -= t;
 
797
 
 
798
      /* c := 2*c - b */
 
799
      t = c << 1;
 
800
      tc = (tc << 1) + (c >> (BITS_PER_MP_LIMB - 1)) - tb - (t < b);
 
801
      c = t - b;
 
802
 
 
803
      /* d := d/3 */
 
804
      d *= MODLIMB_INVERSE_3;
 
805
      td = td - (d >> (BITS_PER_MP_LIMB - 1)) - (d*3 < d);
 
806
      td *= MODLIMB_INVERSE_3;
 
807
 
 
808
      /* b, d := b + d, b - d */
 
809
      t = b + d;
 
810
      tt = tb + td + (t < b);
 
811
      td = tb - td - (b < d);
 
812
      d = b - d;
 
813
      b = t;
 
814
      tb = tt;
 
815
 
 
816
      /* Now:
 
817
       *         b = 4*x1
 
818
       *         c = 2*x2
 
819
       *         d = 4*x3
 
820
       */
 
821
 
 
822
      /* sb has period 2. */
 
823
      b += sb;
 
824
      tb += b < sb;
 
825
      sb &= maskOffHalf;
 
826
      sb |= sb >> (BITS_PER_MP_LIMB >> 1);
 
827
      sb += tb;
 
828
 
 
829
      /* sc has period 1. */
 
830
      c += sc;
 
831
      tc += c < sc;
 
832
      /* TO DO: choose one of the following alternatives. */
 
833
#if 1
 
834
      sc = (mp_limb_signed_t) sc >> (BITS_PER_MP_LIMB - 1);
 
835
      sc += tc;
 
836
#else
 
837
      sc = tc - ((mp_limb_signed_t) sc < 0L);
 
838
#endif
 
839
 
 
840
      /* sd has period 2. */
 
841
      d += sd;
 
842
      td += d < sd;
 
843
      sd &= maskOffHalf;
 
844
      sd |= sd >> (BITS_PER_MP_LIMB >> 1);
 
845
      sd += td;
 
846
 
 
847
      if (i != 0)
 
848
        {
 
849
          B[-1] = ob | b << (BITS_PER_MP_LIMB - 2);
 
850
          C[-1] = oc | c << (BITS_PER_MP_LIMB - 1);
 
851
          D[-1] = od | d << (BITS_PER_MP_LIMB - 2);
 
852
        }
 
853
      ob = b >> 2;
 
854
      oc = c >> 1;
 
855
      od = d >> 2;
 
856
 
 
857
      ++A; ++B; ++C; ++D; ++E;
 
858
    }
 
859
 
 
860
  /* Handle top words. */
 
861
  b = *ptb;
 
862
  c = *ptc;
 
863
  d = *ptd;
 
864
 
 
865
  t = b + d;
 
866
  d = b - d;
 
867
  b = t;
 
868
  b -= c << 3;
 
869
  c = (c << 1) - b;
 
870
  d *= MODLIMB_INVERSE_3;
 
871
  t = b + d;
 
872
  d = b - d;
 
873
  b = t;
 
874
 
 
875
  b += sb;
 
876
  c += sc;
 
877
  d += sd;
 
878
 
 
879
  B[-1] = ob | b << (BITS_PER_MP_LIMB - 2);
 
880
  C[-1] = oc | c << (BITS_PER_MP_LIMB - 1);
 
881
  D[-1] = od | d << (BITS_PER_MP_LIMB - 2);
 
882
 
 
883
  b >>= 2;
 
884
  c >>= 1;
 
885
  d >>= 2;
 
886
 
 
887
#if WANT_ASSERT
 
888
  ASSERT (b < 2);
 
889
  if (l == ls)
 
890
    {
 
891
      ASSERT (c < 3);
 
892
      ASSERT (d < 2);
 
893
    }
 
894
  else
 
895
    {
 
896
      ASSERT (c < 2);
 
897
      ASSERT (!d);
 
898
    }
 
899
#endif
 
900
 
 
901
  *ptb = b;
 
902
  *ptc = c;
 
903
  *ptd = d;
 
904
}
 
905
#endif
 
906
 
 
907
 
 
908
/*-- mpn_toom3_mul_n --------------------------------------------------------------*/
 
909
 
 
910
/* Multiplies using 5 mults of one third size and so on recursively.
 
911
 * p[0..2*n-1] := product of a[0..n-1] and b[0..n-1].
 
912
 * No overlap of p[...] with a[...] or b[...].
 
913
 * ws is workspace.
 
914
 */
 
915
 
 
916
/* TO DO: If TOOM3_MUL_THRESHOLD is much bigger than KARATSUBA_MUL_THRESHOLD then the
 
917
 *        recursion in mpn_toom3_mul_n() will always bottom out with mpn_kara_mul_n()
 
918
 *        because the "n < KARATSUBA_MUL_THRESHOLD" test here will always be false.
 
919
 */
 
920
 
 
921
#define TOOM3_MUL_REC(p, a, b, n, ws) \
 
922
  do {                                                          \
 
923
    if (n < KARATSUBA_MUL_THRESHOLD)                            \
 
924
      mpn_mul_basecase (p, a, n, b, n);                         \
 
925
    else if (n < TOOM3_MUL_THRESHOLD)                           \
 
926
      mpn_kara_mul_n (p, a, b, n, ws);                          \
 
927
    else                                                        \
 
928
      mpn_toom3_mul_n (p, a, b, n, ws);                         \
 
929
  } while (0)
 
930
 
 
931
void
 
932
mpn_toom3_mul_n (mp_ptr p, mp_srcptr a, mp_srcptr b, mp_size_t n, mp_ptr ws)
 
933
{
 
934
  mp_limb_t cB,cC,cD, dB,dC,dD, tB,tC,tD;
 
935
  mp_limb_t *A,*B,*C,*D,*E, *W;
 
936
  mp_size_t l,l2,l3,l4,l5,ls;
 
937
 
 
938
  /* Break n words into chunks of size l, l and ls.
 
939
   * n = 3*k   => l = k,   ls = k
 
940
   * n = 3*k+1 => l = k+1, ls = k-1
 
941
   * n = 3*k+2 => l = k+1, ls = k
 
942
   */
 
943
  {
 
944
    mp_limb_t m;
 
945
 
 
946
    /* this is probably unnecessarily strict */
 
947
    ASSERT (n >= TOOM3_MUL_THRESHOLD);
 
948
 
 
949
    l = ls = n / 3;
 
950
    m = n - l * 3;
 
951
    if (m != 0)
 
952
      ++l;
 
953
    if (m == 1)
 
954
      --ls;
 
955
 
 
956
    l2 = l * 2;
 
957
    l3 = l * 3;
 
958
    l4 = l * 4;
 
959
    l5 = l * 5;
 
960
    A = p;
 
961
    B = ws;
 
962
    C = p + l2;
 
963
    D = ws + l2;
 
964
    E = p + l4;
 
965
    W = ws + l4;
 
966
  }
 
967
 
 
968
  ASSERT (l >= 1);
 
969
  ASSERT (ls >= 1);
 
970
 
 
971
  /** First stage: evaluation at points 0, 1/2, 1, 2, oo. **/
 
972
  evaluate3 (A, B, C, &cB, &cC, &cD, a, a + l, a + l2, l, ls);
 
973
  evaluate3 (A + l, B + l, C + l, &dB, &dC, &dD, b, b + l, b + l2, l, ls);
 
974
 
 
975
  /** Second stage: pointwise multiplies. **/
 
976
  TOOM3_MUL_REC(D, C, C + l, l, W);
 
977
  tD = cD*dD;
 
978
  if (cD) tD += mpn_addmul_1 (D + l, C + l, l, cD);
 
979
  if (dD) tD += mpn_addmul_1 (D + l, C, l, dD);
 
980
  ASSERT (tD < 49);
 
981
  TOOM3_MUL_REC(C, B, B + l, l, W);
 
982
  tC = cC*dC;
 
983
  /* TO DO: choose one of the following alternatives. */
 
984
#if 0
 
985
  if (cC) tC += mpn_addmul_1 (C + l, B + l, l, cC);
 
986
  if (dC) tC += mpn_addmul_1 (C + l, B, l, dC);
 
987
#else
 
988
  if (cC)
 
989
    {
 
990
      if (cC == 1) tC += mpn_add_n (C + l, C + l, B + l, l);
 
991
      else tC += add2Times (C + l, C + l, B + l, l);
 
992
    }
 
993
  if (dC)
 
994
    {
 
995
      if (dC == 1) tC += mpn_add_n (C + l, C + l, B, l);
 
996
      else tC += add2Times (C + l, C + l, B, l);
 
997
    }
 
998
#endif
 
999
  ASSERT (tC < 9);
 
1000
  TOOM3_MUL_REC(B, A, A + l, l, W);
 
1001
  tB = cB*dB;
 
1002
  if (cB) tB += mpn_addmul_1 (B + l, A + l, l, cB);
 
1003
  if (dB) tB += mpn_addmul_1 (B + l, A, l, dB);
 
1004
  ASSERT (tB < 49);
 
1005
  TOOM3_MUL_REC(A, a, b, l, W);
 
1006
  TOOM3_MUL_REC(E, a + l2, b + l2, ls, W);
 
1007
 
 
1008
  /** Third stage: interpolation. **/
 
1009
  interpolate3 (A, B, C, D, E, &tB, &tC, &tD, l2, ls << 1);
 
1010
 
 
1011
  /** Final stage: add up the coefficients. **/
 
1012
  tB += mpn_add_n (p + l, p + l, B, l2);
 
1013
  tD += mpn_add_n (p + l3, p + l3, D, l2);
 
1014
  MPN_INCR_U (p + l3, 2*n - l3, tB);
 
1015
  MPN_INCR_U (p + l4, 2*n - l4, tC);
 
1016
  MPN_INCR_U (p + l5, 2*n - l5, tD);
 
1017
}
 
1018
 
 
1019
/*-- mpn_toom3_sqr_n --------------------------------------------------------------*/
 
1020
 
 
1021
/* Like previous function but for squaring */
 
1022
 
 
1023
/* FIXME: If TOOM3_SQR_THRESHOLD is big enough it might never get into the
 
1024
   basecase range.  Try to arrange those conditonals go dead.  */
 
1025
#define TOOM3_SQR_REC(p, a, n, ws)                              \
 
1026
  do {                                                          \
 
1027
    if (BELOW_THRESHOLD (n, BASECASE_SQR_THRESHOLD))            \
 
1028
      mpn_mul_basecase (p, a, n, a, n);                         \
 
1029
    else if (BELOW_THRESHOLD (n, KARATSUBA_SQR_THRESHOLD))      \
 
1030
      mpn_sqr_basecase (p, a, n);                               \
 
1031
    else if (BELOW_THRESHOLD (n, TOOM3_SQR_THRESHOLD))          \
 
1032
      mpn_kara_sqr_n (p, a, n, ws);                             \
 
1033
    else                                                        \
 
1034
      mpn_toom3_sqr_n (p, a, n, ws);                            \
 
1035
  } while (0)
 
1036
 
 
1037
void
 
1038
mpn_toom3_sqr_n (mp_ptr p, mp_srcptr a, mp_size_t n, mp_ptr ws)
 
1039
{
 
1040
  mp_limb_t cB,cC,cD, tB,tC,tD;
 
1041
  mp_limb_t *A,*B,*C,*D,*E, *W;
 
1042
  mp_size_t l,l2,l3,l4,l5,ls;
 
1043
 
 
1044
  /* Break n words into chunks of size l, l and ls.
 
1045
   * n = 3*k   => l = k,   ls = k
 
1046
   * n = 3*k+1 => l = k+1, ls = k-1
 
1047
   * n = 3*k+2 => l = k+1, ls = k
 
1048
   */
 
1049
  {
 
1050
    mp_limb_t m;
 
1051
 
 
1052
    /* this is probably unnecessarily strict */
 
1053
    ASSERT (n >= TOOM3_SQR_THRESHOLD);
 
1054
 
 
1055
    l = ls = n / 3;
 
1056
    m = n - l * 3;
 
1057
    if (m != 0)
 
1058
      ++l;
 
1059
    if (m == 1)
 
1060
      --ls;
 
1061
 
 
1062
    l2 = l * 2;
 
1063
    l3 = l * 3;
 
1064
    l4 = l * 4;
 
1065
    l5 = l * 5;
 
1066
    A = p;
 
1067
    B = ws;
 
1068
    C = p + l2;
 
1069
    D = ws + l2;
 
1070
    E = p + l4;
 
1071
    W = ws + l4;
 
1072
  }
 
1073
 
 
1074
  ASSERT (l >= 1);
 
1075
  ASSERT (ls >= 1);
 
1076
 
 
1077
  /** First stage: evaluation at points 0, 1/2, 1, 2, oo. **/
 
1078
  evaluate3 (A, B, C, &cB, &cC, &cD, a, a + l, a + l2, l, ls);
 
1079
 
 
1080
  /** Second stage: pointwise multiplies. **/
 
1081
  TOOM3_SQR_REC(D, C, l, W);
 
1082
  tD = cD*cD;
 
1083
  if (cD) tD += mpn_addmul_1 (D + l, C, l, 2*cD);
 
1084
  ASSERT (tD < 49);
 
1085
  TOOM3_SQR_REC(C, B, l, W);
 
1086
  tC = cC*cC;
 
1087
  /* TO DO: choose one of the following alternatives. */
 
1088
#if 0
 
1089
  if (cC) tC += mpn_addmul_1 (C + l, B, l, 2*cC);
 
1090
#else
 
1091
  if (cC >= 1)
 
1092
    {
 
1093
      tC += add2Times (C + l, C + l, B, l);
 
1094
      if (cC == 2)
 
1095
        tC += add2Times (C + l, C + l, B, l);
 
1096
    }
 
1097
#endif
 
1098
  ASSERT (tC < 9);
 
1099
  TOOM3_SQR_REC(B, A, l, W);
 
1100
  tB = cB*cB;
 
1101
  if (cB) tB += mpn_addmul_1 (B + l, A, l, 2*cB);
 
1102
  ASSERT (tB < 49);
 
1103
  TOOM3_SQR_REC(A, a, l, W);
 
1104
  TOOM3_SQR_REC(E, a + l2, ls, W);
 
1105
 
 
1106
  /** Third stage: interpolation. **/
 
1107
  interpolate3 (A, B, C, D, E, &tB, &tC, &tD, l2, ls << 1);
 
1108
 
 
1109
  /** Final stage: add up the coefficients. **/
 
1110
  tB += mpn_add_n (p + l, p + l, B, l2);
 
1111
  tD += mpn_add_n (p + l3, p + l3, D, l2);
 
1112
  MPN_INCR_U (p + l3, 2*n - l3, tB);
 
1113
  MPN_INCR_U (p + l4, 2*n - l4, tC);
 
1114
  MPN_INCR_U (p + l5, 2*n - l5, tD);
 
1115
}
 
1116
 
 
1117
void
 
1118
mpn_mul_n (mp_ptr p, mp_srcptr a, mp_srcptr b, mp_size_t n)
 
1119
{
 
1120
  ASSERT (n >= 1);
 
1121
  ASSERT (! MPN_OVERLAP_P (p, 2*n, a, n));
 
1122
  ASSERT (! MPN_OVERLAP_P (p, 2*n, b, n));
 
1123
 
 
1124
  if (n < KARATSUBA_MUL_THRESHOLD)
 
1125
    mpn_mul_basecase (p, a, n, b, n);
 
1126
  else if (n < TOOM3_MUL_THRESHOLD)
 
1127
    {
 
1128
      /* Allocate workspace of fixed size on stack: fast! */
 
1129
#if TUNE_PROGRAM_BUILD
 
1130
      mp_limb_t ws[MPN_KARA_MUL_N_TSIZE (TOOM3_MUL_THRESHOLD_LIMIT-1)];
 
1131
#else
 
1132
      mp_limb_t ws[MPN_KARA_MUL_N_TSIZE (TOOM3_MUL_THRESHOLD-1)];
 
1133
#endif
 
1134
      mpn_kara_mul_n (p, a, b, n, ws);
 
1135
    }
 
1136
#if WANT_FFT || TUNE_PROGRAM_BUILD
 
1137
  else if (n < FFT_MUL_THRESHOLD)
 
1138
#else
 
1139
  else
 
1140
#endif
 
1141
    {
 
1142
      /* Use workspace of unknown size in heap, as stack space may
 
1143
       * be limited.  Since n is at least TOOM3_MUL_THRESHOLD, the
 
1144
       * multiplication will take much longer than malloc()/free().  */
 
1145
      mp_limb_t wsLen, *ws;
 
1146
      wsLen = MPN_TOOM3_MUL_N_TSIZE (n);
 
1147
#ifdef BAD_ALLOCA
 
1148
      ws = __GMP_ALLOCATE_FUNC_LIMBS ((size_t) wsLen);
 
1149
#else
 
1150
      ws = TMP_ALLOC ((size_t) wsLen * sizeof(mp_limb_t));
 
1151
#endif
 
1152
      mpn_toom3_mul_n (p, a, b, n, ws);
 
1153
#ifdef BAD_ALLOCA
 
1154
      __GMP_FREE_FUNC_LIMBS (ws, (size_t) wsLen);
 
1155
#endif
 
1156
    }
 
1157
#if WANT_FFT || TUNE_PROGRAM_BUILD
 
1158
  else
 
1159
    {
 
1160
      mpn_mul_fft_full (p, a, n, b, n);      
 
1161
    }
 
1162
#endif
 
1163
}