|
| 1 | +//===-- double-precision cospi function ----------------------------------===// |
| 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 | +#include "src/math/cospi.h" |
| 10 | +#include "sincos_eval.h" |
| 11 | +#include "src/__support/FPUtil/BasicOperations.h" |
| 12 | +#include "src/__support/FPUtil/FEnvImpl.h" |
| 13 | +#include "src/__support/FPUtil/FPBits.h" |
| 14 | +#include "src/__support/FPUtil/PolyEval.h" |
| 15 | +#include "src/__support/FPUtil/double_double.h" |
| 16 | +#include "src/__support/FPUtil/multiply_add.h" |
| 17 | +#include "src/__support/common.h" |
| 18 | +#include "src/__support/macros/config.h" |
| 19 | +#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY |
| 20 | +#include "src/math/fmul.h" |
| 21 | +#include "src/math/generic/sincosf_utils.h" |
| 22 | +#include "src/math/pow.h" |
| 23 | + |
| 24 | +namespace LIBC_NAMESPACE_DECL { |
| 25 | + |
| 26 | +static LIBC_INLINE void sincospi_poly_eval(double k, double y, double &sin_k, |
| 27 | + double &cos_k, double &sin_y, |
| 28 | + |
| 29 | + double &cosm1_y) { |
| 30 | + |
| 31 | + // Q3 = fpminimax(sin(x*pi), 7, [|64...|], [-0.0078125, 0.0078125]); |
| 32 | + sin_k = |
| 33 | + k * fputil::polyeval(k, 0x1.59b6a771a45cbab8p-94, 0x1.921fb54442d1846ap1, |
| 34 | + -0x1.8633470ba8bd806cp-76, -0x1.4abbce625be56346p2, |
| 35 | + 0x1.d3e01dfd72e97a92p-61, 0x1.466bc67713dbbfp1, |
| 36 | + -0x1.14c2648595e2ad4p-47, -0x1.32d1cc20b89301fcp-1); |
| 37 | + |
| 38 | + sin_y = |
| 39 | + y * fputil::polyeval(k, 0x1.59b6a771a45cbab8p-94, 0x1.921fb54442d1846ap1, |
| 40 | + -0x1.8633470ba8bd806cp-76, -0x1.4abbce625be56346p2, |
| 41 | + 0x1.d3e01dfd72e97a92p-61, 0x1.466bc67713dbbfp1, |
| 42 | + -0x1.14c2648595e2ad4p-47, -0x1.32d1cc20b89301fcp-1); |
| 43 | + |
| 44 | + // Q1 = fpminimax(cos(x * pi), 7, [|64...|], [-0.0078125, 0.0078125]); |
| 45 | + cos_k = |
| 46 | + k * fputil::polyeval(k, 0x1p0, 0x1.a5b22c564ee1d862p-84, |
| 47 | + -0x1.3bd3cc9be45d30e6p2, -0x1.5c2328fefbe60d3ep-66, |
| 48 | + 0x1.03c1f080a6907a6p2, 0x1.569a4d5c5018eecap-51, |
| 49 | + -0x1.55d1f72455a9848ap0, -0x1.6b18e5f7fc6c39a6p-38); |
| 50 | + |
| 51 | + cosm1_y = |
| 52 | + y * fputil::polyeval(y, 0x1p0, 0x1.a5b22c564ee1d862p-84, |
| 53 | + -0x1.3bd3cc9be45d30e6p2, -0x1.5c2328fefbe60d3ep-66, |
| 54 | + 0x1.03c1f080a6907a6p2, 0x1.569a4d5c5018eecap-51, |
| 55 | + -0x1.55d1f72455a9848ap0, -0x1.6b18e5f7fc6c39a6p-38); |
| 56 | +} |
| 57 | + |
| 58 | +LLVM_LIBC_FUNCTION(double, cospi, (double x)) { |
| 59 | + using FPBits = typename fputil::FPBits<double>; |
| 60 | + FPBits xbits(x); |
| 61 | + |
| 62 | + double x_abs = fputil::abs(x); |
| 63 | + double p = 0x1p52; // 2^p where p is the precision |
| 64 | + |
| 65 | + if (LIBC_UNLIKELY(x_abs == 0U)) |
| 66 | + return x; |
| 67 | + |
| 68 | + if (x_abs >= p) { |
| 69 | + if (xbits.is_nan()) |
| 70 | + return x; |
| 71 | + if (xbits.is_inf()) { |
| 72 | + fputil::set_errno_if_required(EDOM); |
| 73 | + fputil::raise_except_if_required(FE_INVALID); |
| 74 | + return FPBits::quiet_nan().get_val(); |
| 75 | + } |
| 76 | + return FPBits::zero(xbits.sign()).get_val(); |
| 77 | + } |
| 78 | + double n = pow(2, -52); |
| 79 | + double k = fputil::nearest_integer(x * n); |
| 80 | + double y = x - k; |
| 81 | + double sin_k, cos_k, sin_y, cosm1_y; |
| 82 | + |
| 83 | + sincospi_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y); |
| 84 | + |
| 85 | + if (LIBC_UNLIKELY(sin_y == 0 && cos_k == 0)) |
| 86 | + return FPBits::zero(xbits.sign()).get_val(); |
| 87 | + |
| 88 | + return fputil::cast<double>(fputil::multiply_add( |
| 89 | + cos_k, cosm1_y, fputil::multiply_add(-sin_k, sin_y, cos_k))); |
| 90 | +} |
| 91 | +} // namespace LIBC_NAMESPACE_DECL |
0 commit comments