Make function intToSBigNat# preserve sign (fixes #14085)
[ghc.git] / libraries / integer-gmp / cbits / wrappers.c
1 /*
2 * `integer-gmp` GMP FFI wrappers
3 *
4 * Copyright (c) 2014, Herbert Valerio Riedel <hvr@gnu.org>
5 *
6 * BSD3 licensed, see ../LICENSE file for details
7 *
8 */
9
10 #define _ISOC99_SOURCE
11
12 #include "HsFFI.h"
13 #include "MachDeps.h"
14
15 #include <assert.h>
16 #include <stdbool.h>
17 #include <stdlib.h>
18 #include <stdint.h>
19 #include <string.h>
20 #include <math.h>
21 #include <float.h>
22 #include <stdio.h>
23
24 #include <gmp.h>
25
26
27 // GMP 4.x compatibility
28 #if !defined(__GNU_MP_VERSION)
29 # error __GNU_MP_VERSION not defined
30 #elif __GNU_MP_VERSION < 4
31 # error need GMP 4.0 or later
32 #elif __GNU_MP_VERSION < 5
33 typedef unsigned long int mp_bitcnt_t;
34 #endif
35
36 #if (GMP_NUMB_BITS) != (GMP_LIMB_BITS)
37 # error GMP_NUMB_BITS != GMP_LIMB_BITS not supported
38 #endif
39
40 #if (WORD_SIZE_IN_BITS) != (GMP_LIMB_BITS)
41 # error WORD_SIZE_IN_BITS != GMP_LIMB_BITS not supported
42 #endif
43
44 // sanity check
45 #if (SIZEOF_HSWORD*8) != WORD_SIZE_IN_BITS
46 # error (SIZEOF_HSWORD*8) != WORD_SIZE_IN_BITS
47 #endif
48
49 // Turn a (const) {xp,xn} pair into static initializer
50 #define CONST_MPZ_INIT(xp,xn) \
51 {{ ._mp_alloc = 0, ._mp_size = (xn), ._mp_d = (mp_limb_t*)(xp) }}
52
53 // Test if {sp,sn} represents a zero value
54 static inline int
55 mp_limb_zero_p(const mp_limb_t sp[], mp_size_t sn)
56 {
57 return !sn || ((sn == 1 || sn == -1) && !sp[0]);
58 }
59
60 static inline mp_size_t
61 mp_size_abs(const mp_size_t x)
62 {
63 return x>=0 ? x : -x;
64 }
65
66 static inline mp_size_t
67 mp_size_min(const mp_size_t x, const mp_size_t y)
68 {
69 return x<y ? x : y;
70 }
71
72 static inline mp_size_t
73 mp_size_minabs(const mp_size_t x, const mp_size_t y)
74 {
75 return mp_size_min(mp_size_abs(x), mp_size_abs(y));
76 }
77
78 /* Perform arithmetic right shift on MPNs (multi-precision naturals)
79 *
80 * pre-conditions:
81 * - 0 < count < sn*GMP_NUMB_BITS
82 * - rn = sn - floor(count / GMP_NUMB_BITS)
83 * - sn > 0
84 *
85 * write {sp,sn} right-shifted by count bits into {rp,rn}
86 *
87 * return value: most-significant limb stored in {rp,rn} result
88 */
89 mp_limb_t
90 integer_gmp_mpn_rshift (mp_limb_t rp[], const mp_limb_t sp[], mp_size_t sn,
91 mp_bitcnt_t count)
92 {
93 const mp_size_t limb_shift = count / GMP_NUMB_BITS;
94 const unsigned int bit_shift = count % GMP_NUMB_BITS;
95 const mp_size_t rn = sn - limb_shift;
96
97 if (bit_shift)
98 mpn_rshift(rp, &sp[limb_shift], rn, bit_shift);
99 else
100 memcpy(rp, &sp[limb_shift], rn*sizeof(mp_limb_t));
101
102 return rp[rn-1];
103 }
104
105 /* Twos-complement version of 'integer_gmp_mpn_rshift' for performing
106 * arithmetic right shifts on "negative" MPNs.
107 *
108 * pre-conditions:
109 * - 0 < count < sn*GMP_NUMB_BITS
110 * - rn = sn - floor((count - 1) / GMP_NUMB_BITS)
111 * - sn > 0
112 *
113 * This variant is needed to operate on MPNs interpreted as negative
114 * numbers, which require "rounding" towards minus infinity iff a
115 * non-zero bit is shifted out.
116 */
117 mp_limb_t
118 integer_gmp_mpn_rshift_2c (mp_limb_t rp[], const mp_limb_t sp[],
119 const mp_size_t sn, const mp_bitcnt_t count)
120 {
121 const mp_size_t limb_shift = count / GMP_NUMB_BITS;
122 const unsigned int bit_shift = count % GMP_NUMB_BITS;
123 mp_size_t rn = sn - limb_shift;
124
125 // whether non-zero bits were shifted out
126 bool nz_shift_out = false;
127
128 if (bit_shift) {
129 if (mpn_rshift(rp, &sp[limb_shift], rn, bit_shift))
130 nz_shift_out = true;
131 } else {
132 // rp was allocated (rn + 1) limbs, to prevent carry
133 // on mpn_add_1 when all the bits of {rp, rn} are 1.
134 memset(&rp[rn], 0, sizeof(mp_limb_t));
135 memcpy(rp, &sp[limb_shift], rn*sizeof(mp_limb_t));
136 rn++;
137 }
138
139 if (!nz_shift_out)
140 for (unsigned i = 0; i < limb_shift; i++)
141 if (sp[i]) {
142 nz_shift_out = true;
143 break;
144 }
145
146 // round if non-zero bits were shifted out
147 if (nz_shift_out)
148 if (mpn_add_1(rp, rp, rn, 1))
149 abort(); /* should never happen */
150
151 return rp[rn-1];
152 }
153
154 /* Perform left-shift operation on MPN
155 *
156 * pre-conditions:
157 * - 0 < count
158 * - rn = sn + ceil(count / GMP_NUMB_BITS)
159 * - sn > 0
160 *
161 * return value: most-significant limb stored in {rp,rn} result
162 */
163 mp_limb_t
164 integer_gmp_mpn_lshift (mp_limb_t rp[], const mp_limb_t sp[],
165 const mp_size_t sn, const mp_bitcnt_t count)
166 {
167 const mp_size_t limb_shift = count / GMP_NUMB_BITS;
168 const unsigned int bit_shift = count % GMP_NUMB_BITS;
169 const mp_size_t rn0 = sn + limb_shift;
170
171 memset(rp, 0, limb_shift*sizeof(mp_limb_t));
172 if (bit_shift) {
173 const mp_limb_t msl = mpn_lshift(&rp[limb_shift], sp, sn, bit_shift);
174 rp[rn0] = msl;
175 return msl;
176 } else {
177 memcpy(&rp[limb_shift], sp, sn*sizeof(mp_limb_t));
178 return rp[rn0-1];
179 }
180 }
181
182 /* Convert bignum to a `double`, truncating if necessary
183 * (i.e. rounding towards zero).
184 *
185 * sign of mp_size_t argument controls sign of converted double
186 */
187 HsDouble
188 integer_gmp_mpn_get_d (const mp_limb_t sp[], const mp_size_t sn,
189 const HsInt exponent)
190 {
191 if (mp_limb_zero_p(sp, sn))
192 return 0.0;
193
194 const mpz_t mpz = CONST_MPZ_INIT(sp, sn);
195
196 if (!exponent)
197 return mpz_get_d(mpz);
198
199 long e = 0;
200 double d = mpz_get_d_2exp (&e, mpz);
201
202 // TODO: over/underflow handling?
203 return ldexp(d, e+exponent);
204 }
205
206 mp_limb_t
207 integer_gmp_gcd_word(const mp_limb_t x, const mp_limb_t y)
208 {
209 if (!x) return y;
210 if (!y) return x;
211
212 return mpn_gcd_1(&x, 1, y);
213 }
214
215 mp_limb_t
216 integer_gmp_mpn_gcd_1(const mp_limb_t x[], const mp_size_t xn,
217 const mp_limb_t y)
218 {
219 assert (xn > 0);
220 assert (xn == 1 || y != 0);
221
222 if (xn == 1)
223 return integer_gmp_gcd_word(x[0], y);
224
225 return mpn_gcd_1(x, xn, y);
226 }
227
228
229 mp_size_t
230 integer_gmp_mpn_gcd(mp_limb_t r[],
231 const mp_limb_t x0[], const mp_size_t xn,
232 const mp_limb_t y0[], const mp_size_t yn)
233 {
234 assert (xn >= yn);
235 assert (yn > 0);
236 assert (xn == yn || yn > 1 || y0[0] != 0);
237 /* post-condition: rn <= xn */
238
239 if (yn == 1) {
240 if (y0[0]) {
241 r[0] = integer_gmp_mpn_gcd_1(x0, xn, y0[0]);
242 return 1;
243 } else { /* {y0,yn} == 0 */
244 assert (xn==yn); /* NB: redundant assertion */
245 memcpy(r, x0, xn*sizeof(mp_limb_t));
246 return xn;
247 }
248 } else {
249 // mpn_gcd() seems to require non-trivial normalization of its
250 // input arguments (which does not seem to be documented anywhere,
251 // see source of mpz_gcd() for more details), so we resort to just
252 // use mpz_gcd() which does the tiresome normalization for us at
253 // the cost of a few additional temporary buffer allocations in
254 // C-land.
255
256 const mpz_t op1 = CONST_MPZ_INIT(x0, xn);
257 const mpz_t op2 = CONST_MPZ_INIT(y0, yn);
258
259 mpz_t rop;
260 mpz_init (rop);
261
262 mpz_gcd(rop, op1, op2);
263
264 const mp_size_t rn = rop[0]._mp_size;
265 assert(rn > 0);
266 assert(rn <= xn);
267
268 /* the allocation/memcpy of the result can be neglectable since
269 mpz_gcd() already has to allocate other temporary buffers
270 anyway */
271 memcpy(r, rop[0]._mp_d, rn*sizeof(mp_limb_t));
272
273 mpz_clear(rop);
274
275 return rn;
276 }
277 }
278
279 /* wraps mpz_gcdext()
280 *
281 * Set g to the greatest common divisor of x and y, and in addition
282 * set s and t to coefficients satisfying x*s + y*t = g.
283 *
284 * The {gp,gn} array is zero-padded (as otherwise 'gn' can't be
285 * reconstructed).
286 *
287 * g must have space for exactly gn=min(xn,yn) limbs.
288 * s must have space for at least xn limbs.
289 *
290 * return value: signed 'sn' of {sp,sn}
291 */
292 mp_size_t
293 integer_gmp_gcdext(mp_limb_t s0[], mp_limb_t g0[],
294 const mp_limb_t x0[], const mp_size_t xn,
295 const mp_limb_t y0[], const mp_size_t yn)
296 {
297 const mp_size_t gn0 = mp_size_minabs(xn, yn);
298 const mpz_t x = CONST_MPZ_INIT(x0, mp_limb_zero_p(x0,xn) ? 0 : xn);
299 const mpz_t y = CONST_MPZ_INIT(y0, mp_limb_zero_p(y0,yn) ? 0 : yn);
300
301 mpz_t g, s;
302 mpz_init (g);
303 mpz_init (s);
304
305 mpz_gcdext (g, s, NULL, x, y);
306
307 const mp_size_t gn = g[0]._mp_size;
308 assert(0 <= gn && gn <= gn0);
309 memset(g0, 0, gn0*sizeof(mp_limb_t));
310 memcpy(g0, g[0]._mp_d, gn*sizeof(mp_limb_t));
311 mpz_clear (g);
312
313 const mp_size_t ssn = s[0]._mp_size;
314 const mp_size_t sn = mp_size_abs(ssn);
315 assert(sn <= mp_size_abs(xn));
316 memcpy(s0, s[0]._mp_d, sn*sizeof(mp_limb_t));
317 mpz_clear (s);
318
319 if (!sn) {
320 s0[0] = 0;
321 return 1;
322 }
323
324 return ssn;
325 }
326
327 /* Truncating (i.e. rounded towards zero) integer division-quotient of MPN */
328 void
329 integer_gmp_mpn_tdiv_q (mp_limb_t q[],
330 const mp_limb_t n[], const mp_size_t nn,
331 const mp_limb_t d[], const mp_size_t dn)
332 {
333 /* qn = 1+nn-dn; rn = dn */
334 assert(nn>=dn);
335
336 if (dn > 128) {
337 // Use temporary heap allocated throw-away buffer for MPNs larger
338 // than 1KiB for 64bit-sized limbs (larger than 512bytes for
339 // 32bit-sized limbs)
340 mp_limb_t *const r = malloc(dn*sizeof(mp_limb_t));
341 mpn_tdiv_qr(q, r, 0, n, nn, d, dn);
342 free (r);
343 } else { // allocate smaller arrays on the stack
344 mp_limb_t r[dn];
345 mpn_tdiv_qr(q, r, 0, n, nn, d, dn);
346 }
347 }
348
349 /* Truncating (i.e. rounded towards zero) integer division-remainder of MPNs */
350 void
351 integer_gmp_mpn_tdiv_r (mp_limb_t r[],
352 const mp_limb_t n[], const mp_size_t nn,
353 const mp_limb_t d[], const mp_size_t dn)
354 {
355 /* qn = 1+nn-dn; rn = dn */
356 assert(nn>=dn);
357 const mp_size_t qn = 1+nn-dn;
358
359 if (qn > 128) {
360 // Use temporary heap allocated throw-away buffer for MPNs larger
361 // than 1KiB for 64bit-sized limbs (larger than 512bytes for
362 // 32bit-sized limbs)
363 mp_limb_t *const q = malloc(qn*sizeof(mp_limb_t));
364 mpn_tdiv_qr(q, r, 0, n, nn, d, dn);
365 free (q);
366 } else { // allocate smaller arrays on the stack
367 mp_limb_t q[qn];
368 mpn_tdiv_qr(q, r, 0, n, nn, d, dn);
369 }
370 }
371
372
373 /* Wraps GMP's 'mpz_sizeinbase()' function */
374 HsWord
375 integer_gmp_mpn_sizeinbase(const mp_limb_t s[], const mp_size_t sn,
376 const HsInt base)
377 {
378 assert (2 <= base && base <= 256);
379
380 if (mp_limb_zero_p(s,sn)) return 1;
381
382 const mpz_t zs = CONST_MPZ_INIT(s, sn);
383
384 return mpz_sizeinbase(zs, base);
385 }
386
387 /* Single-limb version of 'integer_gmp_mpn_sizeinbase()' */
388 HsWord
389 integer_gmp_mpn_sizeinbase1(const mp_limb_t s, const HsInt base)
390 {
391 return s ? integer_gmp_mpn_sizeinbase(&s, 1, base) : 1;
392 }
393
394 /* Wrapper around GMP's 'mpz_export()' function */
395 HsWord
396 integer_gmp_mpn_export(const mp_limb_t s[], const mp_size_t sn,
397 void *destptr, HsInt destofs, HsInt msbf)
398 {
399 /* TODO: implement w/o GMP, c.f. 'integer_gmp_mpn_import()' */
400 assert (msbf == 0 || msbf == 1);
401
402 if (mp_limb_zero_p(s,sn)) return 0;
403
404 const mpz_t zs = CONST_MPZ_INIT(s, sn);
405
406 size_t written = 0;
407
408 // mpz_export (void *rop, size_t *countp, int order, size_t size, int endian,
409 // size_t nails, const mpz_t op)
410 (void) mpz_export(((char *)destptr)+destofs, &written, !msbf ? -1 : 1,
411 /* size */ 1, /* endian */ 0, /* nails */ 0, zs);
412
413 return written;
414 }
415
416 /* Single-limb version of 'integer_gmp_mpn_export()' */
417 HsWord
418 integer_gmp_mpn_export1(const mp_limb_t s,
419 void *destptr, const HsInt destofs, const HsInt msbf)
420 {
421 /* TODO: implement w/o GMP */
422 return integer_gmp_mpn_export(&s, 1, destptr, destofs, msbf);
423 }
424
425 /* Import single limb from memory location
426 *
427 * We can't use GMP's 'mpz_import()'
428 */
429 HsWord
430 integer_gmp_mpn_import1(const uint8_t *srcptr, const HsWord srcofs,
431 const HsWord srclen, const HsInt msbf)
432 {
433 assert (msbf == 0 || msbf == 1);
434 assert (srclen <= SIZEOF_HSWORD);
435
436 srcptr += srcofs;
437
438 HsWord result = 0;
439
440 if (msbf)
441 for (unsigned i = 0; i < srclen; ++i)
442 result |= (HsWord)srcptr[i] << ((srclen-i-1)*8);
443 else // lsbf
444 for (unsigned i = 0; i < srclen; ++i)
445 result |= (HsWord)srcptr[i] << (i*8);
446
447 return result;
448 }
449
450 /* import into mp_limb_t[] from memory location */
451 void
452 integer_gmp_mpn_import(mp_limb_t * restrict r, const uint8_t * restrict srcptr,
453 const HsWord srcofs, const HsWord srclen,
454 const HsInt msbf)
455 {
456 assert (msbf == 0 || msbf == 1);
457
458 srcptr += srcofs;
459
460 const unsigned limb_cnt_rem = srclen % SIZEOF_HSWORD;
461 const mp_size_t limb_cnt = srclen / SIZEOF_HSWORD;
462
463 if (msbf) {
464 if (limb_cnt_rem) { // partial limb
465 r[limb_cnt] = integer_gmp_mpn_import1(srcptr, 0, limb_cnt_rem, 1);
466 srcptr += limb_cnt_rem;
467 }
468
469 for (unsigned ri = 0; ri < limb_cnt; ++ri) {
470 r[limb_cnt-ri-1] = integer_gmp_mpn_import1(srcptr, 0, SIZEOF_HSWORD, 1);
471 srcptr += SIZEOF_HSWORD;
472 }
473 } else { // lsbf
474 for (unsigned ri = 0; ri < limb_cnt; ++ri) {
475 r[ri] = integer_gmp_mpn_import1(srcptr, 0, SIZEOF_HSWORD, 0);
476 srcptr += SIZEOF_HSWORD;
477 }
478
479 if (limb_cnt_rem) // partial limb
480 r[limb_cnt] = integer_gmp_mpn_import1(srcptr, 0, limb_cnt_rem, 0);
481 }
482 }
483
484 /* Scan for first non-zero byte starting at srcptr[srcofs], ending at
485 * srcptr[srcofs+srclen-1];
486 *
487 * If no non-zero byte found, returns srcofs+srclen; otherwise returns
488 * index of srcptr where first non-zero byte was found.
489 */
490 HsWord
491 integer_gmp_scan_nzbyte(const uint8_t *srcptr,
492 const HsWord srcofs, const HsWord srclen)
493 {
494 // TODO: consider implementing this function in Haskell-land
495 srcptr += srcofs;
496
497 for (unsigned i = 0; i < srclen; ++i)
498 if (srcptr[i])
499 return srcofs+i;
500
501 return srcofs+srclen;
502 }
503
504 /* Reverse scan for non-zero byte
505 * starting at srcptr[srcofs+srclen-1], ending at srcptr[srcofs].
506 *
507 * Returns new length srclen1 such that srcptr[srcofs+i] == 0 for
508 * srclen1 <= i < srclen.
509 */
510 HsWord
511 integer_gmp_rscan_nzbyte(const uint8_t *srcptr,
512 const HsWord srcofs, const HsWord srclen)
513 {
514 // TODO: consider implementing this function in Haskell-land
515 srcptr += srcofs;
516
517 for (unsigned i = srclen; i > 0; --i)
518 if (srcptr[i-1])
519 return i;
520
521 return 0;
522 }
523
524 /* wrapper around mpz_probab_prime_p */
525 HsInt
526 integer_gmp_test_prime(const mp_limb_t s[], const mp_size_t sn, const HsInt rep)
527 {
528 if (mp_limb_zero_p(s,sn)) return 0;
529
530 const mpz_t sz = CONST_MPZ_INIT(s, sn);
531
532 // int mpz_probab_prime_p (const mpz_t n, int reps)
533 return mpz_probab_prime_p(sz, rep);
534 }
535
536 /* wrapper around mpz_probab_prime_p */
537 HsInt
538 integer_gmp_test_prime1(const mp_limb_t limb, const HsInt rep)
539 {
540 if (!limb) return 0;
541
542 return integer_gmp_test_prime(&limb, 1, rep);
543 }
544
545 /* wrapper around mpz_nextprime()
546 *
547 * Stores next prime (relative to {sp,sn}) in {rp,sn}.
548 * Return value is most significant limb of {rp,sn+1}.
549 */
550 mp_limb_t
551 integer_gmp_next_prime(mp_limb_t rp[], const mp_limb_t sp[],
552 const mp_size_t sn)
553 {
554 assert (sn>=0);
555
556 if (!sn) return 2;
557 if (sn == 1 && sp[0] < 2) {
558 rp[0] = 2;
559 return 0;
560 }
561
562 const mpz_t op = CONST_MPZ_INIT(sp, sn);
563
564 mpz_t rop;
565 mpz_init (rop);
566 mpz_nextprime (rop, op);
567
568 const mp_size_t rn = rop[0]._mp_size;
569
570 // copy result into {rp,sn} buffer
571 assert (rn == sn || rn == sn+1);
572 memcpy(rp, rop[0]._mp_d, sn*sizeof(mp_limb_t));
573 const mp_limb_t result = rn>sn ? rop[0]._mp_d[sn] : 0;
574
575 mpz_clear (rop);
576
577 return result;
578 }
579
580 /* wrapper around mpz_nextprime()
581 *
582 * returns next prime modulo 2^GMP_LIMB_BITS
583 */
584 mp_limb_t
585 integer_gmp_next_prime1(const mp_limb_t limb)
586 {
587 if (limb < 2) return 2;
588
589 const mpz_t op = CONST_MPZ_INIT(&limb, 1);
590
591 mpz_t rop;
592 mpz_init (rop);
593 mpz_nextprime (rop, op);
594 assert (rop[0]._mp_size > 0);
595 const mp_limb_t result = rop[0]._mp_d[0];
596 mpz_clear (rop);
597
598 return result;
599 }
600
601 /* wrapper around mpz_powm()
602 *
603 * Store '(B^E) mod M' in {rp,rn}
604 *
605 * rp must have allocated mn limbs; This function's return value is
606 * the actual number rn (0 < rn <= mn) of limbs written to the rp limb-array.
607 *
608 * bn and en are allowed to be negative to denote negative numbers
609 */
610 mp_size_t
611 integer_gmp_powm(mp_limb_t rp[], // result
612 const mp_limb_t bp[], const mp_size_t bn, // base
613 const mp_limb_t ep[], const mp_size_t en, // exponent
614 const mp_limb_t mp[], const mp_size_t mn) // mod
615 {
616 assert(!mp_limb_zero_p(mp,mn));
617
618 if ((mn == 1 || mn == -1) && mp[0] == 1) {
619 rp[0] = 0;
620 return 1;
621 }
622
623 if (mp_limb_zero_p(ep,en)) {
624 rp[0] = 1;
625 return 1;
626 }
627
628 const mpz_t b = CONST_MPZ_INIT(bp, mp_limb_zero_p(bp,bn) ? 0 : bn);
629 const mpz_t e = CONST_MPZ_INIT(ep, mp_limb_zero_p(ep,en) ? 0 : en);
630 const mpz_t m = CONST_MPZ_INIT(mp, mn);
631
632 mpz_t r;
633 mpz_init (r);
634
635 mpz_powm(r, b, e, m);
636
637 const mp_size_t rn = r[0]._mp_size;
638
639 if (rn) {
640 assert(0 < rn && rn <= mn);
641 memcpy(rp, r[0]._mp_d, rn*sizeof(mp_limb_t));
642 }
643
644 mpz_clear (r);
645
646 if (!rn) {
647 rp[0] = 0;
648 return 1;
649 }
650
651 return rn;
652 }
653
654 /* version of integer_gmp_powm() for single-limb moduli */
655 mp_limb_t
656 integer_gmp_powm1(const mp_limb_t bp[], const mp_size_t bn, // base
657 const mp_limb_t ep[], const mp_size_t en, // exponent
658 const mp_limb_t m0) // mod
659 {
660 assert(m0);
661
662 if (m0==1) return 0;
663 if (mp_limb_zero_p(ep,en)) return 1;
664
665 const mpz_t b = CONST_MPZ_INIT(bp, mp_limb_zero_p(bp,bn) ? 0 : bn);
666 const mpz_t e = CONST_MPZ_INIT(ep, en);
667 const mpz_t m = CONST_MPZ_INIT(&m0, !!m0);
668
669 mpz_t r;
670 mpz_init (r);
671 mpz_powm(r, b, e, m);
672
673 assert(r[0]._mp_size == 0 || r[0]._mp_size == 1);
674 const mp_limb_t result = r[0]._mp_size ? r[0]._mp_d[0] : 0;
675
676 mpz_clear (r);
677
678 return result;
679 }
680
681 /* version of integer_gmp_powm() for single-limb arguments */
682 mp_limb_t
683 integer_gmp_powm_word(const mp_limb_t b0, // base
684 const mp_limb_t e0, // exponent
685 const mp_limb_t m0) // mod
686 {
687 return integer_gmp_powm1(&b0, !!b0, &e0, !!e0, m0);
688 }
689
690
691 /* wrapper around mpz_invert()
692 *
693 * Store '(1/X) mod abs(M)' in {rp,rn}
694 *
695 * rp must have allocated mn limbs; This function's return value is
696 * the actual number rn (0 < rn <= mn) of limbs written to the rp limb-array.
697 *
698 * Returns 0 if inverse does not exist.
699 */
700 mp_size_t
701 integer_gmp_invert(mp_limb_t rp[], // result
702 const mp_limb_t xp[], const mp_size_t xn, // base
703 const mp_limb_t mp[], const mp_size_t mn) // mod
704 {
705 if (mp_limb_zero_p(xp,xn)
706 || mp_limb_zero_p(mp,mn)
707 || ((mn == 1 || mn == -1) && mp[0] == 1)) {
708 rp[0] = 0;
709 return 1;
710 }
711
712 const mpz_t x = CONST_MPZ_INIT(xp, xn);
713 const mpz_t m = CONST_MPZ_INIT(mp, mn);
714
715 mpz_t r;
716 mpz_init (r);
717
718 const int inv_exists = mpz_invert(r, x, m);
719
720 const mp_size_t rn = inv_exists ? r[0]._mp_size : 0;
721
722 if (rn) {
723 assert(0 < rn && rn <= mn);
724 memcpy(rp, r[0]._mp_d, rn*sizeof(mp_limb_t));
725 }
726
727 mpz_clear (r);
728
729 if (!rn) {
730 rp[0] = 0;
731 return 1;
732 }
733
734 return rn;
735 }
736
737
738 /* Version of integer_gmp_invert() operating on single limbs */
739 mp_limb_t
740 integer_gmp_invert_word(const mp_limb_t x0, const mp_limb_t m0)
741 {
742 if (!x0 || m0<=1) return 0;
743 if (x0 == 1) return 1;
744
745 const mpz_t x = CONST_MPZ_INIT(&x0, 1);
746 const mpz_t m = CONST_MPZ_INIT(&m0, 1);
747
748 mpz_t r;
749 mpz_init (r);
750
751 const int inv_exists = mpz_invert(r, x, m);
752 const mp_size_t rn = inv_exists ? r[0]._mp_size : 0;
753
754 assert (rn == 0 || rn == 1);
755 const mp_limb_t r0 = rn ? r[0]._mp_d[0] : 0;
756
757 mpz_clear (r);
758
759 return r0;
760 }
761
762
763 /* Wrappers for GMP 4.x compat
764 *
765 * In GMP 5.0 the following operations were added:
766 *
767 * mpn_sqr, mpn_and_n, mpn_ior_n, mpn_xor_n, mpn_nand_n, mpn_nior_n,
768 * mpn_xnor_n, mpn_andn_n, mpn_iorn_n, mpn_com, mpn_neg, mpn_copyi,
769 * mpn_copyd, mpn_zero
770 *
771 * We use some of those, but for GMP 4.x compatibility we need to
772 * emulate those (while incurring some overhead).
773 */
774 #if __GNU_MP_VERSION < 5
775
776 #define MPN_LOGIC_OP_WRAPPER(MPN_WRAPPER, MPZ_OP) \
777 void \
778 MPN_WRAPPER(mp_limb_t *rp, const mp_limb_t *s1p, \
779 const mp_limb_t *s2p, mp_size_t n) \
780 { \
781 assert(n > 0); \
782 \
783 const mpz_t s1 = CONST_MPZ_INIT(s1p, n); \
784 const mpz_t s2 = CONST_MPZ_INIT(s2p, n); \
785 \
786 mpz_t r; \
787 mpz_init (r); \
788 MPZ_OP (r, s1, s2); \
789 \
790 const mp_size_t rn = r[0]._mp_size; \
791 memset (rp, 0, n*sizeof(mp_limb_t)); \
792 memcpy (rp, r[0]._mp_d, mp_size_minabs(rn,n)*sizeof(mp_limb_t)); \
793 \
794 mpz_clear (r); \
795 }
796
797 static void
798 __mpz_andn(mpz_t r, const mpz_t s1, const mpz_t s2)
799 {
800 mpz_t s2c;
801 mpz_init (s2c);
802 mpz_com (s2c, s2);
803 mpz_and (r, s1, s2c);
804 mpz_clear (s2c);
805 }
806
807 MPN_LOGIC_OP_WRAPPER(integer_gmp_mpn_and_n, mpz_and)
808 MPN_LOGIC_OP_WRAPPER(integer_gmp_mpn_andn_n, __mpz_andn)
809 MPN_LOGIC_OP_WRAPPER(integer_gmp_mpn_ior_n, mpz_ior)
810 MPN_LOGIC_OP_WRAPPER(integer_gmp_mpn_xor_n, mpz_xor)
811
812 #else /* __GNU_MP_VERSION >= 5 */
813 void
814 integer_gmp_mpn_and_n(mp_limb_t *rp, const mp_limb_t *s1p,
815 const mp_limb_t *s2p, mp_size_t n)
816 {
817 mpn_and_n(rp, s1p, s2p, n);
818 }
819
820 void
821 integer_gmp_mpn_andn_n(mp_limb_t *rp, const mp_limb_t *s1p,
822 const mp_limb_t *s2p, mp_size_t n)
823 {
824 mpn_andn_n(rp, s1p, s2p, n);
825 }
826
827 void
828 integer_gmp_mpn_ior_n(mp_limb_t *rp, const mp_limb_t *s1p,
829 const mp_limb_t *s2p, mp_size_t n)
830 {
831 mpn_ior_n(rp, s1p, s2p, n);
832 }
833
834 void
835 integer_gmp_mpn_xor_n(mp_limb_t *rp, const mp_limb_t *s1p,
836 const mp_limb_t *s2p, mp_size_t n)
837 {
838 mpn_xor_n(rp, s1p, s2p, n);
839 }
840 #endif