1 /*
2 * Double-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 #define _GNU_SOURCE
9 #include <math.h>
10 #include "mathlib.h"
11 #include "math_config.h"
12 #include "test_sig.h"
13 #include "test_defs.h"
14 #include "poly_scalar_f64.h"
15
16 /* Taylor series coefficents for sin(pi * x).
17 C2 coefficient (orginally ~=5.16771278) has been split into two parts:
18 C2_hi = 4, C2_lo = C2 - C2_hi (~=1.16771278)
19 This change in magnitude reduces floating point rounding errors.
20 C2_hi is then reintroduced after the polynomial approxmation. */
21 static const double poly[]
22 = { 0x1.921fb54442d184p1, -0x1.2aef39896f94bp0, 0x1.466bc6775ab16p1,
23 -0x1.32d2cce62dc33p-1, 0x1.507834891188ep-4, -0x1.e30750a28c88ep-8,
24 0x1.e8f48308acda4p-12, -0x1.6fc0032b3c29fp-16, 0x1.af86ae521260bp-21,
25 -0x1.012a9870eeb7dp-25 };
26
27 #define Shift 0x1.8p+52
28 /* TODO Store constant in structure for more efficient load. */
29 #define Pi 0x1.921fb54442d18p+1
30
31 /* Approximation for scalar double-precision sinpi(x).
32 Maximum error: 3.03 ULP:
33 sinpi(0x1.a90da2818f8b5p+7) got 0x1.fe358f255a4b3p-1
34 want 0x1.fe358f255a4b6p-1. */
35 double
arm_math_sinpi(double x)36 arm_math_sinpi (double x)
37 {
38 if (isinf (x) || isnan (x))
39 return __math_invalid (x);
40
41 double r = asdouble (asuint64 (x) & ~0x8000000000000000);
42 uint64_t sign = asuint64 (x) & 0x8000000000000000;
43
44 /* Edge cases for when sinpif should be exactly 0. (Integers)
45 0x1p53 is the limit for single precision to store any decimal places. */
46 if (r >= 0x1p53)
47 return asdouble (sign);
48
49 /* If x is an integer, return 0. */
50 uint64_t m = (uint64_t) r;
51 if (r == m)
52 return asdouble (sign);
53
54 /* For very small inputs, squaring r causes underflow.
55 Values below this threshold can be approximated via sinpi(x) ≈ pi*x. */
56 if (r < 0x1p-63)
57 return Pi * x;
58
59 /* Any non-integer values >= 0x1x51 will be int + 0.5.
60 These values should return exactly 1 or -1. */
61 if (r >= 0x1p51)
62 {
63 uint64_t iy = ((m & 1) << 63) ^ asuint64 (1.0);
64 return asdouble (sign ^ iy);
65 }
66
67 /* n = rint(|x|). */
68 double n = r + Shift;
69 sign ^= (asuint64 (n) << 63);
70 n = n - Shift;
71
72 /* r = |x| - n (range reduction into -1/2 .. 1/2). */
73 r = r - n;
74
75 /* y = sin(r). */
76 double r2 = r * r;
77 double y = horner_9_f64 (r2, poly);
78 y = y * r;
79
80 /* Reintroduce C2_hi. */
81 y = fma (-4 * r2, r, y);
82
83 /* Copy sign of x to sin(|x|). */
84 return asdouble (asuint64 (y) ^ sign);
85 }
86
87 #if WANT_EXPERIMENTAL_MATH
88 double
sinpi(double x)89 sinpi (double x)
90 {
91 return arm_math_sinpi (x);
92 }
93 #endif
94
95 #if WANT_TRIGPI_TESTS
96 TEST_ULP (arm_math_sinpi, 2.53)
97 TEST_SYM_INTERVAL (arm_math_sinpi, 0, 0x1p-63, 5000)
98 TEST_SYM_INTERVAL (arm_math_sinpi, 0x1p-63, 0.5, 10000)
99 TEST_SYM_INTERVAL (arm_math_sinpi, 0.5, 0x1p51, 10000)
100 TEST_SYM_INTERVAL (arm_math_sinpi, 0x1p51, inf, 10000)
101 #endif
102