131914882SAlex Richardson /* 231914882SAlex Richardson * Single-precision sin/cos function. 331914882SAlex Richardson * 4*072a4ba8SAndrew Turner * Copyright (c) 2018-2021, Arm Limited. 5*072a4ba8SAndrew Turner * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 631914882SAlex Richardson */ 731914882SAlex Richardson 831914882SAlex Richardson #include <stdint.h> 931914882SAlex Richardson #include <math.h> 1031914882SAlex Richardson #include "math_config.h" 1131914882SAlex Richardson #include "sincosf.h" 1231914882SAlex Richardson 1331914882SAlex Richardson /* Fast sincosf implementation. Worst-case ULP is 0.5607, maximum relative 1431914882SAlex Richardson error is 0.5303 * 2^-23. A single-step range reduction is used for 1531914882SAlex Richardson small values. Large inputs have their range reduced using fast integer 1631914882SAlex Richardson arithmetic. */ 1731914882SAlex Richardson void 1831914882SAlex Richardson sincosf (float y, float *sinp, float *cosp) 1931914882SAlex Richardson { 2031914882SAlex Richardson double x = y; 2131914882SAlex Richardson double s; 2231914882SAlex Richardson int n; 2331914882SAlex Richardson const sincos_t *p = &__sincosf_table[0]; 2431914882SAlex Richardson 25d49ad206SAndrew Turner if (abstop12 (y) < abstop12 (pio4f)) 2631914882SAlex Richardson { 2731914882SAlex Richardson double x2 = x * x; 2831914882SAlex Richardson 2931914882SAlex Richardson if (unlikely (abstop12 (y) < abstop12 (0x1p-12f))) 3031914882SAlex Richardson { 3131914882SAlex Richardson if (unlikely (abstop12 (y) < abstop12 (0x1p-126f))) 3231914882SAlex Richardson /* Force underflow for tiny y. */ 3331914882SAlex Richardson force_eval_float (x2); 3431914882SAlex Richardson *sinp = y; 3531914882SAlex Richardson *cosp = 1.0f; 3631914882SAlex Richardson return; 3731914882SAlex Richardson } 3831914882SAlex Richardson 3931914882SAlex Richardson sincosf_poly (x, x2, p, 0, sinp, cosp); 4031914882SAlex Richardson } 4131914882SAlex Richardson else if (abstop12 (y) < abstop12 (120.0f)) 4231914882SAlex Richardson { 4331914882SAlex Richardson x = reduce_fast (x, p, &n); 4431914882SAlex Richardson 4531914882SAlex Richardson /* Setup the signs for sin and cos. */ 4631914882SAlex Richardson s = p->sign[n & 3]; 4731914882SAlex Richardson 4831914882SAlex Richardson if (n & 2) 4931914882SAlex Richardson p = &__sincosf_table[1]; 5031914882SAlex Richardson 5131914882SAlex Richardson sincosf_poly (x * s, x * x, p, n, sinp, cosp); 5231914882SAlex Richardson } 5331914882SAlex Richardson else if (likely (abstop12 (y) < abstop12 (INFINITY))) 5431914882SAlex Richardson { 5531914882SAlex Richardson uint32_t xi = asuint (y); 5631914882SAlex Richardson int sign = xi >> 31; 5731914882SAlex Richardson 5831914882SAlex Richardson x = reduce_large (xi, &n); 5931914882SAlex Richardson 6031914882SAlex Richardson /* Setup signs for sin and cos - include original sign. */ 6131914882SAlex Richardson s = p->sign[(n + sign) & 3]; 6231914882SAlex Richardson 6331914882SAlex Richardson if ((n + sign) & 2) 6431914882SAlex Richardson p = &__sincosf_table[1]; 6531914882SAlex Richardson 6631914882SAlex Richardson sincosf_poly (x * s, x * x, p, n, sinp, cosp); 6731914882SAlex Richardson } 6831914882SAlex Richardson else 6931914882SAlex Richardson { 7031914882SAlex Richardson /* Return NaN if Inf or NaN for both sin and cos. */ 7131914882SAlex Richardson *sinp = *cosp = y - y; 7231914882SAlex Richardson #if WANT_ERRNO 7331914882SAlex Richardson /* Needed to set errno for +-Inf, the add is a hack to work 7431914882SAlex Richardson around a gcc register allocation issue: just passing y 7531914882SAlex Richardson affects code generation in the fast path. */ 7631914882SAlex Richardson __math_invalidf (y + y); 7731914882SAlex Richardson #endif 7831914882SAlex Richardson } 7931914882SAlex Richardson } 80