1 /*
2 * Single-precision scalar sinpi function.
3 *
4 * Copyright (c) 2023-2024, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8 #include "mathlib.h"
9 #include "math_config.h"
10 #include "test_sig.h"
11 #include "test_defs.h"
12
13 /* Taylor series coefficents for sin(pi * x). */
14 #define C0 0x1.921fb6p1f
15 #define C1 -0x1.4abbcep2f
16 #define C2 0x1.466bc6p1f
17 #define C3 -0x1.32d2ccp-1f
18 #define C4 0x1.50783p-4f
19 #define C5 -0x1.e30750p-8f
20
21 #define Shift 0x1.0p+23f
22
23 /* Approximation for scalar single-precision sinpi(x) - sinpif.
24 Maximum error: 2.48 ULP:
25 sinpif(0x1.d062b6p-2) got 0x1.fa8c06p-1
26 want 0x1.fa8c02p-1. */
27 float
arm_math_sinpif(float x)28 arm_math_sinpif (float x)
29 {
30 if (isinf (x) || isnan (x))
31 return __math_invalidf (x);
32
33 float r = asfloat (asuint (x) & ~0x80000000);
34 uint32_t sign = asuint (x) & 0x80000000;
35
36 /* Edge cases for when sinpif should be exactly 0. (Integers)
37 0x1p23 is the limit for single precision to store any decimal places. */
38 if (r >= 0x1p23f)
39 return asfloat (sign);
40
41 int32_t m = roundf (r);
42 if (m == r)
43 return asfloat (sign);
44
45 /* For very small inputs, squaring r causes underflow.
46 Values below this threshold can be approximated via sinpi(x) ~= pi*x. */
47 if (r < 0x1p-31f)
48 return C0 * x;
49
50 /* Any non-integer values >= 0x1p22f will be int + 0.5.
51 These values should return exactly 1 or -1. */
52 if (r >= 0x1p22f)
53 {
54 uint32_t iy = ((m & 1) << 31) ^ asuint (-1.0f);
55 return asfloat (sign ^ iy);
56 }
57
58 /* n = rint(|x|). */
59 float n = r + Shift;
60 sign ^= (asuint (n) << 31);
61 n = n - Shift;
62
63 /* r = |x| - n (range reduction into -1/2 .. 1/2). */
64 r = r - n;
65
66 /* y = sin(pi * r). */
67 float r2 = r * r;
68 float y = fmaf (C5, r2, C4);
69 y = fmaf (y, r2, C3);
70 y = fmaf (y, r2, C2);
71 y = fmaf (y, r2, C1);
72 y = fmaf (y, r2, C0);
73
74 /* Copy sign of x to sin(|x|). */
75 return asfloat (asuint (y * r) ^ sign);
76 }
77
78 #if WANT_EXPERIMENTAL_MATH
79 float
sinpif(float x)80 sinpif (float x)
81 {
82 return arm_math_sinpif (x);
83 }
84 #endif
85
86 #if WANT_TRIGPI_TESTS
87 TEST_ULP (arm_math_sinpif, 1.99)
88 TEST_SYM_INTERVAL (arm_math_sinpif, 0, 0x1p-31, 5000)
89 TEST_SYM_INTERVAL (arm_math_sinpif, 0x1p-31, 0.5, 10000)
90 TEST_SYM_INTERVAL (arm_math_sinpif, 0.5, 0x1p22f, 10000)
91 TEST_SYM_INTERVAL (arm_math_sinpif, 0x1p22f, inf, 10000)
92 #endif
93