Viewing file:  mpfr-mini-gmp.c (8.15 KB) -rw-r--r--
Select action/file-type:
 (+) |  (+) |  (+) | Code (+) | Session (+) |  (+) | SDB (+) |  (+) |  (+) |  (+) |  (+) |  (+) |
/* mpfr-mini-gmp.c -- Interface functions for mini-gmp.
Copyright 2014-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
The GNU MPFR 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 3 of the License, or (at your
option) any later version.
The GNU MPFR 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 MPFR Library; see the file COPYING.LESSER. If not, see
https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
/* The following include will do 2 things: include the config.h
if there is one (as it may define MPFR_USE_MINI_GMP), and avoid
an empty translation unit (see ISO C99, 6.9). */
#include "mpfr-impl.h"
#ifdef MPFR_USE_MINI_GMP
/************************ random generation functions ************************/
#ifdef WANT_gmp_randinit_default
void
gmp_randinit_default (gmp_randstate_t state)
{
}
#endif
#ifdef WANT_gmp_randseed_ui
void
gmp_randseed_ui (gmp_randstate_t state, unsigned long int seed)
{
/* With a portable version of the conversion from unsigned long to long
(at least GCC and Clang optimize this expression to identity). */
srand48 (seed > LONG_MAX ? -1 - (long) ~seed : (long) seed);
}
#endif
#ifdef WANT_gmp_randclear
void
gmp_randclear (gmp_randstate_t state)
{
}
#endif
#ifdef WANT_gmp_randinit_set
void
gmp_randinit_set (gmp_randstate_t s1, gmp_randstate_t s2)
{
}
#endif
static mp_limb_t
random_limb (void)
{
/* lrand48() returns a random number in [0, 2^31-1],
but the low 15 bits do not depend on the random seed,
thus it is safer to use the upper bits */
#if GMP_NUMB_BITS < 32
/* use the upper GMP_NUMB_BITS bits from lrand48 () */
return (mp_limb_t) (lrand48 () >> (31 - GMP_NUMB_BITS));
#elif GMP_NUMB_BITS == 32
/* use the upper 16 bits from two lrand48 calls */
return (lrand48 () >> 15) + ((lrand48 () >> 15) << 16);
#elif GMP_NUMB_BITS == 64
/* use the upper 16 bits from four lrand48 calls */
return (lrand48 () >> 15) + ((((mp_limb_t) lrand48 ()) >> 15) << 16)
+ ((((mp_limb_t) lrand48 ()) >> 15) << 32)
+ ((((mp_limb_t) lrand48 ()) >> 15) << 48);
#else
#error "GMP_NUMB_BITS should be 8, 16, 32 or >= 64"
#endif
}
#ifdef WANT_mpz_urandomb
void
mpz_urandomb (mpz_t rop, gmp_randstate_t state, mp_bitcnt_t nbits)
{
unsigned long n, i;
mpz_realloc2 (rop, nbits);
n = (nbits - 1) / GMP_NUMB_BITS + 1; /* number of limbs */
for (i = n; i-- > 0;)
rop->_mp_d[i] = random_limb ();
i = n * GMP_NUMB_BITS - nbits;
/* mask the upper i bits */
if (i)
rop->_mp_d[n-1] = MPFR_LIMB_LSHIFT(rop->_mp_d[n-1], i) >> i;
while (n > 0 && (rop->_mp_d[n-1] == 0))
n--;
rop->_mp_size = n;
}
#endif
#ifdef WANT_gmp_urandomm_ui
/* generates a random unsigned long */
static unsigned long
random_ulong (void)
{
#ifdef MPFR_LONG_WITHIN_LIMB
/* we assume a limb and an unsigned long have both a number of different
values that is a power of two, thus when we cast a random limb into
an unsigned long, we still get an uniform distribution */
return random_limb ();
#else
/* with the same assumption as above, we need to generate as many random
limbs needed to "fill" an unsigned long */
unsigned long u, v;
v = MPFR_LIMB_MAX;
u = random_limb ();
while (v < ULONG_MAX)
{
v = (v << GMP_NUMB_BITS) + MPFR_LIMB_MAX;
u = (u << GMP_NUMB_BITS) + random_limb ();
}
return u;
#endif
}
unsigned long
gmp_urandomm_ui (gmp_randstate_t state, unsigned long n)
{
unsigned long p, q, r;
MPFR_ASSERTN (n > 0);
p = random_ulong (); /* p is in [0, ULONG_MAX], thus p is uniform among
ULONG_MAX+1 values */
q = n * (ULONG_MAX / n);
r = ULONG_MAX % n;
if (r != n - 1) /* ULONG_MAX+1 is not multiple of n, will happen whenever
n is not a power of two */
while (p >= q)
p = random_ulong ();
return p % n;
}
#endif
#ifdef WANT_gmp_urandomb_ui
unsigned long
gmp_urandomb_ui (gmp_randstate_t state, unsigned long n)
{
#ifdef MPFR_LONG_WITHIN_LIMB
return random_limb () % (1UL << n);
#else
unsigned long res = 0;
int m = n; /* remaining bits to generate */
while (m >= GMP_NUMB_BITS)
{
/* we can generate a full limb */
res = (res << GMP_NUMB_BITS) | (unsigned long) random_limb ();
m -= GMP_NUMB_BITS;
}
/* now m < GMP_NUMB_BITS */
if (m) /* generate m extra bits */
res = (res << m) | (unsigned long) (random_limb () % (1UL << m));
return res;
#endif
}
#endif
/************************* division functions ********************************/
#ifdef WANT_mpn_divrem_1
mp_limb_t
mpn_divrem_1 (mp_limb_t *qp, mp_size_t qxn, mp_limb_t *np, mp_size_t nn,
mp_limb_t d0)
{
mpz_t q, r, n, d;
mp_limb_t ret, dd[1];
d->_mp_d = dd;
d->_mp_d[0] = d0;
d->_mp_size = 1;
mpz_init (q);
mpz_init (r);
if (qxn == 0)
{
n->_mp_d = np;
n->_mp_size = nn;
}
else
{
mpz_init2 (n, (nn + qxn) * GMP_NUMB_BITS);
mpn_copyi (n->_mp_d + qxn, np, nn);
mpn_zero (n->_mp_d, qxn);
n->_mp_size = nn + qxn;
}
mpz_tdiv_qr (q, r, n, d);
if (q->_mp_size > 0)
mpn_copyi (qp, q->_mp_d, q->_mp_size);
if (q->_mp_size < nn + qxn)
mpn_zero (qp + q->_mp_size, nn + qxn - q->_mp_size);
ret = (r->_mp_size == 1) ? r->_mp_d[0] : 0;
mpz_clear (q);
mpz_clear (r);
if (qxn != 0)
mpz_clear (n);
return ret;
}
#endif
#ifdef WANT_mpn_divrem
mp_limb_t
mpn_divrem (mp_limb_t *qp, mp_size_t qn, mp_limb_t *np,
mp_size_t nn, const mp_limb_t *dp, mp_size_t dn)
{
mpz_t q, r, n, d;
mp_limb_t ret;
MPFR_ASSERTN(qn == 0);
qn = nn - dn;
n->_mp_d = np;
n->_mp_size = nn;
d->_mp_d = (mp_limb_t*) dp;
d->_mp_size = dn;
mpz_init (q);
mpz_init (r);
mpz_tdiv_qr (q, r, n, d);
MPFR_ASSERTN(q->_mp_size == qn || q->_mp_size == qn + 1);
mpn_copyi (qp, q->_mp_d, qn);
ret = (q->_mp_size == qn) ? 0 : q->_mp_d[qn];
if (r->_mp_size > 0)
mpn_copyi (np, r->_mp_d, r->_mp_size);
if (r->_mp_size < dn)
mpn_zero (np + r->_mp_size, dn - r->_mp_size);
mpz_clear (q);
mpz_clear (r);
return ret;
}
#endif
#ifdef WANT_mpn_tdiv_qr
void
mpn_tdiv_qr (mp_limb_t *qp, mp_limb_t *rp, mp_size_t qxn,
const mp_limb_t *np, mp_size_t nn,
const mp_limb_t *dp, mp_size_t dn)
{
mpz_t q, r, n, d;
MPFR_ASSERTN(qxn == 0);
n->_mp_d = (mp_limb_t*) np;
n->_mp_size = nn;
d->_mp_d = (mp_limb_t*) dp;
d->_mp_size = dn;
mpz_init (q);
mpz_init (r);
mpz_tdiv_qr (q, r, n, d);
MPFR_ASSERTN(q->_mp_size > 0);
mpn_copyi (qp, q->_mp_d, q->_mp_size);
if (q->_mp_size < nn - dn + 1)
mpn_zero (qp + q->_mp_size, nn - dn + 1 - q->_mp_size);
if (r->_mp_size > 0)
mpn_copyi (rp, r->_mp_d, r->_mp_size);
if (r->_mp_size < dn)
mpn_zero (rp + r->_mp_size, dn - r->_mp_size);
mpz_clear (q);
mpz_clear (r);
}
#endif
#if 0 /* this function is useful for debugging, thus please keep it here */
void
mpz_dump (mpz_t z)
{
mp_size_t n = z->_mp_size;
MPFR_STAT_STATIC_ASSERT ((GMP_NUMB_BITS % 4) == 0);
if (n == 0)
printf ("0");
else
{
int first = 1;
if (n < 0)
{
printf ("-");
n = -n;
}
while (n > 0)
{
if (first)
{
printf ("%lx", (unsigned long) z->_mp_d[n-1]);
first = 0;
}
else
{
char s[17];
int len;
sprintf (s, "%lx", (unsigned long) z->_mp_d[n-1]);
len = strlen (s);
/* one character should be printed for 4 bits */
while (len++ < GMP_NUMB_BITS / 4)
printf ("0");
printf ("%lx", (unsigned long) z->_mp_d[n-1]);
}
n--;
}
}
printf ("\n");
}
#endif
#endif /* MPFR_USE_MINI_GMP */
|
