1*05a0b428SJohn Marino /* @(#)s_sin.c 5.1 93/09/24 */ 2*05a0b428SJohn Marino /* 3*05a0b428SJohn Marino * ==================================================== 4*05a0b428SJohn Marino * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5*05a0b428SJohn Marino * 6*05a0b428SJohn Marino * Developed at SunPro, a Sun Microsystems, Inc. business. 7*05a0b428SJohn Marino * Permission to use, copy, modify, and distribute this 8*05a0b428SJohn Marino * software is freely granted, provided that this notice 9*05a0b428SJohn Marino * is preserved. 10*05a0b428SJohn Marino * ==================================================== 11*05a0b428SJohn Marino */ 12*05a0b428SJohn Marino 13*05a0b428SJohn Marino /* sin(x) 14*05a0b428SJohn Marino * Return sine function of x. 15*05a0b428SJohn Marino * 16*05a0b428SJohn Marino * kernel function: 17*05a0b428SJohn Marino * __kernel_sin ... sine function on [-pi/4,pi/4] 18*05a0b428SJohn Marino * __kernel_cos ... cose function on [-pi/4,pi/4] 19*05a0b428SJohn Marino * __ieee754_rem_pio2 ... argument reduction routine 20*05a0b428SJohn Marino * 21*05a0b428SJohn Marino * Method. 22*05a0b428SJohn Marino * Let S,C and T denote the sin, cos and tan respectively on 23*05a0b428SJohn Marino * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 24*05a0b428SJohn Marino * in [-pi/4 , +pi/4], and let n = k mod 4. 25*05a0b428SJohn Marino * We have 26*05a0b428SJohn Marino * 27*05a0b428SJohn Marino * n sin(x) cos(x) tan(x) 28*05a0b428SJohn Marino * ---------------------------------------------------------- 29*05a0b428SJohn Marino * 0 S C T 30*05a0b428SJohn Marino * 1 C -S -1/T 31*05a0b428SJohn Marino * 2 -S -C T 32*05a0b428SJohn Marino * 3 -C S -1/T 33*05a0b428SJohn Marino * ---------------------------------------------------------- 34*05a0b428SJohn Marino * 35*05a0b428SJohn Marino * Special cases: 36*05a0b428SJohn Marino * Let trig be any of sin, cos, or tan. 37*05a0b428SJohn Marino * trig(+-INF) is NaN, with signals; 38*05a0b428SJohn Marino * trig(NaN) is that NaN; 39*05a0b428SJohn Marino * 40*05a0b428SJohn Marino * Accuracy: 41*05a0b428SJohn Marino * TRIG(x) returns trig(x) nearly rounded 42*05a0b428SJohn Marino */ 43*05a0b428SJohn Marino 44*05a0b428SJohn Marino #include <float.h> 45*05a0b428SJohn Marino #include <math.h> 46*05a0b428SJohn Marino 47*05a0b428SJohn Marino #include "math_private.h" 48*05a0b428SJohn Marino 49*05a0b428SJohn Marino double 50*05a0b428SJohn Marino sin(double x) 51*05a0b428SJohn Marino { 52*05a0b428SJohn Marino double y[2],z=0.0; 53*05a0b428SJohn Marino int32_t n, ix; 54*05a0b428SJohn Marino 55*05a0b428SJohn Marino /* High word of x. */ 56*05a0b428SJohn Marino GET_HIGH_WORD(ix,x); 57*05a0b428SJohn Marino 58*05a0b428SJohn Marino /* |x| ~< pi/4 */ 59*05a0b428SJohn Marino ix &= 0x7fffffff; 60*05a0b428SJohn Marino if(ix <= 0x3fe921fb) return __kernel_sin(x,z,0); 61*05a0b428SJohn Marino 62*05a0b428SJohn Marino /* sin(Inf or NaN) is NaN */ 63*05a0b428SJohn Marino else if (ix>=0x7ff00000) return x-x; 64*05a0b428SJohn Marino 65*05a0b428SJohn Marino /* argument reduction needed */ 66*05a0b428SJohn Marino else { 67*05a0b428SJohn Marino n = __ieee754_rem_pio2(x,y); 68*05a0b428SJohn Marino switch(n&3) { 69*05a0b428SJohn Marino case 0: return __kernel_sin(y[0],y[1],1); 70*05a0b428SJohn Marino case 1: return __kernel_cos(y[0],y[1]); 71*05a0b428SJohn Marino case 2: return -__kernel_sin(y[0],y[1],1); 72*05a0b428SJohn Marino default: 73*05a0b428SJohn Marino return -__kernel_cos(y[0],y[1]); 74*05a0b428SJohn Marino } 75*05a0b428SJohn Marino } 76*05a0b428SJohn Marino } 77*05a0b428SJohn Marino 78*05a0b428SJohn Marino #if LDBL_MANT_DIG == DBL_MANT_DIG 79*05a0b428SJohn Marino __strong_alias(sinl, sin); 80*05a0b428SJohn Marino #endif /* LDBL_MANT_DIG == DBL_MANT_DIG */ 81