1 /* dumbmp mini GMP compatible library.
2
3 Copyright 2001, 2002, 2004 Free Software Foundation, Inc.
4
5 This file is part of the GNU MP Library.
6
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of the GNU Lesser General Public License as published by
9 the Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 The GNU MP Library is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 License for more details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
19
20
21 /* The code here implements a subset (a very limited subset) of the main GMP
22 functions. It's designed for use in a few build-time calculations and
23 will be slow, but highly portable.
24
25 None of the normal GMP configure things are used, nor any of the normal
26 gmp.h or gmp-impl.h. To use this file in a program just #include
27 "dumbmp.c".
28
29 ANSI function definitions can be used here, since ansi2knr is run if
30 necessary. But other ANSI-isms like "const" should be avoided.
31
32 mp_limb_t here is an unsigned long, since that's a sensible type
33 everywhere we know of, with 8*sizeof(unsigned long) giving the number of
34 bits in the type (that not being true for instance with int or short on
35 Cray vector systems.)
36
37 Only the low half of each mp_limb_t is used, so as to make carry handling
38 and limb multiplies easy. GMP_LIMB_BITS is the number of bits used. */
39
40 #include <stdio.h>
41 #include <stdlib.h>
42 #include <string.h>
43
44
45 typedef unsigned long mp_limb_t;
46
47 typedef struct {
48 int _mp_alloc;
49 int _mp_size;
50 mp_limb_t *_mp_d;
51 } mpz_t[1];
52
53 #define GMP_LIMB_BITS (sizeof (mp_limb_t) * 8 / 2)
54
55 #define ABS(x) ((x) >= 0 ? (x) : -(x))
56 #define MIN(l,o) ((l) < (o) ? (l) : (o))
57 #define MAX(h,i) ((h) > (i) ? (h) : (i))
58
59 #define ALLOC(x) ((x)->_mp_alloc)
60 #define PTR(x) ((x)->_mp_d)
61 #define SIZ(x) ((x)->_mp_size)
62 #define ABSIZ(x) ABS (SIZ (x))
63 #define LOMASK ((1L << GMP_LIMB_BITS) - 1)
64 #define LO(x) ((x) & LOMASK)
65 #define HI(x) ((x) >> GMP_LIMB_BITS)
66
67 #define ASSERT(cond) \
68 do { \
69 if (! (cond)) \
70 { \
71 fprintf (stderr, "Assertion failure\n"); \
72 abort (); \
73 } \
74 } while (0)
75
76
77 char *
xmalloc(int n)78 xmalloc (int n)
79 {
80 char *p;
81 p = malloc (n);
82 if (p == NULL)
83 {
84 fprintf (stderr, "Out of memory (alloc %d bytes)\n", n);
85 abort ();
86 }
87 return p;
88 }
89
90 mp_limb_t *
xmalloc_limbs(int n)91 xmalloc_limbs (int n)
92 {
93 return (mp_limb_t *) xmalloc (n * sizeof (mp_limb_t));
94 }
95
96 void
mem_copyi(char * dst,char * src,int size)97 mem_copyi (char *dst, char *src, int size)
98 {
99 int i;
100 for (i = 0; i < size; i++)
101 dst[i] = src[i];
102 }
103
104 static int
isprime(unsigned long int t)105 isprime (unsigned long int t)
106 {
107 unsigned long int q, r, d;
108
109 if (t < 32)
110 return (0xa08a28acUL >> t) & 1;
111 if ((t & 1) == 0)
112 return 0;
113
114 if (t % 3 == 0)
115 return 0;
116 if (t % 5 == 0)
117 return 0;
118 if (t % 7 == 0)
119 return 0;
120
121 for (d = 11;;)
122 {
123 q = t / d;
124 r = t - q * d;
125 if (q < d)
126 return 1;
127 if (r == 0)
128 break;
129 d += 2;
130 q = t / d;
131 r = t - q * d;
132 if (q < d)
133 return 1;
134 if (r == 0)
135 break;
136 d += 4;
137 }
138 return 0;
139 }
140
141 int
log2_ceil(int n)142 log2_ceil (int n)
143 {
144 int e;
145 ASSERT (n >= 1);
146 for (e = 0; ; e++)
147 if ((1 << e) >= n)
148 break;
149 return e;
150 }
151
152 void
mpz_realloc(mpz_t r,int n)153 mpz_realloc (mpz_t r, int n)
154 {
155 if (n <= ALLOC(r))
156 return;
157
158 ALLOC(r) = n;
159 PTR(r) = (mp_limb_t *) realloc (PTR(r), n * sizeof (mp_limb_t));
160 if (PTR(r) == NULL)
161 {
162 fprintf (stderr, "Out of memory (realloc to %d)\n", n);
163 abort ();
164 }
165 }
166
167 void
mpn_normalize(mp_limb_t * rp,int * rnp)168 mpn_normalize (mp_limb_t *rp, int *rnp)
169 {
170 int rn = *rnp;
171 while (rn > 0 && rp[rn-1] == 0)
172 rn--;
173 *rnp = rn;
174 }
175
176 void
mpn_copyi(mp_limb_t * dst,mp_limb_t * src,int n)177 mpn_copyi (mp_limb_t *dst, mp_limb_t *src, int n)
178 {
179 int i;
180 for (i = 0; i < n; i++)
181 dst[i] = src[i];
182 }
183
184 void
mpn_zero(mp_limb_t * rp,int rn)185 mpn_zero (mp_limb_t *rp, int rn)
186 {
187 int i;
188 for (i = 0; i < rn; i++)
189 rp[i] = 0;
190 }
191
192 void
mpz_init(mpz_t r)193 mpz_init (mpz_t r)
194 {
195 ALLOC(r) = 1;
196 PTR(r) = xmalloc_limbs (ALLOC(r));
197 PTR(r)[0] = 0;
198 SIZ(r) = 0;
199 }
200
201 void
mpz_clear(mpz_t r)202 mpz_clear (mpz_t r)
203 {
204 free (PTR (r));
205 ALLOC(r) = -1;
206 SIZ (r) = 0xbadcafeL;
207 PTR (r) = (mp_limb_t *) 0xdeadbeefL;
208 }
209
210 int
mpz_sgn(mpz_t a)211 mpz_sgn (mpz_t a)
212 {
213 return (SIZ(a) > 0 ? 1 : SIZ(a) == 0 ? 0 : -1);
214 }
215
216 int
mpz_odd_p(mpz_t a)217 mpz_odd_p (mpz_t a)
218 {
219 if (SIZ(a) == 0)
220 return 0;
221 else
222 return (PTR(a)[0] & 1) != 0;
223 }
224
225 int
mpz_even_p(mpz_t a)226 mpz_even_p (mpz_t a)
227 {
228 if (SIZ(a) == 0)
229 return 1;
230 else
231 return (PTR(a)[0] & 1) == 0;
232 }
233
234 size_t
mpz_sizeinbase(mpz_t a,int base)235 mpz_sizeinbase (mpz_t a, int base)
236 {
237 int an = ABSIZ (a);
238 mp_limb_t *ap = PTR (a);
239 int cnt;
240 mp_limb_t hi;
241
242 if (base != 2)
243 abort ();
244
245 if (an == 0)
246 return 1;
247
248 cnt = 0;
249 for (hi = ap[an - 1]; hi != 0; hi >>= 1)
250 cnt += 1;
251 return (an - 1) * GMP_LIMB_BITS + cnt;
252 }
253
254 void
mpz_set(mpz_t r,mpz_t a)255 mpz_set (mpz_t r, mpz_t a)
256 {
257 mpz_realloc (r, ABSIZ (a));
258 SIZ(r) = SIZ(a);
259 mpn_copyi (PTR(r), PTR(a), ABSIZ (a));
260 }
261
262 void
mpz_init_set(mpz_t r,mpz_t a)263 mpz_init_set (mpz_t r, mpz_t a)
264 {
265 mpz_init (r);
266 mpz_set (r, a);
267 }
268
269 void
mpz_set_ui(mpz_t r,unsigned long ui)270 mpz_set_ui (mpz_t r, unsigned long ui)
271 {
272 int rn;
273 mpz_realloc (r, 2);
274 PTR(r)[0] = LO(ui);
275 PTR(r)[1] = HI(ui);
276 rn = 2;
277 mpn_normalize (PTR(r), &rn);
278 SIZ(r) = rn;
279 }
280
281 void
mpz_init_set_ui(mpz_t r,unsigned long ui)282 mpz_init_set_ui (mpz_t r, unsigned long ui)
283 {
284 mpz_init (r);
285 mpz_set_ui (r, ui);
286 }
287
288 void
mpz_setbit(mpz_t r,unsigned long bit)289 mpz_setbit (mpz_t r, unsigned long bit)
290 {
291 int limb, rn, extend;
292 mp_limb_t *rp;
293
294 rn = SIZ(r);
295 if (rn < 0)
296 abort (); /* only r>=0 */
297
298 limb = bit / GMP_LIMB_BITS;
299 bit %= GMP_LIMB_BITS;
300
301 mpz_realloc (r, limb+1);
302 rp = PTR(r);
303 extend = (limb+1) - rn;
304 if (extend > 0)
305 mpn_zero (rp + rn, extend);
306
307 rp[limb] |= (mp_limb_t) 1 << bit;
308 SIZ(r) = MAX (rn, limb+1);
309 }
310
311 int
mpz_tstbit(mpz_t r,unsigned long bit)312 mpz_tstbit (mpz_t r, unsigned long bit)
313 {
314 int limb;
315
316 if (SIZ(r) < 0)
317 abort (); /* only r>=0 */
318
319 limb = bit / GMP_LIMB_BITS;
320 if (SIZ(r) <= limb)
321 return 0;
322
323 bit %= GMP_LIMB_BITS;
324 return (PTR(r)[limb] >> bit) & 1;
325 }
326
327 int
popc_limb(mp_limb_t a)328 popc_limb (mp_limb_t a)
329 {
330 int ret = 0;
331 while (a != 0)
332 {
333 ret += (a & 1);
334 a >>= 1;
335 }
336 return ret;
337 }
338
339 unsigned long
mpz_popcount(mpz_t a)340 mpz_popcount (mpz_t a)
341 {
342 unsigned long ret;
343 int i;
344
345 if (SIZ(a) < 0)
346 abort ();
347
348 ret = 0;
349 for (i = 0; i < SIZ(a); i++)
350 ret += popc_limb (PTR(a)[i]);
351 return ret;
352 }
353
354 void
mpz_add(mpz_t r,mpz_t a,mpz_t b)355 mpz_add (mpz_t r, mpz_t a, mpz_t b)
356 {
357 int an = ABSIZ (a), bn = ABSIZ (b), rn;
358 mp_limb_t *rp, *ap, *bp;
359 int i;
360 mp_limb_t t, cy;
361
362 if ((SIZ (a) ^ SIZ (b)) < 0)
363 abort (); /* really subtraction */
364 if (SIZ (a) < 0)
365 abort ();
366
367 mpz_realloc (r, MAX (an, bn) + 1);
368 ap = PTR (a); bp = PTR (b); rp = PTR (r);
369 if (an < bn)
370 {
371 mp_limb_t *tp; int tn;
372 tn = an; an = bn; bn = tn;
373 tp = ap; ap = bp; bp = tp;
374 }
375
376 cy = 0;
377 for (i = 0; i < bn; i++)
378 {
379 t = ap[i] + bp[i] + cy;
380 rp[i] = LO (t);
381 cy = HI (t);
382 }
383 for (i = bn; i < an; i++)
384 {
385 t = ap[i] + cy;
386 rp[i] = LO (t);
387 cy = HI (t);
388 }
389 rp[an] = cy;
390 rn = an + 1;
391
392 mpn_normalize (rp, &rn);
393 SIZ (r) = rn;
394 }
395
396 void
mpz_add_ui(mpz_t r,mpz_t a,unsigned long int ui)397 mpz_add_ui (mpz_t r, mpz_t a, unsigned long int ui)
398 {
399 mpz_t b;
400
401 mpz_init (b);
402 mpz_set_ui (b, ui);
403 mpz_add (r, a, b);
404 mpz_clear (b);
405 }
406
407 void
mpz_sub(mpz_t r,mpz_t a,mpz_t b)408 mpz_sub (mpz_t r, mpz_t a, mpz_t b)
409 {
410 int an = ABSIZ (a), bn = ABSIZ (b), rn;
411 mp_limb_t *rp, *ap, *bp;
412 int i;
413 mp_limb_t t, cy;
414
415 if ((SIZ (a) ^ SIZ (b)) < 0)
416 abort (); /* really addition */
417 if (SIZ (a) < 0)
418 abort ();
419
420 mpz_realloc (r, MAX (an, bn) + 1);
421 ap = PTR (a); bp = PTR (b); rp = PTR (r);
422 if (an < bn)
423 {
424 mp_limb_t *tp; int tn;
425 tn = an; an = bn; bn = tn;
426 tp = ap; ap = bp; bp = tp;
427 }
428
429 cy = 0;
430 for (i = 0; i < bn; i++)
431 {
432 t = ap[i] - bp[i] - cy;
433 rp[i] = LO (t);
434 cy = LO (-HI (t));
435 }
436 for (i = bn; i < an; i++)
437 {
438 t = ap[i] - cy;
439 rp[i] = LO (t);
440 cy = LO (-HI (t));
441 }
442 rp[an] = cy;
443 rn = an + 1;
444
445 if (cy != 0)
446 {
447 cy = 0;
448 for (i = 0; i < rn; i++)
449 {
450 t = -rp[i] - cy;
451 rp[i] = LO (t);
452 cy = LO (-HI (t));
453 }
454 SIZ (r) = -rn;
455 return;
456 }
457
458 mpn_normalize (rp, &rn);
459 SIZ (r) = rn;
460 }
461
462 void
mpz_sub_ui(mpz_t r,mpz_t a,unsigned long int ui)463 mpz_sub_ui (mpz_t r, mpz_t a, unsigned long int ui)
464 {
465 mpz_t b;
466
467 mpz_init (b);
468 mpz_set_ui (b, ui);
469 mpz_sub (r, a, b);
470 mpz_clear (b);
471 }
472
473 void
mpz_mul(mpz_t r,mpz_t a,mpz_t b)474 mpz_mul (mpz_t r, mpz_t a, mpz_t b)
475 {
476 int an = ABSIZ (a), bn = ABSIZ (b), rn;
477 mp_limb_t *scratch, *tmp, *ap = PTR (a), *bp = PTR (b);
478 int ai, bi;
479 mp_limb_t t, cy;
480
481 scratch = xmalloc_limbs (an + bn);
482 tmp = scratch;
483
484 for (bi = 0; bi < bn; bi++)
485 tmp[bi] = 0;
486
487 for (ai = 0; ai < an; ai++)
488 {
489 tmp = scratch + ai;
490 cy = 0;
491 for (bi = 0; bi < bn; bi++)
492 {
493 t = ap[ai] * bp[bi] + tmp[bi] + cy;
494 tmp[bi] = LO (t);
495 cy = HI (t);
496 }
497 tmp[bn] = cy;
498 }
499
500 rn = an + bn;
501 mpn_normalize (scratch, &rn);
502 free (PTR (r));
503 PTR (r) = scratch;
504 SIZ (r) = (SIZ (a) ^ SIZ (b)) >= 0 ? rn : -rn;
505 ALLOC (r) = an + bn;
506 }
507
508 void
mpz_mul_ui(mpz_t r,mpz_t a,unsigned long int ui)509 mpz_mul_ui (mpz_t r, mpz_t a, unsigned long int ui)
510 {
511 mpz_t b;
512
513 mpz_init (b);
514 mpz_set_ui (b, ui);
515 mpz_mul (r, a, b);
516 mpz_clear (b);
517 }
518
519 void
mpz_mul_2exp(mpz_t r,mpz_t a,unsigned long int bcnt)520 mpz_mul_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
521 {
522 mpz_set (r, a);
523 while (bcnt)
524 {
525 mpz_add (r, r, r);
526 bcnt -= 1;
527 }
528 }
529
530 void
mpz_ui_pow_ui(mpz_t r,unsigned long b,unsigned long e)531 mpz_ui_pow_ui (mpz_t r, unsigned long b, unsigned long e)
532 {
533 unsigned long i;
534 mpz_t bz;
535
536 mpz_init (bz);
537 mpz_set_ui (bz, b);
538
539 mpz_set_ui (r, 1L);
540 for (i = 0; i < e; i++)
541 mpz_mul (r, r, bz);
542
543 mpz_clear (bz);
544 }
545
546 void
mpz_tdiv_q_2exp(mpz_t r,mpz_t a,unsigned long int bcnt)547 mpz_tdiv_q_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
548 {
549 int as, rn;
550 int cnt, tnc;
551 int lcnt;
552 mp_limb_t high_limb, low_limb;
553 int i;
554
555 as = SIZ (a);
556 lcnt = bcnt / GMP_LIMB_BITS;
557 rn = ABS (as) - lcnt;
558 if (rn <= 0)
559 SIZ (r) = 0;
560 else
561 {
562 mp_limb_t *rp, *ap;
563
564 mpz_realloc (r, rn);
565
566 rp = PTR (r);
567 ap = PTR (a);
568
569 cnt = bcnt % GMP_LIMB_BITS;
570 if (cnt != 0)
571 {
572 ap += lcnt;
573 tnc = GMP_LIMB_BITS - cnt;
574 high_limb = *ap++;
575 low_limb = high_limb >> cnt;
576
577 for (i = rn - 1; i != 0; i--)
578 {
579 high_limb = *ap++;
580 *rp++ = low_limb | LO (high_limb << tnc);
581 low_limb = high_limb >> cnt;
582 }
583 *rp = low_limb;
584 rn -= low_limb == 0;
585 }
586 else
587 {
588 ap += lcnt;
589 mpn_copyi (rp, ap, rn);
590 }
591
592 SIZ (r) = as >= 0 ? rn : -rn;
593 }
594 }
595
596 void
mpz_tdiv_r_2exp(mpz_t r,mpz_t a,unsigned long int bcnt)597 mpz_tdiv_r_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
598 {
599 int rn, bwhole;
600
601 mpz_set (r, a);
602 rn = ABSIZ(r);
603
604 bwhole = bcnt / GMP_LIMB_BITS;
605 bcnt %= GMP_LIMB_BITS;
606 if (rn > bwhole)
607 {
608 rn = bwhole+1;
609 PTR(r)[rn-1] &= ((mp_limb_t) 1 << bcnt) - 1;
610 mpn_normalize (PTR(r), &rn);
611 SIZ(r) = (SIZ(r) >= 0 ? rn : -rn);
612 }
613 }
614
615 int
mpz_cmp(mpz_t a,mpz_t b)616 mpz_cmp (mpz_t a, mpz_t b)
617 {
618 mp_limb_t *ap, *bp, al, bl;
619 int as = SIZ (a), bs = SIZ (b);
620 int i;
621 int sign;
622
623 if (as != bs)
624 return as > bs ? 1 : -1;
625
626 sign = as > 0 ? 1 : -1;
627
628 ap = PTR (a);
629 bp = PTR (b);
630 for (i = ABS (as) - 1; i >= 0; i--)
631 {
632 al = ap[i];
633 bl = bp[i];
634 if (al != bl)
635 return al > bl ? sign : -sign;
636 }
637 return 0;
638 }
639
640 int
mpz_cmp_ui(mpz_t a,unsigned long b)641 mpz_cmp_ui (mpz_t a, unsigned long b)
642 {
643 mpz_t bz;
644 int ret;
645 mpz_init_set_ui (bz, b);
646 ret = mpz_cmp (a, bz);
647 mpz_clear (bz);
648 return ret;
649 }
650
651 void
mpz_tdiv_qr(mpz_t q,mpz_t r,mpz_t a,mpz_t b)652 mpz_tdiv_qr (mpz_t q, mpz_t r, mpz_t a, mpz_t b)
653 {
654 mpz_t tmpr, tmpb;
655 unsigned long cnt;
656
657 ASSERT (mpz_sgn(a) >= 0);
658 ASSERT (mpz_sgn(b) > 0);
659
660 mpz_init_set (tmpr, a);
661 mpz_init_set (tmpb, b);
662 mpz_set_ui (q, 0L);
663
664 if (mpz_cmp (tmpr, tmpb) > 0)
665 {
666 cnt = mpz_sizeinbase (tmpr, 2) - mpz_sizeinbase (tmpb, 2) + 1;
667 mpz_mul_2exp (tmpb, tmpb, cnt);
668
669 for ( ; cnt > 0; cnt--)
670 {
671 mpz_mul_2exp (q, q, 1);
672 mpz_tdiv_q_2exp (tmpb, tmpb, 1L);
673 if (mpz_cmp (tmpr, tmpb) >= 0)
674 {
675 mpz_sub (tmpr, tmpr, tmpb);
676 mpz_add_ui (q, q, 1L);
677 ASSERT (mpz_cmp (tmpr, tmpb) < 0);
678 }
679 }
680 }
681
682 mpz_set (r, tmpr);
683 mpz_clear (tmpr);
684 mpz_clear (tmpb);
685 }
686
687 void
mpz_tdiv_qr_ui(mpz_t q,mpz_t r,mpz_t a,unsigned long b)688 mpz_tdiv_qr_ui (mpz_t q, mpz_t r, mpz_t a, unsigned long b)
689 {
690 mpz_t bz;
691 mpz_init_set_ui (bz, b);
692 mpz_tdiv_qr (q, r, a, bz);
693 mpz_clear (bz);
694 }
695
696 void
mpz_tdiv_q(mpz_t q,mpz_t a,mpz_t b)697 mpz_tdiv_q (mpz_t q, mpz_t a, mpz_t b)
698 {
699 mpz_t r;
700
701 mpz_init (r);
702 mpz_tdiv_qr (q, r, a, b);
703 mpz_clear (r);
704 }
705
706 void
mpz_tdiv_r(mpz_t r,mpz_t a,mpz_t b)707 mpz_tdiv_r (mpz_t r, mpz_t a, mpz_t b)
708 {
709 mpz_t q;
710
711 mpz_init (q);
712 mpz_tdiv_qr (q, r, a, b);
713 mpz_clear (q);
714 }
715
716 void
mpz_tdiv_q_ui(mpz_t q,mpz_t n,unsigned long d)717 mpz_tdiv_q_ui (mpz_t q, mpz_t n, unsigned long d)
718 {
719 mpz_t dz;
720 mpz_init_set_ui (dz, d);
721 mpz_tdiv_q (q, n, dz);
722 mpz_clear (dz);
723 }
724
725 /* Set inv to the inverse of d, in the style of invert_limb, ie. for
726 udiv_qrnnd_preinv. */
727 void
mpz_preinv_invert(mpz_t inv,mpz_t d,int numb_bits)728 mpz_preinv_invert (mpz_t inv, mpz_t d, int numb_bits)
729 {
730 mpz_t t;
731 int norm;
732 ASSERT (SIZ(d) > 0);
733
734 norm = numb_bits - mpz_sizeinbase (d, 2);
735 ASSERT (norm >= 0);
736 mpz_init_set_ui (t, 1L);
737 mpz_mul_2exp (t, t, 2*numb_bits - norm);
738 mpz_tdiv_q (inv, t, d);
739 mpz_set_ui (t, 1L);
740 mpz_mul_2exp (t, t, numb_bits);
741 mpz_sub (inv, inv, t);
742
743 mpz_clear (t);
744 }
745
746 /* Remove leading '0' characters from the start of a string, by copying the
747 remainder down. */
748 void
strstrip_leading_zeros(char * s)749 strstrip_leading_zeros (char *s)
750 {
751 char c, *p;
752
753 p = s;
754 while (*s == '0')
755 s++;
756
757 do
758 {
759 c = *s++;
760 *p++ = c;
761 }
762 while (c != '\0');
763 }
764
765 char *
mpz_get_str(char * buf,int base,mpz_t a)766 mpz_get_str (char *buf, int base, mpz_t a)
767 {
768 static char tohex[] = "0123456789abcdef";
769
770 mp_limb_t alimb, *ap;
771 int an, bn, i, j;
772 char *bp;
773
774 if (base != 16)
775 abort ();
776 if (SIZ (a) < 0)
777 abort ();
778
779 if (buf == 0)
780 buf = xmalloc (ABSIZ (a) * (GMP_LIMB_BITS / 4) + 3);
781
782 an = ABSIZ (a);
783 if (an == 0)
784 {
785 buf[0] = '0';
786 buf[1] = '\0';
787 return buf;
788 }
789
790 ap = PTR (a);
791 bn = an * (GMP_LIMB_BITS / 4);
792 bp = buf + bn;
793
794 for (i = 0; i < an; i++)
795 {
796 alimb = ap[i];
797 for (j = 0; j < GMP_LIMB_BITS / 4; j++)
798 {
799 bp--;
800 *bp = tohex [alimb & 0xF];
801 alimb >>= 4;
802 }
803 ASSERT (alimb == 0);
804 }
805 ASSERT (bp == buf);
806
807 buf[bn] = '\0';
808
809 strstrip_leading_zeros (buf);
810 return buf;
811 }
812
813 void
mpz_out_str(FILE * file,int base,mpz_t a)814 mpz_out_str (FILE *file, int base, mpz_t a)
815 {
816 char *str;
817
818 if (file == 0)
819 file = stdout;
820
821 str = mpz_get_str (0, 16, a);
822 fputs (str, file);
823 free (str);
824 }
825
826 /* Calculate r satisfying r*d == 1 mod 2^n. */
827 void
mpz_invert_2exp(mpz_t r,mpz_t a,unsigned long n)828 mpz_invert_2exp (mpz_t r, mpz_t a, unsigned long n)
829 {
830 unsigned long i;
831 mpz_t inv, prod;
832
833 ASSERT (mpz_odd_p (a));
834
835 mpz_init_set_ui (inv, 1L);
836 mpz_init (prod);
837
838 for (i = 1; i < n; i++)
839 {
840 mpz_mul (prod, inv, a);
841 if (mpz_tstbit (prod, i) != 0)
842 mpz_setbit (inv, i);
843 }
844
845 mpz_mul (prod, inv, a);
846 mpz_tdiv_r_2exp (prod, prod, n);
847 ASSERT (mpz_cmp_ui (prod, 1L) == 0);
848
849 mpz_set (r, inv);
850
851 mpz_clear (inv);
852 mpz_clear (prod);
853 }
854
855 /* Calculate inv satisfying r*a == 1 mod 2^n. */
856 void
mpz_invert_ui_2exp(mpz_t r,unsigned long a,unsigned long n)857 mpz_invert_ui_2exp (mpz_t r, unsigned long a, unsigned long n)
858 {
859 mpz_t az;
860 mpz_init_set_ui (az, a);
861 mpz_invert_2exp (r, az, n);
862 mpz_clear (az);
863 }
864
865 /* x=y^z */
866 void
mpz_pow_ui(mpz_t x,mpz_t y,unsigned long z)867 mpz_pow_ui (mpz_t x, mpz_t y, unsigned long z)
868 {
869 mpz_t t;
870
871 mpz_init_set_ui (t, 1);
872 for (; z != 0; z--)
873 mpz_mul (t, t, y);
874 mpz_set (x, t);
875 mpz_clear (t);
876 }
877
878 /* x=x+y*z */
879 void
mpz_addmul_ui(mpz_t x,mpz_t y,unsigned long z)880 mpz_addmul_ui (mpz_t x, mpz_t y, unsigned long z)
881 {
882 mpz_t t;
883
884 mpz_init (t);
885 mpz_mul_ui (t, y, z);
886 mpz_add (x, x, t);
887 mpz_clear (t);
888 }
889
890 /* x=floor(y^(1/z)) */
891 void
mpz_root(mpz_t x,mpz_t y,unsigned long z)892 mpz_root (mpz_t x, mpz_t y, unsigned long z)
893 {
894 mpz_t t, u;
895
896 if (mpz_sgn (y) < 0)
897 {
898 fprintf (stderr, "mpz_root does not accept negative values\n");
899 abort ();
900 }
901 if (mpz_cmp_ui (y, 1) <= 0)
902 {
903 mpz_set (x, y);
904 return;
905 }
906 mpz_init (t);
907 mpz_init_set (u, y);
908 do
909 {
910 mpz_pow_ui (t, u, z - 1);
911 mpz_tdiv_q (t, y, t);
912 mpz_addmul_ui (t, u, z - 1);
913 mpz_tdiv_q_ui (t, t, z);
914 if (mpz_cmp (t, u) >= 0)
915 break;
916 mpz_set (u, t);
917 }
918 while (1);
919 mpz_set (x, u);
920 mpz_clear (t);
921 mpz_clear (u);
922 }
923