Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpn/mul.c
Line
Count
Source (jump to first uncovered line)
1
/* mpn_mul -- Multiply two natural numbers.
2
3
   Contributed to the GNU project by Torbjorn Granlund.
4
5
Copyright 1991, 1993, 1994, 1996, 1997, 1999-2003, 2005-2007, 2009, 2010, 2012,
6
2014, 2019 Free Software Foundation, Inc.
7
8
This file is part of the GNU MP Library.
9
10
The GNU MP Library is free software; you can redistribute it and/or modify
11
it under the terms of either:
12
13
  * the GNU Lesser General Public License as published by the Free
14
    Software Foundation; either version 3 of the License, or (at your
15
    option) any later version.
16
17
or
18
19
  * the GNU General Public License as published by the Free Software
20
    Foundation; either version 2 of the License, or (at your option) any
21
    later version.
22
23
or both in parallel, as here.
24
25
The GNU MP Library is distributed in the hope that it will be useful, but
26
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
27
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
28
for more details.
29
30
You should have received copies of the GNU General Public License and the
31
GNU Lesser General Public License along with the GNU MP Library.  If not,
32
see https://www.gnu.org/licenses/.  */
33
34
#include "gmp-impl.h"
35
36
37
#ifndef MUL_BASECASE_MAX_UN
38
46.1k
#define MUL_BASECASE_MAX_UN 500
39
#endif
40
41
/* Areas where the different toom algorithms can be called (extracted
42
   from the t-toom*.c files, and ignoring small constant offsets):
43
44
   1/6  1/5 1/4 4/13 1/3 3/8 2/5 5/11 1/2 3/5 2/3 3/4 4/5   1 vn/un
45
                                        4/7              6/7
46
               6/11
47
                                       |--------------------| toom22 (small)
48
                                                           || toom22 (large)
49
                                                       |xxxx| toom22 called
50
                      |-------------------------------------| toom32
51
                                         |xxxxxxxxxxxxxxxx| | toom32 called
52
                                               |------------| toom33
53
                                                          |x| toom33 called
54
             |---------------------------------|            | toom42
55
                |xxxxxxxxxxxxxxxxxxxxxxxx|            | toom42 called
56
                                       |--------------------| toom43
57
                                               |xxxxxxxxxx|   toom43 called
58
         |-----------------------------|                      toom52 (unused)
59
                                                   |--------| toom44
60
               |xxxxxxxx| toom44 called
61
                              |--------------------|        | toom53
62
                                        |xxxxxx|              toom53 called
63
    |-------------------------|                               toom62 (unused)
64
                                           |----------------| toom54 (unused)
65
                      |--------------------|                  toom63
66
                        |xxxxxxxxx|                   | toom63 called
67
                          |---------------------------------| toom6h
68
               |xxxxxxxx| toom6h called
69
                                  |-------------------------| toom8h (32 bit)
70
                 |------------------------------------------| toom8h (64 bit)
71
               |xxxxxxxx| toom8h called
72
*/
73
74
#define TOOM33_OK(an,bn) (6 + 2 * an < 3 * bn)
75
401
#define TOOM44_OK(an,bn) (12 + 3 * an < 4 * bn)
76
77
/* Multiply the natural numbers u (pointed to by UP, with UN limbs) and v
78
   (pointed to by VP, with VN limbs), and store the result at PRODP.  The
79
   result is UN + VN limbs.  Return the most significant limb of the result.
80
81
   NOTE: The space pointed to by PRODP is overwritten before finished with U
82
   and V, so overlap is an error.
83
84
   Argument constraints:
85
   1. UN >= VN.
86
   2. PRODP != UP and PRODP != VP, i.e. the destination must be distinct from
87
      the multiplier and the multiplicand.  */
