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++)
105 isprime (unsigned long int t)
107 unsigned long int q, r, d;
110 return (0xa08a28acUL >> t) & 1;
153 mpz_realloc (mpz_t r, int n)
159 PTR(r) = (mp_limb_t *) realloc (PTR(r), n * sizeof (mp_limb_t));
162 fprintf (stderr, "Out of memory (realloc to %d)\n", n);
168 mpn_normalize (mp_limb_t *rp, int *rnp)
171 while (rn > 0 && rp[rn-1] == 0)
177 mpn_copyi (mp_limb_t *dst, mp_limb_t *src, int n)
180 for (i = 0; i < n; i++)
185 mpn_zero (mp_limb_t *rp, int rn)
188 for (i = 0; i < rn; i++)
196 PTR(r) = xmalloc_limbs (ALLOC(r));
206 SIZ (r) = 0xbadcafeL;
207 PTR (r) = (mp_limb_t *) 0xdeadbeefL;
213 return (SIZ(a) > 0 ? 1 : SIZ(a) == 0 ? 0 : -1);
222 return (PTR(a)[0] & 1) != 0;
231 return (PTR(a)[0] & 1) == 0;
235 mpz_sizeinbase (mpz_t a, int base)
238 mp_limb_t *ap = PTR (a);
249 for (hi = ap[an - 1]; hi != 0; hi >>= 1)
251 return (an - 1) * GMP_LIMB_BITS + cnt;
255 mpz_set (mpz_t r, mpz_t a)
257 mpz_realloc (r, ABSIZ (a));
259 mpn_copyi (PTR(r), PTR(a), ABSIZ (a));
263 mpz_init_set (mpz_t r, mpz_t a)
270 mpz_set_ui (mpz_t r, unsigned long ui)
277 mpn_normalize (PTR(r), &rn);
282 mpz_init_set_ui (mpz_t r, unsigned long ui)
289 mpz_setbit (mpz_t r, unsigned long bit)
291 int limb, rn, extend;
296 abort (); /* only r>=0 */
298 limb = bit / GMP_LIMB_BITS;
299 bit %= GMP_LIMB_BITS;
301 mpz_realloc (r, limb+1);
303 extend = (limb+1) - rn;
305 mpn_zero (rp + rn, extend);
307 rp[limb] |= (mp_limb_t) 1 << bit;
308 SIZ(r) = MAX (rn, limb+1);
312 mpz_tstbit (mpz_t r, unsigned long bit)
317 abort (); /* only r>=0 */
319 limb = bit / GMP_LIMB_BITS;
323 bit %= GMP_LIMB_BITS;
324 return (PTR(r)[limb] >> bit) & 1;
328 popc_limb (mp_limb_t a)
340 mpz_popcount (mpz_t a)
349 for (i = 0; i < SIZ(a); i++)
350 ret += popc_limb (PTR(a)[i]);
355 mpz_add (mpz_t r, mpz_t a, mpz_t b)
357 int an = ABSIZ (a), bn = ABSIZ (b), rn;
358 mp_limb_t *rp, *ap, *bp;
362 if ((SIZ (a) ^ SIZ (b)) < 0)
363 abort (); /* really subtraction */
367 mpz_realloc (r, MAX (an, bn) + 1);
368 ap = PTR (a); bp = PTR (b); rp = PTR (r);
371 mp_limb_t *tp; int tn;
372 tn = an; an = bn; bn = tn;
373 tp = ap; ap = bp; bp = tp;
377 for (i = 0; i < bn; i++)
379 t = ap[i] + bp[i] + cy;
383 for (i = bn; i < an; i++)
392 mpn_normalize (rp, &rn);
397 mpz_add_ui (mpz_t r, mpz_t a, unsigned long int ui)
408 mpz_sub (mpz_t r, mpz_t a, mpz_t b)
410 int an = ABSIZ (a), bn = ABSIZ (b), rn;
411 mp_limb_t *rp, *ap, *bp;
415 if ((SIZ (a) ^ SIZ (b)) < 0)
416 abort (); /* really addition */
420 mpz_realloc (r, MAX (an, bn) + 1);
421 ap = PTR (a); bp = PTR (b); rp = PTR (r);
424 mp_limb_t *tp; int tn;
425 tn = an; an = bn; bn = tn;
426 tp = ap; ap = bp; bp = tp;
430 for (i = 0; i < bn; i++)
432 t = ap[i] - bp[i] - cy;
436 for (i = bn; i < an; i++)
448 for (i = 0; i < rn; i++)
458 mpn_normalize (rp, &rn);
463 mpz_sub_ui (mpz_t r, mpz_t a, unsigned long int ui)
474 mpz_mul (mpz_t r, mpz_t a, mpz_t b)
476 int an = ABSIZ (a), bn = ABSIZ (b), rn;
477 mp_limb_t *scratch, *tmp, *ap = PTR (a), *bp = PTR (b);
481 scratch = xmalloc_limbs (an + bn);
484 for (bi = 0; bi < bn; bi++)
487 for (ai = 0; ai < an; ai++)
491 for (bi = 0; bi < bn; bi++)
493 t = ap[ai] * bp[bi] + tmp[bi] + cy;
501 mpn_normalize (scratch, &rn);
504 SIZ (r) = (SIZ (a) ^ SIZ (b)) >= 0 ? rn : -rn;
509 mpz_mul_ui (mpz_t r, mpz_t a, unsigned long int ui)
520 mpz_mul_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
531 mpz_ui_pow_ui (mpz_t r, unsigned long b, unsigned long e)
540 for (i = 0; i < e; i++)
547 mpz_tdiv_q_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
552 mp_limb_t high_limb, low_limb;
556 lcnt = bcnt / GMP_LIMB_BITS;
557 rn = ABS (as) - lcnt;
569 cnt = bcnt % GMP_LIMB_BITS;
573 tnc = GMP_LIMB_BITS - cnt;
575 low_limb = high_limb >> cnt;
577 for (i = rn - 1; i != 0; i--)
580 *rp++ = low_limb | LO (high_limb << tnc);
581 low_limb = high_limb >> cnt;
589 mpn_copyi (rp, ap, rn);
592 SIZ (r) = as >= 0 ? rn : -rn;
597 mpz_tdiv_r_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
604 bwhole = bcnt / GMP_LIMB_BITS;
605 bcnt %= GMP_LIMB_BITS;
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);
616 mpz_cmp (mpz_t a, mpz_t b)
618 mp_limb_t *ap, *bp, al, bl;
619 int as = SIZ (a), bs = SIZ (b);
624 return as > bs ? 1 : -1;
626 sign = as > 0 ? 1 : -1;
630 for (i = ABS (as) - 1; i >= 0; i--)
635 return al > bl ? sign : -sign;
641 mpz_cmp_ui (mpz_t a, unsigned long b)
645 mpz_init_set_ui (bz, b);
646 ret = mpz_cmp (a, bz);
652 mpz_tdiv_qr (mpz_t q, mpz_t r, mpz_t a, mpz_t b)
657 ASSERT (mpz_sgn(a) >= 0);
658 ASSERT (mpz_sgn(b) > 0);
660 mpz_init_set (tmpr, a);
661 mpz_init_set (tmpb, b);
664 if (mpz_cmp (tmpr, tmpb) > 0)
666 cnt = mpz_sizeinbase (tmpr, 2) - mpz_sizeinbase (tmpb, 2) + 1;
667 mpz_mul_2exp (tmpb, tmpb, cnt);
669 for ( ; cnt > 0; cnt--)
671 mpz_mul_2exp (q, q, 1);
672 mpz_tdiv_q_2exp (tmpb, tmpb, 1L);
673 if (mpz_cmp (tmpr, tmpb) >= 0)
675 mpz_sub (tmpr, tmpr, tmpb);
676 mpz_add_ui (q, q, 1L);
677 ASSERT (mpz_cmp (tmpr, tmpb) < 0);
688 mpz_tdiv_qr_ui (mpz_t q, mpz_t r, mpz_t a, unsigned long b)
691 mpz_init_set_ui (bz, b);
692 mpz_tdiv_qr (q, r, a, bz);
697 mpz_tdiv_q (mpz_t q, mpz_t a, mpz_t b)
702 mpz_tdiv_qr (q, r, a, b);
707 mpz_tdiv_r (mpz_t r, mpz_t a, mpz_t b)
712 mpz_tdiv_qr (q, r, a, b);
717 mpz_tdiv_q_ui (mpz_t q, mpz_t n, unsigned long d)
720 mpz_init_set_ui (dz, d);
721 mpz_tdiv_q (q, n, dz);
725 /* Set inv to the inverse of d, in the style of invert_limb, ie. for
726 udiv_qrnnd_preinv. */
728 mpz_preinv_invert (mpz_t inv, mpz_t d, int numb_bits)
734 norm = numb_bits - mpz_sizeinbase (d, 2);
736 mpz_init_set_ui (t, 1L);
737 mpz_mul_2exp (t, t, 2*numb_bits - norm);
738 mpz_tdiv_q (inv, t, d);
740 mpz_mul_2exp (t, t, numb_bits);
741 mpz_sub (inv, inv, t);
746 /* Remove leading '0' characters from the start of a string, by copying the
749 strstrip_leading_zeros (char *s)
766 mpz_get_str (char *buf, int base, mpz_t a)
768 static char tohex[] = "0123456789abcdef";
770 mp_limb_t alimb, *ap;
780 buf = xmalloc (ABSIZ (a) * (GMP_LIMB_BITS / 4) + 3);
791 bn = an * (GMP_LIMB_BITS / 4);
794 for (i = 0; i < an; i++)
797 for (j = 0; j < GMP_LIMB_BITS / 4; j++)
800 *bp = tohex [alimb & 0xF];
809 strstrip_leading_zeros (buf);
814 mpz_out_str (FILE *file, int base, mpz_t a)
821 str = mpz_get_str (0, 16, a);
826 /* Calculate r satisfying r*d == 1 mod 2^n. */
828 mpz_invert_2exp (mpz_t r, mpz_t a, unsigned long n)
833 ASSERT (mpz_odd_p (a));
835 mpz_init_set_ui (inv, 1L);
838 for (i = 1; i < n; i++)
840 mpz_mul (prod, inv, a);
841 if (mpz_tstbit (prod, i) != 0)
845 mpz_mul (prod, inv, a);
846 mpz_tdiv_r_2exp (prod, prod, n);
847 ASSERT (mpz_cmp_ui (prod, 1L) == 0);
855 /* Calculate inv satisfying r*a == 1 mod 2^n. */
857 mpz_invert_ui_2exp (mpz_t r, unsigned long a, unsigned long n)
860 mpz_init_set_ui (az, a);
861 mpz_invert_2exp (r, az, n);
867 mpz_pow_ui (mpz_t x, mpz_t y, unsigned long z)
871 mpz_init_set_ui (t, 1);
880 mpz_addmul_ui (mpz_t x, mpz_t y, unsigned long z)
885 mpz_mul_ui (t, y, z);
890 /* x=floor(y^(1/z)) */
892 mpz_root (mpz_t x, mpz_t y, unsigned long z)
898 fprintf (stderr, "mpz_root does not accept negative values\n");
901 if (mpz_cmp_ui (y, 1) <= 0)
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)