Coverage Report

Created: 2025-07-23 06:43

/src/gmp/mpn/mulmod_bknp1.c
Line
Count
Source (jump to first uncovered line)
1
/* Mulptiplication mod B^n+1, for small operands.
2
3
   Contributed to the GNU project by Marco Bodrato.
4
5
   THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
6
   SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
7
   GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
8
9
Copyright 2020-2022 Free Software Foundation, Inc.
10
11
This file is part of the GNU MP Library.
12
13
The GNU MP Library is free software; you can redistribute it and/or modify
14
it under the terms of either:
15
16
  * the GNU Lesser General Public License as published by the Free
17
    Software Foundation; either version 3 of the License, or (at your
18
    option) any later version.
19
20
or
21
22
  * the GNU General Public License as published by the Free Software
23
    Foundation; either version 2 of the License, or (at your option) any
24
    later version.
25
26
or both in parallel, as here.
27
28
The GNU MP Library is distributed in the hope that it will be useful, but
29
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
30
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
31
for more details.
32
33
You should have received copies of the GNU General Public License and the
34
GNU Lesser General Public License along with the GNU MP Library.  If not,
35
see https://www.gnu.org/licenses/.  */
36
37
#include "gmp-impl.h"
38
39
#ifndef MOD_BKNP1_USE11
40
#define MOD_BKNP1_USE11 ((GMP_NUMB_BITS % 8 != 0) && (GMP_NUMB_BITS % 2 == 0))
41
#endif
42
#ifndef MOD_BKNP1_ONLY3
43
#define MOD_BKNP1_ONLY3 0
44
#endif
45
46
/* {rp, (k - 1) * n} = {op, k * n + 1} % (B^{k*n}+1) / (B^n+1) */
47
static void
48
_mpn_modbknp1dbnp1_n (mp_ptr rp, mp_srcptr op, mp_size_t n, unsigned k)
49
6.97M
{
50
6.97M
  mp_limb_t hl;
51
6.97M
  mp_srcptr hp;
52
6.97M
  unsigned i;
53
54
#if MOD_BKNP1_ONLY3
55
  ASSERT (k == 3);
56
  k = 3;
57
#endif
58
6.97M
  ASSERT (k > 2);
59
6.97M
  ASSERT (k % 2 == 1);
60
61
6.97M
  --k;
62
63
6.97M
  rp += k * n;
64
6.97M
  op += k * n;
65
6.97M
  hp = op;
66
6.97M
  hl = hp[n]; /* initial op[k*n]. */
67
6.97M
  ASSERT (hl < GMP_NUMB_MAX - 1);
68
69
6.97M
#if MOD_BKNP1_ONLY3 == 0
70
  /* The first MPN_INCR_U (rp + n, 1, cy); in the loop should be
71
     rp[n] = cy;            */
72
6.97M
  *rp = 0;
73
6.97M
#endif
74
75
6.97M
  i = k >> 1;
76
6.97M
  do
77
10.4M
   {
78
10.4M
     mp_limb_t cy, bw;
79
10.4M
     rp -= n;
80
10.4M
     op -= n;
81
10.4M
     cy = hl + mpn_add_n (rp, op, hp, n);
82
#if MOD_BKNP1_ONLY3
83
     rp[n] = cy;
84
#else
85
10.4M
     MPN_INCR_U (rp + n, (k - i * 2) * n + 1, cy);
86
10.4M
#endif
87
10.4M
     rp -= n;
88
10.4M
     op -= n;
89
10.4M
     bw = hl + mpn_sub_n (rp, op, hp, n);
90
10.4M
     MPN_DECR_U (rp + n, (k - i * 2 + 1) * n + 1, bw);
91
10.4M
   }
92
10.4M
  while (--i != 0);
93
94
10.3M
  for (; (hl = *(rp += k * n)) != 0; ) /* Should run only once... */
95
3.34M
    {
96
3.34M
      *rp = 0;
97
3.34M
      i = k >> 1;
98
3.34M
      do
99
5.10M
  {
100
5.10M
    rp -= n;
101
5.10M
    MPN_INCR_U (rp, (k - i * 2 + 1) * n + 1, hl);
102
5.10M
    rp -= n;
103
5.10M
    MPN_DECR_U (rp, (k - i * 2 + 2) * n + 1, hl);
104
5.10M
  }
105
5.10M
      while (--i != 0);
106
3.34M
    }
107
6.97M
}
108
109
static void
110
_mpn_modbnp1_pn_ip (mp_ptr r, mp_size_t n, mp_limb_t h)
111
5.55M
{
112
5.55M
  ASSERT (r[n] == h);
113
114
  /* Fully normalise */
115
5.55M
  MPN_DECR_U (r, n + 1, h);
116
5.55M
  h -= r[n];
117
5.55M
  r[n] = 0;
118
5.55M
  MPN_INCR_U (r, n + 1, h);
119
5.55M
}
120
121
static void
122
_mpn_modbnp1_neg_ip (mp_ptr r, mp_size_t n, mp_limb_t h)
123
3.64M
{
124
3.64M
  r[n] = 0;
125
3.64M
  MPN_INCR_U (r, n + 1, -h);
126
3.64M
  if (UNLIKELY (r[n] != 0))
127
4.02k
    _mpn_modbnp1_pn_ip (r, n, 1);
128
3.64M
}
129
130
static void
131
_mpn_modbnp1_nc_ip (mp_ptr r, mp_size_t n, mp_limb_t h)
132
10.4M
{
133
10.4M
  if (h & GMP_NUMB_HIGHBIT) /* This means h < 0 */
134
1.96M
    {
135
1.96M
      _mpn_modbnp1_neg_ip (r, n, h);
136
1.96M
    }
137
8.49M
  else
138
8.49M
    {
139
8.49M
      r[n] = h;
140
8.49M
      if (h)
141
2.05M
  _mpn_modbnp1_pn_ip(r, n, h);
142
8.49M
    }
143
10.4M
}
144
145
/* {rp, rn + 1} = {op, on} mod (B^{rn}+1) */
146
/* Used when rn < on < 2*rn. */
147
static void
148
_mpn_modbnp1 (mp_ptr rp, mp_size_t rn, mp_srcptr op, mp_size_t on)
149
3.48M
{
150
3.48M
  mp_limb_t bw;
151
152
#if 0
153
  if (UNLIKELY (on <= rn))
154
    {
155
      MPN_COPY (rp, op, on);
156
      MPN_ZERO (rp + on, rn - on);
157
      return;
158
    }
159
#endif
160
161
3.48M
  ASSERT (on > rn);
162
3.48M
  ASSERT (on <= 2 * rn);
163
164
3.48M
  bw = mpn_sub (rp, op, rn, op + rn, on - rn);
165
3.48M
  rp[rn] = 0;
166
3.48M
  MPN_INCR_U (rp, rn + 1, bw);
167
3.48M
}
168
169
/* {rp, rn + 1} = {op, k * rn + 1} % (B^{rn}+1) */
170
/* With odd k >= 3. */
171
static void
172
_mpn_modbnp1_kn (mp_ptr rp, mp_srcptr op, mp_size_t rn, unsigned k)
173
10.4M
{
174
10.4M
  mp_limb_t cy;
175
176
#if MOD_BKNP1_ONLY3
177
  ASSERT (k == 3);
178
  k = 3;
179
#endif
180
10.4M
  ASSERT (k & 1);
181
10.4M
  k >>= 1;
182
10.4M
  ASSERT (0 < k && k < GMP_NUMB_HIGHBIT - 3);
183
10.4M
  ASSERT (op[(1 + 2 * k) * rn] < GMP_NUMB_HIGHBIT - 2 - k);
184
185
10.4M
  cy = - mpn_sub_n (rp, op, op + rn, rn);
186
15.6M
  for (;;) {
187
15.6M
    op += 2 * rn;
188
15.6M
    cy += mpn_add_n (rp, rp, op, rn);
189
15.6M
    if (--k == 0)
190
10.4M
      break;
191
5.14M
    cy -= mpn_sub_n (rp, rp, op + rn, rn);
192
5.14M
  };
193
194
10.4M
  cy += op[rn];
195
10.4M
  _mpn_modbnp1_nc_ip (rp, rn, cy);
196
10.4M
}
197
198
/* For the various mpn_divexact_byN here, fall back to using either
199
   mpn_pi1_bdiv_q_1 or mpn_divexact_1.  The former has less overhead and is
200
   faster if it is native.  For now, since mpn_divexact_1 is native on
201
   platforms where mpn_pi1_bdiv_q_1 does not yet exist, do not use
202
   mpn_pi1_bdiv_q_1 unconditionally.  FIXME.  */
