diff options
25 files changed, 350 insertions, 303 deletions
diff --git a/ChangeLog b/ChangeLog index 964320481d..0741b0f34b 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,10 +1,26 @@ +2001-03-12 Ulrich Drepper <drepper@redhat.com> + + * sysdeps/ieee754/dbl-64/e_remainder.c: Fix handling of boundary + conditions. + + * sysdeps/ieee754/dbl-64/e_pow.c: Fix handling of boundary + conditions. + + * sysdeps/ieee754/dbl-64/s_sin.c (__sin): Handle Inf and NaN + correctly. + (__cos): Likewise. + + * sysdeps/ieee754/dbl-64/e_asin.c (__ieee754_asin): Handle NaN + correctly. + (__ieee754_acos): Likewise. + 2001-03-12 Andreas Jaeger <aj@suse.de> * sysdeps/unix/sysv/linux/s390/sysdep.h (_LINUX_S390_SYSDEP_H): Fix typo. Patch by Martin Schwidefsky <schwidefsky@de.ibm.com>. * sysdeps/s390/bits/string.h: Protect __STRING_INLINE against - redefinition. + redefinition. 2001-03-11 Roland McGrath <roland@frob.com> @@ -12,6 +28,12 @@ 2001-03-11 Ulrich Drepper <drepper@redhat.com> + * sysdeps/ieee754/dbl-64/endian.h: Define also one of BIG_ENDI and + LITTLE_ENDI. + + * sysdeps/ieee754/dbl-64/MathLib.h (Init_Lib): Use void as + parameter list. + Last-bit accurate math library implementation by IBM Haifa. Contributed by Abraham Ziv <ziv@il.ibm.com>, Moshe Olshansky <olshansk@il.ibm.com>, Ealan Henis <ealan@il.ibm.com>, and diff --git a/sysdeps/ieee754/dbl-64/branred.c b/sysdeps/ieee754/dbl-64/branred.c index 2bceb947c0..3253d34aba 100644 --- a/sysdeps/ieee754/dbl-64/branred.c +++ b/sysdeps/ieee754/dbl-64/branred.c @@ -45,7 +45,7 @@ /* x=n*pi/2+(a+aa), abs(a+aa)<pi/4, n=0,+-1,+-2,.... */ /* Routine return integer (n mod 4) */ /*******************************************************************/ -int branred(double x, double *a, double *aa) +int __branred(double x, double *a, double *aa) { int i,k; #if 0 diff --git a/sysdeps/ieee754/dbl-64/doasin.c b/sysdeps/ieee754/dbl-64/doasin.c index 4ea108b9d7..d3cc88bb19 100644 --- a/sysdeps/ieee754/dbl-64/doasin.c +++ b/sysdeps/ieee754/dbl-64/doasin.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,12 +15,12 @@ * * 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: doasin.c */ /* */ -/* FUNCTION: doasin */ +/* FUNCTION: doasin */ /* */ /* FILES NEEDED:endian.h mydefs.h dla.h doasin.h */ /* mpa.c */ @@ -31,17 +31,17 @@ #include "endian.h" #include "mydefs.h" -#include "dla.h" +#include "dla.h" /********************************************************************/ /* Compute arcsin(x,dx,v) of double-length number (x+dx) the result */ /* stored in v where v= v[0]+v[1] =arcsin(x+dx) */ /********************************************************************/ -void doasin(double x, double dx, double v[]) { - +void __doasin(double x, double dx, double v[]) { + #include "doasin.h" - static const double + static const double d5 = 0.22372159090911789889975459505194491E-01, d6 = 0.17352764422456822913014975683014622E-01, d7 = 0.13964843843786693521653681033981614E-01, @@ -49,10 +49,10 @@ void doasin(double x, double dx, double v[]) { d9 = 0.97622386568166960207425666787248914E-02, d10 = 0.83638737193775788576092749009744976E-02, d11 = 0.79470250400727425881446981833568758E-02; - + double xx,p,pp,u,uu,r,s; double hx,tx,hy,ty,tp,tq,tc,tcc; - + /* Taylor series for arcsin for Double-Length numbers */ xx = x*x+2.0*x*dx; @@ -73,9 +73,3 @@ void doasin(double x, double dx, double v[]) { v[0]=p; v[1]=pp; /* arcsin(x+dx)=v[0]+v[1] */ } - - - - - - diff --git a/sysdeps/ieee754/dbl-64/dosincos.c b/sysdeps/ieee754/dbl-64/dosincos.c index 692fc85889..498cc60929 100644 --- a/sysdeps/ieee754/dbl-64/dosincos.c +++ b/sysdeps/ieee754/dbl-64/dosincos.c @@ -45,7 +45,7 @@ /*(x+dx) between 0 and PI/4 */ /***********************************************************************/ -void dubsin(double x, double dx, double v[]) { +void __dubsin(double x, double dx, double v[]) { double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee, sn,ssn,cs,ccs,ds,dss,dc,dcc; #if 0 @@ -96,7 +96,7 @@ void dubsin(double x, double dx, double v[]) { /*(x+dx) between 0 and PI/4 */ /**********************************************************************/ -void dubcos(double x, double dx, double v[]) { +void __dubcos(double x, double dx, double v[]) { double r,s,p,hx,tx,hy,ty,q,c,cc,d,dd,d2,dd2,e,ee, sn,ssn,cs,ccs,ds,dss,dc,dcc; #if 0 @@ -159,20 +159,20 @@ void dubcos(double x, double dx, double v[]) { /* Routine receive Double-Length number (x+dx) and computes cos(x+dx) */ /* as Double-Length number and store it in array v */ /**********************************************************************/ -void docos(double x, double dx, double v[]) { +void __docos(double x, double dx, double v[]) { double y,yy,p,w[2]; if (x>0) {y=x; yy=dx;} else {y=-x; yy=-dx;} if (y<0.5*hp0.x) /* y< PI/4 */ - {dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];} + {__dubcos(y,yy,w); v[0]=w[0]; v[1]=w[1];} else if (y<1.5*hp0.x) { /* y< 3/4 * PI */ p=hp0.x-y; /* p = PI/2 - y */ yy=hp1.x-yy; y=p+yy; yy=(p-y)+yy; - if (y>0) {dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];} + if (y>0) {__dubsin(y,yy,w); v[0]=w[0]; v[1]=w[1];} /* cos(x) = sin ( 90 - x ) */ - else {dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1]; + else {__dubsin(-y,-yy,w); v[0]=-w[0]; v[1]=-w[1]; } } else { /* y>= 3/4 * PI */ @@ -180,7 +180,7 @@ void docos(double x, double dx, double v[]) { yy=2.0*hp1.x-yy; y=p+yy; yy=(p-y)+yy; - dubcos(y,yy,w); + __dubcos(y,yy,w); v[0]=-w[0]; v[1]=-w[1]; } diff --git a/sysdeps/ieee754/dbl-64/e_asin.c b/sysdeps/ieee754/dbl-64/e_asin.c index 2efa63166e..a2b309e4d5 100644 --- a/sysdeps/ieee754/dbl-64/e_asin.c +++ b/sysdeps/ieee754/dbl-64/e_asin.c @@ -312,6 +312,8 @@ double __ieee754_asin(double x){ } /* else if (k < 0x3ff00000) */ /*---------------------------- |x|>=1 -------------------------------*/ else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?hp0.x:-hp0.x; + else + if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x; else { u.i[HIGH_HALF]=0x7ff00000; v.i[HIGH_HALF]=0x7ff00000; @@ -621,6 +623,8 @@ double __ieee754_acos(double x) /*---------------------------- |x|>=1 -----------------------*/ else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?0:2.0*hp0.x; + else + if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x; else { u.i[HIGH_HALF]=0x7ff00000; v.i[HIGH_HALF]=0x7ff00000; diff --git a/sysdeps/ieee754/dbl-64/e_atan2.c b/sysdeps/ieee754/dbl-64/e_atan2.c index b77505db36..adf7a0d11b 100644 --- a/sysdeps/ieee754/dbl-64/e_atan2.c +++ b/sysdeps/ieee754/dbl-64/e_atan2.c @@ -374,9 +374,9 @@ static double normalized(double ax,double ay,double y, double z) { int p; mp_no mpx,mpy,mpz,mperr,mpz2,mpt1; p=6; - dbl_mp(ax,&mpx,p); dbl_mp(ay,&mpy,p); dvd(&mpy,&mpx,&mpz,p); - dbl_mp(ue.d,&mpt1,p); mul(&mpz,&mpt1,&mperr,p); - sub(&mpz,&mperr,&mpz2,p); __mp_dbl(&mpz2,&z,p); + __dbl_mp(ax,&mpx,p); __dbl_mp(ay,&mpy,p); __dvd(&mpy,&mpx,&mpz,p); + __dbl_mp(ue.d,&mpt1,p); __mul(&mpz,&mpt1,&mperr,p); + __sub(&mpz,&mperr,&mpz2,p); __mp_dbl(&mpz2,&z,p); return signArctan2(y,z); } /* Fix the sign and return after stage 1 or stage 2 */ @@ -392,10 +392,10 @@ static double atan2Mp(double x,double y,const int pr[]) mp_no mpx,mpy,mpz,mpz1,mpz2,mperr,mpt1; for (i=0; i<MM; i++) { p = pr[i]; - dbl_mp(x,&mpx,p); dbl_mp(y,&mpy,p); + __dbl_mp(x,&mpx,p); __dbl_mp(y,&mpy,p); __mpatan2(&mpy,&mpx,&mpz,p); - dbl_mp(ud[i].d,&mpt1,p); mul(&mpz,&mpt1,&mperr,p); - add(&mpz,&mperr,&mpz1,p); sub(&mpz,&mperr,&mpz2,p); + __dbl_mp(ud[i].d,&mpt1,p); __mul(&mpz,&mpt1,&mperr,p); + __add(&mpz,&mperr,&mpz1,p); __sub(&mpz,&mperr,&mpz2,p); __mp_dbl(&mpz1,&z1,p); __mp_dbl(&mpz2,&z2,p); if (z1==z2) return z1; } diff --git a/sysdeps/ieee754/dbl-64/e_lgamma_r.c b/sysdeps/ieee754/dbl-64/e_lgamma_r.c index e5b402a023..863eaa4362 100644 --- a/sysdeps/ieee754/dbl-64/e_lgamma_r.c +++ b/sysdeps/ieee754/dbl-64/e_lgamma_r.c @@ -175,7 +175,7 @@ static double zero= 0.00000000000000000000e+00; GET_HIGH_WORD(ix,x); ix &= 0x7fffffff; - if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0); + if(ix<0x3fd00000) return sin(pi*x); y = -x; /* x is assume negative */ /* diff --git a/sysdeps/ieee754/dbl-64/e_log.c b/sysdeps/ieee754/dbl-64/e_log.c index 814ac13720..1a099726d5 100644 --- a/sysdeps/ieee754/dbl-64/e_log.c +++ b/sysdeps/ieee754/dbl-64/e_log.c @@ -189,10 +189,10 @@ double __ieee754_log(double x) { for (i=0; i<M; i++) { p = pr[i]; - dbl_mp(x,&mpx,p); dbl_mp(y,&mpy,p); + __dbl_mp(x,&mpx,p); __dbl_mp(y,&mpy,p); __mplog(&mpx,&mpy,p); - dbl_mp(e[i].d,&mperr,p); - add(&mpy,&mperr,&mpy1,p); sub(&mpy,&mperr,&mpy2,p); + __dbl_mp(e[i].d,&mperr,p); + __add(&mpy,&mperr,&mpy1,p); __sub(&mpy,&mperr,&mpy2,p); __mp_dbl(&mpy1,&y1,p); __mp_dbl(&mpy2,&y2,p); if (y1==y2) return y1; } diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c index cd6f27f33e..dc92922d18 100644 --- a/sysdeps/ieee754/dbl-64/e_pow.c +++ b/sysdeps/ieee754/dbl-64/e_pow.c @@ -45,7 +45,7 @@ double __exp1(double x, double xx, double error); static double log1(double x, double *delta, double *error); static double log2(double x, double *delta, double *error); -double slowpow(double x, double y,double z); +double __slowpow(double x, double y,double z); static double power1(double x, double y); static int checkint(double x); @@ -69,8 +69,8 @@ double __ieee754_pow(double x, double y) { if (((qx==0x7ff00000)&&(u.i[LOW_HALF]!=0))||(qx>0x7ff00000)) return NaNQ.x; if (y == 1.0) return x; if (y == 2.0) return x*x; - if (y == -1.0) return (x!=0)?1.0/x:NaNQ.x; - if (y == 0) return ((x>0)&&(qx<0x7ff00000))?1.0:NaNQ.x; + if (y == -1.0) return 1.0/x; + if (y == 0) return 1.0; } /* else */ if(((u.i[HIGH_HALF]>0 && u.i[HIGH_HALF]<0x7ff00000)|| /* x>0 and not x->0 */ @@ -94,16 +94,37 @@ double __ieee754_pow(double x, double y) { } if (x == 0) { - if (ABS(y) > 1.0e20) return (y>0)?0:NaNQ.x; + if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0) + || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000) + return y; + if (ABS(y) > 1.0e20) return (y>0)?0:INF.x; k = checkint(y); - if (k == 0 || y < 0) return NaNQ.x; - else return (k==1)?0:x; /* return 0 */ + if (k == -1) + return y < 0 ? 1.0/x : x; + else + return y < 0 ? 1.0/ABS(x) : 0.0; /* return 0 */ } /* if x<0 */ if (u.i[HIGH_HALF] < 0) { k = checkint(y); - if (k==0) return NaNQ.x; /* y not integer and x<0 */ - return (k==1)?upow(-x,y):-upow(-x,y); /* if y even or odd */ + if (k==0) { + if ((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] == 0) { + if (x == -1.0) return 1.0; + else if (x > -1.0) return v.i[HIGH_HALF] < 0 ? INF.x : 0.0; + else return v.i[HIGH_HALF] < 0 ? 0.0 : INF.x; + } + else if (u.i[HIGH_HALF] == 0xfff00000 && u.i[LOW_HALF] == 0) + return y < 0 ? 0.0 : INF.x; + return NaNQ.x; /* y not integer and x<0 */ + } + else if (u.i[HIGH_HALF] == 0xfff00000 && u.i[LOW_HALF] == 0) + { + if (k < 0) + return y < 0 ? nZERO.x : nINF.x; + else + return y < 0 ? 0.0 : INF.x; + } + return (k==1)?__ieee754_pow(-x,y):-__ieee754_pow(-x,y); /* if y even or odd */ } /* x>0 */ qx = u.i[HIGH_HALF]&0x7fffffff; /* no sign */ @@ -111,7 +132,8 @@ double __ieee754_pow(double x, double y) { if (qx > 0x7ff00000 || (qx == 0x7ff00000 && u.i[LOW_HALF] != 0)) return NaNQ.x; /* if 0<x<2^-0x7fe */ - if (qy > 0x7ff00000 || (qy == 0x7ff00000 && v.i[LOW_HALF] != 0)) return NaNQ.x; + if (qy > 0x7ff00000 || (qy == 0x7ff00000 && v.i[LOW_HALF] != 0)) + return x == 1.0 ? 1.0 : NaNQ.x; /* if y<2^-0x7fe */ if (qx == 0x7ff00000) /* x= 2^-0x3ff */ @@ -124,7 +146,7 @@ double __ieee754_pow(double x, double y) { if (y<0) return (x<1.0)?INF.x:0; } - if (x == 1.0) return NaNQ.x; + if (x == 1.0) return 1.0; if (y>0) return (x>1.0)?INF.x:0; if (y<0) return (x<1.0)?INF.x:0; return 0; /* unreachable, to make the compiler happy */ @@ -148,7 +170,7 @@ static double power1(double x, double y) { a2 = (a-a1)+aa; error = error*ABS(y); t = __exp1(a1,a2,1.9e16*error); - return (t >= 0)?t:slowpow(x,y,z); + return (t >= 0)?t:__slowpow(x,y,z); } /****************************************************************************/ diff --git a/sysdeps/ieee754/dbl-64/e_remainder.c b/sysdeps/ieee754/dbl-64/e_remainder.c index f025fdf1b5..c59e5895d8 100644 --- a/sysdeps/ieee754/dbl-64/e_remainder.c +++ b/sysdeps/ieee754/dbl-64/e_remainder.c @@ -103,11 +103,13 @@ double __ieee754_remainder(double x, double y) else { if (kx<0x7ff00000&&ky<0x7ff00000&&(ky>0||t.i[LOW_HALF]!=0)) { y=ABS(y)*t128.x; - z=uremainder(x,y)*t128.x; - z=uremainder(z,y)*tm128.x; + z=__ieee754_remainder(x,y)*t128.x; + z=__ieee754_remainder(z,y)*tm128.x; return z; } else { /* if x is too big */ + if (kx == 0x7ff00000 && u.i[LOW_HALF] == 0 && y == 1.0) + return x / x; if (kx>=0x7ff00000||(ky==0&&t.i[LOW_HALF]==0)||ky>0x7ff00000|| (ky==0x7ff00000&&t.i[LOW_HALF]!=0)) return (u.i[HIGH_HALF]&0x80000000)?nNAN.x:NAN.x; diff --git a/sysdeps/ieee754/dbl-64/halfulp.c b/sysdeps/ieee754/dbl-64/halfulp.c index 7901ec4e36..929ca91c50 100644 --- a/sysdeps/ieee754/dbl-64/halfulp.c +++ b/sysdeps/ieee754/dbl-64/halfulp.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,12 +15,12 @@ * * 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:halfulp.c */ -/* */ +/* MODULE_NAME:halfulp.c */ +/* */ /* FUNCTIONS:halfulp */ /* FILES NEEDED: mydefs.h dla.h endian.h */ /* uroot.c */ @@ -35,11 +35,11 @@ /*3. if x can be represented by x=2**n for some integer n. */ /************************************************************************/ -#include "endian.h" +#include "endian.h" #include "mydefs.h" #include "dla.h" -double usqrt(double x); +double __ieee754_sqrt(double x); int4 tab54[32] = { 262143, 11585, 1782, 511, 210, 107, 63, 42, @@ -48,17 +48,17 @@ int4 tab54[32] = { 3, 3, 3, 3, 3, 3, 3, 3 }; -double halfulp(double x, double y) +double __halfulp(double x, double y) { mynumber v; double z,u,uu,j1,j2,j3,j4,j5; int4 k,l,m,n; if (y <= 0) { /*if power is negative or zero */ v.x = y; - if (v.i[LOW_HALF] != 0) return -10.0; + if (v.i[LOW_HALF] != 0) return -10.0; v.x = x; - if (v.i[LOW_HALF] != 0) return -10.0; - if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10; /* if x =2 ^ n */ + if (v.i[LOW_HALF] != 0) return -10.0; + if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10; /* if x =2 ^ n */ k = ((v.i[HIGH_HALF]&0x7fffffff)>>20)-1023; /* find this n */ z = (double) k; return (z*y == -1075.0)?0: -10.0; @@ -66,19 +66,19 @@ double halfulp(double x, double y) /* if y > 0 */ v.x = y; if (v.i[LOW_HALF] != 0) return -10.0; - + v.x=x; - /* case where x = 2**n for some integer n */ + /* case where x = 2**n for some integer n */ if (((v.i[HIGH_HALF]&0x000fffff)|v.i[LOW_HALF]) == 0) { k=(v.i[HIGH_HALF]>>20)-1023; return (((double) k)*y == -1075.0)?0:-10.0; - } - + } + v.x = y; k = v.i[HIGH_HALF]; m = k<<12; l = 0; - while (m) + while (m) {m = m<<1; l++; } n = (k&0x000fffff)|0x00100000; n = n>>(20-l); /* n is the odd integer of y */ @@ -88,19 +88,19 @@ double halfulp(double x, double y) if (n > 34) return -10.0; k = -k; if (k>5) return -10.0; - + /* now treat x */ while (k>0) { - z = usqrt(x); + z = __ieee754_sqrt(x); EMULV(z,z,u,uu,j1,j2,j3,j4,j5); if (((u-x)+uu) != 0) break; x = z; k--; } - if (k) return -10.0; - + if (k) return -10.0; + /* it is impossible that n == 2, so the mantissa of x must be short */ - + v.x = x; if (v.i[LOW_HALF]) return -10.0; k = v.i[HIGH_HALF]; @@ -109,16 +109,14 @@ double halfulp(double x, double y) while (m) {m = m<<1; l++; } m = (k&0x000fffff)|0x00100000; m = m>>(20-l); /* m is the odd integer of x */ - + /* now check whether the length of m**n is at most 54 bits */ - + if (m > tab54[n-3]) return -10.0; - + /* yes, it is - now compute x**n by simple multiplications */ - + u = x; for (k=1;k<n;k++) u = u*x; return u; } - - diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c index 1586d91bc2..e9840b037f 100644 --- a/sysdeps/ieee754/dbl-64/mpa.c +++ b/sysdeps/ieee754/dbl-64/mpa.c @@ -62,7 +62,7 @@ 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) { @@ -81,20 +81,20 @@ int acr(const mp_no *x, const mp_no *y, int p) { /* cr90 compares the values of two multiple precision numbers */ -int cr(const mp_no *x, const mp_no *y, int p) { +int __cr(const mp_no *x, const mp_no *y, int p) { int i; if (X[0] > Y[0]) i= 1; else if (X[0] < Y[0]) i=-1; - else if (X[0] < ZERO ) i= acr(y,x,p); - else i= acr(x,y,p); + else if (X[0] < ZERO ) i= __acr(y,x,p); + else i= __acr(x,y,p); return i; } /* Copy a multiple precision number. Set *y=*x. x=y is permissible. */ -void cpy(const mp_no *x, mp_no *y, int p) { +void __cpy(const mp_no *x, mp_no *y, int p) { int i; EY = EX; @@ -110,7 +110,7 @@ void cpy(const mp_no *x, mp_no *y, int p) { /* 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) { +void __cpymn(const mp_no *x, int m, mp_no *y, int n) { int i,k; @@ -246,7 +246,7 @@ void __mp_dbl(const mp_no *x, double *y, int p) { /* 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 __dbl_mp(double x, mp_no *y, int p) { int i,n; double u; @@ -286,7 +286,7 @@ static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { i=p; j=p+ EY - EX; k=p+1; if (j<1) - {cpy(x,z,p); return; } + {__cpy(x,z,p); return; } else Z[k] = ZERO; for (; j>0; i--,j--) { @@ -330,7 +330,7 @@ static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { Z[k] = Z[k+1] = ZERO; } else { j= EX - EY; - if (j > p) {cpy(x,z,p); return; } + if (j > p) {__cpy(x,z,p); return; } else { i=p; j=p+1-j; k=p; if (Y[j] > ZERO) { @@ -375,19 +375,19 @@ 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 __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] == 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]; } + 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]; } + 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; } @@ -399,19 +399,19 @@ 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 __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] == 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]; } + 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]; } + 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; } @@ -424,7 +424,7 @@ 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 __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { int i, i1, i2, j, k, k2; double u; @@ -464,7 +464,7 @@ 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. */ -void inv(const mp_no *x, mp_no *y, int p) { +void __inv(const mp_no *x, mp_no *y, int p) { int i; #if 0 int l; @@ -478,14 +478,14 @@ void inv(const mp_no *x, mp_no *y, int p) { 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; + __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); + __cpy(y,&w,p); + __mul(x,&w,y,p); + __sub(&mptwo,y,&z,p); + __mul(&w,&z,y,p); } return; } @@ -496,11 +496,11 @@ 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 __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);} + else {__inv(y,&w,p); __mul(x,&w,z,p);} return; } diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h index bb63e69002..d136283993 100644 --- a/sysdeps/ieee754/dbl-64/mpa.h +++ b/sysdeps/ieee754/dbl-64/mpa.h @@ -65,14 +65,14 @@ typedef union { int i[2]; double d; } number; #define MIN(x,y) ((x) < (y) ? (x) : (y)) #define ABS(x) ((x) < 0 ? -(x) : (x)) -int acr(const mp_no *, const mp_no *, int); -int cr(const mp_no *, const mp_no *, int); -void cpy(const mp_no *, mp_no *, int); -void cpymn(const mp_no *, int, mp_no *, int); +int __acr(const mp_no *, const mp_no *, int); +int __cr(const mp_no *, const mp_no *, int); +void __cpy(const mp_no *, mp_no *, int); +void __cpymn(const mp_no *, int, mp_no *, int); void __mp_dbl(const mp_no *, double *, int); -void dbl_mp(double, mp_no *, int); -void add(const mp_no *, const mp_no *, mp_no *, int); -void sub(const mp_no *, const mp_no *, mp_no *, int); -void mul(const mp_no *, const mp_no *, mp_no *, int); -void inv(const mp_no *, mp_no *, int); -void dvd(const mp_no *, const mp_no *, mp_no *, int); +void __dbl_mp(double, mp_no *, int); +void __add(const mp_no *, const mp_no *, mp_no *, int); +void __sub(const mp_no *, const mp_no *, mp_no *, int); +void __mul(const mp_no *, const mp_no *, mp_no *, int); +void __inv(const mp_no *, mp_no *, int); +void __dvd(const mp_no *, const mp_no *, mp_no *, int); diff --git a/sysdeps/ieee754/dbl-64/mpatan.c b/sysdeps/ieee754/dbl-64/mpatan.c index 416014bab7..0c8556db60 100644 --- a/sysdeps/ieee754/dbl-64/mpatan.c +++ b/sysdeps/ieee754/dbl-64/mpatan.c @@ -33,9 +33,9 @@ #include "endian.h" #include "mpa.h" -void mpsqrt(mp_no *, mp_no *, int); +void __mpsqrt(mp_no *, mp_no *, int); -void mpatan(mp_no *x, mp_no *y, int p) { +void __mpatan(mp_no *x, mp_no *y, int p) { #include "mpatan.h" int i,m,n; @@ -66,36 +66,36 @@ void mpatan(mp_no *x, mp_no *y, int p) { mptwo.d[1] = TWO; /* Reduce x m times */ - mul(x,x,&mpsm,p); - if (m==0) cpy(x,&mps,p); + __mul(x,x,&mpsm,p); + if (m==0) __cpy(x,&mps,p); else { for (i=0; i<m; i++) { - add(&mpone,&mpsm,&mpt1,p); - mpsqrt(&mpt1,&mpt2,p); - add(&mpt2,&mpt2,&mpt1,p); - add(&mptwo,&mpsm,&mpt2,p); - add(&mpt1,&mpt2,&mpt3,p); - dvd(&mpsm,&mpt3,&mpt1,p); - cpy(&mpt1,&mpsm,p); + __add(&mpone,&mpsm,&mpt1,p); + __mpsqrt(&mpt1,&mpt2,p); + __add(&mpt2,&mpt2,&mpt1,p); + __add(&mptwo,&mpsm,&mpt2,p); + __add(&mpt1,&mpt2,&mpt3,p); + __dvd(&mpsm,&mpt3,&mpt1,p); + __cpy(&mpt1,&mpsm,p); } - mpsqrt(&mpsm,&mps,p); mps.d[0] = X[0]; + __mpsqrt(&mpsm,&mps,p); mps.d[0] = X[0]; } /* Evaluate a truncated power series for Atan(s) */ n=np[p]; mptwoim1.d[1] = twonm1[p].d; - dvd(&mpsm,&mptwoim1,&mpt,p); + __dvd(&mpsm,&mptwoim1,&mpt,p); for (i=n-1; i>1; i--) { mptwoim1.d[1] -= TWO; - dvd(&mpsm,&mptwoim1,&mpt1,p); - mul(&mpsm,&mpt,&mpt2,p); - sub(&mpt1,&mpt2,&mpt,p); + __dvd(&mpsm,&mptwoim1,&mpt1,p); + __mul(&mpsm,&mpt,&mpt2,p); + __sub(&mpt1,&mpt2,&mpt,p); } - mul(&mps,&mpt,&mpt1,p); - sub(&mps,&mpt1,&mpt,p); + __mul(&mps,&mpt,&mpt1,p); + __sub(&mps,&mpt1,&mpt,p); /* Compute Atan(x) */ mptwoim1.d[1] = twom[m].d; - mul(&mptwoim1,&mpt,y,p); + __mul(&mptwoim1,&mpt,y,p); return; } diff --git a/sysdeps/ieee754/dbl-64/mpatan2.c b/sysdeps/ieee754/dbl-64/mpatan2.c index e41350e72a..2d1625b823 100644 --- a/sysdeps/ieee754/dbl-64/mpatan2.c +++ b/sysdeps/ieee754/dbl-64/mpatan2.c @@ -37,12 +37,12 @@ #include "mpa.h" -void mpsqrt(mp_no *, mp_no *, int); -void mpatan(mp_no *, mp_no *, int); +void __mpsqrt(mp_no *, mp_no *, int); +void __mpatan(mp_no *, mp_no *, int); /* Multi-Precision Atan2(y,x) function subroutine, for p >= 4. */ /* y=0 is not permitted if x<=0. No error messages are given. */ -void mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) { +void __mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) { static const double ZERO = 0.0, ONE = 1.0; @@ -54,15 +54,15 @@ void mpatan2(mp_no *y, mp_no *x, mp_no *z, int p) { if (X[0] <= ZERO) { mpone.e = 1; mpone.d[0] = mpone.d[1] = ONE; - dvd(x,y,&mpt1,p); mul(&mpt1,&mpt1,&mpt2,p); + __dvd(x,y,&mpt1,p); __mul(&mpt1,&mpt1,&mpt2,p); if (mpt1.d[0] != ZERO) mpt1.d[0] = ONE; - add(&mpt2,&mpone,&mpt3,p); mpsqrt(&mpt3,&mpt2,p); - add(&mpt1,&mpt2,&mpt3,p); mpt3.d[0]=Y[0]; - mpatan(&mpt3,&mpt1,p); add(&mpt1,&mpt1,z,p); + __add(&mpt2,&mpone,&mpt3,p); __mpsqrt(&mpt3,&mpt2,p); + __add(&mpt1,&mpt2,&mpt3,p); mpt3.d[0]=Y[0]; + __mpatan(&mpt3,&mpt1,p); __add(&mpt1,&mpt1,z,p); } else - { dvd(y,x,&mpt1,p); - mpatan(&mpt1,z,p); + { __dvd(y,x,&mpt1,p); + __mpatan(&mpt1,z,p); } return; diff --git a/sysdeps/ieee754/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c index 8b7c1dc4bb..2f0b826605 100644 --- a/sysdeps/ieee754/dbl-64/mpexp.c +++ b/sysdeps/ieee754/dbl-64/mpexp.c @@ -35,7 +35,7 @@ /* Multi-Precision exponential function subroutine (for p >= 4, */ /* 2**(-55) <= abs(x) <= 1024). */ -void mpexp(mp_no *x, mp_no *y, int p) { +void __mpexp(mp_no *x, mp_no *y, int p) { int i,j,k,m,m1,m2,n; double a,b; @@ -75,30 +75,30 @@ void mpexp(mp_no *x, mp_no *y, int p) { } /* Compute s=x*2**(-m). Put result in mps */ - dbl_mp(a,&mpt1,p); - mul(x,&mpt1,&mps,p); + __dbl_mp(a,&mpt1,p); + __mul(x,&mpt1,&mps,p); /* Evaluate the polynomial. Put result in mpt2 */ mpone.e=1; mpone.d[0]=ONE; mpone.d[1]=ONE; mpk.e = 1; mpk.d[0] = ONE; mpk.d[1]=nn[n].d; - dvd(&mps,&mpk,&mpt1,p); - add(&mpone,&mpt1,&mpak,p); + __dvd(&mps,&mpk,&mpt1,p); + __add(&mpone,&mpt1,&mpak,p); for (k=n-1; k>1; k--) { - mul(&mps,&mpak,&mpt1,p); + __mul(&mps,&mpak,&mpt1,p); mpk.d[1]=nn[k].d; - dvd(&mpt1,&mpk,&mpt2,p); - add(&mpone,&mpt2,&mpak,p); + __dvd(&mpt1,&mpk,&mpt2,p); + __add(&mpone,&mpt2,&mpak,p); } - mul(&mps,&mpak,&mpt1,p); - add(&mpone,&mpt1,&mpt2,p); + __mul(&mps,&mpak,&mpt1,p); + __add(&mpone,&mpt1,&mpt2,p); /* Raise polynomial value to the power of 2**m. Put result in y */ for (k=0,j=0; k<m; ) { - mul(&mpt2,&mpt2,&mpt1,p); k++; + __mul(&mpt2,&mpt2,&mpt1,p); k++; if (k==m) { j=1; break; } - mul(&mpt1,&mpt1,&mpt2,p); k++; + __mul(&mpt1,&mpt1,&mpt2,p); k++; } - if (j) cpy(&mpt1,y,p); - else cpy(&mpt2,y,p); + if (j) __cpy(&mpt1,y,p); + else __cpy(&mpt2,y,p); return; } diff --git a/sysdeps/ieee754/dbl-64/mplog.c b/sysdeps/ieee754/dbl-64/mplog.c index 5957c43cc2..cca0c7e4ed 100644 --- a/sysdeps/ieee754/dbl-64/mplog.c +++ b/sysdeps/ieee754/dbl-64/mplog.c @@ -37,9 +37,9 @@ #include "endian.h" #include "mpa.h" -void mpexp(mp_no *, mp_no *, int); +void __mpexp(mp_no *, mp_no *, int); -void mplog(mp_no *x, mp_no *y, int p) { +void __mplog(mp_no *x, mp_no *y, int p) { #include "mplog.h" int i,m; #if 0 @@ -58,14 +58,14 @@ void mplog(mp_no *x, mp_no *y, int p) { /* Perform m newton iterations to solve for y: exp(y)-x=0. */ /* The iterations formula is: y(n+1)=y(n)+(x*exp(-y(n))-1). */ - cpy(y,&mpt1,p); + __cpy(y,&mpt1,p); for (i=0; i<m; i++) { mpt1.d[0]=-mpt1.d[0]; - mpexp(&mpt1,&mpt2,p); - mul(x,&mpt2,&mpt1,p); - sub(&mpt1,&mpone,&mpt2,p); - add(y,&mpt2,&mpt1,p); - cpy(&mpt1,y,p); + __mpexp(&mpt1,&mpt2,p); + __mul(x,&mpt2,&mpt1,p); + __sub(&mpt1,&mpone,&mpt2,p); + __add(y,&mpt2,&mpt1,p); + __cpy(&mpt1,y,p); } return; } diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c index 3a4632bdc9..824c03f746 100644 --- a/sysdeps/ieee754/dbl-64/mpsqrt.c +++ b/sysdeps/ieee754/dbl-64/mpsqrt.c @@ -42,7 +42,7 @@ double fastiroot(double); -void mpsqrt(mp_no *x, mp_no *y, int p) { +void __mpsqrt(mp_no *x, mp_no *y, int p) { #include "mpsqrt.h" int i,m,ex,ey; @@ -60,19 +60,19 @@ void mpsqrt(mp_no *x, mp_no *y, int p) { mphalf.e =0; mphalf.d[0] =ONE; mphalf.d[1] =HALFRAD; mp3halfs.e=1; mp3halfs.d[0]=ONE; mp3halfs.d[1]=ONE; mp3halfs.d[2]=HALFRAD; - ex=EX; ey=EX/2; cpy(x,&mpxn,p); mpxn.e -= (ey+ey); - __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); dbl_mp(dy,&mpu,p); - mul(&mpxn,&mphalf,&mpz,p); + ex=EX; ey=EX/2; __cpy(x,&mpxn,p); mpxn.e -= (ey+ey); + __mp_dbl(&mpxn,&dx,p); dy=fastiroot(dx); __dbl_mp(dy,&mpu,p); + __mul(&mpxn,&mphalf,&mpz,p); m=mp[p]; for (i=0; i<m; i++) { - mul(&mpu,&mpu,&mpt1,p); - mul(&mpt1,&mpz,&mpt2,p); - sub(&mp3halfs,&mpt2,&mpt1,p); - mul(&mpu,&mpt1,&mpt2,p); - cpy(&mpt2,&mpu,p); + __mul(&mpu,&mpu,&mpt1,p); + __mul(&mpt1,&mpz,&mpt2,p); + __sub(&mp3halfs,&mpt2,&mpt1,p); + __mul(&mpu,&mpt1,&mpt2,p); + __cpy(&mpt2,&mpu,p); } - mul(&mpxn,&mpu,y,p); EY += ey; + __mul(&mpxn,&mpu,y,p); EY += ey; return; } diff --git a/sysdeps/ieee754/dbl-64/mptan.c b/sysdeps/ieee754/dbl-64/mptan.c index cbfbd0260b..b51a9c0901 100644 --- a/sysdeps/ieee754/dbl-64/mptan.c +++ b/sysdeps/ieee754/dbl-64/mptan.c @@ -37,23 +37,23 @@ #include "endian.h" #include "mpa.h" -int mpranred(double, mp_no *, int); +int __mpranred(double, mp_no *, int); void __c32(mp_no *, mp_no *, mp_no *, int); -void mptan(double x, mp_no *mpy, int p) { +void __mptan(double x, mp_no *mpy, int p) { static const double MONE = -1.0; int n; mp_no mpw, mpc, mps; - n = mpranred(x, &mpw, p) & 0x00000001; /* negative or positive result */ + n = __mpranred(x, &mpw, p) & 0x00000001; /* negative or positive result */ __c32(&mpw, &mpc, &mps, p); /* computing sin(x) and cos(x) */ if (n) /* second or fourth quarter of unit circle */ - { dvd(&mpc,&mps,mpy,p); + { __dvd(&mpc,&mps,mpy,p); mpy->d[0] *= MONE; } /* tan is negative in this area */ - else dvd(&mps,&mpc,mpy,p); + else __dvd(&mps,&mpc,mpy,p); return; } diff --git a/sysdeps/ieee754/dbl-64/s_atan.c b/sysdeps/ieee754/dbl-64/s_atan.c index 0b71c0d5c8..2872c75db7 100644 --- a/sysdeps/ieee754/dbl-64/s_atan.c +++ b/sysdeps/ieee754/dbl-64/s_atan.c @@ -214,9 +214,9 @@ static double atanMp(double x,const int pr[]){ for (i=0; i<M; i++) { p = pr[i]; - dbl_mp(x,&mpx,p); __mpatan(&mpx,&mpy,p); - dbl_mp(u9[i].d,&mpt1,p); mul(&mpy,&mpt1,&mperr,p); - add(&mpy,&mperr,&mpy1,p); sub(&mpy,&mperr,&mpy2,p); + __dbl_mp(x,&mpx,p); __mpatan(&mpx,&mpy,p); + __dbl_mp(u9[i].d,&mpt1,p); __mul(&mpy,&mpt1,&mperr,p); + __add(&mpy,&mperr,&mpy1,p); __sub(&mpy,&mperr,&mpy2,p); __mp_dbl(&mpy1,&y1,p); __mp_dbl(&mpy2,&y2,p); if (y1==y2) return y1; } diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c index da9ea27d77..ff6cf01bb8 100644 --- a/sysdeps/ieee754/dbl-64/s_sin.c +++ b/sysdeps/ieee754/dbl-64/s_sin.c @@ -60,8 +60,8 @@ static const double cs4 = -4.16666666666664434524222570944589E-02, cs6 = 1.38888874007937613028114285595617E-03; -void dubsin(double x, double dx, double w[]); -void docos(double x, double dx, double w[]); +void __dubsin(double x, double dx, double w[]); +void __docos(double x, double dx, double w[]); double __mpsin(double x, double dx); double __mpcos(double x, double dx); double __mpsin1(double x); @@ -75,7 +75,7 @@ static double sloww2(double x, double dx, double orig, int n); static double bsloww(double x, double dx, double orig, int n); static double bsloww1(double x, double dx, double orig, int n); static double bsloww2(double x, double dx, double orig, int n); -int branred(double x, double *a, double *aa); +int __branred(double x, double *a, double *aa); static double cslow2(double x); static double csloww(double x, double dx, double orig); static double csloww1(double x, double dx, double orig); @@ -84,7 +84,7 @@ static double csloww2(double x, double dx, double orig, int n); /* An ultimate sin routine. Given an IEEE double machine number x */ /* it computes the correctly rounded (to nearest) value of sin(x) */ /*******************************************************************/ -double sin(double x){ +double __sin(double x){ double xx,res,t,cor,y,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2; #if 0 double w[2]; @@ -307,7 +307,7 @@ double sin(double x){ /* -----------------281474976710656 <|x| <2^1024----------------------------*/ else if (k < 0x7ff00000) { - n = branred(x,&a,&da); + n = __branred(x,&a,&da); switch (n) { case 0: if (a*a < 0.01588) return bsloww(a,da,x,n); @@ -327,7 +327,7 @@ double sin(double x){ } /* else if (k < 0x7ff00000 ) */ /*--------------------- |x| > 2^1024 ----------------------------------*/ - else return NAN.x; + else return x / x; return 0; /* unreachable */ } @@ -337,7 +337,7 @@ double sin(double x){ /* it computes the correctly rounded (to nearest) value of cos(x) */ /*******************************************************************/ -double cos(double x) +double __cos(double x) { double y,xx,res,t,cor,s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2; mynumber u,v; @@ -548,7 +548,7 @@ double cos(double x) else if (k < 0x7ff00000) {/* 281474976710656 <|x| <2^1024 */ - n = branred(x,&a,&da); + n = __branred(x,&a,&da); switch (n) { case 1: if (a*a < 0.01588) return bsloww(-a,-da,x,n); @@ -570,7 +570,7 @@ double cos(double x) - else return NAN.x; /* |x| > 2^1024 */ + else return x / x; /* |x| > 2^1024 */ return 0; } @@ -594,7 +594,7 @@ static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ cor = (r-res)+t; if (res == res + 1.0007*cor) return res; else { - dubsin(ABS(x),0,w); + __dubsin(ABS(x),0,w); if (w[0] == w[0]+1.000000001*w[1]) return (x>0)?w[0]:-w[0]; else return (x>0)?__mpsin(x,0):-__mpsin(-x,0); } @@ -631,7 +631,7 @@ static double slow1(double x) { cor=(y-res)+cor; if (res == res+1.0005*cor) return (x>0)?res:-res; else { - dubsin(ABS(x),0,w); + __dubsin(ABS(x),0,w); if (w[0] == w[0]+1.000000005*w[1]) return (x>0)?w[0]:-w[0]; else return (x>0)?__mpsin(x,0):-__mpsin(-x,0); } @@ -679,7 +679,7 @@ static double slow2(double x) { y=ABS(x)-hp0.x; y1=y-hp1.x; y2=(y-y1)-hp1.x; - docos(y1,y2,w); + __docos(y1,y2,w); if (w[0] == w[0]+1.000000005*w[1]) return (x>0)?w[0]:-w[0]; else return (x>0)?__mpsin(x,0):-__mpsin(-x,0); } @@ -709,7 +709,7 @@ static double sloww(double x,double dx, double orig) { cor = (cor>0)? 1.0005*cor+ABS(orig)*3.1e-30 : 1.0005*cor-ABS(orig)*3.1e-30; if (res == res + cor) return res; else { - (x>0)? dubsin(x,dx,w) : dubsin(-x,-dx,w); + (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w); cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-30 : 1.000000001*w[1] - ABS(orig)*1.1e-30; if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0]; else { @@ -725,7 +725,7 @@ static double sloww(double x,double dx, double orig) { a = t - y; da = ((t-a)-y)+da; if (n&2) {a=-a; da=-da;} - (a>0)? dubsin(a,da,w) : dubsin(-a,-da,w); + (a>0)? __dubsin(a,da,w) : __dubsin(-a,-da,w); cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-40 : 1.000000001*w[1] - ABS(orig)*1.1e-40; if (w[0] == w[0]+cor) return (a>0)?w[0]:-w[0]; else return __mpsin1(orig); @@ -768,7 +768,7 @@ static double sloww1(double x, double dx, double orig) { cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig); if (res == res + cor) return (x>0)?res:-res; else { - dubsin(ABS(x),dx,w); + __dubsin(ABS(x),dx,w); cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig); if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0]; else return __mpsin1(orig); @@ -811,7 +811,7 @@ static double sloww2(double x, double dx, double orig, int n) { cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig); if (res == res + cor) return (n&2)?-res:res; else { - docos(ABS(x),dx,w); + __docos(ABS(x),dx,w); cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig); if (w[0] == w[0]+cor) return (n&2)?-w[0]:w[0]; else return __mpsin1(orig); @@ -844,7 +844,7 @@ static double bsloww(double x,double dx, double orig,int n) { cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24; if (res == res + cor) return res; else { - (x>0)? dubsin(x,dx,w) : dubsin(-x,-dx,w); + (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w); cor = (w[1]>0)? 1.000000001*w[1] + 1.1e-24 : 1.000000001*w[1] - 1.1e-24; if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0]; else return (n&1)?__mpcos1(orig):__mpsin1(orig); @@ -887,7 +887,7 @@ mynumber u; cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24; if (res == res + cor) return (x>0)?res:-res; else { - dubsin(ABS(x),dx,w); + __dubsin(ABS(x),dx,w); cor = (w[1]>0)? 1.000000005*w[1]+1.1e-24: 1.000000005*w[1]-1.1e-24; if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0]; else return (n&1)?__mpcos1(orig):__mpsin1(orig); @@ -931,7 +931,7 @@ mynumber u; cor = (cor>0)? 1.0005*cor+1.1e-24 : 1.0005*cor-1.1e-24; if (res == res + cor) return (n&2)?-res:res; else { - docos(ABS(x),dx,w); + __docos(ABS(x),dx,w); cor = (w[1]>0)? 1.000000005*w[1]+1.1e-24 : 1.000000005*w[1]-1.1e-24; if (w[0] == w[0]+cor) return (n&2)?-w[0]:w[0]; else return (n&1)?__mpsin1(orig):__mpcos1(orig); @@ -972,7 +972,7 @@ static double cslow2(double x) { return res; else { y=ABS(x); - docos(y,0,w); + __docos(y,0,w); if (w[0] == w[0]+1.000000005*w[1]) return w[0]; else return __mpcos(x,0); } @@ -1004,7 +1004,7 @@ static double csloww(double x,double dx, double orig) { cor = (cor>0)? 1.0005*cor+ABS(orig)*3.1e-30 : 1.0005*cor-ABS(orig)*3.1e-30; if (res == res + cor) return res; else { - (x>0)? dubsin(x,dx,w) : dubsin(-x,-dx,w); + (x>0)? __dubsin(x,dx,w) : __dubsin(-x,-dx,w); cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-30 : 1.000000001*w[1] - ABS(orig)*1.1e-30; if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0]; else { @@ -1020,7 +1020,7 @@ static double csloww(double x,double dx, double orig) { a = t - y; da = ((t-a)-y)+da; if (n==1) {a=-a; da=-da;} - (a>0)? dubsin(a,da,w) : dubsin(-a,-da,w); + (a>0)? __dubsin(a,da,w) : __dubsin(-a,-da,w); cor = (w[1]>0)? 1.000000001*w[1] + ABS(orig)*1.1e-40 : 1.000000001*w[1] - ABS(orig)*1.1e-40; if (w[0] == w[0]+cor) return (a>0)?w[0]:-w[0]; else return __mpcos1(orig); @@ -1064,7 +1064,7 @@ static double csloww1(double x, double dx, double orig) { cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig); if (res == res + cor) return (x>0)?res:-res; else { - dubsin(ABS(x),dx,w); + __dubsin(ABS(x),dx,w); cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig); if (w[0] == w[0]+cor) return (x>0)?w[0]:-w[0]; else return __mpcos1(orig); @@ -1109,14 +1109,19 @@ static double csloww2(double x, double dx, double orig, int n) { cor = (cor>0)? 1.0005*cor+3.1e-30*ABS(orig) : 1.0005*cor-3.1e-30*ABS(orig); if (res == res + cor) return (n)?-res:res; else { - docos(ABS(x),dx,w); + __docos(ABS(x),dx,w); cor = (w[1]>0)? 1.000000005*w[1]+1.1e-30*ABS(orig) : 1.000000005*w[1]-1.1e-30*ABS(orig); if (w[0] == w[0]+cor) return (n)?-w[0]:w[0]; else return __mpcos1(orig); } } +weak_alias (__cos, cos) +weak_alias (__sin, sin) + #ifdef NO_LONG_DOUBLE -weak_alias (sin, sinl) -weak_alias (cos, cosl) +strong_alias (__sin, __sinl) +weak_alias (__sin, sinl) +strong_alias (__cos, __cosl) +weak_alias (__cos, cosl) #endif diff --git a/sysdeps/ieee754/dbl-64/s_tan.c b/sysdeps/ieee754/dbl-64/s_tan.c index 26d53b1757..7b5dc4f3da 100644 --- a/sysdeps/ieee754/dbl-64/s_tan.c +++ b/sysdeps/ieee754/dbl-64/s_tan.c @@ -53,8 +53,8 @@ double tan(double x) { mp_no mpy; #endif - int branred(double, double *, double *); - int mpranred(double, mp_no *, int); + int __branred(double, double *, double *); + int __mpranred(double, mp_no *, int); /* x=+-INF, x=NaN */ num.d = x; ux = num.i[HIGH_HALF]; @@ -361,7 +361,7 @@ double tan(double x) { /* (---) The case 1e8 < abs(x) < 2**1024 */ /* Range reduction by algorithm iii */ - n = (branred(x,&a,&da)) & 0x00000001; + n = (__branred(x,&a,&da)) & 0x00000001; EADD(a,da,t1,t2) a=t1; da=t2; if (a<ZERO) {ya=-a; yya=-da; sy=MONE;} else {ya= a; yya= da; sy= ONE;} @@ -384,9 +384,9 @@ double tan(double x) { /* Second stage */ /* Reduction by algorithm iv */ - p=10; n = (mpranred(x,&mpa,p)) & 0x00000001; - __mp_dbl(&mpa,&a,p); dbl_mp(a,&mpt1,p); - sub(&mpa,&mpt1,&mpt2,p); __mp_dbl(&mpt2,&da,p); + p=10; n = (__mpranred(x,&mpa,p)) & 0x00000001; + __mp_dbl(&mpa,&a,p); __dbl_mp(a,&mpt1,p); + __sub(&mpa,&mpt1,&mpt2,p); __mp_dbl(&mpt2,&da,p); MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8) c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+ diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c index 6e8c9d5530..0fee643a54 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.c +++ b/sysdeps/ieee754/dbl-64/sincos32.c @@ -61,17 +61,17 @@ static void ss32(mp_no *x, mp_no *y, int p) { #endif for (i=1;i<=p;i++) mpk.d[i]=0; - mul(x,x,&x2,p); - cpy(&oofac27,&gor,p); - cpy(&gor,&sum,p); + __mul(x,x,&x2,p); + __cpy(&oofac27,&gor,p); + __cpy(&gor,&sum,p); for (a=27.0;a>1.0;a-=2.0) { mpk.d[1]=a*(a-1.0); - mul(&gor,&mpk,&mpt1,p); - cpy(&mpt1,&gor,p); - mul(&x2,&sum,&mpt1,p); - sub(&gor,&mpt1,&sum,p); + __mul(&gor,&mpk,&mpt1,p); + __cpy(&mpt1,&gor,p); + __mul(&x2,&sum,&mpt1,p); + __sub(&gor,&mpt1,&sum,p); } - mul(x,&sum,y,p); + __mul(x,&sum,y,p); } /**********************************************************************/ @@ -91,18 +91,18 @@ static void cc32(mp_no *x, mp_no *y, int p) { #endif for (i=1;i<=p;i++) mpk.d[i]=0; - mul(x,x,&x2,p); + __mul(x,x,&x2,p); mpk.d[1]=27.0; - mul(&oofac27,&mpk,&gor,p); - cpy(&gor,&sum,p); + __mul(&oofac27,&mpk,&gor,p); + __cpy(&gor,&sum,p); for (a=26.0;a>2.0;a-=2.0) { mpk.d[1]=a*(a-1.0); - mul(&gor,&mpk,&mpt1,p); - cpy(&mpt1,&gor,p); - mul(&x2,&sum,&mpt1,p); - sub(&gor,&mpt1,&sum,p); + __mul(&gor,&mpk,&mpt1,p); + __cpy(&mpt1,&gor,p); + __mul(&x2,&sum,&mpt1,p); + __sub(&gor,&mpt1,&sum,p); } - mul(&x2,&sum,y,p); + __mul(&x2,&sum,y,p); } /***************************************************************************/ @@ -112,20 +112,20 @@ void __c32(mp_no *x, mp_no *y, mp_no *z, int p) { static const mp_no mpt={1,{1.0,2.0}}, one={1,{1.0,1.0}}; mp_no u,t,t1,t2,c,s; int i; - cpy(x,&u,p); + __cpy(x,&u,p); u.e=u.e-1; cc32(&u,&c,p); ss32(&u,&s,p); for (i=0;i<24;i++) { - mul(&c,&s,&t,p); - sub(&s,&t,&t1,p); - add(&t1,&t1,&s,p); - sub(&mpt,&c,&t1,p); - mul(&t1,&c,&t2,p); - add(&t2,&t2,&c,p); + __mul(&c,&s,&t,p); + __sub(&s,&t,&t1,p); + __add(&t1,&t1,&s,p); + __sub(&mpt,&c,&t1,p); + __mul(&t1,&c,&t2,p); + __add(&t2,&t2,&c,p); } - sub(&one,&c,y,p); - cpy(&s,z,p); + __sub(&one,&c,y,p); + __cpy(&s,z,p); } /************************************************************************/ @@ -137,16 +137,16 @@ double __sin32(double x, double res, double res1) { int p; mp_no a,b,c; p=32; - dbl_mp(res,&a,p); - dbl_mp(0.5*(res1-res),&b,p); - add(&a,&b,&c,p); + __dbl_mp(res,&a,p); + __dbl_mp(0.5*(res1-res),&b,p); + __add(&a,&b,&c,p); if (x>0.8) - { sub(&hp,&c,&a,p); + { __sub(&hp,&c,&a,p); __c32(&a,&b,&c,p); } else __c32(&c,&a,&b,p); /* b=sin(0.5*(res+res1)) */ - dbl_mp(x,&c,p); /* c = x */ - sub(&b,&c,&a,p); + __dbl_mp(x,&c,p); /* c = x */ + __sub(&b,&c,&a,p); /* if a>0 return min(res,res1), otherwise return max(res,res1) */ if (a.d[0]>0) return (res<res1)?res:res1; else return (res>res1)?res:res1; @@ -161,21 +161,21 @@ double __cos32(double x, double res, double res1) { int p; mp_no a,b,c; p=32; - dbl_mp(res,&a,p); - dbl_mp(0.5*(res1-res),&b,p); - add(&a,&b,&c,p); + __dbl_mp(res,&a,p); + __dbl_mp(0.5*(res1-res),&b,p); + __add(&a,&b,&c,p); if (x>2.4) - { sub(&pi,&c,&a,p); + { __sub(&pi,&c,&a,p); __c32(&a,&b,&c,p); b.d[0]=-b.d[0]; } else if (x>0.8) - { sub(&hp,&c,&a,p); + { __sub(&hp,&c,&a,p); __c32(&a,&c,&b,p); } else __c32(&c,&b,&a,p); /* b=cos(0.5*(res+res1)) */ - dbl_mp(x,&c,p); /* c = x */ - sub(&b,&c,&a,p); + __dbl_mp(x,&c,p); /* c = x */ + __sub(&b,&c,&a,p); /* if a>0 return max(res,res1), otherwise return min(res,res1) */ if (a.d[0]>0) return (res>res1)?res:res1; else return (res<res1)?res:res1; @@ -190,10 +190,10 @@ double __mpsin(double x, double dx) { double y; mp_no a,b,c; p=32; - dbl_mp(x,&a,p); - dbl_mp(dx,&b,p); - add(&a,&b,&c,p); - if (x>0.8) { sub(&hp,&c,&a,p); __c32(&a,&b,&c,p); } + __dbl_mp(x,&a,p); + __dbl_mp(dx,&b,p); + __add(&a,&b,&c,p); + if (x>0.8) { __sub(&hp,&c,&a,p); __c32(&a,&b,&c,p); } else __c32(&c,&a,&b,p); /* b = sin(x+dx) */ __mp_dbl(&b,&y,p); return y; @@ -208,11 +208,11 @@ double __mpcos(double x, double dx) { double y; mp_no a,b,c; p=32; - dbl_mp(x,&a,p); - dbl_mp(dx,&b,p); - add(&a,&b,&c,p); + __dbl_mp(x,&a,p); + __dbl_mp(dx,&b,p); + __add(&a,&b,&c,p); if (x>0.8) - { sub(&hp,&c,&b,p); + { __sub(&hp,&c,&b,p); __c32(&b,&a,&c,p); } else __c32(&c,&a,&b,p); /* a = cos(x+dx) */ @@ -226,7 +226,7 @@ double __mpcos(double x, double dx) { /* n=0,+-1,+-2,.... */ /* Return int which indicates in which quarter of circle x is */ /******************************************************************/ -int mpranred(double x, mp_no *y, int p) +int __mpranred(double x, mp_no *y, int p) { number v; double t,xn; @@ -239,31 +239,31 @@ int mpranred(double x, mp_no *y, int p) xn = t - toint.d; v.d = t; n =v.i[LOW_HALF]&3; - dbl_mp(xn,&a,p); - mul(&a,&hp,&b,p); - dbl_mp(x,&c,p); - sub(&c,&b,y,p); + __dbl_mp(xn,&a,p); + __mul(&a,&hp,&b,p); + __dbl_mp(x,&c,p); + __sub(&c,&b,y,p); return n; } else { /* if x is very big more precision required */ - dbl_mp(x,&a,p); + __dbl_mp(x,&a,p); a.d[0]=1.0; k = a.e-5; if (k < 0) k=0; b.e = -k; b.d[0] = 1.0; for (i=0;i<p;i++) b.d[i+1] = toverp[i+k]; - mul(&a,&b,&c,p); + __mul(&a,&b,&c,p); t = c.d[c.e]; for (i=1;i<=p-c.e;i++) c.d[i]=c.d[i+c.e]; for (i=p+1-c.e;i<=p;i++) c.d[i]=0; c.e=0; if (c.d[1] >= 8388608.0) { t +=1.0; - sub(&c,&one,&b,p); - mul(&b,&hp,y,p); + __sub(&c,&one,&b,p); + __mul(&b,&hp,y,p); } - else mul(&c,&hp,y,p); + else __mul(&c,&hp,y,p); n = (int) t; if (x < 0) { y->d[0] = - y->d[0]; n = -n; } return (n&3); @@ -274,14 +274,14 @@ int mpranred(double x, mp_no *y, int p) /* Multi-Precision sin() function subroutine, for p=32. It is */ /* based on the routines mpranred() and c32(). */ /*******************************************************************/ -double mpsin1(double x) +double __mpsin1(double x) { int p; int n; mp_no u,s,c; double y; p=32; - n=mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */ + n=__mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */ __c32(&u,&c,&s,p); switch (n) { /* in which quarter of unit circle y is*/ case 0: @@ -313,7 +313,7 @@ double mpsin1(double x) /* based on the routines mpranred() and c32(). */ /*****************************************************************/ -double mpcos1(double x) +double __mpcos1(double x) { int p; int n; @@ -321,7 +321,7 @@ double mpcos1(double x) double y; p=32; - n=mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */ + n=__mpranred(x,&u,p); /* n is 0, 1, 2 or 3 */ __c32(&u,&c,&s,p); switch (n) { /* in what quarter of unit circle y is*/ diff --git a/sysdeps/ieee754/dbl-64/slowexp.c b/sysdeps/ieee754/dbl-64/slowexp.c index 1c2779b0e9..3a3758b305 100644 --- a/sysdeps/ieee754/dbl-64/slowexp.c +++ b/sysdeps/ieee754/dbl-64/slowexp.c @@ -30,10 +30,10 @@ /**************************************************************************/ #include "mpa.h" -void mpexp(mp_no *x, mp_no *y, int p); +void __mpexp(mp_no *x, mp_no *y, int p); /*Converting from double precision to Multi-precision and calculating e^x */ -double slowexp(double x) { +double __slowexp(double x) { double w,z,res,eps=3.0e-26; #if 0 double y; @@ -45,20 +45,20 @@ double slowexp(double x) { mp_no mpx, mpy, mpz,mpw,mpeps,mpcor; p=6; - dbl_mp(x,&mpx,p); /* Convert a double precision number x */ + __dbl_mp(x,&mpx,p); /* Convert a double precision number x */ /* into a multiple precision number mpx with prec. p. */ - mpexp(&mpx, &mpy, p); /* Multi-Precision exponential function */ - dbl_mp(eps,&mpeps,p); - mul(&mpeps,&mpy,&mpcor,p); - add(&mpy,&mpcor,&mpw,p); - sub(&mpy,&mpcor,&mpz,p); + __mpexp(&mpx, &mpy, p); /* Multi-Precision exponential function */ + __dbl_mp(eps,&mpeps,p); + __mul(&mpeps,&mpy,&mpcor,p); + __add(&mpy,&mpcor,&mpw,p); + __sub(&mpy,&mpcor,&mpz,p); __mp_dbl(&mpw, &w, p); __mp_dbl(&mpz, &z, p); if (w == z) return w; else { /* if calculating is not exactly */ p = 32; - dbl_mp(x,&mpx,p); - mpexp(&mpx, &mpy, p); + __dbl_mp(x,&mpx,p); + __mpexp(&mpx, &mpy, p); __mp_dbl(&mpy, &res, p); return res; } diff --git a/sysdeps/ieee754/dbl-64/slowpow.c b/sysdeps/ieee754/dbl-64/slowpow.c index 11da7e43d9..3412197aef 100644 --- a/sysdeps/ieee754/dbl-64/slowpow.c +++ b/sysdeps/ieee754/dbl-64/slowpow.c @@ -34,40 +34,40 @@ #include "mpa.h" -void mpexp(mp_no *x, mp_no *y, int p); -void mplog(mp_no *x, mp_no *y, int p); +void __mpexp(mp_no *x, mp_no *y, int p); +void __mplog(mp_no *x, mp_no *y, int p); double ulog(double); -double halfulp(double x,double y); +double __halfulp(double x,double y); -double slowpow(double x, double y, double z) { +double __slowpow(double x, double y, double z) { double res,res1; mp_no mpx, mpy, mpz,mpw,mpp,mpr,mpr1; static const mp_no eps = {-3,{1.0,4.0}}; int p; - res = halfulp(x,y); /* halfulp() returns -10 or x^y */ + res = __halfulp(x,y); /* halfulp() returns -10 or x^y */ if (res >= 0) return res; /* if result was really computed by halfulp */ /* else, if result was not really computed by halfulp */ p = 10; /* p=precision */ - dbl_mp(x,&mpx,p); - dbl_mp(y,&mpy,p); - dbl_mp(z,&mpz,p); - mplog(&mpx, &mpz, p); /* log(x) = z */ - mul(&mpy,&mpz,&mpw,p); /* y * z =w */ - mpexp(&mpw, &mpp, p); /* e^w =pp */ - add(&mpp,&eps,&mpr,p); /* pp+eps =r */ + __dbl_mp(x,&mpx,p); + __dbl_mp(y,&mpy,p); + __dbl_mp(z,&mpz,p); + __mplog(&mpx, &mpz, p); /* log(x) = z */ + __mul(&mpy,&mpz,&mpw,p); /* y * z =w */ + __mpexp(&mpw, &mpp, p); /* e^w =pp */ + __add(&mpp,&eps,&mpr,p); /* pp+eps =r */ __mp_dbl(&mpr, &res, p); - sub(&mpp,&eps,&mpr1,p); /* pp -eps =r1 */ + __sub(&mpp,&eps,&mpr1,p); /* pp -eps =r1 */ __mp_dbl(&mpr1, &res1, p); /* converting into double precision */ if (res == res1) return res; p = 32; /* if we get here result wasn't calculated exactly, continue */ - dbl_mp(x,&mpx,p); /* for more exact calculation */ - dbl_mp(y,&mpy,p); - dbl_mp(z,&mpz,p); - mplog(&mpx, &mpz, p); /* log(c)=z */ - mul(&mpy,&mpz,&mpw,p); /* y*z =w */ - mpexp(&mpw, &mpp, p); /* e^w=pp */ + __dbl_mp(x,&mpx,p); /* for more exact calculation */ + __dbl_mp(y,&mpy,p); + __dbl_mp(z,&mpz,p); + __mplog(&mpx, &mpz, p); /* log(c)=z */ + __mul(&mpy,&mpz,&mpw,p); /* y*z =w */ + __mpexp(&mpw, &mpp, p); /* e^w=pp */ __mp_dbl(&mpp, &res, p); /* converting into double precision */ return res; } |