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