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
sin(double x)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