LibM: Add scalbn and associated functions

Also implement ldexp in terms of them, and add ldexpl

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);
+}
 }

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

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)