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 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 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