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