diff options
Diffstat (limited to 'sysdeps/ieee754')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/dla.h | 119 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/doasin.c | 7 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/dosincos.c | 20 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_atan2.c | 7 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_log.c | 5 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_pow.c | 5 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/halfulp.c | 5 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_atan.c | 115 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_tan.c | 59 |
9 files changed, 192 insertions, 150 deletions
diff --git a/sysdeps/ieee754/dbl-64/dla.h b/sysdeps/ieee754/dbl-64/dla.h index bf73fa902e..9f095f9bf5 100644 --- a/sysdeps/ieee754/dbl-64/dla.h +++ b/sysdeps/ieee754/dbl-64/dla.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001 Free Software Foundation, Inc. + * Copyright (C) 2001, 2011 Free Software Foundation, Inc. * * 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 @@ -35,6 +35,18 @@ /* IEEE double. */ /***********************************************************************/ +/* We can use fma instructions if available. */ +#if defined __x86_64__ || (defined __i386__ && defined __SSE2_MATH__) +# ifdef __FMA4__ +# define DLA_FMA(x,y,z) \ + ({ double __zz; \ + asm ("vfmaddsd %3, %2, %1, %0" \ + : "=x" (__zz) : "x" (x), "xm" (y), "x" (-z)); \ + __zz; }) +# endif +#endif + + /* CN = 1+2**27 = '41a0000002000000' IEEE double format */ #define CN 134217729.0 @@ -44,7 +56,7 @@ /* z+zz = x+y exactly. */ #define EADD(x,y,z,zz) \ - z=(x)+(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))+(y)) : (((y)-(z))+(x)); + z=(x)+(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))+(y)) : (((y)-(z))+(x)); /* Exact subtraction of two single-length floating point numbers, Dekker. */ @@ -52,7 +64,7 @@ /* z+zz = x-y exactly. */ #define ESUB(x,y,z,zz) \ - z=(x)-(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))-(y)) : ((x)-((y)+(z))); + z=(x)-(y); zz=(ABS(x)>ABS(y)) ? (((x)-(z))-(y)) : ((x)-((y)+(z))); /* Exact multiplication of two single-length floating point numbers, */ @@ -60,10 +72,15 @@ /* satisfies z+zz = x*y exactly. p,hx,tx,hy,ty are temporary */ /* storage variables of type double. */ -#define EMULV(x,y,z,zz,p,hx,tx,hy,ty) \ - p=CN*(x); hx=((x)-p)+p; tx=(x)-hx; \ - p=CN*(y); hy=((y)-p)+p; ty=(y)-hy; \ - z=(x)*(y); zz=(((hx*hy-z)+hx*ty)+tx*hy)+tx*ty; +#ifdef DLA_FMA +# define EMULV(x,y,z,zz,p,hx,tx,hy,ty) \ + z=x*y; zz=DLA_FMA(x,y,z); +#else +# define EMULV(x,y,z,zz,p,hx,tx,hy,ty) \ + p=CN*(x); hx=((x)-p)+p; tx=(x)-hx; \ + p=CN*(y); hy=((y)-p)+p; ty=(y)-hy; \ + z=(x)*(y); zz=(((hx*hy-z)+hx*ty)+tx*hy)+tx*ty; +#endif /* Exact multiplication of two single-length floating point numbers, Dekker. */ @@ -71,10 +88,15 @@ /* that satisfies z+zz = x*y exactly. p,hx,tx,hy,ty,q are temporary */ /* storage variables of type double. */ -#define MUL12(x,y,z,zz,p,hx,tx,hy,ty,q) \ - p=CN*(x); hx=((x)-p)+p; tx=(x)-hx; \ - p=CN*(y); hy=((y)-p)+p; ty=(y)-hy; \ - p=hx*hy; q=hx*ty+tx*hy; z=p+q; zz=((p-z)+q)+tx*ty; +#ifdef DLA_FMA +# define MUL12(x,y,z,zz,p,hx,tx,hy,ty,q) \ + EMULV(x,y,z,zz,p,hx,tx,hy,ty) +#else +# define MUL12(x,y,z,zz,p,hx,tx,hy,ty,q) \ + p=CN*(x); hx=((x)-p)+p; tx=(x)-hx; \ + p=CN*(y); hy=((y)-p)+p; ty=(y)-hy; \ + p=hx*hy; q=hx*ty+tx*hy; z=p+q; zz=((p-z)+q)+tx*ty; +#endif /* Double-length addition, Dekker. The macro produces a double-length */ @@ -84,10 +106,10 @@ /* storage variables of type double. */ #define ADD2(x,xx,y,yy,z,zz,r,s) \ - r=(x)+(y); s=(ABS(x)>ABS(y)) ? \ - (((((x)-r)+(y))+(yy))+(xx)) : \ - (((((y)-r)+(x))+(xx))+(yy)); \ - z=r+s; zz=(r-z)+s; + r=(x)+(y); s=(ABS(x)>ABS(y)) ? \ + (((((x)-r)+(y))+(yy))+(xx)) : \ + (((((y)-r)+(x))+(xx))+(yy)); \ + z=r+s; zz=(r-z)+s; /* Double-length subtraction, Dekker. The macro produces a double-length */ @@ -97,10 +119,10 @@ /* storage variables of type double. */ #define SUB2(x,xx,y,yy,z,zz,r,s) \ - r=(x)-(y); s=(ABS(x)>ABS(y)) ? \ - (((((x)-r)-(y))-(yy))+(xx)) : \ - ((((x)-((y)+r))+(xx))-(yy)); \ - z=r+s; zz=(r-z)+s; + r=(x)-(y); s=(ABS(x)>ABS(y)) ? \ + (((((x)-r)-(y))-(yy))+(xx)) : \ + ((((x)-((y)+r))+(xx))-(yy)); \ + z=r+s; zz=(r-z)+s; /* Double-length multiplication, Dekker. The macro produces a double-length */ @@ -110,8 +132,8 @@ /* temporary storage variables of type double. */ #define MUL2(x,xx,y,yy,z,zz,p,hx,tx,hy,ty,q,c,cc) \ - MUL12(x,y,c,cc,p,hx,tx,hy,ty,q) \ - cc=((x)*(yy)+(xx)*(y))+cc; z=c+cc; zz=(c-z)+cc; + MUL12(x,y,c,cc,p,hx,tx,hy,ty,q) \ + cc=((x)*(yy)+(xx)*(y))+cc; z=c+cc; zz=(c-z)+cc; /* Double-length division, Dekker. The macro produces a double-length */ @@ -121,8 +143,8 @@ /* are temporary storage variables of type double. */ #define DIV2(x,xx,y,yy,z,zz,p,hx,tx,hy,ty,q,c,cc,u,uu) \ - c=(x)/(y); MUL12(c,y,u,uu,p,hx,tx,hy,ty,q) \ - cc=(((((x)-u)-uu)+(xx))-c*(yy))/(y); z=c+cc; zz=(c-z)+cc; + c=(x)/(y); MUL12(c,y,u,uu,p,hx,tx,hy,ty,q) \ + cc=(((((x)-u)-uu)+(xx))-c*(yy))/(y); z=c+cc; zz=(c-z)+cc; /* Double-length addition, slower but more accurate than ADD2. */ @@ -133,17 +155,17 @@ /* are temporary storage variables of type double. */ #define ADD2A(x,xx,y,yy,z,zz,r,rr,s,ss,u,uu,w) \ - r=(x)+(y); \ - if (ABS(x)>ABS(y)) { rr=((x)-r)+(y); s=(rr+(yy))+(xx); } \ - else { rr=((y)-r)+(x); s=(rr+(xx))+(yy); } \ - if (rr!=0.0) { \ - z=r+s; zz=(r-z)+s; } \ - else { \ - ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)+(yy)) : (((yy)-s)+(xx)); \ - u=r+s; \ - uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \ - w=uu+ss; z=u+w; \ - zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; } + r=(x)+(y); \ + if (ABS(x)>ABS(y)) { rr=((x)-r)+(y); s=(rr+(yy))+(xx); } \ + else { rr=((y)-r)+(x); s=(rr+(xx))+(yy); } \ + if (rr!=0.0) { \ + z=r+s; zz=(r-z)+s; } \ + else { \ + ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)+(yy)) : (((yy)-s)+(xx)); \ + u=r+s; \ + uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \ + w=uu+ss; z=u+w; \ + zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; } /* Double-length subtraction, slower but more accurate than SUB2. */ @@ -154,21 +176,14 @@ /* are temporary storage variables of type double. */ #define SUB2A(x,xx,y,yy,z,zz,r,rr,s,ss,u,uu,w) \ - r=(x)-(y); \ - if (ABS(x)>ABS(y)) { rr=((x)-r)-(y); s=(rr-(yy))+(xx); } \ - else { rr=(x)-((y)+r); s=(rr+(xx))-(yy); } \ - if (rr!=0.0) { \ - z=r+s; zz=(r-z)+s; } \ - else { \ - ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)-(yy)) : ((xx)-((yy)+s)); \ - u=r+s; \ - uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \ - w=uu+ss; z=u+w; \ - zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; } - - - - - - - + r=(x)-(y); \ + if (ABS(x)>ABS(y)) { rr=((x)-r)-(y); s=(rr-(yy))+(xx); } \ + else { rr=(x)-((y)+r); s=(rr+(xx))-(yy); } \ + if (rr!=0.0) { \ + z=r+s; zz=(r-z)+s; } \ + else { \ + ss=(ABS(xx)>ABS(yy)) ? (((xx)-s)-(yy)) : ((xx)-((yy)+s)); \ + u=r+s; \ + uu=(ABS(r)>ABS(s)) ? ((r-u)+s) : ((s-u)+r) ; \ + w=uu+ss; z=u+w; \ + zz=(ABS(u)>ABS(w)) ? ((u-z)+w) : ((w-z)+u) ; } diff --git a/sysdeps/ieee754/dbl-64/doasin.c b/sysdeps/ieee754/dbl-64/doasin.c index 79f344af2d..9ed7609541 100644 --- a/sysdeps/ieee754/dbl-64/doasin.c +++ b/sysdeps/ieee754/dbl-64/doasin.c @@ -1,7 +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 @@ -52,7 +52,10 @@ void __doasin(double x, double dx, double v[]) { d11 = 0.79470250400727425881446981833568758E-02; double xx,p,pp,u,uu,r,s; - double hx,tx,hy,ty,tp,tq,tc,tcc; + double tc,tcc; +#ifndef DLA_FMA + double hx,tx,hy,ty,tp,tq; +#endif /* Taylor series for arcsin for Double-Length numbers */ diff --git a/sysdeps/ieee754/dbl-64/dosincos.c b/sysdeps/ieee754/dbl-64/dosincos.c index 1d347a4bc7..654f3424e8 100644 --- a/sysdeps/ieee754/dbl-64/dosincos.c +++ b/sysdeps/ieee754/dbl-64/dosincos.c @@ -1,7 +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,8 +48,11 @@ /***********************************************************************/ 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, + double r,s,c,cc,d,dd,d2,dd2,e,ee, sn,ssn,cs,ccs,ds,dss,dc,dcc; +#ifndef DLA_FMA + double p,hx,tx,hy,ty,q; +#endif #if 0 double xx,y,yy,z,zz; #endif @@ -61,7 +64,7 @@ void __dubsin(double x, double dx, double v[]) { x=x-(u.x-big.x); d=x+dx; dd=(x-d)+dx; - /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */ + /* sin(x+dx)=sin(Xi+t)=sin(Xi)*cos(t) + cos(Xi)sin(t) where t ->0 */ MUL2(d,dd,d,dd,d2,dd2,p,hx,tx,hy,ty,q,c,cc); sn=sincos.x[k]; /* */ ssn=sincos.x[k+1]; /* sin(Xi) and cos(Xi) */ @@ -99,8 +102,11 @@ void __dubsin(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, + double r,s,c,cc,d,dd,d2,dd2,e,ee, sn,ssn,cs,ccs,ds,dss,dc,dcc; +#ifndef DLA_FMA + double p,hx,tx,hy,ty,q; +#endif #if 0 double xx,y,yy,z,zz; #endif @@ -166,15 +172,15 @@ void __docos(double x, double dx, double v[]) { 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];} - /* cos(x) = sin ( 90 - x ) */ - else {__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 { /* y>= 3/4 * PI */ diff --git a/sysdeps/ieee754/dbl-64/e_atan2.c b/sysdeps/ieee754/dbl-64/e_atan2.c index 784fc5a0c3..4d8c23af95 100644 --- a/sysdeps/ieee754/dbl-64/e_atan2.c +++ b/sysdeps/ieee754/dbl-64/e_atan2.c @@ -1,7 +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 @@ -62,8 +62,11 @@ double __ieee754_atan2(double y,double x) { int p; #endif static const int pr[MM]={6,8,10,20,32}; - double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t4,t5,t6,t7,t8, + double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t7,t8, z,zz,cor,s1,ss1,s2,ss2; +#ifndef DLA_FMA + double t4,t5,t6; +#endif #if 0 double z1,z2; #endif diff --git a/sysdeps/ieee754/dbl-64/e_log.c b/sysdeps/ieee754/dbl-64/e_log.c index 0fc6f920a6..c158c8be30 100644 --- a/sysdeps/ieee754/dbl-64/e_log.c +++ b/sysdeps/ieee754/dbl-64/e_log.c @@ -56,8 +56,11 @@ double __ieee754_log(double x) { #endif double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj, sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb, - t1,t2,t3,t4,t5,t6,t7,t8,t,ra,rb,ww, + t1,t2,t7,t8,t,ra,rb,ww, a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c; +#ifndef DLA_FMA + double t3,t4,t5,t6; +#endif number num; mp_no mpx,mpy,mpy1,mpy2,mperr; diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c index 83a5eff5c2..643e1cb9cf 100644 --- a/sysdeps/ieee754/dbl-64/e_pow.c +++ b/sysdeps/ieee754/dbl-64/e_pow.c @@ -283,7 +283,10 @@ static double my_log2(double x, double *delta, double *error) { double cor; #endif double ou1,ou2,lu1,lu2,ov,lv1,lv2,a,a1,a2; - double y,yy,z,zz,j1,j2,j3,j4,j5,j6,j7,j8; + double y,yy,z,zz,j1,j2,j7,j8; +#ifndef DLA_FMA + double j3,j4,j5,j6; +#endif mynumber u,v; #ifdef BIG_ENDI mynumber diff --git a/sysdeps/ieee754/dbl-64/halfulp.c b/sysdeps/ieee754/dbl-64/halfulp.c index 42b21fb61d..5d27334630 100644 --- a/sysdeps/ieee754/dbl-64/halfulp.c +++ b/sysdeps/ieee754/dbl-64/halfulp.c @@ -50,7 +50,10 @@ static const int4 tab54[32] = { double __halfulp(double x, double y) { mynumber v; - double z,u,uu,j1,j2,j3,j4,j5; + double z,u,uu; +#ifndef DLA_FMA + double j1,j2,j3,j4,j5; +#endif int4 k,l,m,n; if (y <= 0) { /*if power is negative or zero */ v.x = y; diff --git a/sysdeps/ieee754/dbl-64/s_atan.c b/sysdeps/ieee754/dbl-64/s_atan.c index 556f5b216d..b948f503a6 100644 --- a/sysdeps/ieee754/dbl-64/s_atan.c +++ b/sysdeps/ieee754/dbl-64/s_atan.c @@ -1,7 +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 @@ -52,8 +52,11 @@ double __signArctan(double,double); double atan(double x) { - double cor,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,u,u2,u3, - v,vv,w,ww,y,yy,z,zz; + double cor,s1,ss1,s2,ss2,t1,t2,t3,t7,t8,t9,t10,u,u2,u3, + v,vv,w,ww,y,yy,z,zz; +#ifndef DLA_FMA + double t4,t5,t6; +#endif #if 0 double y1,y2; #endif @@ -78,44 +81,44 @@ double atan(double x) { if (u<C) { if (u<B) { if (u<A) { /* u < A */ - return x; } + return x; } else { /* A <= u < B */ - v=x*x; yy=x*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d))))); - if ((y=x+(yy-U1*x)) == x+(yy+U1*x)) return y; - - EMULV(x,x,v,vv,t1,t2,t3,t4,t5) /* v+vv=x^2 */ - s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d)))); - ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - MUL2(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2) - if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1)) return y; - - return atanMp(x,pr); + v=x*x; yy=x*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d))))); + if ((y=x+(yy-U1*x)) == x+(yy+U1*x)) return y; + + EMULV(x,x,v,vv,t1,t2,t3,t4,t5) /* v+vv=x^2 */ + s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d)))); + ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + MUL2(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2) + if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1)) return y; + + return atanMp(x,pr); } } else { /* B <= u < C */ i=(TWO52+TWO8*u)-TWO52; i-=16; z=u-cij[i][0].d; yy=z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+ - z*(cij[i][5].d+z* cij[i][6].d)))); + z*(cij[i][5].d+z* cij[i][6].d)))); t1=cij[i][1].d; if (i<112) { - if (i<48) u2=U21; /* u < 1/4 */ - else u2=U22; } /* 1/4 <= u < 1/2 */ + if (i<48) u2=U21; /* u < 1/4 */ + else u2=U22; } /* 1/4 <= u < 1/2 */ else { - if (i<176) u2=U23; /* 1/2 <= u < 3/4 */ - else u2=U24; } /* 3/4 <= u <= 1 */ + if (i<176) u2=U23; /* 1/2 <= u < 3/4 */ + else u2=U24; } /* 3/4 <= u <= 1 */ if ((y=t1+(yy-u2*t1)) == t1+(yy+u2*t1)) return __signArctan(x,y); z=u-hij[i][0].d; s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+ - z*(hij[i][14].d+z* hij[i][15].d)))); + z*(hij[i][14].d+z* hij[i][15].d)))); ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2) MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2) @@ -138,7 +141,7 @@ double atan(double x) { i=(TWO52+TWO8*w)-TWO52; i-=16; z=(w-cij[i][0].d)+ww; yy=HPI1-z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+ - z*(cij[i][5].d+z* cij[i][6].d)))); + z*(cij[i][5].d+z* cij[i][6].d)))); t1=HPI-cij[i][1].d; if (i<112) u3=U31; /* w < 1/2 */ else u3=U32; /* w >= 1/2 */ @@ -148,7 +151,7 @@ double atan(double x) { t1=w-hij[i][0].d; EADD(t1,ww,z,zz) s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+ - z*(hij[i][14].d+z* hij[i][15].d)))); + z*(hij[i][14].d+z* hij[i][15].d)))); ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2) MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2) @@ -165,36 +168,36 @@ double atan(double x) { } else { if (u<E) { /* D <= u < E */ - w=ONE/u; v=w*w; - EMULV(w,u,t1,t2,t3,t4,t5,t6,t7) - yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d))))); - ww=w*((ONE-t1)-t2); - ESUB(HPI,w,t3,cor) - yy=((HPI1+cor)-ww)-yy; - if ((y=t3+(yy-U4)) == t3+(yy+U4)) return __signArctan(x,y); - - DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10) - MUL2(w,ww,w,ww,v,vv,t1,t2,t3,t4,t5,t6,t7,t8) - s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d)))); - ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2) - MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) - MUL2(w,ww,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) - ADD2(w,ww,s2,ss2,s1,ss1,t1,t2) - SUB2(HPI,HPI1,s1,ss1,s2,ss2,t1,t2) - if ((y=s2+(ss2-U8)) == s2+(ss2+U8)) return __signArctan(x,y); + w=ONE/u; v=w*w; + EMULV(w,u,t1,t2,t3,t4,t5,t6,t7) + yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d))))); + ww=w*((ONE-t1)-t2); + ESUB(HPI,w,t3,cor) + yy=((HPI1+cor)-ww)-yy; + if ((y=t3+(yy-U4)) == t3+(yy+U4)) return __signArctan(x,y); + + DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10) + MUL2(w,ww,w,ww,v,vv,t1,t2,t3,t4,t5,t6,t7,t8) + s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d)))); + ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2) + MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8) + MUL2(w,ww,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(w,ww,s2,ss2,s1,ss1,t1,t2) + SUB2(HPI,HPI1,s1,ss1,s2,ss2,t1,t2) + if ((y=s2+(ss2-U8)) == s2+(ss2+U8)) return __signArctan(x,y); return atanMp(x,pr); } else { - /* u >= E */ - if (x>0) return HPI; - else return MHPI; } + /* u >= E */ + if (x>0) return HPI; + else return MHPI; } } } diff --git a/sysdeps/ieee754/dbl-64/s_tan.c b/sysdeps/ieee754/dbl-64/s_tan.c index 4e26d90ae1..015b027dc4 100644 --- a/sysdeps/ieee754/dbl-64/s_tan.c +++ b/sysdeps/ieee754/dbl-64/s_tan.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001, 2009 Free Software Foundation + * Copyright (C) 2001, 2009, 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 @@ -50,7 +50,10 @@ double tan(double x) { int ux,i,n; double a,da,a2,b,db,c,dc,c1,cc1,c2,cc2,c3,cc3,fi,ffi,gi,pz,s,sy, - t,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,w,x2,xn,xx2,y,ya,yya,z0,z,zz,z2,zz2; + t,t1,t2,t3,t4,t7,t8,t9,t10,w,x2,xn,xx2,y,ya,yya,z0,z,zz,z2,zz2; +#ifndef DLA_FMA + double t5,t6; +#endif int p; number num,v; mp_no mpa,mpt1,mpt2; @@ -84,7 +87,7 @@ double tan(double x) { /* Second stage */ c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+ - x2*a27.d)))))); + x2*a27.d)))))); EMULV(x,x,x2,xx2,t1,t2,t3,t4,t5) ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2) MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8) @@ -159,13 +162,13 @@ double tan(double x) { a2 = a*a; t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d)))); if (n) { - /* First stage -cot */ - EADD(a,t2,b,db) - DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10) - if ((y=c+(dc-u6.d*c))==c+(dc+u6.d*c)) return (-y); } + /* First stage -cot */ + EADD(a,t2,b,db) + DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10) + if ((y=c+(dc-u6.d*c))==c+(dc+u6.d*c)) return (-y); } else { - /* First stage tan */ - if ((y=a+(t2-u5.d*a))==a+(t2+u5.d*a)) return y; } + /* First stage tan */ + if ((y=a+(t2-u5.d*a))==a+(t2+u5.d*a)) return y; } /* Second stage */ /* Range reduction by algorithm ii */ t = (x*hpinv.d + toint.d); @@ -184,7 +187,7 @@ double tan(double x) { EADD(a,da,t1,t2) a=t1; da=t2; 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+ - x2*a27.d)))))); + x2*a27.d)))))); ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2) MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8) ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2) @@ -201,12 +204,12 @@ double tan(double x) { ADD2(a ,da ,c2,cc2,c1,cc1,t1,t2) if (n) { - /* Second stage -cot */ - DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10) - if ((y=c2+(cc2-u8.d*c2)) == c2+(cc2+u8.d*c2)) return (-y); } + /* Second stage -cot */ + DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10) + if ((y=c2+(cc2-u8.d*c2)) == c2+(cc2+u8.d*c2)) return (-y); } else { - /* Second stage tan */ - if ((y=c1+(cc1-u7.d*c1)) == c1+(cc1+u7.d*c1)) return y; } + /* Second stage tan */ + if ((y=c1+(cc1-u7.d*c1)) == c1+(cc1+u7.d*c1)) return y; } return tanMp(x); } @@ -287,18 +290,18 @@ double tan(double x) { a2 = a*a; t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d)))); if (n) { - /* First stage -cot */ - EADD(a,t2,b,db) - DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10) - if ((y=c+(dc-u14.d*c))==c+(dc+u14.d*c)) return (-y); } + /* First stage -cot */ + EADD(a,t2,b,db) + DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10) + if ((y=c+(dc-u14.d*c))==c+(dc+u14.d*c)) return (-y); } else { - /* First stage tan */ - if ((y=a+(t2-u13.d*a))==a+(t2+u13.d*a)) return y; } + /* First stage tan */ + if ((y=a+(t2-u13.d*a))==a+(t2+u13.d*a)) return y; } /* Second stage */ 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+ - x2*a27.d)))))); + x2*a27.d)))))); ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2) MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8) ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2) @@ -315,12 +318,12 @@ double tan(double x) { ADD2(a ,da ,c2,cc2,c1,cc1,t1,t2) if (n) { - /* Second stage -cot */ - DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10) - if ((y=c2+(cc2-u16.d*c2)) == c2+(cc2+u16.d*c2)) return (-y); } + /* Second stage -cot */ + DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10) + if ((y=c2+(cc2-u16.d*c2)) == c2+(cc2+u16.d*c2)) return (-y); } else { - /* Second stage tan */ - if ((y=c1+(cc1-u15.d*c1)) == c1+(cc1+u15.d*c1)) return (y); } + /* Second stage tan */ + if ((y=c1+(cc1-u15.d*c1)) == c1+(cc1+u15.d*c1)) return (y); } return tanMp(x); } @@ -404,7 +407,7 @@ double tan(double x) { 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+ - x2*a27.d)))))); + x2*a27.d)))))); ADD2(a13.d,aa13.d,c1,zero.d,c2,cc2,t1,t2) MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8) ADD2(a11.d,aa11.d,c1,cc1,c2,cc2,t1,t2) |