1*31931Szliu /* 2*31931Szliu * Copyright (c) 1987 Regents of the University of California. 3*31931Szliu * 4*31931Szliu * Use and reproduction of this software are granted in accordance with 5*31931Szliu * the terms and conditions specified in the Berkeley Software License 6*31931Szliu * Agreement (in particular, this entails acknowledgement of the programs' 7*31931Szliu * source, and inclusion of this notice) with the additional understanding 8*31931Szliu * that all recipients should regard themselves as participants in an 9*31931Szliu * ongoing research project and hence should feel obligated to report 10*31931Szliu * their experiences (good or bad) with these elementary function codes, 11*31931Szliu * using "sendbug 4bsd-bugs@BERKELEY", to the authors. 12*31931Szliu */ 13*31931Szliu 14*31931Szliu #ifndef lint 15*31931Szliu static char sccsid[] = "@(#)sincos.c 1.1 1.1 (ucb.elefunt) 07/24/87"; 16*31931Szliu #endif /* not lint */ 17*31931Szliu 18*31931Szliu #include "trig.h" 19*31931Szliu double 20*31931Szliu sin(x) 21*31931Szliu double x; 22*31931Szliu { 23*31931Szliu double a,c,z; 24*31931Szliu 25*31931Szliu if(!finite(x)) /* sin(NaN) and sin(INF) must be NaN */ 26*31931Szliu return x-x; 27*31931Szliu x=drem(x,PI2); /* reduce x into [-PI,PI] */ 28*31931Szliu a=copysign(x,one); 29*31931Szliu if (a >= PIo4) { 30*31931Szliu if(a >= PI3o4) /* ... in [3PI/4,PI] */ 31*31931Szliu x = copysign((a = PI-a),x); 32*31931Szliu else { /* ... in [PI/4,3PI/4] */ 33*31931Szliu a = PIo2-a; /* rtn. sign(x)*C(PI/2-|x|) */ 34*31931Szliu z = a*a; 35*31931Szliu c = cos__C(z); 36*31931Szliu z *= half; 37*31931Szliu a = (z >= thresh ? half-((z-half)-c) : one-(z-c)); 38*31931Szliu return copysign(a,x); 39*31931Szliu } 40*31931Szliu } 41*31931Szliu 42*31931Szliu if (a < small) { /* rtn. S(x) */ 43*31931Szliu big+a; 44*31931Szliu return x; 45*31931Szliu } 46*31931Szliu return x+x*sin__S(x*x); 47*31931Szliu } 48*31931Szliu 49*31931Szliu double 50*31931Szliu cos(x) 51*31931Szliu double x; 52*31931Szliu { 53*31931Szliu double a,c,z,s = 1.0; 54*31931Szliu 55*31931Szliu if(!finite(x)) /* cos(NaN) and cos(INF) must be NaN */ 56*31931Szliu return x-x; 57*31931Szliu x=drem(x,PI2); /* reduce x into [-PI,PI] */ 58*31931Szliu a=copysign(x,one); 59*31931Szliu if (a >= PIo4) { 60*31931Szliu if (a >= PI3o4) { /* ... in [3PI/4,PI] */ 61*31931Szliu a = PI-a; 62*31931Szliu s = negone; 63*31931Szliu } 64*31931Szliu else { /* ... in [PI/4,3PI/4] */ 65*31931Szliu a = PIo2-a; 66*31931Szliu return a+a*sin__S(a*a); /* rtn. S(PI/2-|x|) */ 67*31931Szliu } 68*31931Szliu } 69*31931Szliu if (a < small) { 70*31931Szliu big+a; 71*31931Szliu return s; /* rtn. s*C(a) */ 72*31931Szliu } 73*31931Szliu z = a*a; 74*31931Szliu c = cos__C(z); 75*31931Szliu z *= half; 76*31931Szliu a = (z >= thresh ? half-((z-half)-c) : one-(z-c)); 77*31931Szliu return copysign(a,s); 78*31931Szliu } 79