1 /* dumbmp mini GMP compatible library.
3 Copyright 2001, 2002, 2004 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
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.
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.
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/. */
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.
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
29 ANSI function definitions can be used here, since ansi2knr is run if
30 necessary. But other ANSI-isms like "const" should be avoided.
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
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. */
45 typedef unsigned long mp_limb_t;
53 #define GMP_LIMB_BITS (sizeof (mp_limb_t) * 8 / 2)
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))
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)
67 #define ASSERT(cond) \
71 fprintf (stderr, "Assertion failure\n"); \
84 fprintf (stderr, "Out of memory (alloc %d bytes)\n", n);
93 return (mp_limb_t *) xmalloc (n * sizeof (mp_limb_t));
97 mem_copyi (char *dst, char *src, int size)
100 for (i = 0; i < size; i++)
110 for (i = 2; i < n; i++)
128 mpz_realloc (mpz_t r, int n)
134 PTR(r) = (mp_limb_t *) realloc (PTR(r), n * sizeof (mp_limb_t));
137 fprintf (stderr, "Out of memory (realloc to %d)\n", n);
143 mpn_normalize (mp_limb_t *rp, int *rnp)
146 while (rn > 0 && rp[rn-1] == 0)
152 mpn_copyi (mp_limb_t *dst, mp_limb_t *src, int n)
155 for (i = 0; i < n; i++)
160 mpn_zero (mp_limb_t *rp, int rn)
163 for (i = 0; i < rn; i++)
171 PTR(r) = xmalloc_limbs (ALLOC(r));
181 SIZ (r) = 0xbadcafeL;
182 PTR (r) = (mp_limb_t *) 0xdeadbeefL;
188 return (SIZ(a) > 0 ? 1 : SIZ(a) == 0 ? 0 : -1);
197 return (PTR(a)[0] & 1) != 0;
206 return (PTR(a)[0] & 1) == 0;
210 mpz_sizeinbase (mpz_t a, int base)
213 mp_limb_t *ap = PTR (a);
224 for (hi = ap[an - 1]; hi != 0; hi >>= 1)
226 return (an - 1) * GMP_LIMB_BITS + cnt;
230 mpz_set (mpz_t r, mpz_t a)
232 mpz_realloc (r, ABSIZ (a));
234 mpn_copyi (PTR(r), PTR(a), ABSIZ (a));
238 mpz_init_set (mpz_t r, mpz_t a)
245 mpz_set_ui (mpz_t r, unsigned long ui)
252 mpn_normalize (PTR(r), &rn);
257 mpz_init_set_ui (mpz_t r, unsigned long ui)
264 mpz_setbit (mpz_t r, unsigned long bit)
266 int limb, rn, extend;
271 abort (); /* only r>=0 */
273 limb = bit / GMP_LIMB_BITS;
274 bit %= GMP_LIMB_BITS;
276 mpz_realloc (r, limb+1);
278 extend = (limb+1) - rn;
280 mpn_zero (rp + rn, extend);
282 rp[limb] |= (mp_limb_t) 1 << bit;
283 SIZ(r) = MAX (rn, limb+1);
287 mpz_tstbit (mpz_t r, unsigned long bit)
292 abort (); /* only r>=0 */
294 limb = bit / GMP_LIMB_BITS;
298 bit %= GMP_LIMB_BITS;
299 return (PTR(r)[limb] >> bit) & 1;
303 popc_limb (mp_limb_t a)
315 mpz_popcount (mpz_t a)
324 for (i = 0; i < SIZ(a); i++)
325 ret += popc_limb (PTR(a)[i]);
330 mpz_add (mpz_t r, mpz_t a, mpz_t b)
332 int an = ABSIZ (a), bn = ABSIZ (b), rn;
333 mp_limb_t *rp, *ap, *bp;
337 if ((SIZ (a) ^ SIZ (b)) < 0)
338 abort (); /* really subtraction */
342 mpz_realloc (r, MAX (an, bn) + 1);
343 ap = PTR (a); bp = PTR (b); rp = PTR (r);
346 mp_limb_t *tp; int tn;
347 tn = an; an = bn; bn = tn;
348 tp = ap; ap = bp; bp = tp;
352 for (i = 0; i < bn; i++)
354 t = ap[i] + bp[i] + cy;
358 for (i = bn; i < an; i++)
367 mpn_normalize (rp, &rn);
372 mpz_add_ui (mpz_t r, mpz_t a, unsigned long int ui)
383 mpz_sub (mpz_t r, mpz_t a, mpz_t b)
385 int an = ABSIZ (a), bn = ABSIZ (b), rn;
386 mp_limb_t *rp, *ap, *bp;
390 if ((SIZ (a) ^ SIZ (b)) < 0)
391 abort (); /* really addition */
395 mpz_realloc (r, MAX (an, bn) + 1);
396 ap = PTR (a); bp = PTR (b); rp = PTR (r);
399 mp_limb_t *tp; int tn;
400 tn = an; an = bn; bn = tn;
401 tp = ap; ap = bp; bp = tp;
405 for (i = 0; i < bn; i++)
407 t = ap[i] - bp[i] - cy;
411 for (i = bn; i < an; i++)
423 for (i = 0; i < rn; i++)
433 mpn_normalize (rp, &rn);
438 mpz_sub_ui (mpz_t r, mpz_t a, unsigned long int ui)
449 mpz_mul (mpz_t r, mpz_t a, mpz_t b)
451 int an = ABSIZ (a), bn = ABSIZ (b), rn;
452 mp_limb_t *scratch, *tmp, *ap = PTR (a), *bp = PTR (b);
456 scratch = xmalloc_limbs (an + bn);
459 for (bi = 0; bi < bn; bi++)
462 for (ai = 0; ai < an; ai++)
466 for (bi = 0; bi < bn; bi++)
468 t = ap[ai] * bp[bi] + tmp[bi] + cy;
476 mpn_normalize (scratch, &rn);
479 SIZ (r) = (SIZ (a) ^ SIZ (b)) >= 0 ? rn : -rn;
484 mpz_mul_ui (mpz_t r, mpz_t a, unsigned long int ui)
495 mpz_mul_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
506 mpz_ui_pow_ui (mpz_t r, unsigned long b, unsigned long e)
515 for (i = 0; i < e; i++)
522 mpz_tdiv_q_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
527 mp_limb_t high_limb, low_limb;
531 lcnt = bcnt / GMP_LIMB_BITS;
532 rn = ABS (as) - lcnt;
544 cnt = bcnt % GMP_LIMB_BITS;
548 tnc = GMP_LIMB_BITS - cnt;
550 low_limb = high_limb >> cnt;
552 for (i = rn - 1; i != 0; i--)
555 *rp++ = low_limb | LO (high_limb << tnc);
556 low_limb = high_limb >> cnt;
564 mpn_copyi (rp, ap, rn);
567 SIZ (r) = as >= 0 ? rn : -rn;
572 mpz_tdiv_r_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
579 bwhole = bcnt / GMP_LIMB_BITS;
580 bcnt %= GMP_LIMB_BITS;
584 PTR(r)[rn-1] &= ((mp_limb_t) 1 << bcnt) - 1;
585 mpn_normalize (PTR(r), &rn);
586 SIZ(r) = (SIZ(r) >= 0 ? rn : -rn);
591 mpz_cmp (mpz_t a, mpz_t b)
593 mp_limb_t *ap, *bp, al, bl;
594 int as = SIZ (a), bs = SIZ (b);
599 return as > bs ? 1 : -1;
601 sign = as > 0 ? 1 : -1;
605 for (i = ABS (as) - 1; i >= 0; i--)
610 return al > bl ? sign : -sign;
616 mpz_cmp_ui (mpz_t a, unsigned long b)
620 mpz_init_set_ui (bz, b);
621 ret = mpz_cmp (a, bz);
627 mpz_tdiv_qr (mpz_t q, mpz_t r, mpz_t a, mpz_t b)
632 ASSERT (mpz_sgn(a) >= 0);
633 ASSERT (mpz_sgn(b) > 0);
635 mpz_init_set (tmpr, a);
636 mpz_init_set (tmpb, b);
639 if (mpz_cmp (tmpr, tmpb) > 0)
641 cnt = mpz_sizeinbase (tmpr, 2) - mpz_sizeinbase (tmpb, 2) + 1;
642 mpz_mul_2exp (tmpb, tmpb, cnt);
644 for ( ; cnt > 0; cnt--)
646 mpz_mul_2exp (q, q, 1);
647 mpz_tdiv_q_2exp (tmpb, tmpb, 1L);
648 if (mpz_cmp (tmpr, tmpb) >= 0)
650 mpz_sub (tmpr, tmpr, tmpb);
651 mpz_add_ui (q, q, 1L);
652 ASSERT (mpz_cmp (tmpr, tmpb) < 0);
663 mpz_tdiv_qr_ui (mpz_t q, mpz_t r, mpz_t a, unsigned long b)
666 mpz_init_set_ui (bz, b);
667 mpz_tdiv_qr (q, r, a, bz);
672 mpz_tdiv_q (mpz_t q, mpz_t a, mpz_t b)
677 mpz_tdiv_qr (q, r, a, b);
682 mpz_tdiv_q_ui (mpz_t q, mpz_t n, unsigned long d)
685 mpz_init_set_ui (dz, d);
686 mpz_tdiv_q (q, n, dz);
690 /* Set inv to the inverse of d, in the style of invert_limb, ie. for
691 udiv_qrnnd_preinv. */
693 mpz_preinv_invert (mpz_t inv, mpz_t d, int numb_bits)
699 norm = numb_bits - mpz_sizeinbase (d, 2);
701 mpz_init_set_ui (t, 1L);
702 mpz_mul_2exp (t, t, 2*numb_bits - norm);
703 mpz_tdiv_q (inv, t, d);
705 mpz_mul_2exp (t, t, numb_bits);
706 mpz_sub (inv, inv, t);
711 /* Remove leading '0' characters from the start of a string, by copying the
714 strstrip_leading_zeros (char *s)
731 mpz_get_str (char *buf, int base, mpz_t a)
733 static char tohex[] = "0123456789abcdef";
735 mp_limb_t alimb, *ap;
745 buf = xmalloc (ABSIZ (a) * (GMP_LIMB_BITS / 4) + 3);
756 bn = an * (GMP_LIMB_BITS / 4);
759 for (i = 0; i < an; i++)
762 for (j = 0; j < GMP_LIMB_BITS / 4; j++)
765 *bp = tohex [alimb & 0xF];
774 strstrip_leading_zeros (buf);
779 mpz_out_str (FILE *file, int base, mpz_t a)
786 str = mpz_get_str (0, 16, a);
791 /* Calculate r satisfying r*d == 1 mod 2^n. */
793 mpz_invert_2exp (mpz_t r, mpz_t a, unsigned long n)
798 ASSERT (mpz_odd_p (a));
800 mpz_init_set_ui (inv, 1L);
803 for (i = 1; i < n; i++)
805 mpz_mul (prod, inv, a);
806 if (mpz_tstbit (prod, i) != 0)
810 mpz_mul (prod, inv, a);
811 mpz_tdiv_r_2exp (prod, prod, n);
812 ASSERT (mpz_cmp_ui (prod, 1L) == 0);
820 /* Calculate inv satisfying r*a == 1 mod 2^n. */
822 mpz_invert_ui_2exp (mpz_t r, unsigned long a, unsigned long n)
825 mpz_init_set_ui (az, a);
826 mpz_invert_2exp (r, az, n);
832 mpz_pow_ui (mpz_t x, mpz_t y, unsigned long z)
836 mpz_init_set_ui (t, 1);
845 mpz_addmul_ui (mpz_t x, mpz_t y, unsigned long z)
850 mpz_mul_ui (t, y, z);
855 /* x=floor(y^(1/z)) */
857 mpz_root (mpz_t x, mpz_t y, unsigned long z)
863 fprintf (stderr, "mpz_root does not accept negative values\n");
866 if (mpz_cmp_ui (y, 1) <= 0)
875 mpz_pow_ui (t, u, z - 1);
876 mpz_tdiv_q (t, y, t);
877 mpz_addmul_ui (t, u, z - 1);
878 mpz_tdiv_q_ui (t, t, z);
879 if (mpz_cmp (t, u) >= 0)