1*31914882SAlex Richardson /* 2*31914882SAlex Richardson * Single-precision sin/cos function. 3*31914882SAlex Richardson * 4*31914882SAlex Richardson * Copyright (c) 2018-2019, Arm Limited. 5*31914882SAlex Richardson * SPDX-License-Identifier: MIT 6*31914882SAlex Richardson */ 7*31914882SAlex Richardson 8*31914882SAlex Richardson #include <stdint.h> 9*31914882SAlex Richardson #include <math.h> 10*31914882SAlex Richardson #include "math_config.h" 11*31914882SAlex Richardson #include "sincosf.h" 12*31914882SAlex Richardson 13*31914882SAlex Richardson /* Fast sincosf implementation. Worst-case ULP is 0.5607, maximum relative 14*31914882SAlex Richardson error is 0.5303 * 2^-23. A single-step range reduction is used for 15*31914882SAlex Richardson small values. Large inputs have their range reduced using fast integer 16*31914882SAlex Richardson arithmetic. */ 17*31914882SAlex Richardson void 18*31914882SAlex Richardson sincosf (float y, float *sinp, float *cosp) 19*31914882SAlex Richardson { 20*31914882SAlex Richardson double x = y; 21*31914882SAlex Richardson double s; 22*31914882SAlex Richardson int n; 23*31914882SAlex Richardson const sincos_t *p = &__sincosf_table[0]; 24*31914882SAlex Richardson 25*31914882SAlex Richardson if (abstop12 (y) < abstop12 (pio4)) 26*31914882SAlex Richardson { 27*31914882SAlex Richardson double x2 = x * x; 28*31914882SAlex Richardson 29*31914882SAlex Richardson if (unlikely (abstop12 (y) < abstop12 (0x1p-12f))) 30*31914882SAlex Richardson { 31*31914882SAlex Richardson if (unlikely (abstop12 (y) < abstop12 (0x1p-126f))) 32*31914882SAlex Richardson /* Force underflow for tiny y. */ 33*31914882SAlex Richardson force_eval_float (x2); 34*31914882SAlex Richardson *sinp = y; 35*31914882SAlex Richardson *cosp = 1.0f; 36*31914882SAlex Richardson return; 37*31914882SAlex Richardson } 38*31914882SAlex Richardson 39*31914882SAlex Richardson sincosf_poly (x, x2, p, 0, sinp, cosp); 40*31914882SAlex Richardson } 41*31914882SAlex Richardson else if (abstop12 (y) < abstop12 (120.0f)) 42*31914882SAlex Richardson { 43*31914882SAlex Richardson x = reduce_fast (x, p, &n); 44*31914882SAlex Richardson 45*31914882SAlex Richardson /* Setup the signs for sin and cos. */ 46*31914882SAlex Richardson s = p->sign[n & 3]; 47*31914882SAlex Richardson 48*31914882SAlex Richardson if (n & 2) 49*31914882SAlex Richardson p = &__sincosf_table[1]; 50*31914882SAlex Richardson 51*31914882SAlex Richardson sincosf_poly (x * s, x * x, p, n, sinp, cosp); 52*31914882SAlex Richardson } 53*31914882SAlex Richardson else if (likely (abstop12 (y) < abstop12 (INFINITY))) 54*31914882SAlex Richardson { 55*31914882SAlex Richardson uint32_t xi = asuint (y); 56*31914882SAlex Richardson int sign = xi >> 31; 57*31914882SAlex Richardson 58*31914882SAlex Richardson x = reduce_large (xi, &n); 59*31914882SAlex Richardson 60*31914882SAlex Richardson /* Setup signs for sin and cos - include original sign. */ 61*31914882SAlex Richardson s = p->sign[(n + sign) & 3]; 62*31914882SAlex Richardson 63*31914882SAlex Richardson if ((n + sign) & 2) 64*31914882SAlex Richardson p = &__sincosf_table[1]; 65*31914882SAlex Richardson 66*31914882SAlex Richardson sincosf_poly (x * s, x * x, p, n, sinp, cosp); 67*31914882SAlex Richardson } 68*31914882SAlex Richardson else 69*31914882SAlex Richardson { 70*31914882SAlex Richardson /* Return NaN if Inf or NaN for both sin and cos. */ 71*31914882SAlex Richardson *sinp = *cosp = y - y; 72*31914882SAlex Richardson #if WANT_ERRNO 73*31914882SAlex Richardson /* Needed to set errno for +-Inf, the add is a hack to work 74*31914882SAlex Richardson around a gcc register allocation issue: just passing y 75*31914882SAlex Richardson affects code generation in the fast path. */ 76*31914882SAlex Richardson __math_invalidf (y + y); 77*31914882SAlex Richardson #endif 78*31914882SAlex Richardson } 79*31914882SAlex Richardson } 80