6c188a3684d4ce396f74913226d70dc19843d26f
[ghc.git] / libraries / integer-gmp2 / cbits / wrappers.c
1 #define _ISOC99_SOURCE
2
3 #include <assert.h>
4 #include <stdbool.h>
5 #include <stdlib.h>
6 #include <stdint.h>
7 #include <string.h>
8 #include <math.h>
9 #include <float.h>
10 #include <stdio.h>
11
12 #include <gmp.h>
13
14 #include "HsFFI.h"
15 #include "MachDeps.h"
16
17 // GMP 4.x compatibility
18 #if !defined(__GNU_MP_VERSION)
19 # error __GNU_MP_VERSION not defined
20 #elif __GNU_MP_VERSION < 4
21 # error need GMP 4.0 or later
22 #elif __GNU_MP_VERSION < 5
23 typedef unsigned long int mp_bitcnt_t;
24 #endif
25
26 #if (GMP_NUMB_BITS) != (GMP_LIMB_BITS)
27 # error GMP_NUMB_BITS != GMP_LIMB_BITS not supported
28 #endif
29
30 #if (WORD_SIZE_IN_BITS) != (GMP_LIMB_BITS)
31 # error WORD_SIZE_IN_BITS != GMP_LIMB_BITS not supported
32 #endif
33
34 // sanity check
35 #if (SIZEOF_HSWORD*8) != WORD_SIZE_IN_BITS
36 # error (SIZEOF_HSWORD*8) != WORD_SIZE_IN_BITS
37 #endif
38
39 /* Perform arithmetic right shift on MPNs (multi-precision naturals)
40 *
41 * pre-conditions:
42 * - 0 < count < sn*GMP_NUMB_BITS
43 * - rn = sn - floor(count / GMP_NUMB_BITS)
44 * - sn > 0
45 *
46 * write {sp,sn} right-shifted by count bits into {rp,rn}
47 *
48 * return value: most-significant limb stored in {rp,rn} result
49 */
50 mp_limb_t
51 integer_gmp_mpn_rshift (mp_limb_t rp[], const mp_limb_t sp[], mp_size_t sn,
52 mp_bitcnt_t count)
53 {
54 const mp_size_t limb_shift = count / GMP_NUMB_BITS;
55 const unsigned int bit_shift = count % GMP_NUMB_BITS;
56 const mp_size_t rn = sn - limb_shift;
57
58 if (bit_shift)
59 mpn_rshift(rp, &sp[limb_shift], rn, bit_shift);
60 else
61 memcpy(rp, &sp[limb_shift], rn*sizeof(mp_limb_t));
62
63 return rp[rn-1];
64 }
65
66 /* Twos-complement version of 'integer_gmp_mpn_rshift' for performing
67 * arithmetic right shifts on "negative" MPNs.
68 *
69 * Same pre-conditions as 'integer_gmp_mpn_rshift'
70 *
71 * This variant is needed to operate on MPNs interpreted as negative
72 * numbers, which require "rounding" towards minus infinity iff a
73 * non-zero bit is shifted out.
74 */
75 mp_limb_t
76 integer_gmp_mpn_rshift_2c (mp_limb_t rp[], const mp_limb_t sp[],
77 const mp_size_t sn, const mp_bitcnt_t count)
78 {
79 const mp_size_t limb_shift = count / GMP_NUMB_BITS;
80 const unsigned int bit_shift = count % GMP_NUMB_BITS;
81 const mp_size_t rn = sn - limb_shift;
82
83 // whether non-zero bits were shifted out
84 bool nz_shift_out = false;
85
86 if (bit_shift) {
87 if (mpn_rshift(rp, &sp[limb_shift], rn, bit_shift))
88 nz_shift_out = true;
89 } else
90 memcpy(rp, &sp[limb_shift], rn*sizeof(mp_limb_t));
91
92 if (!nz_shift_out)
93 for (unsigned i = 0; i < limb_shift; i++)
94 if (sp[i]) {
95 nz_shift_out = true;
96 break;
97 }
98
99 // round if non-zero bits were shifted out
100 if (nz_shift_out)
101 if (mpn_add_1(rp, rp, rn, 1))
102 abort(); /* should never happen */
103
104 return rp[rn-1];
105 }
106
107 /* Perform left-shift operation on MPN
108 *
109 * pre-conditions:
110 * - 0 < count
111 * - rn = sn + ceil(count / GMP_NUMB_BITS)
112 * - sn > 0
113 *
114 * return value: most-significant limb stored in {rp,rn} result
115 */
116 mp_limb_t
117 integer_gmp_mpn_lshift (mp_limb_t rp[], const mp_limb_t sp[],
118 const mp_size_t sn, const mp_bitcnt_t count)
119 {
120 const mp_size_t limb_shift = count / GMP_NUMB_BITS;
121 const unsigned int bit_shift = count % GMP_NUMB_BITS;
122 const mp_size_t rn0 = sn + limb_shift;
123
124 memset(rp, 0, limb_shift*sizeof(mp_limb_t));
125 if (bit_shift) {
126 const mp_limb_t msl = mpn_lshift(&rp[limb_shift], sp, sn, bit_shift);
127 rp[rn0] = msl;
128 return msl;
129 } else {
130 memcpy(&rp[limb_shift], sp, sn*sizeof(mp_limb_t));
131 return rp[rn0-1];
132 }
133 }
134
135 /*
136 *
137 * sign of mp_size_t argument controls sign of converted double
138 */
139 HsDouble
140 integer_gmp_mpn_get_d (const mp_limb_t sp[], const mp_size_t sn,
141 const HsInt exponent)
142 {
143 if (sn == 0)
144 return 0.0; // should not happen
145
146 if (sn == 1 && sp[0] == 0)
147 return 0.0;
148
149 __mpz_struct const mpz = {
150 ._mp_alloc = abs(sn),
151 ._mp_size = sn,
152 ._mp_d = (mp_limb_t*)sp
153 };
154
155 if (!exponent)
156 return mpz_get_d(&mpz);
157
158 long e = 0;
159 double d = mpz_get_d_2exp (&e, &mpz);
160
161 // TODO: over/underflow handling?
162 return ldexp(d, e+exponent);
163 }
164
165 mp_limb_t
166 integer_gmp_gcd_word(const mp_limb_t x, const mp_limb_t y)
167 {
168 if (!x) return y;
169 if (!y) return x;
170
171 return mpn_gcd_1(&x, 1, y);
172 }
173
174 mp_limb_t
175 integer_gmp_mpn_gcd_1(const mp_limb_t x[], const mp_size_t xn,
176 const mp_limb_t y)
177 {
178 assert (xn > 0);
179 assert (xn == 1 || y != 0);
180
181 if (xn == 1)
182 return integer_gmp_gcd_word(x[0], y);
183
184 return mpn_gcd_1(x, xn, y);
185 }
186
187
188 mp_size_t
189 integer_gmp_mpn_gcd(mp_limb_t r[],
190 const mp_limb_t x0[], const mp_size_t xn,
191 const mp_limb_t y0[], const mp_size_t yn)
192 {
193 assert (xn >= yn);
194 assert (yn > 0);
195 assert (xn == yn || yn > 1 || y0[0] != 0);
196 /* post-condition: rn <= xn */
197
198 if (yn == 1) {
199 if (y0[0]) {
200 r[0] = integer_gmp_mpn_gcd_1(x0, xn, y0[0]);
201 return 1;
202 } else { /* {y0,yn} == 0 */
203 assert (xn==yn); /* NB: redundant assertion */
204 memcpy(r, x0, xn*sizeof(mp_limb_t));
205 return xn;
206 }
207 } else {
208 // mpn_gcd() seems to require non-trivial normalization of its
209 // input arguments (which does not seem to be documented anywhere,
210 // see source of mpz_gcd() for more details), so we resort to just
211 // use mpz_gcd() which does the tiresome normalization for us at
212 // the cost of a few additional temporary buffer allocations in
213 // C-land.
214
215 const mpz_t op1 = {{
216 ._mp_alloc = xn,
217 ._mp_size = xn,
218 ._mp_d = (mp_limb_t*)x0
219 }};
220
221 const mpz_t op2 = {{
222 ._mp_alloc = yn,
223 ._mp_size = yn,
224 ._mp_d = (mp_limb_t*)y0
225 }};
226
227 mpz_t rop;
228 mpz_init (rop);
229
230 mpz_gcd(rop, op1, op2);
231
232 const mp_size_t rn = rop[0]._mp_size;
233 assert(rn > 0);
234 assert(rn <= xn);
235
236 /* the allocation/memcpy of the result can be neglectable since
237 mpz_gcd() already has to allocate other temporary buffers
238 anyway */
239 memcpy(r, rop[0]._mp_d, rn*sizeof(mp_limb_t));
240
241 mpz_clear(rop);
242
243 return rn;
244 }
245 }
246
247 /* Truncating (i.e. rounded towards zero) integer division-quotient of MPN */
248 void
249 integer_gmp_mpn_tdiv_q (mp_limb_t q[],
250 const mp_limb_t n[], const mp_size_t nn,
251 const mp_limb_t d[], const mp_size_t dn)
252 {
253 /* qn = 1+nn-dn; rn = dn */
254 assert(nn>=dn);
255
256 if (dn > 128) {
257 // Use temporary heap allocated throw-away buffer for MPNs larger
258 // than 1KiB for 64bit-sized limbs (larger than 512bytes for
259 // 32bit-sized limbs)
260 mp_limb_t *const r = malloc(dn*sizeof(mp_limb_t));
261 mpn_tdiv_qr(q, r, 0, n, nn, d, dn);
262 free (r);
263 } else { // allocate smaller arrays on the stack
264 mp_limb_t r[dn];
265 mpn_tdiv_qr(q, r, 0, n, nn, d, dn);
266 }
267 }
268
269 /* Truncating (i.e. rounded towards zero) integer division-remainder of MPNs */
270 void
271 integer_gmp_mpn_tdiv_r (mp_limb_t r[],
272 const mp_limb_t n[], const mp_size_t nn,
273 const mp_limb_t d[], const mp_size_t dn)
274 {
275 /* qn = 1+nn-dn; rn = dn */
276 assert(nn>=dn);
277 const mp_size_t qn = 1+nn-dn;
278
279 if (qn > 128) {
280 // Use temporary heap allocated throw-away buffer for MPNs larger
281 // than 1KiB for 64bit-sized limbs (larger than 512bytes for
282 // 32bit-sized limbs)
283 mp_limb_t *const q = malloc(qn*sizeof(mp_limb_t));
284 mpn_tdiv_qr(q, r, 0, n, nn, d, dn);
285 free (q);
286 } else { // allocate smaller arrays on the stack
287 mp_limb_t q[qn];
288 mpn_tdiv_qr(q, r, 0, n, nn, d, dn);
289 }
290 }
291
292
293 /* Wraps GMP's 'mpz_sizeinbase()' function */
294 HsWord
295 integer_gmp_mpn_sizeinbase(const mp_limb_t s[], const mp_size_t sn,
296 const HsInt base)
297 {
298 assert (2 <= base && base <= 256);
299
300 if (!sn) return 1;
301
302 const mpz_t zs = {{
303 ._mp_alloc = sn,
304 ._mp_size = sn,
305 ._mp_d = (mp_limb_t*)s
306 }};
307
308 return mpz_sizeinbase(zs, base);
309 }
310
311 /* Single-limb version of 'integer_gmp_mpn_sizeinbase()' */
312 HsWord
313 integer_gmp_mpn_sizeinbase1(const mp_limb_t s, const HsInt base)
314 {
315 return s ? integer_gmp_mpn_sizeinbase(&s, 1, base) : 1;
316 }
317
318 /* Wrapper around GMP's 'mpz_export()' function */
319 HsWord
320 integer_gmp_mpn_export(const mp_limb_t s[], const mp_size_t sn,
321 void *destptr, HsInt destofs, HsInt msbf)
322 {
323 /* TODO: implement w/o GMP, c.f. 'integer_gmp_mpn_import()' */
324 assert (msbf == 0 || msbf == 1);
325
326 if (!sn || (sn == 1 && !s[0]))
327 return 0;
328
329 const mpz_t zs = {{
330 ._mp_alloc = sn,
331 ._mp_size = sn,
332 ._mp_d = (mp_limb_t*)s
333 }};
334
335 size_t written = 0;
336
337 // mpz_export (void *rop, size_t *countp, int order, size_t size, int endian,
338 // size_t nails, const mpz_t op)
339 (void) mpz_export(((char *)destptr)+destofs, &written, !msbf ? -1 : 1,
340 /* size */ 1, /* endian */ 0, /* nails */ 0, zs);
341
342 return written;
343 }
344
345 /* Single-limb version of 'integer_gmp_mpn_export()' */
346 HsWord
347 integer_gmp_mpn_export1(const mp_limb_t s,
348 void *destptr, const HsInt destofs, const HsInt msbf)
349 {
350 /* TODO: implement w/o GMP */
351 return integer_gmp_mpn_export(&s, 1, destptr, destofs, msbf);
352 }
353
354 /* Import single limb from memory location
355 *
356 * We can't use GMP's 'mpz_import()'
357 */
358 inline HsWord
359 integer_gmp_mpn_import1(const uint8_t *srcptr, const HsWord srcofs,
360 const HsWord srclen, const HsInt msbf)
361 {
362 assert (msbf == 0 || msbf == 1);
363 assert (srclen <= SIZEOF_HSWORD);
364
365 srcptr += srcofs;
366
367 HsWord result = 0;
368
369 if (msbf)
370 for (unsigned i = 0; i < srclen; ++i)
371 result |= (HsWord)srcptr[i] << ((srclen-i-1)*8);
372 else // lsbf
373 for (unsigned i = 0; i < srclen; ++i)
374 result |= (HsWord)srcptr[i] << (i*8);
375
376 return result;
377 }
378
379 /* import into mp_limb_t[] from memory location */
380 void
381 integer_gmp_mpn_import(mp_limb_t * restrict r, const uint8_t * restrict srcptr,
382 const HsWord srcofs, const HsWord srclen,
383 const HsInt msbf)
384 {
385 assert (msbf == 0 || msbf == 1);
386
387 srcptr += srcofs;
388
389 const unsigned limb_cnt_rem = srclen % SIZEOF_HSWORD;
390 const mp_size_t limb_cnt = srclen / SIZEOF_HSWORD;
391
392 if (msbf) {
393 if (limb_cnt_rem) { // partial limb
394 r[limb_cnt] = integer_gmp_mpn_import1(srcptr, 0, limb_cnt_rem, 1);
395 srcptr += limb_cnt_rem;
396 }
397
398 for (unsigned ri = 0; ri < limb_cnt; ++ri) {
399 r[limb_cnt-ri-1] = integer_gmp_mpn_import1(srcptr, 0, SIZEOF_HSWORD, 1);
400 srcptr += SIZEOF_HSWORD;
401 }
402 } else { // lsbf
403 for (unsigned ri = 0; ri < limb_cnt; ++ri) {
404 r[ri] = integer_gmp_mpn_import1(srcptr, 0, SIZEOF_HSWORD, 0);
405 srcptr += SIZEOF_HSWORD;
406 }
407
408 if (limb_cnt_rem) // partial limb
409 r[limb_cnt] = integer_gmp_mpn_import1(srcptr, 0, limb_cnt_rem, 0);
410 }
411 }
412
413 /* Scan for first non-zero byte starting at srcptr[srcofs], ending at
414 * srcptr[srcofs+srclen-1];
415 *
416 * If no non-zero byte found, returns srcofs+srclen; otherwise returns
417 * index of srcptr where first non-zero byte was found.
418 */
419 HsWord
420 integer_gmp_scan_nzbyte(const uint8_t *srcptr,
421 const HsWord srcofs, const HsWord srclen)
422 {
423 // TODO: consider implementing this function in Haskell-land
424 srcptr += srcofs;
425
426 for (unsigned i = 0; i < srclen; ++i)
427 if (srcptr[i])
428 return srcofs+i;
429
430 return srcofs+srclen;
431 }
432
433 /* Reverse scan for non-zero byte
434 * starting at srcptr[srcofs+srclen-1], ending at srcptr[srcofs].
435 *
436 * Returns new length srclen1 such that srcptr[srcofs+i] == 0 for
437 * srclen1 <= i < srclen.
438 */
439 HsWord
440 integer_gmp_rscan_nzbyte(const uint8_t *srcptr,
441 const HsWord srcofs, const HsWord srclen)
442 {
443 // TODO: consider implementing this function in Haskell-land
444 srcptr += srcofs;
445
446 for (unsigned i = srclen; i > 0; --i)
447 if (srcptr[i-1])
448 return i;
449
450 return 0;
451 }
452
453 /* wrapper around mpz_probab_prime_p */
454 HsInt
455 integer_gmp_test_prime(const mp_limb_t s[], const mp_size_t sn, const HsInt rep)
456 {
457 if (!sn) return 0;
458
459 const mpz_t sz = {{
460 ._mp_alloc = sn,
461 ._mp_size = sn,
462 ._mp_d = (mp_limb_t*)s
463 }};
464
465 // int mpz_probab_prime_p (const mpz_t n, int reps)
466 return mpz_probab_prime_p(sz, rep);
467 }
468
469 /* wrapper around mpz_probab_prime_p */
470 HsInt
471 integer_gmp_test_prime1(const mp_limb_t limb, const HsInt rep)
472 {
473 if (!limb) return 0;
474
475 return integer_gmp_test_prime(&limb, 1, rep);
476 }