128 lines · cpp
1//===-- Half-precision asinpif16(x) 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/asinpif16.h"10#include "hdr/errno_macros.h"11#include "hdr/fenv_macros.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/cast.h"16#include "src/__support/FPUtil/except_value_utils.h"17#include "src/__support/FPUtil/multiply_add.h"18#include "src/__support/FPUtil/sqrt.h"19#include "src/__support/macros/optimization.h"20 21namespace LIBC_NAMESPACE_DECL {22 23LLVM_LIBC_FUNCTION(float16, asinpif16, (float16 x)) {24 using FPBits = fputil::FPBits<float16>;25 26 FPBits xbits(x);27 bool is_neg = xbits.is_neg();28 double x_abs = fputil::cast<double>(xbits.abs().get_val());29 30 auto signed_result = [is_neg](auto r) -> auto { return is_neg ? -r : r; };31 32 if (LIBC_UNLIKELY(x_abs > 1.0)) {33 // aspinf16(NaN) = NaN34 if (xbits.is_nan()) {35 if (xbits.is_signaling_nan()) {36 fputil::raise_except_if_required(FE_INVALID);37 return FPBits::quiet_nan().get_val();38 }39 return x;40 }41 42 // 1 < |x| <= +/-inf43 fputil::raise_except_if_required(FE_INVALID);44 fputil::set_errno_if_required(EDOM);45 46 return FPBits::quiet_nan().get_val();47 }48 49 // the coefficients for the polynomial approximation of asin(x)/pi in the50 // range [0, 0.5] extracted using python-sympy51 //52 // Python code to generate the coefficients:53 // > from sympy import *54 // > import math55 // > x = symbols('x')56 // > print(series(asin(x)/math.pi, x, 0, 21))57 //58 // OUTPUT:59 //60 // 0.318309886183791*x + 0.0530516476972984*x**3 + 0.0238732414637843*x**5 +61 // 0.0142102627760621*x**7 + 0.00967087327815336*x**9 +62 // 0.00712127941391293*x**11 + 0.00552355646848375*x**13 +63 // 0.00444514782463692*x**15 + 0.00367705242846804*x**17 +64 // 0.00310721681820837*x**19 + O(x**21)65 //66 // it's very accurate in the range [0, 0.5] and has a maximum error of67 // 0.0000000000000001 in the range [0, 0.5].68 constexpr double POLY_COEFFS[] = {69 0x1.45f306dc9c889p-2, // x^170 0x1.b2995e7b7b5fdp-5, // x^371 0x1.8723a1d588a36p-6, // x^572 0x1.d1a452f20430dp-7, // x^773 0x1.3ce52a3a09f61p-7, // x^974 0x1.d2b33e303d375p-8, // x^1175 0x1.69fde663c674fp-8, // x^1376 0x1.235134885f19bp-8, // x^1577 };78 // polynomial evaluation using horner's method79 // work only for |x| in [0, 0.5]80 auto asinpi_polyeval = [&](double x) -> double {81 return x * fputil::polyeval(x * x, POLY_COEFFS[0], POLY_COEFFS[1],82 POLY_COEFFS[2], POLY_COEFFS[3], POLY_COEFFS[4],83 POLY_COEFFS[5], POLY_COEFFS[6], POLY_COEFFS[7]);84 };85 86 // if |x| <= 0.5:87 if (LIBC_UNLIKELY(x_abs <= 0.5)) {88 // Use polynomial approximation of asin(x)/pi in the range [0, 0.5]89 double result = asinpi_polyeval(fputil::cast<double>(x));90 return fputil::cast<float16>(result);91 }92 93 // If |x| > 0.5, we need to use the range reduction method:94 // y = asin(x) => x = sin(y)95 // because: sin(a) = cos(pi/2 - a)96 // therefore:97 // x = cos(pi/2 - y)98 // let z = pi/2 - y,99 // x = cos(z)100 // because: cos(2a) = 1 - 2 * sin^2(a), z = 2a, a = z/2101 // therefore:102 // cos(z) = 1 - 2 * sin^2(z/2)103 // sin(z/2) = sqrt((1 - cos(z))/2)104 // sin(z/2) = sqrt((1 - x)/2)105 // let u = (1 - x)/2106 // then:107 // sin(z/2) = sqrt(u)108 // z/2 = asin(sqrt(u))109 // z = 2 * asin(sqrt(u))110 // pi/2 - y = 2 * asin(sqrt(u))111 // y = pi/2 - 2 * asin(sqrt(u))112 // y/pi = 1/2 - 2 * asin(sqrt(u))/pi113 //114 // Finally, we can write:115 // asinpi(x) = 1/2 - 2 * asinpi(sqrt(u))116 // where u = (1 - x) /2117 // = 0.5 - 0.5 * x118 // = multiply_add(-0.5, x, 0.5)119 120 double u = fputil::multiply_add(-0.5, x_abs, 0.5);121 double asinpi_sqrt_u = asinpi_polyeval(fputil::sqrt<double>(u));122 double result = fputil::multiply_add(-2.0, asinpi_sqrt_u, 0.5);123 124 return fputil::cast<float16>(signed_result(result));125}126 127} // namespace LIBC_NAMESPACE_DECL128