1*412f47f9SXin Li /* 2*412f47f9SXin Li * Extended precision scalar reference functions for trigpi. 3*412f47f9SXin Li * 4*412f47f9SXin Li * Copyright (c) 2023, Arm Limited. 5*412f47f9SXin Li * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6*412f47f9SXin Li */ 7*412f47f9SXin Li 8*412f47f9SXin Li #define _GNU_SOURCE 9*412f47f9SXin Li #include "math_config.h" 10*412f47f9SXin Li #include "mathlib.h" 11*412f47f9SXin Li 12*412f47f9SXin Li long double sinpil(long double x)13*412f47f9SXin Lisinpil (long double x) 14*412f47f9SXin Li { 15*412f47f9SXin Li /* sin(inf) should return nan, as defined by C23. */ 16*412f47f9SXin Li if (isinf (x)) 17*412f47f9SXin Li return __math_invalid (x); 18*412f47f9SXin Li 19*412f47f9SXin Li long double ax = fabsl (x); 20*412f47f9SXin Li 21*412f47f9SXin Li /* Return 0 for all values above 2^64 to prevent 22*412f47f9SXin Li overflow when casting to uint64_t. */ 23*412f47f9SXin Li if (ax >= 0x1p64) 24*412f47f9SXin Li return 0; 25*412f47f9SXin Li 26*412f47f9SXin Li /* All integer cases should return 0. */ 27*412f47f9SXin Li if (ax == (uint64_t) ax) 28*412f47f9SXin Li return 0; 29*412f47f9SXin Li 30*412f47f9SXin Li return sinl (x * M_PIl); 31*412f47f9SXin Li } 32*412f47f9SXin Li 33*412f47f9SXin Li long double cospil(long double x)34*412f47f9SXin Licospil (long double x) 35*412f47f9SXin Li { 36*412f47f9SXin Li /* cos(inf) should return nan, as defined by C23. */ 37*412f47f9SXin Li if (isinf (x)) 38*412f47f9SXin Li return __math_invalid (x); 39*412f47f9SXin Li 40*412f47f9SXin Li long double ax = fabsl (x); 41*412f47f9SXin Li 42*412f47f9SXin Li if (ax >= 0x1p64) 43*412f47f9SXin Li return 1; 44*412f47f9SXin Li 45*412f47f9SXin Li uint64_t m = (uint64_t) ax; 46*412f47f9SXin Li 47*412f47f9SXin Li /* Integer values of cospi(x) should return +/-1. 48*412f47f9SXin Li The sign depends on if x is odd or even. */ 49*412f47f9SXin Li if (m == ax) 50*412f47f9SXin Li return (m & 1) ? -1 : 1; 51*412f47f9SXin Li 52*412f47f9SXin Li /* Values of Integer + 0.5 should always return 0. */ 53*412f47f9SXin Li if (ax - 0.5 == m || ax + 0.5 == m) 54*412f47f9SXin Li return 0; 55*412f47f9SXin Li 56*412f47f9SXin Li return cosl (ax * M_PIl); 57*412f47f9SXin Li }