aboutsummaryrefslogtreecommitdiff
path: root/sysdeps/ieee754/dbl-64/s_sin.c
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754/dbl-64/s_sin.c')
-rw-r--r--sysdeps/ieee754/dbl-64/s_sin.c1178
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
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 02:56:46 +0000