[libc-commits] [libc] [libc][math] Improve performance of double precision trig functions. (PR #111793)
Nick Desaulniers via libc-commits
libc-commits at lists.llvm.org
Thu Oct 10 08:52:09 PDT 2024
================
@@ -15,174 +15,62 @@
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/common.h"
#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/math/generic/range_reduction_double_common.h"
namespace LIBC_NAMESPACE_DECL {
-namespace fma {
-
using LIBC_NAMESPACE::fputil::DoubleDouble;
-LIBC_INLINE constexpr int FAST_PASS_EXPONENT = 32;
+LIBC_INLINE unsigned LargeRangeReduction::fast(double x, DoubleDouble &u) {
+ using FPBits = typename fputil::FPBits<double>;
+ FPBits xbits(x);
-// Digits of 2^(16*i) / pi, generated by Sollya with:
-// For [2..62]:
-// > for i from 3 to 63 do {
-// pi_inv = 2^(16*(i - 3)) / pi;
-// pn = nearestint(pi_inv);
-// pi_frac = pi_inv - pn;
-// a = round(pi_frac, D, RN);
-// b = round(pi_frac - a, D, RN);
-// c = round(pi_frac - a - b, D, RN);
-// d = round(pi_frac - a - b - c, D, RN);
-// print("{", 2^7 * a, ",", 2^7 * b, ",", 2^7 * c, ",", 2^7 * d, "},");
-// };
-// For [0..1]:
-// The leading bit of 2^(16*(i - 3)) / pi is very small, so we add 0.25 so that
-// the conditions for the algorithms are still satisfied, and one of those
-// conditions guarantees that ulp(0.25 * x_reduced) >= 2, and will safely be
-// discarded.
-// for i from 0 to 2 do {
-// pi_frac = 0.25 + 2^(16*(i - 3)) / pi;
-// a = round(pi_frac, D, RN);
-// b = round(pi_frac - a, D, RN);
-// c = round(pi_frac - a - b, D, RN);
-// d = round(pi_frac - a - b - c, D, RN);
-// print("{", 2^7 * a, ",", 2^7 * b, ",", 2^7 * c, ",", 2^7 * d, "},");
-// };
-// For The fast pass using double-double, we only need 3 parts (a, b, c), but
-// for the accurate pass using Float128, instead of using another table of
-// Float128s, we simply add the fourth path (a, b, c, d), which simplify the
-// implementation a bit and saving some memory.
-LIBC_INLINE constexpr double ONE_TWENTY_EIGHT_OVER_PI[64][4] = {
- {0x1.0000000000014p5, 0x1.7cc1b727220a9p-49, 0x1.3f84eafa3ea6ap-103,
- -0x1.11f924eb53362p-157},
- {0x1.0000000145f3p5, 0x1.b727220a94fe1p-49, 0x1.d5f47d4d37703p-104,
- 0x1.b6295993c439p-158},
- {0x1.000145f306dcap5, -0x1.bbead603d8a83p-50, 0x1.f534ddc0db629p-106,
- 0x1.664f10e4107f9p-160},
- {0x1.45f306dc9c883p5, -0x1.6b01ec5417056p-49, -0x1.6447e493ad4cep-103,
- 0x1.e21c820ff28b2p-157},
- {-0x1.f246c6efab581p4, 0x1.3abe8fa9a6eep-53, 0x1.b6c52b3278872p-107,
- 0x1.07f9458eaf7afp-164},
- {0x1.391054a7f09d6p4, -0x1.70565911f924fp-53, 0x1.2b3278872084p-107,
- -0x1.ae9c5421443aap-162},
- {0x1.529fc2757d1f5p2, 0x1.a6ee06db14acdp-53, -0x1.8778df7c035d4p-107,
- 0x1.d5ef5de2b0db9p-161},
- {-0x1.ec54170565912p-1, 0x1.b6c52b3278872p-59, 0x1.07f9458eaf7afp-116,
- -0x1.d4f246dc8e2dfp-173},
- {-0x1.505c1596447e5p5, 0x1.b14acc9e21c82p-49, 0x1.fe5163abdebbcp-106,
- 0x1.586dc91b8e909p-160},
- {-0x1.596447e493ad5p1, 0x1.93c439041fe51p-54, 0x1.8eaf7aef1586ep-108,
- -0x1.b7238b7b645a4p-163},
- {0x1.bb81b6c52b328p5, -0x1.de37df00d74e3p-49, 0x1.7bd778ac36e49p-103,
- -0x1.1c5bdb22d1ffap-158},
- {0x1.b6c52b3278872p5, 0x1.07f9458eaf7afp-52, -0x1.d4f246dc8e2dfp-109,
- 0x1.374b801924bbbp-164},
- {0x1.2b3278872084p5, -0x1.ae9c5421443aap-50, 0x1.b7246e3a424ddp-106,
- 0x1.700324977504fp-161},
- {-0x1.8778df7c035d4p5, 0x1.d5ef5de2b0db9p-49, 0x1.1b8e909374b8p-104,
- 0x1.924bba8274648p-160},
- {-0x1.bef806ba71508p4, -0x1.443a9e48db91cp-50, -0x1.6f6c8b47fe6dbp-104,
- -0x1.115f62e6de302p-158},
- {-0x1.ae9c5421443aap-2, 0x1.b7246e3a424ddp-58, 0x1.700324977504fp-113,
- -0x1.cdbc603c429c7p-167},
- {-0x1.38a84288753c9p5, -0x1.b7238b7b645a4p-51, 0x1.924bba8274648p-112,
- 0x1.cfe1deb1cb12ap-166},
- {-0x1.0a21d4f246dc9p3, 0x1.d2126e9700325p-53, -0x1.a22bec5cdbc6p-107,
- -0x1.e214e34ed658cp-162},
- {-0x1.d4f246dc8e2dfp3, 0x1.374b801924bbbp-52, -0x1.f62e6de301e21p-106,
- -0x1.38d3b5963045ep-160},
- {-0x1.236e4716f6c8bp4, -0x1.1ff9b6d115f63p-50, 0x1.921cfe1deb1cbp-106,
- 0x1.29a73ee88235fp-162},
- {0x1.b8e909374b802p4, -0x1.b6d115f62e6dep-50, -0x1.80f10a71a76b3p-105,
- 0x1.cfba208d7d4bbp-160},
- {0x1.09374b801924cp4, -0x1.15f62e6de301ep-50, -0x1.0a71a76b2c609p-105,
- 0x1.1046bea5d7689p-159},
- {-0x1.68ffcdb688afbp3, -0x1.736f180f10a72p-53, 0x1.62534e7dd1047p-107,
- -0x1.0568a25dbd8b3p-161},
- {0x1.924bba8274648p0, 0x1.cfe1deb1cb12ap-54, -0x1.63045df7282b4p-108,
- -0x1.44bb7b16638fep-162},
- {-0x1.a22bec5cdbc6p5, -0x1.e214e34ed658cp-50, -0x1.177dca0ad144cp-106,
- 0x1.213a671c09ad1p-160},
- {0x1.3a32439fc3bd6p1, 0x1.cb129a73ee882p-54, 0x1.afa975da24275p-109,
- -0x1.8e3f652e8207p-164},
- {-0x1.b78c0788538d4p4, 0x1.29a73ee88235fp-50, 0x1.4baed1213a672p-104,
- -0x1.fb29741037d8dp-159},
- {0x1.fc3bd63962535p5, -0x1.822efb9415a29p-51, 0x1.a24274ce38136p-105,
- -0x1.741037d8cdc54p-159},
- {-0x1.4e34ed658c117p2, -0x1.f7282b4512edfp-52, 0x1.d338e04d68bfp-107,
- -0x1.bec66e29c67cbp-162},
- {0x1.62534e7dd1047p5, -0x1.0568a25dbd8b3p-49, -0x1.c7eca5d040df6p-105,
- -0x1.9b8a719f2b318p-160},
- {-0x1.63045df7282b4p4, -0x1.44bb7b16638fep-50, 0x1.ad17df904e647p-104,
- 0x1.639835339f49dp-158},
- {0x1.d1046bea5d769p5, -0x1.bd8b31c7eca5dp-49, -0x1.037d8cdc538dp-107,
- 0x1.a99cfa4e422fcp-161},
- {0x1.afa975da24275p3, -0x1.8e3f652e8207p-52, 0x1.3991d63983534p-106,
- -0x1.82d8dee81d108p-160},
- {-0x1.a28976f62cc72p5, 0x1.35a2fbf209cc9p-53, -0x1.4e33e566305b2p-109,
- 0x1.08bf177bf2507p-163},
- {-0x1.76f62cc71fb29p5, -0x1.d040df633714ep-49, -0x1.9f2b3182d8defp-104,
- 0x1.f8bbdf9283b2p-158},
- {0x1.d338e04d68bfp5, -0x1.bec66e29c67cbp-50, 0x1.9cfa4e422fc5ep-105,
- -0x1.036be27003b4p-161},
- {0x1.c09ad17df904ep4, 0x1.91d639835339fp-50, 0x1.272117e2ef7e5p-104,
- -0x1.7c4e007680022p-158},
- {0x1.68befc827323bp5, -0x1.c67cacc60b638p-50, 0x1.17e2ef7e4a0ecp-104,
- 0x1.ff897ffde0598p-158},
- {-0x1.037d8cdc538dp5, 0x1.a99cfa4e422fcp-49, 0x1.77bf250763ff1p-103,
- 0x1.7ffde05980fefp-158},
- {-0x1.8cdc538cf9599p5, 0x1.f49c845f8bbep-50, -0x1.b5f13801da001p-104,
- 0x1.e05980fef2f12p-158},
- {-0x1.4e33e566305b2p3, 0x1.08bf177bf2507p-51, 0x1.8ffc4bffef02dp-105,
- -0x1.fc04343b9d298p-160},
- {-0x1.f2b3182d8dee8p4, -0x1.d1081b5f13802p-52, 0x1.2fffbc0b301fep-107,
- -0x1.a1dce94beb25cp-163},
- {-0x1.8c16c6f740e88p5, -0x1.036be27003b4p-49, -0x1.0fd33f8086877p-109,
- -0x1.d297d64b824b2p-164},
- {0x1.3908bf177bf25p5, 0x1.d8ffc4bffef03p-53, -0x1.9fc04343b9d29p-108,
- -0x1.f592e092c9813p-162},
- {0x1.7e2ef7e4a0ec8p4, -0x1.da00087e99fcp-56, -0x1.0d0ee74a5f593p-110,
- 0x1.f6d367ecf27cbp-166},
- {-0x1.081b5f13801dap4, -0x1.0fd33f8086877p-61, -0x1.d297d64b824b2p-116,
- -0x1.8130d834f648bp-170},
- {-0x1.af89c00ed0004p5, -0x1.fa67f010d0ee7p-50, -0x1.297d64b824b26p-104,
- -0x1.30d834f648b0cp-162},
- {-0x1.c00ed00043f4dp5, 0x1.fde5e2316b415p-55, -0x1.2e092c98130d8p-110,
- -0x1.a7b24585ce04dp-165},
- {0x1.2fffbc0b301fep5, -0x1.a1dce94beb25cp-51, -0x1.25930261b069fp-107,
- 0x1.b74f463f669e6p-162},
- {-0x1.0fd33f8086877p3, -0x1.d297d64b824b2p-52, -0x1.8130d834f648bp-106,
- -0x1.738132c3402bap-163},
- {-0x1.9fc04343b9d29p4, -0x1.f592e092c9813p-50, -0x1.b069ec9161738p-107,
- -0x1.32c3402ba515bp-163},
- {-0x1.0d0ee74a5f593p2, 0x1.f6d367ecf27cbp-54, 0x1.36e9e8c7ecd3dp-111,
- -0x1.00ae9456c229cp-165},
- {-0x1.dce94beb25c12p5, -0x1.64c0986c1a7b2p-49, -0x1.161738132c34p-103,
- -0x1.5d28ad8453814p-158},
- {-0x1.4beb25c12593p5, -0x1.30d834f648b0cp-50, 0x1.8fd9a797fa8b6p-104,
- -0x1.5b08a7028341dp-159},
- {0x1.b47db4d9fb3cap4, -0x1.a7b24585ce04dp-53, 0x1.3cbfd45aea4f7p-107,
- 0x1.63f5f2f8bd9e8p-161},
- {-0x1.25930261b069fp5, 0x1.b74f463f669e6p-50, -0x1.5d28ad8453814p-110,
- -0x1.a0e84c2f8c608p-166},
- {0x1.fb3c9f2c26dd4p4, -0x1.738132c3402bap-51, -0x1.456c229c0a0dp-105,
- -0x1.d0985f18c10ebp-159},
- {-0x1.b069ec9161738p5, -0x1.32c3402ba515bp-51, -0x1.14e050683a131p-108,
- 0x1.0739f78a5292fp-162},
- {-0x1.ec9161738132cp5, -0x1.a015d28ad8454p-50, 0x1.faf97c5ecf41dp-104,
- -0x1.821d6b5b4565p-160},
- {-0x1.61738132c3403p5, 0x1.16ba93dd63f5fp-49, 0x1.7c5ecf41ce7dep-104,
- 0x1.4a525d4d7f6bfp-159},
- {0x1.fb34f2ff516bbp3, -0x1.b08a7028341d1p-51, 0x1.9e839cfbc5295p-105,
- -0x1.a2b2809409dc1p-159},
- {0x1.3cbfd45aea4f7p5, 0x1.63f5f2f8bd9e8p-49, 0x1.ce7de294a4baap-104,
- -0x1.404a04ee072a3p-158},
- {-0x1.5d28ad8453814p2, -0x1.a0e84c2f8c608p-54, -0x1.d6b5b45650128p-108,
- -0x1.3b81ca8bdea7fp-164},
- {-0x1.15b08a7028342p5, 0x1.7b3d0739f78a5p-50, 0x1.497535fdafd89p-105,
- -0x1.ca8bdea7f33eep-164},
-};
+ int x_e_m62 = xbits.get_biased_exponent() - (FPBits::EXP_BIAS + 62);
+ idx = static_cast<unsigned>((x_e_m62 >> 4) + 3);
+ // Scale x down by 2^(-(16 * (idx - 3))
+ xbits.set_biased_exponent((x_e_m62 & 15) + FPBits::EXP_BIAS + 62);
+ // 2^62 <= |x_reduced| < 2^(62 + 16) = 2^78
+ x_reduced = xbits.get_val();
+ // x * c_hi = ph.hi + ph.lo exactly.
+ DoubleDouble ph =
+ fputil::exact_mult<NO_FMA>(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][0]);
+ // x * c_mid = pm.hi + pm.lo exactly.
+ DoubleDouble pm =
+ fputil::exact_mult<NO_FMA>(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][1]);
+ // x * c_lo = pl.hi + pl.lo exactly.
+ DoubleDouble pl =
+ fputil::exact_mult<NO_FMA>(x_reduced, ONE_TWENTY_EIGHT_OVER_PI[idx][2]);
----------------
nickdesaulniers wrote:
Should these invocations of `exact_mult` (and on line 70) be using FMA rather than NO_FMA?
It's hard to tell what's different between this definition and the nofma.h version. Perhaps some code can be shared between the two, rather than duplicated?
https://github.com/llvm/llvm-project/pull/111793
More information about the libc-commits
mailing list