203
204
#ifndef mpn_divexact_by5
205
#if HAVE_NATIVE_mpn_pi1_bdiv_q_1
206
#define BINVERT_5 \
207
  ((((GMP_NUMB_MAX >> (GMP_NUMB_BITS % 4)) / 5 * 3 << 3) + 5) & GMP_NUMB_MAX)
208
#define mpn_divexact_by5(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,5,BINVERT_5,0)
209
#else
210
#define mpn_divexact_by5(dst,src,size) mpn_divexact_1(dst,src,size,5)
211
#endif
212
#endif
213
214
#ifndef mpn_divexact_by7
215
#if HAVE_NATIVE_mpn_pi1_bdiv_q_1
216
#define BINVERT_7 \
217
  ((((GMP_NUMB_MAX >> (GMP_NUMB_BITS % 3)) / 7 * 3 << 4) + 7) & GMP_NUMB_MAX)
218
#define mpn_divexact_by7(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,7,BINVERT_7,0)
219
#else
220
435k
#define mpn_divexact_by7(dst,src,size) mpn_divexact_1(dst,src,size,7)
221
#endif
222
#endif
223
224
#ifndef mpn_divexact_by11
225
#if HAVE_NATIVE_mpn_pi1_bdiv_q_1
226
#define BINVERT_11 \
227
  ((((GMP_NUMB_MAX >> (GMP_NUMB_BITS % 10)) / 11 << 5) + 3) & GMP_NUMB_MAX)
