1/*-
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#if defined(LIBM_SCCS) && !defined(lint)
17__RCSID("$NetBSD: s_sincos.c,v 1.5 2022/08/29 01:48:34 riastradh Exp $");
18#endif
19#if 0
20__FBSDID("$FreeBSD: head/lib/msun/src/s_sincos.c 319047 2017-05-28 06:13:38Z mmel $");
21#endif
22
23#include "namespace.h"
24#include <float.h>
25
26#include "math.h"
27// #define INLINE_REM_PIO2
28#include "math_private.h"
29// #include "e_rem_pio2.c"
30#include "k_sincos.h"
31
32#ifdef __weak_alias
33__weak_alias(sincos,_sincos)
34#endif
35
36void
37sincos(double x, double *sn, double *cs)
38{
39	double y[2];
40	int32_t n, ix;
41
42	/* High word of x. */
43	GET_HIGH_WORD(ix, x);
44
45	/* |x| ~< pi/4 */
46	ix &= 0x7fffffff;
47	if (ix <= 0x3fe921fb) {
48		if (ix < 0x3e400000) {		/* |x| < 2**-27 */
49			if ((int)x == 0) {	/* Generate inexact. */
50				*sn = x;
51				*cs = 1;
52				return;
53			}
54		}
55		__kernel_sincos(x, 0, 0, sn, cs);
56		return;
57	}
58
59	/* If x = Inf or NaN, then sin(x) = NaN and cos(x) = NaN. */
60	if (ix >= 0x7ff00000) {
61		*sn = x - x;
62		*cs = x - x;
63		return;
64	}
65
66	/* Argument reduction. */
67	n = __ieee754_rem_pio2(x, y);
68
69	switch(n & 3) {
70	case 0:
71		__kernel_sincos(y[0], y[1], 1, sn, cs);
72		break;
73	case 1:
74		__kernel_sincos(y[0], y[1], 1, cs, sn);
75		*cs = -*cs;
76		break;
77	case 2:
78		__kernel_sincos(y[0], y[1], 1, sn, cs);
79		*sn = -*sn;
80		*cs = -*cs;
81		break;
82	default:
83		__kernel_sincos(y[0], y[1], 1, cs, sn);
84		*sn = -*sn;
85	}
86}
87
88#if !defined(__HAVE_LONG_DOUBLE)
89__weak_alias(sincosl, sincos);
90#endif
91