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