1 /* 2 * Extended precision scalar reference functions for trigpi. 3 * 4 * Copyright (c) 2023-2024, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 8 #include "math_config.h" 9 10 #ifndef M_PIl 11 # define M_PIl 3.141592653589793238462643383279502884l 12 #endif 13 14 long double arm_math_sinpil(long double x)15arm_math_sinpil (long double x) 16 { 17 /* sin(inf) should return nan, as defined by C23. */ 18 if (isinf (x)) 19 return __math_invalid (x); 20 21 long double ax = fabsl (x); 22 23 /* Return 0 for all values above 2^64 to prevent 24 overflow when casting to uint64_t. */ 25 if (ax >= 0x1p64) 26 return x < 0 ? -0.0l : 0.0l; 27 28 /* All integer cases should return 0, with unchanged sign for zero. */ 29 if (x == 0.0l) 30 return x; 31 if (ax == (uint64_t) ax) 32 return x < 0 ? -0.0l : 0.0l; 33 34 return sinl (x * M_PIl); 35 } 36 37 long double arm_math_cospil(long double x)38arm_math_cospil (long double x) 39 { 40 /* cos(inf) should return nan, as defined by C23. */ 41 if (isinf (x)) 42 return __math_invalid (x); 43 44 long double ax = fabsl (x); 45 46 if (ax >= 0x1p64) 47 return 1; 48 49 uint64_t m = (uint64_t) ax; 50 51 /* Integer values of cospi(x) should return +/-1. 52 The sign depends on if x is odd or even. */ 53 if (m == ax) 54 return (m & 1) ? -1 : 1; 55 56 /* Values of Integer + 0.5 should always return 0. */ 57 if (ax - 0.5 == m || ax + 0.5 == m) 58 return 0; 59 60 return cosl (ax * M_PIl); 61 } 62 63 long double arm_math_tanpil(long double x)64arm_math_tanpil (long double x) 65 { 66 /* inf and x = n + 0.5 for any integral n should return nan. */ 67 if (fabsl (x) >= 0x1p54l) 68 { 69 if (isinf (x)) 70 return __math_invalid (x); 71 return x < 0 ? -0.0l : 0.0l; 72 } 73 74 long double i = roundl (x); 75 long double f = x - i; 76 int64_t m = (int64_t) i; 77 78 if (x == 0) 79 { 80 return x; 81 } 82 else if (x == i) 83 { 84 if (x < 0) 85 { 86 return m & 1 ? 0.0l : -0.0l; 87 } 88 else 89 { 90 return m & 1 ? -0.0l : 0.0l; 91 } 92 } 93 else if (fabsl (f) == 0.5l) 94 { 95 if (x < 0) 96 { 97 return m & 1 ? -1.0l / 0.0l : 1.0l / 0.0l; 98 } 99 else 100 { 101 return m & 1 ? 1.0l / 0.0l : -1.0l / 0.0l; 102 } 103 } 104 105 return tanl (f * M_PIl); 106 } 107