1 /* mini-gmp, a minimalistic implementation of a GNU GMP subset.
2
3 Contributed to the GNU project by Niels Möller
4
5 Copyright 1991-1997, 1999-2016 Free Software Foundation, Inc.
6
7 This file is part of the GNU MP Library.
8
9 The GNU MP Library is free software; you can redistribute it and/or modify
10 it under the terms of either:
11
12 * the GNU Lesser General Public License as published by the Free
13 Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
15
16 or
17
18 * the GNU General Public License as published by the Free Software
19 Foundation; either version 2 of the License, or (at your option) any
20 later version.
21
22 or both in parallel, as here.
23
24 The GNU MP Library is distributed in the hope that it will be useful, but
25 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
26 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
27 for more details.
28
29 You should have received copies of the GNU General Public License and the
30 GNU Lesser General Public License along with the GNU MP Library. If not,
31 see https://www.gnu.org/licenses/. */
32
33 /* NOTE: All functions in this file which are not declared in
34 mini-gmp.h are internal, and are not intended to be compatible
35 neither with GMP nor with future versions of mini-gmp. */
36
37 /* Much of the material copied from GMP files, including: gmp-impl.h,
38 longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c,
39 mpn/generic/lshift.c, mpn/generic/mul_1.c,
40 mpn/generic/mul_basecase.c, mpn/generic/rshift.c,
41 mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c,
42 mpn/generic/submul_1.c. */
43
44 #include <assert.h>
45 #include <ctype.h>
46 #include <limits.h>
47 #include <stdio.h>
48 #include <stdlib.h>
49 #include <string.h>
50
51 #include "mini-gmp.h"
52
53
54 /* Macros */
55 #define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
56
57 #define GMP_LIMB_MAX (~ (mp_limb_t) 0)
58 #define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1))
59
60 #define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2))
61 #define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1)
62
63 #define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT)
64 #define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1))
65
66 #define GMP_ABS(x) ((x) >= 0 ? (x) : -(x))
67 #define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1))
68
69 #define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
70 #define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))
71
72 #define GMP_CMP(a,b) (((a) > (b)) - ((a) < (b)))
73
74 #define gmp_assert_nocarry(x) do { \
75 mp_limb_t __cy = (x); \
76 assert (__cy == 0); \
77 } while (0)
78
79 #define gmp_clz(count, x) do { \
80 mp_limb_t __clz_x = (x); \
81 unsigned __clz_c; \
82 for (__clz_c = 0; \
83 (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \
84 __clz_c += 8) \
85 __clz_x <<= 8; \
86 for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++) \
87 __clz_x <<= 1; \
88 (count) = __clz_c; \
89 } while (0)
90
91 #define gmp_ctz(count, x) do { \
92 mp_limb_t __ctz_x = (x); \
93 unsigned __ctz_c = 0; \
94 gmp_clz (__ctz_c, __ctz_x & - __ctz_x); \
95 (count) = GMP_LIMB_BITS - 1 - __ctz_c; \
96 } while (0)
97
98 #define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \
99 do { \
100 mp_limb_t __x; \
101 __x = (al) + (bl); \
102 (sh) = (ah) + (bh) + (__x < (al)); \
103 (sl) = __x; \
104 } while (0)
105
106 #define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \
107 do { \
108 mp_limb_t __x; \
109 __x = (al) - (bl); \
110 (sh) = (ah) - (bh) - ((al) < (bl)); \
111 (sl) = __x; \
112 } while (0)
113
114 #define gmp_umul_ppmm(w1, w0, u, v) \
115 do { \
116 mp_limb_t __x0, __x1, __x2, __x3; \
117 unsigned __ul, __vl, __uh, __vh; \
118 mp_limb_t __u = (u), __v = (v); \
119 \
120 __ul = __u & GMP_LLIMB_MASK; \
121 __uh = __u >> (GMP_LIMB_BITS / 2); \
122 __vl = __v & GMP_LLIMB_MASK; \
123 __vh = __v >> (GMP_LIMB_BITS / 2); \
124 \
125 __x0 = (mp_limb_t) __ul * __vl; \
126 __x1 = (mp_limb_t) __ul * __vh; \
127 __x2 = (mp_limb_t) __uh * __vl; \
128 __x3 = (mp_limb_t) __uh * __vh; \
129 \
130 __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \
131 __x1 += __x2; /* but this indeed can */ \
132 if (__x1 < __x2) /* did we get it? */ \
133 __x3 += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */ \
134 \
135 (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2)); \
136 (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK); \
137 } while (0)
138
139 #define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
140 do { \
141 mp_limb_t _qh, _ql, _r, _mask; \
142 gmp_umul_ppmm (_qh, _ql, (nh), (di)); \
143 gmp_add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl)); \
144 _r = (nl) - _qh * (d); \
145 _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \
146 _qh += _mask; \
147 _r += _mask & (d); \
148 if (_r >= (d)) \
149 { \
150 _r -= (d); \
151 _qh++; \
152 } \
153 \
154 (r) = _r; \
155 (q) = _qh; \
156 } while (0)
157
158 #define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv) \
159 do { \
160 mp_limb_t _q0, _t1, _t0, _mask; \
161 gmp_umul_ppmm ((q), _q0, (n2), (dinv)); \
162 gmp_add_ssaaaa ((q), _q0, (q), _q0, (n2), (n1)); \
163 \
164 /* Compute the two most significant limbs of n - q'd */ \
165 (r1) = (n1) - (d1) * (q); \
166 gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0)); \
167 gmp_umul_ppmm (_t1, _t0, (d0), (q)); \
168 gmp_sub_ddmmss ((r1), (r0), (r1), (r0), _t1, _t0); \
169 (q)++; \
170 \
171 /* Conditionally adjust q and the remainders */ \
172 _mask = - (mp_limb_t) ((r1) >= _q0); \
173 (q) += _mask; \
174 gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \
175 if ((r1) >= (d1)) \
176 { \
177 if ((r1) > (d1) || (r0) >= (d0)) \
178 { \
179 (q)++; \
180 gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0)); \
181 } \
182 } \
183 } while (0)
184
185 /* Swap macros. */
186 #define MP_LIMB_T_SWAP(x, y) \
187 do { \
188 mp_limb_t __mp_limb_t_swap__tmp = (x); \
189 (x) = (y); \
190 (y) = __mp_limb_t_swap__tmp; \
191 } while (0)
192 #define MP_SIZE_T_SWAP(x, y) \
193 do { \
194 mp_size_t __mp_size_t_swap__tmp = (x); \
195 (x) = (y); \
196 (y) = __mp_size_t_swap__tmp; \
197 } while (0)
198 #define MP_BITCNT_T_SWAP(x,y) \
199 do { \
200 mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x); \
201 (x) = (y); \
202 (y) = __mp_bitcnt_t_swap__tmp; \
203 } while (0)
204 #define MP_PTR_SWAP(x, y) \
205 do { \
206 mp_ptr __mp_ptr_swap__tmp = (x); \
207 (x) = (y); \
208 (y) = __mp_ptr_swap__tmp; \
209 } while (0)
210 #define MP_SRCPTR_SWAP(x, y) \
211 do { \
212 mp_srcptr __mp_srcptr_swap__tmp = (x); \
213 (x) = (y); \
214 (y) = __mp_srcptr_swap__tmp; \
215 } while (0)
216
217 #define MPN_PTR_SWAP(xp,xs, yp,ys) \
218 do { \
219 MP_PTR_SWAP (xp, yp); \
220 MP_SIZE_T_SWAP (xs, ys); \
221 } while(0)
222 #define MPN_SRCPTR_SWAP(xp,xs, yp,ys) \
223 do { \
224 MP_SRCPTR_SWAP (xp, yp); \
225 MP_SIZE_T_SWAP (xs, ys); \
226 } while(0)
227
228 #define MPZ_PTR_SWAP(x, y) \
229 do { \
230 mpz_ptr __mpz_ptr_swap__tmp = (x); \
231 (x) = (y); \
232 (y) = __mpz_ptr_swap__tmp; \
233 } while (0)
234 #define MPZ_SRCPTR_SWAP(x, y) \
235 do { \
236 mpz_srcptr __mpz_srcptr_swap__tmp = (x); \
237 (x) = (y); \
238 (y) = __mpz_srcptr_swap__tmp; \
239 } while (0)
240
241 const int mp_bits_per_limb = GMP_LIMB_BITS;
242
243
244 /* Memory allocation and other helper functions. */
245 static void
gmp_die(const char * msg)246 gmp_die (const char *msg)
247 {
248 fprintf (stderr, "%s\n", msg);
249 abort();
250 }
251
252 static void *
gmp_default_alloc(size_t size)253 gmp_default_alloc (size_t size)
254 {
255 void *p;
256
257 assert (size > 0);
258
259 p = malloc (size);
260 if (!p)
261 gmp_die("gmp_default_alloc: Virtual memory exhausted.");
262
263 return p;
264 }
265
266 static void *
gmp_default_realloc(void * old,size_t old_size,size_t new_size)267 gmp_default_realloc (void *old, size_t old_size, size_t new_size)
268 {
269 void * p;
270
271 p = realloc (old, new_size);
272
273 if (!p)
274 gmp_die("gmp_default_realloc: Virtual memory exhausted.");
275
276 return p;
277 }
278
279 static void
gmp_default_free(void * p,size_t size)280 gmp_default_free (void *p, size_t size)
281 {
282 free (p);
283 }
284
285 static void * (*gmp_allocate_func) (size_t) = gmp_default_alloc;
286 static void * (*gmp_reallocate_func) (void *, size_t, size_t) = gmp_default_realloc;
287 static void (*gmp_free_func) (void *, size_t) = gmp_default_free;
288
289 void
mp_get_memory_functions(void * (** alloc_func)(size_t),void * (** realloc_func)(void *,size_t,size_t),void (** free_func)(void *,size_t))290 mp_get_memory_functions (void *(**alloc_func) (size_t),
291 void *(**realloc_func) (void *, size_t, size_t),
292 void (**free_func) (void *, size_t))
293 {
294 if (alloc_func)
295 *alloc_func = gmp_allocate_func;
296
297 if (realloc_func)
298 *realloc_func = gmp_reallocate_func;
299
300 if (free_func)
301 *free_func = gmp_free_func;
302 }
303
304 void
mp_set_memory_functions(void * (* alloc_func)(size_t),void * (* realloc_func)(void *,size_t,size_t),void (* free_func)(void *,size_t))305 mp_set_memory_functions (void *(*alloc_func) (size_t),
306 void *(*realloc_func) (void *, size_t, size_t),
307 void (*free_func) (void *, size_t))
308 {
309 if (!alloc_func)
310 alloc_func = gmp_default_alloc;
311 if (!realloc_func)
312 realloc_func = gmp_default_realloc;
313 if (!free_func)
314 free_func = gmp_default_free;
315
316 gmp_allocate_func = alloc_func;
317 gmp_reallocate_func = realloc_func;
318 gmp_free_func = free_func;
319 }
320
321 #define gmp_xalloc(size) ((*gmp_allocate_func)((size)))
322 #define gmp_free(p) ((*gmp_free_func) ((p), 0))
323
324 static mp_ptr
gmp_xalloc_limbs(mp_size_t size)325 gmp_xalloc_limbs (mp_size_t size)
326 {
327 return (mp_ptr) gmp_xalloc (size * sizeof (mp_limb_t));
328 }
329
330 static mp_ptr
gmp_xrealloc_limbs(mp_ptr old,mp_size_t size)331 gmp_xrealloc_limbs (mp_ptr old, mp_size_t size)
332 {
333 assert (size > 0);
334 return (mp_ptr) (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t));
335 }
336
337
338 /* MPN interface */
339
340 void
mpn_copyi(mp_ptr d,mp_srcptr s,mp_size_t n)341 mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n)
342 {
343 mp_size_t i;
344 for (i = 0; i < n; i++)
345 d[i] = s[i];
346 }
347
348 void
mpn_copyd(mp_ptr d,mp_srcptr s,mp_size_t n)349 mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
350 {
351 while (--n >= 0)
352 d[n] = s[n];
353 }
354
355 int
mpn_cmp(mp_srcptr ap,mp_srcptr bp,mp_size_t n)356 mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
357 {
358 while (--n >= 0)
359 {
360 if (ap[n] != bp[n])
361 return ap[n] > bp[n] ? 1 : -1;
362 }
363 return 0;
364 }
365
366 static int
mpn_cmp4(mp_srcptr ap,mp_size_t an,mp_srcptr bp,mp_size_t bn)367 mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
368 {
369 if (an != bn)
370 return an < bn ? -1 : 1;
371 else
372 return mpn_cmp (ap, bp, an);
373 }
374
375 static mp_size_t
mpn_normalized_size(mp_srcptr xp,mp_size_t n)376 mpn_normalized_size (mp_srcptr xp, mp_size_t n)
377 {
378 while (n > 0 && xp[n-1] == 0)
379 --n;
380 return n;
381 }
382
383 int
mpn_zero_p(mp_srcptr rp,mp_size_t n)384 mpn_zero_p(mp_srcptr rp, mp_size_t n)
385 {
386 return mpn_normalized_size (rp, n) == 0;
387 }
388
389 void
mpn_zero(mp_ptr rp,mp_size_t n)390 mpn_zero (mp_ptr rp, mp_size_t n)
391 {
392 while (--n >= 0)
393 rp[n] = 0;
394 }
395
396 mp_limb_t
mpn_add_1(mp_ptr rp,mp_srcptr ap,mp_size_t n,mp_limb_t b)397 mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
398 {
399 mp_size_t i;
400
401 assert (n > 0);
402 i = 0;
403 do
404 {
405 mp_limb_t r = ap[i] + b;
406 /* Carry out */
407 b = (r < b);
408 rp[i] = r;
409 }
410 while (++i < n);
411
412 return b;
413 }
414
415 mp_limb_t
mpn_add_n(mp_ptr rp,mp_srcptr ap,mp_srcptr bp,mp_size_t n)416 mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
417 {
418 mp_size_t i;
419 mp_limb_t cy;
420
421 for (i = 0, cy = 0; i < n; i++)
422 {
423 mp_limb_t a, b, r;
424 a = ap[i]; b = bp[i];
425 r = a + cy;
426 cy = (r < cy);
427 r += b;
428 cy += (r < b);
429 rp[i] = r;
430 }
431 return cy;
432 }
433
434 mp_limb_t
mpn_add(mp_ptr rp,mp_srcptr ap,mp_size_t an,mp_srcptr bp,mp_size_t bn)435 mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
436 {
437 mp_limb_t cy;
438
439 assert (an >= bn);
440
441 cy = mpn_add_n (rp, ap, bp, bn);
442 if (an > bn)
443 cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy);
444 return cy;
445 }
446
447 mp_limb_t
mpn_sub_1(mp_ptr rp,mp_srcptr ap,mp_size_t n,mp_limb_t b)448 mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
449 {
450 mp_size_t i;
451
452 assert (n > 0);
453
454 i = 0;
455 do
456 {
457 mp_limb_t a = ap[i];
458 /* Carry out */
459 mp_limb_t cy = a < b;
460 rp[i] = a - b;
461 b = cy;
462 }
463 while (++i < n);
464
465 return b;
466 }
467
468 mp_limb_t
mpn_sub_n(mp_ptr rp,mp_srcptr ap,mp_srcptr bp,mp_size_t n)469 mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
470 {
471 mp_size_t i;
472 mp_limb_t cy;
473
474 for (i = 0, cy = 0; i < n; i++)
475 {
476 mp_limb_t a, b;
477 a = ap[i]; b = bp[i];
478 b += cy;
479 cy = (b < cy);
480 cy += (a < b);
481 rp[i] = a - b;
482 }
483 return cy;
484 }
485
486 mp_limb_t
mpn_sub(mp_ptr rp,mp_srcptr ap,mp_size_t an,mp_srcptr bp,mp_size_t bn)487 mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
488 {
489 mp_limb_t cy;
490
491 assert (an >= bn);
492
493 cy = mpn_sub_n (rp, ap, bp, bn);
494 if (an > bn)
495 cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy);
496 return cy;
497 }
498
499 mp_limb_t
mpn_mul_1(mp_ptr rp,mp_srcptr up,mp_size_t n,mp_limb_t vl)500 mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
501 {
502 mp_limb_t ul, cl, hpl, lpl;
503
504 assert (n >= 1);
505
506 cl = 0;
507 do
508 {
509 ul = *up++;
510 gmp_umul_ppmm (hpl, lpl, ul, vl);
511
512 lpl += cl;
513 cl = (lpl < cl) + hpl;
514
515 *rp++ = lpl;
516 }
517 while (--n != 0);
518
519 return cl;
520 }
521
522 mp_limb_t
mpn_addmul_1(mp_ptr rp,mp_srcptr up,mp_size_t n,mp_limb_t vl)523 mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
524 {
525 mp_limb_t ul, cl, hpl, lpl, rl;
526
527 assert (n >= 1);
528
529 cl = 0;
530 do
531 {
532 ul = *up++;
533 gmp_umul_ppmm (hpl, lpl, ul, vl);
534
535 lpl += cl;
536 cl = (lpl < cl) + hpl;
537
538 rl = *rp;
539 lpl = rl + lpl;
540 cl += lpl < rl;
541 *rp++ = lpl;
542 }
543 while (--n != 0);
544
545 return cl;
546 }
547
548 mp_limb_t
mpn_submul_1(mp_ptr rp,mp_srcptr up,mp_size_t n,mp_limb_t vl)549 mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
550 {
551 mp_limb_t ul, cl, hpl, lpl, rl;
552
553 assert (n >= 1);
554
555 cl = 0;
556 do
557 {
558 ul = *up++;
559 gmp_umul_ppmm (hpl, lpl, ul, vl);
560
561 lpl += cl;
562 cl = (lpl < cl) + hpl;
563
564 rl = *rp;
565 lpl = rl - lpl;
566 cl += lpl > rl;
567 *rp++ = lpl;
568 }
569 while (--n != 0);
570
571 return cl;
572 }
573
574 mp_limb_t
mpn_mul(mp_ptr rp,mp_srcptr up,mp_size_t un,mp_srcptr vp,mp_size_t vn)575 mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
576 {
577 assert (un >= vn);
578 assert (vn >= 1);
579
580 /* We first multiply by the low order limb. This result can be
581 stored, not added, to rp. We also avoid a loop for zeroing this
582 way. */
583
584 rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
585
586 /* Now accumulate the product of up[] and the next higher limb from
587 vp[]. */
588
589 while (--vn >= 1)
590 {
591 rp += 1, vp += 1;
592 rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
593 }
594 return rp[un];
595 }
596
597 void
mpn_mul_n(mp_ptr rp,mp_srcptr ap,mp_srcptr bp,mp_size_t n)598 mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
599 {
600 mpn_mul (rp, ap, n, bp, n);
601 }
602
603 void
mpn_sqr(mp_ptr rp,mp_srcptr ap,mp_size_t n)604 mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n)
605 {
606 mpn_mul (rp, ap, n, ap, n);
607 }
608
609 mp_limb_t
mpn_lshift(mp_ptr rp,mp_srcptr up,mp_size_t n,unsigned int cnt)610 mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
611 {
612 mp_limb_t high_limb, low_limb;
613 unsigned int tnc;
614 mp_limb_t retval;
615
616 assert (n >= 1);
617 assert (cnt >= 1);
618 assert (cnt < GMP_LIMB_BITS);
619
620 up += n;
621 rp += n;
622
623 tnc = GMP_LIMB_BITS - cnt;
624 low_limb = *--up;
625 retval = low_limb >> tnc;
626 high_limb = (low_limb << cnt);
627
628 while (--n != 0)
629 {
630 low_limb = *--up;
631 *--rp = high_limb | (low_limb >> tnc);
632 high_limb = (low_limb << cnt);
633 }
634 *--rp = high_limb;
635
636 return retval;
637 }
638
639 mp_limb_t
mpn_rshift(mp_ptr rp,mp_srcptr up,mp_size_t n,unsigned int cnt)640 mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
641 {
642 mp_limb_t high_limb, low_limb;
643 unsigned int tnc;
644 mp_limb_t retval;
645
646 assert (n >= 1);
647 assert (cnt >= 1);
648 assert (cnt < GMP_LIMB_BITS);
649
650 tnc = GMP_LIMB_BITS - cnt;
651 high_limb = *up++;
652 retval = (high_limb << tnc);
653 low_limb = high_limb >> cnt;
654
655 while (--n != 0)
656 {
657 high_limb = *up++;
658 *rp++ = low_limb | (high_limb << tnc);
659 low_limb = high_limb >> cnt;
660 }
661 *rp = low_limb;
662
663 return retval;
664 }
665
666 static mp_bitcnt_t
mpn_common_scan(mp_limb_t limb,mp_size_t i,mp_srcptr up,mp_size_t un,mp_limb_t ux)667 mpn_common_scan (mp_limb_t limb, mp_size_t i, mp_srcptr up, mp_size_t un,
668 mp_limb_t ux)
669 {
670 unsigned cnt;
671
672 assert (ux == 0 || ux == GMP_LIMB_MAX);
673 assert (0 <= i && i <= un );
674
675 while (limb == 0)
676 {
677 i++;
678 if (i == un)
679 return (ux == 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
680 limb = ux ^ up[i];
681 }
682 gmp_ctz (cnt, limb);
683 return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
684 }
685
686 mp_bitcnt_t
mpn_scan1(mp_srcptr ptr,mp_bitcnt_t bit)687 mpn_scan1 (mp_srcptr ptr, mp_bitcnt_t bit)
688 {
689 mp_size_t i;
690 i = bit / GMP_LIMB_BITS;
691
692 return mpn_common_scan ( ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
693 i, ptr, i, 0);
694 }
695
696 mp_bitcnt_t
mpn_scan0(mp_srcptr ptr,mp_bitcnt_t bit)697 mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit)
698 {
699 mp_size_t i;
700 i = bit / GMP_LIMB_BITS;
701
702 return mpn_common_scan (~ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
703 i, ptr, i, GMP_LIMB_MAX);
704 }
705
706 void
mpn_com(mp_ptr rp,mp_srcptr up,mp_size_t n)707 mpn_com (mp_ptr rp, mp_srcptr up, mp_size_t n)
708 {
709 while (--n >= 0)
710 *rp++ = ~ *up++;
711 }
712
713 mp_limb_t
mpn_neg(mp_ptr rp,mp_srcptr up,mp_size_t n)714 mpn_neg (mp_ptr rp, mp_srcptr up, mp_size_t n)
715 {
716 while (*up == 0)
717 {
718 *rp = 0;
719 if (!--n)
720 return 0;
721 ++up; ++rp;
722 }
723 *rp = - *up;
724 mpn_com (++rp, ++up, --n);
725 return 1;
726 }
727
728
729 /* MPN division interface. */
730
731 /* The 3/2 inverse is defined as
732
733 m = floor( (B^3-1) / (B u1 + u0)) - B
734 */
735 mp_limb_t
mpn_invert_3by2(mp_limb_t u1,mp_limb_t u0)736 mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
737 {
738 mp_limb_t r, p, m, ql;
739 unsigned ul, uh, qh;
740
741 assert (u1 >= GMP_LIMB_HIGHBIT);
742
743 /* For notation, let b denote the half-limb base, so that B = b^2.
744 Split u1 = b uh + ul. */
745 ul = u1 & GMP_LLIMB_MASK;
746 uh = u1 >> (GMP_LIMB_BITS / 2);
747
748 /* Approximation of the high half of quotient. Differs from the 2/1
749 inverse of the half limb uh, since we have already subtracted
750 u0. */
751 qh = ~u1 / uh;
752
753 /* Adjust to get a half-limb 3/2 inverse, i.e., we want
754
755 qh' = floor( (b^3 - 1) / u) - b = floor ((b^3 - b u - 1) / u
756 = floor( (b (~u) + b-1) / u),
757
758 and the remainder
759
760 r = b (~u) + b-1 - qh (b uh + ul)
761 = b (~u - qh uh) + b-1 - qh ul
762
763 Subtraction of qh ul may underflow, which implies adjustments.
764 But by normalization, 2 u >= B > qh ul, so we need to adjust by
765 at most 2.
766 */
767
768 r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK;
769
770 p = (mp_limb_t) qh * ul;
771 /* Adjustment steps taken from udiv_qrnnd_c */
772 if (r < p)
773 {
774 qh--;
775 r += u1;
776 if (r >= u1) /* i.e. we didn't get carry when adding to r */
777 if (r < p)
778 {
779 qh--;
780 r += u1;
781 }
782 }
783 r -= p;
784
785 /* Low half of the quotient is
786
787 ql = floor ( (b r + b-1) / u1).
788
789 This is a 3/2 division (on half-limbs), for which qh is a
790 suitable inverse. */
791
792 p = (r >> (GMP_LIMB_BITS / 2)) * qh + r;
793 /* Unlike full-limb 3/2, we can add 1 without overflow. For this to
794 work, it is essential that ql is a full mp_limb_t. */
795 ql = (p >> (GMP_LIMB_BITS / 2)) + 1;
796
797 /* By the 3/2 trick, we don't need the high half limb. */
798 r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1;
799
800 if (r >= (p << (GMP_LIMB_BITS / 2)))
801 {
802 ql--;
803 r += u1;
804 }
805 m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql;
806 if (r >= u1)
807 {
808 m++;
809 r -= u1;
810 }
811
812 /* Now m is the 2/1 invers of u1. If u0 > 0, adjust it to become a
813 3/2 inverse. */
814 if (u0 > 0)
815 {
816 mp_limb_t th, tl;
817 r = ~r;
818 r += u0;
819 if (r < u0)
820 {
821 m--;
822 if (r >= u1)
823 {
824 m--;
825 r -= u1;
826 }
827 r -= u1;
828 }
829 gmp_umul_ppmm (th, tl, u0, m);
830 r += th;
831 if (r < th)
832 {
833 m--;
834 m -= ((r > u1) | ((r == u1) & (tl > u0)));
835 }
836 }
837
838 return m;
839 }
840
841 struct gmp_div_inverse
842 {
843 /* Normalization shift count. */
844 unsigned shift;
845 /* Normalized divisor (d0 unused for mpn_div_qr_1) */
846 mp_limb_t d1, d0;
847 /* Inverse, for 2/1 or 3/2. */
848 mp_limb_t di;
849 };
850
851 static void
mpn_div_qr_1_invert(struct gmp_div_inverse * inv,mp_limb_t d)852 mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d)
853 {
854 unsigned shift;
855
856 assert (d > 0);
857 gmp_clz (shift, d);
858 inv->shift = shift;
859 inv->d1 = d << shift;
860 inv->di = mpn_invert_limb (inv->d1);
861 }
862
863 static void
mpn_div_qr_2_invert(struct gmp_div_inverse * inv,mp_limb_t d1,mp_limb_t d0)864 mpn_div_qr_2_invert (struct gmp_div_inverse *inv,
865 mp_limb_t d1, mp_limb_t d0)
866 {
867 unsigned shift;
868
869 assert (d1 > 0);
870 gmp_clz (shift, d1);
871 inv->shift = shift;
872 if (shift > 0)
873 {
874 d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
875 d0 <<= shift;
876 }
877 inv->d1 = d1;
878 inv->d0 = d0;
879 inv->di = mpn_invert_3by2 (d1, d0);
880 }
881
882 static void
mpn_div_qr_invert(struct gmp_div_inverse * inv,mp_srcptr dp,mp_size_t dn)883 mpn_div_qr_invert (struct gmp_div_inverse *inv,
884 mp_srcptr dp, mp_size_t dn)
885 {
886 assert (dn > 0);
887
888 if (dn == 1)
889 mpn_div_qr_1_invert (inv, dp[0]);
890 else if (dn == 2)
891 mpn_div_qr_2_invert (inv, dp[1], dp[0]);
892 else
893 {
894 unsigned shift;
895 mp_limb_t d1, d0;
896
897 d1 = dp[dn-1];
898 d0 = dp[dn-2];
899 assert (d1 > 0);
900 gmp_clz (shift, d1);
901 inv->shift = shift;
902 if (shift > 0)
903 {
904 d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
905 d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift));
906 }
907 inv->d1 = d1;
908 inv->d0 = d0;
909 inv->di = mpn_invert_3by2 (d1, d0);
910 }
911 }
912
913 /* Not matching current public gmp interface, rather corresponding to
914 the sbpi1_div_* functions. */
915 static mp_limb_t
mpn_div_qr_1_preinv(mp_ptr qp,mp_srcptr np,mp_size_t nn,const struct gmp_div_inverse * inv)916 mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
917 const struct gmp_div_inverse *inv)
918 {
919 mp_limb_t d, di;
920 mp_limb_t r;
921 mp_ptr tp = NULL;
922
923 if (inv->shift > 0)
924 {
925 tp = gmp_xalloc_limbs (nn);
926 r = mpn_lshift (tp, np, nn, inv->shift);
927 np = tp;
928 }
929 else
930 r = 0;
931
932 d = inv->d1;
933 di = inv->di;
934 while (--nn >= 0)
935 {
936 mp_limb_t q;
937
938 gmp_udiv_qrnnd_preinv (q, r, r, np[nn], d, di);
939 if (qp)
940 qp[nn] = q;
941 }
942 if (inv->shift > 0)
943 gmp_free (tp);
944
945 return r >> inv->shift;
946 }
947
948 static mp_limb_t
mpn_div_qr_1(mp_ptr qp,mp_srcptr np,mp_size_t nn,mp_limb_t d)949 mpn_div_qr_1 (mp_ptr qp, mp_srcptr np, mp_size_t nn, mp_limb_t d)
950 {
951 assert (d > 0);
952
953 /* Special case for powers of two. */
954 if ((d & (d-1)) == 0)
955 {
956 mp_limb_t r = np[0] & (d-1);
957 if (qp)
958 {
959 if (d <= 1)
960 mpn_copyi (qp, np, nn);
961 else
962 {
963 unsigned shift;
964 gmp_ctz (shift, d);
965 mpn_rshift (qp, np, nn, shift);
966 }
967 }
968 return r;
969 }
970 else
971 {
972 struct gmp_div_inverse inv;
973 mpn_div_qr_1_invert (&inv, d);
974 return mpn_div_qr_1_preinv (qp, np, nn, &inv);
975 }
976 }
977
978 static void
mpn_div_qr_2_preinv(mp_ptr qp,mp_ptr rp,mp_srcptr np,mp_size_t nn,const struct gmp_div_inverse * inv)979 mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
980 const struct gmp_div_inverse *inv)
981 {
982 unsigned shift;
983 mp_size_t i;
984 mp_limb_t d1, d0, di, r1, r0;
985 mp_ptr tp;
986
987 assert (nn >= 2);
988 shift = inv->shift;
989 d1 = inv->d1;
990 d0 = inv->d0;
991 di = inv->di;
992
993 if (shift > 0)
994 {
995 tp = gmp_xalloc_limbs (nn);
996 r1 = mpn_lshift (tp, np, nn, shift);
997 np = tp;
998 }
999 else
1000 r1 = 0;
1001
1002 r0 = np[nn - 1];
1003
1004 i = nn - 2;
1005 do
1006 {
1007 mp_limb_t n0, q;
1008 n0 = np[i];
1009 gmp_udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di);
1010
1011 if (qp)
1012 qp[i] = q;
1013 }
1014 while (--i >= 0);
1015
1016 if (shift > 0)
1017 {
1018 assert ((r0 << (GMP_LIMB_BITS - shift)) == 0);
1019 r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift));
1020 r1 >>= shift;
1021
1022 gmp_free (tp);
1023 }
1024
1025 rp[1] = r1;
1026 rp[0] = r0;
1027 }
1028
1029 #if 0
1030 static void
1031 mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
1032 mp_limb_t d1, mp_limb_t d0)
1033 {
1034 struct gmp_div_inverse inv;
1035 assert (nn >= 2);
1036
1037 mpn_div_qr_2_invert (&inv, d1, d0);
1038 mpn_div_qr_2_preinv (qp, rp, np, nn, &inv);
1039 }
1040 #endif
1041
1042 static void
mpn_div_qr_pi1(mp_ptr qp,mp_ptr np,mp_size_t nn,mp_limb_t n1,mp_srcptr dp,mp_size_t dn,mp_limb_t dinv)1043 mpn_div_qr_pi1 (mp_ptr qp,
1044 mp_ptr np, mp_size_t nn, mp_limb_t n1,
1045 mp_srcptr dp, mp_size_t dn,
1046 mp_limb_t dinv)
1047 {
1048 mp_size_t i;
1049
1050 mp_limb_t d1, d0;
1051 mp_limb_t cy, cy1;
1052 mp_limb_t q;
1053
1054 assert (dn > 2);
1055 assert (nn >= dn);
1056
1057 d1 = dp[dn - 1];
1058 d0 = dp[dn - 2];
1059
1060 assert ((d1 & GMP_LIMB_HIGHBIT) != 0);
1061 /* Iteration variable is the index of the q limb.
1062 *
1063 * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]>
1064 * by <d1, d0, dp[dn-3], ..., dp[0] >
1065 */
1066
1067 i = nn - dn;
1068 do
1069 {
1070 mp_limb_t n0 = np[dn-1+i];
1071
1072 if (n1 == d1 && n0 == d0)
1073 {
1074 q = GMP_LIMB_MAX;
1075 mpn_submul_1 (np+i, dp, dn, q);
1076 n1 = np[dn-1+i]; /* update n1, last loop's value will now be invalid */
1077 }
1078 else
1079 {
1080 gmp_udiv_qr_3by2 (q, n1, n0, n1, n0, np[dn-2+i], d1, d0, dinv);
1081
1082 cy = mpn_submul_1 (np + i, dp, dn-2, q);
1083
1084 cy1 = n0 < cy;
1085 n0 = n0 - cy;
1086 cy = n1 < cy1;
1087 n1 = n1 - cy1;
1088 np[dn-2+i] = n0;
1089
1090 if (cy != 0)
1091 {
1092 n1 += d1 + mpn_add_n (np + i, np + i, dp, dn - 1);
1093 q--;
1094 }
1095 }
1096
1097 if (qp)
1098 qp[i] = q;
1099 }
1100 while (--i >= 0);
1101
1102 np[dn - 1] = n1;
1103 }
1104
1105 static void
mpn_div_qr_preinv(mp_ptr qp,mp_ptr np,mp_size_t nn,mp_srcptr dp,mp_size_t dn,const struct gmp_div_inverse * inv)1106 mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
1107 mp_srcptr dp, mp_size_t dn,
1108 const struct gmp_div_inverse *inv)
1109 {
1110 assert (dn > 0);
1111 assert (nn >= dn);
1112
1113 if (dn == 1)
1114 np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv);
1115 else if (dn == 2)
1116 mpn_div_qr_2_preinv (qp, np, np, nn, inv);
1117 else
1118 {
1119 mp_limb_t nh;
1120 unsigned shift;
1121
1122 assert (inv->d1 == dp[dn-1]);
1123 assert (inv->d0 == dp[dn-2]);
1124 assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0);
1125
1126 shift = inv->shift;
1127 if (shift > 0)
1128 nh = mpn_lshift (np, np, nn, shift);
1129 else
1130 nh = 0;
1131
1132 mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);
1133
1134 if (shift > 0)
1135 gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
1136 }
1137 }
1138
1139 static void
mpn_div_qr(mp_ptr qp,mp_ptr np,mp_size_t nn,mp_srcptr dp,mp_size_t dn)1140 mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
1141 {
1142 struct gmp_div_inverse inv;
1143 mp_ptr tp = NULL;
1144
1145 assert (dn > 0);
1146 assert (nn >= dn);
1147
1148 mpn_div_qr_invert (&inv, dp, dn);
1149 if (dn > 2 && inv.shift > 0)
1150 {
1151 tp = gmp_xalloc_limbs (dn);
1152 gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift));
1153 dp = tp;
1154 }
1155 mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv);
1156 if (tp)
1157 gmp_free (tp);
1158 }
1159
1160
1161 /* MPN base conversion. */
1162 static unsigned
mpn_base_power_of_two_p(unsigned b)1163 mpn_base_power_of_two_p (unsigned b)
1164 {
1165 switch (b)
1166 {
1167 case 2: return 1;
1168 case 4: return 2;
1169 case 8: return 3;
1170 case 16: return 4;
1171 case 32: return 5;
1172 case 64: return 6;
1173 case 128: return 7;
1174 case 256: return 8;
1175 default: return 0;
1176 }
1177 }
1178
1179 struct mpn_base_info
1180 {
1181 /* bb is the largest power of the base which fits in one limb, and
1182 exp is the corresponding exponent. */
1183 unsigned exp;
1184 mp_limb_t bb;
1185 };
1186
1187 static void
mpn_get_base_info(struct mpn_base_info * info,mp_limb_t b)1188 mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b)
1189 {
1190 mp_limb_t m;
1191 mp_limb_t p;
1192 unsigned exp;
1193
1194 m = GMP_LIMB_MAX / b;
1195 for (exp = 1, p = b; p <= m; exp++)
1196 p *= b;
1197
1198 info->exp = exp;
1199 info->bb = p;
1200 }
1201
1202 static mp_bitcnt_t
mpn_limb_size_in_base_2(mp_limb_t u)1203 mpn_limb_size_in_base_2 (mp_limb_t u)
1204 {
1205 unsigned shift;
1206
1207 assert (u > 0);
1208 gmp_clz (shift, u);
1209 return GMP_LIMB_BITS - shift;
1210 }
1211
1212 static size_t
mpn_get_str_bits(unsigned char * sp,unsigned bits,mp_srcptr up,mp_size_t un)1213 mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un)
1214 {
1215 unsigned char mask;
1216 size_t sn, j;
1217 mp_size_t i;
1218 unsigned shift;
1219
1220 sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1])
1221 + bits - 1) / bits;
1222
1223 mask = (1U << bits) - 1;
1224
1225 for (i = 0, j = sn, shift = 0; j-- > 0;)
1226 {
1227 unsigned char digit = up[i] >> shift;
1228
1229 shift += bits;
1230
1231 if (shift >= GMP_LIMB_BITS && ++i < un)
1232 {
1233 shift -= GMP_LIMB_BITS;
1234 digit |= up[i] << (bits - shift);
1235 }
1236 sp[j] = digit & mask;
1237 }
1238 return sn;
1239 }
1240
1241 /* We generate digits from the least significant end, and reverse at
1242 the end. */
1243 static size_t
mpn_limb_get_str(unsigned char * sp,mp_limb_t w,const struct gmp_div_inverse * binv)1244 mpn_limb_get_str (unsigned char *sp, mp_limb_t w,
1245 const struct gmp_div_inverse *binv)
1246 {
1247 mp_size_t i;
1248 for (i = 0; w > 0; i++)
1249 {
1250 mp_limb_t h, l, r;
1251
1252 h = w >> (GMP_LIMB_BITS - binv->shift);
1253 l = w << binv->shift;
1254
1255 gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di);
1256 assert ( (r << (GMP_LIMB_BITS - binv->shift)) == 0);
1257 r >>= binv->shift;
1258
1259 sp[i] = r;
1260 }
1261 return i;
1262 }
1263
1264 static size_t
mpn_get_str_other(unsigned char * sp,int base,const struct mpn_base_info * info,mp_ptr up,mp_size_t un)1265 mpn_get_str_other (unsigned char *sp,
1266 int base, const struct mpn_base_info *info,
1267 mp_ptr up, mp_size_t un)
1268 {
1269 struct gmp_div_inverse binv;
1270 size_t sn;
1271 size_t i;
1272
1273 mpn_div_qr_1_invert (&binv, base);
1274
1275 sn = 0;
1276
1277 if (un > 1)
1278 {
1279 struct gmp_div_inverse bbinv;
1280 mpn_div_qr_1_invert (&bbinv, info->bb);
1281
1282 do
1283 {
1284 mp_limb_t w;
1285 size_t done;
1286 w = mpn_div_qr_1_preinv (up, up, un, &bbinv);
1287 un -= (up[un-1] == 0);
1288 done = mpn_limb_get_str (sp + sn, w, &binv);
1289
1290 for (sn += done; done < info->exp; done++)
1291 sp[sn++] = 0;
1292 }
1293 while (un > 1);
1294 }
1295 sn += mpn_limb_get_str (sp + sn, up[0], &binv);
1296
1297 /* Reverse order */
1298 for (i = 0; 2*i + 1 < sn; i++)
1299 {
1300 unsigned char t = sp[i];
1301 sp[i] = sp[sn - i - 1];
1302 sp[sn - i - 1] = t;
1303 }
1304
1305 return sn;
1306 }
1307
1308 size_t
mpn_get_str(unsigned char * sp,int base,mp_ptr up,mp_size_t un)1309 mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un)
1310 {
1311 unsigned bits;
1312
1313 assert (un > 0);
1314 assert (up[un-1] > 0);
1315
1316 bits = mpn_base_power_of_two_p (base);
1317 if (bits)
1318 return mpn_get_str_bits (sp, bits, up, un);
1319 else
1320 {
1321 struct mpn_base_info info;
1322
1323 mpn_get_base_info (&info, base);
1324 return mpn_get_str_other (sp, base, &info, up, un);
1325 }
1326 }
1327
1328 static mp_size_t
mpn_set_str_bits(mp_ptr rp,const unsigned char * sp,size_t sn,unsigned bits)1329 mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn,
1330 unsigned bits)
1331 {
1332 mp_size_t rn;
1333 size_t j;
1334 unsigned shift;
1335
1336 for (j = sn, rn = 0, shift = 0; j-- > 0; )
1337 {
1338 if (shift == 0)
1339 {
1340 rp[rn++] = sp[j];
1341 shift += bits;
1342 }
1343 else
1344 {
1345 rp[rn-1] |= (mp_limb_t) sp[j] << shift;
1346 shift += bits;
1347 if (shift >= GMP_LIMB_BITS)
1348 {
1349 shift -= GMP_LIMB_BITS;
1350 if (shift > 0)
1351 rp[rn++] = (mp_limb_t) sp[j] >> (bits - shift);
1352 }
1353 }
1354 }
1355 rn = mpn_normalized_size (rp, rn);
1356 return rn;
1357 }
1358
1359 /* Result is usually normalized, except for all-zero input, in which
1360 case a single zero limb is written at *RP, and 1 is returned. */
1361 static mp_size_t
mpn_set_str_other(mp_ptr rp,const unsigned char * sp,size_t sn,mp_limb_t b,const struct mpn_base_info * info)1362 mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
1363 mp_limb_t b, const struct mpn_base_info *info)
1364 {
1365 mp_size_t rn;
1366 mp_limb_t w;
1367 unsigned k;
1368 size_t j;
1369
1370 assert (sn > 0);
1371
1372 k = 1 + (sn - 1) % info->exp;
1373
1374 j = 0;
1375 w = sp[j++];
1376 while (--k != 0)
1377 w = w * b + sp[j++];
1378
1379 rp[0] = w;
1380
1381 for (rn = 1; j < sn;)
1382 {
1383 mp_limb_t cy;
1384
1385 w = sp[j++];
1386 for (k = 1; k < info->exp; k++)
1387 w = w * b + sp[j++];
1388
1389 cy = mpn_mul_1 (rp, rp, rn, info->bb);
1390 cy += mpn_add_1 (rp, rp, rn, w);
1391 if (cy > 0)
1392 rp[rn++] = cy;
1393 }
1394 assert (j == sn);
1395
1396 return rn;
1397 }
1398
1399 mp_size_t
mpn_set_str(mp_ptr rp,const unsigned char * sp,size_t sn,int base)1400 mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base)
1401 {
1402 unsigned bits;
1403
1404 if (sn == 0)
1405 return 0;
1406
1407 bits = mpn_base_power_of_two_p (base);
1408 if (bits)
1409 return mpn_set_str_bits (rp, sp, sn, bits);
1410 else
1411 {
1412 struct mpn_base_info info;
1413
1414 mpn_get_base_info (&info, base);
1415 return mpn_set_str_other (rp, sp, sn, base, &info);
1416 }
1417 }
1418
1419
1420 /* MPZ interface */
1421 void
mpz_init(mpz_t r)1422 mpz_init (mpz_t r)
1423 {
1424 static const mp_limb_t dummy_limb = 0xc1a0;
1425
1426 r->_mp_alloc = 0;
1427 r->_mp_size = 0;
1428 r->_mp_d = (mp_ptr) &dummy_limb;
1429 }
1430
1431 /* The utility of this function is a bit limited, since many functions
1432 assigns the result variable using mpz_swap. */
1433 void
mpz_init2(mpz_t r,mp_bitcnt_t bits)1434 mpz_init2 (mpz_t r, mp_bitcnt_t bits)
1435 {
1436 mp_size_t rn;
1437
1438 bits -= (bits != 0); /* Round down, except if 0 */
1439 rn = 1 + bits / GMP_LIMB_BITS;
1440
1441 r->_mp_alloc = rn;
1442 r->_mp_size = 0;
1443 r->_mp_d = gmp_xalloc_limbs (rn);
1444 }
1445
1446 void
mpz_clear(mpz_t r)1447 mpz_clear (mpz_t r)
1448 {
1449 if (r->_mp_alloc)
1450 gmp_free (r->_mp_d);
1451 }
1452
1453 static mp_ptr
mpz_realloc(mpz_t r,mp_size_t size)1454 mpz_realloc (mpz_t r, mp_size_t size)
1455 {
1456 size = GMP_MAX (size, 1);
1457
1458 if (r->_mp_alloc)
1459 r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size);
1460 else
1461 r->_mp_d = gmp_xalloc_limbs (size);
1462 r->_mp_alloc = size;
1463
1464 if (GMP_ABS (r->_mp_size) > size)
1465 r->_mp_size = 0;
1466
1467 return r->_mp_d;
1468 }
1469
1470 /* Realloc for an mpz_t WHAT if it has less than NEEDED limbs. */
1471 #define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc \
1472 ? mpz_realloc(z,n) \
1473 : (z)->_mp_d)
1474
1475 /* MPZ assignment and basic conversions. */
1476 void
mpz_set_si(mpz_t r,signed long int x)1477 mpz_set_si (mpz_t r, signed long int x)
1478 {
1479 if (x >= 0)
1480 mpz_set_ui (r, x);
1481 else /* (x < 0) */
1482 {
1483 r->_mp_size = -1;
1484 MPZ_REALLOC (r, 1)[0] = GMP_NEG_CAST (unsigned long int, x);
1485 }
1486 }
1487
1488 void
mpz_set_ui(mpz_t r,unsigned long int x)1489 mpz_set_ui (mpz_t r, unsigned long int x)
1490 {
1491 if (x > 0)
1492 {
1493 r->_mp_size = 1;
1494 MPZ_REALLOC (r, 1)[0] = x;
1495 }
1496 else
1497 r->_mp_size = 0;
1498 }
1499
1500 void
mpz_set(mpz_t r,const mpz_t x)1501 mpz_set (mpz_t r, const mpz_t x)
1502 {
1503 /* Allow the NOP r == x */
1504 if (r != x)
1505 {
1506 mp_size_t n;
1507 mp_ptr rp;
1508
1509 n = GMP_ABS (x->_mp_size);
1510 rp = MPZ_REALLOC (r, n);
1511
1512 mpn_copyi (rp, x->_mp_d, n);
1513 r->_mp_size = x->_mp_size;
1514 }
1515 }
1516
1517 void
mpz_init_set_si(mpz_t r,signed long int x)1518 mpz_init_set_si (mpz_t r, signed long int x)
1519 {
1520 mpz_init (r);
1521 mpz_set_si (r, x);
1522 }
1523
1524 void
mpz_init_set_ui(mpz_t r,unsigned long int x)1525 mpz_init_set_ui (mpz_t r, unsigned long int x)
1526 {
1527 mpz_init (r);
1528 mpz_set_ui (r, x);
1529 }
1530
1531 void
mpz_init_set(mpz_t r,const mpz_t x)1532 mpz_init_set (mpz_t r, const mpz_t x)
1533 {
1534 mpz_init (r);
1535 mpz_set (r, x);
1536 }
1537
1538 int
mpz_fits_slong_p(const mpz_t u)1539 mpz_fits_slong_p (const mpz_t u)
1540 {
1541 mp_size_t us = u->_mp_size;
1542
1543 if (us == 1)
1544 return u->_mp_d[0] < GMP_LIMB_HIGHBIT;
1545 else if (us == -1)
1546 return u->_mp_d[0] <= GMP_LIMB_HIGHBIT;
1547 else
1548 return (us == 0);
1549 }
1550
1551 int
mpz_fits_ulong_p(const mpz_t u)1552 mpz_fits_ulong_p (const mpz_t u)
1553 {
1554 mp_size_t us = u->_mp_size;
1555
1556 return (us == (us > 0));
1557 }
1558
1559 long int
mpz_get_si(const mpz_t u)1560 mpz_get_si (const mpz_t u)
1561 {
1562 if (u->_mp_size < 0)
1563 /* This expression is necessary to properly handle 0x80000000 */
1564 return -1 - (long) ((u->_mp_d[0] - 1) & ~GMP_LIMB_HIGHBIT);
1565 else
1566 return (long) (mpz_get_ui (u) & ~GMP_LIMB_HIGHBIT);
1567 }
1568
1569 unsigned long int
mpz_get_ui(const mpz_t u)1570 mpz_get_ui (const mpz_t u)
1571 {
1572 return u->_mp_size == 0 ? 0 : u->_mp_d[0];
1573 }
1574
1575 size_t
mpz_size(const mpz_t u)1576 mpz_size (const mpz_t u)
1577 {
1578 return GMP_ABS (u->_mp_size);
1579 }
1580
1581 mp_limb_t
mpz_getlimbn(const mpz_t u,mp_size_t n)1582 mpz_getlimbn (const mpz_t u, mp_size_t n)
1583 {
1584 if (n >= 0 && n < GMP_ABS (u->_mp_size))
1585 return u->_mp_d[n];
1586 else
1587 return 0;
1588 }
1589
1590 void
mpz_realloc2(mpz_t x,mp_bitcnt_t n)1591 mpz_realloc2 (mpz_t x, mp_bitcnt_t n)
1592 {
1593 mpz_realloc (x, 1 + (n - (n != 0)) / GMP_LIMB_BITS);
1594 }
1595
1596 mp_srcptr
mpz_limbs_read(mpz_srcptr x)1597 mpz_limbs_read (mpz_srcptr x)
1598 {
1599 return x->_mp_d;
1600 }
1601
1602 mp_ptr
mpz_limbs_modify(mpz_t x,mp_size_t n)1603 mpz_limbs_modify (mpz_t x, mp_size_t n)
1604 {
1605 assert (n > 0);
1606 return MPZ_REALLOC (x, n);
1607 }
1608
1609 mp_ptr
mpz_limbs_write(mpz_t x,mp_size_t n)1610 mpz_limbs_write (mpz_t x, mp_size_t n)
1611 {
1612 return mpz_limbs_modify (x, n);
1613 }
1614
1615 void
mpz_limbs_finish(mpz_t x,mp_size_t xs)1616 mpz_limbs_finish (mpz_t x, mp_size_t xs)
1617 {
1618 mp_size_t xn;
1619 xn = mpn_normalized_size (x->_mp_d, GMP_ABS (xs));
1620 x->_mp_size = xs < 0 ? -xn : xn;
1621 }
1622
1623 mpz_srcptr
mpz_roinit_n(mpz_t x,mp_srcptr xp,mp_size_t xs)1624 mpz_roinit_n (mpz_t x, mp_srcptr xp, mp_size_t xs)
1625 {
1626 x->_mp_alloc = 0;
1627 x->_mp_d = (mp_ptr) xp;
1628 mpz_limbs_finish (x, xs);
1629 return x;
1630 }
1631
1632
1633 /* Conversions and comparison to double. */
1634 void
mpz_set_d(mpz_t r,double x)1635 mpz_set_d (mpz_t r, double x)
1636 {
1637 int sign;
1638 mp_ptr rp;
1639 mp_size_t rn, i;
1640 double B;
1641 double Bi;
1642 mp_limb_t f;
1643
1644 /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
1645 zero or infinity. */
1646 if (x != x || x == x * 0.5)
1647 {
1648 r->_mp_size = 0;
1649 return;
1650 }
1651
1652 sign = x < 0.0 ;
1653 if (sign)
1654 x = - x;
1655
1656 if (x < 1.0)
1657 {
1658 r->_mp_size = 0;
1659 return;
1660 }
1661 B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1662 Bi = 1.0 / B;
1663 for (rn = 1; x >= B; rn++)
1664 x *= Bi;
1665
1666 rp = MPZ_REALLOC (r, rn);
1667
1668 f = (mp_limb_t) x;
1669 x -= f;
1670 assert (x < 1.0);
1671 i = rn-1;
1672 rp[i] = f;
1673 while (--i >= 0)
1674 {
1675 x = B * x;
1676 f = (mp_limb_t) x;
1677 x -= f;
1678 assert (x < 1.0);
1679 rp[i] = f;
1680 }
1681
1682 r->_mp_size = sign ? - rn : rn;
1683 }
1684
1685 void
mpz_init_set_d(mpz_t r,double x)1686 mpz_init_set_d (mpz_t r, double x)
1687 {
1688 mpz_init (r);
1689 mpz_set_d (r, x);
1690 }
1691
1692 double
mpz_get_d(const mpz_t u)1693 mpz_get_d (const mpz_t u)
1694 {
1695 mp_size_t un;
1696 double x;
1697 double B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1698
1699 un = GMP_ABS (u->_mp_size);
1700
1701 if (un == 0)
1702 return 0.0;
1703
1704 x = u->_mp_d[--un];
1705 while (un > 0)
1706 x = B*x + u->_mp_d[--un];
1707
1708 if (u->_mp_size < 0)
1709 x = -x;
1710
1711 return x;
1712 }
1713
1714 int
mpz_cmpabs_d(const mpz_t x,double d)1715 mpz_cmpabs_d (const mpz_t x, double d)
1716 {
1717 mp_size_t xn;
1718 double B, Bi;
1719 mp_size_t i;
1720
1721 xn = x->_mp_size;
1722 d = GMP_ABS (d);
1723
1724 if (xn != 0)
1725 {
1726 xn = GMP_ABS (xn);
1727
1728 B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1729 Bi = 1.0 / B;
1730
1731 /* Scale d so it can be compared with the top limb. */
1732 for (i = 1; i < xn; i++)
1733 d *= Bi;
1734
1735 if (d >= B)
1736 return -1;
1737
1738 /* Compare floor(d) to top limb, subtract and cancel when equal. */
1739 for (i = xn; i-- > 0;)
1740 {
1741 mp_limb_t f, xl;
1742
1743 f = (mp_limb_t) d;
1744 xl = x->_mp_d[i];
1745 if (xl > f)
1746 return 1;
1747 else if (xl < f)
1748 return -1;
1749 d = B * (d - f);
1750 }
1751 }
1752 return - (d > 0.0);
1753 }
1754
1755 int
mpz_cmp_d(const mpz_t x,double d)1756 mpz_cmp_d (const mpz_t x, double d)
1757 {
1758 if (x->_mp_size < 0)
1759 {
1760 if (d >= 0.0)
1761 return -1;
1762 else
1763 return -mpz_cmpabs_d (x, d);
1764 }
1765 else
1766 {
1767 if (d < 0.0)
1768 return 1;
1769 else
1770 return mpz_cmpabs_d (x, d);
1771 }
1772 }
1773
1774
1775 /* MPZ comparisons and the like. */
1776 int
mpz_sgn(const mpz_t u)1777 mpz_sgn (const mpz_t u)
1778 {
1779 return GMP_CMP (u->_mp_size, 0);
1780 }
1781
1782 int
mpz_cmp_si(const mpz_t u,long v)1783 mpz_cmp_si (const mpz_t u, long v)
1784 {
1785 mp_size_t usize = u->_mp_size;
1786
1787 if (usize < -1)
1788 return -1;
1789 else if (v >= 0)
1790 return mpz_cmp_ui (u, v);
1791 else if (usize >= 0)
1792 return 1;
1793 else /* usize == -1 */
1794 return GMP_CMP (GMP_NEG_CAST (mp_limb_t, v), u->_mp_d[0]);
1795 }
1796
1797 int
mpz_cmp_ui(const mpz_t u,unsigned long v)1798 mpz_cmp_ui (const mpz_t u, unsigned long v)
1799 {
1800 mp_size_t usize = u->_mp_size;
1801
1802 if (usize > 1)
1803 return 1;
1804 else if (usize < 0)
1805 return -1;
1806 else
1807 return GMP_CMP (mpz_get_ui (u), v);
1808 }
1809
1810 int
mpz_cmp(const mpz_t a,const mpz_t b)1811 mpz_cmp (const mpz_t a, const mpz_t b)
1812 {
1813 mp_size_t asize = a->_mp_size;
1814 mp_size_t bsize = b->_mp_size;
1815
1816 if (asize != bsize)
1817 return (asize < bsize) ? -1 : 1;
1818 else if (asize >= 0)
1819 return mpn_cmp (a->_mp_d, b->_mp_d, asize);
1820 else
1821 return mpn_cmp (b->_mp_d, a->_mp_d, -asize);
1822 }
1823
1824 int
mpz_cmpabs_ui(const mpz_t u,unsigned long v)1825 mpz_cmpabs_ui (const mpz_t u, unsigned long v)
1826 {
1827 if (GMP_ABS (u->_mp_size) > 1)
1828 return 1;
1829 else
1830 return GMP_CMP (mpz_get_ui (u), v);
1831 }
1832
1833 int
mpz_cmpabs(const mpz_t u,const mpz_t v)1834 mpz_cmpabs (const mpz_t u, const mpz_t v)
1835 {
1836 return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size),
1837 v->_mp_d, GMP_ABS (v->_mp_size));
1838 }
1839
1840 void
mpz_abs(mpz_t r,const mpz_t u)1841 mpz_abs (mpz_t r, const mpz_t u)
1842 {
1843 mpz_set (r, u);
1844 r->_mp_size = GMP_ABS (r->_mp_size);
1845 }
1846
1847 void
mpz_neg(mpz_t r,const mpz_t u)1848 mpz_neg (mpz_t r, const mpz_t u)
1849 {
1850 mpz_set (r, u);
1851 r->_mp_size = -r->_mp_size;
1852 }
1853
1854 void
mpz_swap(mpz_t u,mpz_t v)1855 mpz_swap (mpz_t u, mpz_t v)
1856 {
1857 MP_SIZE_T_SWAP (u->_mp_size, v->_mp_size);
1858 MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc);
1859 MP_PTR_SWAP (u->_mp_d, v->_mp_d);
1860 }
1861
1862
1863 /* MPZ addition and subtraction */
1864
1865 /* Adds to the absolute value. Returns new size, but doesn't store it. */
1866 static mp_size_t
mpz_abs_add_ui(mpz_t r,const mpz_t a,unsigned long b)1867 mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1868 {
1869 mp_size_t an;
1870 mp_ptr rp;
1871 mp_limb_t cy;
1872
1873 an = GMP_ABS (a->_mp_size);
1874 if (an == 0)
1875 {
1876 MPZ_REALLOC (r, 1)[0] = b;
1877 return b > 0;
1878 }
1879
1880 rp = MPZ_REALLOC (r, an + 1);
1881
1882 cy = mpn_add_1 (rp, a->_mp_d, an, b);
1883 rp[an] = cy;
1884 an += cy;
1885
1886 return an;
1887 }
1888
1889 /* Subtract from the absolute value. Returns new size, (or -1 on underflow),
1890 but doesn't store it. */
1891 static mp_size_t
mpz_abs_sub_ui(mpz_t r,const mpz_t a,unsigned long b)1892 mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1893 {
1894 mp_size_t an = GMP_ABS (a->_mp_size);
1895 mp_ptr rp;
1896
1897 if (an == 0)
1898 {
1899 MPZ_REALLOC (r, 1)[0] = b;
1900 return -(b > 0);
1901 }
1902 rp = MPZ_REALLOC (r, an);
1903 if (an == 1 && a->_mp_d[0] < b)
1904 {
1905 rp[0] = b - a->_mp_d[0];
1906 return -1;
1907 }
1908 else
1909 {
1910 gmp_assert_nocarry (mpn_sub_1 (rp, a->_mp_d, an, b));
1911 return mpn_normalized_size (rp, an);
1912 }
1913 }
1914
1915 void
mpz_add_ui(mpz_t r,const mpz_t a,unsigned long b)1916 mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1917 {
1918 if (a->_mp_size >= 0)
1919 r->_mp_size = mpz_abs_add_ui (r, a, b);
1920 else
1921 r->_mp_size = -mpz_abs_sub_ui (r, a, b);
1922 }
1923
1924 void
mpz_sub_ui(mpz_t r,const mpz_t a,unsigned long b)1925 mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1926 {
1927 if (a->_mp_size < 0)
1928 r->_mp_size = -mpz_abs_add_ui (r, a, b);
1929 else
1930 r->_mp_size = mpz_abs_sub_ui (r, a, b);
1931 }
1932
1933 void
mpz_ui_sub(mpz_t r,unsigned long a,const mpz_t b)1934 mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b)
1935 {
1936 if (b->_mp_size < 0)
1937 r->_mp_size = mpz_abs_add_ui (r, b, a);
1938 else
1939 r->_mp_size = -mpz_abs_sub_ui (r, b, a);
1940 }
1941
1942 static mp_size_t
mpz_abs_add(mpz_t r,const mpz_t a,const mpz_t b)1943 mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b)
1944 {
1945 mp_size_t an = GMP_ABS (a->_mp_size);
1946 mp_size_t bn = GMP_ABS (b->_mp_size);
1947 mp_ptr rp;
1948 mp_limb_t cy;
1949
1950 if (an < bn)
1951 {
1952 MPZ_SRCPTR_SWAP (a, b);
1953 MP_SIZE_T_SWAP (an, bn);
1954 }
1955
1956 rp = MPZ_REALLOC (r, an + 1);
1957 cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
1958
1959 rp[an] = cy;
1960
1961 return an + cy;
1962 }
1963
1964 static mp_size_t
mpz_abs_sub(mpz_t r,const mpz_t a,const mpz_t b)1965 mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b)
1966 {
1967 mp_size_t an = GMP_ABS (a->_mp_size);
1968 mp_size_t bn = GMP_ABS (b->_mp_size);
1969 int cmp;
1970 mp_ptr rp;
1971
1972 cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn);
1973 if (cmp > 0)
1974 {
1975 rp = MPZ_REALLOC (r, an);
1976 gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn));
1977 return mpn_normalized_size (rp, an);
1978 }
1979 else if (cmp < 0)
1980 {
1981 rp = MPZ_REALLOC (r, bn);
1982 gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an));
1983 return -mpn_normalized_size (rp, bn);
1984 }
1985 else
1986 return 0;
1987 }
1988
1989 void
mpz_add(mpz_t r,const mpz_t a,const mpz_t b)1990 mpz_add (mpz_t r, const mpz_t a, const mpz_t b)
1991 {
1992 mp_size_t rn;
1993
1994 if ( (a->_mp_size ^ b->_mp_size) >= 0)
1995 rn = mpz_abs_add (r, a, b);
1996 else
1997 rn = mpz_abs_sub (r, a, b);
1998
1999 r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
2000 }
2001
2002 void
mpz_sub(mpz_t r,const mpz_t a,const mpz_t b)2003 mpz_sub (mpz_t r, const mpz_t a, const mpz_t b)
2004 {
2005 mp_size_t rn;
2006
2007 if ( (a->_mp_size ^ b->_mp_size) >= 0)
2008 rn = mpz_abs_sub (r, a, b);
2009 else
2010 rn = mpz_abs_add (r, a, b);
2011
2012 r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
2013 }
2014
2015
2016 /* MPZ multiplication */
2017 void
mpz_mul_si(mpz_t r,const mpz_t u,long int v)2018 mpz_mul_si (mpz_t r, const mpz_t u, long int v)
2019 {
2020 if (v < 0)
2021 {
2022 mpz_mul_ui (r, u, GMP_NEG_CAST (unsigned long int, v));
2023 mpz_neg (r, r);
2024 }
2025 else
2026 mpz_mul_ui (r, u, (unsigned long int) v);
2027 }
2028
2029 void
mpz_mul_ui(mpz_t r,const mpz_t u,unsigned long int v)2030 mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2031 {
2032 mp_size_t un, us;
2033 mp_ptr tp;
2034 mp_limb_t cy;
2035
2036 us = u->_mp_size;
2037
2038 if (us == 0 || v == 0)
2039 {
2040 r->_mp_size = 0;
2041 return;
2042 }
2043
2044 un = GMP_ABS (us);
2045
2046 tp = MPZ_REALLOC (r, un + 1);
2047 cy = mpn_mul_1 (tp, u->_mp_d, un, v);
2048 tp[un] = cy;
2049
2050 un += (cy > 0);
2051 r->_mp_size = (us < 0) ? - un : un;
2052 }
2053
2054 void
mpz_mul(mpz_t r,const mpz_t u,const mpz_t v)2055 mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
2056 {
2057 int sign;
2058 mp_size_t un, vn, rn;
2059 mpz_t t;
2060 mp_ptr tp;
2061
2062 un = u->_mp_size;
2063 vn = v->_mp_size;
2064
2065 if (un == 0 || vn == 0)
2066 {
2067 r->_mp_size = 0;
2068 return;
2069 }
2070
2071 sign = (un ^ vn) < 0;
2072
2073 un = GMP_ABS (un);
2074 vn = GMP_ABS (vn);
2075
2076 mpz_init2 (t, (un + vn) * GMP_LIMB_BITS);
2077
2078 tp = t->_mp_d;
2079 if (un >= vn)
2080 mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn);
2081 else
2082 mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un);
2083
2084 rn = un + vn;
2085 rn -= tp[rn-1] == 0;
2086
2087 t->_mp_size = sign ? - rn : rn;
2088 mpz_swap (r, t);
2089 mpz_clear (t);
2090 }
2091
2092 void
mpz_mul_2exp(mpz_t r,const mpz_t u,mp_bitcnt_t bits)2093 mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
2094 {
2095 mp_size_t un, rn;
2096 mp_size_t limbs;
2097 unsigned shift;
2098 mp_ptr rp;
2099
2100 un = GMP_ABS (u->_mp_size);
2101 if (un == 0)
2102 {
2103 r->_mp_size = 0;
2104 return;
2105 }
2106
2107 limbs = bits / GMP_LIMB_BITS;
2108 shift = bits % GMP_LIMB_BITS;
2109
2110 rn = un + limbs + (shift > 0);
2111 rp = MPZ_REALLOC (r, rn);
2112 if (shift > 0)
2113 {
2114 mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift);
2115 rp[rn-1] = cy;
2116 rn -= (cy == 0);
2117 }
2118 else
2119 mpn_copyd (rp + limbs, u->_mp_d, un);
2120
2121 mpn_zero (rp, limbs);
2122
2123 r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
2124 }
2125
2126 void
mpz_addmul_ui(mpz_t r,const mpz_t u,unsigned long int v)2127 mpz_addmul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2128 {
2129 mpz_t t;
2130 mpz_init (t);
2131 mpz_mul_ui (t, u, v);
2132 mpz_add (r, r, t);
2133 mpz_clear (t);
2134 }
2135
2136 void
mpz_submul_ui(mpz_t r,const mpz_t u,unsigned long int v)2137 mpz_submul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2138 {
2139 mpz_t t;
2140 mpz_init (t);
2141 mpz_mul_ui (t, u, v);
2142 mpz_sub (r, r, t);
2143 mpz_clear (t);
2144 }
2145
2146 void
mpz_addmul(mpz_t r,const mpz_t u,const mpz_t v)2147 mpz_addmul (mpz_t r, const mpz_t u, const mpz_t v)
2148 {
2149 mpz_t t;
2150 mpz_init (t);
2151 mpz_mul (t, u, v);
2152 mpz_add (r, r, t);
2153 mpz_clear (t);
2154 }
2155
2156 void
mpz_submul(mpz_t r,const mpz_t u,const mpz_t v)2157 mpz_submul (mpz_t r, const mpz_t u, const mpz_t v)
2158 {
2159 mpz_t t;
2160 mpz_init (t);
2161 mpz_mul (t, u, v);
2162 mpz_sub (r, r, t);
2163 mpz_clear (t);
2164 }
2165
2166
2167 /* MPZ division */
2168 enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };
2169
2170 /* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */
2171 static int
mpz_div_qr(mpz_t q,mpz_t r,const mpz_t n,const mpz_t d,enum mpz_div_round_mode mode)2172 mpz_div_qr (mpz_t q, mpz_t r,
2173 const mpz_t n, const mpz_t d, enum mpz_div_round_mode mode)
2174 {
2175 mp_size_t ns, ds, nn, dn, qs;
2176 ns = n->_mp_size;
2177 ds = d->_mp_size;
2178
2179 if (ds == 0)
2180 gmp_die("mpz_div_qr: Divide by zero.");
2181
2182 if (ns == 0)
2183 {
2184 if (q)
2185 q->_mp_size = 0;
2186 if (r)
2187 r->_mp_size = 0;
2188 return 0;
2189 }
2190
2191 nn = GMP_ABS (ns);
2192 dn = GMP_ABS (ds);
2193
2194 qs = ds ^ ns;
2195
2196 if (nn < dn)
2197 {
2198 if (mode == GMP_DIV_CEIL && qs >= 0)
2199 {
2200 /* q = 1, r = n - d */
2201 if (r)
2202 mpz_sub (r, n, d);
2203 if (q)
2204 mpz_set_ui (q, 1);
2205 }
2206 else if (mode == GMP_DIV_FLOOR && qs < 0)
2207 {
2208 /* q = -1, r = n + d */
2209 if (r)
2210 mpz_add (r, n, d);
2211 if (q)
2212 mpz_set_si (q, -1);
2213 }
2214 else
2215 {
2216 /* q = 0, r = d */
2217 if (r)
2218 mpz_set (r, n);
2219 if (q)
2220 q->_mp_size = 0;
2221 }
2222 return 1;
2223 }
2224 else
2225 {
2226 mp_ptr np, qp;
2227 mp_size_t qn, rn;
2228 mpz_t tq, tr;
2229
2230 mpz_init_set (tr, n);
2231 np = tr->_mp_d;
2232
2233 qn = nn - dn + 1;
2234
2235 if (q)
2236 {
2237 mpz_init2 (tq, qn * GMP_LIMB_BITS);
2238 qp = tq->_mp_d;
2239 }
2240 else
2241 qp = NULL;
2242
2243 mpn_div_qr (qp, np, nn, d->_mp_d, dn);
2244
2245 if (qp)
2246 {
2247 qn -= (qp[qn-1] == 0);
2248
2249 tq->_mp_size = qs < 0 ? -qn : qn;
2250 }
2251 rn = mpn_normalized_size (np, dn);
2252 tr->_mp_size = ns < 0 ? - rn : rn;
2253
2254 if (mode == GMP_DIV_FLOOR && qs < 0 && rn != 0)
2255 {
2256 if (q)
2257 mpz_sub_ui (tq, tq, 1);
2258 if (r)
2259 mpz_add (tr, tr, d);
2260 }
2261 else if (mode == GMP_DIV_CEIL && qs >= 0 && rn != 0)
2262 {
2263 if (q)
2264 mpz_add_ui (tq, tq, 1);
2265 if (r)
2266 mpz_sub (tr, tr, d);
2267 }
2268
2269 if (q)
2270 {
2271 mpz_swap (tq, q);
2272 mpz_clear (tq);
2273 }
2274 if (r)
2275 mpz_swap (tr, r);
2276
2277 mpz_clear (tr);
2278
2279 return rn != 0;
2280 }
2281 }
2282
2283 void
mpz_cdiv_qr(mpz_t q,mpz_t r,const mpz_t n,const mpz_t d)2284 mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2285 {
2286 mpz_div_qr (q, r, n, d, GMP_DIV_CEIL);
2287 }
2288
2289 void
mpz_fdiv_qr(mpz_t q,mpz_t r,const mpz_t n,const mpz_t d)2290 mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2291 {
2292 mpz_div_qr (q, r, n, d, GMP_DIV_FLOOR);
2293 }
2294
2295 void
mpz_tdiv_qr(mpz_t q,mpz_t r,const mpz_t n,const mpz_t d)2296 mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2297 {
2298 mpz_div_qr (q, r, n, d, GMP_DIV_TRUNC);
2299 }
2300
2301 void
mpz_cdiv_q(mpz_t q,const mpz_t n,const mpz_t d)2302 mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2303 {
2304 mpz_div_qr (q, NULL, n, d, GMP_DIV_CEIL);
2305 }
2306
2307 void
mpz_fdiv_q(mpz_t q,const mpz_t n,const mpz_t d)2308 mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2309 {
2310 mpz_div_qr (q, NULL, n, d, GMP_DIV_FLOOR);
2311 }
2312
2313 void
mpz_tdiv_q(mpz_t q,const mpz_t n,const mpz_t d)2314 mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2315 {
2316 mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC);
2317 }
2318
2319 void
mpz_cdiv_r(mpz_t r,const mpz_t n,const mpz_t d)2320 mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2321 {
2322 mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
2323 }
2324
2325 void
mpz_fdiv_r(mpz_t r,const mpz_t n,const mpz_t d)2326 mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2327 {
2328 mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
2329 }
2330
2331 void
mpz_tdiv_r(mpz_t r,const mpz_t n,const mpz_t d)2332 mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2333 {
2334 mpz_div_qr (NULL, r, n, d, GMP_DIV_TRUNC);
2335 }
2336
2337 void
mpz_mod(mpz_t r,const mpz_t n,const mpz_t d)2338 mpz_mod (mpz_t r, const mpz_t n, const mpz_t d)
2339 {
2340 mpz_div_qr (NULL, r, n, d, d->_mp_size >= 0 ? GMP_DIV_FLOOR : GMP_DIV_CEIL);
2341 }
2342
2343 static void
mpz_div_q_2exp(mpz_t q,const mpz_t u,mp_bitcnt_t bit_index,enum mpz_div_round_mode mode)2344 mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
2345 enum mpz_div_round_mode mode)
2346 {
2347 mp_size_t un, qn;
2348 mp_size_t limb_cnt;
2349 mp_ptr qp;
2350 int adjust;
2351
2352 un = u->_mp_size;
2353 if (un == 0)
2354 {
2355 q->_mp_size = 0;
2356 return;
2357 }
2358 limb_cnt = bit_index / GMP_LIMB_BITS;
2359 qn = GMP_ABS (un) - limb_cnt;
2360 bit_index %= GMP_LIMB_BITS;
2361
2362 if (mode == ((un > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* un != 0 here. */
2363 /* Note: Below, the final indexing at limb_cnt is valid because at
2364 that point we have qn > 0. */
2365 adjust = (qn <= 0
2366 || !mpn_zero_p (u->_mp_d, limb_cnt)
2367 || (u->_mp_d[limb_cnt]
2368 & (((mp_limb_t) 1 << bit_index) - 1)));
2369 else
2370 adjust = 0;
2371
2372 if (qn <= 0)
2373 qn = 0;
2374 else
2375 {
2376 qp = MPZ_REALLOC (q, qn);
2377
2378 if (bit_index != 0)
2379 {
2380 mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index);
2381 qn -= qp[qn - 1] == 0;
2382 }
2383 else
2384 {
2385 mpn_copyi (qp, u->_mp_d + limb_cnt, qn);
2386 }
2387 }
2388
2389 q->_mp_size = qn;
2390
2391 if (adjust)
2392 mpz_add_ui (q, q, 1);
2393 if (un < 0)
2394 mpz_neg (q, q);
2395 }
2396
2397 static void
mpz_div_r_2exp(mpz_t r,const mpz_t u,mp_bitcnt_t bit_index,enum mpz_div_round_mode mode)2398 mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
2399 enum mpz_div_round_mode mode)
2400 {
2401 mp_size_t us, un, rn;
2402 mp_ptr rp;
2403 mp_limb_t mask;
2404
2405 us = u->_mp_size;
2406 if (us == 0 || bit_index == 0)
2407 {
2408 r->_mp_size = 0;
2409 return;
2410 }
2411 rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
2412 assert (rn > 0);
2413
2414 rp = MPZ_REALLOC (r, rn);
2415 un = GMP_ABS (us);
2416
2417 mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index);
2418
2419 if (rn > un)
2420 {
2421 /* Quotient (with truncation) is zero, and remainder is
2422 non-zero */
2423 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2424 {
2425 /* Have to negate and sign extend. */
2426 mp_size_t i;
2427
2428 gmp_assert_nocarry (! mpn_neg (rp, u->_mp_d, un));
2429 for (i = un; i < rn - 1; i++)
2430 rp[i] = GMP_LIMB_MAX;
2431
2432 rp[rn-1] = mask;
2433 us = -us;
2434 }
2435 else
2436 {
2437 /* Just copy */
2438 if (r != u)
2439 mpn_copyi (rp, u->_mp_d, un);
2440
2441 rn = un;
2442 }
2443 }
2444 else
2445 {
2446 if (r != u)
2447 mpn_copyi (rp, u->_mp_d, rn - 1);
2448
2449 rp[rn-1] = u->_mp_d[rn-1] & mask;
2450
2451 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2452 {
2453 /* If r != 0, compute 2^{bit_count} - r. */
2454 mpn_neg (rp, rp, rn);
2455
2456 rp[rn-1] &= mask;
2457
2458 /* us is not used for anything else, so we can modify it
2459 here to indicate flipped sign. */
2460 us = -us;
2461 }
2462 }
2463 rn = mpn_normalized_size (rp, rn);
2464 r->_mp_size = us < 0 ? -rn : rn;
2465 }
2466
2467 void
mpz_cdiv_q_2exp(mpz_t r,const mpz_t u,mp_bitcnt_t cnt)2468 mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2469 {
2470 mpz_div_q_2exp (r, u, cnt, GMP_DIV_CEIL);
2471 }
2472
2473 void
mpz_fdiv_q_2exp(mpz_t r,const mpz_t u,mp_bitcnt_t cnt)2474 mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2475 {
2476 mpz_div_q_2exp (r, u, cnt, GMP_DIV_FLOOR);
2477 }
2478
2479 void
mpz_tdiv_q_2exp(mpz_t r,const mpz_t u,mp_bitcnt_t cnt)2480 mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2481 {
2482 mpz_div_q_2exp (r, u, cnt, GMP_DIV_TRUNC);
2483 }
2484
2485 void
mpz_cdiv_r_2exp(mpz_t r,const mpz_t u,mp_bitcnt_t cnt)2486 mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2487 {
2488 mpz_div_r_2exp (r, u, cnt, GMP_DIV_CEIL);
2489 }
2490
2491 void
mpz_fdiv_r_2exp(mpz_t r,const mpz_t u,mp_bitcnt_t cnt)2492 mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2493 {
2494 mpz_div_r_2exp (r, u, cnt, GMP_DIV_FLOOR);
2495 }
2496
2497 void
mpz_tdiv_r_2exp(mpz_t r,const mpz_t u,mp_bitcnt_t cnt)2498 mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2499 {
2500 mpz_div_r_2exp (r, u, cnt, GMP_DIV_TRUNC);
2501 }
2502
2503 void
mpz_divexact(mpz_t q,const mpz_t n,const mpz_t d)2504 mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d)
2505 {
2506 gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC));
2507 }
2508
2509 int
mpz_divisible_p(const mpz_t n,const mpz_t d)2510 mpz_divisible_p (const mpz_t n, const mpz_t d)
2511 {
2512 return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2513 }
2514
2515 int
mpz_congruent_p(const mpz_t a,const mpz_t b,const mpz_t m)2516 mpz_congruent_p (const mpz_t a, const mpz_t b, const mpz_t m)
2517 {
2518 mpz_t t;
2519 int res;
2520
2521 /* a == b (mod 0) iff a == b */
2522 if (mpz_sgn (m) == 0)
2523 return (mpz_cmp (a, b) == 0);
2524
2525 mpz_init (t);
2526 mpz_sub (t, a, b);
2527 res = mpz_divisible_p (t, m);
2528 mpz_clear (t);
2529
2530 return res;
2531 }
2532
2533 static unsigned long
mpz_div_qr_ui(mpz_t q,mpz_t r,const mpz_t n,unsigned long d,enum mpz_div_round_mode mode)2534 mpz_div_qr_ui (mpz_t q, mpz_t r,
2535 const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
2536 {
2537 mp_size_t ns, qn;
2538 mp_ptr qp;
2539 mp_limb_t rl;
2540 mp_size_t rs;
2541
2542 ns = n->_mp_size;
2543 if (ns == 0)
2544 {
2545 if (q)
2546 q->_mp_size = 0;
2547 if (r)
2548 r->_mp_size = 0;
2549 return 0;
2550 }
2551
2552 qn = GMP_ABS (ns);
2553 if (q)
2554 qp = MPZ_REALLOC (q, qn);
2555 else
2556 qp = NULL;
2557
2558 rl = mpn_div_qr_1 (qp, n->_mp_d, qn, d);
2559 assert (rl < d);
2560
2561 rs = rl > 0;
2562 rs = (ns < 0) ? -rs : rs;
2563
2564 if (rl > 0 && ( (mode == GMP_DIV_FLOOR && ns < 0)
2565 || (mode == GMP_DIV_CEIL && ns >= 0)))
2566 {
2567 if (q)
2568 gmp_assert_nocarry (mpn_add_1 (qp, qp, qn, 1));
2569 rl = d - rl;
2570 rs = -rs;
2571 }
2572
2573 if (r)
2574 {
2575 MPZ_REALLOC (r, 1)[0] = rl;
2576 r->_mp_size = rs;
2577 }
2578 if (q)
2579 {
2580 qn -= (qp[qn-1] == 0);
2581 assert (qn == 0 || qp[qn-1] > 0);
2582
2583 q->_mp_size = (ns < 0) ? - qn : qn;
2584 }
2585
2586 return rl;
2587 }
2588
2589 unsigned long
mpz_cdiv_qr_ui(mpz_t q,mpz_t r,const mpz_t n,unsigned long d)2590 mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2591 {
2592 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_CEIL);
2593 }
2594
2595 unsigned long
mpz_fdiv_qr_ui(mpz_t q,mpz_t r,const mpz_t n,unsigned long d)2596 mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2597 {
2598 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_FLOOR);
2599 }
2600
2601 unsigned long
mpz_tdiv_qr_ui(mpz_t q,mpz_t r,const mpz_t n,unsigned long d)2602 mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2603 {
2604 return mpz_div_qr_ui (q, r, n, d, GMP_DIV_TRUNC);
2605 }
2606
2607 unsigned long
mpz_cdiv_q_ui(mpz_t q,const mpz_t n,unsigned long d)2608 mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2609 {
2610 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_CEIL);
2611 }
2612
2613 unsigned long
mpz_fdiv_q_ui(mpz_t q,const mpz_t n,unsigned long d)2614 mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2615 {
2616 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_FLOOR);
2617 }
2618
2619 unsigned long
mpz_tdiv_q_ui(mpz_t q,const mpz_t n,unsigned long d)2620 mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2621 {
2622 return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC);
2623 }
2624
2625 unsigned long
mpz_cdiv_r_ui(mpz_t r,const mpz_t n,unsigned long d)2626 mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2627 {
2628 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_CEIL);
2629 }
2630 unsigned long
mpz_fdiv_r_ui(mpz_t r,const mpz_t n,unsigned long d)2631 mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2632 {
2633 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2634 }
2635 unsigned long
mpz_tdiv_r_ui(mpz_t r,const mpz_t n,unsigned long d)2636 mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2637 {
2638 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_TRUNC);
2639 }
2640
2641 unsigned long
mpz_cdiv_ui(const mpz_t n,unsigned long d)2642 mpz_cdiv_ui (const mpz_t n, unsigned long d)
2643 {
2644 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_CEIL);
2645 }
2646
2647 unsigned long
mpz_fdiv_ui(const mpz_t n,unsigned long d)2648 mpz_fdiv_ui (const mpz_t n, unsigned long d)
2649 {
2650 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_FLOOR);
2651 }
2652
2653 unsigned long
mpz_tdiv_ui(const mpz_t n,unsigned long d)2654 mpz_tdiv_ui (const mpz_t n, unsigned long d)
2655 {
2656 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC);
2657 }
2658
2659 unsigned long
mpz_mod_ui(mpz_t r,const mpz_t n,unsigned long d)2660 mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long d)
2661 {
2662 return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2663 }
2664
2665 void
mpz_divexact_ui(mpz_t q,const mpz_t n,unsigned long d)2666 mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d)
2667 {
2668 gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC));
2669 }
2670
2671 int
mpz_divisible_ui_p(const mpz_t n,unsigned long d)2672 mpz_divisible_ui_p (const mpz_t n, unsigned long d)
2673 {
2674 return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2675 }
2676
2677
2678 /* GCD */
2679 static mp_limb_t
mpn_gcd_11(mp_limb_t u,mp_limb_t v)2680 mpn_gcd_11 (mp_limb_t u, mp_limb_t v)
2681 {
2682 unsigned shift;
2683
2684 assert ( (u | v) > 0);
2685
2686 if (u == 0)
2687 return v;
2688 else if (v == 0)
2689 return u;
2690
2691 gmp_ctz (shift, u | v);
2692
2693 u >>= shift;
2694 v >>= shift;
2695
2696 if ( (u & 1) == 0)
2697 MP_LIMB_T_SWAP (u, v);
2698
2699 while ( (v & 1) == 0)
2700 v >>= 1;
2701
2702 while (u != v)
2703 {
2704 if (u > v)
2705 {
2706 u -= v;
2707 do
2708 u >>= 1;
2709 while ( (u & 1) == 0);
2710 }
2711 else
2712 {
2713 v -= u;
2714 do
2715 v >>= 1;
2716 while ( (v & 1) == 0);
2717 }
2718 }
2719 return u << shift;
2720 }
2721
2722 unsigned long
mpz_gcd_ui(mpz_t g,const mpz_t u,unsigned long v)2723 mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v)
2724 {
2725 mp_size_t un;
2726
2727 if (v == 0)
2728 {
2729 if (g)
2730 mpz_abs (g, u);
2731 }
2732 else
2733 {
2734 un = GMP_ABS (u->_mp_size);
2735 if (un != 0)
2736 v = mpn_gcd_11 (mpn_div_qr_1 (NULL, u->_mp_d, un, v), v);
2737
2738 if (g)
2739 mpz_set_ui (g, v);
2740 }
2741
2742 return v;
2743 }
2744
2745 static mp_bitcnt_t
mpz_make_odd(mpz_t r)2746 mpz_make_odd (mpz_t r)
2747 {
2748 mp_bitcnt_t shift;
2749
2750 assert (r->_mp_size > 0);
2751 /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */
2752 shift = mpn_common_scan (r->_mp_d[0], 0, r->_mp_d, 0, 0);
2753 mpz_tdiv_q_2exp (r, r, shift);
2754
2755 return shift;
2756 }
2757
2758 void
mpz_gcd(mpz_t g,const mpz_t u,const mpz_t v)2759 mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
2760 {
2761 mpz_t tu, tv;
2762 mp_bitcnt_t uz, vz, gz;
2763
2764 if (u->_mp_size == 0)
2765 {
2766 mpz_abs (g, v);
2767 return;
2768 }
2769 if (v->_mp_size == 0)
2770 {
2771 mpz_abs (g, u);
2772 return;
2773 }
2774
2775 mpz_init (tu);
2776 mpz_init (tv);
2777
2778 mpz_abs (tu, u);
2779 uz = mpz_make_odd (tu);
2780 mpz_abs (tv, v);
2781 vz = mpz_make_odd (tv);
2782 gz = GMP_MIN (uz, vz);
2783
2784 if (tu->_mp_size < tv->_mp_size)
2785 mpz_swap (tu, tv);
2786
2787 mpz_tdiv_r (tu, tu, tv);
2788 if (tu->_mp_size == 0)
2789 {
2790 mpz_swap (g, tv);
2791 }
2792 else
2793 for (;;)
2794 {
2795 int c;
2796
2797 mpz_make_odd (tu);
2798 c = mpz_cmp (tu, tv);
2799 if (c == 0)
2800 {
2801 mpz_swap (g, tu);
2802 break;
2803 }
2804 if (c < 0)
2805 mpz_swap (tu, tv);
2806
2807 if (tv->_mp_size == 1)
2808 {
2809 mp_limb_t vl = tv->_mp_d[0];
2810 mp_limb_t ul = mpz_tdiv_ui (tu, vl);
2811 mpz_set_ui (g, mpn_gcd_11 (ul, vl));
2812 break;
2813 }
2814 mpz_sub (tu, tu, tv);
2815 }
2816 mpz_clear (tu);
2817 mpz_clear (tv);
2818 mpz_mul_2exp (g, g, gz);
2819 }
2820
2821 void
mpz_gcdext(mpz_t g,mpz_t s,mpz_t t,const mpz_t u,const mpz_t v)2822 mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
2823 {
2824 mpz_t tu, tv, s0, s1, t0, t1;
2825 mp_bitcnt_t uz, vz, gz;
2826 mp_bitcnt_t power;
2827
2828 if (u->_mp_size == 0)
2829 {
2830 /* g = 0 u + sgn(v) v */
2831 signed long sign = mpz_sgn (v);
2832 mpz_abs (g, v);
2833 if (s)
2834 mpz_set_ui (s, 0);
2835 if (t)
2836 mpz_set_si (t, sign);
2837 return;
2838 }
2839
2840 if (v->_mp_size == 0)
2841 {
2842 /* g = sgn(u) u + 0 v */
2843 signed long sign = mpz_sgn (u);
2844 mpz_abs (g, u);
2845 if (s)
2846 mpz_set_si (s, sign);
2847 if (t)
2848 mpz_set_ui (t, 0);
2849 return;
2850 }
2851
2852 mpz_init (tu);
2853 mpz_init (tv);
2854 mpz_init (s0);
2855 mpz_init (s1);
2856 mpz_init (t0);
2857 mpz_init (t1);
2858
2859 mpz_abs (tu, u);
2860 uz = mpz_make_odd (tu);
2861 mpz_abs (tv, v);
2862 vz = mpz_make_odd (tv);
2863 gz = GMP_MIN (uz, vz);
2864
2865 uz -= gz;
2866 vz -= gz;
2867
2868 /* Cofactors corresponding to odd gcd. gz handled later. */
2869 if (tu->_mp_size < tv->_mp_size)
2870 {
2871 mpz_swap (tu, tv);
2872 MPZ_SRCPTR_SWAP (u, v);
2873 MPZ_PTR_SWAP (s, t);
2874 MP_BITCNT_T_SWAP (uz, vz);
2875 }
2876
2877 /* Maintain
2878 *
2879 * u = t0 tu + t1 tv
2880 * v = s0 tu + s1 tv
2881 *
2882 * where u and v denote the inputs with common factors of two
2883 * eliminated, and det (s0, t0; s1, t1) = 2^p. Then
2884 *
2885 * 2^p tu = s1 u - t1 v
2886 * 2^p tv = -s0 u + t0 v
2887 */
2888
2889 /* After initial division, tu = q tv + tu', we have
2890 *
2891 * u = 2^uz (tu' + q tv)
2892 * v = 2^vz tv
2893 *
2894 * or
2895 *
2896 * t0 = 2^uz, t1 = 2^uz q
2897 * s0 = 0, s1 = 2^vz
2898 */
2899
2900 mpz_setbit (t0, uz);
2901 mpz_tdiv_qr (t1, tu, tu, tv);
2902 mpz_mul_2exp (t1, t1, uz);
2903
2904 mpz_setbit (s1, vz);
2905 power = uz + vz;
2906
2907 if (tu->_mp_size > 0)
2908 {
2909 mp_bitcnt_t shift;
2910 shift = mpz_make_odd (tu);
2911 mpz_mul_2exp (t0, t0, shift);
2912 mpz_mul_2exp (s0, s0, shift);
2913 power += shift;
2914
2915 for (;;)
2916 {
2917 int c;
2918 c = mpz_cmp (tu, tv);
2919 if (c == 0)
2920 break;
2921
2922 if (c < 0)
2923 {
2924 /* tv = tv' + tu
2925 *
2926 * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv'
2927 * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */
2928
2929 mpz_sub (tv, tv, tu);
2930 mpz_add (t0, t0, t1);
2931 mpz_add (s0, s0, s1);
2932
2933 shift = mpz_make_odd (tv);
2934 mpz_mul_2exp (t1, t1, shift);
2935 mpz_mul_2exp (s1, s1, shift);
2936 }
2937 else
2938 {
2939 mpz_sub (tu, tu, tv);
2940 mpz_add (t1, t0, t1);
2941 mpz_add (s1, s0, s1);
2942
2943 shift = mpz_make_odd (tu);
2944 mpz_mul_2exp (t0, t0, shift);
2945 mpz_mul_2exp (s0, s0, shift);
2946 }
2947 power += shift;
2948 }
2949 }
2950
2951 /* Now tv = odd part of gcd, and -s0 and t0 are corresponding
2952 cofactors. */
2953
2954 mpz_mul_2exp (tv, tv, gz);
2955 mpz_neg (s0, s0);
2956
2957 /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To
2958 adjust cofactors, we need u / g and v / g */
2959
2960 mpz_divexact (s1, v, tv);
2961 mpz_abs (s1, s1);
2962 mpz_divexact (t1, u, tv);
2963 mpz_abs (t1, t1);
2964
2965 while (power-- > 0)
2966 {
2967 /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */
2968 if (mpz_odd_p (s0) || mpz_odd_p (t0))
2969 {
2970 mpz_sub (s0, s0, s1);
2971 mpz_add (t0, t0, t1);
2972 }
2973 mpz_divexact_ui (s0, s0, 2);
2974 mpz_divexact_ui (t0, t0, 2);
2975 }
2976
2977 /* Arrange so that |s| < |u| / 2g */
2978 mpz_add (s1, s0, s1);
2979 if (mpz_cmpabs (s0, s1) > 0)
2980 {
2981 mpz_swap (s0, s1);
2982 mpz_sub (t0, t0, t1);
2983 }
2984 if (u->_mp_size < 0)
2985 mpz_neg (s0, s0);
2986 if (v->_mp_size < 0)
2987 mpz_neg (t0, t0);
2988
2989 mpz_swap (g, tv);
2990 if (s)
2991 mpz_swap (s, s0);
2992 if (t)
2993 mpz_swap (t, t0);
2994
2995 mpz_clear (tu);
2996 mpz_clear (tv);
2997 mpz_clear (s0);
2998 mpz_clear (s1);
2999 mpz_clear (t0);
3000 mpz_clear (t1);
3001 }
3002
3003 void
mpz_lcm(mpz_t r,const mpz_t u,const mpz_t v)3004 mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v)
3005 {
3006 mpz_t g;
3007
3008 if (u->_mp_size == 0 || v->_mp_size == 0)
3009 {
3010 r->_mp_size = 0;
3011 return;
3012 }
3013
3014 mpz_init (g);
3015
3016 mpz_gcd (g, u, v);
3017 mpz_divexact (g, u, g);
3018 mpz_mul (r, g, v);
3019
3020 mpz_clear (g);
3021 mpz_abs (r, r);
3022 }
3023
3024 void
mpz_lcm_ui(mpz_t r,const mpz_t u,unsigned long v)3025 mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v)
3026 {
3027 if (v == 0 || u->_mp_size == 0)
3028 {
3029 r->_mp_size = 0;
3030 return;
3031 }
3032
3033 v /= mpz_gcd_ui (NULL, u, v);
3034 mpz_mul_ui (r, u, v);
3035
3036 mpz_abs (r, r);
3037 }
3038
3039 int
mpz_invert(mpz_t r,const mpz_t u,const mpz_t m)3040 mpz_invert (mpz_t r, const mpz_t u, const mpz_t m)
3041 {
3042 mpz_t g, tr;
3043 int invertible;
3044
3045 if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0)
3046 return 0;
3047
3048 mpz_init (g);
3049 mpz_init (tr);
3050
3051 mpz_gcdext (g, tr, NULL, u, m);
3052 invertible = (mpz_cmp_ui (g, 1) == 0);
3053
3054 if (invertible)
3055 {
3056 if (tr->_mp_size < 0)
3057 {
3058 if (m->_mp_size >= 0)
3059 mpz_add (tr, tr, m);
3060 else
3061 mpz_sub (tr, tr, m);
3062 }
3063 mpz_swap (r, tr);
3064 }
3065
3066 mpz_clear (g);
3067 mpz_clear (tr);
3068 return invertible;
3069 }
3070
3071
3072 /* Higher level operations (sqrt, pow and root) */
3073
3074 void
mpz_pow_ui(mpz_t r,const mpz_t b,unsigned long e)3075 mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e)
3076 {
3077 unsigned long bit;
3078 mpz_t tr;
3079 mpz_init_set_ui (tr, 1);
3080
3081 bit = GMP_ULONG_HIGHBIT;
3082 do
3083 {
3084 mpz_mul (tr, tr, tr);
3085 if (e & bit)
3086 mpz_mul (tr, tr, b);
3087 bit >>= 1;
3088 }
3089 while (bit > 0);
3090
3091 mpz_swap (r, tr);
3092 mpz_clear (tr);
3093 }
3094
3095 void
mpz_ui_pow_ui(mpz_t r,unsigned long blimb,unsigned long e)3096 mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e)
3097 {
3098 mpz_t b;
3099 mpz_pow_ui (r, mpz_roinit_n (b, &blimb, 1), e);
3100 }
3101
3102 void
mpz_powm(mpz_t r,const mpz_t b,const mpz_t e,const mpz_t m)3103 mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
3104 {
3105 mpz_t tr;
3106 mpz_t base;
3107 mp_size_t en, mn;
3108 mp_srcptr mp;
3109 struct gmp_div_inverse minv;
3110 unsigned shift;
3111 mp_ptr tp = NULL;
3112
3113 en = GMP_ABS (e->_mp_size);
3114 mn = GMP_ABS (m->_mp_size);
3115 if (mn == 0)
3116 gmp_die ("mpz_powm: Zero modulo.");
3117
3118 if (en == 0)
3119 {
3120 mpz_set_ui (r, 1);
3121 return;
3122 }
3123
3124 mp = m->_mp_d;
3125 mpn_div_qr_invert (&minv, mp, mn);
3126 shift = minv.shift;
3127
3128 if (shift > 0)
3129 {
3130 /* To avoid shifts, we do all our reductions, except the final
3131 one, using a *normalized* m. */
3132 minv.shift = 0;
3133
3134 tp = gmp_xalloc_limbs (mn);
3135 gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift));
3136 mp = tp;
3137 }
3138
3139 mpz_init (base);
3140
3141 if (e->_mp_size < 0)
3142 {
3143 if (!mpz_invert (base, b, m))
3144 gmp_die ("mpz_powm: Negative exponent and non-invertible base.");
3145 }
3146 else
3147 {
3148 mp_size_t bn;
3149 mpz_abs (base, b);
3150
3151 bn = base->_mp_size;
3152 if (bn >= mn)
3153 {
3154 mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv);
3155 bn = mn;
3156 }
3157
3158 /* We have reduced the absolute value. Now take care of the
3159 sign. Note that we get zero represented non-canonically as
3160 m. */
3161 if (b->_mp_size < 0)
3162 {
3163 mp_ptr bp = MPZ_REALLOC (base, mn);
3164 gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn));
3165 bn = mn;
3166 }
3167 base->_mp_size = mpn_normalized_size (base->_mp_d, bn);
3168 }
3169 mpz_init_set_ui (tr, 1);
3170
3171 while (--en >= 0)
3172 {
3173 mp_limb_t w = e->_mp_d[en];
3174 mp_limb_t bit;
3175
3176 bit = GMP_LIMB_HIGHBIT;
3177 do
3178 {
3179 mpz_mul (tr, tr, tr);
3180 if (w & bit)
3181 mpz_mul (tr, tr, base);
3182 if (tr->_mp_size > mn)
3183 {
3184 mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3185 tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3186 }
3187 bit >>= 1;
3188 }
3189 while (bit > 0);
3190 }
3191
3192 /* Final reduction */
3193 if (tr->_mp_size >= mn)
3194 {
3195 minv.shift = shift;
3196 mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3197 tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3198 }
3199 if (tp)
3200 gmp_free (tp);
3201
3202 mpz_swap (r, tr);
3203 mpz_clear (tr);
3204 mpz_clear (base);
3205 }
3206
3207 void
mpz_powm_ui(mpz_t r,const mpz_t b,unsigned long elimb,const mpz_t m)3208 mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
3209 {
3210 mpz_t e;
3211 mpz_powm (r, b, mpz_roinit_n (e, &elimb, 1), m);
3212 }
3213
3214 /* x=trunc(y^(1/z)), r=y-x^z */
3215 void
mpz_rootrem(mpz_t x,mpz_t r,const mpz_t y,unsigned long z)3216 mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
3217 {
3218 int sgn;
3219 mpz_t t, u;
3220
3221 sgn = y->_mp_size < 0;
3222 if ((~z & sgn) != 0)
3223 gmp_die ("mpz_rootrem: Negative argument, with even root.");
3224 if (z == 0)
3225 gmp_die ("mpz_rootrem: Zeroth root.");
3226
3227 if (mpz_cmpabs_ui (y, 1) <= 0) {
3228 if (x)
3229 mpz_set (x, y);
3230 if (r)
3231 r->_mp_size = 0;
3232 return;
3233 }
3234
3235 mpz_init (u);
3236 mpz_init (t);
3237 mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1);
3238
3239 if (z == 2) /* simplify sqrt loop: z-1 == 1 */
3240 do {
3241 mpz_swap (u, t); /* u = x */
3242 mpz_tdiv_q (t, y, u); /* t = y/x */
3243 mpz_add (t, t, u); /* t = y/x + x */
3244 mpz_tdiv_q_2exp (t, t, 1); /* x'= (y/x + x)/2 */
3245 } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */
3246 else /* z != 2 */ {
3247 mpz_t v;
3248
3249 mpz_init (v);
3250 if (sgn)
3251 mpz_neg (t, t);
3252
3253 do {
3254 mpz_swap (u, t); /* u = x */
3255 mpz_pow_ui (t, u, z - 1); /* t = x^(z-1) */
3256 mpz_tdiv_q (t, y, t); /* t = y/x^(z-1) */
3257 mpz_mul_ui (v, u, z - 1); /* v = x*(z-1) */
3258 mpz_add (t, t, v); /* t = y/x^(z-1) + x*(z-1) */
3259 mpz_tdiv_q_ui (t, t, z); /* x'=(y/x^(z-1) + x*(z-1))/z */
3260 } while (mpz_cmpabs (t, u) < 0); /* |x'| < |x| */
3261
3262 mpz_clear (v);
3263 }
3264
3265 if (r) {
3266 mpz_pow_ui (t, u, z);
3267 mpz_sub (r, y, t);
3268 }
3269 if (x)
3270 mpz_swap (x, u);
3271 mpz_clear (u);
3272 mpz_clear (t);
3273 }
3274
3275 int
mpz_root(mpz_t x,const mpz_t y,unsigned long z)3276 mpz_root (mpz_t x, const mpz_t y, unsigned long z)
3277 {
3278 int res;
3279 mpz_t r;
3280
3281 mpz_init (r);
3282 mpz_rootrem (x, r, y, z);
3283 res = r->_mp_size == 0;
3284 mpz_clear (r);
3285
3286 return res;
3287 }
3288
3289 /* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */
3290 void
mpz_sqrtrem(mpz_t s,mpz_t r,const mpz_t u)3291 mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u)
3292 {
3293 mpz_rootrem (s, r, u, 2);
3294 }
3295
3296 void
mpz_sqrt(mpz_t s,const mpz_t u)3297 mpz_sqrt (mpz_t s, const mpz_t u)
3298 {
3299 mpz_rootrem (s, NULL, u, 2);
3300 }
3301
3302 int
mpz_perfect_square_p(const mpz_t u)3303 mpz_perfect_square_p (const mpz_t u)
3304 {
3305 if (u->_mp_size <= 0)
3306 return (u->_mp_size == 0);
3307 else
3308 return mpz_root (NULL, u, 2);
3309 }
3310
3311 int
mpn_perfect_square_p(mp_srcptr p,mp_size_t n)3312 mpn_perfect_square_p (mp_srcptr p, mp_size_t n)
3313 {
3314 mpz_t t;
3315
3316 assert (n > 0);
3317 assert (p [n-1] != 0);
3318 return mpz_root (NULL, mpz_roinit_n (t, p, n), 2);
3319 }
3320
3321 mp_size_t
mpn_sqrtrem(mp_ptr sp,mp_ptr rp,mp_srcptr p,mp_size_t n)3322 mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n)
3323 {
3324 mpz_t s, r, u;
3325 mp_size_t res;
3326
3327 assert (n > 0);
3328 assert (p [n-1] != 0);
3329
3330 mpz_init (r);
3331 mpz_init (s);
3332 mpz_rootrem (s, r, mpz_roinit_n (u, p, n), 2);
3333
3334 assert (s->_mp_size == (n+1)/2);
3335 mpn_copyd (sp, s->_mp_d, s->_mp_size);
3336 mpz_clear (s);
3337 res = r->_mp_size;
3338 if (rp)
3339 mpn_copyd (rp, r->_mp_d, res);
3340 mpz_clear (r);
3341 return res;
3342 }
3343
3344 /* Combinatorics */
3345
3346 void
mpz_fac_ui(mpz_t x,unsigned long n)3347 mpz_fac_ui (mpz_t x, unsigned long n)
3348 {
3349 mpz_set_ui (x, n + (n == 0));
3350 while (n > 2)
3351 mpz_mul_ui (x, x, --n);
3352 }
3353
3354 void
mpz_bin_uiui(mpz_t r,unsigned long n,unsigned long k)3355 mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
3356 {
3357 mpz_t t;
3358
3359 mpz_set_ui (r, k <= n);
3360
3361 if (k > (n >> 1))
3362 k = (k <= n) ? n - k : 0;
3363
3364 mpz_init (t);
3365 mpz_fac_ui (t, k);
3366
3367 for (; k > 0; k--)
3368 mpz_mul_ui (r, r, n--);
3369
3370 mpz_divexact (r, r, t);
3371 mpz_clear (t);
3372 }
3373
3374
3375 /* Primality testing */
3376 static int
gmp_millerrabin(const mpz_t n,const mpz_t nm1,mpz_t y,const mpz_t q,mp_bitcnt_t k)3377 gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y,
3378 const mpz_t q, mp_bitcnt_t k)
3379 {
3380 assert (k > 0);
3381
3382 /* Caller must initialize y to the base. */
3383 mpz_powm (y, y, q, n);
3384
3385 if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
3386 return 1;
3387
3388 while (--k > 0)
3389 {
3390 mpz_powm_ui (y, y, 2, n);
3391 if (mpz_cmp (y, nm1) == 0)
3392 return 1;
3393 /* y == 1 means that the previous y was a non-trivial square root
3394 of 1 (mod n). y == 0 means that n is a power of the base.
3395 In either case, n is not prime. */
3396 if (mpz_cmp_ui (y, 1) <= 0)
3397 return 0;
3398 }
3399 return 0;
3400 }
3401
3402 /* This product is 0xc0cfd797, and fits in 32 bits. */
3403 #define GMP_PRIME_PRODUCT \
3404 (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL)
3405
3406 /* Bit (p+1)/2 is set, for each odd prime <= 61 */
3407 #define GMP_PRIME_MASK 0xc96996dcUL
3408
3409 int
mpz_probab_prime_p(const mpz_t n,int reps)3410 mpz_probab_prime_p (const mpz_t n, int reps)
3411 {
3412 mpz_t nm1;
3413 mpz_t q;
3414 mpz_t y;
3415 mp_bitcnt_t k;
3416 int is_prime;
3417 int j;
3418
3419 /* Note that we use the absolute value of n only, for compatibility
3420 with the real GMP. */
3421 if (mpz_even_p (n))
3422 return (mpz_cmpabs_ui (n, 2) == 0) ? 2 : 0;
3423
3424 /* Above test excludes n == 0 */
3425 assert (n->_mp_size != 0);
3426
3427 if (mpz_cmpabs_ui (n, 64) < 0)
3428 return (GMP_PRIME_MASK >> (n->_mp_d[0] >> 1)) & 2;
3429
3430 if (mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT) != 1)
3431 return 0;
3432
3433 /* All prime factors are >= 31. */
3434 if (mpz_cmpabs_ui (n, 31*31) < 0)
3435 return 2;
3436
3437 /* Use Miller-Rabin, with a deterministic sequence of bases, a[j] =
3438 j^2 + j + 41 using Euler's polynomial. We potentially stop early,
3439 if a[j] >= n - 1. Since n >= 31*31, this can happen only if reps >
3440 30 (a[30] == 971 > 31*31 == 961). */
3441
3442 mpz_init (nm1);
3443 mpz_init (q);
3444 mpz_init (y);
3445
3446 /* Find q and k, where q is odd and n = 1 + 2**k * q. */
3447 nm1->_mp_size = mpz_abs_sub_ui (nm1, n, 1);
3448 k = mpz_scan1 (nm1, 0);
3449 mpz_tdiv_q_2exp (q, nm1, k);
3450
3451 for (j = 0, is_prime = 1; is_prime & (j < reps); j++)
3452 {
3453 mpz_set_ui (y, (unsigned long) j*j+j+41);
3454 if (mpz_cmp (y, nm1) >= 0)
3455 {
3456 /* Don't try any further bases. This "early" break does not affect
3457 the result for any reasonable reps value (<=5000 was tested) */
3458 assert (j >= 30);
3459 break;
3460 }
3461 is_prime = gmp_millerrabin (n, nm1, y, q, k);
3462 }
3463 mpz_clear (nm1);
3464 mpz_clear (q);
3465 mpz_clear (y);
3466
3467 return is_prime;
3468 }
3469
3470
3471 /* Logical operations and bit manipulation. */
3472
3473 /* Numbers are treated as if represented in two's complement (and
3474 infinitely sign extended). For a negative values we get the two's
3475 complement from -x = ~x + 1, where ~ is bitwise complement.
3476 Negation transforms
3477
3478 xxxx10...0
3479
3480 into
3481
3482 yyyy10...0
3483
3484 where yyyy is the bitwise complement of xxxx. So least significant
3485 bits, up to and including the first one bit, are unchanged, and
3486 the more significant bits are all complemented.
3487
3488 To change a bit from zero to one in a negative number, subtract the
3489 corresponding power of two from the absolute value. This can never
3490 underflow. To change a bit from one to zero, add the corresponding
3491 power of two, and this might overflow. E.g., if x = -001111, the
3492 two's complement is 110001. Clearing the least significant bit, we
3493 get two's complement 110000, and -010000. */
3494
3495 int
mpz_tstbit(const mpz_t d,mp_bitcnt_t bit_index)3496 mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index)
3497 {
3498 mp_size_t limb_index;
3499 unsigned shift;
3500 mp_size_t ds;
3501 mp_size_t dn;
3502 mp_limb_t w;
3503 int bit;
3504
3505 ds = d->_mp_size;
3506 dn = GMP_ABS (ds);
3507 limb_index = bit_index / GMP_LIMB_BITS;
3508 if (limb_index >= dn)
3509 return ds < 0;
3510
3511 shift = bit_index % GMP_LIMB_BITS;
3512 w = d->_mp_d[limb_index];
3513 bit = (w >> shift) & 1;
3514
3515 if (ds < 0)
3516 {
3517 /* d < 0. Check if any of the bits below is set: If so, our bit
3518 must be complemented. */
3519 if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0)
3520 return bit ^ 1;
3521 while (--limb_index >= 0)
3522 if (d->_mp_d[limb_index] > 0)
3523 return bit ^ 1;
3524 }
3525 return bit;
3526 }
3527
3528 static void
mpz_abs_add_bit(mpz_t d,mp_bitcnt_t bit_index)3529 mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index)
3530 {
3531 mp_size_t dn, limb_index;
3532 mp_limb_t bit;
3533 mp_ptr dp;
3534
3535 dn = GMP_ABS (d->_mp_size);
3536
3537 limb_index = bit_index / GMP_LIMB_BITS;
3538 bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3539
3540 if (limb_index >= dn)
3541 {
3542 mp_size_t i;
3543 /* The bit should be set outside of the end of the number.
3544 We have to increase the size of the number. */
3545 dp = MPZ_REALLOC (d, limb_index + 1);
3546
3547 dp[limb_index] = bit;
3548 for (i = dn; i < limb_index; i++)
3549 dp[i] = 0;
3550 dn = limb_index + 1;
3551 }
3552 else
3553 {
3554 mp_limb_t cy;
3555
3556 dp = d->_mp_d;
3557
3558 cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit);
3559 if (cy > 0)
3560 {
3561 dp = MPZ_REALLOC (d, dn + 1);
3562 dp[dn++] = cy;
3563 }
3564 }
3565
3566 d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3567 }
3568
3569 static void
mpz_abs_sub_bit(mpz_t d,mp_bitcnt_t bit_index)3570 mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index)
3571 {
3572 mp_size_t dn, limb_index;
3573 mp_ptr dp;
3574 mp_limb_t bit;
3575
3576 dn = GMP_ABS (d->_mp_size);
3577 dp = d->_mp_d;
3578
3579 limb_index = bit_index / GMP_LIMB_BITS;
3580 bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3581
3582 assert (limb_index < dn);
3583
3584 gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index,
3585 dn - limb_index, bit));
3586 dn = mpn_normalized_size (dp, dn);
3587 d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3588 }
3589
3590 void
mpz_setbit(mpz_t d,mp_bitcnt_t bit_index)3591 mpz_setbit (mpz_t d, mp_bitcnt_t bit_index)
3592 {
3593 if (!mpz_tstbit (d, bit_index))
3594 {
3595 if (d->_mp_size >= 0)
3596 mpz_abs_add_bit (d, bit_index);
3597 else
3598 mpz_abs_sub_bit (d, bit_index);
3599 }
3600 }
3601
3602 void
mpz_clrbit(mpz_t d,mp_bitcnt_t bit_index)3603 mpz_clrbit (mpz_t d, mp_bitcnt_t bit_index)
3604 {
3605 if (mpz_tstbit (d, bit_index))
3606 {
3607 if (d->_mp_size >= 0)
3608 mpz_abs_sub_bit (d, bit_index);
3609 else
3610 mpz_abs_add_bit (d, bit_index);
3611 }
3612 }
3613
3614 void
mpz_combit(mpz_t d,mp_bitcnt_t bit_index)3615 mpz_combit (mpz_t d, mp_bitcnt_t bit_index)
3616 {
3617 if (mpz_tstbit (d, bit_index) ^ (d->_mp_size < 0))
3618 mpz_abs_sub_bit (d, bit_index);
3619 else
3620 mpz_abs_add_bit (d, bit_index);
3621 }
3622
3623 void
mpz_com(mpz_t r,const mpz_t u)3624 mpz_com (mpz_t r, const mpz_t u)
3625 {
3626 mpz_neg (r, u);
3627 mpz_sub_ui (r, r, 1);
3628 }
3629
3630 void
mpz_and(mpz_t r,const mpz_t u,const mpz_t v)3631 mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
3632 {
3633 mp_size_t un, vn, rn, i;
3634 mp_ptr up, vp, rp;
3635
3636 mp_limb_t ux, vx, rx;
3637 mp_limb_t uc, vc, rc;
3638 mp_limb_t ul, vl, rl;
3639
3640 un = GMP_ABS (u->_mp_size);
3641 vn = GMP_ABS (v->_mp_size);
3642 if (un < vn)
3643 {
3644 MPZ_SRCPTR_SWAP (u, v);
3645 MP_SIZE_T_SWAP (un, vn);
3646 }
3647 if (vn == 0)
3648 {
3649 r->_mp_size = 0;
3650 return;
3651 }
3652
3653 uc = u->_mp_size < 0;
3654 vc = v->_mp_size < 0;
3655 rc = uc & vc;
3656
3657 ux = -uc;
3658 vx = -vc;
3659 rx = -rc;
3660
3661 /* If the smaller input is positive, higher limbs don't matter. */
3662 rn = vx ? un : vn;
3663
3664 rp = MPZ_REALLOC (r, rn + (mp_size_t) rc);
3665
3666 up = u->_mp_d;
3667 vp = v->_mp_d;
3668
3669 i = 0;
3670 do
3671 {
3672 ul = (up[i] ^ ux) + uc;
3673 uc = ul < uc;
3674
3675 vl = (vp[i] ^ vx) + vc;
3676 vc = vl < vc;
3677
3678 rl = ( (ul & vl) ^ rx) + rc;
3679 rc = rl < rc;
3680 rp[i] = rl;
3681 }
3682 while (++i < vn);
3683 assert (vc == 0);
3684
3685 for (; i < rn; i++)
3686 {
3687 ul = (up[i] ^ ux) + uc;
3688 uc = ul < uc;
3689
3690 rl = ( (ul & vx) ^ rx) + rc;
3691 rc = rl < rc;
3692 rp[i] = rl;
3693 }
3694 if (rc)
3695 rp[rn++] = rc;
3696 else
3697 rn = mpn_normalized_size (rp, rn);
3698
3699 r->_mp_size = rx ? -rn : rn;
3700 }
3701
3702 void
mpz_ior(mpz_t r,const mpz_t u,const mpz_t v)3703 mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
3704 {
3705 mp_size_t un, vn, rn, i;
3706 mp_ptr up, vp, rp;
3707
3708 mp_limb_t ux, vx, rx;
3709 mp_limb_t uc, vc, rc;
3710 mp_limb_t ul, vl, rl;
3711
3712 un = GMP_ABS (u->_mp_size);
3713 vn = GMP_ABS (v->_mp_size);
3714 if (un < vn)
3715 {
3716 MPZ_SRCPTR_SWAP (u, v);
3717 MP_SIZE_T_SWAP (un, vn);
3718 }
3719 if (vn == 0)
3720 {
3721 mpz_set (r, u);
3722 return;
3723 }
3724
3725 uc = u->_mp_size < 0;
3726 vc = v->_mp_size < 0;
3727 rc = uc | vc;
3728
3729 ux = -uc;
3730 vx = -vc;
3731 rx = -rc;
3732
3733 /* If the smaller input is negative, by sign extension higher limbs
3734 don't matter. */
3735 rn = vx ? vn : un;
3736
3737 rp = MPZ_REALLOC (r, rn + (mp_size_t) rc);
3738
3739 up = u->_mp_d;
3740 vp = v->_mp_d;
3741
3742 i = 0;
3743 do
3744 {
3745 ul = (up[i] ^ ux) + uc;
3746 uc = ul < uc;
3747
3748 vl = (vp[i] ^ vx) + vc;
3749 vc = vl < vc;
3750
3751 rl = ( (ul | vl) ^ rx) + rc;
3752 rc = rl < rc;
3753 rp[i] = rl;
3754 }
3755 while (++i < vn);
3756 assert (vc == 0);
3757
3758 for (; i < rn; i++)
3759 {
3760 ul = (up[i] ^ ux) + uc;
3761 uc = ul < uc;
3762
3763 rl = ( (ul | vx) ^ rx) + rc;
3764 rc = rl < rc;
3765 rp[i] = rl;
3766 }
3767 if (rc)
3768 rp[rn++] = rc;
3769 else
3770 rn = mpn_normalized_size (rp, rn);
3771
3772 r->_mp_size = rx ? -rn : rn;
3773 }
3774
3775 void
mpz_xor(mpz_t r,const mpz_t u,const mpz_t v)3776 mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
3777 {
3778 mp_size_t un, vn, i;
3779 mp_ptr up, vp, rp;
3780
3781 mp_limb_t ux, vx, rx;
3782 mp_limb_t uc, vc, rc;
3783 mp_limb_t ul, vl, rl;
3784
3785 un = GMP_ABS (u->_mp_size);
3786 vn = GMP_ABS (v->_mp_size);
3787 if (un < vn)
3788 {
3789 MPZ_SRCPTR_SWAP (u, v);
3790 MP_SIZE_T_SWAP (un, vn);
3791 }
3792 if (vn == 0)
3793 {
3794 mpz_set (r, u);
3795 return;
3796 }
3797
3798 uc = u->_mp_size < 0;
3799 vc = v->_mp_size < 0;
3800 rc = uc ^ vc;
3801
3802 ux = -uc;
3803 vx = -vc;
3804 rx = -rc;
3805
3806 rp = MPZ_REALLOC (r, un + (mp_size_t) rc);
3807
3808 up = u->_mp_d;
3809 vp = v->_mp_d;
3810
3811 i = 0;
3812 do
3813 {
3814 ul = (up[i] ^ ux) + uc;
3815 uc = ul < uc;
3816
3817 vl = (vp[i] ^ vx) + vc;
3818 vc = vl < vc;
3819
3820 rl = (ul ^ vl ^ rx) + rc;
3821 rc = rl < rc;
3822 rp[i] = rl;
3823 }
3824 while (++i < vn);
3825 assert (vc == 0);
3826
3827 for (; i < un; i++)
3828 {
3829 ul = (up[i] ^ ux) + uc;
3830 uc = ul < uc;
3831
3832 rl = (ul ^ ux) + rc;
3833 rc = rl < rc;
3834 rp[i] = rl;
3835 }
3836 if (rc)
3837 rp[un++] = rc;
3838 else
3839 un = mpn_normalized_size (rp, un);
3840
3841 r->_mp_size = rx ? -un : un;
3842 }
3843
3844 static unsigned
gmp_popcount_limb(mp_limb_t x)3845 gmp_popcount_limb (mp_limb_t x)
3846 {
3847 unsigned c;
3848
3849 /* Do 16 bits at a time, to avoid limb-sized constants. */
3850 for (c = 0; x > 0; x >>= 16)
3851 {
3852 unsigned w = ((x >> 1) & 0x5555) + (x & 0x5555);
3853 w = ((w >> 2) & 0x3333) + (w & 0x3333);
3854 w = ((w >> 4) & 0x0f0f) + (w & 0x0f0f);
3855 w = (w >> 8) + (w & 0x00ff);
3856 c += w;
3857 }
3858 return c;
3859 }
3860
3861 mp_bitcnt_t
mpn_popcount(mp_srcptr p,mp_size_t n)3862 mpn_popcount (mp_srcptr p, mp_size_t n)
3863 {
3864 mp_size_t i;
3865 mp_bitcnt_t c;
3866
3867 for (c = 0, i = 0; i < n; i++)
3868 c += gmp_popcount_limb (p[i]);
3869
3870 return c;
3871 }
3872
3873 mp_bitcnt_t
mpz_popcount(const mpz_t u)3874 mpz_popcount (const mpz_t u)
3875 {
3876 mp_size_t un;
3877
3878 un = u->_mp_size;
3879
3880 if (un < 0)
3881 return ~(mp_bitcnt_t) 0;
3882
3883 return mpn_popcount (u->_mp_d, un);
3884 }
3885
3886 mp_bitcnt_t
mpz_hamdist(const mpz_t u,const mpz_t v)3887 mpz_hamdist (const mpz_t u, const mpz_t v)
3888 {
3889 mp_size_t un, vn, i;
3890 mp_limb_t uc, vc, ul, vl, comp;
3891 mp_srcptr up, vp;
3892 mp_bitcnt_t c;
3893
3894 un = u->_mp_size;
3895 vn = v->_mp_size;
3896
3897 if ( (un ^ vn) < 0)
3898 return ~(mp_bitcnt_t) 0;
3899
3900 comp = - (uc = vc = (un < 0));
3901 if (uc)
3902 {
3903 assert (vn < 0);
3904 un = -un;
3905 vn = -vn;
3906 }
3907
3908 up = u->_mp_d;
3909 vp = v->_mp_d;
3910
3911 if (un < vn)
3912 MPN_SRCPTR_SWAP (up, un, vp, vn);
3913
3914 for (i = 0, c = 0; i < vn; i++)
3915 {
3916 ul = (up[i] ^ comp) + uc;
3917 uc = ul < uc;
3918
3919 vl = (vp[i] ^ comp) + vc;
3920 vc = vl < vc;
3921
3922 c += gmp_popcount_limb (ul ^ vl);
3923 }
3924 assert (vc == 0);
3925
3926 for (; i < un; i++)
3927 {
3928 ul = (up[i] ^ comp) + uc;
3929 uc = ul < uc;
3930
3931 c += gmp_popcount_limb (ul ^ comp);
3932 }
3933
3934 return c;
3935 }
3936
3937 mp_bitcnt_t
mpz_scan1(const mpz_t u,mp_bitcnt_t starting_bit)3938 mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit)
3939 {
3940 mp_ptr up;
3941 mp_size_t us, un, i;
3942 mp_limb_t limb, ux;
3943
3944 us = u->_mp_size;
3945 un = GMP_ABS (us);
3946 i = starting_bit / GMP_LIMB_BITS;
3947
3948 /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit
3949 for u<0. Notice this test picks up any u==0 too. */
3950 if (i >= un)
3951 return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit);
3952
3953 up = u->_mp_d;
3954 ux = 0;
3955 limb = up[i];
3956
3957 if (starting_bit != 0)
3958 {
3959 if (us < 0)
3960 {
3961 ux = mpn_zero_p (up, i);
3962 limb = ~ limb + ux;
3963 ux = - (mp_limb_t) (limb >= ux);
3964 }
3965
3966 /* Mask to 0 all bits before starting_bit, thus ignoring them. */
3967 limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
3968 }
3969
3970 return mpn_common_scan (limb, i, up, un, ux);
3971 }
3972
3973 mp_bitcnt_t
mpz_scan0(const mpz_t u,mp_bitcnt_t starting_bit)3974 mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit)
3975 {
3976 mp_ptr up;
3977 mp_size_t us, un, i;
3978 mp_limb_t limb, ux;
3979
3980 us = u->_mp_size;
3981 ux = - (mp_limb_t) (us >= 0);
3982 un = GMP_ABS (us);
3983 i = starting_bit / GMP_LIMB_BITS;
3984
3985 /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for
3986 u<0. Notice this test picks up all cases of u==0 too. */
3987 if (i >= un)
3988 return (ux ? starting_bit : ~(mp_bitcnt_t) 0);
3989
3990 up = u->_mp_d;
3991 limb = up[i] ^ ux;
3992
3993 if (ux == 0)
3994 limb -= mpn_zero_p (up, i); /* limb = ~(~limb + zero_p) */
3995
3996 /* Mask all bits before starting_bit, thus ignoring them. */
3997 limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
3998
3999 return mpn_common_scan (limb, i, up, un, ux);
4000 }
4001
4002
4003 /* MPZ base conversion. */
4004
4005 size_t
mpz_sizeinbase(const mpz_t u,int base)4006 mpz_sizeinbase (const mpz_t u, int base)
4007 {
4008 mp_size_t un;
4009 mp_srcptr up;
4010 mp_ptr tp;
4011 mp_bitcnt_t bits;
4012 struct gmp_div_inverse bi;
4013 size_t ndigits;
4014
4015 assert (base >= 2);
4016 assert (base <= 36);
4017
4018 un = GMP_ABS (u->_mp_size);
4019 if (un == 0)
4020 return 1;
4021
4022 up = u->_mp_d;
4023
4024 bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]);
4025 switch (base)
4026 {
4027 case 2:
4028 return bits;
4029 case 4:
4030 return (bits + 1) / 2;
4031 case 8:
4032 return (bits + 2) / 3;
4033 case 16:
4034 return (bits + 3) / 4;
4035 case 32:
4036 return (bits + 4) / 5;
4037 /* FIXME: Do something more clever for the common case of base
4038 10. */
4039 }
4040
4041 tp = gmp_xalloc_limbs (un);
4042 mpn_copyi (tp, up, un);
4043 mpn_div_qr_1_invert (&bi, base);
4044
4045 ndigits = 0;
4046 do
4047 {
4048 ndigits++;
4049 mpn_div_qr_1_preinv (tp, tp, un, &bi);
4050 un -= (tp[un-1] == 0);
4051 }
4052 while (un > 0);
4053
4054 gmp_free (tp);
4055 return ndigits;
4056 }
4057
4058 char *
mpz_get_str(char * sp,int base,const mpz_t u)4059 mpz_get_str (char *sp, int base, const mpz_t u)
4060 {
4061 unsigned bits;
4062 const char *digits;
4063 mp_size_t un;
4064 size_t i, sn;
4065
4066 if (base >= 0)
4067 {
4068 digits = "0123456789abcdefghijklmnopqrstuvwxyz";
4069 }
4070 else
4071 {
4072 base = -base;
4073 digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
4074 }
4075 if (base <= 1)
4076 base = 10;
4077 if (base > 36)
4078 return NULL;
4079
4080 sn = 1 + mpz_sizeinbase (u, base);
4081 if (!sp)
4082 sp = (char *) gmp_xalloc (1 + sn);
4083
4084 un = GMP_ABS (u->_mp_size);
4085
4086 if (un == 0)
4087 {
4088 sp[0] = '0';
4089 sp[1] = '\0';
4090 return sp;
4091 }
4092
4093 i = 0;
4094
4095 if (u->_mp_size < 0)
4096 sp[i++] = '-';
4097
4098 bits = mpn_base_power_of_two_p (base);
4099
4100 if (bits)
4101 /* Not modified in this case. */
4102 sn = i + mpn_get_str_bits ((unsigned char *) sp + i, bits, u->_mp_d, un);
4103 else
4104 {
4105 struct mpn_base_info info;
4106 mp_ptr tp;
4107
4108 mpn_get_base_info (&info, base);
4109 tp = gmp_xalloc_limbs (un);
4110 mpn_copyi (tp, u->_mp_d, un);
4111
4112 sn = i + mpn_get_str_other ((unsigned char *) sp + i, base, &info, tp, un);
4113 gmp_free (tp);
4114 }
4115
4116 for (; i < sn; i++)
4117 sp[i] = digits[(unsigned char) sp[i]];
4118
4119 sp[sn] = '\0';
4120 return sp;
4121 }
4122
4123 int
mpz_set_str(mpz_t r,const char * sp,int base)4124 mpz_set_str (mpz_t r, const char *sp, int base)
4125 {
4126 unsigned bits;
4127 mp_size_t rn, alloc;
4128 mp_ptr rp;
4129 size_t dn;
4130 int sign;
4131 unsigned char *dp;
4132
4133 assert (base == 0 || (base >= 2 && base <= 36));
4134
4135 while (isspace( (unsigned char) *sp))
4136 sp++;
4137
4138 sign = (*sp == '-');
4139 sp += sign;
4140
4141 if (base == 0)
4142 {
4143 if (sp[0] == '0')
4144 {
4145 if (sp[1] == 'x' || sp[1] == 'X')
4146 {
4147 base = 16;
4148 sp += 2;
4149 }
4150 else if (sp[1] == 'b' || sp[1] == 'B')
4151 {
4152 base = 2;
4153 sp += 2;
4154 }
4155 else
4156 base = 8;
4157 }
4158 else
4159 base = 10;
4160 }
4161
4162 if (!*sp)
4163 {
4164 r->_mp_size = 0;
4165 return -1;
4166 }
4167 dp = (unsigned char *) gmp_xalloc (strlen (sp));
4168
4169 for (dn = 0; *sp; sp++)
4170 {
4171 unsigned digit;
4172
4173 if (isspace ((unsigned char) *sp))
4174 continue;
4175 else if (*sp >= '0' && *sp <= '9')
4176 digit = *sp - '0';
4177 else if (*sp >= 'a' && *sp <= 'z')
4178 digit = *sp - 'a' + 10;
4179 else if (*sp >= 'A' && *sp <= 'Z')
4180 digit = *sp - 'A' + 10;
4181 else
4182 digit = base; /* fail */
4183
4184 if (digit >= (unsigned) base)
4185 {
4186 gmp_free (dp);
4187 r->_mp_size = 0;
4188 return -1;
4189 }
4190
4191 dp[dn++] = digit;
4192 }
4193
4194 if (!dn)
4195 {
4196 gmp_free (dp);
4197 r->_mp_size = 0;
4198 return -1;
4199 }
4200 bits = mpn_base_power_of_two_p (base);
4201
4202 if (bits > 0)
4203 {
4204 alloc = (dn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
4205 rp = MPZ_REALLOC (r, alloc);
4206 rn = mpn_set_str_bits (rp, dp, dn, bits);
4207 }
4208 else
4209 {
4210 struct mpn_base_info info;
4211 mpn_get_base_info (&info, base);
4212 alloc = (dn + info.exp - 1) / info.exp;
4213 rp = MPZ_REALLOC (r, alloc);
4214 rn = mpn_set_str_other (rp, dp, dn, base, &info);
4215 /* Normalization, needed for all-zero input. */
4216 assert (rn > 0);
4217 rn -= rp[rn-1] == 0;
4218 }
4219 assert (rn <= alloc);
4220 gmp_free (dp);
4221
4222 r->_mp_size = sign ? - rn : rn;
4223
4224 return 0;
4225 }
4226
4227 int
mpz_init_set_str(mpz_t r,const char * sp,int base)4228 mpz_init_set_str (mpz_t r, const char *sp, int base)
4229 {
4230 mpz_init (r);
4231 return mpz_set_str (r, sp, base);
4232 }
4233
4234 size_t
mpz_out_str(FILE * stream,int base,const mpz_t x)4235 mpz_out_str (FILE *stream, int base, const mpz_t x)
4236 {
4237 char *str;
4238 size_t len;
4239
4240 str = mpz_get_str (NULL, base, x);
4241 len = strlen (str);
4242 len = fwrite (str, 1, len, stream);
4243 gmp_free (str);
4244 return len;
4245 }
4246
4247
4248 static int
gmp_detect_endian(void)4249 gmp_detect_endian (void)
4250 {
4251 static const int i = 2;
4252 const unsigned char *p = (const unsigned char *) &i;
4253 return 1 - *p;
4254 }
4255
4256 /* Import and export. Does not support nails. */
4257 void
mpz_import(mpz_t r,size_t count,int order,size_t size,int endian,size_t nails,const void * src)4258 mpz_import (mpz_t r, size_t count, int order, size_t size, int endian,
4259 size_t nails, const void *src)
4260 {
4261 const unsigned char *p;
4262 ptrdiff_t word_step;
4263 mp_ptr rp;
4264 mp_size_t rn;
4265
4266 /* The current (partial) limb. */
4267 mp_limb_t limb;
4268 /* The number of bytes already copied to this limb (starting from
4269 the low end). */
4270 size_t bytes;
4271 /* The index where the limb should be stored, when completed. */
4272 mp_size_t i;
4273
4274 if (nails != 0)
4275 gmp_die ("mpz_import: Nails not supported.");
4276
4277 assert (order == 1 || order == -1);
4278 assert (endian >= -1 && endian <= 1);
4279
4280 if (endian == 0)
4281 endian = gmp_detect_endian ();
4282
4283 p = (unsigned char *) src;
4284
4285 word_step = (order != endian) ? 2 * size : 0;
4286
4287 /* Process bytes from the least significant end, so point p at the
4288 least significant word. */
4289 if (order == 1)
4290 {
4291 p += size * (count - 1);
4292 word_step = - word_step;
4293 }
4294
4295 /* And at least significant byte of that word. */
4296 if (endian == 1)
4297 p += (size - 1);
4298
4299 rn = (size * count + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t);
4300 rp = MPZ_REALLOC (r, rn);
4301
4302 for (limb = 0, bytes = 0, i = 0; count > 0; count--, p += word_step)
4303 {
4304 size_t j;
4305 for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4306 {
4307 limb |= (mp_limb_t) *p << (bytes++ * CHAR_BIT);
4308 if (bytes == sizeof(mp_limb_t))
4309 {
4310 rp[i++] = limb;
4311 bytes = 0;
4312 limb = 0;
4313 }
4314 }
4315 }
4316 assert (i + (bytes > 0) == rn);
4317 if (limb != 0)
4318 rp[i++] = limb;
4319 else
4320 i = mpn_normalized_size (rp, i);
4321
4322 r->_mp_size = i;
4323 }
4324
4325 void *
mpz_export(void * r,size_t * countp,int order,size_t size,int endian,size_t nails,const mpz_t u)4326 mpz_export (void *r, size_t *countp, int order, size_t size, int endian,
4327 size_t nails, const mpz_t u)
4328 {
4329 size_t count;
4330 mp_size_t un;
4331
4332 if (nails != 0)
4333 gmp_die ("mpz_import: Nails not supported.");
4334
4335 assert (order == 1 || order == -1);
4336 assert (endian >= -1 && endian <= 1);
4337 assert (size > 0 || u->_mp_size == 0);
4338
4339 un = u->_mp_size;
4340 count = 0;
4341 if (un != 0)
4342 {
4343 size_t k;
4344 unsigned char *p;
4345 ptrdiff_t word_step;
4346 /* The current (partial) limb. */
4347 mp_limb_t limb;
4348 /* The number of bytes left to to in this limb. */
4349 size_t bytes;
4350 /* The index where the limb was read. */
4351 mp_size_t i;
4352
4353 un = GMP_ABS (un);
4354
4355 /* Count bytes in top limb. */
4356 limb = u->_mp_d[un-1];
4357 assert (limb != 0);
4358
4359 k = 0;
4360 do {
4361 k++; limb >>= CHAR_BIT;
4362 } while (limb != 0);
4363
4364 count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size;
4365
4366 if (!r)
4367 r = gmp_xalloc (count * size);
4368
4369 if (endian == 0)
4370 endian = gmp_detect_endian ();
4371
4372 p = (unsigned char *) r;
4373
4374 word_step = (order != endian) ? 2 * size : 0;
4375
4376 /* Process bytes from the least significant end, so point p at the
4377 least significant word. */
4378 if (order == 1)
4379 {
4380 p += size * (count - 1);
4381 word_step = - word_step;
4382 }
4383
4384 /* And at least significant byte of that word. */
4385 if (endian == 1)
4386 p += (size - 1);
4387
4388 for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step)
4389 {
4390 size_t j;
4391 for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4392 {
4393 if (bytes == 0)
4394 {
4395 if (i < un)
4396 limb = u->_mp_d[i++];
4397 bytes = sizeof (mp_limb_t);
4398 }
4399 *p = limb;
4400 limb >>= CHAR_BIT;
4401 bytes--;
4402 }
4403 }
4404 assert (i == un);
4405 assert (k == count);
4406 }
4407
4408 if (countp)
4409 *countp = count;
4410
4411 return r;
4412 }
4413