1*412f47f9SXin Li /*
2*412f47f9SXin Li * Single-precision sin/cos function.
3*412f47f9SXin Li *
4*412f47f9SXin Li * Copyright (c) 2018-2021, Arm Limited.
5*412f47f9SXin Li * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6*412f47f9SXin Li */
7*412f47f9SXin Li
8*412f47f9SXin Li #include <stdint.h>
9*412f47f9SXin Li #include <math.h>
10*412f47f9SXin Li #include "math_config.h"
11*412f47f9SXin Li #include "sincosf.h"
12*412f47f9SXin Li
13*412f47f9SXin Li /* Fast sincosf implementation. Worst-case ULP is 0.5607, maximum relative
14*412f47f9SXin Li error is 0.5303 * 2^-23. A single-step range reduction is used for
15*412f47f9SXin Li small values. Large inputs have their range reduced using fast integer
16*412f47f9SXin Li arithmetic. */
17*412f47f9SXin Li void
sincosf(float y,float * sinp,float * cosp)18*412f47f9SXin Li sincosf (float y, float *sinp, float *cosp)
19*412f47f9SXin Li {
20*412f47f9SXin Li double x = y;
21*412f47f9SXin Li double s;
22*412f47f9SXin Li int n;
23*412f47f9SXin Li const sincos_t *p = &__sincosf_table[0];
24*412f47f9SXin Li
25*412f47f9SXin Li if (abstop12 (y) < abstop12 (pio4f))
26*412f47f9SXin Li {
27*412f47f9SXin Li double x2 = x * x;
28*412f47f9SXin Li
29*412f47f9SXin Li if (unlikely (abstop12 (y) < abstop12 (0x1p-12f)))
30*412f47f9SXin Li {
31*412f47f9SXin Li if (unlikely (abstop12 (y) < abstop12 (0x1p-126f)))
32*412f47f9SXin Li /* Force underflow for tiny y. */
33*412f47f9SXin Li force_eval_float (x2);
34*412f47f9SXin Li *sinp = y;
35*412f47f9SXin Li *cosp = 1.0f;
36*412f47f9SXin Li return;
37*412f47f9SXin Li }
38*412f47f9SXin Li
39*412f47f9SXin Li sincosf_poly (x, x2, p, 0, sinp, cosp);
40*412f47f9SXin Li }
41*412f47f9SXin Li else if (abstop12 (y) < abstop12 (120.0f))
42*412f47f9SXin Li {
43*412f47f9SXin Li x = reduce_fast (x, p, &n);
44*412f47f9SXin Li
45*412f47f9SXin Li /* Setup the signs for sin and cos. */
46*412f47f9SXin Li s = p->sign[n & 3];
47*412f47f9SXin Li
48*412f47f9SXin Li if (n & 2)
49*412f47f9SXin Li p = &__sincosf_table[1];
50*412f47f9SXin Li
51*412f47f9SXin Li sincosf_poly (x * s, x * x, p, n, sinp, cosp);
52*412f47f9SXin Li }
53*412f47f9SXin Li else if (likely (abstop12 (y) < abstop12 (INFINITY)))
54*412f47f9SXin Li {
55*412f47f9SXin Li uint32_t xi = asuint (y);
56*412f47f9SXin Li int sign = xi >> 31;
57*412f47f9SXin Li
58*412f47f9SXin Li x = reduce_large (xi, &n);
59*412f47f9SXin Li
60*412f47f9SXin Li /* Setup signs for sin and cos - include original sign. */
61*412f47f9SXin Li s = p->sign[(n + sign) & 3];
62*412f47f9SXin Li
63*412f47f9SXin Li if ((n + sign) & 2)
64*412f47f9SXin Li p = &__sincosf_table[1];
65*412f47f9SXin Li
66*412f47f9SXin Li sincosf_poly (x * s, x * x, p, n, sinp, cosp);
67*412f47f9SXin Li }
68*412f47f9SXin Li else
69*412f47f9SXin Li {
70*412f47f9SXin Li /* Return NaN if Inf or NaN for both sin and cos. */
71*412f47f9SXin Li *sinp = *cosp = y - y;
72*412f47f9SXin Li #if WANT_ERRNO
73*412f47f9SXin Li /* Needed to set errno for +-Inf, the add is a hack to work
74*412f47f9SXin Li around a gcc register allocation issue: just passing y
75*412f47f9SXin Li affects code generation in the fast path. */
76*412f47f9SXin Li __math_invalidf (y + y);
77*412f47f9SXin Li #endif
78*412f47f9SXin Li }
79*412f47f9SXin Li }
80