1 /* 2 * Double-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_f64.h" 12 13 #define SIGN_MASK 0x8000000000000000 14 15 const static struct tanpi_data 16 { 17 double tan_poly[14], cot_poly[9], pi, invpi; 18 } tanpi_data = { 19 /* Coefficents for tan(pi * x). */ 20 .tan_poly = { 21 0x1.4abbce625be52p3, 22 0x1.466bc6775b0f9p5, 23 0x1.45fff9b426f5ep7, 24 0x1.45f4730dbca5cp9, 25 0x1.45f3265994f85p11, 26 0x1.45f4234b330cap13, 27 0x1.45dca11be79ebp15, 28 0x1.47283fc5eea69p17, 29 0x1.3a6d958cdefaep19, 30 0x1.927896baee627p21, 31 -0x1.89333f6acd922p19, 32 0x1.5d4e912bb8456p27, 33 -0x1.a854d53ab6874p29, 34 0x1.1b76de7681424p32, 35 }, 36 /* Coefficents for cot(pi * x). */ 37 .cot_poly = { 38 -0x1.0c152382d7366p0, 39 -0x1.60c8539c1d316p-1, 40 -0x1.4b9a2f3516354p-1, 41 -0x1.47474060b6ba8p-1, 42 -0x1.464633ad9dcb1p-1, 43 -0x1.45ff229d7edd6p-1, 44 -0x1.46d8dbf492923p-1, 45 -0x1.3873892311c6bp-1, 46 -0x1.b2f3d0ff96d73p-1, 47 }, 48 .pi = 0x1.921fb54442d18p1, 49 .invpi = 0x1.45f306dc9c883p-2, 50 }; 51 52 /* Double-precision scalar tanpi(x) implementation. 53 Maximum error 2.19 ULP: 54 tanpi(0x1.68847e177a855p-2) got 0x1.fe9a0ff9bb9d7p+0 55 want 0x1.fe9a0ff9bb9d5p+0. */ 56 double 57 arm_math_tanpi (double x) 58 { 59 uint64_t xabs_12 = asuint64 (x) >> 52 & 0x7ff; 60 61 /* x >= 0x1p54. */ 62 if (unlikely (xabs_12 >= 0x434)) 63 { 64 /* tanpi(+/-inf) and tanpi(+/-nan) = nan. */ 65 if (unlikely (xabs_12 == 0x7ff)) 66 { 67 return __math_invalid (x); 68 } 69 70 uint64_t x_sign = asuint64 (x) & SIGN_MASK; 71 return asdouble (x_sign); 72 } 73 74 const struct tanpi_data *d = ptr_barrier (&tanpi_data); 75 76 double rounded = round (x); 77 if (unlikely (rounded == x)) 78 { 79 /* If x == 0, return with sign. */ 80 if (x == 0) 81 { 82 return x; 83 } 84 /* Otherwise, return zero with alternating sign. */ 85 int64_t m = (int64_t) rounded; 86 if (x < 0) 87 { 88 return m & 1 ? 0.0 : -0.0; 89 } 90 else 91 { 92 return m & 1 ? -0.0 : 0.0; 93 } 94 } 95 96 double x_reduced = x - rounded; 97 double abs_x_reduced = 0.5 - fabs (x_reduced); 98 99 /* Prevent underflow exceptions. x <= 0x1p-63. */ 100 if (unlikely (xabs_12 < 0x3c0)) 101 { 102 return d->pi * x; 103 } 104 105 double result, offset, scale; 106 107 /* Test 0.25 < abs_x < 0.5 independent from abs_x_reduced. */ 108 double x2 = x + x; 109 int64_t rounded_x2 = (int64_t) round (x2); 110 if (rounded_x2 & 1) 111 { 112 double r_x = abs_x_reduced; 113 114 double r_x2 = r_x * r_x; 115 double r_x4 = r_x2 * r_x2; 116 117 uint64_t sign = asuint64 (x_reduced) & SIGN_MASK; 118 r_x = asdouble (asuint64 (r_x) ^ sign); 119 120 // calculate sign for half-fractional inf values 121 uint64_t is_finite = asuint64 (abs_x_reduced); 122 uint64_t is_odd = (rounded_x2 & 2) << 62; 123 uint64_t is_neg = rounded_x2 & SIGN_MASK; 124 uint64_t keep_sign = is_finite | (is_odd ^ is_neg); 125 offset = d->invpi / (keep_sign ? r_x : -r_x); 126 scale = r_x; 127 128 result = pw_horner_8_f64 (r_x2, r_x4, d->cot_poly); 129 } 130 else 131 { 132 double r_x2 = x_reduced * x_reduced; 133 double r_x4 = r_x2 * r_x2; 134 135 offset = d->pi * x_reduced; 136 scale = x_reduced * r_x2; 137 138 result = pw_horner_13_f64 (r_x2, r_x4, d->tan_poly); 139 } 140 141 return fma (scale, result, offset); 142 } 143 144 #if WANT_EXPERIMENTAL_MATH 145 double 146 tanpi (double x) 147 { 148 return arm_math_tanpi (x); 149 } 150 #endif 151 152 #if WANT_TRIGPI_TESTS 153 TEST_ULP (arm_math_tanpi, 1.69) 154 TEST_SYM_INTERVAL (arm_math_tanpi, 0, 0x1p-63, 50000) 155 TEST_SYM_INTERVAL (arm_math_tanpi, 0x1p-63, 0.5, 100000) 156 TEST_SYM_INTERVAL (arm_math_tanpi, 0.5, 0x1p53, 100000) 157 TEST_SYM_INTERVAL (arm_math_tanpi, 0x1p53, inf, 100000) 158 #endif 159