1 /*
2 * Single-precision scalar tanpi(x) function.
3 *
4 * Copyright (c) 2024, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7 #include "mathlib.h"
8 #include "math_config.h"
9 #include "test_sig.h"
10 #include "test_defs.h"
11 #include "poly_scalar_f32.h"
12
13 const static struct tanpif_data
14 {
15 float tan_poly[6], cot_poly[4], pi, invpi;
16 } tanpif_data = {
17 /* Coefficents for tan(pi * x). */
18 .tan_poly = {
19 0x1.4abbc8p3,
20 0x1.467284p5,
21 0x1.44cf12p7,
22 0x1.596b5p9,
23 0x1.753858p10,
24 0x1.76ff52p14,
25 },
26 /* Coefficents for cot(pi * x). */
27 .cot_poly = {
28 -0x1.0c1522p0,
29 -0x1.60ce32p-1,
30 -0x1.49cd42p-1,
31 -0x1.73f786p-1,
32 },
33 .pi = 0x1.921fb6p1f,
34 .invpi = 0x1.45f308p-2f,
35 };
36
37 /* Single-precision scalar tanpi(x) implementation.
38 Maximum error 2.56 ULP:
39 tanpif(0x1.4bf948p-1) got -0x1.fcc9ep+0
40 want -0x1.fcc9e6p+0. */
41 float
arm_math_tanpif(float x)42 arm_math_tanpif (float x)
43 {
44 uint32_t xabs_12 = asuint (x) >> 20 & 0x7f8;
45
46 /* x >= 0x1p24f. */
47 if (unlikely (xabs_12 >= 0x4b1))
48 {
49 /* tanpif(+/-inf) and tanpif(+/-nan) = nan. */
50 if (unlikely (xabs_12 == 0x7f8))
51 {
52 return __math_invalidf (x);
53 }
54
55 uint32_t x_sign = asuint (x) & 0x80000000;
56 return asfloat (x_sign);
57 }
58
59 const struct tanpif_data *d = ptr_barrier (&tanpif_data);
60
61 /* Prevent underflow exceptions. x <= 0x1p-31. */
62 if (unlikely (xabs_12 < 0x300))
63 {
64 return d->pi * x;
65 }
66
67 float rounded = roundf (x);
68 if (unlikely (rounded == x))
69 {
70 /* If x == 0, return with sign. */
71 if (x == 0)
72 {
73 return x;
74 }
75 /* Otherwise, return zero with alternating sign. */
76 int32_t m = (int32_t) rounded;
77 if (x < 0)
78 {
79 return m & 1 ? 0.0f : -0.0f;
80 }
81 else
82 {
83 return m & 1 ? -0.0f : 0.0f;
84 }
85 }
86
87 float x_reduced = x - rounded;
88 float abs_x_reduced = 0.5f - asfloat (asuint (x_reduced) & 0x7fffffff);
89
90 float result, offset, scale;
91
92 /* Test 0.25 < abs_x < 0.5 independent from abs_x_reduced. */
93 float x2 = x + x;
94 int32_t rounded_x2 = (int32_t) roundf (x2);
95 if (rounded_x2 & 1)
96 {
97 float r_x = abs_x_reduced;
98
99 float r_x2 = r_x * r_x;
100 float r_x4 = r_x2 * r_x2;
101
102 uint32_t sign = asuint (x_reduced) & 0x80000000;
103 r_x = asfloat (asuint (r_x) ^ sign);
104
105 // calculate sign for half-fractional inf values
106 uint32_t is_finite = asuint (abs_x_reduced);
107 uint32_t is_odd = (rounded_x2 & 2) << 30;
108 uint32_t is_neg = rounded_x2 & 0x80000000;
109 uint32_t keep_sign = is_finite | (is_odd ^ is_neg);
110 offset = d->invpi / (keep_sign ? r_x : -r_x);
111 scale = r_x;
112
113 result = pairwise_poly_3_f32 (r_x2, r_x4, d->cot_poly);
114 }
115 else
116 {
117 float r_x = x_reduced;
118
119 float r_x2 = r_x * r_x;
120 float r_x4 = r_x2 * r_x2;
121
122 offset = d->pi * r_x;
123 scale = r_x * r_x2;
124
125 result = pw_horner_5_f32 (r_x2, r_x4, d->tan_poly);
126 }
127
128 return fmaf (scale, result, offset);
129 }
130
131 #if WANT_EXPERIMENTAL_MATH
132 float
tanpif(float x)133 tanpif (float x)
134 {
135 return arm_math_tanpif (x);
136 }
137 #endif
138
139 #if WANT_TRIGPI_TESTS
140 TEST_ULP (arm_math_tanpif, 2.57)
141 TEST_SYM_INTERVAL (arm_math_tanpif, 0, 0x1p-31f, 50000)
142 TEST_SYM_INTERVAL (arm_math_tanpif, 0x1p-31f, 0.5, 100000)
143 TEST_SYM_INTERVAL (arm_math_tanpif, 0.5, 0x1p23f, 100000)
144 TEST_SYM_INTERVAL (arm_math_tanpif, 0x1p23f, inf, 100000)
145 #endif
146