285 lines · c
1//===-- Implementation header for acos --------------------------*- 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_ACOS_H10#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOS_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 23namespace LIBC_NAMESPACE_DECL {24 25namespace math {26 27LIBC_INLINE static constexpr double acos(double x) {28 using DoubleDouble = fputil::DoubleDouble;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^-55.38 if (LIBC_UNLIKELY(x_exp < FPBits::EXP_BIAS - 55)) {39 // When |x| < 2^-55, acos(x) = pi/240#if defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS)41 return PI_OVER_TWO.hi;42#else43 // Force the evaluation and prevent constant propagation so that it44 // is rounded correctly for FE_UPWARD rounding mode.45 return (xbits.abs().get_val() + 0x1.0p-160) + PI_OVER_TWO.hi;46#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS47 }48 49#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS50 // acos(x) = pi/2 - asin(x)51 // = pi/2 - x * P(x^2)52 double p = asin_eval(x * x);53 return PI_OVER_TWO.hi + fputil::multiply_add(-x, p, PI_OVER_TWO.lo);54#else55 unsigned idx = 0;56 DoubleDouble x_sq = fputil::exact_mult(x, x);57 double err = xbits.abs().get_val() * 0x1.0p-51;58 // Polynomial approximation:59 // p ~ asin(x)/x60 DoubleDouble p = asin_eval(x_sq, idx, err);61 // asin(x) ~ x * p62 DoubleDouble r0 = fputil::exact_mult(x, p.hi);63 // acos(x) = pi/2 - asin(x)64 // ~ pi/2 - x * p65 // = pi/2 - x * (p.hi + p.lo)66 double r_hi = fputil::multiply_add(-x, p.hi, PI_OVER_TWO.hi);67 // Use Dekker's 2SUM algorithm to compute the lower part.68 double r_lo = ((PI_OVER_TWO.hi - r_hi) - r0.hi) - r0.lo;69 r_lo = fputil::multiply_add(-x, p.lo, r_lo + PI_OVER_TWO.lo);70 71 // Ziv's accuracy test.72 73 double r_upper = r_hi + (r_lo + err);74 double r_lower = r_hi + (r_lo - err);75 76 if (LIBC_LIKELY(r_upper == r_lower))77 return r_upper;78 79 // Ziv's accuracy test failed, perform 128-bit calculation.80 81 // Recalculate mod 1/64.82 idx = static_cast<unsigned>(fputil::nearest_integer(x_sq.hi * 0x1.0p6));83 84 // Get x^2 - idx/64 exactly. When FMA is available, double-double85 // multiplication will be correct for all rounding modes. Otherwise we use86 // Float128 directly.87 Float128 x_f128(x);88 89#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE90 // u = x^2 - idx/6491 Float128 u_hi(92 fputil::multiply_add(static_cast<double>(idx), -0x1.0p-6, x_sq.hi));93 Float128 u = fputil::quick_add(u_hi, Float128(x_sq.lo));94#else95 Float128 x_sq_f128 = fputil::quick_mul(x_f128, x_f128);96 Float128 u = fputil::quick_add(97 x_sq_f128, Float128(static_cast<double>(idx) * (-0x1.0p-6)));98#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE99 100 Float128 p_f128 = asin_eval(u, idx);101 // Flip the sign of x_f128 to perform subtraction.102 x_f128.sign = x_f128.sign.negate();103 Float128 r =104 fputil::quick_add(PI_OVER_TWO_F128, fputil::quick_mul(x_f128, p_f128));105 106 return static_cast<double>(r);107#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS108 }109 // |x| >= 0.5110 111 double x_abs = xbits.abs().get_val();112 113 // Maintaining the sign:114 constexpr double SIGN[2] = {1.0, -1.0};115 double x_sign = SIGN[xbits.is_neg()];116 // |x| >= 1117 if (LIBC_UNLIKELY(x_exp >= FPBits::EXP_BIAS)) {118 // x = +-1, asin(x) = +- pi/2119 if (x_abs == 1.0) {120 // x = 1, acos(x) = 0,121 // x = -1, acos(x) = pi122 return x == 1.0 ? 0.0 : fputil::multiply_add(-x_sign, PI.hi, PI.lo);123 }124 // |x| > 1, return NaN.125 if (xbits.is_quiet_nan())126 return x;127 128 // Set domain error for non-NaN input.129 if (!xbits.is_nan())130 fputil::set_errno_if_required(EDOM);131 132 fputil::raise_except_if_required(FE_INVALID);133 return FPBits::quiet_nan().get_val();134 }135 136 // When |x| >= 0.5, we perform range reduction as follow:137 //138 // When 0.5 <= x < 1, let:139 // y = acos(x)140 // We will use the double angle formula:141 // cos(2y) = 1 - 2 sin^2(y)142 // and the complement angle identity:143 // x = cos(y) = 1 - 2 sin^2 (y/2)144 // So:145 // sin(y/2) = sqrt( (1 - x)/2 )146 // And hence:147 // y/2 = asin( sqrt( (1 - x)/2 ) )148 // Equivalently:149 // acos(x) = y = 2 * asin( sqrt( (1 - x)/2 ) )150 // Let u = (1 - x)/2, then:151 // acos(x) = 2 * asin( sqrt(u) )152 // Moreover, since 0.5 <= x < 1:153 // 0 < u <= 1/4, and 0 < sqrt(u) <= 0.5,154 // And hence we can reuse the same polynomial approximation of asin(x) when155 // |x| <= 0.5:156 // acos(x) ~ 2 * sqrt(u) * P(u).157 //158 // When -1 < x <= -0.5, we reduce to the previous case using the formula:159 // acos(x) = pi - acos(-x)160 // = pi - 2 * asin ( sqrt( (1 + x)/2 ) )161 // ~ pi - 2 * sqrt(u) * P(u),162 // where u = (1 - |x|)/2.163 164 // u = (1 - |x|)/2165 double u = fputil::multiply_add(x_abs, -0.5, 0.5);166 // v_hi + v_lo ~ sqrt(u).167 // Let:168 // h = u - v_hi^2 = (sqrt(u) - v_hi) * (sqrt(u) + v_hi)169 // Then:170 // sqrt(u) = v_hi + h / (sqrt(u) + v_hi)171 // ~ v_hi + h / (2 * v_hi)172 // So we can use:173 // v_lo = h / (2 * v_hi).174 double v_hi = fputil::sqrt<double>(u);175 176#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS177 constexpr DoubleDouble CONST_TERM[2] = {{0.0, 0.0}, PI};178 DoubleDouble const_term = CONST_TERM[xbits.is_neg()];179 180 double p = asin_eval(u);181 double scale = x_sign * 2.0 * v_hi;182 double r = const_term.hi + fputil::multiply_add(scale, p, const_term.lo);183 return r;184#else185 186#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE187 double h = fputil::multiply_add(v_hi, -v_hi, u);188#else189 DoubleDouble v_hi_sq = fputil::exact_mult(v_hi, v_hi);190 double h = (u - v_hi_sq.hi) - v_hi_sq.lo;191#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE192 193 // Scale v_lo and v_hi by 2 from the formula:194 // vh = v_hi * 2195 // vl = 2*v_lo = h / v_hi.196 double vh = v_hi * 2.0;197 double vl = h / v_hi;198 199 // Polynomial approximation:200 // p ~ asin(sqrt(u))/sqrt(u)201 unsigned idx = 0;202 double err = vh * 0x1.0p-51;203 204 DoubleDouble p = asin_eval(DoubleDouble{0.0, u}, idx, err);205 206 // Perform computations in double-double arithmetic:207 // asin(x) = pi/2 - (v_hi + v_lo) * (ASIN_COEFFS[idx][0] + p)208 DoubleDouble r0 = fputil::quick_mult(DoubleDouble{vl, vh}, p);209 210 double r_hi = 0, r_lo = 0;211 if (xbits.is_pos()) {212 r_hi = r0.hi;213 r_lo = r0.lo;214 } else {215 DoubleDouble r = fputil::exact_add(PI.hi, -r0.hi);216 r_hi = r.hi;217 r_lo = (PI.lo - r0.lo) + r.lo;218 }219 220 // Ziv's accuracy test.221 222 double r_upper = r_hi + (r_lo + err);223 double r_lower = r_hi + (r_lo - err);224 225 if (LIBC_LIKELY(r_upper == r_lower))226 return r_upper;227 228 // Ziv's accuracy test failed, we redo the computations in Float128.229 // Recalculate mod 1/64.230 idx = static_cast<unsigned>(fputil::nearest_integer(u * 0x1.0p6));231 232 // After the first step of Newton-Raphson approximating v = sqrt(u), we have233 // that:234 // sqrt(u) = v_hi + h / (sqrt(u) + v_hi)235 // v_lo = h / (2 * v_hi)236 // With error:237 // sqrt(u) - (v_hi + v_lo) = h * ( 1/(sqrt(u) + v_hi) - 1/(2*v_hi) )238 // = -h^2 / (2*v * (sqrt(u) + v)^2).239 // Since:240 // (sqrt(u) + v_hi)^2 ~ (2sqrt(u))^2 = 4u,241 // we can add another correction term to (v_hi + v_lo) that is:242 // v_ll = -h^2 / (2*v_hi * 4u)243 // = -v_lo * (h / 4u)244 // = -vl * (h / 8u),245 // making the errors:246 // sqrt(u) - (v_hi + v_lo + v_ll) = O(h^3)247 // well beyond 128-bit precision needed.248 249 // Get the rounding error of vl = 2 * v_lo ~ h / vh250 // Get full product of vh * vl251#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE252 double vl_lo = fputil::multiply_add(-v_hi, vl, h) / v_hi;253#else254 DoubleDouble vh_vl = fputil::exact_mult(v_hi, vl);255 double vl_lo = ((h - vh_vl.hi) - vh_vl.lo) / v_hi;256#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE257 // vll = 2*v_ll = -vl * (h / (4u)).258 double t = h * (-0.25) / u;259 double vll = fputil::multiply_add(vl, t, vl_lo);260 // m_v = -(v_hi + v_lo + v_ll).261 Float128 m_v = fputil::quick_add(262 Float128(vh), fputil::quick_add(Float128(vl), Float128(vll)));263 m_v.sign = xbits.sign();264 265 // Perform computations in Float128:266 // acos(x) = (v_hi + v_lo + vll) * P(u) , when 0.5 <= x < 1,267 // = pi - (v_hi + v_lo + vll) * P(u) , when -1 < x <= -0.5.268 Float128 y_f128(fputil::multiply_add(static_cast<double>(idx), -0x1.0p-6, u));269 270 Float128 p_f128 = asin_eval(y_f128, idx);271 Float128 r_f128 = fputil::quick_mul(m_v, p_f128);272 273 if (xbits.is_neg())274 r_f128 = fputil::quick_add(PI_F128, r_f128);275 276 return static_cast<double>(r_f128);277#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS278}279 280} // namespace math281 282} // namespace LIBC_NAMESPACE_DECL283 284#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOS_H285