1*31914882SAlex Richardson /* 2*31914882SAlex Richardson * Single-precision sin 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 <math.h> 9*31914882SAlex Richardson #include "math_config.h" 10*31914882SAlex Richardson #include "sincosf.h" 11*31914882SAlex Richardson 12*31914882SAlex Richardson /* Fast sinf implementation. Worst-case ULP is 0.5607, maximum relative 13*31914882SAlex Richardson error is 0.5303 * 2^-23. A single-step range reduction is used for 14*31914882SAlex Richardson small values. Large inputs have their range reduced using fast integer 15*31914882SAlex Richardson arithmetic. */ 16*31914882SAlex Richardson float 17*31914882SAlex Richardson sinf (float y) 18*31914882SAlex Richardson { 19*31914882SAlex Richardson double x = y; 20*31914882SAlex Richardson double s; 21*31914882SAlex Richardson int n; 22*31914882SAlex Richardson const sincos_t *p = &__sincosf_table[0]; 23*31914882SAlex Richardson 24*31914882SAlex Richardson if (abstop12 (y) < abstop12 (pio4)) 25*31914882SAlex Richardson { 26*31914882SAlex Richardson s = x * x; 27*31914882SAlex Richardson 28*31914882SAlex Richardson if (unlikely (abstop12 (y) < abstop12 (0x1p-12f))) 29*31914882SAlex Richardson { 30*31914882SAlex Richardson if (unlikely (abstop12 (y) < abstop12 (0x1p-126f))) 31*31914882SAlex Richardson /* Force underflow for tiny y. */ 32*31914882SAlex Richardson force_eval_float (s); 33*31914882SAlex Richardson return y; 34*31914882SAlex Richardson } 35*31914882SAlex Richardson 36*31914882SAlex Richardson return sinf_poly (x, s, p, 0); 37*31914882SAlex Richardson } 38*31914882SAlex Richardson else if (likely (abstop12 (y) < abstop12 (120.0f))) 39*31914882SAlex Richardson { 40*31914882SAlex Richardson x = reduce_fast (x, p, &n); 41*31914882SAlex Richardson 42*31914882SAlex Richardson /* Setup the signs for sin and cos. */ 43*31914882SAlex Richardson s = p->sign[n & 3]; 44*31914882SAlex Richardson 45*31914882SAlex Richardson if (n & 2) 46*31914882SAlex Richardson p = &__sincosf_table[1]; 47*31914882SAlex Richardson 48*31914882SAlex Richardson return sinf_poly (x * s, x * x, p, n); 49*31914882SAlex Richardson } 50*31914882SAlex Richardson else if (abstop12 (y) < abstop12 (INFINITY)) 51*31914882SAlex Richardson { 52*31914882SAlex Richardson uint32_t xi = asuint (y); 53*31914882SAlex Richardson int sign = xi >> 31; 54*31914882SAlex Richardson 55*31914882SAlex Richardson x = reduce_large (xi, &n); 56*31914882SAlex Richardson 57*31914882SAlex Richardson /* Setup signs for sin and cos - include original sign. */ 58*31914882SAlex Richardson s = p->sign[(n + sign) & 3]; 59*31914882SAlex Richardson 60*31914882SAlex Richardson if ((n + sign) & 2) 61*31914882SAlex Richardson p = &__sincosf_table[1]; 62*31914882SAlex Richardson 63*31914882SAlex Richardson return sinf_poly (x * s, x * x, p, n); 64*31914882SAlex Richardson } 65*31914882SAlex Richardson else 66*31914882SAlex Richardson return __math_invalidf (y); 67*31914882SAlex Richardson } 68