Coverage Report

Created: 2024-11-21 06:47

/src/libgmp/mpn/toom42_mul.c
Line
Count
Source
1
/* mpn_toom42_mul -- Multiply {ap,an} and {bp,bn} where an is nominally twice
2
   as large as bn.  Or more accurately, (3/2)bn < an < 4bn.
3
4
   Contributed to the GNU project by Torbjorn Granlund.
5
   Additional improvements by Marco Bodrato.
6
7
   The idea of applying toom to unbalanced multiplication is due to Marco
8
   Bodrato and Alberto Zanoni.
9
10
   THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
11
   SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
12
   GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
13
14
Copyright 2006-2008, 2012, 2014 Free Software Foundation, Inc.
15
16
This file is part of the GNU MP Library.
17
18
The GNU MP Library is free software; you can redistribute it and/or modify
19
it under the terms of either:
20
21
  * the GNU Lesser General Public License as published by the Free
22
    Software Foundation; either version 3 of the License, or (at your
23
    option) any later version.
24
25
or
26
27
  * the GNU General Public License as published by the Free Software
28
    Foundation; either version 2 of the License, or (at your option) any
29
    later version.
30
31
or both in parallel, as here.
32
33
The GNU MP Library is distributed in the hope that it will be useful, but
34
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
35
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
36
for more details.
37
38
You should have received copies of the GNU General Public License and the
39
GNU Lesser General Public License along with the GNU MP Library.  If not,
40
see https://www.gnu.org/licenses/.  */
41
42
43
#include "gmp-impl.h"
44
45
/* Evaluate in: -1, 0, +1, +2, +inf
46
47
  <-s-><--n--><--n--><--n-->
48
   ___ ______ ______ ______
49
  |a3_|___a2_|___a1_|___a0_|
50
         |_b1_|___b0_|
51
         <-t--><--n-->
52
53
  v0  =  a0             * b0      #   A(0)*B(0)
54
  v1  = (a0+ a1+ a2+ a3)*(b0+ b1) #   A(1)*B(1)      ah  <= 3  bh <= 1
55
  vm1 = (a0- a1+ a2- a3)*(b0- b1) #  A(-1)*B(-1)    |ah| <= 1  bh  = 0
56
  v2  = (a0+2a1+4a2+8a3)*(b0+2b1) #   A(2)*B(2)      ah  <= 14 bh <= 2
57
  vinf=              a3 *     b1  # A(inf)*B(inf)
58
*/
59
60
#define TOOM42_MUL_N_REC(p, a, b, n, ws)        \
61
3.03k
  do {                 \
62
3.03k
    mpn_mul_n (p, a, b, n);            \
63
3.03k
  } while (0)
64
65
void
66
mpn_toom42_mul (mp_ptr pp,
67
    mp_srcptr ap, mp_size_t an,
68
    mp_srcptr bp, mp_size_t bn,
69
    mp_ptr scratch)