228
#define mpn_divexact_by11(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,11,BINVERT_11,0)
229
#else
230
#define mpn_divexact_by11(dst,src,size) mpn_divexact_1(dst,src,size,11)
231
#endif
232
#endif
233
234
#ifndef mpn_divexact_by13
235
#if HAVE_NATIVE_mpn_pi1_bdiv_q_1
236
#define BINVERT_13 \
237
  ((((GMP_NUMB_MAX >> (GMP_NUMB_BITS % 12)) / 13 * 3 << 14) + 3781) & GMP_NUMB_MAX)
238
#define mpn_divexact_by13(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,13,BINVERT_13,0)
239
#else
240
52.1k
#define mpn_divexact_by13(dst,src,size) mpn_divexact_1(dst,src,size,13)
241
#endif
242
#endif
243
244
#ifndef mpn_divexact_by17
245
#if HAVE_NATIVE_mpn_pi1_bdiv_q_1
246
#define BINVERT_17 \
247
  ((((GMP_NUMB_MAX >> (GMP_NUMB_BITS % 8)) / 17 * 15 << 7) + 113) & GMP_NUMB_MAX)
248
#define mpn_divexact_by17(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,17,BINVERT_17,0)
249
#else
250
#define mpn_divexact_by17(dst,src,size) mpn_divexact_1(dst,src,size,17)
251
#endif
252
#endif
253
254
/* Thanks to Chinese remainder theorem, store
255
   in {rp, k*n+1} the value mod (B^(k*n)+1), given
256
   {ap, k*n+1} mod ((B^(k*n)+1)/(B^n+1)) and
257
   {bp, n+1} mod (B^n+1) .
258
   {tp, n+1} is a scratch area.
259
   tp == rp or rp == ap are possible.
260
*/
261
static void
262
_mpn_crt (mp_ptr rp, mp_srcptr ap, mp_srcptr bp,
263
    mp_size_t n, unsigned k, mp_ptr tp)