88
89
/*
90
  * The cutoff lines in the toomX2 and toomX3 code are now exactly between the
91
    ideal lines of the surrounding algorithms.  Is that optimal?
92
93
  * The toomX3 code now uses a structure similar to the one of toomX2, except
94
    that it loops longer in the unbalanced case.  The result is that the
95
    remaining area might have un < vn.  Should we fix the toomX2 code in a
96
    similar way?
97
98
  * The toomX3 code is used for the largest non-FFT unbalanced operands.  It
99
    therefore calls mpn_mul recursively for certain cases.
100
101
  * Allocate static temp space using THRESHOLD variables (except for toom44
102
    when !WANT_FFT).  That way, we can typically have no TMP_ALLOC at all.
103
104
  * We sort ToomX2 algorithms together, assuming the toom22, toom32, toom42
105
    have the same vn threshold.  This is not true, we should actually use
106
    mul_basecase for slightly larger operands for toom32 than for toom22, and
107
    even larger for toom42.
108
109
  * That problem is even more prevalent for toomX3.  We therefore use special
110
    THRESHOLD variables there.
111
*/
112
113
mp_limb_t
114
mpn_mul (mp_ptr prodp,
115
   mp_srcptr up, mp_size_t un,
116
   mp_srcptr vp, mp_size_t vn)
