1//===-- Compute sin + cos for small angles ----------------------*- C++ -*-===//
2//
3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4// See https://llvm.org/LICENSE.txt for license information.
5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6//
7//===----------------------------------------------------------------------===//
8
9#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_SINCOS_EVAL_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_SINCOS_EVAL_H
11
12#include "src/__support/FPUtil/PolyEval.h"
13#include "src/__support/FPUtil/double_double.h"
14#include "src/__support/FPUtil/dyadic_float.h"
15#include "src/__support/FPUtil/multiply_add.h"
16#include "src/__support/integer_literals.h"
17#include "src/__support/macros/config.h"
18
19namespace LIBC_NAMESPACE_DECL {
20
21namespace math {
22
23namespace sincos_eval_internal {
24
25using fputil::DoubleDouble;
26using DFloat128 = fputil::DyadicFloat<128>;
27
28LIBC_INLINE double sincos_eval(const DoubleDouble &u, DoubleDouble &sin_u,
29 DoubleDouble &cos_u) {
30 // Evaluate sin(y) = sin(x - k * (pi/128))
31 // We use the degree-7 Taylor approximation:
32 // sin(y) ~ y - y^3/3! + y^5/5! - y^7/7!
33 // Then the error is bounded by:
34 // |sin(y) - (y - y^3/3! + y^5/5! - y^7/7!)| < |y|^9/9! < 2^-54/9! < 2^-72.
35 // For y ~ u_hi + u_lo, fully expanding the polynomial and drop any terms
36 // < ulp(u_hi^3) gives us:
37 // y - y^3/3! + y^5/5! - y^7/7! = ...
38 // ~ u_hi + u_hi^3 * (-1/6 + u_hi^2 * (1/120 - u_hi^2 * 1/5040)) +
39 // + u_lo (1 + u_hi^2 * (-1/2 + u_hi^2 / 24))
40 double u_hi_sq = u.hi * u.hi; // Error < ulp(u_hi^2) < 2^(-6 - 52) = 2^-58.
41 // p1 ~ 1/120 + u_hi^2 / 5040.
42 double p1 = fputil::multiply_add(x: u_hi_sq, y: -0x1.a01a01a01a01ap-13,
43 z: 0x1.1111111111111p-7);
44 // q1 ~ -1/2 + u_hi^2 / 24.
45 double q1 = fputil::multiply_add(x: u_hi_sq, y: 0x1.5555555555555p-5, z: -0x1.0p-1);
46 double u_hi_3 = u_hi_sq * u.hi;
47 // p2 ~ -1/6 + u_hi^2 (1/120 - u_hi^2 * 1/5040)
48 double p2 = fputil::multiply_add(x: u_hi_sq, y: p1, z: -0x1.5555555555555p-3);
49 // q2 ~ 1 + u_hi^2 (-1/2 + u_hi^2 / 24)
50 double q2 = fputil::multiply_add(x: u_hi_sq, y: q1, z: 1.0);
51 double sin_lo = fputil::multiply_add(x: u_hi_3, y: p2, z: u.lo * q2);
52 // Overall, |sin(y) - (u_hi + sin_lo)| < 2*ulp(u_hi^3) < 2^-69.
53
54 // Evaluate cos(y) = cos(x - k * (pi/128))
55 // We use the degree-8 Taylor approximation:
56 // cos(y) ~ 1 - y^2/2 + y^4/4! - y^6/6! + y^8/8!
57 // Then the error is bounded by:
58 // |cos(y) - (...)| < |y|^10/10! < 2^-81
59 // For y ~ u_hi + u_lo, fully expanding the polynomial and drop any terms
60 // < ulp(u_hi^3) gives us:
61 // 1 - y^2/2 + y^4/4! - y^6/6! + y^8/8! = ...
62 // ~ 1 - u_hi^2/2 + u_hi^4(1/24 + u_hi^2 (-1/720 + u_hi^2/40320)) +
63 // + u_hi u_lo (-1 + u_hi^2/6)
64 // We compute 1 - u_hi^2 accurately:
65 // v_hi + v_lo ~ 1 - u_hi^2/2
66 // with error <= 2^-105.
67 double u_hi_neg_half = (-0.5) * u.hi;
68 DoubleDouble v;
69
70#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
71 v.hi = fputil::multiply_add(u.hi, u_hi_neg_half, 1.0);
72 v.lo = 1.0 - v.hi; // Exact
73 v.lo = fputil::multiply_add(u.hi, u_hi_neg_half, v.lo);
74#else
75 DoubleDouble u_hi_sq_neg_half = fputil::exact_mult(a: u.hi, b: u_hi_neg_half);
76 v = fputil::exact_add(a: 1.0, b: u_hi_sq_neg_half.hi);
77 v.lo += u_hi_sq_neg_half.lo;
78#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
79
80 // r1 ~ -1/720 + u_hi^2 / 40320
81 double r1 = fputil::multiply_add(x: u_hi_sq, y: 0x1.a01a01a01a01ap-16,
82 z: -0x1.6c16c16c16c17p-10);
83 // s1 ~ -1 + u_hi^2 / 6
84 double s1 = fputil::multiply_add(x: u_hi_sq, y: 0x1.5555555555555p-3, z: -1.0);
85 double u_hi_4 = u_hi_sq * u_hi_sq;
86 double u_hi_u_lo = u.hi * u.lo;
87 // r2 ~ 1/24 + u_hi^2 (-1/720 + u_hi^2 / 40320)
88 double r2 = fputil::multiply_add(x: u_hi_sq, y: r1, z: 0x1.5555555555555p-5);
89 // s2 ~ v_lo + u_hi * u_lo * (-1 + u_hi^2 / 6)
90 double s2 = fputil::multiply_add(x: u_hi_u_lo, y: s1, z: v.lo);
91 double cos_lo = fputil::multiply_add(x: u_hi_4, y: r2, z: s2);
92 // Overall, |cos(y) - (v_hi + cos_lo)| < 2*ulp(u_hi^4) < 2^-75.
93
94 sin_u = fputil::exact_add(a: u.hi, b: sin_lo);
95 cos_u = fputil::exact_add(a: v.hi, b: cos_lo);
96
97 return fputil::multiply_add(x: fputil::FPBits<double>(u_hi_3).abs().get_val(),
98 y: 0x1.0p-51, z: 0x1.0p-105);
99}
100
101LIBC_INLINE void sincos_eval(const DFloat128 &u, DFloat128 &sin_u,
102 DFloat128 &cos_u) {
103 DFloat128 u_sq = fputil::quick_mul(a: u, b: u);
104
105 // sin(u) ~ x - x^3/3! + x^5/5! - x^7/7! + x^9/9! - x^11/11! + x^13/13!
106 constexpr DFloat128 SIN_COEFFS[] = {
107 {Sign::POS, -127, 0x80000000'00000000'00000000'00000000_u128}, // 1
108 {Sign::NEG, -130, 0xaaaaaaaa'aaaaaaaa'aaaaaaaa'aaaaaaab_u128}, // -1/3!
109 {Sign::POS, -134, 0x88888888'88888888'88888888'88888889_u128}, // 1/5!
110 {Sign::NEG, -140, 0xd00d00d0'0d00d00d'00d00d00'd00d00d0_u128}, // -1/7!
111 {Sign::POS, -146, 0xb8ef1d2a'b6399c7d'560e4472'800b8ef2_u128}, // 1/9!
112 {Sign::NEG, -153, 0xd7322b3f'aa271c7f'3a3f25c1'bee38f10_u128}, // -1/11!
113 {Sign::POS, -160, 0xb092309d'43684be5'1c198e91'd7b4269e_u128}, // 1/13!
114 };
115
116 // cos(u) ~ 1 - x^2/2 + x^4/4! - x^6/6! + x^8/8! - x^10/10! + x^12/12!
117 constexpr DFloat128 COS_COEFFS[] = {
118 {Sign::POS, -127, 0x80000000'00000000'00000000'00000000_u128}, // 1.0
119 {Sign::NEG, -128, 0x80000000'00000000'00000000'00000000_u128}, // 1/2
120 {Sign::POS, -132, 0xaaaaaaaa'aaaaaaaa'aaaaaaaa'aaaaaaab_u128}, // 1/4!
121 {Sign::NEG, -137, 0xb60b60b6'0b60b60b'60b60b60'b60b60b6_u128}, // 1/6!
122 {Sign::POS, -143, 0xd00d00d0'0d00d00d'00d00d00'd00d00d0_u128}, // 1/8!
123 {Sign::NEG, -149, 0x93f27dbb'c4fae397'780b69f5'333c725b_u128}, // 1/10!
124 {Sign::POS, -156, 0x8f76c77f'c6c4bdaa'26d4c3d6'7f425f60_u128}, // 1/12!
125 };
126
127 sin_u = fputil::quick_mul(a: u, b: fputil::polyeval(x: u_sq, a0: SIN_COEFFS[0],
128 a: SIN_COEFFS[1], a: SIN_COEFFS[2],
129 a: SIN_COEFFS[3], a: SIN_COEFFS[4],
130 a: SIN_COEFFS[5], a: SIN_COEFFS[6]));
131 cos_u = fputil::polyeval(x: u_sq, a0: COS_COEFFS[0], a: COS_COEFFS[1], a: COS_COEFFS[2],
132 a: COS_COEFFS[3], a: COS_COEFFS[4], a: COS_COEFFS[5],
133 a: COS_COEFFS[6]);
134}
135
136} // namespace sincos_eval_internal
137
138} // namespace math
139
140} // namespace LIBC_NAMESPACE_DECL
141
142#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINCOS_EVAL_H
143