1bbb75554SSiva Chandra //===-- Single-precision sincos function ----------------------------------===// 2bbb75554SSiva Chandra // 3bbb75554SSiva Chandra // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. 4bbb75554SSiva Chandra // See https://llvm.org/LICENSE.txt for license information. 5bbb75554SSiva Chandra // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception 6bbb75554SSiva Chandra // 7bbb75554SSiva Chandra //===----------------------------------------------------------------------===// 8bbb75554SSiva Chandra 9bbb75554SSiva Chandra #include "src/math/sincosf.h" 10bbb75554SSiva Chandra #include "math_utils.h" 11bbb75554SSiva Chandra #include "sincosf_utils.h" 12bbb75554SSiva Chandra 13bbb75554SSiva Chandra #include "src/__support/common.h" 14bbb75554SSiva Chandra #include <math.h> 15bbb75554SSiva Chandra 16bbb75554SSiva Chandra #include <stdint.h> 17bbb75554SSiva Chandra 18bbb75554SSiva Chandra namespace __llvm_libc { 19bbb75554SSiva Chandra 20bbb75554SSiva Chandra // Fast sincosf implementation. Worst-case ULP is 0.5607, maximum relative 21bbb75554SSiva Chandra // error is 0.5303 * 2^-23. A single-step range reduction is used for 22bbb75554SSiva Chandra // small values. Large inputs have their range reduced using fast integer 23bbb75554SSiva Chandra // arithmetic. 24bbb75554SSiva Chandra LLVM_LIBC_FUNCTION(void, sincosf, (float y, float *sinp, float *cosp)) { 25bbb75554SSiva Chandra double x = y; 26bbb75554SSiva Chandra double s; 27bbb75554SSiva Chandra int n; 28*75d2fcb0SSiva Chandra Reddy const sincos_t *p = &SINCOSF_TABLE[0]; 29bbb75554SSiva Chandra 30*75d2fcb0SSiva Chandra Reddy if (abstop12(y) < abstop12(PIO4)) { 31bbb75554SSiva Chandra double x2 = x * x; 32bbb75554SSiva Chandra 33bbb75554SSiva Chandra if (unlikely(abstop12(y) < abstop12(as_float(0x39800000)))) { 34bbb75554SSiva Chandra if (unlikely(abstop12(y) < abstop12(as_float(0x800000)))) 35bbb75554SSiva Chandra // Force underflow for tiny y. 36bbb75554SSiva Chandra force_eval<float>(x2); 37bbb75554SSiva Chandra *sinp = y; 38bbb75554SSiva Chandra *cosp = 1.0f; 39bbb75554SSiva Chandra return; 40bbb75554SSiva Chandra } 41bbb75554SSiva Chandra 42bbb75554SSiva Chandra sincosf_poly(x, x2, p, 0, sinp, cosp); 43bbb75554SSiva Chandra } else if (abstop12(y) < abstop12(120.0f)) { 44bbb75554SSiva Chandra x = reduce_fast(x, p, &n); 45bbb75554SSiva Chandra 46bbb75554SSiva Chandra // Setup the signs for sin and cos. 47bbb75554SSiva Chandra s = p->sign[n & 3]; 48bbb75554SSiva Chandra 49bbb75554SSiva Chandra if (n & 2) 50*75d2fcb0SSiva Chandra Reddy p = &SINCOSF_TABLE[1]; 51bbb75554SSiva Chandra 52bbb75554SSiva Chandra sincosf_poly(x * s, x * x, p, n, sinp, cosp); 53bbb75554SSiva Chandra } else if (likely(abstop12(y) < abstop12(INFINITY))) { 54bbb75554SSiva Chandra uint32_t xi = as_uint32_bits(y); 55bbb75554SSiva Chandra int sign = xi >> 31; 56bbb75554SSiva Chandra 57bbb75554SSiva Chandra x = reduce_large(xi, &n); 58bbb75554SSiva Chandra 59bbb75554SSiva Chandra // Setup signs for sin and cos - include original sign. 60bbb75554SSiva Chandra s = p->sign[(n + sign) & 3]; 61bbb75554SSiva Chandra 62bbb75554SSiva Chandra if ((n + sign) & 2) 63*75d2fcb0SSiva Chandra Reddy p = &SINCOSF_TABLE[1]; 64bbb75554SSiva Chandra 65bbb75554SSiva Chandra sincosf_poly(x * s, x * x, p, n, sinp, cosp); 66bbb75554SSiva Chandra } else { 67bbb75554SSiva Chandra // Return NaN if Inf or NaN for both sin and cos. 68bbb75554SSiva Chandra *sinp = *cosp = y - y; 69bbb75554SSiva Chandra 70bbb75554SSiva Chandra // Needed to set errno for +-Inf, the add is a hack to work 71bbb75554SSiva Chandra // around a gcc register allocation issue: just passing y 72bbb75554SSiva Chandra // affects code generation in the fast path. 73bbb75554SSiva Chandra invalid(y + y); 74bbb75554SSiva Chandra } 75bbb75554SSiva Chandra } 76bbb75554SSiva Chandra 77bbb75554SSiva Chandra } // namespace __llvm_libc 78