diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/mpa.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpa.c | 43 |
1 files changed, 26 insertions, 17 deletions
diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c index cdd15a1..1586d91 100644 --- a/sysdeps/ieee754/dbl-64/mpa.c +++ b/sysdeps/ieee754/dbl-64/mpa.c @@ -5,9 +5,9 @@ * * 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 - * the Free Software Foundation; either version 2 of the License, or + * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. - * + * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the @@ -15,7 +15,7 @@ * * You should have received a copy of the GNU Lesser General Public License * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ /************************************************************************/ /* MODULE_NAME: mpa.c */ @@ -123,10 +123,13 @@ void cpymn(const mp_no *x, int m, mp_no *y, int n) { /* 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.d - int i,k; + int i; +#if 0 + int k; +#endif double a,c,u,v,z[5]; if (p<5) { if (p==1) c = X[1]; @@ -135,7 +138,7 @@ static void norm(const mp_no *x, double *y, int p) 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; ) + for (a=ONE, z[1]=X[1]; z[1] < TWO23; ) {a *= TWO; z[1] *= TWO; } for (i=2; i<5; i++) { @@ -175,10 +178,13 @@ static void norm(const mp_no *x, double *y, int p) /* Convert a multiple precision number *x into a double precision */ /* number *y, denormalized 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 a,c,u,v,z[5]; + double c,u,z[5]; +#if 0 + double a,v; +#endif #define R radixi.d if (EX<-44 || (EX==-44 && X[1]<TWO5)) @@ -222,10 +228,11 @@ static void denorm(const mp_no *x, double *y, int p) /* Convert a multiple precision number *x into a double precision number *y. */ /* The result is correctly rounded to the nearest/even. *x is left unchanged */ -void mp_dbl(const mp_no *x, double *y, int p) { - +void __mp_dbl(const mp_no *x, double *y, int p) { +#if 0 int i,k; double a,c,u,v,z[5]; +#endif if (X[0] == ZERO) {*y = ZERO; return; } @@ -458,17 +465,20 @@ void mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { /* *x=0 is not permissible. *x is left unchanged. */ void inv(const mp_no *x, mp_no *y, int p) { - int i,l; + int i; +#if 0 + int l; +#endif 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}; - 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, + 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}}; - cpy(x,&z,p); z.e=0; mp_dbl(&z,&t,p); + 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++) { @@ -494,4 +504,3 @@ void dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) { else {inv(y,&w,p); mul(x,&w,z,p);} return; } - |