298 lines · c
1//===-- Implementation header for asin --------------------------*- 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-exception6//7//===----------------------------------------------------------------------===//8 9#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_ASIN_H10#define LLVM_LIBC_SRC___SUPPORT_MATH_ASIN_H11 12#include "asin_utils.h"13#include "src/__support/FPUtil/FEnvImpl.h"14#include "src/__support/FPUtil/FPBits.h"15#include "src/__support/FPUtil/double_double.h"16#include "src/__support/FPUtil/dyadic_float.h"17#include "src/__support/FPUtil/multiply_add.h"18#include "src/__support/FPUtil/sqrt.h"19#include "src/__support/macros/config.h"20#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY21#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA22#include "src/__support/math/asin_utils.h"23 24namespace LIBC_NAMESPACE_DECL {25 26namespace math {27 28LIBC_INLINE static constexpr double asin(double x) {29 using namespace asin_internal;30 using FPBits = fputil::FPBits<double>;31 32 FPBits xbits(x);33 int x_exp = xbits.get_biased_exponent();34 35 // |x| < 0.5.36 if (x_exp < FPBits::EXP_BIAS - 1) {37 // |x| < 2^-26.38 if (LIBC_UNLIKELY(x_exp < FPBits::EXP_BIAS - 26)) {39 // When |x| < 2^-26, the relative error of the approximation asin(x) ~ x40 // is:41 // |asin(x) - x| / |asin(x)| < |x^3| / (6|x|)42 // = x^2 / 643 // < 2^-5444 // < epsilon(1)/2.45 // So the correctly rounded values of asin(x) are:46 // = x + sign(x)*eps(x) if rounding mode = FE_TOWARDZERO,47 // or (rounding mode = FE_UPWARD and x is48 // negative),49 // = x otherwise.50 // To simplify the rounding decision and make it more efficient, we use51 // fma(x, 2^-54, x) instead.52 // Note: to use the formula x + 2^-54*x to decide the correct rounding, we53 // do need fma(x, 2^-54, x) to prevent underflow caused by 2^-54*x when54 // |x| < 2^-1022. For targets without FMA instructions, when x is close to55 // denormal range, we normalize x,56#if defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS)57 return x;58#elif defined(LIBC_TARGET_CPU_HAS_FMA_DOUBLE)59 return fputil::multiply_add(x, 0x1.0p-54, x);60#else61 if (xbits.abs().uintval() == 0)62 return x;63 // Get sign(x) * min_normal.64 FPBits eps_bits = FPBits::min_normal();65 eps_bits.set_sign(xbits.sign());66 double eps = eps_bits.get_val();67 double normalize_const = (x_exp == 0) ? eps : 0.0;68 double scaled_normal =69 fputil::multiply_add(x + normalize_const, 0x1.0p54, eps);70 return fputil::multiply_add(scaled_normal, 0x1.0p-54, -normalize_const);71#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS72 }73 74#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS75 return x * asin_eval(x * x);76#else77 using Float128 = fputil::DyadicFloat<128>;78 using DoubleDouble = fputil::DoubleDouble;79 80 unsigned idx = 0;81 DoubleDouble x_sq = fputil::exact_mult(x, x);82 double err = xbits.abs().get_val() * 0x1.0p-51;83 // Polynomial approximation:84 // p ~ asin(x)/x85 86 DoubleDouble p = asin_eval(x_sq, idx, err);87 // asin(x) ~ x * (ASIN_COEFFS[idx][0] + p)88 DoubleDouble r0 = fputil::exact_mult(x, p.hi);89 double r_lo = fputil::multiply_add(x, p.lo, r0.lo);90 91 // Ziv's accuracy test.92 93 double r_upper = r0.hi + (r_lo + err);94 double r_lower = r0.hi + (r_lo - err);95 96 if (LIBC_LIKELY(r_upper == r_lower))97 return r_upper;98 99 // Ziv's accuracy test failed, perform 128-bit calculation.100 101 // Recalculate mod 1/64.102 idx = static_cast<unsigned>(fputil::nearest_integer(x_sq.hi * 0x1.0p6));103 104 // Get x^2 - idx/64 exactly. When FMA is available, double-double105 // multiplication will be correct for all rounding modes. Otherwise we use106 // Float128 directly.107 Float128 x_f128(x);108 109#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE110 // u = x^2 - idx/64111 Float128 u_hi(112 fputil::multiply_add(static_cast<double>(idx), -0x1.0p-6, x_sq.hi));113 Float128 u = fputil::quick_add(u_hi, Float128(x_sq.lo));114#else115 Float128 x_sq_f128 = fputil::quick_mul(x_f128, x_f128);116 Float128 u = fputil::quick_add(117 x_sq_f128, Float128(static_cast<double>(idx) * (-0x1.0p-6)));118#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE119 120 Float128 p_f128 = asin_eval(u, idx);121 Float128 r = fputil::quick_mul(x_f128, p_f128);122 123 return static_cast<double>(r);124#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS125 }126 // |x| >= 0.5127 128 double x_abs = xbits.abs().get_val();129 130 // Maintaining the sign:131 constexpr double SIGN[2] = {1.0, -1.0};132 double x_sign = SIGN[xbits.is_neg()];133 134 // |x| >= 1135 if (LIBC_UNLIKELY(x_exp >= FPBits::EXP_BIAS)) {136 // x = +-1, asin(x) = +- pi/2137 if (x_abs == 1.0) {138 // return +- pi/2139 return fputil::multiply_add(x_sign, PI_OVER_TWO.hi,140 x_sign * PI_OVER_TWO.lo);141 }142 // |x| > 1, return NaN.143 if (xbits.is_quiet_nan())144 return x;145 146 // Set domain error for non-NaN input.147 if (!xbits.is_nan())148 fputil::set_errno_if_required(EDOM);149 150 fputil::raise_except_if_required(FE_INVALID);151 return FPBits::quiet_nan().get_val();152 }153 154 // When |x| >= 0.5, we perform range reduction as follow:155 //156 // Assume further that 0.5 <= x < 1, and let:157 // y = asin(x)158 // We will use the double angle formula:159 // cos(2y) = 1 - 2 sin^2(y)160 // and the complement angle identity:161 // x = sin(y) = cos(pi/2 - y)162 // = 1 - 2 sin^2 (pi/4 - y/2)163 // So:164 // sin(pi/4 - y/2) = sqrt( (1 - x)/2 )165 // And hence:166 // pi/4 - y/2 = asin( sqrt( (1 - x)/2 ) )167 // Equivalently:168 // asin(x) = y = pi/2 - 2 * asin( sqrt( (1 - x)/2 ) )169 // Let u = (1 - x)/2, then:170 // asin(x) = pi/2 - 2 * asin( sqrt(u) )171 // Moreover, since 0.5 <= x < 1:172 // 0 < u <= 1/4, and 0 < sqrt(u) <= 0.5,173 // And hence we can reuse the same polynomial approximation of asin(x) when174 // |x| <= 0.5:175 // asin(x) ~ pi/2 - 2 * sqrt(u) * P(u),176 177 // u = (1 - |x|)/2178 double u = fputil::multiply_add(x_abs, -0.5, 0.5);179 // v_hi + v_lo ~ sqrt(u).180 // Let:181 // h = u - v_hi^2 = (sqrt(u) - v_hi) * (sqrt(u) + v_hi)182 // Then:183 // sqrt(u) = v_hi + h / (sqrt(u) + v_hi)184 // ~ v_hi + h / (2 * v_hi)185 // So we can use:186 // v_lo = h / (2 * v_hi).187 // Then,188 // asin(x) ~ pi/2 - 2*(v_hi + v_lo) * P(u)189 double v_hi = fputil::sqrt<double>(u);190 191#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS192 double p = asin_eval(u);193 double r = x_sign * fputil::multiply_add(-2.0 * v_hi, p, PI_OVER_TWO.hi);194 return r;195#else196 197#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE198 double h = fputil::multiply_add(v_hi, -v_hi, u);199#else200 DoubleDouble v_hi_sq = fputil::exact_mult(v_hi, v_hi);201 double h = (u - v_hi_sq.hi) - v_hi_sq.lo;202#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE203 204 // Scale v_lo and v_hi by 2 from the formula:205 // vh = v_hi * 2206 // vl = 2*v_lo = h / v_hi.207 double vh = v_hi * 2.0;208 double vl = h / v_hi;209 210 // Polynomial approximation:211 // p ~ asin(sqrt(u))/sqrt(u)212 unsigned idx = 0;213 double err = vh * 0x1.0p-51;214 215 DoubleDouble p = asin_eval(DoubleDouble{0.0, u}, idx, err);216 217 // Perform computations in double-double arithmetic:218 // asin(x) = pi/2 - (v_hi + v_lo) * (ASIN_COEFFS[idx][0] + p)219 DoubleDouble r0 = fputil::quick_mult(DoubleDouble{vl, vh}, p);220 DoubleDouble r = fputil::exact_add(PI_OVER_TWO.hi, -r0.hi);221 222 double r_lo = PI_OVER_TWO.lo - r0.lo + r.lo;223 224 // Ziv's accuracy test.225 226#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE227 double r_upper = fputil::multiply_add(228 r.hi, x_sign, fputil::multiply_add(r_lo, x_sign, err));229 double r_lower = fputil::multiply_add(230 r.hi, x_sign, fputil::multiply_add(r_lo, x_sign, -err));231#else232 r_lo *= x_sign;233 r.hi *= x_sign;234 double r_upper = r.hi + (r_lo + err);235 double r_lower = r.hi + (r_lo - err);236#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE237 238 if (LIBC_LIKELY(r_upper == r_lower))239 return r_upper;240 241 // Ziv's accuracy test failed, we redo the computations in Float128.242 // Recalculate mod 1/64.243 idx = static_cast<unsigned>(fputil::nearest_integer(u * 0x1.0p6));244 245 // After the first step of Newton-Raphson approximating v = sqrt(u), we have246 // that:247 // sqrt(u) = v_hi + h / (sqrt(u) + v_hi)248 // v_lo = h / (2 * v_hi)249 // With error:250 // sqrt(u) - (v_hi + v_lo) = h * ( 1/(sqrt(u) + v_hi) - 1/(2*v_hi) )251 // = -h^2 / (2*v * (sqrt(u) + v)^2).252 // Since:253 // (sqrt(u) + v_hi)^2 ~ (2sqrt(u))^2 = 4u,254 // we can add another correction term to (v_hi + v_lo) that is:255 // v_ll = -h^2 / (2*v_hi * 4u)256 // = -v_lo * (h / 4u)257 // = -vl * (h / 8u),258 // making the errors:259 // sqrt(u) - (v_hi + v_lo + v_ll) = O(h^3)260 // well beyond 128-bit precision needed.261 262 // Get the rounding error of vl = 2 * v_lo ~ h / vh263 // Get full product of vh * vl264#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE265 double vl_lo = fputil::multiply_add(-v_hi, vl, h) / v_hi;266#else267 DoubleDouble vh_vl = fputil::exact_mult(v_hi, vl);268 double vl_lo = ((h - vh_vl.hi) - vh_vl.lo) / v_hi;269#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE270 // vll = 2*v_ll = -vl * (h / (4u)).271 double t = h * (-0.25) / u;272 double vll = fputil::multiply_add(vl, t, vl_lo);273 // m_v = -(v_hi + v_lo + v_ll).274 Float128 m_v = fputil::quick_add(275 Float128(vh), fputil::quick_add(Float128(vl), Float128(vll)));276 m_v.sign = Sign::NEG;277 278 // Perform computations in Float128:279 // asin(x) = pi/2 - (v_hi + v_lo + vll) * P(u).280 Float128 y_f128(fputil::multiply_add(static_cast<double>(idx), -0x1.0p-6, u));281 282 Float128 p_f128 = asin_eval(y_f128, idx);283 Float128 r0_f128 = fputil::quick_mul(m_v, p_f128);284 Float128 r_f128 = fputil::quick_add(PI_OVER_TWO_F128, r0_f128);285 286 if (xbits.is_neg())287 r_f128.sign = Sign::NEG;288 289 return static_cast<double>(r_f128);290#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS291}292 293} // namespace math294 295} // namespace LIBC_NAMESPACE_DECL296 297#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASIN_H298