From 0269578d3eadaa117dc0c030f739fce09e59d8d2 Mon Sep 17 00:00:00 2001 From: Andreas Kling Date: Fri, 5 Feb 2021 12:27:00 +0100 Subject: [PATCH] LibM: Implement nextafter() and nexttoward() Patch from Anonymous. --- Userland/Libraries/LibM/math.cpp | 99 +++++++++++++++++++++++++++---- Userland/Tests/LibM/test-math.cpp | 97 ++++++++++++++++++++++++++++++ 2 files changed, 186 insertions(+), 10 deletions(-) diff --git a/Userland/Libraries/LibM/math.cpp b/Userland/Libraries/LibM/math.cpp index d10df79c44..9df74829b1 100644 --- a/Userland/Libraries/LibM/math.cpp +++ b/Userland/Libraries/LibM/math.cpp @@ -70,11 +70,14 @@ union FloatExtractor; template<> union FloatExtractor { static const int mantissa_bits = 52; + static const unsigned long long mantissa_max = (1ull << 52) - 1; static const int exponent_bias = 1023; + static const int exponent_bits = 11; + static const unsigned exponent_max = 2047; struct { unsigned long long mantissa : 52; unsigned exponent : 11; - int sign : 1; + unsigned sign : 1; }; double d; }; @@ -82,11 +85,14 @@ union FloatExtractor { template<> union FloatExtractor { static const int mantissa_bits = 23; + static const unsigned mantissa_max = (1 << 23) - 1; static const int exponent_bias = 127; + static const int exponent_bits = 8; + static const unsigned exponent_max = 255; struct { unsigned long long mantissa : 23; unsigned exponent : 8; - int sign : 1; + unsigned sign : 1; }; float d; }; @@ -152,6 +158,71 @@ static FloatType internal_to_integer(FloatType x, RoundingMode rounding_mode) return extractor.d; } +// This is much branchier than it really needs to be +template +static FloatType internal_nextafter(FloatType x, bool up) +{ + if (!isfinite(x)) + return x; + using Extractor = FloatExtractor; + Extractor extractor; + extractor.d = x; + if (x == 0) { + if (!extractor.sign) { + extractor.mantissa = 1; + extractor.sign = !up; + return extractor.d; + } + if (up) { + extractor.sign = false; + extractor.mantissa = 1; + return extractor.d; + } + extractor.mantissa = 1; + extractor.sign = up != extractor.sign; + return extractor.d; + } + if (up != extractor.sign) { + extractor.mantissa++; + if (!extractor.mantissa) { + // no need to normalize the mantissa as we just hit a power + // of two. + extractor.exponent++; + if (extractor.exponent == Extractor::exponent_max) { + extractor.exponent = Extractor::exponent_max - 1; + extractor.mantissa = Extractor::mantissa_max; + } + } + return extractor.d; + } + + if (!extractor.mantissa) { + if (extractor.exponent) { + extractor.exponent--; + extractor.mantissa = Extractor::mantissa_max; + } else { + extractor.d = 0; + } + return extractor.d; + } + + extractor.mantissa--; + if (extractor.mantissa != Extractor::mantissa_max) + return extractor.d; + if (extractor.exponent) { + extractor.exponent--; + // normalize + extractor.mantissa <<= 1; + } else { + if (extractor.sign) { + // Negative infinity + extractor.mantissa = 0; + extractor.exponent = Extractor::exponent_max; + } + } + return extractor.d; +} + extern "C" { double trunc(double x) NOEXCEPT @@ -634,14 +705,18 @@ double erfc(double x) NOEXCEPT return 1 - erf(x); } -double nextafter(double, double) NOEXCEPT +double nextafter(double x, double target) NOEXCEPT { - TODO(); + if (x == target) + return target; + return internal_nextafter(x, target >= x); } -float nextafterf(float, float) NOEXCEPT +float nextafterf(float x, float target) NOEXCEPT { - TODO(); + if (x == target) + return target; + return internal_nextafter(x, target >= x); } long double nextafterl(long double, long double) NOEXCEPT @@ -649,14 +724,18 @@ long double nextafterl(long double, long double) NOEXCEPT TODO(); } -double nexttoward(double, long double) NOEXCEPT +double nexttoward(double x, long double target) NOEXCEPT { - TODO(); + if (x == target) + return target; + return internal_nextafter(x, target >= x); } -float nexttowardf(float, long double) NOEXCEPT +float nexttowardf(float x, long double target) NOEXCEPT { - TODO(); + if (x == target) + return target; + return internal_nextafter(x, target >= x); } long double nexttowardl(long double, long double) NOEXCEPT diff --git a/Userland/Tests/LibM/test-math.cpp b/Userland/Tests/LibM/test-math.cpp index 0a82f318a0..f70d1f6c99 100644 --- a/Userland/Tests/LibM/test-math.cpp +++ b/Userland/Tests/LibM/test-math.cpp @@ -111,4 +111,101 @@ TEST_CASE(logarithms) EXPECT_CLOSE(log10(5), 0.698988) } +union Extractor { + explicit Extractor(double d) + : d(d) + { + } + Extractor(unsigned sign, unsigned exponent, unsigned long long mantissa) + : mantissa(mantissa) + , exponent(exponent) + , sign(sign) + { + } + struct { + unsigned long long mantissa : 52; + unsigned exponent : 11; + unsigned sign : 1; + }; + double d; + + bool operator==(const Extractor& other) const + { + return other.sign == sign && other.exponent == exponent && other.mantissa == mantissa; + } +}; +namespace AK { +template<> +struct Formatter : StandardFormatter { + void format(FormatBuilder& builder, const Extractor& value) + { + builder.put_literal("{"); + builder.put_u64(value.sign); + builder.put_literal(", "); + builder.put_u64(value.exponent, 16, true); + builder.put_literal(", "); + builder.put_u64(value.mantissa, 16, true); + builder.put_literal("}"); + } +}; +} + +static Extractor nextafter_translator(Extractor x, Extractor target) +{ + return Extractor(nextafter(x.d, target.d)); +} + +TEST_CASE(nextafter) +{ + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x7fe, 0xfffffffffffff), Extractor(0x0, 0x7fe, 0xfffffffffffff)), Extractor(0x0, 0x7fe, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x1, 0x0), Extractor(0x0, 0x412, 0xe848000000000)), Extractor(0x0, 0x1, 0x1)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x3ff, 0x0), Extractor(0x0, 0x412, 0xe848200000000)), Extractor(0x0, 0x3ff, 0x1)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x0, 0x0), Extractor(0x0, 0x412, 0xe848000000000)), Extractor(0x0, 0x0, 0x1)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x0, 0x0), Extractor(0x0, 0x412, 0xe848000000000)), Extractor(0x0, 0x0, 0x1)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x3ff, 0x0), Extractor(0x0, 0x412, 0xe847e00000000)), Extractor(0x1, 0x3fe, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x0, 0x1), Extractor(0x0, 0x412, 0xe848000000000)), Extractor(0x0, 0x0, 0x2)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x7fe, 0xfffffffffffff), Extractor(0x0, 0x7fe, 0xfffffffffffff)), Extractor(0x0, 0x7fe, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x412, 0xe848000000000), Extractor(0x0, 0x1, 0x0)), Extractor(0x0, 0x412, 0xe847fffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x412, 0xe848200000000), Extractor(0x0, 0x3ff, 0x0)), Extractor(0x0, 0x412, 0xe8481ffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x412, 0xe848000000000), Extractor(0x1, 0x0, 0x0)), Extractor(0x0, 0x412, 0xe847fffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x412, 0xe848000000000), Extractor(0x0, 0x0, 0x0)), Extractor(0x0, 0x412, 0xe847fffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x412, 0xe847e00000000), Extractor(0x1, 0x3ff, 0x0)), Extractor(0x0, 0x412, 0xe847dffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x412, 0xe848000000000), Extractor(0x0, 0x0, 0x1)), Extractor(0x0, 0x412, 0xe847fffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x7fe, 0xfffffffffffff), Extractor(0x0, 0x7fe, 0xfffffffffffff)), Extractor(0x0, 0x7fe, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x1, 0x0), Extractor(0x0, 0x1, 0x0)), Extractor(0x0, 0x1, 0x0)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x3ff, 0x0), Extractor(0x0, 0x3ff, 0x0)), Extractor(0x0, 0x3ff, 0x0)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x0, 0x0), Extractor(0x1, 0x0, 0x0)), Extractor(0x1, 0x0, 0x0)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x0, 0x0), Extractor(0x0, 0x0, 0x0)), Extractor(0x0, 0x0, 0x0)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x3ff, 0x0), Extractor(0x1, 0x3ff, 0x0)), Extractor(0x1, 0x3ff, 0x0)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x0, 0x1), Extractor(0x0, 0x0, 0x1)), Extractor(0x0, 0x0, 0x1)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x7fe, 0xfffffffffffff), Extractor(0x0, 0x7fe, 0xfffffffffffff)), Extractor(0x1, 0x7fe, 0xffffffffffffe)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x1, 0x0), Extractor(0x0, 0x1, 0x0)), Extractor(0x1, 0x0, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x3ff, 0x0), Extractor(0x0, 0x3ff, 0x0)), Extractor(0x1, 0x3fe, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x0, 0x0), Extractor(0x1, 0x0, 0x0)), Extractor(0x1, 0x0, 0x0)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x0, 0x0), Extractor(0x0, 0x0, 0x0)), Extractor(0x0, 0x0, 0x0)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x3ff, 0x0), Extractor(0x1, 0x3ff, 0x0)), Extractor(0x0, 0x3fe, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x0, 0x1), Extractor(0x0, 0x0, 0x1)), Extractor(0x1, 0x0, 0x0)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x7fe, 0xfffffffffffff), Extractor(0x1, 0x7fe, 0xfffffffffffff)), Extractor(0x0, 0x7fe, 0xffffffffffffe)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x1, 0x0), Extractor(0x1, 0x1, 0x0)), Extractor(0x0, 0x0, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x3ff, 0x0), Extractor(0x1, 0x3ff, 0x0)), Extractor(0x0, 0x3fe, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x0, 0x0), Extractor(0x0, 0x0, 0x0)), Extractor(0x0, 0x0, 0x0)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x0, 0x0), Extractor(0x1, 0x0, 0x0)), Extractor(0x1, 0x0, 0x0)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x3ff, 0x0), Extractor(0x0, 0x3ff, 0x0)), Extractor(0x1, 0x3fe, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x0, 0x1), Extractor(0x1, 0x0, 0x1)), Extractor(0x0, 0x0, 0x0)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x7fe, 0xfffffffffffff), Extractor(0x0, 0x7fe, 0xfffffffffffff)), Extractor(0x0, 0x7fe, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x1, 0x0), Extractor(0x1, 0x419, 0x7d78400000000)), Extractor(0x0, 0x0, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x3ff, 0x0), Extractor(0x1, 0x419, 0x7d783fc000000)), Extractor(0x0, 0x3fe, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x0, 0x0), Extractor(0x1, 0x419, 0x7d78400000000)), Extractor(0x1, 0x0, 0x1)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x0, 0x0), Extractor(0x1, 0x419, 0x7d78400000000)), Extractor(0x1, 0x0, 0x1)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x3ff, 0x0), Extractor(0x1, 0x419, 0x7d78404000000)), Extractor(0x1, 0x3ff, 0x1)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x0, 0x1), Extractor(0x1, 0x419, 0x7d78400000000)), Extractor(0x0, 0x0, 0x0)); + EXPECT_EQ(nextafter_translator(Extractor(0x0, 0x7fe, 0xfffffffffffff), Extractor(0x0, 0x7fe, 0xfffffffffffff)), Extractor(0x0, 0x7fe, 0xfffffffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x419, 0x7d78400000000), Extractor(0x0, 0x1, 0x0)), Extractor(0x1, 0x419, 0x7d783ffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x419, 0x7d783fc000000), Extractor(0x0, 0x3ff, 0x0)), Extractor(0x1, 0x419, 0x7d783fbffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x419, 0x7d78400000000), Extractor(0x1, 0x0, 0x0)), Extractor(0x1, 0x419, 0x7d783ffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x419, 0x7d78400000000), Extractor(0x0, 0x0, 0x0)), Extractor(0x1, 0x419, 0x7d783ffffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x419, 0x7d78404000000), Extractor(0x1, 0x3ff, 0x0)), Extractor(0x1, 0x419, 0x7d78403ffffff)); + EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x419, 0x7d78400000000), Extractor(0x0, 0x0, 0x1)), Extractor(0x1, 0x419, 0x7d783ffffffff)); +} + TEST_MAIN(Math)