264
3.48M
{
265
3.48M
  mp_limb_t mod;
266
3.48M
  unsigned i;
267
268
#if MOD_BKNP1_ONLY3
269
  ASSERT (k == 3);
270
  k = 3;
271
#endif
272
3.48M
  _mpn_modbnp1_kn (tp, ap, n, k);
273
3.48M
  if (mpn_sub_n (tp, bp, tp, n + 1))
274
1.67M
    _mpn_modbnp1_neg_ip (tp, n, tp[n]);
275
276
#if MOD_BKNP1_USE11
277
  if (UNLIKELY (k == 11))
278
    {
279
      ASSERT (GMP_NUMB_BITS % 2 == 0);
280
      /* mod <- -Mod(B^n+1,11)^-1 */
281
      mod = n * (GMP_NUMB_BITS % 5) % 5;
282
      if ((mod > 2) || UNLIKELY (mod == 0))
283
  mod += 5;
284
285
      mod *= mpn_mod_1 (tp, n + 1, 11);
286
    }
287
  else
288
#endif
289
3.48M
    {
290
3.48M
#if GMP_NUMB_BITS % 8 == 0
291
  /* (2^6 - 1) | (2^{GMP_NUMB_BITS*3/4} - 1)  */
292
  /* (2^6 - 1) = 3^2 * 7      */
293
3.48M
  mod = mpn_mod_34lsub1 (tp, n + 1);
294
3.48M
  ASSERT ((GMP_NUMB_MAX >> (GMP_NUMB_BITS >> 2)) % k == 0);
295
  /* (2^12 - 1) = 3^2 * 5 * 7 * 13    */
296
  /* (2^24 - 1) = 3^2 * 5 * 7 * 13 * 17 * 241 */
297
3.48M
  ASSERT (k == 3 || k == 5 || k == 7 || k == 13 || k == 17);
298
299
3.48M
#if GMP_NUMB_BITS % 3 != 0
300
3.48M
  if (UNLIKELY (k != 3))
301
1.06M
    {
302
1.06M
      ASSERT ((GMP_NUMB_MAX % k == 0) || (n % 3 != 0));
303
1.06M
      if ((GMP_NUMB_BITS % 16 == 0) && LIKELY (k == 5))
304
582k
  mod <<= 1; /* k >> 1 = 1 << 1 */
305
487k
      else if ((GMP_NUMB_BITS % 16 != 0) || LIKELY (k == 7))
306
435k
  mod <<= (n << (GMP_NUMB_BITS % 3 >> 1)) % 3;
307
52.1k
      else if ((GMP_NUMB_BITS % 32 != 0) || LIKELY (k == 13))
308
52.1k
  mod *= ((n << (GMP_NUMB_BITS % 3 >> 1)) % 3 == 1) ? 3 : 9;
309
0
      else /* k == 17 */
310
0
  mod <<= 3; /* k >> 1 = 1 << 3 */
311
#if 0
312
      if ((GMP_NUMB_BITS == 8) /* && (k == 7) */ ||
313
    (GMP_NUMB_BITS == 16) && (k == 13))
314
  mod = ((mod & (GMP_NUMB_MAX >> (GMP_NUMB_BITS >> 2))) +
315
         (mod >> (3 * GMP_NUMB_BITS >> 2)));
316
#endif
317
1.06M
    }
318
#else
319
  ASSERT (GMP_NUMB_MAX % k == 0);
320
  /* 2^{GMP_NUMB_BITS} - 1  = 0 (mod k) */
321
  /* 2^{GMP_NUMB_BITS}    = 1 (mod k) */
322
  /* 2^{n*GMP_NUMB_BITS} + 1  = 2 (mod k) */
323
  /* -2^{-1}  = k >> 1 (mod k) */
324
  mod *= k >> 1;
325
#endif
326
#else
327
  ASSERT_ALWAYS (k == 0); /* Not implemented, should not be used. */
328
#endif
329
3.48M
    }
330
331
3.48M
  MPN_INCR_U (tp, n + 1, mod);
332
3.48M
  tp[n] += mod;
333
334
3.48M
  if (LIKELY (k == 3))
335
2.41M
    ASSERT_NOCARRY (mpn_divexact_by3 (tp, tp, n + 1));
336
1.06M
  else if ((GMP_NUMB_BITS % 16 == 0) && LIKELY (k == 5))
337
582k
    mpn_divexact_by5 (tp, tp, n + 1);
338
487k
  else if (((! MOD_BKNP1_USE11) && (GMP_NUMB_BITS % 16 != 0))
339
487k
     || LIKELY (k == 7))
340
487k
    mpn_divexact_by7 (tp, tp, n + 1);
341
#if MOD_BKNP1_USE11
342
  else if (k == 11)
343
    mpn_divexact_by11 (tp, tp, n + 1);
344
#endif
345
52.1k
  else if ((GMP_NUMB_BITS % 32 != 0) || LIKELY (k == 13))
346
52.1k
    mpn_divexact_by13 (tp, tp, n + 1);
347
0
  else /* (k == 17) */
348
0
    mpn_divexact_by17 (tp, tp, n + 1);
349
350
3.48M
  rp += k * n;
351
3.48M
  ap += k * n; /* tp - 1 */
352
353
3.48M
  rp -= n;
354
3.48M
  ap -= n;
355
3.48M
  ASSERT_NOCARRY (mpn_add_n (rp, ap, tp, n + 1));
356
357
3.48M
  i = k >> 1;
358
3.48M
  do
359
5.20M
   {
360
5.20M
      mp_limb_t cy, bw;
361
5.20M
      rp -= n;
362
5.20M
      ap -= n;
363
5.20M
      bw = mpn_sub_n (rp, ap, tp, n) + tp[n];
364
5.20M
      MPN_DECR_U (rp + n, (k - i * 2) * n + 1, bw);
365
5.20M
      rp -= n;
366
5.20M
      ap -= n;
367
5.20M
      cy = mpn_add_n (rp, ap, tp, n) + tp[n];
368
5.20M
      MPN_INCR_U (rp + n, (k - i * 2 + 1) * n + 1, cy);
369
5.20M
    }
370
5.20M
  while (--i != 0);
371
372
  /* if (LIKELY (rp[k * n])) */
373
3.48M
    _mpn_modbnp1_pn_ip (rp, k * n, rp[k * n]);
374
3.48M
}
375
376
377
static void
378
_mpn_mulmod_bnp1_tp (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t rn,
379
        mp_ptr tp)
