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