Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpn/toom44_mul.c
Line
Count
Source (jump to first uncovered line)
1
/* mpn_toom44_mul -- Multiply {ap,an} and {bp,bn} where an and bn are close in
2
   size.  Or more accurately, bn <= an < (4/3)bn.
3
4
   Contributed to the GNU project by Torbjorn Granlund and Marco Bodrato.
5
6
   THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
7
   SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
8
   GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
9
10
Copyright 2006-2008, 2013 Free Software 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 either:
16
17
  * the GNU Lesser General Public License as published by the Free
18
    Software Foundation; either version 3 of the License, or (at your
19
    option) any later version.
20
21
or
22
23
  * the GNU General Public License as published by the Free Software
24
    Foundation; either version 2 of the License, or (at your option) any
25
    later version.
26
27
or both in parallel, as here.
28
29
The GNU MP Library is distributed in the hope that it will be useful, but
30
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
31
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
32
for more details.
33
34
You should have received copies of the GNU General Public License and the
35
GNU Lesser General Public License along with the GNU MP Library.  If not,
36
see https://www.gnu.org/licenses/.  */
37
38
39
#include "gmp-impl.h"
40
41
/* Evaluate in: 0, +1, -1, +2, -2, 1/2, +inf
42
43
  <-s--><--n--><--n--><--n-->
44
   ____ ______ ______ ______
45
  |_a3_|___a2_|___a1_|___a0_|
46
   |b3_|___b2_|___b1_|___b0_|
47
   <-t-><--n--><--n--><--n-->
48
49
  v0  =   a0             *  b0              #    A(0)*B(0)
50
  v1  = ( a0+ a1+ a2+ a3)*( b0+ b1+ b2+ b3) #    A(1)*B(1)      ah  <= 3   bh  <= 3
51
  vm1 = ( a0- a1+ a2- a3)*( b0- b1+ b2- b3) #   A(-1)*B(-1)    |ah| <= 1  |bh| <= 1
52
  v2  = ( a0+2a1+4a2+8a3)*( b0+2b1+4b2+8b3) #    A(2)*B(2)      ah  <= 14  bh  <= 14
53
  vm2 = ( a0-2a1+4a2-8a3)*( b0-2b1+4b2-8b3) #   A(-2)*B(-2)    |ah| <= 9  |bh| <= 9
54
  vh  = (8a0+4a1+2a2+ a3)*(8b0+4b1+2b2+ b3) #  A(1/2)*B(1/2)    ah  <= 14  bh  <= 14
55
  vinf=               a3 *          b2      #  A(inf)*B(inf)
56
*/
57
58
#if TUNE_PROGRAM_BUILD
59
#define MAYBE_mul_basecase 1
60
#define MAYBE_mul_toom22   1
61
#define MAYBE_mul_toom44   1
62
#else
63
#define MAYBE_mul_basecase            \
64
8.89k
  (MUL_TOOM44_THRESHOLD < 4 * MUL_TOOM22_THRESHOLD)
65
#define MAYBE_mul_toom22            \
66
8.89k
  (MUL_TOOM44_THRESHOLD < 4 * MUL_TOOM33_THRESHOLD)
67
#define MAYBE_mul_toom44            \
68
0
  (MUL_TOOM6H_THRESHOLD >= 4 * MUL_TOOM44_THRESHOLD)
69
#endif
70
71
#define TOOM44_MUL_N_REC(p, a, b, n, ws)        \
72
4.44k
  do {                 \
73
4.44k
    if (MAYBE_mul_basecase            \
74
4.44k
  && BELOW_THRESHOLD (n, MUL_TOOM22_THRESHOLD))     \
75
4.44k
      mpn_mul_basecase (p, a, n, b, n);         \
76
4.44k
    else if (MAYBE_mul_toom22            \
77
4.44k
       && BELOW_THRESHOLD (n, MUL_TOOM33_THRESHOLD))   \
78
4.44k
      mpn_toom22_mul (p, a, n, b, n, ws);       \
79
4.44k
    else if (! MAYBE_mul_toom44            \
80
0
       || BELOW_THRESHOLD (n, MUL_TOOM44_THRESHOLD))   \
81
0
      mpn_toom33_mul (p, a, n, b, n, ws);       \
82
0
    else                \
83
0
      mpn_toom44_mul (p, a, n, b, n, ws);       \
84
4.44k
  } while (0)
