1/*
2 * Double-precision scalar sinpi function.
3 *
4 * Copyright (c) 2023, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8#define _GNU_SOURCE
9#include <math.h>
10#include "mathlib.h"
11#include "math_config.h"
12#include "pl_sig.h"
13#include "pl_test.h"
14#include "poly_scalar_f64.h"
15
16/* Taylor series coefficents for sin(pi * x).
17   C2 coefficient (orginally ~=5.16771278) has been split into two parts:
18   C2_hi = 4, C2_lo = C2 - C2_hi (~=1.16771278)
19   This change in magnitude reduces floating point rounding errors.
20   C2_hi is then reintroduced after the polynomial approxmation.  */
21static const double poly[]
22    = { 0x1.921fb54442d184p1,  -0x1.2aef39896f94bp0,   0x1.466bc6775ab16p1,
23	-0x1.32d2cce62dc33p-1, 0x1.507834891188ep-4,   -0x1.e30750a28c88ep-8,
24	0x1.e8f48308acda4p-12, -0x1.6fc0032b3c29fp-16, 0x1.af86ae521260bp-21,
25	-0x1.012a9870eeb7dp-25 };
26
27#define Shift 0x1.8p+52
28
29/* Approximation for scalar double-precision sinpi(x).
30   Maximum error: 3.03 ULP:
31   sinpi(0x1.a90da2818f8b5p+7) got 0x1.fe358f255a4b3p-1
32			      want 0x1.fe358f255a4b6p-1.  */
33double
34sinpi (double x)
35{
36  if (isinf (x))
37    return __math_invalid (x);
38
39  double r = asdouble (asuint64 (x) & ~0x8000000000000000);
40  uint64_t sign = asuint64 (x) & 0x8000000000000000;
41
42  /* Edge cases for when sinpif should be exactly 0. (Integers)
43     0x1p53 is the limit for single precision to store any decimal places.  */
44  if (r >= 0x1p53)
45    return 0;
46
47  /* If x is an integer, return 0.  */
48  uint64_t m = (uint64_t) r;
49  if (r == m)
50    return 0;
51
52  /* For very small inputs, squaring r causes underflow.
53     Values below this threshold can be approximated via sinpi(x) ��� pi*x.  */
54  if (r < 0x1p-63)
55    return M_PI * x;
56
57  /* Any non-integer values >= 0x1x51 will be int + 0.5.
58     These values should return exactly 1 or -1.  */
59  if (r >= 0x1p51)
60    {
61      uint64_t iy = ((m & 1) << 63) ^ asuint64 (1.0);
62      return asdouble (sign ^ iy);
63    }
64
65  /* n = rint(|x|).  */
66  double n = r + Shift;
67  sign ^= (asuint64 (n) << 63);
68  n = n - Shift;
69
70  /* r = |x| - n (range reduction into -1/2 .. 1/2).  */
71  r = r - n;
72
73  /* y = sin(r).  */
74  double r2 = r * r;
75  double y = horner_9_f64 (r2, poly);
76  y = y * r;
77
78  /* Reintroduce C2_hi.  */
79  y = fma (-4 * r2, r, y);
80
81  /* Copy sign of x to sin(|x|).  */
82  return asdouble (asuint64 (y) ^ sign);
83}
84
85PL_SIG (S, D, 1, sinpi, -0.9, 0.9)
86PL_TEST_ULP (sinpi, 2.53)
87PL_TEST_SYM_INTERVAL (sinpi, 0, 0x1p-63, 5000)
88PL_TEST_SYM_INTERVAL (sinpi, 0x1p-63, 0.5, 10000)
89PL_TEST_SYM_INTERVAL (sinpi, 0.5, 0x1p51, 10000)
90PL_TEST_SYM_INTERVAL (sinpi, 0x1p51, inf, 10000)
91