70
758
{
71
758
  mp_size_t n, s, t;
72
758
  int vm1_neg;
73
758
  mp_limb_t cy, vinf0;
74
758
  mp_ptr a0_a2;
75
758
  mp_ptr as1, asm1, as2;
76
758
  mp_ptr bs1, bsm1, bs2;
77
758
  mp_ptr tmp;
78
758
  TMP_DECL;
79
80
758
#define a0  ap
81
758
#define a1  (ap + n)
82
1.48k
#define a2  (ap + 2*n)
83
1.51k
#define a3  (ap + 3*n)
84
2.49k
#define b0  bp
85
3.41k
#define b1  (bp + n)
86
87
758
  n = an >= 2 * bn ? (an + 3) >> 2 : (bn + 1) >> 1;
88
89
758
  s = an - 3 * n;
90
758
  t = bn - n;
91
92
758
  ASSERT (0 < s && s <= n);
93
758
  ASSERT (0 < t && t <= n);
94
95
758
  TMP_MARK;
96
97
758
  tmp = TMP_ALLOC_LIMBS (6 * n + 5);
98
758
  as1  = tmp; tmp += n + 1;
99
758
  asm1 = tmp; tmp += n + 1;
100
758
  as2  = tmp; tmp += n + 1;
101
758
  bs1  = tmp; tmp += n + 1;
102
758
  bsm1 = tmp; tmp += n;
103
758
  bs2  = tmp; tmp += n + 1;
104
105
758
  a0_a2 = pp;
106
107
  /* Compute as1 and asm1.  */
108
758
  vm1_neg = mpn_toom_eval_dgr3_pm1 (as1, asm1, ap, n, s, a0_a2) & 1;
109
110
  /* Compute as2.  */
111
758
#if HAVE_NATIVE_mpn_addlsh1_n
112
758
  cy  = mpn_addlsh1_n (as2, a2, a3, s);
113
758
  if (s != n)
114
730
    cy = mpn_add_1 (as2 + s, a2 + s, n - s, cy);
115
758
  cy = 2 * cy + mpn_addlsh1_n (as2, a1, as2, n);
116
758
  cy = 2 * cy + mpn_addlsh1_n (as2, a0, as2, n);
117
#else
118
  cy  = mpn_lshift (as2, a3, s, 1);
119
  cy += mpn_add_n (as2, a2, as2, s);
120
  if (s != n)
121
    cy = mpn_add_1 (as2 + s, a2 + s, n - s, cy);
122
  cy = 2 * cy + mpn_lshift (as2, as2, n, 1);
123
  cy += mpn_add_n (as2, a1, as2, n);
124
  cy = 2 * cy + mpn_lshift (as2, as2, n, 1);
125
  cy += mpn_add_n (as2, a0, as2, n);
126
#endif
127
758
  as2[n] = cy;
128
129
  /* Compute bs1 and bsm1.  */
130
758
  if (t == n)
131
169
    {
132
#if HAVE_NATIVE_mpn_add_n_sub_n
133
      if (mpn_cmp (b0, b1, n) < 0)
134
  {
135
    cy = mpn_add_n_sub_n (bs1, bsm1, b1, b0, n);
136
    vm1_neg ^= 1;
137
  }
138
      else
139
  {
140
    cy = mpn_add_n_sub_n (bs1, bsm1, b0, b1, n);
141
  }
142
      bs1[n] = cy >> 1;
143
#else
144
169
      bs1[n] = mpn_add_n (bs1, b0, b1, n);
145
146
169
      if (mpn_cmp (b0, b1, n) < 0)
147
14
  {
148
14
    mpn_sub_n (bsm1, b1, b0, n);
149
14
    vm1_neg ^= 1;
150
14
  }
151
155
      else
152
155
  {
153
155
    mpn_sub_n (bsm1, b0, b1, n);
154
155
  }
155
169
#endif
156
169
    }
157
589
  else
158
589
    {
159
589
      bs1[n] = mpn_add (bs1, b0, n, b1, t);
160
161
589
      if (mpn_zero_p (b0 + t, n - t) && mpn_cmp (b0, b1, t) < 0)
162
52
  {
163
52
    mpn_sub_n (bsm1, b1, b0, t);
164
52
    MPN_ZERO (bsm1 + t, n - t);
165
52
    vm1_neg ^= 1;
166
52
  }
167
537
      else
168
537
  {
169
537
    mpn_sub (bsm1, b0, n, b1, t);
170
537
  }
171
589
    }
172
173
  /* Compute bs2, recycling bs1. bs2=bs1+b1  */
174
758
  mpn_add (bs2, bs1, n + 1, b1, t);
175
176
758
  ASSERT (as1[n] <= 3);
177
758
  ASSERT (bs1[n] <= 1);
178
758
  ASSERT (asm1[n] <= 1);
179
/*ASSERT (bsm1[n] == 0);*/
180
758
  ASSERT (as2[n] <= 14);
181
758
  ASSERT (bs2[n] <= 2);
182
183
758
#define v0    pp        /* 2n */
184
1.75k
#define v1    (pp + 2 * n)      /* 2n+1 */
185
1.51k
#define vinf  (pp + 4 * n)      /* s+t */
186
1.75k
#define vm1   scratch        /* 2n+1 */
187
758
#define v2    (scratch + 2 * n + 1)    /* 2n+2 */
188
758
#define scratch_out scratch + 4 * n + 4 /* Currently unused. */
189
190
  /* vm1, 2n+1 limbs */
191
758
  TOOM42_MUL_N_REC (vm1, asm1, bsm1, n, scratch_out);
192
758
  cy = 0;
193
758
  if (asm1[n] != 0)
194
118
    cy = mpn_add_n (vm1 + n, vm1 + n, bsm1, n);
195
758
  vm1[2 * n] = cy;
196
197
758
  TOOM42_MUL_N_REC (v2, as2, bs2, n + 1, scratch_out);  /* v2, 2n+1 limbs */
198
199
  /* vinf, s+t limbs */
200
758
  if (s > t)  mpn_mul (vinf, a3, s, b1, t);
201
622
  else        mpn_mul (vinf, b1, t, a3, s);
202
203
758
  vinf0 = vinf[0];        /* v1 overlaps with this */
204
205
  /* v1, 2n+1 limbs */
206
758
  TOOM42_MUL_N_REC (v1, as1, bs1, n, scratch_out);
207
758
  if (as1[n] == 1)
208
478
    {
209
478
      cy = mpn_add_n (v1 + n, v1 + n, bs1, n);
210
478
    }
211
280
  else if (as1[n] == 2)
212
122
    {
213
122
#if HAVE_NATIVE_mpn_addlsh1_n_ip1
214
122
      cy = mpn_addlsh1_n_ip1 (v1 + n, bs1, n);
215
#else
216
      cy = mpn_addmul_1 (v1 + n, bs1, n, CNST_LIMB(2));
217
#endif
218
122
    }
219
158
  else if (as1[n] == 3)
220
3
    {
221
3
      cy = mpn_addmul_1 (v1 + n, bs1, n, CNST_LIMB(3));
222
3
    }
223
155
  else
224
155
    cy = 0;
225
758
  if (bs1[n] != 0)
226
21
    cy += as1[n] + mpn_add_n (v1 + n, v1 + n, as1, n);
227
758
  v1[2 * n] = cy;
228
229
758
  TOOM42_MUL_N_REC (v0, ap, bp, n, scratch_out);  /* v0, 2n limbs */
230
231
758
  mpn_toom_interpolate_5pts (pp, v2, vm1, n, s + t, vm1_neg, vinf0);
232
233
758
  TMP_FREE;
234
758
}