#include <stdio.h>
#include <stdlib.h>
#include "gmp.h"
#include "gmp-impl.h"
#include "tests.h"
void
refmpf_add (mpf_ptr w, mpf_srcptr u, mpf_srcptr v)
{
mp_size_t hi, lo, size;
mp_ptr ut, vt, wt;
int neg;
mp_exp_t exp;
mp_limb_t cy;
TMP_DECL;
TMP_MARK;
if (SIZ (u) == 0)
{
size = ABSIZ (v);
wt = TMP_ALLOC_LIMBS (size + 1);
MPN_COPY (wt, PTR (v), size);
exp = EXP (v);
neg = SIZ (v) < 0;
goto done;
}
if (SIZ (v) == 0)
{
size = ABSIZ (u);
wt = TMP_ALLOC_LIMBS (size + 1);
MPN_COPY (wt, PTR (u), size);
exp = EXP (u);
neg = SIZ (u) < 0;
goto done;
}
if ((SIZ (u) ^ SIZ (v)) < 0)
{
mpf_t tmp;
SIZ (tmp) = -SIZ (v);
EXP (tmp) = EXP (v);
PTR (tmp) = PTR (v);
refmpf_sub (w, u, tmp);
return;
}
neg = SIZ (u) < 0;
hi = MAX (EXP (u), EXP (v));
lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v));
size = hi - lo;
ut = TMP_ALLOC_LIMBS (size + 1);
vt = TMP_ALLOC_LIMBS (size + 1);
wt = TMP_ALLOC_LIMBS (size + 1);
MPN_ZERO (ut, size);
MPN_ZERO (vt, size);
{int off;
off = size + (EXP (u) - hi) - ABSIZ (u);
MPN_COPY (ut + off, PTR (u), ABSIZ (u));
off = size + (EXP (v) - hi) - ABSIZ (v);
MPN_COPY (vt + off, PTR (v), ABSIZ (v));
}
cy = mpn_add_n (wt, ut, vt, size);
wt[size] = cy;
size += cy;
exp = hi + cy;
done:
if (size > PREC (w))
{
wt += size - PREC (w);
size = PREC (w);
}
MPN_COPY (PTR (w), wt, size);
SIZ (w) = neg == 0 ? size : -size;
EXP (w) = exp;
TMP_FREE;
}
void
refmpf_add_ulp (mpf_ptr f)
{
mp_ptr fp = PTR(f);
mp_size_t fsize = SIZ(f);
mp_size_t abs_fsize = ABSIZ(f);
mp_limb_t c;
if (fsize == 0)
{
printf ("Oops, refmpf_add_ulp called with f==0\n");
abort ();
}
c = refmpn_add_1 (fp, fp, abs_fsize, CNST_LIMB(1));
if (c != 0)
{
if (abs_fsize >= PREC(f) + 1)
{
printf ("Oops, refmpf_add_ulp carried out of prec+1 limbs\n");
abort ();
}
fp[abs_fsize] = c;
abs_fsize++;
SIZ(f) = (fsize > 0 ? abs_fsize : - abs_fsize);
EXP(f)++;
}
}
void
refmpf_fill (mpf_ptr f, mp_size_t size, mp_limb_t value)
{
ASSERT (size >= 0);
size = MIN (PREC(f) + 1, size);
SIZ(f) = size;
EXP(f) = size;
refmpn_fill (PTR(f), size, value);
}
void
refmpf_normalize (mpf_ptr f)
{
while (SIZ(f) != 0 && PTR(f)[ABSIZ(f)-1] == 0)
{
SIZ(f) = (SIZ(f) >= 0 ? SIZ(f)-1 : SIZ(f)+1);
EXP(f) --;
}
if (SIZ(f) == 0)
EXP(f) = 0;
}
unsigned long
refmpf_set_overlap (mpf_ptr dst, mpf_srcptr src)
{
mp_size_t dprec = PREC(dst);
mp_size_t ssize = ABSIZ(src);
unsigned long ret;
refmpf_set_prec_limbs (dst, (unsigned long) MAX (dprec, ssize));
mpf_set (dst, src);
ret = mpf_get_prec (dst);
PREC(dst) = dprec;
return ret;
}
void
refmpf_set_prec_limbs (mpf_ptr f, unsigned long prec)
{
mpf_set_prec (f, __GMPF_PREC_TO_BITS (prec));
}
void
refmpf_sub (mpf_ptr w, mpf_srcptr u, mpf_srcptr v)
{
mp_size_t hi, lo, size;
mp_ptr ut, vt, wt;
int neg;
mp_exp_t exp;
TMP_DECL;
TMP_MARK;
if (SIZ (u) == 0)
{
size = ABSIZ (v);
wt = TMP_ALLOC_LIMBS (size + 1);
MPN_COPY (wt, PTR (v), size);
exp = EXP (v);
neg = SIZ (v) > 0;
goto done;
}
if (SIZ (v) == 0)
{
size = ABSIZ (u);
wt = TMP_ALLOC_LIMBS (size + 1);
MPN_COPY (wt, PTR (u), size);
exp = EXP (u);
neg = SIZ (u) < 0;
goto done;
}
if ((SIZ (u) ^ SIZ (v)) < 0)
{
mpf_t tmp;
SIZ (tmp) = -SIZ (v);
EXP (tmp) = EXP (v);
PTR (tmp) = PTR (v);
refmpf_add (w, u, tmp);
if (SIZ (u) < 0)
mpf_neg (w, w);
return;
}
neg = SIZ (u) < 0;
hi = MAX (EXP (u), EXP (v));
lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v));
size = hi - lo;
ut = TMP_ALLOC_LIMBS (size + 1);
vt = TMP_ALLOC_LIMBS (size + 1);
wt = TMP_ALLOC_LIMBS (size + 1);
MPN_ZERO (ut, size);
MPN_ZERO (vt, size);
{int off;
off = size + (EXP (u) - hi) - ABSIZ (u);
MPN_COPY (ut + off, PTR (u), ABSIZ (u));
off = size + (EXP (v) - hi) - ABSIZ (v);
MPN_COPY (vt + off, PTR (v), ABSIZ (v));
}
if (mpn_cmp (ut, vt, size) >= 0)
mpn_sub_n (wt, ut, vt, size);
else
{
mpn_sub_n (wt, vt, ut, size);
neg ^= 1;
}
exp = hi;
while (size != 0 && wt[size - 1] == 0)
{
size--;
exp--;
}
done:
if (size > PREC (w))
{
wt += size - PREC (w);
size = PREC (w);
}
MPN_COPY (PTR (w), wt, size);
SIZ (w) = neg == 0 ? size : -size;
EXP (w) = exp;
TMP_FREE;
}
int
refmpf_validate (const char *name, mpf_srcptr got, mpf_srcptr want)
{
int bad = 0;
mp_size_t gsize, wsize, cmpsize, i;
mp_srcptr gp, wp;
mp_limb_t glimb, wlimb;
MPF_CHECK_FORMAT (got);
if (EXP (got) != EXP (want))
{
printf ("%s: wrong exponent\n", name);
bad = 1;
}
gsize = SIZ (got);
wsize = SIZ (want);
if ((gsize < 0 && wsize > 0) || (gsize > 0 && wsize < 0))
{
printf ("%s: wrong sign\n", name);
bad = 1;
}
gsize = ABS (gsize);
wsize = ABS (wsize);
gp = PTR (got) + gsize - 1;
wp = PTR (want) + wsize - 1;
cmpsize = MAX (PREC (got), gsize);
for (i = 0; i < cmpsize; i++)
{
glimb = (i < gsize ? gp[-i] : 0);
wlimb = (i < wsize ? wp[-i] : 0);
if (glimb != wlimb)
{
printf ("%s: wrong data starting at index %ld from top\n",
name, (long) i);
bad = 1;
break;
}
}
if (bad)
{
printf (" prec %d\n", PREC(got));
printf (" exp got %ld\n", (long) EXP(got));
printf (" exp want %ld\n", (long) EXP(want));
printf (" size got %d\n", SIZ(got));
printf (" size want %d\n", SIZ(want));
printf (" limbs (high to low)\n");
printf (" got ");
for (i = ABSIZ(got)-1; i >= 0; i--)
{
gmp_printf ("%MX", PTR(got)[i]);
if (i != 0)
printf (",");
}
printf ("\n");
printf (" want ");
for (i = ABSIZ(want)-1; i >= 0; i--)
{
gmp_printf ("%MX", PTR(want)[i]);
if (i != 0)
printf (",");
}
printf ("\n");
return 0;
}
return 1;
}
int
refmpf_validate_division (const char *name, mpf_srcptr got,
mpf_srcptr n, mpf_srcptr d)
{
mp_size_t nsize, dsize, sign, prec, qsize, tsize;
mp_srcptr np, dp;
mp_ptr tp, qp, rp;
mpf_t want;
int ret;
nsize = SIZ (n);
dsize = SIZ (d);
ASSERT_ALWAYS (dsize != 0);
sign = nsize ^ dsize;
nsize = ABS (nsize);
dsize = ABS (dsize);
np = PTR (n);
dp = PTR (d);
prec = PREC (got);
EXP (want) = EXP (n) - EXP (d) + 1;
qsize = prec + 2;
tsize = qsize + dsize - 1;
tp = refmpn_malloc_limbs (tsize);
refmpn_copy_extend (tp, tsize, np, nsize);
qp = refmpn_malloc_limbs (qsize);
rp = refmpn_malloc_limbs (dsize);
ASSERT_ALWAYS (qsize == tsize - dsize + 1);
refmpn_tdiv_qr (qp, rp, (mp_size_t) 0, tp, tsize, dp, dsize);
PTR (want) = qp;
SIZ (want) = (sign >= 0 ? qsize : -qsize);
refmpf_normalize (want);
ret = refmpf_validate (name, got, want);
free (tp);
free (qp);
free (rp);
return ret;
}