117
202k
{
118
202k
  ASSERT (un >= vn);
119
202k
  ASSERT (vn >= 1);
120
202k
  ASSERT (! MPN_OVERLAP_P (prodp, un+vn, up, un));
121
202k
  ASSERT (! MPN_OVERLAP_P (prodp, un+vn, vp, vn));
122
123
202k
  if (BELOW_THRESHOLD (un, MUL_TOOM22_THRESHOLD))
124
134k
    {
125
      /* When un (and thus vn) is below the toom22 range, do mul_basecase.
126
   Test un and not vn here not to thwart the un >> vn code below.
127
   This special case is not necessary, but cuts the overhead for the
128
   smallest operands. */
129
134k
      mpn_mul_basecase (prodp, up, un, vp, vn);
130
134k
    }
131
68.8k
  else if (un == vn)
132
5.47k
    {
133
5.47k
      mpn_mul_n (prodp, up, vp, un);
134
5.47k
    }
135
63.3k
  else if (vn < MUL_TOOM22_THRESHOLD)
136
23.0k
    { /* plain schoolbook multiplication */
137
138
      /* Unless un is very large, or else if have an applicable mpn_mul_N,
139
   perform basecase multiply directly.  */
140
23.0k
      if (un <= MUL_BASECASE_MAX_UN
141
23.0k
#if HAVE_NATIVE_mpn_mul_2
142
23.0k
    || vn <= 2
143
#else
144
    || vn == 1
145
#endif
146
23.0k
    )
147
23.0k
  mpn_mul_basecase (prodp, up, un, vp, vn);
148
0
      else
149
0
  {
150
    /* We have un >> MUL_BASECASE_MAX_UN > vn.  For better memory
151
       locality, split up[] into MUL_BASECASE_MAX_UN pieces and multiply
152
       these pieces with the vp[] operand.  After each such partial
153
       multiplication (but the last) we copy the most significant vn
154
       limbs into a temporary buffer since that part would otherwise be
155
       overwritten by the next multiplication.  After the next
156
       multiplication, we add it back.  This illustrates the situation:
157
158
                                                    -->vn<--
159
                                                      |  |<------- un ------->|
160
                                                         _____________________|
161
                                                        X                    /|
162
                                                      /XX__________________/  |
163
                                    _____________________                     |
164
                                   X                    /                     |
165
                                 /XX__________________/                       |
166
               _____________________                                          |
167
              /                    /                                          |
168
            /____________________/                                            |
169
      ==================================================================
170
171
      The parts marked with X are the parts whose sums are copied into
172
      the temporary buffer.  */
173
174
0
    mp_limb_t tp[MUL_TOOM22_THRESHOLD_LIMIT];
175
0
    mp_limb_t cy;
176
0
    ASSERT (MUL_TOOM22_THRESHOLD <= MUL_TOOM22_THRESHOLD_LIMIT);
177
178
0
    mpn_mul_basecase (prodp, up, MUL_BASECASE_MAX_UN, vp, vn);
179
0
    prodp += MUL_BASECASE_MAX_UN;
180
0
    MPN_COPY (tp, prodp, vn);   /* preserve high triangle */
181
0
    up += MUL_BASECASE_MAX_UN;
182
0
    un -= MUL_BASECASE_MAX_UN;
183
0
    while (un > MUL_BASECASE_MAX_UN)
184
0
      {
185
0
        mpn_mul_basecase (prodp, up, MUL_BASECASE_MAX_UN, vp, vn);
186
0
        cy = mpn_add_n (prodp, prodp, tp, vn); /* add back preserved triangle */
187
0
        mpn_incr_u (prodp + vn, cy);
188
0
        prodp += MUL_BASECASE_MAX_UN;
189
0
        MPN_COPY (tp, prodp, vn);   /* preserve high triangle */
190
0
        up += MUL_BASECASE_MAX_UN;
191
0
        un -= MUL_BASECASE_MAX_UN;
192
0
      }
193
0
    if (un > vn)
194
0
      {
195
0
        mpn_mul_basecase (prodp, up, un, vp, vn);
196
0
      }
197
0
    else
198
0
      {
199
0
        ASSERT (un > 0);
200
0
        mpn_mul_basecase (prodp, vp, vn, up, un);
201
0
      }
202
0
    cy = mpn_add_n (prodp, prodp, tp, vn); /* add back preserved triangle */
203
0
    mpn_incr_u (prodp + vn, cy);
204
0
  }
205
23.0k
    }
206
40.2k
  else if (BELOW_THRESHOLD (vn, MUL_TOOM33_THRESHOLD))
207
38.8k
    {
208
      /* Use ToomX2 variants */
209
38.8k
      mp_ptr scratch;
210
38.8k
      TMP_SDECL; TMP_SMARK;
211
212
38.8k
#define ITCH_TOOMX2 (9 * vn / 2 + GMP_NUMB_BITS * 2)
213
38.8k
      scratch = TMP_SALLOC_LIMBS (ITCH_TOOMX2);
214
38.8k
      ASSERT (mpn_toom22_mul_itch ((5*vn-1)/4, vn) <= ITCH_TOOMX2); /* 5vn/2+ */
215
38.8k
      ASSERT (mpn_toom32_mul_itch ((7*vn-1)/4, vn) <= ITCH_TOOMX2); /* 7vn/6+ */
216
38.8k
      ASSERT (mpn_toom42_mul_itch (3 * vn - 1, vn) <= ITCH_TOOMX2); /* 9vn/2+ */
217
38.8k
#undef ITCH_TOOMX2
218
219
      /* FIXME: This condition (repeated in the loop below) leaves from a vn*vn
220
   square to a (3vn-1)*vn rectangle.  Leaving such a rectangle is hardly
221
   wise; we would get better balance by slightly moving the bound.  We
222
   will sometimes end up with un < vn, like in the X3 arm below.  */
223
38.8k
      if (un >= 3 * vn)
224
257
  {
225
257
    mp_limb_t cy;
226
257
    mp_ptr ws;
227
228
    /* The maximum ws usage is for the mpn_mul result.  */
229
257
    ws = TMP_SALLOC_LIMBS (4 * vn);
230
231
257
    mpn_toom42_mul (prodp, up, 2 * vn, vp, vn, scratch);
232
257
    un -= 2 * vn;
233
257
    up += 2 * vn;
234
257
    prodp += 2 * vn;
235
236
339
    while (un >= 3 * vn)
237
82
      {
238
82
        mpn_toom42_mul (ws, up, 2 * vn, vp, vn, scratch);
239
82
        un -= 2 * vn;
240
82
        up += 2 * vn;
241
82
        cy = mpn_add_n (prodp, prodp, ws, vn);
242
82
        MPN_COPY (prodp + vn, ws + vn, 2 * vn);
243
82
        mpn_incr_u (prodp + vn, cy);
244
82
        prodp += 2 * vn;
245
82
      }
246
247
    /* vn <= un < 3vn */
248
249
257
    if (4 * un < 5 * vn)
250
24
      mpn_toom22_mul (ws, up, un, vp, vn, scratch);
251
233
    else if (4 * un < 7 * vn)
252
93
      mpn_toom32_mul (ws, up, un, vp, vn, scratch);
253
140
    else
254
140
      mpn_toom42_mul (ws, up, un, vp, vn, scratch);
255
256
257
    cy = mpn_add_n (prodp, prodp, ws, vn);
257
257
    MPN_COPY (prodp + vn, ws + vn, un);
258
257
    mpn_incr_u (prodp + vn, cy);
259
257
  }
260
38.5k
      else
261
38.5k
  {
262
38.5k
    if (4 * un < 5 * vn)
263
33.4k
      mpn_toom22_mul (prodp, up, un, vp, vn, scratch);
264
5.12k
    else if (4 * un < 7 * vn)
265
4.93k
      mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
266
193
    else
267
193
      mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
268
38.5k
  }
269
38.8k
      TMP_SFREE;
270
38.8k
    }
271
1.47k
  else if (BELOW_THRESHOLD ((un + vn) >> 1, MUL_FFT_THRESHOLD) ||
272
1.47k
     BELOW_THRESHOLD (3 * vn, MUL_FFT_THRESHOLD))
273
1.47k
    {
274
      /* Handle the largest operands that are not in the FFT range.  The 2nd
275
   condition makes very unbalanced operands avoid the FFT code (except
276
   perhaps as coefficient products of the Toom code.  */
277
278
1.47k
      if (BELOW_THRESHOLD (vn, MUL_TOOM44_THRESHOLD) || !TOOM44_OK (un, vn))
279
1.14k
  {
280
    /* Use ToomX3 variants */
281
1.14k
    mp_ptr scratch;
282
1.14k
    TMP_DECL; TMP_MARK;
283
284
1.14k
#define ITCH_TOOMX3 (4 * vn + GMP_NUMB_BITS)
285
1.14k
    scratch = TMP_ALLOC_LIMBS (ITCH_TOOMX3);
286
1.14k
    ASSERT (mpn_toom33_mul_itch ((7*vn-1)/6, vn) <= ITCH_TOOMX3); /* 7vn/2+ */
287
1.14k
    ASSERT (mpn_toom43_mul_itch ((3*vn-1)/2, vn) <= ITCH_TOOMX3); /* 9vn/4+ */
288
1.14k
    ASSERT (mpn_toom32_mul_itch ((7*vn-1)/4, vn) <= ITCH_TOOMX3); /* 7vn/6+ */
289
1.14k
    ASSERT (mpn_toom53_mul_itch ((11*vn-1)/6, vn) <= ITCH_TOOMX3); /* 11vn/3+ */
290
1.14k
    ASSERT (mpn_toom42_mul_itch ((5*vn-1)/2, vn) <= ITCH_TOOMX3); /* 15vn/4+ */
291
1.14k
    ASSERT (mpn_toom63_mul_itch ((5*vn-1)/2, vn) <= ITCH_TOOMX3); /* 15vn/4+ */
292
1.14k
#undef ITCH_TOOMX3
293
294
1.14k
    if (2 * un >= 5 * vn)
295
56
      {
296
56
        mp_limb_t cy;
297
56
        mp_ptr ws;
298
299
        /* The maximum ws usage is for the mpn_mul result.  */
300
56
        ws = TMP_ALLOC_LIMBS (7 * vn >> 1);
301
302
56
        if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
303
14
    mpn_toom42_mul (prodp, up, 2 * vn, vp, vn, scratch);
304
42
        else
305
42
    mpn_toom63_mul (prodp, up, 2 * vn, vp, vn, scratch);
306
56
        un -= 2 * vn;
307
56
        up += 2 * vn;
308
56
        prodp += 2 * vn;
309
310
209
        while (2 * un >= 5 * vn)  /* un >= 2.5vn */
311
153
    {
312
153
      if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
313
42
        mpn_toom42_mul (ws, up, 2 * vn, vp, vn, scratch);
314
111
      else
315
111
        mpn_toom63_mul (ws, up, 2 * vn, vp, vn, scratch);
316
153
      un -= 2 * vn;
317
153
      up += 2 * vn;
318
153
      cy = mpn_add_n (prodp, prodp, ws, vn);
319
153
      MPN_COPY (prodp + vn, ws + vn, 2 * vn);
320
153
      mpn_incr_u (prodp + vn, cy);
321
153
      prodp += 2 * vn;
322
153
    }
323
324
        /* vn / 2 <= un < 2.5vn */
325
326
56
        if (un < vn)
327
13
    mpn_mul (ws, vp, vn, up, un);
328
43
        else
329
43
    mpn_mul (ws, up, un, vp, vn);
330
331
56
        cy = mpn_add_n (prodp, prodp, ws, vn);
332
56
        MPN_COPY (prodp + vn, ws + vn, un);
333
56
        mpn_incr_u (prodp + vn, cy);
334
56
      }
335
1.08k
    else
336
1.08k
      {
337
1.08k
        if (6 * un < 7 * vn)
338
914
    mpn_toom33_mul (prodp, up, un, vp, vn, scratch);
339
174
        else if (2 * un < 3 * vn)
340
105
    {
341
105
      if (BELOW_THRESHOLD (vn, MUL_TOOM32_TO_TOOM43_THRESHOLD))
342
1
        mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
343
104
      else
344
104
        mpn_toom43_mul (prodp, up, un, vp, vn, scratch);
345
105
    }
346
69
        else if (6 * un < 11 * vn)
347
67
    {
348
67
      if (4 * un < 7 * vn)
349
55
        {
350
55
          if (BELOW_THRESHOLD (vn, MUL_TOOM32_TO_TOOM53_THRESHOLD))
351
32
      mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
352
23
          else
353
23
      mpn_toom53_mul (prodp, up, un, vp, vn, scratch);
354
55
        }
355
12
      else
356
12
        {
357
12
          if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM53_THRESHOLD))
358
3
      mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
359
9
          else
360
9
      mpn_toom53_mul (prodp, up, un, vp, vn, scratch);
361
12
        }
362
67
    }
363
2
        else
364
2
    {
365
2
      if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
366
1
        mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
367
1
      else
368
1
        mpn_toom63_mul (prodp, up, un, vp, vn, scratch);
369
2
    }
370
1.08k
      }
371
1.14k
    TMP_FREE;
372
1.14k
  }
