math/detail/Sincospi.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2023-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file corecel/math/detail/Sincospi.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011 #include <cstdint>
0012
0013 namespace celeritas
0014 {
0015 namespace detail
0016 {
0017 //---------------------------------------------------------------------------//
0018 /*!
0019  * Calculate sin and cosine of a fraction of pi.
0020  * From
0021  * https://stackoverflow.com/questions/42792939/implementation-of-sinpi-and-cospi-using-standard-c-math-library
0022  * from njuffa, (c) CC BY-SA
0023  *
0024  */
0025 inline void sincospi_impl(double a, double* sptr, double* cptr)
0026 {
0027     double c, r, s, t, az;
0028     std::int64_t i;
0029
0030     az = a * 0.0;  // must be evaluated with IEEE-754 semantics
0031     // for |a| >= 2**53, cospi(a) = 1.0, but cospi(Inf) = NaN
0032     a = (std::fabs(a) < 9.0071992547409920e+15) ? a : az;  // 0x1.0p53
0033     // reduce argument to primary approximation interval (-0.25, 0.25)
0034     r = std::nearbyint(a + a);  // must use IEEE-754 "to nearest" rounding
0035     i = static_cast<std::int64_t>(r);
0036     t = std::fma(-0.5, r, a);
0037     // compute core approximations
0038     s = t * t;
0039     // Approximate cos(pi*x) for x in [-0.25,0.25]
0040     r = -1.0369917389758117e-4;
0041     r = std::fma(r, s, 1.9294935641298806e-3);
0042     r = std::fma(r, s, -2.5806887942825395e-2);
0043     r = std::fma(r, s, 2.3533063028328211e-1);
0044     r = std::fma(r, s, -1.3352627688538006e+0);
0045     r = std::fma(r, s, 4.0587121264167623e+0);
0046     r = std::fma(r, s, -4.9348022005446790e+0);
0047     c = std::fma(r, s, 1.0000000000000000e+0);
0048     // Approximate sin(pi*x) for x in [-0.25,0.25]
0049     r = 4.6151442520157035e-4;
0050     r = std::fma(r, s, -7.3700183130883555e-3);
0051     r = std::fma(r, s, 8.2145868949323936e-2);
0052     r = std::fma(r, s, -5.9926452893214921e-1);
0053     r = std::fma(r, s, 2.5501640398732688e+0);
0054     r = std::fma(r, s, -5.1677127800499516e+0);
0055     s = s * t;
0056     r = r * s;
0057     s = std::fma(t, 3.1415926535897931e+0, r);
0058     // Map results according to quadrant
0059     if (i & 2)
0060     {
0061         s = 0.0 - s;  // must be evaluated with IEEE-754 semantics
0062         c = 0.0 - c;  // must be evaluated with IEEE-754 semantics
0063     }
0064     if (i & 1)
0065     {
0066         t = 0.0 - s;  // must be evaluated with IEEE-754 semantics
0067         s = c;
0068         c = t;
0069     }
0070     // IEEE-754: sinPi(+n) is +0 and sinPi(-n) is -0 for positive integers n
0071     if (a == std::floor(a))
0072         s = az;
0073     *sptr = s;
0074     *cptr = c;
0075 }
0076
0077 //---------------------------------------------------------------------------//
0078 /*!
0079  * Calculate sin and cosine of a fraction of pi.
0080  * From
0081  * https://stackoverflow.com/questions/42792939/implementation-of-sinpi-and-cospi-using-standard-c-math-library
0082  * from njuffa, (c) CC BY-SA
0083  */
0084 inline void sincospi_impl(float a, float* sptr, float* cptr)
0085 {
0086     // NOTE: some C++ library implementations (GCC?) don't define
0087     // single-precision functions in std namespace
0088     using namespace std;
0089
0090     float az, t, c, r, s;
0091     int32_t i;
0092
0093     az = a * 0.0f;  // Must be evaluated with IEEE-754 semantics
0094     // For |a| > 2**24, cospi(a) = 1.0f, but cospi(Inf) = NaN
0095     a = (fabsf(a) < 0x1.0p24f) ? a : az;
0096     r = nearbyintf(a + a);  // Must use IEEE-754 "to nearest" rounding
0097     i = static_cast<int32_t>(r);
0098     t = fmaf(-0.5f, r, a);
0099     // Compute core approximations
0100     s = t * t;
0101     // Approximate cos(pi*x) for x in [-0.25,0.25] */
0102     r = 0x1.d9e000p-3f;
0103     r = fmaf(r, s, -0x1.55c400p+0f);
0104     r = fmaf(r, s, 0x1.03c1cep+2f);
0105     r = fmaf(r, s, -0x1.3bd3ccp+2f);
0106     c = fmaf(r, s, 0x1.000000p+0f);
0107     // Approximate sin(pi*x) for x in [-0.25,0.25] */
0108     r = -0x1.310000p-1f;
0109     r = fmaf(r, s, 0x1.46737ep+1f);
0110     r = fmaf(r, s, -0x1.4abbfep+2f);
0111     r = (t * s) * r;
0112     s = fmaf(t, 0x1.921fb6p+1f, r);
0113     if (i & 2)
0114     {
0115         s = 0.0f - s;  // Must be evaluated with IEEE-754 semantics
0116         c = 0.0f - c;  // Must be evaluated with IEEE-754 semantics
0117     }
0118     if (i & 1)
0119     {
0120         t = 0.0f - s;  // Must be evaluated with IEEE-754 semantics
0121         s = c;
0122         c = t;
0123     }
0124     // IEEE-754: sinPi(+n) is +0 and sinPi(-n) is -0 for positive integers n
0125     if (a == floorf(a))
0126         s = az;
0127     *sptr = s;
0128     *cptr = c;
0129 }
0130
0131 //---------------------------------------------------------------------------//
0132 //! Lazy implementation of sin(x * pi)
0133 template<class T>
0134 CELER_FORCEINLINE T sinpi_impl(T v)
0135 {
0136     T result;
0137     T unused;
0138     sincospi_impl(v, &result, &unused);
0139     return result;
0140 }
0141
0142 //---------------------------------------------------------------------------//
0143 //! Lazy implementation of cos(x * pi)
0144 template<class T>
0145 CELER_FORCEINLINE T cospi_impl(T v)
0146 {
0147     T unused;
0148     T result;
0149     sincospi_impl(v, &unused, &result);
0150     return result;
0151 }
0152
0153 //---------------------------------------------------------------------------//
0154 //! Lazy implementation of sincos
0155 template<class T>
0156 CELER_FORCEINLINE void sincos_impl(T x, T* sptr, T* cptr)
0157 {
0158     *sptr = std::sin(x);
0159     *cptr = std::cos(x);
0160 }
0161
0162 //---------------------------------------------------------------------------//
0163
0164 CELER_FORCEINLINE void sincospif(float x, float* sptr, float* cptr)
0165 {
0166     return sincospi_impl(x, sptr, cptr);
0167 }
0168
0169 CELER_FORCEINLINE void sincosf(float x, float* sptr, float* cptr)
0170 {
0171     return sincos_impl(x, sptr, cptr);
0172 }
0173
0174 CELER_FORCEINLINE float sinpif(float x)
0175 {
0176     return sinpi_impl(x);
0177 }
0178
0179 CELER_FORCEINLINE float cospif(float x)
0180 {
0181     return cospi_impl(x);
0182 }
0183
0184 CELER_FORCEINLINE void sincospi(double x, double* sptr, double* cptr)
0185 {
0186     return sincospi_impl(x, sptr, cptr);
0187 }
0188
0189 CELER_FORCEINLINE void sincos(double x, double* sptr, double* cptr)
0190 {
0191     return sincos_impl(x, sptr, cptr);
0192 }
0193
0194 CELER_FORCEINLINE double sinpi(double x)
0195 {
0196     return sinpi_impl(x);
0197 }
0198
0199 CELER_FORCEINLINE double cospi(double x)
0200 {
0201     return cospi_impl(x);
0202 }
0203
0204 //---------------------------------------------------------------------------//
0205 }  // namespace detail
0206 }  // namespace celeritas