1*31930Szliu /* 2*31930Szliu * Copyright (c) 1987 Regents of the University of California. 3*31930Szliu * 4*31930Szliu * Use and reproduction of this software are granted in accordance with 5*31930Szliu * the terms and conditions specified in the Berkeley Software License 6*31930Szliu * Agreement (in particular, this entails acknowledgement of the programs' 7*31930Szliu * source, and inclusion of this notice) with the additional understanding 8*31930Szliu * that all recipients should regard themselves as participants in an 9*31930Szliu * ongoing research project and hence should feel obligated to report 10*31930Szliu * their experiences (good or bad) with these elementary function codes, 11*31930Szliu * using "sendbug 4bsd-bugs@BERKELEY", to the authors. 12*31930Szliu */ 13*31930Szliu 14*31930Szliu #ifndef lint 15*31930Szliu static char sccsid[] = "@(#)tan.c 1.1 1.1 (ucb.elefunt) 07/24/87"; 16*31930Szliu #endif /* not lint */ 17*31930Szliu 18*31930Szliu #include "trig.h" 19*31930Szliu double 20*31930Szliu tan(x) 21*31930Szliu double x; 22*31930Szliu { 23*31930Szliu double a,z,ss,cc,c; 24*31930Szliu int k; 25*31930Szliu 26*31930Szliu if(!finite(x)) /* tan(NaN) and tan(INF) must be NaN */ 27*31930Szliu return x-x; 28*31930Szliu x = drem(x,PI); /* reduce x into [-PI/2, PI/2] */ 29*31930Szliu a = copysign(x,one); /* ... = abs(x) */ 30*31930Szliu if (a >= PIo4) { 31*31930Szliu k = 1; 32*31930Szliu x = copysign(PIo2-a,x); 33*31930Szliu } 34*31930Szliu else { 35*31930Szliu k = 0; 36*31930Szliu if (a < small) { 37*31930Szliu big+a; 38*31930Szliu return x; 39*31930Szliu } 40*31930Szliu } 41*31930Szliu z = x*x; 42*31930Szliu cc = cos__C(z); 43*31930Szliu ss = sin__S(z); 44*31930Szliu z *= half; /* Next get c = cos(x) accurately */ 45*31930Szliu c = (z >= thresh ? half-((z-half)-cc) : one-(z-cc)); 46*31930Szliu if (k == 0) 47*31930Szliu return x+(x*(z-(cc-ss)))/c; /* ... sin/cos */ 48*31930Szliu #ifdef national 49*31930Szliu else if (x == zero) 50*31930Szliu return copysign(fmax,x); /* no inf on 32k */ 51*31930Szliu #endif /* national */ 52*31930Szliu else 53*31930Szliu return c/(x+x*ss); /* ... cos/sin */ 54*31930Szliu } 55