85
86
/* Use of scratch space. In the product area, we store
87
88
      ___________________
89
     |vinf|____|_v1_|_v0_|
90
      s+t  2n-1 2n+1  2n
91
92
   The other recursive products, vm1, v2, vm2, vh are stored in the
93
   scratch area. When computing them, we use the product area for
94
   intermediate values.
95
96
   Next, we compute v1. We can store the intermediate factors at v0
97
   and at vh + 2n + 2.
98
99
   Finally, for v0 and vinf, factors are parts of the input operands,
100
   and we need scratch space only for the recursive multiplication.
101
102
   In all, if S(an) is the scratch need, the needed space is bounded by
103
104
     S(an) <= 4 (2*ceil(an/4) + 1) + 1 + S(ceil(an/4) + 1)
105
106
   which should give S(n) = 8 n/3 + c log(n) for some constant c.
107
*/
108
109
void
110
mpn_toom44_mul (mp_ptr pp,
111
    mp_srcptr ap, mp_size_t an,
112
    mp_srcptr bp, mp_size_t bn,
113
    mp_ptr scratch)
114
651
{
115
651
  mp_size_t n, s, t;
116
651
  mp_limb_t cy;
117
651
  enum toom7_flags flags;
118
119
651
#define a0  ap
120
651
#define a1  (ap + n)
121
651
#define a2  (ap + 2*n)
122
761
#define a3  (ap + 3*n)
123
651
#define b0  bp
124
651
#define b1  (bp + n)
125
651
#define b2  (bp + 2*n)
126
761
#define b3  (bp + 3*n)
127
128
651
  ASSERT (an >= bn);
129
130
651
  n = (an + 3) >> 2;
131
132
651
  s = an - 3 * n;
133
651
  t = bn - 3 * n;
134
135
651
  ASSERT (0 < s && s <= n);
136
651
  ASSERT (0 < t && t <= n);
137
651
  ASSERT (s >= t);
138
139
  /* NOTE: The multiplications to v2, vm2, vh and vm1 overwrites the
140
   * following limb, so these must be computed in order, and we need a
141
   * one limb gap to tp. */
142
651
#define v0    pp        /* 2n */
143
651
#define v1    (pp + 2 * n)      /* 2n+1 */
144
651
#define vinf  (pp + 6 * n)      /* s+t */
145
651
#define v2    scratch        /* 2n+1 */
146
651
#define vm2   (scratch + 2 * n + 1)    /* 2n+1 */
147
651
#define vh    (scratch + 4 * n + 2)    /* 2n+1 */
148
1.30k
#define vm1   (scratch + 6 * n + 3)    /* 2n+1 */
149
3.25k
#define tp (scratch + 8*n + 5)
150
151
  /* apx and bpx must not overlap with v1 */
152
6.39k
#define apx   pp        /* n+1 */
153
1.30k
#define amx   (pp + n + 1)      /* n+1 */
154
1.30k
#define bmx   (pp + 2*n + 2)      /* n+1 */
155
6.40k
#define bpx   (pp + 4*n + 2)      /* n+1 */
156
157
  /* Total scratch need: 8*n + 5 + scratch for recursive calls. This
158
     gives roughly 32 n/3 + log term. */
159
160
  /* Compute apx = a0 + 2 a1 + 4 a2 + 8 a3 and amx = a0 - 2 a1 + 4 a2 - 8 a3.  */
161
651
  flags = (enum toom7_flags) (toom7_w1_neg & mpn_toom_eval_dgr3_pm2 (apx, amx, ap, n, s, tp));
162
163
  /* Compute bpx = b0 + 2 b1 + 4 b2 + 8 b3 and bmx = b0 - 2 b1 + 4 b2 - 8 b3.  */
164
651
  flags = (enum toom7_flags) (flags ^ (toom7_w1_neg & mpn_toom_eval_dgr3_pm2 (bpx, bmx, bp, n, t, tp)));
165
166
651
  TOOM44_MUL_N_REC (v2, apx, bpx, n + 1, tp); /* v2,  2n+1 limbs */
167
651
  TOOM44_MUL_N_REC (vm2, amx, bmx, n + 1, tp);  /* vm2,  2n+1 limbs */
168
169
  /* Compute apx = 8 a0 + 4 a1 + 2 a2 + a3 = (((2*a0 + a1) * 2 + a2) * 2 + a3 */
170
651
#if HAVE_NATIVE_mpn_addlsh1_n
171
651
  cy = mpn_addlsh1_n (apx, a1, a0, n);
172
651
  cy = 2*cy + mpn_addlsh1_n (apx, a2, apx, n);
173
651
  if (s < n)
174
595
    {
175
595
      mp_limb_t cy2;
176
595
      cy2 = mpn_addlsh1_n (apx, a3, apx, s);
177
595
      apx[n] = 2*cy + mpn_lshift (apx + s, apx + s, n - s, 1);
178
595
      MPN_INCR_U (apx + s, n+1-s, cy2);
179
595
    }
180
56
  else
181
56
    apx[n] = 2*cy + mpn_addlsh1_n (apx, a3, apx, n);
182
#else
183
  cy = mpn_lshift (apx, a0, n, 1);
184
  cy += mpn_add_n (apx, apx, a1, n);
185
  cy = 2*cy + mpn_lshift (apx, apx, n, 1);
186
  cy += mpn_add_n (apx, apx, a2, n);
187
  cy = 2*cy + mpn_lshift (apx, apx, n, 1);
188
  apx[n] = cy + mpn_add (apx, apx, n, a3, s);
189
#endif
190
191
  /* Compute bpx = 8 b0 + 4 b1 + 2 b2 + b3 = (((2*b0 + b1) * 2 + b2) * 2 + b3 */
192
651
#if HAVE_NATIVE_mpn_addlsh1_n
193
651
  cy = mpn_addlsh1_n (bpx, b1, b0, n);
194
651
  cy = 2*cy + mpn_addlsh1_n (bpx, b2, bpx, n);
195
651
  if (t < n)
196
599
    {
197
599
      mp_limb_t cy2;
198
599
      cy2 = mpn_addlsh1_n (bpx, b3, bpx, t);
199
599
      bpx[n] = 2*cy + mpn_lshift (bpx + t, bpx + t, n - t, 1);
200
599
      MPN_INCR_U (bpx + t, n+1-t, cy2);
201
599
    }
202
52
  else
203
52
    bpx[n] = 2*cy + mpn_addlsh1_n (bpx, b3, bpx, n);
204
#else
205
  cy = mpn_lshift (bpx, b0, n, 1);
206
  cy += mpn_add_n (bpx, bpx, b1, n);
207
  cy = 2*cy + mpn_lshift (bpx, bpx, n, 1);
208
  cy += mpn_add_n (bpx, bpx, b2, n);
209
  cy = 2*cy + mpn_lshift (bpx, bpx, n, 1);
210
  bpx[n] = cy + mpn_add (bpx, bpx, n, b3, t);
211
#endif
212
213
651
  ASSERT (apx[n] < 15);
214
651
  ASSERT (bpx[n] < 15);
215
216
651
  TOOM44_MUL_N_REC (vh, apx, bpx, n + 1, tp); /* vh,  2n+1 limbs */
217
218
  /* Compute apx = a0 + a1 + a2 + a3 and amx = a0 - a1 + a2 - a3.  */
219
651
  flags = (enum toom7_flags) (flags | (toom7_w3_neg & mpn_toom_eval_dgr3_pm1 (apx, amx, ap, n, s, tp)));
220
221
  /* Compute bpx = b0 + b1 + b2 + b3 and bmx = b0 - b1 + b2 - b3.  */
222
651
  flags = (enum toom7_flags) (flags ^ (toom7_w3_neg & mpn_toom_eval_dgr3_pm1 (bpx, bmx, bp, n, t, tp)));
223
224
651
  ASSERT (amx[n] <= 1);
225
651
  ASSERT (bmx[n] <= 1);
226
227
651
  vm1 [2 * n] = 0;
228
651
  TOOM44_MUL_N_REC (vm1, amx, bmx, n + (bmx[n] | amx[n]), tp);  /* vm1,  2n+1 limbs */
229
  /* Clobbers amx, bmx. */
230
651
  TOOM44_MUL_N_REC (v1, apx, bpx, n + 1, tp); /* v1,  2n+1 limbs */
231
232
651
  TOOM44_MUL_N_REC (v0, a0, b0, n, tp);
233
651
  if (s > t)
234
110
    mpn_mul (vinf, a3, s, b3, t);
235
541
  else
236
541
    TOOM44_MUL_N_REC (vinf, a3, b3, s, tp); /* vinf, s+t limbs */
237
238
651
  mpn_toom_interpolate_7pts (pp, n, flags, vm2, vm1, v2, vh, s + t, tp);
239
651
}