-
Notifications
You must be signed in to change notification settings - Fork 14.5k
[libc][math][c23] implement C23 math function asinpif16 #146226
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: main
Are you sure you want to change the base?
Conversation
@llvm/pr-subscribers-libc Author: Mohamed Emad (hulxv) ChangesThe function is implemented using the following Taylor series that's generated using python-sympy, and it is very accurate for x $$ Geogebra graphCloses #132210 Full diff: https://github.com/llvm/llvm-project/pull/146226.diff 10 Files Affected:
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index be5f5a66016b5..240d48f5bcdb8 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -651,6 +651,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
list(APPEND TARGET_LIBM_ENTRYPOINTS
# math.h C23 _Float16 entrypoints
# libc.src.math.acoshf16
+ libc.src.math.asinpif16
libc.src.math.canonicalizef16
libc.src.math.ceilf16
libc.src.math.copysignf16
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 73dfeae1a2c94..c9e5b0096099f 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -664,6 +664,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.acoshf16
libc.src.math.asinf16
libc.src.math.asinhf16
+ libc.src.math.asinpif16
libc.src.math.canonicalizef16
libc.src.math.ceilf16
libc.src.math.copysignf16
diff --git a/libc/docs/headers/math/index.rst b/libc/docs/headers/math/index.rst
index 947bd4b60b391..aaacc3ca92832 100644
--- a/libc/docs/headers/math/index.rst
+++ b/libc/docs/headers/math/index.rst
@@ -259,7 +259,7 @@ Higher Math Functions
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| asinh | |check| | | | |check| | | 7.12.5.2 | F.10.2.2 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| asinpi | | | | | | 7.12.4.9 | F.10.1.9 |
+| asinpi | | | | |check| | | 7.12.4.9 | F.10.1.9 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| atan | |check| | 1 ULP | | | | 7.12.4.3 | F.10.1.3 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index fef829422244d..b49c17b871a93 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -71,7 +71,14 @@ functions:
- stdc
return_type: double
arguments:
- - type: double
+ - type: double
+ - name: asinpif16
+ standards:
+ - stdc
+ return_type: _Float16
+ arguments:
+ - type: _Float16
+ guard: LIBC_TYPES_HAS_FLOAT16
- name: atan2
standards:
- stdc
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index d177ff79141c0..e7ccb3726c54c 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -56,6 +56,8 @@ add_math_entrypoint_object(asinh)
add_math_entrypoint_object(asinhf)
add_math_entrypoint_object(asinhf16)
+add_math_entrypoint_object(asinpif16)
+
add_math_entrypoint_object(atan)
add_math_entrypoint_object(atanf)
diff --git a/libc/src/math/asinpif16.h b/libc/src/math/asinpif16.h
new file mode 100644
index 0000000000000..67ccb4ff4ac3d
--- /dev/null
+++ b/libc/src/math/asinpif16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for asinpif16 ---------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception.
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_ASINFPI16_H
+#define LLVM_LIBC_SRC_MATH_ASINFPI16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 asinpif16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_ASINPIF16_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index adbed5b2de48c..76accf72058d3 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4017,6 +4017,22 @@ add_entrypoint_object(
libc.src.__support.macros.properties.types
)
+add_entrypoint_object(
+ asinpif16
+ SRCS
+ asinpif16.cpp
+ HDRS
+ ../asinpif16.h
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.__support.FPUtil.cast
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.macros.optimization
+)
+
add_entrypoint_object(
atanhf
SRCS
diff --git a/libc/src/math/generic/asinpif16.cpp b/libc/src/math/generic/asinpif16.cpp
new file mode 100644
index 0000000000000..34b0b7e7779fc
--- /dev/null
+++ b/libc/src/math/generic/asinpif16.cpp
@@ -0,0 +1,177 @@
+//===-- Half-precision asinf16(x) function --------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception.
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/asinpif16.h"
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/except_value_utils.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/sqrt.h"
+#include "src/__support/macros/optimization.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+static constexpr float16 ONE_OVER_TWO = 0.5f16;
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+static constexpr size_t N_ASINFPI_EXCEPTS = 3;
+static constexpr float16 ONE_OVER_SIX = 0.166748046875f16;
+
+static constexpr fputil::ExceptValues<float16, N_ASINFPI_EXCEPTS>
+ ASINFPI_EXCEPTS{{
+ // (input_hex, RZ_output_hex, RU_offset, RD_offset, RN_offset)
+ // x = 0.0, asinfpi(0.0) = 0.0
+ {0x0000, 0x0000, 0, 0, 0},
+
+ // x = 1.0, asinfpi(1.0) = 0.5
+ {(fputil::FPBits<float16>(1.0f16)).uintval(),
+ (fputil::FPBits<float16>(ONE_OVER_TWO)).uintval(), 0, 0, 0},
+
+ // x = 0.5, asinfpi(0.5) = 1/6
+ {(fputil::FPBits<float16>(0.5f16)).uintval(),
+ (fputil::FPBits<float16>(ONE_OVER_SIX)).uintval(), 0, 0, 0},
+
+ }};
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+LLVM_LIBC_FUNCTION(float16, asinpif16, (float16 x)) {
+ using FPBits = fputil::FPBits<float16>;
+
+ FPBits xbits(x);
+ uint16_t x_uint = xbits.uintval();
+ bool is_neg = static_cast<bool>(x_uint >> 15);
+ float16 x_abs = is_neg ? -x : x;
+
+ auto __signed_result = [is_neg](float16 r) -> float16 {
+ return is_neg ? -r : r;
+ };
+
+ if (LIBC_UNLIKELY(x_abs > 1.0f16)) {
+ // aspinf16(NaN) = NaN
+ if (xbits.is_nan()) {
+ if (xbits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+ return x;
+ }
+
+ // 1 < |x| <= +/-inf
+ fputil::raise_except_if_required(FE_INVALID);
+ fputil::set_errno_if_required(EDOM);
+
+ return FPBits::quiet_nan().get_val();
+ }
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // exceptional values
+ if (auto r = ASINFPI_EXCEPTS.lookup(x_uint); LIBC_UNLIKELY(r.has_value())) {
+ return (r.value());
+ }
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+ // zero
+ if (LIBC_UNLIKELY(x_abs == 0.0f16)) {
+ return 0.0f16;
+ }
+
+ // +/-1.0
+ // if x is +/-1.0, return +/-0.5
+ if (LIBC_UNLIKELY(x_abs == 1.0f16)) {
+ return fputil::cast<float16>(__signed_result(ONE_OVER_TWO));
+ }
+
+ // the coefficients for the polynomial approximation of asin(x)/pi in the
+ // range [0, 0.5] extracted using python-sympy
+ //
+ // Python code to generate the coefficients:
+ // > from sympy import *
+ // > import math
+ // > x = symbols('x')
+ // > print(series(asin(x)/math.pi, x, 0, 21))
+ //
+ // OUTPUT:
+ //
+ // 0.318309886183791*x + 0.0530516476972984*x**3 + 0.0238732414637843*x**5 +
+ // 0.0142102627760621*x**7 + 0.00967087327815336*x**9 +
+ // 0.00712127941391293*x**11 + 0.00552355646848375*x**13 +
+ // 0.00444514782463692*x**15 + 0.00367705242846804*x**17 +
+ // 0.00310721681820837*x**19 + O(x**21)
+ //
+ // it's very accurate in the range [0, 0.5] and has a maximum error of
+ // 0.0000000000000001 in the range [0, 0.5].
+ static constexpr float16 POLY_COEFFS[10] = {
+ 0.318309886183791f16, // x^1
+ 0.0530516476972984f16, // x^3
+ 0.0238732414637843f16, // x^5
+ 0.0142102627760621f16, // x^7
+ 0.00967087327815336f16, // x^9
+ 0.00712127941391293f16, // x^11
+ 0.00552355646848375f16, // x^13
+ 0.00444514782463692f16, // x^15
+ 0.00367705242846804f16, // x^17
+ 0.00310721681820837f16 // x^19
+ };
+
+ // polynomial evaluation using horner's method
+ // work only for |x| in [0, 0.5]
+ auto __asinpi_polyeval = [](float16 xsq) -> float16 {
+ return fputil::polyeval(xsq, POLY_COEFFS[0], POLY_COEFFS[1], POLY_COEFFS[2],
+ POLY_COEFFS[3], POLY_COEFFS[4], POLY_COEFFS[5],
+ POLY_COEFFS[6], POLY_COEFFS[7], POLY_COEFFS[8],
+ POLY_COEFFS[9]);
+ };
+
+ // if |x| <= 0.5:
+ if (LIBC_UNLIKELY(x_abs <= ONE_OVER_TWO)) {
+ // Use polynomial approximation of asin(x)/pi in the range [0, 0.5]
+ float16 xsq = x * x;
+ float16 result = x * __asinpi_polyeval(xsq);
+ return fputil::cast<float16>(__signed_result(result));
+ }
+
+ // If |x| > 0.5, we need to use the range reduction method:
+ // y = asin(x) => x = sin(y)
+ // because: sin(a) = cos(pi/2 - a)
+ // therefore:
+ // x = cos(pi/2 - y)
+ // let z = pi/2 - y,
+ // x = cos(z)
+ // becuase: cos(2a) = 1 - 2 * sin^2(a), z = 2a, a = z/2
+ // therefore:
+ // cos(z) = 1 - 2 * sin^2(z/2)
+ // sin(z/2) = sqrt((1 - cos(z))/2)
+ // sin(z/2) = sqrt((1 - x)/2)
+ // let u = (1 - x)/2
+ // then:
+ // sin(z/2) = sqrt(u)
+ // z/2 = asin(sqrt(u))
+ // z = 2 * asin(sqrt(u))
+ // pi/2 - y = 2 * asin(sqrt(u))
+ // y = pi/2 - 2 * asin(sqrt(u))
+ // y/pi = 1/2 - 2 * asin(sqrt(u))/pi
+ //
+ // Finally, we can write:
+ // asinpi(x) = 1/2 - 2 * asinpi(sqrt(u))
+ // where u = (1 - x) /2
+ // = 0.5 - 0.5 * x
+ // = multiply_add(-0.5, x, 0.5)
+
+ float16 u = fputil::multiply_add(-ONE_OVER_TWO, x_abs, ONE_OVER_TWO);
+ float16 u_sqrt = fputil::sqrt<float16>(u);
+ float16 asinpi_sqrt_u = u_sqrt * __asinpi_polyeval(u);
+ float16 result = fputil::multiply_add(-2.0f16, asinpi_sqrt_u, ONE_OVER_TWO);
+
+ return fputil::cast<float16>(__signed_result(result));
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 7ee8b86135557..bd62b76817ceb 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2245,6 +2245,19 @@ add_fp_unittest(
libc.src.math.asinf16
)
+add_fp_unittest(
+ asinpif16_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ asinpif16_test.cpp
+ DEPENDS
+ libc.src.math.fabs
+ libc.src.math.asinpif16
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_fp_unittest(
acosf_test
NEED_MPFR
diff --git a/libc/test/src/math/asinpif16_test.cpp b/libc/test/src/math/asinpif16_test.cpp
new file mode 100644
index 0000000000000..6d8f76568902e
--- /dev/null
+++ b/libc/test/src/math/asinpif16_test.cpp
@@ -0,0 +1,110 @@
+//===-- Unittests for asinpif16 -------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/errno/libc_errno.h"
+#include "src/math/asinpif16.h"
+#include "src/math/fabs.h"
+#include "test/UnitTest/FPMatcher.h"
+
+#include <errno.h>
+#include <stdint.h>
+
+using LlvmLibcAsinpif16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcAsinpif16Test, SpecialNumbers) {
+ using FPBits = LIBC_NAMESPACE::fputil::FPBits<float16>;
+
+ // zero
+ EXPECT_FP_EQ(0.0f16, LIBC_NAMESPACE::asinpif16(0.0f16));
+
+ // +/-1
+ EXPECT_FP_EQ(float16(0.5), LIBC_NAMESPACE::asinpif16(float16(1.0)));
+ EXPECT_FP_EQ(float16(-0.5), LIBC_NAMESPACE::asinpif16(float16(-1.0)));
+
+ // NaN inputs
+ EXPECT_FP_EQ(FPBits::quiet_nan().get_val(),
+ LIBC_NAMESPACE::asinpif16(FPBits::quiet_nan().get_val()));
+
+ EXPECT_FP_EQ(FPBits::quiet_nan().get_val(),
+ LIBC_NAMESPACE::asinpif16(FPBits::signaling_nan().get_val()));
+
+ // infinity inputs -> should return NaN
+ errno = 0;
+ EXPECT_FP_EQ(FPBits::quiet_nan().get_val(), LIBC_NAMESPACE::asinpif16(inf));
+ EXPECT_MATH_ERRNO(EDOM);
+
+ errno = 0;
+ EXPECT_FP_EQ(FPBits::quiet_nan().get_val(),
+ LIBC_NAMESPACE::asinpif16(neg_inf));
+ EXPECT_MATH_ERRNO(EDOM);
+}
+
+TEST_F(LlvmLibcAsinpif16Test, OutOfRange) {
+ using FPBits = LIBC_NAMESPACE::fputil::FPBits<float16>;
+
+ // Test values > 1
+ errno = 0;
+ EXPECT_FP_EQ(FPBits::quiet_nan().get_val(),
+ LIBC_NAMESPACE::asinpif16(float16(1.5)));
+ EXPECT_MATH_ERRNO(EDOM);
+
+ errno = 0;
+ EXPECT_FP_EQ(FPBits::quiet_nan().get_val(),
+ LIBC_NAMESPACE::asinpif16(float16(2.0)));
+ EXPECT_MATH_ERRNO(EDOM);
+
+ // Test values < -1
+ errno = 0;
+ EXPECT_FP_EQ(FPBits::quiet_nan().get_val(),
+ LIBC_NAMESPACE::asinpif16(float16(-1.5)));
+ EXPECT_MATH_ERRNO(EDOM);
+
+ errno = 0;
+ EXPECT_FP_EQ(FPBits::quiet_nan().get_val(),
+ LIBC_NAMESPACE::asinpif16(float16(-2.0)));
+ EXPECT_MATH_ERRNO(EDOM);
+
+ // Test maximum normal value (should be > 1 for float16)
+ errno = 0;
+ EXPECT_FP_EQ(FPBits::quiet_nan().get_val(),
+ LIBC_NAMESPACE::asinpif16(FPBits::max_normal().get_val()));
+ EXPECT_MATH_ERRNO(EDOM);
+}
+
+TEST_F(LlvmLibcAsinpif16Test, SymmetryProperty) {
+ // Test that asinpi(-x) = -asinpi(x)
+ constexpr float16 test_vals[] = {0.1f16, 0.25f16, 0.5f16, 0.75f16,
+ 0.9f16, 0.99f16, 1.0f16};
+
+ for (float16 x : test_vals) {
+ if (x <= 1.0) {
+ float16 pos_result = LIBC_NAMESPACE::asinpif16(x);
+ float16 neg_result = LIBC_NAMESPACE::asinpif16(-x);
+
+ EXPECT_FP_EQ(pos_result, LIBC_NAMESPACE::fabs(neg_result));
+ }
+ }
+}
+
+TEST_F(LlvmLibcAsinpif16Test, RangeValidation) {
+ // Test that output is always in [-0.5, 0.5] for valid inputs
+ constexpr int num_tests = 1000;
+
+ for (int i = 0; i <= num_tests; ++i) {
+ float16 t = -1.0f16 + (2.0f16 * static_cast<float16>(i)) /
+ static_cast<float16>(num_tests);
+
+ if (LIBC_NAMESPACE::fabs(t) <= 1.0) {
+ float16 result = LIBC_NAMESPACE::asinpif16(t);
+
+ // should be in [-0.5, 0.5]
+ EXPECT_TRUE(result >= -0.5f16);
+ EXPECT_TRUE(result <= 0.5f16);
+ }
+ }
+}
|
✅ With the latest revision this PR passed the C/C++ code formatter. |
all tests work now |
libc/src/math/generic/asinpif16.cpp
Outdated
0.318309886183791, // x^1 | ||
0.0530516476972984, // x^3 | ||
0.0238732414637843, // x^5 | ||
0.0142102627760621, // x^7 | ||
0.00967087327815336, // x^9 | ||
0.00712127941391293, // x^11 | ||
0.00552355646848375, // x^13 | ||
0.00444514782463692, // x^15 | ||
0.00367705242846804, // x^17 | ||
0.00310721681820837 // x^19 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Use hexadecimal float constants.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Also can float16's asinpi get away with using float
as intermediate computational type? What should the degree of the polynomial be, and/or how many more exceptional points needed? You should be able to iterate quickly with the added tests to figure these question out. One of the main reasons is that sqrt
in double is significantly slower than sqrt
in float.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
When using float
instead of double
, the exhaustive test fails at many points with ULP error = 1
/home/mohamed/dev/contrib/llvm/llvm-project/libc/test/src/math/asinpif16_test.cpp:38: FAILURE
Failed to match __llvm_libc_21_0_0_git::asinpif16(x) against LIBC_NAMESPACE::testing::mpfr::get_mpfr_matcher<mpfr::Operation::Asinpi>( x, __llvm_libc_21_0_0_git::asinpif16(x), 0.5, mpfr::RoundingMode::Upward).
Match value not within tolerance value of MPFR result:
Input decimal: -0.99365234375000000000000000000000000000000000000000
Input bits: 0xBBF3 = (S: 1, E: 0x000E, M: 0x03F3)
Match decimal: -0.46386718750000000000000000000000000000000000000000
Match bits: 0xB76C = (S: 1, E: 0x000D, M: 0x036C)
MPFR result: -0.46411132812500000000000000000000000000000000000000
MPFR rounded: 0xB76D = (S: 1, E: 0x000D, M: 0x036D)
ULP error: 1.00000000000000000000000000000000000000000000000000
/home/mohamed/dev/contrib/llvm/llvm-project/libc/test/src/math/asinpif16_test.cpp:38: FAILURE
Failed to match __llvm_libc_21_0_0_git::asinpif16(x) against LIBC_NAMESPACE::testing::mpfr::get_mpfr_matcher<mpfr::Operation::Asinpi>( x, __llvm_libc_21_0_0_git::asinpif16(x), 0.5, mpfr::RoundingMode::Downward).
Match value not within tolerance value of MPFR result:
Input decimal: -0.99365234375000000000000000000000000000000000000000
Input bits: 0xBBF3 = (S: 1, E: 0x000E, M: 0x03F3)
Match decimal: -0.46411132812500000000000000000000000000000000000000
Match bits: 0xB76D = (S: 1, E: 0x000D, M: 0x036D)
MPFR result: -0.46435546875000000000000000000000000000000000000000
MPFR rounded: 0xB76E = (S: 1, E: 0x000D, M: 0x036E)
ULP error: 1.00000000000000000000000000000000000000000000000000
/home/mohamed/dev/contrib/llvm/llvm-project/libc/test/src/math/asinpif16_test.cpp:38: FAILURE
Failed to match __llvm_libc_21_0_0_git::asinpif16(x) against LIBC_NAMESPACE::testing::mpfr::get_mpfr_matcher<mpfr::Operation::Asinpi>( x, __llvm_libc_21_0_0_git::asinpif16(x), 0.5, mpfr::RoundingMode::TowardZero).
Match value not within tolerance value of MPFR result:
Input decimal: -0.99365234375000000000000000000000000000000000000000
Input bits: 0xBBF3 = (S: 1, E: 0x000E, M: 0x03F3)
Match decimal: -0.46386718750000000000000000000000000000000000000000
Match bits: 0xB76C = (S: 1, E: 0x000D, M: 0x036C)
MPFR result: -0.46411132812500000000000000000000000000000000000000
MPFR rounded: 0xB76D = (S: 1, E: 0x000D, M: 0x036D)
ULP error: 1.00000000000000000000000000000000000000000000000000
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Then that's ok, we can fine-tune it later.
// 0.0142102627760621*x**7 + 0.00967087327815336*x**9 + | ||
// 0.00712127941391293*x**11 + 0.00552355646848375*x**13 + | ||
// 0.00444514782463692*x**15 + 0.00367705242846804*x**17 + | ||
// 0.00310721681820837*x**19 + O(x**21) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
nit: Update comments to match with the polynomial that is used below.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This output comes from the python script above, As a clarification
9d18512
to
fcc4f8a
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
LGTM. Let wait for @overmighty whether he has any further comments.
The function is implemented using the following Taylor series that's generated using python-sympy, and it is very accurate for |x|$$\in [0, 0.5]$$ and has been verified using Geogebra. The range reduction is used for the rest range (0.5, 1].
Geogebra graph
Closes #132210