/* mini-gmp, a minimalistic implementation of a GNU GMP subset.
Contributed to the GNU project by Niels Möller
Copyright 1991-1997, 1999-2019 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 3 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. If not, see http://www.gnu.org/licenses/. */
/* NOTE: All functions in this file which are not declared in
mini-gmp.h are internal, and are not intended to be compatible
neither with GMP nor with future versions of mini-gmp. */
/* Much of the material copied from GMP files, including: gmp-impl.h,
longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c,
mpn/generic/lshift.c, mpn/generic/mul_1.c,
mpn/generic/mul_basecase.c, mpn/generic/rshift.c,
mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c,
mpn/generic/submul_1.c. */
#include <assert.h>
#include <ctype.h>
#include <limits.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "mini-gmp.h"
#if !defined(MINI_GMP_DONT_USE_FLOAT_H)
#include <float.h>
#endif
/* Macros */
#define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
#define GMP_LIMB_MAX ((mp_limb_t) ~ (mp_limb_t) 0)
#define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1))
#define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2))
#define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1)
#define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT)
#define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1))
#define GMP_ABS(x) ((x) >= 0 ? (x) : -(x))
#define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1))
#define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
#define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))
#define GMP_CMP(a,b) (((a) > (b)) - ((a) < (b)))
#if defined(DBL_MANT_DIG) && FLT_RADIX == 2
#define GMP_DBL_MANT_BITS DBL_MANT_DIG
#else
#define GMP_DBL_MANT_BITS (53)
#endif
/* Return non-zero if xp,xsize and yp,ysize overlap.
If xp+xsize<=yp there's no overlap, or if yp+ysize<=xp there's no
overlap. If both these are false, there's an overlap. */
#define GMP_MPN_OVERLAP_P(xp, xsize, yp, ysize) \
((xp) + (xsize) > (yp) && (yp) + (ysize) > (xp))
#define gmp_assert_nocarry(x) do { \
mp_limb_t __cy = (x); \
assert (__cy == 0); \
} while (0)
#define gmp_clz(count, x) do { \
mp_limb_t __clz_x = (x); \
unsigned __clz_c = 0; \
int LOCAL_SHIFT_BITS = 8; \
if (GMP_LIMB_BITS > LOCAL_SHIFT_BITS) \
for (; \
(__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0; \
__clz_c += 8) \
{ __clz_x <<= LOCAL_SHIFT_BITS; } \
for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++) \
__clz_x <<= 1; \
(count) = __clz_c; \
} while (0)
#define gmp_ctz(count, x) do { \
mp_limb_t __ctz_x = (x); \
unsigned __ctz_c = 0; \
gmp_clz (__ctz_c, __ctz_x & - __ctz_x); \
(count) = GMP_LIMB_BITS - 1 - __ctz_c; \
} while (0)
#define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \
do { \
mp_limb_t __x; \
__x = (al) + (bl); \
(sh) = (ah) + (bh) + (__x < (al)); \
(sl) = __x; \
} while (0)
#define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \
do { \
mp_limb_t __x; \
__x = (al) - (bl); \
(sh) = (ah) - (bh) - ((al) < (bl)); \
(sl) = __x; \
} while (0)
#define gmp_umul_ppmm(w1, w0, u, v) \
do { \
int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS; \
if (sizeof(unsigned int) * CHAR_BIT >= 2 * GMP_LIMB_BITS) \
{ \
unsigned int __ww = (unsigned int) (u) * (v); \
w0 = (mp_limb_t) __ww; \
w1 = (mp_limb_t) (__ww >> LOCAL_GMP_LIMB_BITS); \
} \
else if (GMP_ULONG_BITS >= 2 * GMP_LIMB_BITS) \
{ \
unsigned long int __ww = (unsigned long int) (u) * (v); \
w0 = (mp_limb_t) __ww; \
w1 = (mp_limb_t) (__ww >> LOCAL_GMP_LIMB_BITS); \
} \
else { \
mp_limb_t __x0, __x1, __x2, __x3; \
unsigned __ul, __vl, __uh, __vh; \
mp_limb_t __u = (u), __v = (v); \
\
__ul = __u & GMP_LLIMB_MASK; \
__uh = __u >> (GMP_LIMB_BITS / 2); \
__vl = __v & GMP_LLIMB_MASK; \
__vh = __v >> (GMP_LIMB_BITS / 2); \
\
__x0 = (mp_limb_t) __ul * __vl; \
__x1 = (mp_limb_t) __ul * __vh; \
__x2 = (mp_limb_t) __uh * __vl; \
__x3 = (mp_limb_t) __uh * __vh; \
\
__x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */ \
__x1 += __x2; /* but this indeed can */ \
if (__x1 < __x2) /* did we get it? */ \
__x3 += GMP_HLIMB_BIT; /* yes, add it in the proper pos. */ \
\
(w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2)); \
(w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK); \
} \
} while (0)
#define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
do { \
mp_limb_t _qh, _ql, _r, _mask; \
gmp_umul_ppmm (_qh, _ql, (nh), (di)); \
gmp_add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl)); \
_r = (nl) - _qh * (d); \
_mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */ \
_qh += _mask; \
_r += _mask & (d); \
if (_r >= (d)) \
{ \
_r -= (d); \
_qh++; \
} \
\
(r) = _r; \
(q) = _qh; \
} while (0)
#define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv) \
do { \
mp_limb_t _q0, _t1, _t0, _mask; \
gmp_umul_ppmm ((q), _q0, (n2), (dinv)); \
gmp_add_ssaaaa ((q), _q0, (q), _q0, (n2), (n1)); \
\
/* Compute the two most significant limbs of n - q'd */ \
(r1) = (n1) - (d1) * (q); \
gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0)); \
gmp_umul_ppmm (_t1, _t0, (d0), (q)); \
gmp_sub_ddmmss ((r1), (r0), (r1), (r0), _t1, _t0); \
(q)++; \
\
/* Conditionally adjust q and the remainders */ \
_mask = - (mp_limb_t) ((r1) >= _q0); \
(q) += _mask; \
gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \
if ((r1) >= (d1)) \
{ \
if ((r1) > (d1) || (r0) >= (d0)) \
{ \
(q)++; \
gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0)); \
} \
} \
} while (0)
/* Swap macros. */
#define MP_LIMB_T_SWAP(x, y) \
do { \
mp_limb_t __mp_limb_t_swap__tmp = (x); \
(x) = (y); \
(y) = __mp_limb_t_swap__tmp; \
} while (0)
#define MP_SIZE_T_SWAP(x, y) \
do { \
mp_size_t __mp_size_t_swap__tmp = (x); \
(x) = (y); \
(y) = __mp_size_t_swap__tmp; \
} while (0)
#define MP_BITCNT_T_SWAP(x,y) \
do { \
mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x); \
(x) = (y); \
(y) = __mp_bitcnt_t_swap__tmp; \
} while (0)
#define MP_PTR_SWAP(x, y) \
do { \
mp_ptr __mp_ptr_swap__tmp = (x); \
(x) = (y); \
(y) = __mp_ptr_swap__tmp; \
} while (0)
#define MP_SRCPTR_SWAP(x, y) \
do { \
mp_srcptr __mp_srcptr_swap__tmp = (x); \
(x) = (y); \
(y) = __mp_srcptr_swap__tmp; \
} while (0)
#define MPN_PTR_SWAP(xp,xs, yp,ys) \
do { \
MP_PTR_SWAP (xp, yp); \
MP_SIZE_T_SWAP (xs, ys); \
} while(0)
#define MPN_SRCPTR_SWAP(xp,xs, yp,ys) \
do { \
MP_SRCPTR_SWAP (xp, yp); \
MP_SIZE_T_SWAP (xs, ys); \
} while(0)
#define MPZ_PTR_SWAP(x, y) \
do { \
mpz_ptr __mpz_ptr_swap__tmp = (x); \
(x) = (y); \
(y) = __mpz_ptr_swap__tmp; \
} while (0)
#define MPZ_SRCPTR_SWAP(x, y) \
do { \
mpz_srcptr __mpz_srcptr_swap__tmp = (x); \
(x) = (y); \
(y) = __mpz_srcptr_swap__tmp; \
} while (0)
const int mp_bits_per_limb = GMP_LIMB_BITS;
/* Memory allocation and other helper functions. */
static void
gmp_die (const char *msg)
{
fprintf (stderr, "%s\n", msg);
abort();
}
static void *
gmp_default_alloc (size_t size)
{
void *p;
assert (size > 0);
p = malloc (size);
if (!p)
gmp_die("gmp_default_alloc: Virtual memory exhausted.");
return p;
}
static void *
gmp_default_realloc (void *old, size_t unused_old_size, size_t new_size)
{
void * p;
p = realloc (old, new_size);
if (!p)
gmp_die("gmp_default_realloc: Virtual memory exhausted.");
return p;
}
static void
gmp_default_free (void *p, size_t unused_size)
{
free (p);
}
static void * (*gmp_allocate_func) (size_t) = gmp_default_alloc;
static void * (*gmp_reallocate_func) (void *, size_t, size_t) = gmp_default_realloc;
static void (*gmp_free_func) (void *, size_t) = gmp_default_free;
void
mp_get_memory_functions (void *(**alloc_func) (size_t),
void *(**realloc_func) (void *, size_t, size_t),
void (**free_func) (void *, size_t))
{
if (alloc_func)
*alloc_func = gmp_allocate_func;
if (realloc_func)
*realloc_func = gmp_reallocate_func;
if (free_func)
*free_func = gmp_free_func;
}
void
mp_set_memory_functions (void *(*alloc_func) (size_t),
void *(*realloc_func) (void *, size_t, size_t),
void (*free_func) (void *, size_t))
{
if (!alloc_func)
alloc_func = gmp_default_alloc;
if (!realloc_func)
realloc_func = gmp_default_realloc;
if (!free_func)
free_func = gmp_default_free;
gmp_allocate_func = alloc_func;
gmp_reallocate_func = realloc_func;
gmp_free_func = free_func;
}
#define gmp_xalloc(size) ((*gmp_allocate_func)((size)))
#define gmp_free(p) ((*gmp_free_func) ((p), 0))
static mp_ptr
gmp_xalloc_limbs (mp_size_t size)
{
return (mp_ptr) gmp_xalloc (size * sizeof (mp_limb_t));
}
static mp_ptr
gmp_xrealloc_limbs (mp_ptr old, mp_size_t size)
{
assert (size > 0);
return (mp_ptr) (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t));
}
/* MPN interface */
void
mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n)
{
mp_size_t i;
for (i = 0; i < n; i++)
d[i] = s[i];
}
void
mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
{
while (--n >= 0)
d[n] = s[n];
}
int
mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
{
while (--n >= 0)
{
if (ap[n] != bp[n])
return ap[n] > bp[n] ? 1 : -1;
}
return 0;
}
static int
mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
{
if (an != bn)
return an < bn ? -1 : 1;
else
return mpn_cmp (ap, bp, an);
}
static mp_size_t
mpn_normalized_size (mp_srcptr xp, mp_size_t n)
{
while (n > 0 && xp[n-1] == 0)
--n;
return n;
}
int
mpn_zero_p(mp_srcptr rp, mp_size_t n)
{
return mpn_normalized_size (rp, n) == 0;
}
void
mpn_zero (mp_ptr rp, mp_size_t n)
{
while (--n >= 0)
rp[n] = 0;
}
mp_limb_t
mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
{
mp_size_t i;
assert (n > 0);
i = 0;
do
{
mp_limb_t r = ap[i] + b;
/* Carry out */
b = (r < b);
rp[i] = r;
}
while (++i < n);
return b;
}
mp_limb_t
mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
{
mp_size_t i;
mp_limb_t cy;
for (i = 0, cy = 0; i < n; i++)
{
mp_limb_t a, b, r;
a = ap[i]; b = bp[i];
r = a + cy;
cy = (r < cy);
r += b;
cy += (r < b);
rp[i] = r;
}
return cy;
}
mp_limb_t
mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
{
mp_limb_t cy;
assert (an >= bn);
cy = mpn_add_n (rp, ap, bp, bn);
if (an > bn)
cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy);
return cy;
}
mp_limb_t
mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
{
mp_size_t i;
assert (n > 0);
i = 0;
do
{
mp_limb_t a = ap[i];
/* Carry out */
mp_limb_t cy = a < b;
rp[i] = a - b;
b = cy;
}
while (++i < n);
return b;
}
mp_limb_t
mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
{
mp_size_t i;
mp_limb_t cy;
for (i = 0, cy = 0; i < n; i++)
{
mp_limb_t a, b;
a = ap[i]; b = bp[i];
b += cy;
cy = (b < cy);
cy += (a < b);
rp[i] = a - b;
}
return cy;
}
mp_limb_t
mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
{
mp_limb_t cy;
assert (an >= bn);
cy = mpn_sub_n (rp, ap, bp, bn);
if (an > bn)
cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy);
return cy;
}
mp_limb_t
mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
{
mp_limb_t ul, cl, hpl, lpl;
assert (n >= 1);
cl = 0;
do
{
ul = *up++;
gmp_umul_ppmm (hpl, lpl, ul, vl);
lpl += cl;
cl = (lpl < cl) + hpl;
*rp++ = lpl;
}
while (--n != 0);
return cl;
}
mp_limb_t
mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
{
mp_limb_t ul, cl, hpl, lpl, rl;
assert (n >= 1);
cl = 0;
do
{
ul = *up++;
gmp_umul_ppmm (hpl, lpl, ul, vl);
lpl += cl;
cl = (lpl < cl) + hpl;
rl = *rp;
lpl = rl + lpl;
cl += lpl < rl;
*rp++ = lpl;
}
while (--n != 0);
return cl;
}
mp_limb_t
mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
{
mp_limb_t ul, cl, hpl, lpl, rl;
assert (n >= 1);
cl = 0;
do
{
ul = *up++;
gmp_umul_ppmm (hpl, lpl, ul, vl);
lpl += cl;
cl = (lpl < cl) + hpl;
rl = *rp;
lpl = rl - lpl;
cl += lpl > rl;
*rp++ = lpl;
}
while (--n != 0);
return cl;
}
mp_limb_t
mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
{
assert (un >= vn);
assert (vn >= 1);
assert (!GMP_MPN_OVERLAP_P(rp, un + vn, up, un));
assert (!GMP_MPN_OVERLAP_P(rp, un + vn, vp, vn));
/* We first multiply by the low order limb. This result can be
stored, not added, to rp. We also avoid a loop for zeroing this
way. */
rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
/* Now accumulate the product of up[] and the next higher limb from
vp[]. */
while (--vn >= 1)
{
rp += 1, vp += 1;
rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
}
return rp[un];
}
void
mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
{
mpn_mul (rp, ap, n, bp, n);
}
void
mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n)
{
mpn_mul (rp, ap, n, ap, n);
}
mp_limb_t
mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
{
mp_limb_t high_limb, low_limb;
unsigned int tnc;
mp_limb_t retval;
assert (n >= 1);
assert (cnt >= 1);
assert (cnt < GMP_LIMB_BITS);
up += n;
rp += n;
tnc = GMP_LIMB_BITS - cnt;
low_limb = *--up;
retval = low_limb >> tnc;
high_limb = (low_limb << cnt);
while (--n != 0)
{
low_limb = *--up;
*--rp = high_limb | (low_limb >> tnc);
high_limb = (low_limb << cnt);
}
*--rp = high_limb;
return retval;
}
mp_limb_t
mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
{
mp_limb_t high_limb, low_limb;
unsigned int tnc;
mp_limb_t retval;
assert (n >= 1);
assert (cnt >= 1);
assert (cnt < GMP_LIMB_BITS);
tnc = GMP_LIMB_BITS - cnt;
high_limb = *up++;
retval = (high_limb << tnc);
low_limb = high_limb >> cnt;
while (--n != 0)
{
high_limb = *up++;
*rp++ = low_limb | (high_limb << tnc);
low_limb = high_limb >> cnt;
}
*rp = low_limb;
return retval;
}
static mp_bitcnt_t
mpn_common_scan (mp_limb_t limb, mp_size_t i, mp_srcptr up, mp_size_t un,
mp_limb_t ux)
{
unsigned cnt;
assert (ux == 0 || ux == GMP_LIMB_MAX);
assert (0 <= i && i <= un );
while (limb == 0)
{
i++;
if (i == un)
return (ux == 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
limb = ux ^ up[i];
}
gmp_ctz (cnt, limb);
|