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