Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpn/mulmod_bnm1.c
Line
Count
Source (jump to first uncovered line)
1
/* mulmod_bnm1.c -- multiplication mod B^n-1.
2
3
   Contributed to the GNU project by Niels Möller, Torbjorn Granlund and
4
   Marco Bodrato.
5
6
   THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
7
   SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
8
   GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
9
10
Copyright 2009, 2010, 2012, 2013, 2020, 2022 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
#include "longlong.h"
41
42
/* Inputs are {ap,rn} and {bp,rn}; output is {rp,rn}, computation is
43
   mod B^rn - 1, and values are semi-normalised; zero is represented
44
   as either 0 or B^n - 1.  Needs a scratch of 2rn limbs at tp.
45
   tp==rp is allowed. */
46
void
47
mpn_bc_mulmod_bnm1 (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t rn,
48
        mp_ptr tp)
49
2.22M
{
50
2.22M
  mp_limb_t cy;
51
52
2.22M
  ASSERT (0 < rn);
53
54
2.22M
  mpn_mul_n (tp, ap, bp, rn);
55
2.22M
  cy = mpn_add_n (rp, tp, tp + rn, rn);
56
  /* If cy == 1, then the value of rp is at most B^rn - 2, so there can
57
   * be no overflow when adding in the carry. */
58
2.22M
  MPN_INCR_U (rp, rn, cy);
59
2.22M
}
60
61
62
/* Inputs are {ap,rn+1} and {bp,rn+1}; output is {rp,rn+1}, in
63
   normalised representation, computation is mod B^rn + 1. Needs
64
   a scratch area of 2rn limbs at tp; tp == rp is allowed.
65
   Output is normalised. */
66
static void
67
mpn_bc_mulmod_bnp1 (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t rn,
68
        mp_ptr tp)
69
6.42M
{
70
6.42M
  mp_limb_t cy;
71
6.42M
  unsigned k;
72
73
6.42M
  ASSERT (0 < rn);
74
75
6.42M
  if (UNLIKELY (ap[rn] | bp [rn]))
76
0
    {
77
0
      if (ap[rn])
78
0
  cy = bp [rn] + mpn_neg (rp, bp, rn);
79
0
      else /* ap[rn] == 0 */
80
0
  cy = mpn_neg (rp, ap, rn);
81
0
    }
82
6.42M
  else if (MPN_MULMOD_BKNP1_USABLE (rn, k, MUL_FFT_MODF_THRESHOLD))
83
2.11M
    {
84
2.11M
      mp_size_t n_k = rn / k;
85
2.11M
      TMP_DECL;
86
87
2.11M
      TMP_MARK;
88
2.11M
      mpn_mulmod_bknp1 (rp, ap, bp, n_k, k,
89
2.11M
                       TMP_ALLOC_LIMBS (mpn_mulmod_bknp1_itch (rn)));
90
2.11M
      TMP_FREE;
91
2.11M
      return;
92
2.11M
    }
93
4.30M
  else
94
4.30M
    {
95
4.30M
      mpn_mul_n (tp, ap, bp, rn);
96
4.30M
      cy = mpn_sub_n (rp, tp, tp + rn, rn);
97
4.30M
    }
98
4.30M
  rp[rn] = 0;
99
4.30M
  MPN_INCR_U (rp, rn + 1, cy);
100
4.30M
}
101
102
103
/* Computes {rp,MIN(rn,an+bn)} <- {ap,an}*{bp,bn} Mod(B^rn-1)
104
 *
105
 * The result is expected to be ZERO if and only if one of the operand
106
 * already is. Otherwise the class [0] Mod(B^rn-1) is represented by
107
 * B^rn-1. This should not be a problem if mulmod_bnm1 is used to
108
 * combine results and obtain a natural number when one knows in
109
 * advance that the final value is less than (B^rn-1).
110
 * Moreover it should not be a problem if mulmod_bnm1 is used to
111
 * compute the full product with an+bn <= rn, because this condition
112
 * implies (B^an-1)(B^bn-1) < (B^rn-1) .
113
 *
114
 * Requires 0 < bn <= an <= rn and an + bn > rn/2
115
 * Scratch need: rn + (need for recursive call OR rn + 4). This gives
116
 *
117
 * S(n) <= rn + MAX (rn + 4, S(n/2)) <= 2rn + 4
118
 */
