diff options
| -rw-r--r-- | Userland/Libraries/LibM/math.cpp | 66 | ||||
| -rw-r--r-- | Userland/Libraries/LibM/math.h | 7 | ||||
| -rw-r--r-- | Userland/Tests/LibM/test-math.cpp | 17 |
3 files changed, 88 insertions, 2 deletions
diff --git a/Userland/Libraries/LibM/math.cpp b/Userland/Libraries/LibM/math.cpp index f0d8724f26..46e506d138 100644 --- a/Userland/Libraries/LibM/math.cpp +++ b/Userland/Libraries/LibM/math.cpp @@ -1,5 +1,6 @@ /* * Copyright (c) 2018-2020, Andreas Kling <kling@serenityos.org> + * Copyright (c) 2021, Mițca Dumitru <dumitru0mitca@gmail.com> * All rights reserved. * * Redistribution and use in source and binary forms, with or without @@ -263,6 +264,31 @@ static int internal_ilogb(FloatT x) NOEXCEPT return (int)extractor.exponent - Extractor::exponent_bias; } +template<typename FloatT> +static FloatT internal_scalbn(FloatT x, int exponent) NOEXCEPT +{ + if (x == 0 || !isfinite(x) || isnan(x) || exponent == 0) + return x; + + using Extractor = FloatExtractor<FloatT>; + Extractor extractor; + extractor.d = x; + + if (extractor.exponent != 0) { + extractor.exponent = clamp((int)extractor.exponent + exponent, 0, (int)Extractor::exponent_max); + return extractor.d; + } + + unsigned leading_mantissa_zeroes = extractor.mantissa == 0 ? 32 : __builtin_clz(extractor.mantissa); + int shift = min((int)leading_mantissa_zeroes, exponent); + exponent = max(exponent - shift, 0); + + extractor.exponent <<= shift; + extractor.exponent = exponent + 1; + + return extractor.d; +} + extern "C" { double trunc(double x) NOEXCEPT @@ -333,14 +359,20 @@ float powf(float x, float y) NOEXCEPT return (float)pow(x, y); } +// On systems where FLT_RADIX == 2, ldexp is equivalent to scalbn +long double ldexpl(long double x, int exp) NOEXCEPT +{ + return internal_scalbn(x, exp); +} + double ldexp(double x, int exp) NOEXCEPT { - return x * exp2(exp); + return internal_scalbn(x, exp); } float ldexpf(float x, int exp) NOEXCEPT { - return x * exp2f(exp); + return internal_scalbn(x, exp); } double tanh(double x) NOEXCEPT @@ -822,4 +854,34 @@ double copysign(double x, double y) ex.sign = ey.sign; return ex.d; } + +float scalbnf(float x, int exponent) NOEXCEPT +{ + return internal_scalbn(x, exponent); +} + +double scalbn(double x, int exponent) NOEXCEPT +{ + return internal_scalbn(x, exponent); +} + +long double scalbnl(long double x, int exponent) NOEXCEPT +{ + return internal_scalbn(x, exponent); +} + +float scalbnlf(float x, long exponent) NOEXCEPT +{ + return internal_scalbn(x, exponent); +} + +double scalbln(double x, long exponent) NOEXCEPT +{ + return internal_scalbn(x, exponent); +} + +long double scalblnl(long double x, long exponent) NOEXCEPT +{ + return internal_scalbn(x, exponent); +} } diff --git a/Userland/Libraries/LibM/math.h b/Userland/Libraries/LibM/math.h index 5b68470376..a0019cbdc3 100644 --- a/Userland/Libraries/LibM/math.h +++ b/Userland/Libraries/LibM/math.h @@ -163,6 +163,13 @@ double nexttoward(double, long double) NOEXCEPT; float nexttowardf(float, long double) NOEXCEPT; long double nexttowardl(long double, long double) NOEXCEPT; +float scalbnf(float, int) NOEXCEPT; +double scalbn(double, int) NOEXCEPT; +long double scalbnl(long double, int) NOEXCEPT; +float scalbnlf(float, long) NOEXCEPT; +double scalbln(double, long) NOEXCEPT; +long double scalblnl(long double, long) NOEXCEPT; + double copysign(double x, double y); __END_DECLS diff --git a/Userland/Tests/LibM/test-math.cpp b/Userland/Tests/LibM/test-math.cpp index 555e05df30..63a979d4e9 100644 --- a/Userland/Tests/LibM/test-math.cpp +++ b/Userland/Tests/LibM/test-math.cpp @@ -26,6 +26,7 @@ #include <AK/TestSuite.h> +#include <float.h> #include <math.h> TEST_CASE(trig) @@ -203,4 +204,20 @@ TEST_CASE(nextafter) EXPECT_EQ(nextafter_translator(Extractor(0x1, 0x419, 0x7d78400000000), Extractor(0x0, 0x0, 0x1)), Extractor(0x1, 0x419, 0x7d783ffffffff)); } +TEST_CASE(scalbn) +{ + EXPECT(isnan(scalbn(NAN, 3))); + EXPECT(!isfinite(scalbn(INFINITY, 5))); + EXPECT_EQ(scalbn(0, 3), 0); + EXPECT_EQ(scalbn(15.3, 0), 15.3); + + EXPECT_EQ(scalbn(0x0.0000000000008p-1022, 16), 0x0.0000000000008p-1006); + static constexpr auto biggest_subnormal = DBL_MIN - DBL_TRUE_MIN; + auto smallest_normal = scalbn(biggest_subnormal, 1); + Extractor ex(smallest_normal); + EXPECT(ex.exponent != 0); + + EXPECT_EQ(scalbn(2.0, 4), 32.0); +} + TEST_MAIN(Math) |