diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/s_sin.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_sin.c | 1178 |
1 files changed, 1102 insertions, 76 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c index 376c69ed00..b40c47f389 100644 --- a/sysdeps/ieee754/dbl-64/s_sin.c +++ b/sysdeps/ieee754/dbl-64/s_sin.c @@ -1,87 +1,1113 @@ -/* @(#)s_sin.c 5.1 93/09/24 */ /* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * IBM Accurate Mathematical Library + * Copyright (c) International Business Machines Corp., 2001 * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_sin.c,v 1.7 1995/05/10 20:48:15 jtc Exp $"; -#endif - -/* sin(x) - * Return sine function of x. - * - * kernel function: - * __kernel_sin ... sine function on [-pi/4,pi/4] - * __kernel_cos ... cose function on [-pi/4,pi/4] - * __ieee754_rem_pio2 ... argument reduction routine - * - * Method. - * Let S,C and T denote the sin, cos and tan respectively on - * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 - * in [-pi/4 , +pi/4], and let n = k mod 4. - * We have - * - * n sin(x) cos(x) tan(x) - * ---------------------------------------------------------- - * 0 S C T - * 1 C -S -1/T - * 2 -S -C T - * 3 -C S -1/T - * ---------------------------------------------------------- + * 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 + * (at your option) any later version. * - * Special cases: - * Let trig be any of sin, cos, or tan. - * trig(+-INF) is NaN, with signals; - * trig(NaN) is that NaN; + * 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 + * GNU General Public License for more details. * - * Accuracy: - * TRIG(x) returns trig(x) nearly rounded + * 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. */ +/****************************************************************************/ +/* */ +/* MODULE_NAME:usncs.c */ +/* */ +/* FUNCTIONS: usin */ +/* ucos */ +/* slow */ +/* slow1 */ +/* slow2 */ +/* sloww */ +/* sloww1 */ +/* sloww2 */ +/* bsloww */ +/* bsloww1 */ +/* bsloww2 */ +/* cslow2 */ +/* csloww */ +/* csloww1 */ +/* csloww2 */ +/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */ +/* branred.c sincos32.c dosincos.c mpa.c */ +/* sincos.tbl */ +/* */ +/* An ultimate sin and routine. Given an IEEE double machine number x */ +/* it computes the correctly rounded (to nearest) value of sin(x) or cos(x) */ +/* Assumption: Machine arithmetic operations are performed in */ +/* round to nearest mode of IEEE 754 standard. */ +/* */ +/****************************************************************************/ -#include "math.h" -#include "math_private.h" -#ifdef __STDC__ - double __sin(double x) -#else - double __sin(x) - double x; -#endif -{ - double y[2],z=0.0; - int32_t n, ix; - - /* High word of x. */ - GET_HIGH_WORD(ix,x); - - /* |x| ~< pi/4 */ - ix &= 0x7fffffff; - if(ix <= 0x3fe921fb) return __kernel_sin(x,z,0); - - /* sin(Inf or NaN) is NaN */ - else if (ix>=0x7ff00000) return x-x; - - /* argument reduction needed */ - else { - n = __ieee754_rem_pio2(x,y); - switch(n&3) { - case 0: return __kernel_sin(y[0],y[1],1); - case 1: return __kernel_cos(y[0],y[1]); - case 2: return -__kernel_sin(y[0],y[1],1); - default: - return -__kernel_cos(y[0],y[1]); +#include "endian.h" +#include "mydefs.h" +#include "usncs.h" +#include "MathLib.h" +#include "sincos.tbl" + +static const double + sn3 = -1.66666666666664880952546298448555E-01, + sn5 = 8.33333214285722277379541354343671E-03, + cs2 = 4.99999999999999999999950396842453E-01, + cs4 = -4.16666666666664434524222570944589E-02, + cs6 = 1.38888874007937613028114285595617E-03; + +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); +double mpcos1(double x); +static double slow(double x); +static double slow1(double x); +static double slow2(double x); +static double sloww(double x, double dx, double orig); +static double sloww1(double x, double dx, double orig); +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); +static double cslow2(double x); +static double csloww(double x, double dx, double orig); +static double csloww1(double x, double dx, double orig); +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 xx,res,t,cor,y,w[2],s,c,sn,ssn,cs,ccs,xn,a,da,db,eps,xn1,xn2; + mynumber u,v; + int4 k,m,n,nn; + + u.x = x; + m = u.i[HIGH_HALF]; + k = 0x7fffffff&m; /* no sign */ + if (k < 0x3e500000) /* if x->0 =>sin(x)=x */ + return x; + /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/ + else if (k < 0x3fd00000){ + xx = x*x; + /*Taylor series */ + t = ((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*(xx*x); + res = x+t; + cor = (x-res)+t; + return (res == res + 1.07*cor)? res : slow(x); + } /* else if (k < 0x3fd00000) */ +/*---------------------------- 0.25<|x|< 0.855469---------------------- */ + else if (k < 0x3feb6000) { + u.x=(m>0)?big.x+x:big.x-x; + y=(m>0)?x-(u.x-big.x):x+(u.x-big.x); + xx=y*y; + s = y + y*xx*(sn3 +xx*sn5); + c = xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=(m>0)?sincos.x[k]:-sincos.x[k]; + ssn=(m>0)?sincos.x[k+1]:-sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + cor=(ssn+s*ccs-sn*c)+cs*s; + res=sn+cor; + cor=(sn-res)+cor; + return (res==res+1.025*cor)? res : slow1(x); + } /* else if (k < 0x3feb6000) */ + +/*----------------------- 0.855469 <|x|<2.426265 ----------------------*/ + else if (k < 0x400368fd ) { + + y = (m>0)? hp0.x-x:hp0.x+x; + if (y>=0) { + u.x = big.x+y; + y = (y-(u.x-big.x))+hp1.x; + } + else { + u.x = big.x-y; + y = (-hp1.x) - (y+(u.x-big.x)); + } + xx=y*y; + s = y + y*xx*(sn3 +xx*sn5); + c = xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + cor=(ccs-s*ssn-cs*c)-sn*s; + res=cs+cor; + cor=(cs-res)+cor; + return (res==res+1.020*cor)? ((m>0)?res:-res) : slow2(x); + } /* else if (k < 0x400368fd) */ + +/*-------------------------- 2.426265<|x|< 105414350 ----------------------*/ + else if (k < 0x419921FB ) { + t = (x*hpinv.x + toint.x); + xn = t - toint.x; + v.x = t; + y = (x - xn*mp1.x) - xn*mp2.x; + n =v.i[LOW_HALF]&3; + da = xn*mp3.x; + a=y-da; + da = (y-a)-da; + eps = ABS(x)*1.2e-30; + + switch (n) { /* quarter of unit circle */ + case 0: + case 2: + xx = a*a; + if (n) {a=-a;da=-da;} + if (xx < 0.01588) { + /*Taylor series */ + t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da; + res = a+t; + cor = (a-res)+t; + cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps; + return (res == res + cor)? res : sloww(a,da,x); } - } + else { + if (a>0) + {m=1;t=a;db=da;} + else + {m=0;t=-a;db=-da;} + u.x=big.x+t; + y=t-(u.x-big.x); + xx=y*y; + s = y + (db+y*xx*(sn3 +xx*sn5)); + c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + cor=(ssn+s*ccs-sn*c)+cs*s; + res=sn+cor; + cor=(sn-res)+cor; + cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps; + return (res==res+cor)? ((m)?res:-res) : sloww1(a,da,x); + } + break; + + case 1: + case 3: + if (a<0) + {a=-a;da=-da;} + u.x=big.x+a; + y=a-(u.x-big.x)+da; + xx=y*y; + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + s = y + y*xx*(sn3 +xx*sn5); + c = xx*(cs2 +xx*(cs4 + xx*cs6)); + cor=(ccs-s*ssn-cs*c)-sn*s; + res=cs+cor; + cor=(cs-res)+cor; + cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps; + return (res==res+cor)? ((n&2)?-res:res) : sloww2(a,da,x,n); + + break; + + } + + } /* else if (k < 0x419921FB ) */ + +/*---------------------105414350 <|x|< 281474976710656 --------------------*/ + else if (k < 0x42F00000 ) { + t = (x*hpinv.x + toint.x); + xn = t - toint.x; + v.x = t; + xn1 = (xn+8.0e22)-8.0e22; + xn2 = xn - xn1; + y = ((((x - xn1*mp1.x) - xn1*mp2.x)-xn2*mp1.x)-xn2*mp2.x); + n =v.i[LOW_HALF]&3; + da = xn1*pp3.x; + t=y-da; + da = (y-t)-da; + da = (da - xn2*pp3.x) -xn*pp4.x; + a = t+da; + da = (t-a)+da; + eps = 1.0e-24; + + switch (n) { + case 0: + case 2: + xx = a*a; + if (n) {a=-a;da=-da;} + if (xx < 0.01588) { + /* Taylor series */ + t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da; + res = a+t; + cor = (a-res)+t; + cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps; + return (res == res + cor)? res : bsloww(a,da,x,n); + } + else { + if (a>0) {m=1;t=a;db=da;} + else {m=0;t=-a;db=-da;} + u.x=big.x+t; + y=t-(u.x-big.x); + xx=y*y; + s = y + (db+y*xx*(sn3 +xx*sn5)); + c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + cor=(ssn+s*ccs-sn*c)+cs*s; + res=sn+cor; + cor=(sn-res)+cor; + cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps; + return (res==res+cor)? ((m)?res:-res) : bsloww1(a,da,x,n); + } + break; + + case 1: + case 3: + if (a<0) + {a=-a;da=-da;} + u.x=big.x+a; + y=a-(u.x-big.x)+da; + xx=y*y; + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + s = y + y*xx*(sn3 +xx*sn5); + c = xx*(cs2 +xx*(cs4 + xx*cs6)); + cor=(ccs-s*ssn-cs*c)-sn*s; + res=cs+cor; + cor=(cs-res)+cor; + cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps; + return (res==res+cor)? ((n&2)?-res:res) : bsloww2(a,da,x,n); + + break; + + } + + } /* else if (k < 0x42F00000 ) */ + +/* -----------------281474976710656 <|x| <2^1024----------------------------*/ + else if (k < 0x7ff00000) { + + n = branred(x,&a,&da); + switch (n) { + case 0: + if (a*a < 0.01588) return bsloww(a,da,x,n); + else return bsloww1(a,da,x,n); + break; + case 2: + if (a*a < 0.01588) return bsloww(-a,-da,x,n); + else return bsloww1(-a,-da,x,n); + break; + + case 1: + case 3: + return bsloww2(a,da,x,n); + break; + } + + } /* else if (k < 0x7ff00000 ) */ + +/*--------------------- |x| > 2^1024 ----------------------------------*/ + else return NAN.x; + return 0; /* unreachable */ +} + + +/*******************************************************************/ +/* An ultimate cos routine. Given an IEEE double machine number x */ +/* it computes the correctly rounded (to nearest) value of cos(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; + int4 k,m,n; + + u.x = x; + m = u.i[HIGH_HALF]; + k = 0x7fffffff&m; + + if (k < 0x3e400000 ) return 1.0; /* |x|<2^-27 => cos(x)=1 */ + + else if (k < 0x3feb6000 ) {/* 2^-27 < |x| < 0.855469 */ + y=ABS(x); + u.x = big.x+y; + y = y-(u.x-big.x); + xx=y*y; + s = y + y*xx*(sn3 +xx*sn5); + c = xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + cor=(ccs-s*ssn-cs*c)-sn*s; + res=cs+cor; + cor=(cs-res)+cor; + return (res==res+1.020*cor)? res : cslow2(x); + +} /* else if (k < 0x3feb6000) */ + + else if (k < 0x400368fd ) {/* 0.855469 <|x|<2.426265 */; + y=hp0.x-ABS(x); + a=y+hp1.x; + da=(y-a)+hp1.x; + xx=a*a; + if (xx < 0.01588) { + t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da; + res = a+t; + cor = (a-res)+t; + cor = (cor>0)? 1.02*cor+1.0e-31 : 1.02*cor -1.0e-31; + return (res == res + cor)? res : csloww(a,da,x); + } + else { + if (a>0) {m=1;t=a;db=da;} + else {m=0;t=-a;db=-da;} + u.x=big.x+t; + y=t-(u.x-big.x); + xx=y*y; + s = y + (db+y*xx*(sn3 +xx*sn5)); + c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + cor=(ssn+s*ccs-sn*c)+cs*s; + res=sn+cor; + cor=(sn-res)+cor; + cor = (cor>0)? 1.035*cor+1.0e-31 : 1.035*cor-1.0e-31; + return (res==res+cor)? ((m)?res:-res) : csloww1(a,da,x); } -weak_alias (__sin, sin) + +} /* else if (k < 0x400368fd) */ + + + else if (k < 0x419921FB ) {/* 2.426265<|x|< 105414350 */ + t = (x*hpinv.x + toint.x); + xn = t - toint.x; + v.x = t; + y = (x - xn*mp1.x) - xn*mp2.x; + n =v.i[LOW_HALF]&3; + da = xn*mp3.x; + a=y-da; + da = (y-a)-da; + eps = ABS(x)*1.2e-30; + + switch (n) { + case 1: + case 3: + xx = a*a; + if (n == 1) {a=-a;da=-da;} + if (xx < 0.01588) { + t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da; + res = a+t; + cor = (a-res)+t; + cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps; + return (res == res + cor)? res : csloww(a,da,x); + } + else { + if (a>0) {m=1;t=a;db=da;} + else {m=0;t=-a;db=-da;} + u.x=big.x+t; + y=t-(u.x-big.x); + xx=y*y; + s = y + (db+y*xx*(sn3 +xx*sn5)); + c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + cor=(ssn+s*ccs-sn*c)+cs*s; + res=sn+cor; + cor=(sn-res)+cor; + cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps; + return (res==res+cor)? ((m)?res:-res) : csloww1(a,da,x); + } + break; + + case 0: + case 2: + if (a<0) {a=-a;da=-da;} + u.x=big.x+a; + y=a-(u.x-big.x)+da; + xx=y*y; + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + s = y + y*xx*(sn3 +xx*sn5); + c = xx*(cs2 +xx*(cs4 + xx*cs6)); + cor=(ccs-s*ssn-cs*c)-sn*s; + res=cs+cor; + cor=(cs-res)+cor; + cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps; + return (res==res+cor)? ((n)?-res:res) : csloww2(a,da,x,n); + + break; + + } + + } /* else if (k < 0x419921FB ) */ + + + else if (k < 0x42F00000 ) { + t = (x*hpinv.x + toint.x); + xn = t - toint.x; + v.x = t; + xn1 = (xn+8.0e22)-8.0e22; + xn2 = xn - xn1; + y = ((((x - xn1*mp1.x) - xn1*mp2.x)-xn2*mp1.x)-xn2*mp2.x); + n =v.i[LOW_HALF]&3; + da = xn1*pp3.x; + t=y-da; + da = (y-t)-da; + da = (da - xn2*pp3.x) -xn*pp4.x; + a = t+da; + da = (t-a)+da; + eps = 1.0e-24; + + switch (n) { + case 1: + case 3: + xx = a*a; + if (n==1) {a=-a;da=-da;} + if (xx < 0.01588) { + t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*a - 0.5*da)*xx+da; + res = a+t; + cor = (a-res)+t; + cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps; + return (res == res + cor)? res : bsloww(a,da,x,n); + } + else { + if (a>0) {m=1;t=a;db=da;} + else {m=0;t=-a;db=-da;} + u.x=big.x+t; + y=t-(u.x-big.x); + xx=y*y; + s = y + (db+y*xx*(sn3 +xx*sn5)); + c = y*db+xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + cor=(ssn+s*ccs-sn*c)+cs*s; + res=sn+cor; + cor=(sn-res)+cor; + cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps; + return (res==res+cor)? ((m)?res:-res) : bsloww1(a,da,x,n); + } + break; + + case 0: + case 2: + if (a<0) {a=-a;da=-da;} + u.x=big.x+a; + y=a-(u.x-big.x)+da; + xx=y*y; + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + s = y + y*xx*(sn3 +xx*sn5); + c = xx*(cs2 +xx*(cs4 + xx*cs6)); + cor=(ccs-s*ssn-cs*c)-sn*s; + res=cs+cor; + cor=(cs-res)+cor; + cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps; + return (res==res+cor)? ((n)?-res:res) : bsloww2(a,da,x,n); + break; + + } + + } /* else if (k < 0x42F00000 ) */ + + else if (k < 0x7ff00000) {/* 281474976710656 <|x| <2^1024 */ + + n = branred(x,&a,&da); + switch (n) { + case 1: + if (a*a < 0.01588) return bsloww(-a,-da,x,n); + else return bsloww1(-a,-da,x,n); + break; + case 3: + if (a*a < 0.01588) return bsloww(a,da,x,n); + else return bsloww1(a,da,x,n); + break; + + case 0: + case 2: + return bsloww2(a,da,x,n); + break; + } + + } /* else if (k < 0x7ff00000 ) */ + + + + + else return NAN.x; /* |x| > 2^1024 */ + return 0; + +} + +/************************************************************************/ +/* Routine compute sin(x) for 2^-26 < |x|< 0.25 by Taylor with more */ +/* precision and if still doesn't accurate enough by mpsin or dubsin */ +/************************************************************************/ + +static double slow(double x) { +static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ + double y,x1,x2,xx,r,t,res,cor,w[2]; + x1=(x+th2_36)-th2_36; + y = aa.x*x1*x1*x1; + r=x+y; + x2=x-x1; + xx=x*x; + t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + bb.x)*xx + 3.0*aa.x*x1*x2)*x +aa.x*x2*x2*x2; + t=((x-r)+y)+t; + res=r+t; + cor = (r-res)+t; + if (res == res + 1.0007*cor) return res; + else { + 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); + } +} +/*******************************************************************************/ +/* Routine compute sin(x) for 0.25<|x|< 0.855469 by sincos.tbl and Taylor */ +/* and if result still doesn't accurate enough by mpsin or dubsin */ +/*******************************************************************************/ + +static double slow1(double x) { + mynumber u; + double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; + static const double t22 = 6291456.0; + int4 k; + y=ABS(x); + u.x=big.x+y; + y=y-(u.x-big.x); + xx=y*y; + s = y*xx*(sn3 +xx*sn5); + c = xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; /* Data */ + ssn=sincos.x[k+1]; /* from */ + cs=sincos.x[k+2]; /* tables */ + ccs=sincos.x[k+3]; /* sincos.tbl */ + y1 = (y+t22)-t22; + y2 = y - y1; + c1 = (cs+t22)-t22; + c2=(cs-c1)+ccs; + cor=(ssn+s*ccs+cs*s+c2*y+c1*y2)-sn*c; + y=sn+c1*y1; + cor = cor+((sn-y)+c1*y1); + res=y+cor; + cor=(y-res)+cor; + if (res == res+1.0005*cor) return (x>0)?res:-res; + else { + 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); + } +} +/**************************************************************************/ +/* Routine compute sin(x) for 0.855469 <|x|<2.426265 by sincos.tbl */ +/* and if result still doesn't accurate enough by mpsin or dubsin */ +/**************************************************************************/ +static double slow2(double x) { + mynumber u; + double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res,del; + static const double t22 = 6291456.0; + int4 k; + y=ABS(x); + y = hp0.x-y; + if (y>=0) { + u.x = big.x+y; + y = y-(u.x-big.x); + del = hp1.x; + } + else { + u.x = big.x-y; + y = -(y+(u.x-big.x)); + del = -hp1.x; + } + xx=y*y; + s = y*xx*(sn3 +xx*sn5); + c = y*del+xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + y1 = (y+t22)-t22; + y2 = (y - y1)+del; + e1 = (sn+t22)-t22; + e2=(sn-e1)+ssn; + cor=(ccs-cs*c-e1*y2-e2*y)-sn*s; + y=cs-e1*y1; + cor = cor+((cs-y)-e1*y1); + res=y+cor; + cor=(y-res)+cor; + if (res == res+1.0005*cor) return (x>0)?res:-res; + else { + y=ABS(x)-hp0.x; + y1=y-hp1.x; + y2=(y-y1)-hp1.x; + 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); + } +} +/***************************************************************************/ +/* Routine compute sin(x+dx) (Double-Length number) where x is small enough*/ +/* to use Taylor series around zero and (x+dx) */ +/* in first or third quarter of unit circle.Routine receive also */ +/* (right argument) the original value of x for computing error of */ +/* result.And if result not accurate enough routine calls mpsin1 or dubsin */ +/***************************************************************************/ + +static double sloww(double x,double dx, double orig) { + static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ + double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn; + union {int4 i[2]; double x;} v; + int4 n; + x1=(x+th2_36)-th2_36; + y = aa.x*x1*x1*x1; + r=x+y; + x2=(x-x1)+dx; + xx=x*x; + t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + bb.x)*xx + 3.0*aa.x*x1*x2)*x +aa.x*x2*x2*x2+dx; + t=((x-r)+y)+t; + res=r+t; + cor = (r-res)+t; + 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); + 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 { + t = (orig*hpinv.x + toint.x); + xn = t - toint.x; + v.x = t; + y = (orig - xn*mp1.x) - xn*mp2.x; + n =v.i[LOW_HALF]&3; + da = xn*pp3.x; + t=y-da; + da = (y-t)-da; + y = xn*pp4.x; + 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); + 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); + } + } +} +/***************************************************************************/ +/* Routine compute sin(x+dx) (Double-Length number) where x in first or */ +/* third quarter of unit circle.Routine receive also (right argument) the */ +/* original value of x for computing error of result.And if result not */ +/* accurate enough routine calls mpsin1 or dubsin */ +/***************************************************************************/ + +static double sloww1(double x, double dx, double orig) { + mynumber u; + double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; + static const double t22 = 6291456.0; + int4 k; + y=ABS(x); + u.x=big.x+y; + y=y-(u.x-big.x); + dx=(x>0)?dx:-dx; + xx=y*y; + s = y*xx*(sn3 +xx*sn5); + c = xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + y1 = (y+t22)-t22; + y2 = (y - y1)+dx; + c1 = (cs+t22)-t22; + c2=(cs-c1)+ccs; + cor=(ssn+s*ccs+cs*s+c2*y+c1*y2-sn*y*dx)-sn*c; + y=sn+c1*y1; + cor = cor+((sn-y)+c1*y1); + res=y+cor; + cor=(y-res)+cor; + 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); + 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); + } +} +/***************************************************************************/ +/* Routine compute sin(x+dx) (Double-Length number) where x in second or */ +/* fourth quarter of unit circle.Routine receive also the original value */ +/* and quarter(n= 1or 3)of x for computing error of result.And if result not*/ +/* accurate enough routine calls mpsin1 or dubsin */ +/***************************************************************************/ + +static double sloww2(double x, double dx, double orig, int n) { + mynumber u; + double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; + static const double t22 = 6291456.0; + int4 k; + y=ABS(x); + u.x=big.x+y; + y=y-(u.x-big.x); + dx=(x>0)?dx:-dx; + xx=y*y; + s = y*xx*(sn3 +xx*sn5); + c = y*dx+xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + + y1 = (y+t22)-t22; + y2 = (y - y1)+dx; + e1 = (sn+t22)-t22; + e2=(sn-e1)+ssn; + cor=(ccs-cs*c-e1*y2-e2*y)-sn*s; + y=cs-e1*y1; + cor = cor+((cs-y)-e1*y1); + res=y+cor; + cor=(y-res)+cor; + 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); + 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); + } +} +/***************************************************************************/ +/* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */ +/* is small enough to use Taylor series around zero and (x+dx) */ +/* in first or third quarter of unit circle.Routine receive also */ +/* (right argument) the original value of x for computing error of */ +/* result.And if result not accurate enough routine calls other routines */ +/***************************************************************************/ + +static double bsloww(double x,double dx, double orig,int n) { + static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ + double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn; + union {int4 i[2]; double x;} v; + x1=(x+th2_36)-th2_36; + y = aa.x*x1*x1*x1; + r=x+y; + x2=(x-x1)+dx; + xx=x*x; + t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + bb.x)*xx + 3.0*aa.x*x1*x2)*x +aa.x*x2*x2*x2+dx; + t=((x-r)+y)+t; + res=r+t; + cor = (r-res)+t; + 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); + 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); + } +} + +/***************************************************************************/ +/* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */ +/* in first or third quarter of unit circle.Routine receive also */ +/* (right argument) the original value of x for computing error of result.*/ +/* And if result not accurate enough routine calls other routines */ +/***************************************************************************/ + +static double bsloww1(double x, double dx, double orig,int n) { +mynumber u; + double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; + static const double t22 = 6291456.0; + int4 k; + y=ABS(x); + u.x=big.x+y; + y=y-(u.x-big.x); + dx=(x>0)?dx:-dx; + xx=y*y; + s = y*xx*(sn3 +xx*sn5); + c = xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + y1 = (y+t22)-t22; + y2 = (y - y1)+dx; + c1 = (cs+t22)-t22; + c2=(cs-c1)+ccs; + cor=(ssn+s*ccs+cs*s+c2*y+c1*y2-sn*y*dx)-sn*c; + y=sn+c1*y1; + cor = cor+((sn-y)+c1*y1); + res=y+cor; + cor=(y-res)+cor; + 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); + 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); + } +} + +/***************************************************************************/ +/* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */ +/* in second or fourth quarter of unit circle.Routine receive also the */ +/* original value and quarter(n= 1or 3)of x for computing error of result. */ +/* And if result not accurate enough routine calls other routines */ +/***************************************************************************/ + +static double bsloww2(double x, double dx, double orig, int n) { +mynumber u; + double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; + static const double t22 = 6291456.0; + int4 k; + y=ABS(x); + u.x=big.x+y; + y=y-(u.x-big.x); + dx=(x>0)?dx:-dx; + xx=y*y; + s = y*xx*(sn3 +xx*sn5); + c = y*dx+xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + + y1 = (y+t22)-t22; + y2 = (y - y1)+dx; + e1 = (sn+t22)-t22; + e2=(sn-e1)+ssn; + cor=(ccs-cs*c-e1*y2-e2*y)-sn*s; + y=cs-e1*y1; + cor = cor+((cs-y)-e1*y1); + res=y+cor; + cor=(y-res)+cor; + 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); + 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); + } +} + +/************************************************************************/ +/* Routine compute cos(x) for 2^-27 < |x|< 0.25 by Taylor with more */ +/* precision and if still doesn't accurate enough by mpcos or docos */ +/************************************************************************/ + +static double cslow2(double x) { + mynumber u; + double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; + static const double t22 = 6291456.0; + int4 k; + y=ABS(x); + u.x = big.x+y; + y = y-(u.x-big.x); + xx=y*y; + s = y*xx*(sn3 +xx*sn5); + c = xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + y1 = (y+t22)-t22; + y2 = y - y1; + e1 = (sn+t22)-t22; + e2=(sn-e1)+ssn; + cor=(ccs-cs*c-e1*y2-e2*y)-sn*s; + y=cs-e1*y1; + cor = cor+((cs-y)-e1*y1); + res=y+cor; + cor=(y-res)+cor; + if (res == res+1.0005*cor) + return res; + else { + y=ABS(x); + docos(y,0,w); + if (w[0] == w[0]+1.000000005*w[1]) return w[0]; + else return mpcos(x,0); + } +} + +/***************************************************************************/ +/* Routine compute cos(x+dx) (Double-Length number) where x is small enough*/ +/* to use Taylor series around zero and (x+dx) .Routine receive also */ +/* (right argument) the original value of x for computing error of */ +/* result.And if result not accurate enough routine calls other routines */ +/***************************************************************************/ + + +static double csloww(double x,double dx, double orig) { + static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ + double y,x1,x2,xx,r,t,res,cor,w[2],a,da,xn; + union {int4 i[2]; double x;} v; + int4 n; + x1=(x+th2_36)-th2_36; + y = aa.x*x1*x1*x1; + r=x+y; + x2=(x-x1)+dx; + xx=x*x; + /* Taylor series */ + t = (((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + bb.x)*xx + 3.0*aa.x*x1*x2)*x +aa.x*x2*x2*x2+dx; + t=((x-r)+y)+t; + res=r+t; + cor = (r-res)+t; + 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); + 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 { + t = (orig*hpinv.x + toint.x); + xn = t - toint.x; + v.x = t; + y = (orig - xn*mp1.x) - xn*mp2.x; + n =v.i[LOW_HALF]&3; + da = xn*pp3.x; + t=y-da; + da = (y-t)-da; + y = xn*pp4.x; + 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); + 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); + } + } +} + +/***************************************************************************/ +/* Routine compute sin(x+dx) (Double-Length number) where x in first or */ +/* third quarter of unit circle.Routine receive also (right argument) the */ +/* original value of x for computing error of result.And if result not */ +/* accurate enough routine calls other routines */ +/***************************************************************************/ + +static double csloww1(double x, double dx, double orig) { + mynumber u; + double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,c1,c2,xx,cor,res; + static const double t22 = 6291456.0; + int4 k; + y=ABS(x); + u.x=big.x+y; + y=y-(u.x-big.x); + dx=(x>0)?dx:-dx; + xx=y*y; + s = y*xx*(sn3 +xx*sn5); + c = xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + y1 = (y+t22)-t22; + y2 = (y - y1)+dx; + c1 = (cs+t22)-t22; + c2=(cs-c1)+ccs; + cor=(ssn+s*ccs+cs*s+c2*y+c1*y2-sn*y*dx)-sn*c; + y=sn+c1*y1; + cor = cor+((sn-y)+c1*y1); + res=y+cor; + cor=(y-res)+cor; + 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); + 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); + } +} + + +/***************************************************************************/ +/* Routine compute sin(x+dx) (Double-Length number) where x in second or */ +/* fourth quarter of unit circle.Routine receive also the original value */ +/* and quarter(n= 1or 3)of x for computing error of result.And if result not*/ +/* accurate enough routine calls other routines */ +/***************************************************************************/ + +static double csloww2(double x, double dx, double orig, int n) { + mynumber u; + double sn,ssn,cs,ccs,s,c,w[2],y,y1,y2,e1,e2,xx,cor,res; + static const double t22 = 6291456.0; + int4 k; + y=ABS(x); + u.x=big.x+y; + y=y-(u.x-big.x); + dx=(x>0)?dx:-dx; + xx=y*y; + s = y*xx*(sn3 +xx*sn5); + c = y*dx+xx*(cs2 +xx*(cs4 + xx*cs6)); + k=u.i[LOW_HALF]<<2; + sn=sincos.x[k]; + ssn=sincos.x[k+1]; + cs=sincos.x[k+2]; + ccs=sincos.x[k+3]; + + y1 = (y+t22)-t22; + y2 = (y - y1)+dx; + e1 = (sn+t22)-t22; + e2=(sn-e1)+ssn; + cor=(ccs-cs*c-e1*y2-e2*y)-sn*s; + y=cs-e1*y1; + cor = cor+((cs-y)-e1*y1); + res=y+cor; + cor=(y-res)+cor; + 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); + 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); + } +} + #ifdef NO_LONG_DOUBLE -strong_alias (__sin, __sinl) -weak_alias (__sin, sinl) +weak_alias (sin, sinl) +weak_alias (cos, cosl) #endif |