373
333
      else
374
333
  {
375
333
    mp_ptr scratch;
376
333
    TMP_DECL; TMP_MARK;
377
378
333
    if (BELOW_THRESHOLD (vn, MUL_TOOM6H_THRESHOLD))
379
110
      {
380
110
        scratch = TMP_SALLOC_LIMBS (mpn_toom44_mul_itch (un, vn));
381
110
        mpn_toom44_mul (prodp, up, un, vp, vn, scratch);
382
110
      }
383
223
    else if (BELOW_THRESHOLD (vn, MUL_TOOM8H_THRESHOLD))
384
77
      {
385
77
        scratch = TMP_SALLOC_LIMBS (mpn_toom6h_mul_itch (un, vn));
386
77
        mpn_toom6h_mul (prodp, up, un, vp, vn, scratch);
387
77
      }
388
146
    else
389
146
      {
390
146
        scratch = TMP_ALLOC_LIMBS (mpn_toom8h_mul_itch (un, vn));
391
146
        mpn_toom8h_mul (prodp, up, un, vp, vn, scratch);
392
146
      }
393
333
    TMP_FREE;
394
333
  }
395
1.47k
    }
396
0
  else
397
0
    {
398
0
      if (un >= 8 * vn)
399
0
  {
400
0
    mp_limb_t cy;
401
0
    mp_ptr ws;
402
0
    TMP_DECL; TMP_MARK;
403
404
    /* The maximum ws usage is for the mpn_mul result.  */
405
0
    ws = TMP_BALLOC_LIMBS (9 * vn >> 1);
406
407
0
    mpn_fft_mul (prodp, up, 3 * vn, vp, vn);
408
0
    un -= 3 * vn;
409
0
    up += 3 * vn;
410
0
    prodp += 3 * vn;
411
412
0
    while (2 * un >= 7 * vn) /* un >= 3.5vn  */
413
0
      {
414
0
        mpn_fft_mul (ws, up, 3 * vn, vp, vn);
415
0
        un -= 3 * vn;
416
0
        up += 3 * vn;
417
0
        cy = mpn_add_n (prodp, prodp, ws, vn);
418
0
        MPN_COPY (prodp + vn, ws + vn, 3 * vn);
419
0
        mpn_incr_u (prodp + vn, cy);
420
0
        prodp += 3 * vn;
421
0
      }
422
423
    /* vn / 2 <= un < 3.5vn */
424
425
0
    if (un < vn)
426
0
      mpn_mul (ws, vp, vn, up, un);
427
0
    else
428
0
      mpn_mul (ws, up, un, vp, vn);
429
430
0
    cy = mpn_add_n (prodp, prodp, ws, vn);
431
0
    MPN_COPY (prodp + vn, ws + vn, un);
432
0
    mpn_incr_u (prodp + vn, cy);
433
434
0
    TMP_FREE;
435
0
  }
436
0
      else
437
0
  mpn_fft_mul (prodp, up, un, vp, vn);
438
0
    }
439
440
202k
  return prodp[un + vn - 1]; /* historic */
441
202k
}