131914882SAlex Richardson /*
231914882SAlex Richardson * Single-precision sin/cos function.
331914882SAlex Richardson *
4*f3087befSAndrew Turner * Copyright (c) 2018-2024, Arm Limited.
5072a4ba8SAndrew 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"
12*f3087befSAndrew Turner #include "test_defs.h"
1331914882SAlex Richardson
1431914882SAlex Richardson /* Fast sincosf implementation. Worst-case ULP is 0.5607, maximum relative
1531914882SAlex Richardson error is 0.5303 * 2^-23. A single-step range reduction is used for
1631914882SAlex Richardson small values. Large inputs have their range reduced using fast integer
1731914882SAlex Richardson arithmetic. */
1831914882SAlex Richardson void
sincosf(float y,float * sinp,float * cosp)1931914882SAlex Richardson sincosf (float y, float *sinp, float *cosp)
2031914882SAlex Richardson {
2131914882SAlex Richardson double x = y;
2231914882SAlex Richardson double s;
2331914882SAlex Richardson int n;
2431914882SAlex Richardson const sincos_t *p = &__sincosf_table[0];
2531914882SAlex Richardson
26d49ad206SAndrew Turner if (abstop12 (y) < abstop12 (pio4f))
2731914882SAlex Richardson {
2831914882SAlex Richardson double x2 = x * x;
2931914882SAlex Richardson
3031914882SAlex Richardson if (unlikely (abstop12 (y) < abstop12 (0x1p-12f)))
3131914882SAlex Richardson {
3231914882SAlex Richardson if (unlikely (abstop12 (y) < abstop12 (0x1p-126f)))
3331914882SAlex Richardson /* Force underflow for tiny y. */
3431914882SAlex Richardson force_eval_float (x2);
3531914882SAlex Richardson *sinp = y;
3631914882SAlex Richardson *cosp = 1.0f;
3731914882SAlex Richardson return;
3831914882SAlex Richardson }
3931914882SAlex Richardson
4031914882SAlex Richardson sincosf_poly (x, x2, p, 0, sinp, cosp);
4131914882SAlex Richardson }
4231914882SAlex Richardson else if (abstop12 (y) < abstop12 (120.0f))
4331914882SAlex Richardson {
4431914882SAlex Richardson x = reduce_fast (x, p, &n);
4531914882SAlex Richardson
4631914882SAlex Richardson /* Setup the signs for sin and cos. */
4731914882SAlex Richardson s = p->sign[n & 3];
4831914882SAlex Richardson
4931914882SAlex Richardson if (n & 2)
5031914882SAlex Richardson p = &__sincosf_table[1];
5131914882SAlex Richardson
5231914882SAlex Richardson sincosf_poly (x * s, x * x, p, n, sinp, cosp);
5331914882SAlex Richardson }
5431914882SAlex Richardson else if (likely (abstop12 (y) < abstop12 (INFINITY)))
5531914882SAlex Richardson {
5631914882SAlex Richardson uint32_t xi = asuint (y);
5731914882SAlex Richardson int sign = xi >> 31;
5831914882SAlex Richardson
5931914882SAlex Richardson x = reduce_large (xi, &n);
6031914882SAlex Richardson
6131914882SAlex Richardson /* Setup signs for sin and cos - include original sign. */
6231914882SAlex Richardson s = p->sign[(n + sign) & 3];
6331914882SAlex Richardson
6431914882SAlex Richardson if ((n + sign) & 2)
6531914882SAlex Richardson p = &__sincosf_table[1];
6631914882SAlex Richardson
6731914882SAlex Richardson sincosf_poly (x * s, x * x, p, n, sinp, cosp);
6831914882SAlex Richardson }
6931914882SAlex Richardson else
7031914882SAlex Richardson {
7131914882SAlex Richardson /* Return NaN if Inf or NaN for both sin and cos. */
7231914882SAlex Richardson *sinp = *cosp = y - y;
7331914882SAlex Richardson #if WANT_ERRNO
7431914882SAlex Richardson /* Needed to set errno for +-Inf, the add is a hack to work
7531914882SAlex Richardson around a gcc register allocation issue: just passing y
7631914882SAlex Richardson affects code generation in the fast path. */
7731914882SAlex Richardson __math_invalidf (y + y);
7831914882SAlex Richardson #endif
7931914882SAlex Richardson }
8031914882SAlex Richardson }
81*f3087befSAndrew Turner
82*f3087befSAndrew Turner TEST_ULP (sincosf_sinf, 0.06)
83*f3087befSAndrew Turner TEST_ULP (sincosf_cosf, 0.06)
84*f3087befSAndrew Turner TEST_ULP_NONNEAREST (sincosf_sinf, 0.5)
85*f3087befSAndrew Turner TEST_ULP_NONNEAREST (sincosf_cosf, 0.5)
86*f3087befSAndrew Turner TEST_INTERVAL (sincosf_sinf, 0, 0xffff0000, 10000)
87*f3087befSAndrew Turner TEST_SYM_INTERVAL (sincosf_sinf, 0x1p-14, 0x1p54, 50000)
88*f3087befSAndrew Turner TEST_INTERVAL (sincosf_cosf, 0, 0xffff0000, 10000)
89*f3087befSAndrew Turner TEST_SYM_INTERVAL (sincosf_cosf, 0x1p-14, 0x1p54, 50000)
90