blob: 72b2d64da329233d58294f9a5cd330a5b72c2a90 [file] [log] [blame]
Elliott Hughese332f652018-05-08 15:07:43 -07001/*-
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 *
11 * s_sin.c and s_cos.c merged by Steven G. Kargl. Descriptions of the
12 * algorithms are contained in the original files.
13 */
14
15#include <sys/cdefs.h>
16__FBSDID("$FreeBSD: head/lib/msun/src/s_sincos.c 319047 2017-05-28 06:13:38Z mmel $");
17
18#include <float.h>
19
20#include "math.h"
21#define INLINE_REM_PIO2
22#include "math_private.h"
23#include "e_rem_pio2.c"
24#include "k_sincos.h"
25
26void
27sincos(double x, double *sn, double *cs)
28{
29 double y[2];
30 int32_t n, ix;
31
32 /* High word of x. */
33 GET_HIGH_WORD(ix, x);
34
35 /* |x| ~< pi/4 */
36 ix &= 0x7fffffff;
37 if (ix <= 0x3fe921fb) {
38 if (ix < 0x3e400000) { /* |x| < 2**-27 */
39 if ((int)x == 0) { /* Generate inexact. */
40 *sn = x;
41 *cs = 1;
42 return;
43 }
44 }
45 __kernel_sincos(x, 0, 0, sn, cs);
46 return;
47 }
48
49 /* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
50 if (ix >= 0x7ff00000) {
51 *sn = x - x;
52 *cs = x - x;
53 return;
54 }
55
56 /* Argument reduction. */
57 n = __ieee754_rem_pio2(x, y);
58
59 switch(n & 3) {
60 case 0:
61 __kernel_sincos(y[0], y[1], 1, sn, cs);
62 break;
63 case 1:
64 __kernel_sincos(y[0], y[1], 1, cs, sn);
65 *cs = -*cs;
66 break;
67 case 2:
68 __kernel_sincos(y[0], y[1], 1, sn, cs);
69 *sn = -*sn;
70 *cs = -*cs;
71 break;
72 default:
73 __kernel_sincos(y[0], y[1], 1, cs, sn);
74 *sn = -*sn;
75 }
76}
77
78#if (LDBL_MANT_DIG == 53)
79__weak_reference(sincos, sincosl);
80#endif