380
3.48M
{
381
3.48M
  mp_limb_t cy;
382
3.48M
  unsigned k;
383
384
3.48M
  ASSERT (0 < rn);
385
3.48M
  ASSERT ((ap[rn] | bp[rn]) <= 1);
386
387
3.48M
  if (UNLIKELY (ap[rn] | bp[rn]))
388
12.8k
    {
389
12.8k
      if (ap[rn])
390
329
  cy = bp[rn] + mpn_neg (rp, bp, rn);
391
12.4k
      else /* ap[rn] == 0 */
392
12.4k
  cy = mpn_neg (rp, ap, rn);
393
12.8k
    }
394
3.47M
  else if (MPN_MULMOD_BKNP1_USABLE(rn, k, MUL_FFT_MODF_THRESHOLD / 3))
395
639k
    {
396
639k
      rn /= k;
397
639k
      mpn_mulmod_bknp1 (rp, ap, bp, rn, k, tp);
398
639k
      return;
399
639k
    }
400
2.83M
  else
401
2.83M
    {
402
2.83M
      mpn_mul_n (tp, ap, bp, rn);
403
2.83M
      cy = mpn_sub_n (rp, tp, tp + rn, rn);
404
2.83M
    }
405
2.84M
  rp[rn] = 0;
406
2.84M
  MPN_INCR_U (rp, rn + 1, cy);
407
2.84M
}
408
409
/* {rp, kn + 1} = {ap, kn + 1} * {bp, kn + 1} % (B^kn + 1) */
410
/* tp must point to at least 4*(k-1)*n+1 limbs*/
411
void
412
mpn_mulmod_bknp1 (mp_ptr rp, mp_srcptr ap, mp_srcptr bp,
413
      mp_size_t n, unsigned k, mp_ptr tp)
