cheng/mpir
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
1516ccfbf2
mpir/dumbmp.c
jasonmoxham f405d0fe68 remove ansitoknr stuff
2009-08-13 03:07:39 +00:00

890 lines
16 KiB
C
Raw Blame History

/* dumbmp mini GMP compatible library.
Copyright 2001, 2002, 2004 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library; see the file COPYING.LIB. If not, write to
the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
MA 02110-1301, USA. */
/* The code here implements a subset (a very limited subset) of the main GMP
functions. It's designed for use in a few build-time calculations and
will be slow, but highly portable.
None of the normal GMP configure things are used, nor any of the normal
mpir.h or gmp-impl.h. To use this file in a program just #include
"dumbmp.c".
ANSI function definitions can be used here but other ANSI-isms like
"const" should be avoided.
mp_limb_t here is an unsigned long, since that's a sensible type
everywhere we know of, with 8*sizeof(unsigned long) giving the number of
bits in the type (that not being true for instance with int or short on
Cray vector systems.)
Only the low half of each mp_limb_t is used, so as to make carry handling
and limb multiplies easy. GMP_LIMB_BITS is the number of bits used. */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
typedef unsigned long mp_limb_t;
typedef struct {
int _mp_alloc;
int _mp_size;
mp_limb_t *_mp_d;
} mpz_t[1];
#define GMP_LIMB_BITS (sizeof (mp_limb_t) * 8 / 2)
#define ABS(x) ((x) >= 0 ? (x) : -(x))
#define MIN(l,o) ((l) < (o) ? (l) : (o))
#define MAX(h,i) ((h) > (i) ? (h) : (i))
#define ALLOC(x) ((x)->_mp_alloc)
#define PTR(x) ((x)->_mp_d)
#define SIZ(x) ((x)->_mp_size)
#define ABSIZ(x) ABS (SIZ (x))
#define LOMASK ((1L << GMP_LIMB_BITS) - 1)
#define LO(x) ((x) & LOMASK)
#define HI(x) ((x) >> GMP_LIMB_BITS)
#define ASSERT(cond) \
do { \
if (! (cond)) \
{ \
fprintf (stderr, "Assertion failure\n"); \
abort (); \
} \
} while (0)
char *
xmalloc (int n)
{
char *p;
p = malloc (n);
if (p == NULL)
{
fprintf (stderr, "Out of memory (alloc %d bytes)\n", n);
abort ();
}
return p;
}
mp_limb_t *
xmalloc_limbs (int n)
{
return (mp_limb_t *) xmalloc (n * sizeof (mp_limb_t));
}
void
mem_copyi (char *dst, char *src, int size)
{
int i;
for (i = 0; i < size; i++)
dst[i] = src[i];
}
int
isprime (int n)
{
int i;
if (n < 2)
return 0;
for (i = 2; i < n; i++)
if ((n % i) == 0)
return 0;
return 1;
}
int
log2_ceil (int n)
{
int e;
ASSERT (n >= 1);
for (e = 0; ; e++)
if ((1 << e) >= n)
break;
return e;
}
void
mpz_realloc (mpz_t r, int n)
{
if (n <= ALLOC(r))
return;
ALLOC(r) = n;
PTR(r) = (mp_limb_t *) realloc (PTR(r), n * sizeof (mp_limb_t));
if (PTR(r) == NULL)
{
fprintf (stderr, "Out of memory (realloc to %d)\n", n);
abort ();
}
}
void
mpn_normalize (mp_limb_t *rp, int *rnp)
{
int rn = *rnp;
while (rn > 0 && rp[rn-1] == 0)
rn--;
*rnp = rn;
}
void
mpn_copyi (mp_limb_t *dst, mp_limb_t *src, int n)
{
int i;
for (i = 0; i < n; i++)
dst[i] = src[i];
}
void
mpn_zero (mp_limb_t *rp, int rn)
{
int i;
for (i = 0; i < rn; i++)
rp[i] = 0;
}
void
mpz_init (mpz_t r)
{
ALLOC(r) = 1;
PTR(r) = xmalloc_limbs (ALLOC(r));
PTR(r)[0] = 0;
SIZ(r) = 0;
}
void
mpz_clear (mpz_t r)
{
free (PTR (r));
ALLOC(r) = -1;
SIZ (r) = 0xbadcafeL;
PTR (r) = (mp_limb_t *) 0xdeadbeefL;
}
int
mpz_sgn (mpz_t a)
{
return (SIZ(a) > 0 ? 1 : SIZ(a) == 0 ? 0 : -1);
}
int
mpz_odd_p (mpz_t a)
{
if (SIZ(a) == 0)
return 0;
else
return (PTR(a)[0] & 1) != 0;
}
int
mpz_even_p (mpz_t a)
{
if (SIZ(a) == 0)
return 1;
else
return (PTR(a)[0] & 1) == 0;
}
size_t
mpz_sizeinbase (mpz_t a, int base)
{
int an = ABSIZ (a);
mp_limb_t *ap = PTR (a);
int cnt;
mp_limb_t hi;
if (base != 2)
abort ();
if (an == 0)
return 1;
cnt = 0;
for (hi = ap[an - 1]; hi != 0; hi >>= 1)
cnt += 1;
return (an - 1) * GMP_LIMB_BITS + cnt;
}
void
mpz_set (mpz_t r, mpz_t a)
{
mpz_realloc (r, ABSIZ (a));
SIZ(r) = SIZ(a);
mpn_copyi (PTR(r), PTR(a), ABSIZ (a));
}
void
mpz_init_set (mpz_t r, mpz_t a)
{
mpz_init (r);
mpz_set (r, a);
}
void
mpz_set_ui (mpz_t r, unsigned long ui)
{
int rn;
mpz_realloc (r, 2);
PTR(r)[0] = LO(ui);
PTR(r)[1] = HI(ui);
rn = 2;
mpn_normalize (PTR(r), &rn);
SIZ(r) = rn;
}
void
mpz_init_set_ui (mpz_t r, unsigned long ui)
{
mpz_init (r);
mpz_set_ui (r, ui);
}
void
mpz_setbit (mpz_t r, unsigned long bit)
{
int limb, rn, extend;
mp_limb_t *rp;
rn = SIZ(r);
if (rn < 0)
abort (); /* only r>=0 */
limb = bit / GMP_LIMB_BITS;
bit %= GMP_LIMB_BITS;
mpz_realloc (r, limb+1);
rp = PTR(r);
extend = (limb+1) - rn;
if (extend > 0)
mpn_zero (rp + rn, extend);
rp[limb] |= (mp_limb_t) 1 << bit;
SIZ(r) = MAX (rn, limb+1);
}
int
mpz_tstbit (mpz_t r, unsigned long bit)
{
int limb;
if (SIZ(r) < 0)
abort (); /* only r>=0 */
limb = bit / GMP_LIMB_BITS;
if (SIZ(r) <= limb)
return 0;
bit %= GMP_LIMB_BITS;
return (PTR(r)[limb] >> bit) & 1;
}
int
popc_limb (mp_limb_t a)
{
int ret = 0;
while (a != 0)
{
ret += (a & 1);
a >>= 1;
}
return ret;
}
unsigned long
mpz_popcount (mpz_t a)
{
unsigned long ret;
int i;
if (SIZ(a) < 0)
abort ();
ret = 0;
for (i = 0; i < SIZ(a); i++)
ret += popc_limb (PTR(a)[i]);
return ret;
}
void
mpz_add (mpz_t r, mpz_t a, mpz_t b)
{
int an = ABSIZ (a), bn = ABSIZ (b), rn;
mp_limb_t *rp, *ap, *bp;
int i;
mp_limb_t t, cy;
if ((SIZ (a) ^ SIZ (b)) < 0)
abort (); /* really subtraction */
if (SIZ (a) < 0)
abort ();
mpz_realloc (r, MAX (an, bn) + 1);
ap = PTR (a); bp = PTR (b); rp = PTR (r);
if (an < bn)
{
mp_limb_t *tp; int tn;
tn = an; an = bn; bn = tn;
tp = ap; ap = bp; bp = tp;
}
cy = 0;
for (i = 0; i < bn; i++)
{
t = ap[i] + bp[i] + cy;
rp[i] = LO (t);
cy = HI (t);
}
for (i = bn; i < an; i++)
{
t = ap[i] + cy;
rp[i] = LO (t);
cy = HI (t);
}
rp[an] = cy;
rn = an + 1;
mpn_normalize (rp, &rn);
SIZ (r) = rn;
}
void
mpz_add_ui (mpz_t r, mpz_t a, unsigned long int ui)
{
mpz_t b;
mpz_init (b);
mpz_set_ui (b, ui);
mpz_add (r, a, b);
mpz_clear (b);
}
void
mpz_sub (mpz_t r, mpz_t a, mpz_t b)
{
int an = ABSIZ (a), bn = ABSIZ (b), rn;
mp_limb_t *rp, *ap, *bp;
int i;
mp_limb_t t, cy;
if ((SIZ (a) ^ SIZ (b)) < 0)
abort (); /* really addition */
if (SIZ (a) < 0)
abort ();
mpz_realloc (r, MAX (an, bn) + 1);
ap = PTR (a); bp = PTR (b); rp = PTR (r);
if (an < bn)
{
mp_limb_t *tp; int tn;
tn = an; an = bn; bn = tn;
tp = ap; ap = bp; bp = tp;
}
cy = 0;
for (i = 0; i < bn; i++)
{
t = ap[i] - bp[i] - cy;
rp[i] = LO (t);
cy = LO (-HI (t));
}
for (i = bn; i < an; i++)
{
t = ap[i] - cy;
rp[i] = LO (t);
cy = LO (-HI (t));
}
rp[an] = cy;
rn = an + 1;
if (cy != 0)
{
cy = 0;
for (i = 0; i < rn; i++)
{
t = -rp[i] - cy;
rp[i] = LO (t);
cy = LO (-HI (t));
}
SIZ (r) = -rn;
return;
}
mpn_normalize (rp, &rn);
SIZ (r) = rn;
}
void
mpz_sub_ui (mpz_t r, mpz_t a, unsigned long int ui)
{
mpz_t b;
mpz_init (b);
mpz_set_ui (b, ui);
mpz_sub (r, a, b);
mpz_clear (b);
}
void
mpz_mul (mpz_t r, mpz_t a, mpz_t b)
{
int an = ABSIZ (a), bn = ABSIZ (b), rn;
mp_limb_t *scratch, *tmp, *ap = PTR (a), *bp = PTR (b);
int ai, bi;
mp_limb_t t, cy;
scratch = xmalloc_limbs (an + bn);
tmp = scratch;
for (bi = 0; bi < bn; bi++)
tmp[bi] = 0;
for (ai = 0; ai < an; ai++)
{
tmp = scratch + ai;
cy = 0;
for (bi = 0; bi < bn; bi++)
{
t = ap[ai] * bp[bi] + tmp[bi] + cy;
tmp[bi] = LO (t);
cy = HI (t);
}
tmp[bn] = cy;
}
rn = an + bn;
mpn_normalize (scratch, &rn);
free (PTR (r));
PTR (r) = scratch;
SIZ (r) = (SIZ (a) ^ SIZ (b)) >= 0 ? rn : -rn;
ALLOC (r) = an + bn;
}
void
mpz_mul_ui (mpz_t r, mpz_t a, unsigned long int ui)
{
mpz_t b;
mpz_init (b);
mpz_set_ui (b, ui);
mpz_mul (r, a, b);
mpz_clear (b);
}
void
mpz_mul_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
{
mpz_set (r, a);
while (bcnt)
{
mpz_add (r, r, r);
bcnt -= 1;
}
}
void
mpz_ui_pow_ui (mpz_t r, unsigned long b, unsigned long e)
{
unsigned long i;
mpz_t bz;
mpz_init (bz);
mpz_set_ui (bz, b);
mpz_set_ui (r, 1L);
for (i = 0; i < e; i++)
mpz_mul (r, r, bz);
mpz_clear (bz);
}
void
mpz_tdiv_q_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
{
int as, rn;
int cnt, tnc;
int lcnt;
mp_limb_t high_limb, low_limb;
int i;
as = SIZ (a);
lcnt = bcnt / GMP_LIMB_BITS;
rn = ABS (as) - lcnt;
if (rn <= 0)
SIZ (r) = 0;
else
{
mp_limb_t *rp, *ap;
mpz_realloc (r, rn);
rp = PTR (r);
ap = PTR (a);
cnt = bcnt % GMP_LIMB_BITS;
if (cnt != 0)
{
ap += lcnt;
tnc = GMP_LIMB_BITS - cnt;
high_limb = *ap++;
low_limb = high_limb >> cnt;
for (i = rn - 1; i != 0; i--)
{
high_limb = *ap++;
*rp++ = low_limb | LO (high_limb << tnc);
low_limb = high_limb >> cnt;
}
*rp = low_limb;
rn -= low_limb == 0;
}
else
{
ap += lcnt;
mpn_copyi (rp, ap, rn);
}
SIZ (r) = as >= 0 ? rn : -rn;
}
}
void
mpz_tdiv_r_2exp (mpz_t r, mpz_t a, unsigned long int bcnt)
{
int rn, bwhole;
mpz_set (r, a);
rn = ABSIZ(r);
bwhole = bcnt / GMP_LIMB_BITS;
bcnt %= GMP_LIMB_BITS;
if (rn > bwhole)
{
rn = bwhole+1;
PTR(r)[rn-1] &= ((mp_limb_t) 1 << bcnt) - 1;
mpn_normalize (PTR(r), &rn);
SIZ(r) = (SIZ(r) >= 0 ? rn : -rn);
}
}
int
mpz_cmp (mpz_t a, mpz_t b)
{
mp_limb_t *ap, *bp, al, bl;
int as = SIZ (a), bs = SIZ (b);
int i;
int sign;
if (as != bs)
return as > bs ? 1 : -1;
sign = as > 0 ? 1 : -1;
ap = PTR (a);
bp = PTR (b);
for (i = ABS (as) - 1; i >= 0; i--)
{
al = ap[i];
bl = bp[i];
if (al != bl)
return al > bl ? sign : -sign;
}
return 0;
}
int
mpz_cmp_ui (mpz_t a, unsigned long b)
{
mpz_t bz;
int ret;
mpz_init_set_ui (bz, b);
ret = mpz_cmp (a, bz);
mpz_clear (bz);
return ret;
}
void
mpz_tdiv_qr (mpz_t q, mpz_t r, mpz_t a, mpz_t b)
{
mpz_t tmpr, tmpb;
unsigned long cnt;
ASSERT (mpz_sgn(a) >= 0);
ASSERT (mpz_sgn(b) > 0);
mpz_init_set (tmpr, a);
mpz_init_set (tmpb, b);
mpz_set_ui (q, 0L);
if (mpz_cmp (tmpr, tmpb) > 0)
{
cnt = mpz_sizeinbase (tmpr, 2) - mpz_sizeinbase (tmpb, 2) + 1;
mpz_mul_2exp (tmpb, tmpb, cnt);
for ( ; cnt > 0; cnt--)
{
mpz_mul_2exp (q, q, 1);
mpz_tdiv_q_2exp (tmpb, tmpb, 1L);
if (mpz_cmp (tmpr, tmpb) >= 0)
{
mpz_sub (tmpr, tmpr, tmpb);
mpz_add_ui (q, q, 1L);
ASSERT (mpz_cmp (tmpr, tmpb) < 0);
}
}
}
mpz_set (r, tmpr);
mpz_clear (tmpr);
mpz_clear (tmpb);
}
void
mpz_tdiv_qr_ui (mpz_t q, mpz_t r, mpz_t a, unsigned long b)
{
mpz_t bz;
mpz_init_set_ui (bz, b);
mpz_tdiv_qr (q, r, a, bz);
mpz_clear (bz);
}
void
mpz_tdiv_q (mpz_t q, mpz_t a, mpz_t b)
{
mpz_t r;
mpz_init (r);
mpz_tdiv_qr (q, r, a, b);
mpz_clear (r);
}
void
mpz_tdiv_q_ui (mpz_t q, mpz_t n, unsigned long d)
{
mpz_t dz;
mpz_init_set_ui (dz, d);
mpz_tdiv_q (q, n, dz);
mpz_clear (dz);
}
/* Set inv to the inverse of d, in the style of invert_limb, ie. for
udiv_qrnnd_preinv. */
void
mpz_preinv_invert (mpz_t inv, mpz_t d, int numb_bits)
{
mpz_t t;
int norm;
ASSERT (SIZ(d) > 0);
norm = numb_bits - mpz_sizeinbase (d, 2);
ASSERT (norm >= 0);
mpz_init_set_ui (t, 1L);
mpz_mul_2exp (t, t, 2*numb_bits - norm);
mpz_tdiv_q (inv, t, d);
mpz_set_ui (t, 1L);
mpz_mul_2exp (t, t, numb_bits);
mpz_sub (inv, inv, t);
mpz_clear (t);
}
/* Remove leading '0' characters from the start of a string, by copying the
remainder down. */
void
strstrip_leading_zeros (char *s)
{
char c, *p;
p = s;
while (*s == '0')
s++;
do
{
c = *s++;
*p++ = c;
}
while (c != '\0');
}
char *
mpz_get_str (char *buf, int base, mpz_t a)
{
static char tohex[] = "0123456789abcdef";
mp_limb_t alimb, *ap;
int an, bn, i, j;
char *bp;
if (base != 16)
abort ();
if (SIZ (a) < 0)
abort ();
if (buf == 0)
buf = xmalloc (ABSIZ (a) * (GMP_LIMB_BITS / 4) + 3);
an = ABSIZ (a);
if (an == 0)
{
buf[0] = '0';
buf[1] = '\0';
return buf;
}
ap = PTR (a);
bn = an * (GMP_LIMB_BITS / 4);
bp = buf + bn;
for (i = 0; i < an; i++)
{
alimb = ap[i];
for (j = 0; j < GMP_LIMB_BITS / 4; j++)
{
bp--;
*bp = tohex [alimb & 0xF];
alimb >>= 4;
}
ASSERT (alimb == 0);
}
ASSERT (bp == buf);
buf[bn] = '\0';
strstrip_leading_zeros (buf);
return buf;
}
void
mpz_out_str (FILE *file, int base, mpz_t a)
{
char *str;
if (file == 0)
file = stdout;
str = mpz_get_str (0, 16, a);
fputs (str, file);
free (str);
}
/* Calculate r satisfying r*d == 1 mod 2^n. */
void
mpz_invert_2exp (mpz_t r, mpz_t a, unsigned long n)
{
unsigned long i;
mpz_t inv, prod;
ASSERT (mpz_odd_p (a));
mpz_init_set_ui (inv, 1L);
mpz_init (prod);
for (i = 1; i < n; i++)
{
mpz_mul (prod, inv, a);
if (mpz_tstbit (prod, i) != 0)
mpz_setbit (inv, i);
}
mpz_mul (prod, inv, a);
mpz_tdiv_r_2exp (prod, prod, n);
ASSERT (mpz_cmp_ui (prod, 1L) == 0);
mpz_set (r, inv);
mpz_clear (inv);
mpz_clear (prod);
}
/* Calculate inv satisfying r*a == 1 mod 2^n. */
void
mpz_invert_ui_2exp (mpz_t r, unsigned long a, unsigned long n)
{
mpz_t az;
mpz_init_set_ui (az, a);
mpz_invert_2exp (r, az, n);
mpz_clear (az);
}
/* x=y^z */
void
mpz_pow_ui (mpz_t x, mpz_t y, unsigned long z)
{
mpz_t t;
mpz_init_set_ui (t, 1);
for (; z != 0; z--)
mpz_mul (t, t, y);
mpz_set (x, t);
mpz_clear (t);
}
/* x=x+y*z */
void
mpz_addmul_ui (mpz_t x, mpz_t y, unsigned long z)
{
mpz_t t;
mpz_init (t);
mpz_mul_ui (t, y, z);
mpz_add (x, x, t);
mpz_clear (t);
}
/* x=floor(y^(1/z)) */
void
mpz_root (mpz_t x, mpz_t y, unsigned long z)
{
mpz_t t, u;
if (mpz_sgn (y) < 0)
{
fprintf (stderr, "mpz_root does not accept negative values\n");
abort ();
}
if (mpz_cmp_ui (y, 1) <= 0)
{
mpz_set (x, y);
return;
}
mpz_init (t);
mpz_init_set (u, y);
do
{
mpz_pow_ui (t, u, z - 1);
mpz_tdiv_q (t, y, t);
mpz_addmul_ui (t, u, z - 1);
mpz_tdiv_q_ui (t, t, z);
if (mpz_cmp (t, u) >= 0)
break;
mpz_set (u, t);
}
while (1);
mpz_set (x, u);
mpz_clear (t);
mpz_clear (u);
}
Reference in New Issue View Git Blame Copy Permalink