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