Use generic mpa.c code for everything except __mul and __sqr
This commit is contained in:
Siddhesh Poyarekar
parent
adbb8027be
commit
82a9811d29
|
|
@ -1,5 +1,12 @@
|
|||
2013-03-07 Siddhesh Poyarekar <siddhesh@redhat.com>
|
||||
|
||||
* sysdeps/ieee754/dbl-64/mpa.c [!NO__MUL]: Define __mul.
|
||||
[!NO__SQR]: Define __sqr.
|
||||
* sysdeps/powerpc/powerpc32/power4/fpu/mpa.c: define NO__MUL
|
||||
and NO__SQR. Remove all code except __mul and __sqr. Include
|
||||
sysdeps/ieee754/dbl-64/mpa.c.
|
||||
* sysdeps/powerpc/powerpc64/power4/fpu/mpa.c: Likewise.
|
||||
|
||||
[BZ #12723]
|
||||
* posix/Makefile (tests): Add tst-pathconf.
|
||||
* posix/tst-pathconf.c: New test case.
|
||||
|
|
|
|||
|
|
@ -611,6 +611,7 @@ __sub (const mp_no *x, const mp_no *y, mp_no *z, int p)
|
|||
}
|
||||
}
|
||||
|
||||
#ifndef NO__MUL
|
||||
/* Multiply *X and *Y and store result in *Z. X and Y may overlap but not X
|
||||
and Z or Y and Z. For P in [1, 2, 3], the exact result is truncated to P
|
||||
digits. In case P > 3 the error is bounded by 1.001 ULP. */
|
||||
|
|
@ -761,7 +762,9 @@ __mul (const mp_no *x, const mp_no *y, mp_no *z, int p)
|
|||
EZ = e;
|
||||
Z[0] = X[0] * Y[0];
|
||||
}
|
||||
#endif
|
||||
|
||||
#ifndef NO__SQR
|
||||
/* Square *X and store result in *Y. X and Y may not overlap. For P in
|
||||
[1, 2, 3], the exact result is truncated to P digits. In case P > 3 the
|
||||
error is bounded by 1.001 ULP. This is a faster special case of
|
||||
|
|
@ -862,6 +865,7 @@ __sqr (const mp_no *x, mp_no *y, int p)
|
|||
|
||||
EY = e;
|
||||
}
|
||||
#endif
|
||||
|
||||
/* Invert *X and store in *Y. Relative error bound:
|
||||
- For P = 2: 1.001 * R ^ (1 - P)
|
||||
|
|
|
|||
|
|
@ -17,582 +17,12 @@
|
|||
* You should have received a copy of the GNU Lesser General Public License
|
||||
* along with this program; if not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
/************************************************************************/
|
||||
/* MODULE_NAME: mpa.c */
|
||||
/* */
|
||||
/* FUNCTIONS: */
|
||||
/* mcr */
|
||||
/* acr */
|
||||
/* cpy */
|
||||
/* norm */
|
||||
/* denorm */
|
||||
/* mp_dbl */
|
||||
/* dbl_mp */
|
||||
/* add_magnitudes */
|
||||
/* sub_magnitudes */
|
||||
/* add */
|
||||
/* sub */
|
||||
/* mul */
|
||||
/* inv */
|
||||
/* dvd */
|
||||
/* */
|
||||
/* Arithmetic functions for multiple precision numbers. */
|
||||
/* Relative errors are bounded */
|
||||
/************************************************************************/
|
||||
|
||||
/* Define __mul and __sqr and use the rest from generic code. */
|
||||
#define NO__MUL
|
||||
#define NO__SQR
|
||||
|
||||
#include "endian.h"
|
||||
#include "mpa.h"
|
||||
#include <sys/param.h>
|
||||
|
||||
const mp_no mpone = {1, {1.0, 1.0}};
|
||||
const mp_no mptwo = {1, {1.0, 2.0}};
|
||||
|
||||
/* Compare mantissa of two multiple precision numbers regardless of the sign
|
||||
and exponent of the numbers. */
|
||||
static int
|
||||
mcr (const mp_no *x, const mp_no *y, int p)
|
||||
{
|
||||
long i;
|
||||
long p2 = p;
|
||||
for (i = 1; i <= p2; i++)
|
||||
{
|
||||
if (X[i] == Y[i])
|
||||
continue;
|
||||
else if (X[i] > Y[i])
|
||||
return 1;
|
||||
else
|
||||
return -1;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Compare the absolute values of two multiple precision numbers. */
|
||||
int
|
||||
__acr (const mp_no *x, const mp_no *y, int p)
|
||||
{
|
||||
long i;
|
||||
|
||||
if (X[0] == ZERO)
|
||||
{
|
||||
if (Y[0] == ZERO)
|
||||
i = 0;
|
||||
else
|
||||
i = -1;
|
||||
}
|
||||
else if (Y[0] == ZERO)
|
||||
i = 1;
|
||||
else
|
||||
{
|
||||
if (EX > EY)
|
||||
i = 1;
|
||||
else if (EX < EY)
|
||||
i = -1;
|
||||
else
|
||||
i = mcr (x, y, p);
|
||||
}
|
||||
|
||||
return i;
|
||||
}
|
||||
|
||||
/* Copy multiple precision number X into Y. They could be the same
|
||||
number. */
|
||||
void
|
||||
__cpy (const mp_no *x, mp_no *y, int p)
|
||||
{
|
||||
long i;
|
||||
|
||||
EY = EX;
|
||||
for (i = 0; i <= p; i++)
|
||||
Y[i] = X[i];
|
||||
}
|
||||
|
||||
/* Convert a multiple precision number *X into a double precision
|
||||
number *Y, normalized case (|x| >= 2**(-1022))). */
|
||||
static void
|
||||
norm (const mp_no *x, double *y, int p)
|
||||
{
|
||||
#define R RADIXI
|
||||
long i;
|
||||
double a, c, u, v, z[5];
|
||||
if (p < 5)
|
||||
{
|
||||
if (p == 1)
|
||||
c = X[1];
|
||||
else if (p == 2)
|
||||
c = X[1] + R * X[2];
|
||||
else if (p == 3)
|
||||
c = X[1] + R * (X[2] + R * X[3]);
|
||||
else if (p == 4)
|
||||
c = (X[1] + R * X[2]) + R * R * (X[3] + R * X[4]);
|
||||
}
|
||||
else
|
||||
{
|
||||
for (a = ONE, z[1] = X[1]; z[1] < TWO23;)
|
||||
{
|
||||
a *= TWO;
|
||||
z[1] *= TWO;
|
||||
}
|
||||
|
||||
for (i = 2; i < 5; i++)
|
||||
{
|
||||
z[i] = X[i] * a;
|
||||
u = (z[i] + CUTTER) - CUTTER;
|
||||
if (u > z[i])
|
||||
u -= RADIX;
|
||||
z[i] -= u;
|
||||
z[i - 1] += u * RADIXI;
|
||||
}
|
||||
|
||||
u = (z[3] + TWO71) - TWO71;
|
||||
if (u > z[3])
|
||||
u -= TWO19;
|
||||
v = z[3] - u;
|
||||
|
||||
if (v == TWO18)
|
||||
{
|
||||
if (z[4] == ZERO)
|
||||
{
|
||||
for (i = 5; i <= p; i++)
|
||||
{
|
||||
if (X[i] == ZERO)
|
||||
continue;
|
||||
else
|
||||
{
|
||||
z[3] += ONE;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
z[3] += ONE;
|
||||
}
|
||||
|
||||
c = (z[1] + R * (z[2] + R * z[3])) / a;
|
||||
}
|
||||
|
||||
c *= X[0];
|
||||
|
||||
for (i = 1; i < EX; i++)
|
||||
c *= RADIX;
|
||||
for (i = 1; i > EX; i--)
|
||||
c *= RADIXI;
|
||||
|
||||
*y = c;
|
||||
#undef R
|
||||
}
|
||||
|
||||
/* Convert a multiple precision number *X into a double precision
|
||||
number *Y, Denormal case (|x| < 2**(-1022))). */
|
||||
static void
|
||||
denorm (const mp_no *x, double *y, int p)
|
||||
{
|
||||
long i, k;
|
||||
long p2 = p;
|
||||
double c, u, z[5];
|
||||
|
||||
#define R RADIXI
|
||||
if (EX < -44 || (EX == -44 && X[1] < TWO5))
|
||||
{
|
||||
*y = ZERO;
|
||||
return;
|
||||
}
|
||||
|
||||
if (p2 == 1)
|
||||
{
|
||||
if (EX == -42)
|
||||
{
|
||||
z[1] = X[1] + TWO10;
|
||||
z[2] = ZERO;
|
||||
z[3] = ZERO;
|
||||
k = 3;
|
||||
}
|
||||
else if (EX == -43)
|
||||
{
|
||||
z[1] = TWO10;
|
||||
z[2] = X[1];
|
||||
z[3] = ZERO;
|
||||
k = 2;
|
||||
}
|
||||
else
|
||||
{
|
||||
z[1] = TWO10;
|
||||
z[2] = ZERO;
|
||||
z[3] = X[1];
|
||||
k = 1;
|
||||
}
|
||||
}
|
||||
else if (p2 == 2)
|
||||
{
|
||||
if (EX == -42)
|
||||
{
|
||||
z[1] = X[1] + TWO10;
|
||||
z[2] = X[2];
|
||||
z[3] = ZERO;
|
||||
k = 3;
|
||||
}
|
||||
else if (EX == -43)
|
||||
{
|
||||
z[1] = TWO10;
|
||||
z[2] = X[1];
|
||||
z[3] = X[2];
|
||||
k = 2;
|
||||
}
|
||||
else
|
||||
{
|
||||
z[1] = TWO10;
|
||||
z[2] = ZERO;
|
||||
z[3] = X[1];
|
||||
k = 1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (EX == -42)
|
||||
{
|
||||
z[1] = X[1] + TWO10;
|
||||
z[2] = X[2];
|
||||
k = 3;
|
||||
}
|
||||
else if (EX == -43)
|
||||
{
|
||||
z[1] = TWO10;
|
||||
z[2] = X[1];
|
||||
k = 2;
|
||||
}
|
||||
else
|
||||
{
|
||||
z[1] = TWO10;
|
||||
z[2] = ZERO;
|
||||
k = 1;
|
||||
}
|
||||
z[3] = X[k];
|
||||
}
|
||||
|
||||
u = (z[3] + TWO57) - TWO57;
|
||||
if (u > z[3])
|
||||
u -= TWO5;
|
||||
|
||||
if (u == z[3])
|
||||
{
|
||||
for (i = k + 1; i <= p2; i++)
|
||||
{
|
||||
if (X[i] == ZERO)
|
||||
continue;
|
||||
else
|
||||
{
|
||||
z[3] += ONE;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
c = X[0] * ((z[1] + R * (z[2] + R * z[3])) - TWO10);
|
||||
|
||||
*y = c * TWOM1032;
|
||||
#undef R
|
||||
}
|
||||
|
||||
/* Convert multiple precision number *X into double precision number *Y. The
|
||||
result is correctly rounded to the nearest/even. */
|
||||
void
|
||||
__mp_dbl (const mp_no *x, double *y, int p)
|
||||
{
|
||||
if (X[0] == ZERO)
|
||||
{
|
||||
*y = ZERO;
|
||||
return;
|
||||
}
|
||||
|
||||
if (__glibc_likely (EX > -42 || (EX == -42 && X[1] >= TWO10)))
|
||||
norm (x, y, p);
|
||||
else
|
||||
denorm (x, y, p);
|
||||
}
|
||||
|
||||
/* Get the multiple precision equivalent of X into *Y. If the precision is too
|
||||
small, the result is truncated. */
|
||||
void
|
||||
__dbl_mp (double x, mp_no *y, int p)
|
||||
{
|
||||
long i, n;
|
||||
long p2 = p;
|
||||
double u;
|
||||
|
||||
/* Sign. */
|
||||
if (x == ZERO)
|
||||
{
|
||||
Y[0] = ZERO;
|
||||
return;
|
||||
}
|
||||
else if (x > ZERO)
|
||||
Y[0] = ONE;
|
||||
else
|
||||
{
|
||||
Y[0] = MONE;
|
||||
x = -x;
|
||||
}
|
||||
|
||||
/* Exponent. */
|
||||
for (EY = ONE; x >= RADIX; EY += ONE)
|
||||
x *= RADIXI;
|
||||
for (; x < ONE; EY -= ONE)
|
||||
x *= RADIX;
|
||||
|
||||
/* Digits. */
|
||||
n = MIN (p2, 4);
|
||||
for (i = 1; i <= n; i++)
|
||||
{
|
||||
u = (x + TWO52) - TWO52;
|
||||
if (u > x)
|
||||
u -= ONE;
|
||||
Y[i] = u;
|
||||
x -= u;
|
||||
x *= RADIX;
|
||||
}
|
||||
for (; i <= p2; i++)
|
||||
Y[i] = ZERO;
|
||||
}
|
||||
|
||||
/* Add magnitudes of *X and *Y assuming that abs (*X) >= abs (*Y) > 0. The
|
||||
sign of the sum *Z is not changed. X and Y may overlap but not X and Z or
|
||||
Y and Z. No guard digit is used. The result equals the exact sum,
|
||||
truncated. */
|
||||
static void
|
||||
add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
|
||||
{
|
||||
long i, j, k;
|
||||
long p2 = p;
|
||||
double zk;
|
||||
|
||||
EZ = EX;
|
||||
|
||||
i = p2;
|
||||
j = p2 + EY - EX;
|
||||
k = p2 + 1;
|
||||
|
||||
if (__glibc_unlikely (j < 1))
|
||||
{
|
||||
__cpy (x, z, p);
|
||||
return;
|
||||
}
|
||||
|
||||
zk = ZERO;
|
||||
|
||||
for (; j > 0; i--, j--)
|
||||
{
|
||||
zk += X[i] + Y[j];
|
||||
if (zk >= RADIX)
|
||||
{
|
||||
Z[k--] = zk - RADIX;
|
||||
zk = ONE;
|
||||
}
|
||||
else
|
||||
{
|
||||
Z[k--] = zk;
|
||||
zk = ZERO;
|
||||
}
|
||||
}
|
||||
|
||||
for (; i > 0; i--)
|
||||
{
|
||||
zk += X[i];
|
||||
if (zk >= RADIX)
|
||||
{
|
||||
Z[k--] = zk - RADIX;
|
||||
zk = ONE;
|
||||
}
|
||||
else
|
||||
{
|
||||
Z[k--] = zk;
|
||||
zk = ZERO;
|
||||
}
|
||||
}
|
||||
|
||||
if (zk == ZERO)
|
||||
{
|
||||
for (i = 1; i <= p2; i++)
|
||||
Z[i] = Z[i + 1];
|
||||
}
|
||||
else
|
||||
{
|
||||
Z[1] = zk;
|
||||
EZ += ONE;
|
||||
}
|
||||
}
|
||||
|
||||
/* Subtract the magnitudes of *X and *Y assuming that abs (*x) > abs (*y) > 0.
|
||||
The sign of the difference *Z is not changed. X and Y may overlap but not X
|
||||
and Z or Y and Z. One guard digit is used. The error is less than one
|
||||
ULP. */
|
||||
static void
|
||||
sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
|
||||
{
|
||||
long i, j, k;
|
||||
long p2 = p;
|
||||
double zk;
|
||||
|
||||
EZ = EX;
|
||||
i = p2;
|
||||
j = p2 + EY - EX;
|
||||
k = p2;
|
||||
|
||||
/* Y is too small compared to X, copy X over to the result. */
|
||||
if (__glibc_unlikely (j < 1))
|
||||
{
|
||||
__cpy (x, z, p);
|
||||
return;
|
||||
}
|
||||
|
||||
/* The relevant least significant digit in Y is non-zero, so we factor it in
|
||||
to enhance accuracy. */
|
||||
if (j < p2 && Y[j + 1] > ZERO)
|
||||
{
|
||||
Z[k + 1] = RADIX - Y[j + 1];
|
||||
zk = MONE;
|
||||
}
|
||||
else
|
||||
zk = Z[k + 1] = ZERO;
|
||||
|
||||
/* Subtract and borrow. */
|
||||
for (; j > 0; i--, j--)
|
||||
{
|
||||
zk += (X[i] - Y[j]);
|
||||
if (zk < ZERO)
|
||||
{
|
||||
Z[k--] = zk + RADIX;
|
||||
zk = MONE;
|
||||
}
|
||||
else
|
||||
{
|
||||
Z[k--] = zk;
|
||||
zk = ZERO;
|
||||
}
|
||||
}
|
||||
|
||||
/* We're done with digits from Y, so it's just digits in X. */
|
||||
for (; i > 0; i--)
|
||||
{
|
||||
zk += X[i];
|
||||
if (zk < ZERO)
|
||||
{
|
||||
Z[k--] = zk + RADIX;
|
||||
zk = MONE;
|
||||
}
|
||||
else
|
||||
{
|
||||
Z[k--] = zk;
|
||||
zk = ZERO;
|
||||
}
|
||||
}
|
||||
|
||||
/* Normalize. */
|
||||
for (i = 1; Z[i] == ZERO; i++);
|
||||
EZ = EZ - i + 1;
|
||||
for (k = 1; i <= p2 + 1;)
|
||||
Z[k++] = Z[i++];
|
||||
for (; k <= p2;)
|
||||
Z[k++] = ZERO;
|
||||
}
|
||||
|
||||
/* Add *X and *Y and store the result in *Z. X and Y may overlap, but not X
|
||||
and Z or Y and Z. One guard digit is used. The error is less than one
|
||||
ULP. */
|
||||
void
|
||||
__add (const mp_no *x, const mp_no *y, mp_no *z, int p)
|
||||
{
|
||||
int n;
|
||||
|
||||
if (X[0] == ZERO)
|
||||
{
|
||||
__cpy (y, z, p);
|
||||
return;
|
||||
}
|
||||
else if (Y[0] == ZERO)
|
||||
{
|
||||
__cpy (x, z, p);
|
||||
return;
|
||||
}
|
||||
|
||||
if (X[0] == Y[0])
|
||||
{
|
||||
if (__acr (x, y, p) > 0)
|
||||
{
|
||||
add_magnitudes (x, y, z, p);
|
||||
Z[0] = X[0];
|
||||
}
|
||||
else
|
||||
{
|
||||
add_magnitudes (y, x, z, p);
|
||||
Z[0] = Y[0];
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if ((n = __acr (x, y, p)) == 1)
|
||||
{
|
||||
sub_magnitudes (x, y, z, p);
|
||||
Z[0] = X[0];
|
||||
}
|
||||
else if (n == -1)
|
||||
{
|
||||
sub_magnitudes (y, x, z, p);
|
||||
Z[0] = Y[0];
|
||||
}
|
||||
else
|
||||
Z[0] = ZERO;
|
||||
}
|
||||
}
|
||||
|
||||
/* Subtract *Y from *X and return the result in *Z. X and Y may overlap but
|
||||
not X and Z or Y and Z. One guard digit is used. The error is less than
|
||||
one ULP. */
|
||||
void
|
||||
__sub (const mp_no *x, const mp_no *y, mp_no *z, int p)
|
||||
{
|
||||
int n;
|
||||
|
||||
if (X[0] == ZERO)
|
||||
{
|
||||
__cpy (y, z, p);
|
||||
Z[0] = -Z[0];
|
||||
return;
|
||||
}
|
||||
else if (Y[0] == ZERO)
|
||||
{
|
||||
__cpy (x, z, p);
|
||||
return;
|
||||
}
|
||||
|
||||
if (X[0] != Y[0])
|
||||
{
|
||||
if (__acr (x, y, p) > 0)
|
||||
{
|
||||
add_magnitudes (x, y, z, p);
|
||||
Z[0] = X[0];
|
||||
}
|
||||
else
|
||||
{
|
||||
add_magnitudes (y, x, z, p);
|
||||
Z[0] = -Y[0];
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if ((n = __acr (x, y, p)) == 1)
|
||||
{
|
||||
sub_magnitudes (x, y, z, p);
|
||||
Z[0] = X[0];
|
||||
}
|
||||
else if (n == -1)
|
||||
{
|
||||
sub_magnitudes (y, x, z, p);
|
||||
Z[0] = -Y[0];
|
||||
}
|
||||
else
|
||||
Z[0] = ZERO;
|
||||
}
|
||||
}
|
||||
#include <sysdeps/ieee754/dbl-64/mpa.c>
|
||||
|
||||
/* Multiply *X and *Y and store result in *Z. X and Y may overlap but not X
|
||||
and Z or Y and Z. For P in [1, 2, 3], the exact result is truncated to P
|
||||
|
|
@ -781,57 +211,3 @@ __sqr (const mp_no *x, mp_no *y, int p)
|
|||
EY--;
|
||||
}
|
||||
}
|
||||
|
||||
/* Invert *X and store in *Y. Relative error bound:
|
||||
- For P = 2: 1.001 * R ^ (1 - P)
|
||||
- For P = 3: 1.063 * R ^ (1 - P)
|
||||
- For P > 3: 2.001 * R ^ (1 - P)
|
||||
|
||||
*X = 0 is not permissible. */
|
||||
static void
|
||||
__inv (const mp_no *x, mp_no *y, int p)
|
||||
{
|
||||
long i;
|
||||
double t;
|
||||
mp_no z, w;
|
||||
static const int np1[] =
|
||||
{ 0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3,
|
||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4
|
||||
};
|
||||
|
||||
__cpy (x, &z, p);
|
||||
z.e = 0;
|
||||
__mp_dbl (&z, &t, p);
|
||||
t = ONE / t;
|
||||
__dbl_mp (t, y, p);
|
||||
EY -= EX;
|
||||
|
||||
for (i = 0; i < np1[p]; i++)
|
||||
{
|
||||
__cpy (y, &w, p);
|
||||
__mul (x, &w, y, p);
|
||||
__sub (&mptwo, y, &z, p);
|
||||
__mul (&w, &z, y, p);
|
||||
}
|
||||
}
|
||||
|
||||
/* Divide *X by *Y and store result in *Z. X and Y may overlap but not X and Z
|
||||
or Y and Z. Relative error bound:
|
||||
- For P = 2: 2.001 * R ^ (1 - P)
|
||||
- For P = 3: 2.063 * R ^ (1 - P)
|
||||
- For P > 3: 3.001 * R ^ (1 - P)
|
||||
|
||||
*X = 0 is not permissible. */
|
||||
void
|
||||
__dvd (const mp_no *x, const mp_no *y, mp_no *z, int p)
|
||||
{
|
||||
mp_no w;
|
||||
|
||||
if (X[0] == ZERO)
|
||||
Z[0] = ZERO;
|
||||
else
|
||||
{
|
||||
__inv (y, &w, p);
|
||||
__mul (x, &w, z, p);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -17,582 +17,12 @@
|
|||
* You should have received a copy of the GNU Lesser General Public License
|
||||
* along with this program; if not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
/************************************************************************/
|
||||
/* MODULE_NAME: mpa.c */
|
||||
/* */
|
||||
/* FUNCTIONS: */
|
||||
/* mcr */
|
||||
/* acr */
|
||||
/* cpy */
|
||||
/* norm */
|
||||
/* denorm */
|
||||
/* mp_dbl */
|
||||
/* dbl_mp */
|
||||
/* add_magnitudes */
|
||||
/* sub_magnitudes */
|
||||
/* add */
|
||||
/* sub */
|
||||
/* mul */
|
||||
/* inv */
|
||||
/* dvd */
|
||||
/* */
|
||||
/* Arithmetic functions for multiple precision numbers. */
|
||||
/* Relative errors are bounded */
|
||||
/************************************************************************/
|
||||
|
||||
/* Define __mul and __sqr and use the rest from generic code. */
|
||||
#define NO__MUL
|
||||
#define NO__SQR
|
||||
|
||||
#include "endian.h"
|
||||
#include "mpa.h"
|
||||
#include <sys/param.h>
|
||||
|
||||
const mp_no mpone = {1, {1.0, 1.0}};
|
||||
const mp_no mptwo = {1, {1.0, 2.0}};
|
||||
|
||||
/* Compare mantissa of two multiple precision numbers regardless of the sign
|
||||
and exponent of the numbers. */
|
||||
static int
|
||||
mcr (const mp_no *x, const mp_no *y, int p)
|
||||
{
|
||||
long i;
|
||||
long p2 = p;
|
||||
for (i = 1; i <= p2; i++)
|
||||
{
|
||||
if (X[i] == Y[i])
|
||||
continue;
|
||||
else if (X[i] > Y[i])
|
||||
return 1;
|
||||
else
|
||||
return -1;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Compare the absolute values of two multiple precision numbers. */
|
||||
int
|
||||
__acr (const mp_no *x, const mp_no *y, int p)
|
||||
{
|
||||
long i;
|
||||
|
||||
if (X[0] == ZERO)
|
||||
{
|
||||
if (Y[0] == ZERO)
|
||||
i = 0;
|
||||
else
|
||||
i = -1;
|
||||
}
|
||||
else if (Y[0] == ZERO)
|
||||
i = 1;
|
||||
else
|
||||
{
|
||||
if (EX > EY)
|
||||
i = 1;
|
||||
else if (EX < EY)
|
||||
i = -1;
|
||||
else
|
||||
i = mcr (x, y, p);
|
||||
}
|
||||
|
||||
return i;
|
||||
}
|
||||
|
||||
/* Copy multiple precision number X into Y. They could be the same
|
||||
number. */
|
||||
void
|
||||
__cpy (const mp_no *x, mp_no *y, int p)
|
||||
{
|
||||
long i;
|
||||
|
||||
EY = EX;
|
||||
for (i = 0; i <= p; i++)
|
||||
Y[i] = X[i];
|
||||
}
|
||||
|
||||
/* Convert a multiple precision number *X into a double precision
|
||||
number *Y, normalized case (|x| >= 2**(-1022))). */
|
||||
static void
|
||||
norm (const mp_no *x, double *y, int p)
|
||||
{
|
||||
#define R RADIXI
|
||||
long i;
|
||||
double a, c, u, v, z[5];
|
||||
if (p < 5)
|
||||
{
|
||||
if (p == 1)
|
||||
c = X[1];
|
||||
else if (p == 2)
|
||||
c = X[1] + R * X[2];
|
||||
else if (p == 3)
|
||||
c = X[1] + R * (X[2] + R * X[3]);
|
||||
else if (p == 4)
|
||||
c = (X[1] + R * X[2]) + R * R * (X[3] + R * X[4]);
|
||||
}
|
||||
else
|
||||
{
|
||||
for (a = ONE, z[1] = X[1]; z[1] < TWO23;)
|
||||
{
|
||||
a *= TWO;
|
||||
z[1] *= TWO;
|
||||
}
|
||||
|
||||
for (i = 2; i < 5; i++)
|
||||
{
|
||||
z[i] = X[i] * a;
|
||||
u = (z[i] + CUTTER) - CUTTER;
|
||||
if (u > z[i])
|
||||
u -= RADIX;
|
||||
z[i] -= u;
|
||||
z[i - 1] += u * RADIXI;
|
||||
}
|
||||
|
||||
u = (z[3] + TWO71) - TWO71;
|
||||
if (u > z[3])
|
||||
u -= TWO19;
|
||||
v = z[3] - u;
|
||||
|
||||
if (v == TWO18)
|
||||
{
|
||||
if (z[4] == ZERO)
|
||||
{
|
||||
for (i = 5; i <= p; i++)
|
||||
{
|
||||
if (X[i] == ZERO)
|
||||
continue;
|
||||
else
|
||||
{
|
||||
z[3] += ONE;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
z[3] += ONE;
|
||||
}
|
||||
|
||||
c = (z[1] + R * (z[2] + R * z[3])) / a;
|
||||
}
|
||||
|
||||
c *= X[0];
|
||||
|
||||
for (i = 1; i < EX; i++)
|
||||
c *= RADIX;
|
||||
for (i = 1; i > EX; i--)
|
||||
c *= RADIXI;
|
||||
|
||||
*y = c;
|
||||
#undef R
|
||||
}
|
||||
|
||||
/* Convert a multiple precision number *X into a double precision
|
||||
number *Y, Denormal case (|x| < 2**(-1022))). */
|
||||
static void
|
||||
denorm (const mp_no *x, double *y, int p)
|
||||
{
|
||||
long i, k;
|
||||
long p2 = p;
|
||||
double c, u, z[5];
|
||||
|
||||
#define R RADIXI
|
||||
if (EX < -44 || (EX == -44 && X[1] < TWO5))
|
||||
{
|
||||
*y = ZERO;
|
||||
return;
|
||||
}
|
||||
|
||||
if (p2 == 1)
|
||||
{
|
||||
if (EX == -42)
|
||||
{
|
||||
z[1] = X[1] + TWO10;
|
||||
z[2] = ZERO;
|
||||
z[3] = ZERO;
|
||||
k = 3;
|
||||
}
|
||||
else if (EX == -43)
|
||||
{
|
||||
z[1] = TWO10;
|
||||
z[2] = X[1];
|
||||
z[3] = ZERO;
|
||||
k = 2;
|
||||
}
|
||||
else
|
||||
{
|
||||
z[1] = TWO10;
|
||||
z[2] = ZERO;
|
||||
z[3] = X[1];
|
||||
k = 1;
|
||||
}
|
||||
}
|
||||
else if (p2 == 2)
|
||||
{
|
||||
if (EX == -42)
|
||||
{
|
||||
z[1] = X[1] + TWO10;
|
||||
z[2] = X[2];
|
||||
z[3] = ZERO;
|
||||
k = 3;
|
||||
}
|
||||
else if (EX == -43)
|
||||
{
|
||||
z[1] = TWO10;
|
||||
z[2] = X[1];
|
||||
z[3] = X[2];
|
||||
k = 2;
|
||||
}
|
||||
else
|
||||
{
|
||||
z[1] = TWO10;
|
||||
z[2] = ZERO;
|
||||
z[3] = X[1];
|
||||
k = 1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (EX == -42)
|
||||
{
|
||||
z[1] = X[1] + TWO10;
|
||||
z[2] = X[2];
|
||||
k = 3;
|
||||
}
|
||||
else if (EX == -43)
|
||||
{
|
||||
z[1] = TWO10;
|
||||
z[2] = X[1];
|
||||
k = 2;
|
||||
}
|
||||
else
|
||||
{
|
||||
z[1] = TWO10;
|
||||
z[2] = ZERO;
|
||||
k = 1;
|
||||
}
|
||||
z[3] = X[k];
|
||||
}
|
||||
|
||||
u = (z[3] + TWO57) - TWO57;
|
||||
if (u > z[3])
|
||||
u -= TWO5;
|
||||
|
||||
if (u == z[3])
|
||||
{
|
||||
for (i = k + 1; i <= p2; i++)
|
||||
{
|
||||
if (X[i] == ZERO)
|
||||
continue;
|
||||
else
|
||||
{
|
||||
z[3] += ONE;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
c = X[0] * ((z[1] + R * (z[2] + R * z[3])) - TWO10);
|
||||
|
||||
*y = c * TWOM1032;
|
||||
#undef R
|
||||
}
|
||||
|
||||
/* Convert multiple precision number *X into double precision number *Y. The
|
||||
result is correctly rounded to the nearest/even. */
|
||||
void
|
||||
__mp_dbl (const mp_no *x, double *y, int p)
|
||||
{
|
||||
if (X[0] == ZERO)
|
||||
{
|
||||
*y = ZERO;
|
||||
return;
|
||||
}
|
||||
|
||||
if (__glibc_likely (EX > -42 || (EX == -42 && X[1] >= TWO10)))
|
||||
norm (x, y, p);
|
||||
else
|
||||
denorm (x, y, p);
|
||||
}
|
||||
|
||||
/* Get the multiple precision equivalent of X into *Y. If the precision is too
|
||||
small, the result is truncated. */
|
||||
void
|
||||
__dbl_mp (double x, mp_no *y, int p)
|
||||
{
|
||||
long i, n;
|
||||
long p2 = p;
|
||||
double u;
|
||||
|
||||
/* Sign. */
|
||||
if (x == ZERO)
|
||||
{
|
||||
Y[0] = ZERO;
|
||||
return;
|
||||
}
|
||||
else if (x > ZERO)
|
||||
Y[0] = ONE;
|
||||
else
|
||||
{
|
||||
Y[0] = MONE;
|
||||
x = -x;
|
||||
}
|
||||
|
||||
/* Exponent. */
|
||||
for (EY = ONE; x >= RADIX; EY += ONE)
|
||||
x *= RADIXI;
|
||||
for (; x < ONE; EY -= ONE)
|
||||
x *= RADIX;
|
||||
|
||||
/* Digits. */
|
||||
n = MIN (p2, 4);
|
||||
for (i = 1; i <= n; i++)
|
||||
{
|
||||
u = (x + TWO52) - TWO52;
|
||||
if (u > x)
|
||||
u -= ONE;
|
||||
Y[i] = u;
|
||||
x -= u;
|
||||
x *= RADIX;
|
||||
}
|
||||
for (; i <= p2; i++)
|
||||
Y[i] = ZERO;
|
||||
}
|
||||
|
||||
/* Add magnitudes of *X and *Y assuming that abs (*X) >= abs (*Y) > 0. The
|
||||
sign of the sum *Z is not changed. X and Y may overlap but not X and Z or
|
||||
Y and Z. No guard digit is used. The result equals the exact sum,
|
||||
truncated. */
|
||||
static void
|
||||
add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
|
||||
{
|
||||
long i, j, k;
|
||||
long p2 = p;
|
||||
double zk;
|
||||
|
||||
EZ = EX;
|
||||
|
||||
i = p2;
|
||||
j = p2 + EY - EX;
|
||||
k = p2 + 1;
|
||||
|
||||
if (__glibc_unlikely (j < 1))
|
||||
{
|
||||
__cpy (x, z, p);
|
||||
return;
|
||||
}
|
||||
|
||||
zk = ZERO;
|
||||
|
||||
for (; j > 0; i--, j--)
|
||||
{
|
||||
zk += X[i] + Y[j];
|
||||
if (zk >= RADIX)
|
||||
{
|
||||
Z[k--] = zk - RADIX;
|
||||
zk = ONE;
|
||||
}
|
||||
else
|
||||
{
|
||||
Z[k--] = zk;
|
||||
zk = ZERO;
|
||||
}
|
||||
}
|
||||
|
||||
for (; i > 0; i--)
|
||||
{
|
||||
zk += X[i];
|
||||
if (zk >= RADIX)
|
||||
{
|
||||
Z[k--] = zk - RADIX;
|
||||
zk = ONE;
|
||||
}
|
||||
else
|
||||
{
|
||||
Z[k--] = zk;
|
||||
zk = ZERO;
|
||||
}
|
||||
}
|
||||
|
||||
if (zk == ZERO)
|
||||
{
|
||||
for (i = 1; i <= p2; i++)
|
||||
Z[i] = Z[i + 1];
|
||||
}
|
||||
else
|
||||
{
|
||||
Z[1] = zk;
|
||||
EZ += ONE;
|
||||
}
|
||||
}
|
||||
|
||||
/* Subtract the magnitudes of *X and *Y assuming that abs (*x) > abs (*y) > 0.
|
||||
The sign of the difference *Z is not changed. X and Y may overlap but not X
|
||||
and Z or Y and Z. One guard digit is used. The error is less than one
|
||||
ULP. */
|
||||
static void
|
||||
sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p)
|
||||
{
|
||||
long i, j, k;
|
||||
long p2 = p;
|
||||
double zk;
|
||||
|
||||
EZ = EX;
|
||||
i = p2;
|
||||
j = p2 + EY - EX;
|
||||
k = p2;
|
||||
|
||||
/* Y is too small compared to X, copy X over to the result. */
|
||||
if (__glibc_unlikely (j < 1))
|
||||
{
|
||||
__cpy (x, z, p);
|
||||
return;
|
||||
}
|
||||
|
||||
/* The relevant least significant digit in Y is non-zero, so we factor it in
|
||||
to enhance accuracy. */
|
||||
if (j < p2 && Y[j + 1] > ZERO)
|
||||
{
|
||||
Z[k + 1] = RADIX - Y[j + 1];
|
||||
zk = MONE;
|
||||
}
|
||||
else
|
||||
zk = Z[k + 1] = ZERO;
|
||||
|
||||
/* Subtract and borrow. */
|
||||
for (; j > 0; i--, j--)
|
||||
{
|
||||
zk += (X[i] - Y[j]);
|
||||
if (zk < ZERO)
|
||||
{
|
||||
Z[k--] = zk + RADIX;
|
||||
zk = MONE;
|
||||
}
|
||||
else
|
||||
{
|
||||
Z[k--] = zk;
|
||||
zk = ZERO;
|
||||
}
|
||||
}
|
||||
|
||||
/* We're done with digits from Y, so it's just digits in X. */
|
||||
for (; i > 0; i--)
|
||||
{
|
||||
zk += X[i];
|
||||
if (zk < ZERO)
|
||||
{
|
||||
Z[k--] = zk + RADIX;
|
||||
zk = MONE;
|
||||
}
|
||||
else
|
||||
{
|
||||
Z[k--] = zk;
|
||||
zk = ZERO;
|
||||
}
|
||||
}
|
||||
|
||||
/* Normalize. */
|
||||
for (i = 1; Z[i] == ZERO; i++);
|
||||
EZ = EZ - i + 1;
|
||||
for (k = 1; i <= p2 + 1;)
|
||||
Z[k++] = Z[i++];
|
||||
for (; k <= p2;)
|
||||
Z[k++] = ZERO;
|
||||
}
|
||||
|
||||
/* Add *X and *Y and store the result in *Z. X and Y may overlap, but not X
|
||||
and Z or Y and Z. One guard digit is used. The error is less than one
|
||||
ULP. */
|
||||
void
|
||||
__add (const mp_no *x, const mp_no *y, mp_no *z, int p)
|
||||
{
|
||||
int n;
|
||||
|
||||
if (X[0] == ZERO)
|
||||
{
|
||||
__cpy (y, z, p);
|
||||
return;
|
||||
}
|
||||
else if (Y[0] == ZERO)
|
||||
{
|
||||
__cpy (x, z, p);
|
||||
return;
|
||||
}
|
||||
|
||||
if (X[0] == Y[0])
|
||||
{
|
||||
if (__acr (x, y, p) > 0)
|
||||
{
|
||||
add_magnitudes (x, y, z, p);
|
||||
Z[0] = X[0];
|
||||
}
|
||||
else
|
||||
{
|
||||
add_magnitudes (y, x, z, p);
|
||||
Z[0] = Y[0];
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if ((n = __acr (x, y, p)) == 1)
|
||||
{
|
||||
sub_magnitudes (x, y, z, p);
|
||||
Z[0] = X[0];
|
||||
}
|
||||
else if (n == -1)
|
||||
{
|
||||
sub_magnitudes (y, x, z, p);
|
||||
Z[0] = Y[0];
|
||||
}
|
||||
else
|
||||
Z[0] = ZERO;
|
||||
}
|
||||
}
|
||||
|
||||
/* Subtract *Y from *X and return the result in *Z. X and Y may overlap but
|
||||
not X and Z or Y and Z. One guard digit is used. The error is less than
|
||||
one ULP. */
|
||||
void
|
||||
__sub (const mp_no *x, const mp_no *y, mp_no *z, int p)
|
||||
{
|
||||
int n;
|
||||
|
||||
if (X[0] == ZERO)
|
||||
{
|
||||
__cpy (y, z, p);
|
||||
Z[0] = -Z[0];
|
||||
return;
|
||||
}
|
||||
else if (Y[0] == ZERO)
|
||||
{
|
||||
__cpy (x, z, p);
|
||||
return;
|
||||
}
|
||||
|
||||
if (X[0] != Y[0])
|
||||
{
|
||||
if (__acr (x, y, p) > 0)
|
||||
{
|
||||
add_magnitudes (x, y, z, p);
|
||||
Z[0] = X[0];
|
||||
}
|
||||
else
|
||||
{
|
||||
add_magnitudes (y, x, z, p);
|
||||
Z[0] = -Y[0];
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if ((n = __acr (x, y, p)) == 1)
|
||||
{
|
||||
sub_magnitudes (x, y, z, p);
|
||||
Z[0] = X[0];
|
||||
}
|
||||
else if (n == -1)
|
||||
{
|
||||
sub_magnitudes (y, x, z, p);
|
||||
Z[0] = -Y[0];
|
||||
}
|
||||
else
|
||||
Z[0] = ZERO;
|
||||
}
|
||||
}
|
||||
#include <sysdeps/ieee754/dbl-64/mpa.c>
|
||||
|
||||
/* Multiply *X and *Y and store result in *Z. X and Y may overlap but not X
|
||||
and Z or Y and Z. For P in [1, 2, 3], the exact result is truncated to P
|
||||
|
|
@ -781,57 +211,3 @@ __sqr (const mp_no *x, mp_no *y, int p)
|
|||
EY--;
|
||||
}
|
||||
}
|
||||
|
||||
/* Invert *X and store in *Y. Relative error bound:
|
||||
- For P = 2: 1.001 * R ^ (1 - P)
|
||||
- For P = 3: 1.063 * R ^ (1 - P)
|
||||
- For P > 3: 2.001 * R ^ (1 - P)
|
||||
|
||||
*X = 0 is not permissible. */
|
||||
static void
|
||||
__inv (const mp_no *x, mp_no *y, int p)
|
||||
{
|
||||
long i;
|
||||
double t;
|
||||
mp_no z, w;
|
||||
static const int np1[] =
|
||||
{ 0, 0, 0, 0, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3,
|
||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4
|
||||
};
|
||||
|
||||
__cpy (x, &z, p);
|
||||
z.e = 0;
|
||||
__mp_dbl (&z, &t, p);
|
||||
t = ONE / t;
|
||||
__dbl_mp (t, y, p);
|
||||
EY -= EX;
|
||||
|
||||
for (i = 0; i < np1[p]; i++)
|
||||
{
|
||||
__cpy (y, &w, p);
|
||||
__mul (x, &w, y, p);
|
||||
__sub (&mptwo, y, &z, p);
|
||||
__mul (&w, &z, y, p);
|
||||
}
|
||||
}
|
||||
|
||||
/* Divide *X by *Y and store result in *Z. X and Y may overlap but not X and Z
|
||||
or Y and Z. Relative error bound:
|
||||
- For P = 2: 2.001 * R ^ (1 - P)
|
||||
- For P = 3: 2.063 * R ^ (1 - P)
|
||||
- For P > 3: 3.001 * R ^ (1 - P)
|
||||
|
||||
*X = 0 is not permissible. */
|
||||
void
|
||||
__dvd (const mp_no *x, const mp_no *y, mp_no *z, int p)
|
||||
{
|
||||
mp_no w;
|
||||
|
||||
if (X[0] == ZERO)
|
||||
Z[0] = ZERO;
|
||||
else
|
||||
{
|
||||
__inv (y, &w, p);
|
||||
__mul (x, &w, z, p);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue