gonzui

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