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