summaryrefslogtreecommitdiff
path: root/libc/sysdeps/ieee754/dbl-64/mpa.c
diff options
context:
space:
mode:
Diffstat (limited to 'libc/sysdeps/ieee754/dbl-64/mpa.c')
-rw-r--r--libc/sysdeps/ieee754/dbl-64/mpa.c113
1 files changed, 64 insertions, 49 deletions
diff --git a/libc/sysdeps/ieee754/dbl-64/mpa.c b/libc/sysdeps/ieee754/dbl-64/mpa.c
index 68647ba33..39c640882 100644
--- a/libc/sysdeps/ieee754/dbl-64/mpa.c
+++ b/libc/sysdeps/ieee754/dbl-64/mpa.c
@@ -1,8 +1,7 @@
-
/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
- * Copyright (C) 2001 Free Software Foundation
+ * Copyright (C) 2001, 2011 Free Software Foundation
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
@@ -48,11 +47,18 @@
#include "mpa.h"
#include "mpa2.h"
#include <sys/param.h> /* For MIN() */
+
+#ifndef SECTION
+# define SECTION
+#endif
+
+#ifndef NO___ACR
/* mcr() compares the sizes of the mantissas of two multiple precision */
/* numbers. Mantissas are compared regardless of the signs of the */
/* numbers, even if x->d[0] or y->d[0] are zero. Exponents are also */
/* disregarded. */
-static int mcr(const mp_no *x, const mp_no *y, int p) {
+static int
+mcr(const mp_no *x, const mp_no *y, int p) {
int i;
for (i=1; i<=p; i++) {
if (X[i] == Y[i]) continue;
@@ -62,9 +68,9 @@ static int mcr(const mp_no *x, const mp_no *y, int p) {
}
-
/* acr() compares the absolute values of two multiple precision numbers */
-int __acr(const mp_no *x, const mp_no *y, int p) {
+int
+__acr(const mp_no *x, const mp_no *y, int p) {
int i;
if (X[0] == ZERO) {
@@ -80,10 +86,12 @@ int __acr(const mp_no *x, const mp_no *y, int p) {
return i;
}
+#endif
-/* cr90 compares the values of two multiple precision numbers */
-int __cr(const mp_no *x, const mp_no *y, int p) {
+#if 0
+/* cr() compares the values of two multiple precision numbers */
+static int __cr(const mp_no *x, const mp_no *y, int p) {
int i;
if (X[0] > Y[0]) i= 1;
@@ -93,36 +101,37 @@ int __cr(const mp_no *x, const mp_no *y, int p) {
return i;
}
+#endif
+#ifndef NO___CPY
/* Copy a multiple precision number. Set *y=*x. x=y is permissible. */
void __cpy(const mp_no *x, mp_no *y, int p) {
- int i;
-
EY = EX;
- for (i=0; i <= p; i++) Y[i] = X[i];
-
- return;
+ for (int i=0; i <= p; i++) Y[i] = X[i];
}
+#endif
+#if 0
/* Copy a multiple precision number x of precision m into a */
/* multiple precision number y of precision n. In case n>m, */
/* the digits of y beyond the m'th are set to zero. In case */
/* n<m, the digits of x beyond the n'th are ignored. */
/* x=y is permissible. */
-void __cpymn(const mp_no *x, int m, mp_no *y, int n) {
+static void __cpymn(const mp_no *x, int m, mp_no *y, int n) {
int i,k;
EY = EX; k=MIN(m,n);
for (i=0; i <= k; i++) Y[i] = X[i];
for ( ; i <= n; i++) Y[i] = ZERO;
-
- return;
}
+#endif
+
+#ifndef NO___MP_DBL
/* 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)
@@ -141,7 +150,7 @@ static void norm(const mp_no *x, double *y, int p)
}
else {
for (a=ONE, z[1]=X[1]; z[1] < TWO23; )
- {a *= TWO; z[1] *= TWO; }
+ {a *= TWO; z[1] *= TWO; }
for (i=2; i<5; i++) {
z[i] = X[i]*a;
@@ -157,10 +166,10 @@ static void norm(const mp_no *x, double *y, int p)
if (v == TWO18) {
if (z[4] == ZERO) {
- for (i=5; i <= p; i++) {
- if (X[i] == ZERO) continue;
- else {z[3] += ONE; break; }
- }
+ for (i=5; i <= p; i++) {
+ if (X[i] == ZERO) continue;
+ else {z[3] += ONE; break; }
+ }
}
else z[3] += ONE;
}
@@ -174,7 +183,6 @@ static void norm(const mp_no *x, double *y, int p)
for (i=1; i>EX; i--) c *= RADIXI;
*y = c;
- return;
#undef R
}
@@ -222,8 +230,6 @@ static void denorm(const mp_no *x, double *y, int p)
c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10);
*y = c*TWOM1032;
- return;
-
#undef R
}
@@ -242,13 +248,16 @@ void __mp_dbl(const mp_no *x, double *y, int p) {
else if (EX==-42 && X[1]>=TWO10) norm(x,y,p);
else denorm(x,y,p);
}
+#endif
/* dbl_mp() converts a double precision number x into a multiple precision */
/* number *y. If the precision p is too small the result is truncated. x is */
/* left unchanged. */
-void __dbl_mp(double x, mp_no *y, int p) {
+void
+SECTION
+__dbl_mp(double x, mp_no *y, int p) {
int i,n;
double u;
@@ -269,7 +278,6 @@ void __dbl_mp(double x, mp_no *y, int p) {
if (u>x) u -= ONE;
Y[i] = u; x -= u; x *= RADIX; }
for ( ; i<=p; i++) Y[i] = ZERO;
- return;
}
@@ -279,7 +287,9 @@ void __dbl_mp(double x, mp_no *y, int p) {
/* No guard digit is used. The result equals the exact sum, truncated. */
/* *x & *y are left unchanged. */
-static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+static void
+SECTION
+add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int i,j,k;
@@ -321,7 +331,9 @@ static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* or y&z. One guard digit is used. The error is less than one ulp. */
/* *x & *y are left unchanged. */
-static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+static void
+SECTION
+sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int i,j,k;
@@ -336,11 +348,11 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
else {
i=p; j=p+1-j; k=p;
if (Y[j] > ZERO) {
- Z[k+1] = RADIX - Y[j--];
- Z[k] = MONE; }
+ Z[k+1] = RADIX - Y[j--];
+ Z[k] = MONE; }
else {
- Z[k+1] = ZERO;
- Z[k] = ZERO; j--;}
+ Z[k+1] = ZERO;
+ Z[k] = ZERO; j--;}
}
}
@@ -368,8 +380,6 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
Z[k++] = Z[i++];
for (; k <= p; )
Z[k++] = ZERO;
-
- return;
}
@@ -377,7 +387,9 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* but not x&z or y&z. One guard digit is used. The error is less than */
/* one ulp. *x & *y are left unchanged. */
-void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+void
+SECTION
+__add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int n;
@@ -393,7 +405,6 @@ void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = Y[0]; }
else Z[0] = ZERO;
}
- return;
}
@@ -401,7 +412,9 @@ void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* overlap but not x&z or y&z. One guard digit is used. The error is */
/* less than one ulp. *x & *y are left unchanged. */
-void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+void
+SECTION
+__sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int n;
@@ -417,7 +430,6 @@ void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = -Y[0]; }
else Z[0] = ZERO;
}
- return;
}
@@ -426,16 +438,18 @@ void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp. */
/* *x & *y are left unchanged. */
-void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+void
+SECTION
+__mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
int i, i1, i2, j, k, k2;
double u;
- /* Is z=0? */
+ /* Is z=0? */
if (X[0]*Y[0]==ZERO)
{ Z[0]=ZERO; return; }
- /* Multiply, add and carry */
+ /* Multiply, add and carry */
k2 = (p<3) ? p+p : p+3;
Z[k2]=ZERO;
for (k=k2; k>1; ) {
@@ -449,7 +463,7 @@ void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
Z[--k] = u*RADIXI;
}
- /* Is there a carry beyond the most significant digit? */
+ /* Is there a carry beyond the most significant digit? */
if (Z[1] == ZERO) {
for (i=1; i<=p; i++) Z[i]=Z[i+1];
EZ = EX + EY - 1; }
@@ -457,7 +471,6 @@ void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
EZ = EX + EY;
Z[0] = X[0] * Y[0];
- return;
}
@@ -466,6 +479,8 @@ void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
/* 2.001*r**(1-p) for p>3. */
/* *x=0 is not permissible. *x is left unchanged. */
+static
+SECTION
void __inv(const mp_no *x, mp_no *y, int p) {
int i;
#if 0
@@ -474,11 +489,11 @@ void __inv(const mp_no *x, mp_no *y, int p) {
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};
+ 4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};
const mp_no mptwo = {1,{1.0,2.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
- 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
+ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
+ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
+ 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
__cpy(x,&z,p); z.e=0; __mp_dbl(&z,&t,p);
t=ONE/t; __dbl_mp(t,y,p); EY -= EX;
@@ -489,7 +504,6 @@ void __inv(const mp_no *x, mp_no *y, int p) {
__sub(&mptwo,y,&z,p);
__mul(&w,&z,y,p);
}
- return;
}
@@ -498,11 +512,12 @@ void __inv(const mp_no *x, mp_no *y, int p) {
/* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3 */
/* and 3.001*r**(1-p) for p>3. *y=0 is not permissible. */
-void __dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
+void
+SECTION
+__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);}
- return;
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-23 00:30:33 +0000