414
3.48M
{
415
3.48M
  mp_ptr hp;
416
417
#if MOD_BKNP1_ONLY3
418
  ASSERT (k == 3);
419
  k = 3;
420
#endif
421
3.48M
  ASSERT (k > 2);
422
3.48M
  ASSERT (k % 2 == 1);
423
424
  /* a % (B^{nn}+1)/(B^{nn/k}+1) */
425
3.48M
  _mpn_modbknp1dbnp1_n (tp + (k - 1) * n * 2, ap, n, k);
426
  /* b % (B^{nn}+1)/(B^{nn/k}+1) */
427
3.48M
  _mpn_modbknp1dbnp1_n (tp + (k - 1) * n * 3, bp, n, k);
428
3.48M
  mpn_mul_n (tp, tp + (k - 1) * n * 2, tp + (k - 1) * n * 3, (k - 1) * n);
429
3.48M
  _mpn_modbnp1 (tp, k * n, tp, (k - 1) * n * 2);
430
431
3.48M
  hp = tp + k * n + 1;
432
  /* a % (B^{nn/k}+1) */
433
3.48M
  ASSERT (ap[k * n] <= 1);
434
3.48M
  _mpn_modbnp1_kn (hp, ap, n, k);
435
  /* b % (B^{nn/k}+1) */
436
3.48M
  ASSERT (bp[k * n] <= 1);
437
3.48M
  _mpn_modbnp1_kn (hp + n + 1, bp, n, k);
438
3.48M
  _mpn_mulmod_bnp1_tp (hp + (n + 1) * 2, hp, hp + n + 1, n, hp + (n + 1) * 2);
439
440
3.48M
  _mpn_crt (rp, tp, hp + (n + 1) * 2, n, k, hp);
441
3.48M
}
442
443
444
static void
445
_mpn_sqrmod_bnp1_tp (mp_ptr rp, mp_srcptr ap, mp_size_t rn,
446
        mp_ptr tp)
447
0
{
448
0
  mp_limb_t cy;
449
0
  unsigned k;
450
451
0
  ASSERT (0 < rn);
452
453
0
  if (UNLIKELY (ap[rn]))
454
0
    {
455
0
      ASSERT (ap[rn] == 1);
456
0
      *rp = 1;
457
0
      MPN_FILL (rp + 1, rn, 0);
458
0
      return;
459
0
    }
460
0
  else if (MPN_SQRMOD_BKNP1_USABLE(rn, k, MUL_FFT_MODF_THRESHOLD / 3))
461
0
    {
462
0
      rn /= k;
463
0
      mpn_sqrmod_bknp1 (rp, ap, rn, k, tp);
464
0
      return;
465
0
    }
466
0
  else
467
0
    {
468
0
      mpn_sqr (tp, ap, rn);
469
0
      cy = mpn_sub_n (rp, tp, tp + rn, rn);
470
0
    }
471
0
  rp[rn] = 0;
472
0
  MPN_INCR_U (rp, rn + 1, cy);
473
0
}
474
475
/* {rp, kn + 1} = {ap, kn + 1}^2 % (B^kn + 1) */
476
/* tp must point to at least 3*(k-1)*n+1 limbs*/
477
void
478
mpn_sqrmod_bknp1 (mp_ptr rp, mp_srcptr ap,
479
      mp_size_t n, unsigned k, mp_ptr tp)
480
0
{
481
0
  mp_ptr hp;
482
483
#if MOD_BKNP1_ONLY3
484
  ASSERT (k == 3);
485
  k = 3;
486
#endif
487
0
  ASSERT (k > 2);
488
0
  ASSERT (k % 2 == 1);
489
490
  /* a % (B^{nn}+1)/(B^{nn/k}+1) */
491
0
  _mpn_modbknp1dbnp1_n (tp + (k - 1) * n * 2, ap, n, k);
492
0
  mpn_sqr (tp, tp + (k - 1) * n * 2, (k - 1) * n);
493
0
  _mpn_modbnp1 (tp, k * n, tp, (k - 1) * n * 2);
494
495
0
  hp = tp + k * n + 1;
496
  /* a % (B^{nn/k}+1) */
497
0
  ASSERT (ap[k * n] <= 1);
498
0
  _mpn_modbnp1_kn (hp, ap, n, k);
499
0
  _mpn_sqrmod_bnp1_tp (hp + (n + 1), hp, n, hp + (n + 1));
500
501
0
  _mpn_crt (rp, tp, hp + (n + 1), n, k, hp);
502
0
}