|
Packit |
5c3484 |
/* mpf_sub -- Subtract two floats.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Copyright 1993-1996, 1999-2002, 2004, 2005, 2011, 2014 Free Software Foundation, Inc.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
This file is part of the GNU MP Library.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The GNU MP Library is free software; you can redistribute it and/or modify
|
|
Packit |
5c3484 |
it under the terms of either:
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
* the GNU Lesser General Public License as published by the Free
|
|
Packit |
5c3484 |
Software Foundation; either version 3 of the License, or (at your
|
|
Packit |
5c3484 |
option) any later version.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
or
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
* the GNU General Public License as published by the Free Software
|
|
Packit |
5c3484 |
Foundation; either version 2 of the License, or (at your option) any
|
|
Packit |
5c3484 |
later version.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
or both in parallel, as here.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The GNU MP Library is distributed in the hope that it will be useful, but
|
|
Packit |
5c3484 |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|
Packit |
5c3484 |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
Packit |
5c3484 |
for more details.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
You should have received copies of the GNU General Public License and the
|
|
Packit |
5c3484 |
GNU Lesser General Public License along with the GNU MP Library. If not,
|
|
Packit |
5c3484 |
see https://www.gnu.org/licenses/. */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#include "gmp.h"
|
|
Packit |
5c3484 |
#include "gmp-impl.h"
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
void
|
|
Packit |
5c3484 |
mpf_sub (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mp_srcptr up, vp;
|
|
Packit |
5c3484 |
mp_ptr rp, tp;
|
|
Packit |
5c3484 |
mp_size_t usize, vsize, rsize;
|
|
Packit |
5c3484 |
mp_size_t prec;
|
|
Packit |
5c3484 |
mp_exp_t exp;
|
|
Packit |
5c3484 |
mp_size_t ediff;
|
|
Packit |
5c3484 |
int negate;
|
|
Packit |
5c3484 |
TMP_DECL;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
usize = SIZ (u);
|
|
Packit |
5c3484 |
vsize = SIZ (v);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Handle special cases that don't work in generic code below. */
|
|
Packit |
5c3484 |
if (usize == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mpf_neg (r, v);
|
|
Packit |
5c3484 |
return;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
if (vsize == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
if (r != u)
|
|
Packit |
5c3484 |
mpf_set (r, u);
|
|
Packit |
5c3484 |
return;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* If signs of U and V are different, perform addition. */
|
|
Packit |
5c3484 |
if ((usize ^ vsize) < 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
__mpf_struct v_negated;
|
|
Packit |
5c3484 |
v_negated._mp_size = -vsize;
|
|
Packit |
5c3484 |
v_negated._mp_exp = EXP (v);
|
|
Packit |
5c3484 |
v_negated._mp_d = PTR (v);
|
|
Packit |
5c3484 |
mpf_add (r, u, &v_negated);
|
|
Packit |
5c3484 |
return;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
TMP_MARK;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Signs are now known to be the same. */
|
|
Packit |
5c3484 |
negate = usize < 0;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Make U be the operand with the largest exponent. */
|
|
Packit |
5c3484 |
if (EXP (u) < EXP (v))
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mpf_srcptr t;
|
|
Packit |
5c3484 |
t = u; u = v; v = t;
|
|
Packit |
5c3484 |
negate ^= 1;
|
|
Packit |
5c3484 |
usize = SIZ (u);
|
|
Packit |
5c3484 |
vsize = SIZ (v);
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
usize = ABS (usize);
|
|
Packit |
5c3484 |
vsize = ABS (vsize);
|
|
Packit |
5c3484 |
up = PTR (u);
|
|
Packit |
5c3484 |
vp = PTR (v);
|
|
Packit |
5c3484 |
rp = PTR (r);
|
|
Packit |
5c3484 |
prec = PREC (r) + 1;
|
|
Packit |
5c3484 |
exp = EXP (u);
|
|
Packit |
5c3484 |
ediff = exp - EXP (v);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* If ediff is 0 or 1, we might have a situation where the operands are
|
|
Packit |
5c3484 |
extremely close. We need to scan the operands from the most significant
|
|
Packit |
5c3484 |
end ignore the initial parts that are equal. */
|
|
Packit |
5c3484 |
if (ediff <= 1)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
if (ediff == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* Skip leading limbs in U and V that are equal. */
|
|
Packit |
5c3484 |
/* This loop normally exits immediately. Optimize for that. */
|
|
Packit |
5c3484 |
while (up[usize - 1] == vp[vsize - 1])
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
usize--;
|
|
Packit |
5c3484 |
vsize--;
|
|
Packit |
5c3484 |
exp--;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (usize == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* u cancels high limbs of v, result is rest of v */
|
|
Packit |
5c3484 |
negate ^= 1;
|
|
Packit |
5c3484 |
cancellation:
|
|
Packit |
5c3484 |
/* strip high zeros before truncating to prec */
|
|
Packit |
5c3484 |
while (vsize != 0 && vp[vsize - 1] == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
vsize--;
|
|
Packit |
5c3484 |
exp--;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
if (vsize > prec)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
vp += vsize - prec;
|
|
Packit |
5c3484 |
vsize = prec;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
MPN_COPY_INCR (rp, vp, vsize);
|
|
Packit |
5c3484 |
rsize = vsize;
|
|
Packit |
5c3484 |
goto done;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
if (vsize == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
vp = up;
|
|
Packit |
5c3484 |
vsize = usize;
|
|
Packit |
5c3484 |
goto cancellation;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (up[usize - 1] < vp[vsize - 1])
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* For simplicity, swap U and V. Note that since the loop above
|
|
Packit |
5c3484 |
wouldn't have exited unless up[usize - 1] and vp[vsize - 1]
|
|
Packit |
5c3484 |
were non-equal, this if-statement catches all cases where U
|
|
Packit |
5c3484 |
is smaller than V. */
|
|
Packit |
5c3484 |
MPN_SRCPTR_SWAP (up,usize, vp,vsize);
|
|
Packit |
5c3484 |
negate ^= 1;
|
|
Packit |
5c3484 |
/* negating ediff not necessary since it is 0. */
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Check for
|
|
Packit |
5c3484 |
x+1 00000000 ...
|
|
Packit |
5c3484 |
x ffffffff ... */
|
|
Packit |
5c3484 |
if (up[usize - 1] != vp[vsize - 1] + 1)
|
|
Packit |
5c3484 |
goto general_case;
|
|
Packit |
5c3484 |
usize--;
|
|
Packit |
5c3484 |
vsize--;
|
|
Packit |
5c3484 |
exp--;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else /* ediff == 1 */
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* Check for
|
|
Packit |
5c3484 |
1 00000000 ...
|
|
Packit |
5c3484 |
0 ffffffff ... */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (up[usize - 1] != 1 || vp[vsize - 1] != GMP_NUMB_MAX
|
|
Packit |
5c3484 |
|| (usize >= 2 && up[usize - 2] != 0))
|
|
Packit |
5c3484 |
goto general_case;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
usize--;
|
|
Packit |
5c3484 |
exp--;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Skip sequences of 00000000/ffffffff */
|
|
Packit |
5c3484 |
while (vsize != 0 && usize != 0 && up[usize - 1] == 0
|
|
Packit |
5c3484 |
&& vp[vsize - 1] == GMP_NUMB_MAX)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
usize--;
|
|
Packit |
5c3484 |
vsize--;
|
|
Packit |
5c3484 |
exp--;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (usize == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
while (vsize != 0 && vp[vsize - 1] == GMP_NUMB_MAX)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
vsize--;
|
|
Packit |
5c3484 |
exp--;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else if (usize > prec - 1)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
up += usize - (prec - 1);
|
|
Packit |
5c3484 |
usize = prec - 1;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
if (vsize > prec - 1)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
vp += vsize - (prec - 1);
|
|
Packit |
5c3484 |
vsize = prec - 1;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
tp = TMP_ALLOC_LIMBS (prec);
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mp_limb_t cy_limb;
|
|
Packit |
5c3484 |
if (vsize == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
MPN_COPY (tp, up, usize);
|
|
Packit |
5c3484 |
tp[usize] = 1;
|
|
Packit |
5c3484 |
rsize = usize + 1;
|
|
Packit |
5c3484 |
exp++;
|
|
Packit |
5c3484 |
goto normalized;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
if (usize == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
cy_limb = mpn_neg (tp, vp, vsize);
|
|
Packit |
5c3484 |
rsize = vsize;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else if (usize >= vsize)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* uuuu */
|
|
Packit |
5c3484 |
/* vv */
|
|
Packit |
5c3484 |
mp_size_t size;
|
|
Packit |
5c3484 |
size = usize - vsize;
|
|
Packit |
5c3484 |
MPN_COPY (tp, up, size);
|
|
Packit |
5c3484 |
cy_limb = mpn_sub_n (tp + size, up + size, vp, vsize);
|
|
Packit |
5c3484 |
rsize = usize;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else /* (usize < vsize) */
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* uuuu */
|
|
Packit |
5c3484 |
/* vvvvvvv */
|
|
Packit |
5c3484 |
mp_size_t size;
|
|
Packit |
5c3484 |
size = vsize - usize;
|
|
Packit |
5c3484 |
cy_limb = mpn_neg (tp, vp, size);
|
|
Packit |
5c3484 |
cy_limb = mpn_sub_nc (tp + size, up, vp + size, usize, cy_limb);
|
|
Packit |
5c3484 |
rsize = vsize;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
if (cy_limb == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
tp[rsize] = 1;
|
|
Packit |
5c3484 |
rsize++;
|
|
Packit |
5c3484 |
exp++;
|
|
Packit |
5c3484 |
goto normalized;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
goto normalize;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
general_case:
|
|
Packit |
5c3484 |
/* If U extends beyond PREC, ignore the part that does. */
|
|
Packit |
5c3484 |
if (usize > prec)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
up += usize - prec;
|
|
Packit |
5c3484 |
usize = prec;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* If V extends beyond PREC, ignore the part that does.
|
|
Packit |
5c3484 |
Note that this may make vsize negative. */
|
|
Packit |
5c3484 |
if (vsize + ediff > prec)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
vp += vsize + ediff - prec;
|
|
Packit |
5c3484 |
vsize = prec - ediff;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (ediff >= prec)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* V completely cancelled. */
|
|
Packit |
5c3484 |
if (rp != up)
|
|
Packit |
5c3484 |
MPN_COPY (rp, up, usize);
|
|
Packit |
5c3484 |
rsize = usize;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* Allocate temp space for the result. Allocate
|
|
Packit |
5c3484 |
just vsize + ediff later??? */
|
|
Packit |
5c3484 |
tp = TMP_ALLOC_LIMBS (prec);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Locate the least significant non-zero limb in (the needed
|
|
Packit |
5c3484 |
parts of) U and V, to simplify the code below. */
|
|
Packit |
5c3484 |
for (;;)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
if (vsize == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
MPN_COPY (rp, up, usize);
|
|
Packit |
5c3484 |
rsize = usize;
|
|
Packit |
5c3484 |
goto done;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
if (vp[0] != 0)
|
|
Packit |
5c3484 |
break;
|
|
Packit |
5c3484 |
vp++, vsize--;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
for (;;)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
if (usize == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
MPN_COPY (rp, vp, vsize);
|
|
Packit |
5c3484 |
rsize = vsize;
|
|
Packit |
5c3484 |
negate ^= 1;
|
|
Packit |
5c3484 |
goto done;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
if (up[0] != 0)
|
|
Packit |
5c3484 |
break;
|
|
Packit |
5c3484 |
up++, usize--;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* uuuu | uuuu | uuuu | uuuu | uuuu */
|
|
Packit |
5c3484 |
/* vvvvvvv | vv | vvvvv | v | vv */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (usize > ediff)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* U and V partially overlaps. */
|
|
Packit |
5c3484 |
if (ediff == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* Have to compare the leading limbs of u and v
|
|
Packit |
5c3484 |
to determine whether to compute u - v or v - u. */
|
|
Packit |
5c3484 |
if (usize >= vsize)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* uuuu */
|
|
Packit |
5c3484 |
/* vv */
|
|
Packit |
5c3484 |
mp_size_t size;
|
|
Packit |
5c3484 |
size = usize - vsize;
|
|
Packit |
5c3484 |
MPN_COPY (tp, up, size);
|
|
Packit |
5c3484 |
mpn_sub_n (tp + size, up + size, vp, vsize);
|
|
Packit |
5c3484 |
rsize = usize;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else /* (usize < vsize) */
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* uuuu */
|
|
Packit |
5c3484 |
/* vvvvvvv */
|
|
Packit |
5c3484 |
mp_size_t size;
|
|
Packit |
5c3484 |
size = vsize - usize;
|
|
Packit |
5c3484 |
ASSERT_CARRY (mpn_neg (tp, vp, size));
|
|
Packit |
5c3484 |
mpn_sub_nc (tp + size, up, vp + size, usize, CNST_LIMB (1));
|
|
Packit |
5c3484 |
rsize = vsize;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
if (vsize + ediff <= usize)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* uuuu */
|
|
Packit |
5c3484 |
/* v */
|
|
Packit |
5c3484 |
mp_size_t size;
|
|
Packit |
5c3484 |
size = usize - ediff - vsize;
|
|
Packit |
5c3484 |
MPN_COPY (tp, up, size);
|
|
Packit |
5c3484 |
mpn_sub (tp + size, up + size, usize - size, vp, vsize);
|
|
Packit |
5c3484 |
rsize = usize;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* uuuu */
|
|
Packit |
5c3484 |
/* vvvvv */
|
|
Packit |
5c3484 |
mp_size_t size;
|
|
Packit |
5c3484 |
rsize = vsize + ediff;
|
|
Packit |
5c3484 |
size = rsize - usize;
|
|
Packit |
5c3484 |
ASSERT_CARRY (mpn_neg (tp, vp, size));
|
|
Packit |
5c3484 |
mpn_sub (tp + size, up, usize, vp + size, usize - ediff);
|
|
Packit |
5c3484 |
/* Should we use sub_nc then sub_1? */
|
|
Packit |
5c3484 |
MPN_DECR_U (tp + size, usize, CNST_LIMB (1));
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
/* uuuu */
|
|
Packit |
5c3484 |
/* vv */
|
|
Packit |
5c3484 |
mp_size_t size, i;
|
|
Packit |
5c3484 |
size = vsize + ediff - usize;
|
|
Packit |
5c3484 |
ASSERT_CARRY (mpn_neg (tp, vp, vsize));
|
|
Packit |
5c3484 |
for (i = vsize; i < size; i++)
|
|
Packit |
5c3484 |
tp[i] = GMP_NUMB_MAX;
|
|
Packit |
5c3484 |
mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1);
|
|
Packit |
5c3484 |
rsize = size + usize;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
normalize:
|
|
Packit |
5c3484 |
/* Full normalize. Optimize later. */
|
|
Packit |
5c3484 |
while (rsize != 0 && tp[rsize - 1] == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
rsize--;
|
|
Packit |
5c3484 |
exp--;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
normalized:
|
|
Packit |
5c3484 |
MPN_COPY (rp, tp, rsize);
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
done:
|
|
Packit |
5c3484 |
TMP_FREE;
|
|
Packit |
5c3484 |
if (rsize == 0) {
|
|
Packit |
5c3484 |
SIZ (r) = 0;
|
|
Packit |
5c3484 |
EXP (r) = 0;
|
|
Packit |
5c3484 |
} else {
|
|
Packit |
5c3484 |
SIZ (r) = negate ? -rsize : rsize;
|
|
Packit |
5c3484 |
EXP (r) = exp;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|