1 /*
2 * Single-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
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 cospi(x) - cospif.
24 Maximum error: 2.64 ULP:
25 cospif(0x1.37e844p-4) got 0x1.f16b3p-1
26 want 0x1.f16b2ap-1. */
27 float
arm_math_cospif(float x)28 arm_math_cospif (float x)
29 {
30 if (isinf (x) || isnan (x))
31 return __math_invalidf (x);
32
33 float ax = asfloat (asuint (x) & ~0x80000000);
34
35 /* Edge cases for when cospif should be exactly +/- 1. (Integers)
36 0x1p23 is the limit for single precision to store any decimal places. */
37 if (ax >= 0x1p24f)
38 return 1;
39
40 uint32_t m = roundf (ax);
41 if (m == ax)
42 return (m & 1) ? -1 : 1;
43
44 /* Any non-integer values >= 0x1p22f will be int +0.5.
45 These values should return exactly 0. */
46 if (ax >= 0x1p22f)
47 return 0;
48
49 /* For very small inputs, squaring r causes underflow.
50 Values below this threshold can be approximated via cospi(x) ~= 1 -
51 (pi*x). */
52 if (ax < 0x1p-31f)
53 return 1 - (C0 * x);
54
55 /* n = rint(|x|). */
56 float n = ax + Shift;
57 uint32_t sign = asuint (n) << 31;
58 n = n - Shift;
59
60 /* We know that cospi(x) = sinpi(0.5 - x)
61 range reduction and offset into sinpi range -1/2 .. 1/2
62 r = 0.5 - |x - rint(x)|. */
63 float r = 0.5f - fabs (ax - n);
64
65 /* y = sin(pi * r). */
66 float r2 = r * r;
67 float y = fmaf (C5, r2, C4);
68 y = fmaf (y, r2, C3);
69 y = fmaf (y, r2, C2);
70 y = fmaf (y, r2, C1);
71 y = fmaf (y, r2, C0);
72
73 /* As all values are reduced to -1/2 .. 1/2, the result of cos(x) always be
74 positive, therefore, the sign must be introduced based upon if x rounds to
75 odd or even. */
76 return asfloat (asuint (y * r) ^ sign);
77 }
78
79 #if WANT_EXPERIMENTAL_MATH
80 float
cospif(float x)81 cospif (float x)
82 {
83 return arm_math_cospif (x);
84 }
85 #endif
86
87 #if WANT_TRIGPI_TESTS
88 TEST_ULP (arm_math_cospif, 2.15)
89 TEST_SYM_INTERVAL (arm_math_cospif, 0, 0x1p-31, 5000)
90 TEST_SYM_INTERVAL (arm_math_cospif, 0x1p-31, 0.5, 10000)
91 TEST_SYM_INTERVAL (arm_math_cospif, 0.5, 0x1p22f, 10000)
92 TEST_SYM_INTERVAL (arm_math_cospif, 0x1p22f, inf, 10000)
93 #endif
94