diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpa.c | 691 |
1 files changed, 463 insertions, 228 deletions
diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c index 7abad67823..b3bfa6c0fe 100644 --- a/sysdeps/ieee754/dbl-64/mpa.c +++ b/sysdeps/ieee754/dbl-64/mpa.c @@ -60,30 +60,45 @@ 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) { +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; - else if (X[i] > Y[i]) return 1; - else return -1; } + for (i = 1; i <= p; 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) { +__acr (const mp_no *x, const mp_no *y, int p) +{ int 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); - } + 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; } @@ -92,59 +107,86 @@ __acr(const mp_no *x, const mp_no *y, int p) { #ifndef NO___CPY /* Copy multiple precision number X into Y. They could be the same number. */ -void __cpy(const mp_no *x, mp_no *y, int p) { +void +__cpy (const mp_no *x, mp_no *y, int p) +{ EY = EX; - for (int i=0; i <= p; i++) Y[i] = X[i]; + for (int i = 0; i <= p; i++) + Y[i] = X[i]; } #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) +static void +norm (const mp_no *x, double *y, int p) { - #define R RADIXI +#define R RADIXI int 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; + 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; + } - c = (z[1] + R *(z[2] + R * z[3]))/a; - } + 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; + for (i = 1; i < EX; i++) + c *= RADIX; + for (i = 1; i > EX; i--) + c *= RADIXI; *y = c; #undef R @@ -152,58 +194,129 @@ static void norm(const mp_no *x, double *y, int p) /* 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) +static void +denorm (const mp_no *x, double *y, int p) { - int i,k; - double c,u,z[5]; + int i, k; + double c, u, z[5]; #define R RADIXI - if (EX<-44 || (EX==-44 && X[1]<TWO5)) - { *y=ZERO; return; } - - if (p==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 (p==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]; - } + if (EX < -44 || (EX == -44 && X[1] < TWO5)) + { + *y = ZERO; + return; + } + + if (p == 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 (p == 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]) + u -= TWO5; - if (u==z[3]) { - for (i=k+1; i <= p; i++) { - if (X[i] == ZERO) continue; - else {z[3] += ONE; break; } + if (u == z[3]) + { + for (i = k + 1; i <= p; i++) + { + if (X[i] == ZERO) + continue; + else + { + z[3] += ONE; + break; + } + } } - } - c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10); + c = X[0] * ((z[1] + R * (z[2] + R * z[3])) - TWO10); - *y = c*TWOM1032; + *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; } +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); + norm (x, y, p); else - denorm(x,y,p); + denorm (x, y, p); } #endif @@ -211,27 +324,44 @@ void __mp_dbl(const mp_no *x, double *y, int p) { small, the result is truncated. */ void SECTION -__dbl_mp(double x, mp_no *y, int p) { - - int i,n; +__dbl_mp (double x, mp_no *y, int p) +{ + int i, n; double u; /* Sign. */ - if (x == ZERO) {Y[0] = ZERO; return; } - else if (x > ZERO) Y[0] = ONE; - else {Y[0] = MONE; x=-x; } + 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; + for (EY = ONE; x >= RADIX; EY += ONE) + x *= RADIXI; + for (; x < ONE; EY -= ONE) + x *= RADIX; /* Digits. */ - n=MIN(p,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<=p; i++) Y[i] = ZERO; + n = MIN (p, 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 <= p; i++) + Y[i] = ZERO; } /* Add magnitudes of *X and *Y assuming that abs (*X) >= abs (*Y) > 0. The @@ -240,39 +370,55 @@ __dbl_mp(double x, mp_no *y, int p) { truncated. */ static void SECTION -add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { - - int i,j,k; +add_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p) +{ + int i, j, k; EZ = EX; - i=p; j=p+ EY - EX; k=p+1; - - if (j<1) - {__cpy(x,z,p); return; } - else Z[k] = ZERO; - - for (; j>0; i--,j--) { - Z[k] += X[i] + Y[j]; - if (Z[k] >= RADIX) { - Z[k] -= RADIX; - Z[--k] = ONE; } - else - Z[--k] = ZERO; - } - - for (; i>0; i--) { - Z[k] += X[i]; - if (Z[k] >= RADIX) { - Z[k] -= RADIX; - Z[--k] = ONE; } - else - Z[--k] = ZERO; - } - - if (Z[1] == ZERO) { - for (i=1; i<=p; i++) Z[i] = Z[i+1]; } - else EZ += ONE; + i = p; + j = p + EY - EX; + k = p + 1; + + if (j < 1) + { + __cpy (x, z, p); + return; + } + else + Z[k] = ZERO; + + for (; j > 0; i--, j--) + { + Z[k] += X[i] + Y[j]; + if (Z[k] >= RADIX) + { + Z[k] -= RADIX; + Z[--k] = ONE; + } + else + Z[--k] = ZERO; + } + + for (; i > 0; i--) + { + Z[k] += X[i]; + if (Z[k] >= RADIX) + { + Z[k] -= RADIX; + Z[--k] = ONE; + } + else + Z[--k] = ZERO; + } + + if (Z[1] == ZERO) + { + for (i = 1; i <= p; i++) + Z[i] = Z[i + 1]; + } + else + EZ += ONE; } /* Subtract the magnitudes of *X and *Y assuming that abs (*x) > abs (*y) > 0. @@ -281,52 +427,73 @@ add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { ULP. */ static void SECTION -sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { - - int i,j,k; +sub_magnitudes (const mp_no *x, const mp_no *y, mp_no *z, int p) +{ + int i, j, k; EZ = EX; - if (EX == EY) { - i=j=k=p; - Z[k] = Z[k+1] = ZERO; } - else { - j= EX - EY; - if (j > p) {__cpy(x,z,p); return; } - else { - i=p; j=p+1-j; k=p; - if (Y[j] > ZERO) { - Z[k+1] = RADIX - Y[j--]; - Z[k] = MONE; } - else { - Z[k+1] = ZERO; - Z[k] = ZERO; j--;} - } - } - - for (; j>0; i--,j--) { - Z[k] += (X[i] - Y[j]); - if (Z[k] < ZERO) { - Z[k] += RADIX; - Z[--k] = MONE; } - else - Z[--k] = ZERO; - } - - for (; i>0; i--) { - Z[k] += X[i]; - if (Z[k] < ZERO) { - Z[k] += RADIX; - Z[--k] = MONE; } - else - Z[--k] = ZERO; - } - - for (i=1; Z[i] == ZERO; i++) ; + if (EX == EY) + { + i = j = k = p; + Z[k] = Z[k + 1] = ZERO; + } + else + { + j = EX - EY; + if (j > p) + { + __cpy (x, z, p); + return; + } + else + { + i = p; + j = p + 1 - j; + k = p; + if (Y[j] > ZERO) + { + Z[k + 1] = RADIX - Y[j--]; + Z[k] = MONE; + } + else + { + Z[k + 1] = ZERO; + Z[k] = ZERO; + j--; + } + } + } + + for (; j > 0; i--, j--) + { + Z[k] += (X[i] - Y[j]); + if (Z[k] < ZERO) + { + Z[k] += RADIX; + Z[--k] = MONE; + } + else + Z[--k] = ZERO; + } + + for (; i > 0; i--) + { + Z[k] += X[i]; + if (Z[k] < ZERO) + { + Z[k] += RADIX; + Z[--k] = MONE; + } + else + Z[--k] = ZERO; + } + + for (i = 1; Z[i] == ZERO; i++); EZ = EZ - i + 1; - for (k=1; i <= p+1; ) + for (k = 1; i <= p + 1;) Z[k++] = Z[i++]; - for (; k <= p; ) + for (; k <= p;) Z[k++] = ZERO; } @@ -335,22 +502,49 @@ sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { ULP. */ void SECTION -__add(const mp_no *x, const mp_no *y, mp_no *z, int p) { - +__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; - } + 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 @@ -358,22 +552,50 @@ __add(const mp_no *x, const mp_no *y, mp_no *z, int p) { one ULP. */ void SECTION -__sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { - +__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; - } + 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; + } } /* Multiply *X and *Y and store result in *Z. X and Y may overlap but not X @@ -381,15 +603,15 @@ __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { digits. In case P > 3 the error is bounded by 1.001 ULP. */ void SECTION -__mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { - +__mul (const mp_no *x, const mp_no *y, mp_no *z, int p) +{ int i, j, k, k2; double u; /* Is z=0? */ if (__glibc_unlikely (X[0] * Y[0] == ZERO)) { - Z[0]=ZERO; + Z[0] = ZERO; return; } @@ -397,7 +619,7 @@ __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { k2 = (__glibc_unlikely (p < 3)) ? p + p : p + 3; Z[k2] = ZERO; - for (k = k2; k > p; ) + for (k = k2; k > p;) { for (i = k - p, j = p; i < p + 1; i++, j--) Z[k] += X[i] * Y[j]; @@ -411,7 +633,7 @@ __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { while (k > 1) { - for (i = 1,j = k - 1; i < k; i++, j--) + for (i = 1, j = k - 1; i < k; i++, j--) Z[k] += X[i] * Y[j]; u = (Z[k] + CUTTER) - CUTTER; @@ -426,7 +648,7 @@ __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { if (__glibc_unlikely (Z[1] == ZERO)) { for (i = 1; i <= p; i++) - Z[i] = Z[i+1]; + Z[i] = Z[i + 1]; EZ--; } @@ -439,24 +661,32 @@ __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { - For P > 3: 2.001 * R ^ (1 - P) *X = 0 is not permissible. */ -static +static void SECTION -void __inv(const mp_no *x, mp_no *y, int p) { +__inv (const mp_no *x, mp_no *y, int p) +{ int 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); - } + 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 @@ -468,10 +698,15 @@ void __inv(const mp_no *x, mp_no *y, int p) { *X = 0 is not permissible. */ void SECTION -__dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) { - +__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);} + if (X[0] == ZERO) + Z[0] = ZERO; + else + { + __inv (y, &w, p); + __mul (x, &w, z, p); + } } |