119
void
120
mpn_mulmod_bnm1 (mp_ptr rp, mp_size_t rn, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn, mp_ptr tp)
121
8.64M
{
122
8.64M
  ASSERT (0 < bn);
123
8.64M
  ASSERT (bn <= an);
124
8.64M
  ASSERT (an <= rn);
125
126
8.64M
  if ((rn & 1) != 0 || BELOW_THRESHOLD (rn, MULMOD_BNM1_THRESHOLD))
127
2.22M
    {
128
2.22M
      if (UNLIKELY (bn < rn))
129
0
  {
130
0
    if (UNLIKELY (an + bn <= rn))
131
0
      {
132
0
        mpn_mul (rp, ap, an, bp, bn);
133
0
      }
134
0
    else
135
0
      {
136
0
        mp_limb_t cy;
137
0
        mpn_mul (tp, ap, an, bp, bn);
138
0
        cy = mpn_add (rp, tp, rn, tp + rn, an + bn - rn);
139
0
        MPN_INCR_U (rp, rn, cy);
140
0
      }
141
0
  }
142
2.22M
      else
143
2.22M
  mpn_bc_mulmod_bnm1 (rp, ap, bp, rn, tp);
144
2.22M
    }
145
6.42M
  else
146
6.42M
    {
147
6.42M
      mp_size_t n;
148
6.42M
      mp_limb_t cy;
149
6.42M
      mp_limb_t hi;
150
151
6.42M
      n = rn >> 1;
152
153
      /* We need at least an + bn >= n, to be able to fit one of the
154
   recursive products at rp. Requiring strict inequality makes
155
   the code slightly simpler. If desired, we could avoid this
156
   restriction by initially halving rn as long as rn is even and
157
   an + bn <= rn/2. */
158
159
6.42M
      ASSERT (an + bn > n);
160
161
      /* Compute xm = a*b mod (B^n - 1), xp = a*b mod (B^n + 1)
162
   and crt together as
163
164
   x = -xp * B^n + (B^n + 1) * [ (xp + xm)/2 mod (B^n-1)]
165
      */
166
167
12.8M
#define a0 ap
168
12.8M
#define a1 (ap + n)
169
25.6M
#define b0 bp
170
12.8M
#define b1 (bp + n)
171
172
57.7M
#define xp  tp  /* 2n + 2 */
173
      /* am1  maybe in {xp, n} */
174
      /* bm1  maybe in {xp + n, n} */
175
38.5M
#define sp1 (tp + 2*n + 2)
176
      /* ap1  maybe in {sp1, n + 1} */
177
      /* bp1  maybe in {sp1 + n + 1, n + 1} */
178
179
6.42M
      {
180
6.42M
  mp_srcptr am1, bm1;
181
6.42M
  mp_size_t anm, bnm;
182
6.42M
  mp_ptr so;
183
184
6.42M
  bm1 = b0;
185
6.42M
  bnm = bn;
186
6.42M
  if (LIKELY (an > n))
187
6.42M
    {
188
6.42M
      am1 = xp;
189
6.42M
      cy = mpn_add (xp, a0, n, a1, an - n);
190
6.42M
      MPN_INCR_U (xp, n, cy);
191
6.42M
      anm = n;
192
6.42M
      so = xp + n;
193
6.42M
      if (LIKELY (bn > n))
194
6.42M
        {
195
6.42M
    bm1 = so;
196
6.42M
    cy = mpn_add (so, b0, n, b1, bn - n);
197
6.42M
    MPN_INCR_U (so, n, cy);
198
6.42M
    bnm = n;
199
6.42M
    so += n;
200
6.42M
        }
201
6.42M
    }
202
0
  else
203
0
    {
204
0
      so = xp;
205
0
      am1 = a0;
206
0
      anm = an;
207
0
    }
208
209
6.42M
  mpn_mulmod_bnm1 (rp, n, am1, anm, bm1, bnm, so);
210
6.42M
      }
211
212
0
      {
213
6.42M
  int       k;
214
6.42M
  mp_srcptr ap1, bp1;
215
6.42M
  mp_size_t anp, bnp;
216
217
6.42M
  bp1 = b0;
218
6.42M
  bnp = bn;
219
6.42M
  if (LIKELY (an > n)) {
220
6.42M
    ap1 = sp1;
221
6.42M
    cy = mpn_sub (sp1, a0, n, a1, an - n);
222
6.42M
    sp1[n] = 0;
223
6.42M
    MPN_INCR_U (sp1, n + 1, cy);
224
6.42M
    anp = n + ap1[n];
225
6.42M
    if (LIKELY (bn > n)) {
226
6.42M
      bp1 = sp1 + n + 1;
227
6.42M
      cy = mpn_sub (sp1 + n + 1, b0, n, b1, bn - n);
228
6.42M
      sp1[2*n+1] = 0;
229
6.42M
      MPN_INCR_U (sp1 + n + 1, n + 1, cy);
230
6.42M
      bnp = n + bp1[n];
231
6.42M
    }
232
6.42M
  } else {
233
0
    ap1 = a0;
234
0
    anp = an;
235
0
  }
236
237
6.42M
  if (BELOW_THRESHOLD (n, MUL_FFT_MODF_THRESHOLD))
238
6.42M
    k=0;
239
0
  else
240
0
    {
241
0
      int mask;
242
0
      k = mpn_fft_best_k (n, 0);
243
0
      mask = (1<<k) - 1;
244
0
      while (n & mask) {k--; mask >>=1;};
245
0
    }
246
6.42M
  if (k >= FFT_FIRST_K)
247
0
    xp[n] = mpn_mul_fft (xp, n, ap1, anp, bp1, bnp, k);
248
6.42M
  else if (UNLIKELY (bp1 == b0))
249
0
    {
250
0
      ASSERT (anp + bnp <= 2*n+1);
251
0
      ASSERT (anp + bnp > n);
252
0
      ASSERT (anp >= bnp);
253
0
      mpn_mul (xp, ap1, anp, bp1, bnp);
254
0
      anp = anp + bnp - n;
255
0
      ASSERT (anp <= n || xp[2*n]==0);
256
0
      anp-= anp > n;
257
0
      cy = mpn_sub (xp, xp, n, xp + n, anp);
258
0
      xp[n] = 0;
259
0
      MPN_INCR_U (xp, n+1, cy);
260
0
    }
261
6.42M
  else
262
6.42M
    mpn_bc_mulmod_bnp1 (xp, ap1, bp1, n, xp);
263
6.42M
      }
264
265
      /* Here the CRT recomposition begins.
266
267
   xm <- (xp + xm)/2 = (xp + xm)B^n/2 mod (B^n-1)
268
   Division by 2 is a bitwise rotation.
269
270
   Assumes xp normalised mod (B^n+1).
271
272
   The residue class [0] is represented by [B^n-1]; except when
273
   both input are ZERO.
274
      */
275
276
6.42M
#if HAVE_NATIVE_mpn_rsh1add_n || HAVE_NATIVE_mpn_rsh1add_nc
277
6.42M
#if HAVE_NATIVE_mpn_rsh1add_nc
278
6.42M
      cy = mpn_rsh1add_nc(rp, rp, xp, n, xp[n]); /* B^n = 1 */
279
6.42M
      hi = cy << (GMP_NUMB_BITS - 1);
280
6.42M
      cy = 0;
281
      /* next update of rp[n-1] will set cy = 1 only if rp[n-1]+=hi
282
   overflows, i.e. a further increment will not overflow again. */
283
#else /* ! _nc */
284
      cy = xp[n] + mpn_rsh1add_n(rp, rp, xp, n); /* B^n = 1 */
285
      hi = (cy<<(GMP_NUMB_BITS-1))&GMP_NUMB_MASK; /* (cy&1) << ... */
286
      cy >>= 1;
287
      /* cy = 1 only if xp[n] = 1 i.e. {xp,n} = ZERO, this implies that
288
   the rsh1add was a simple rshift: the top bit is 0. cy=1 => hi=0. */
289
#endif
290
6.42M
#if GMP_NAIL_BITS == 0
291
6.42M
      add_ssaaaa(cy, rp[n-1], cy, rp[n-1], 0, hi);
292
#else
293
      cy += (hi & rp[n-1]) >> (GMP_NUMB_BITS-1);
294
      rp[n-1] ^= hi;
295
#endif
296
#else /* ! HAVE_NATIVE_mpn_rsh1add_n */
297
#if HAVE_NATIVE_mpn_add_nc
298
      cy = mpn_add_nc(rp, rp, xp, n, xp[n]);
299
#else /* ! _nc */
300
      cy = xp[n] + mpn_add_n(rp, rp, xp, n); /* xp[n] == 1 implies {xp,n} == ZERO */
301
#endif
302
      cy += (rp[0]&1);
303
      mpn_rshift(rp, rp, n, 1);
304
      ASSERT (cy <= 2);
305
      hi = (cy<<(GMP_NUMB_BITS-1))&GMP_NUMB_MASK; /* (cy&1) << ... */
306
      cy >>= 1;
307
      /* We can have cy != 0 only if hi = 0... */
308
      ASSERT ((rp[n-1] & GMP_NUMB_HIGHBIT) == 0);
309
      rp[n-1] |= hi;
310
      /* ... rp[n-1] + cy can not overflow, the following INCR is correct. */
311
#endif
312
6.42M
      ASSERT (cy <= 1);
313
      /* Next increment can not overflow, read the previous comments about cy. */
314
6.42M
      ASSERT ((cy == 0) || ((rp[n-1] & GMP_NUMB_HIGHBIT) == 0));
315
6.42M
      MPN_INCR_U(rp, n, cy);
316
317
      /* Compute the highest half:
318
   ([(xp + xm)/2 mod (B^n-1)] - xp ) * B^n
319
       */
320
6.42M
      if (UNLIKELY (an + bn < rn))
321
0
  {
322
    /* Note that in this case, the only way the result can equal
323
       zero mod B^{rn} - 1 is if one of the inputs is zero, and
324
       then the output of both the recursive calls and this CRT
325
       reconstruction is zero, not B^{rn} - 1. Which is good,
326
       since the latter representation doesn't fit in the output
327
       area.*/
328
0
    cy = mpn_sub_n (rp + n, rp, xp, an + bn - n);
329
330
    /* FIXME: This subtraction of the high parts is not really
331
       necessary, we do it to get the carry out, and for sanity
332
       checking. */
333
0
    cy = xp[n] + mpn_sub_nc (xp + an + bn - n, rp + an + bn - n,
334
0
           xp + an + bn - n, rn - (an + bn), cy);
335
0
    ASSERT (an + bn == rn - 1 ||
336
0
      mpn_zero_p (xp + an + bn - n + 1, rn - 1 - (an + bn)));
337
0
    cy = mpn_sub_1 (rp, rp, an + bn, cy);
338
0
    ASSERT (cy == (xp + an + bn - n)[0]);
339
0
  }
340
6.42M
      else
341
6.42M
  {
342
6.42M
    cy = xp[n] + mpn_sub_n (rp + n, rp, xp, n);
343
    /* cy = 1 only if {xp,n+1} is not ZERO, i.e. {rp,n} is not ZERO.
344
       DECR will affect _at most_ the lowest n limbs. */
345
6.42M
    MPN_DECR_U (rp, 2*n, cy);
346
6.42M
  }
347
6.42M
#undef a0
348
6.42M
#undef a1
349
6.42M
#undef b0
350
6.42M
#undef b1
351
6.42M
#undef xp
352
6.42M
#undef sp1
353
6.42M
    }
354
8.64M
}
355
356
mp_size_t
357
mpn_mulmod_bnm1_next_size (mp_size_t n)
358
2.25M
{
359
2.25M
  mp_size_t nh;
360
361
2.25M
  if (BELOW_THRESHOLD (n,     MULMOD_BNM1_THRESHOLD))
362
25.7k
    return n;
363
2.23M
  if (BELOW_THRESHOLD (n, 4 * (MULMOD_BNM1_THRESHOLD - 1) + 1))
364
1.22k
    return (n + (2-1)) & (-2);
365
2.23M
  if (BELOW_THRESHOLD (n, 8 * (MULMOD_BNM1_THRESHOLD - 1) + 1))
366
1.60M
    return (n + (4-1)) & (-4);
367
368
623k
  nh = (n + 1) >> 1;
369
370
623k
  if (BELOW_THRESHOLD (nh, MUL_FFT_MODF_THRESHOLD))
371
623k
    return (n + (8-1)) & (-8);
372
373
0
  return 2 * mpn_fft_next_size (nh, mpn_fft_best_k (nh, 0));
374
623k
}