POK
sin.c
1 /*
2  * POK header
3  *
4  * The following file is a part of the POK project. Any modification should
5  * made according to the POK licence. You CANNOT use this file or a part of
6  * this file is this part of a file for your own project
7  *
8  * For more information on the POK licence, please see our LICENCE FILE
9  *
10  * Please follow the coding guidelines described in doc/CODING_GUIDELINES
11  *
12  * Copyright (c) 2007-2009 POK team
13  *
14  * Created by julien on Fri Jan 30 14:41:34 2009
15  */
16 
17 /* @(#)s_sin.c 5.1 93/09/24 */
18 /*
19  * ====================================================
20  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
21  *
22  * Developed at SunPro, a Sun Microsystems, Inc. business.
23  * Permission to use, copy, modify, and distribute this
24  * software is freely granted, provided that this notice
25  * is preserved.
26  * ====================================================
27  */
28 
29 /* sin(x)
30  * Return sine function of x.
31  *
32  * kernel function:
33  * __kernel_sin ... sine function on [-pi/4,pi/4]
34  * __kernel_cos ... cose function on [-pi/4,pi/4]
35  * __ieee754_rem_pio2 ... argument reduction routine
36  *
37  * Method.
38  * Let S,C and T denote the sin, cos and tan respectively on
39  * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
40  * in [-pi/4 , +pi/4], and let n = k mod 4.
41  * We have
42  *
43  * n sin(x) cos(x) tan(x)
44  * ----------------------------------------------------------
45  * 0 S C T
46  * 1 C -S -1/T
47  * 2 -S -C T
48  * 3 -C S -1/T
49  * ----------------------------------------------------------
50  *
51  * Special cases:
52  * Let trig be any of sin, cos, or tan.
53  * trig(+-INF) is NaN, with signals;
54  * trig(NaN) is that NaN;
55  *
56  * Accuracy:
57  * TRIG(x) returns trig(x) nearly rounded
58  */
59 
60 #ifdef POK_NEEDS_LIBMATH
61 
62 #include <libm.h>
63 #include "namespace.h"
64 #include "math_private.h"
65 
66 #if 0 /* notyet */
67 #ifdef __weak_alias
68 __weak_alias(sin, _sin)
69 #endif
70 #endif
71 
72 double
73 sin(double x)
74 {
75  double y[2],z=0.0;
76  int32_t n, ix;
77 
78  /* High word of x. */
79  GET_HIGH_WORD(ix,x);
80 
81  /* |x| ~< pi/4 */
82  ix &= 0x7fffffff;
83  if(ix <= 0x3fe921fb) return __kernel_sin(x,z,0);
84 
85  /* sin(Inf or NaN) is NaN */
86  else if (ix>=0x7ff00000) return x-x;
87 
88  /* argument reduction needed */
89  else {
90  n = __ieee754_rem_pio2(x,y);
91  switch(n&3) {
92  case 0: return __kernel_sin(y[0],y[1],1);
93  case 1: return __kernel_cos(y[0],y[1]);
94  case 2: return -__kernel_sin(y[0],y[1],1);
95  default:
96  return -__kernel_cos(y[0],y[1]);
97  }
98  }
99 }
100 
101 #endif