d: Merge upstream dmd 56589f0f4, druntime 651389b5, phobos 1516ecad9.

D front-end changes:

    - Import latest bug fixes to mainline.

D runtime changes:

    - Import latest bug fixes to mainline.

Phobos changes:

    - Import latest bug fixes to mainline.

gcc/d/ChangeLog:

	* dmd/MERGE: Merge upstream dmd 56589f0f4.

libphobos/ChangeLog:

	* libdruntime/MERGE: Merge upstream druntime 651389b5.
	* src/MERGE: Merge upstream phobos 1516ecad9.

4 files changed, 422 insertions, 267 deletions

diff --git a/libphobos/src/MERGE b/libphobos/src/MERGE
index a4daa8419bf..744e5ad5e78 100644
--- a/libphobos/src/MERGE
+++ b/libphobos/src/MERGE
@@ -1,4 +1,4 @@
-a4a18d21c4ea7930f80309f85e38c571c5f6d4b8
+1516ecad932d88a1618163384e6f69009d125391
 
 The first line of this file holds the git revision number of the last
 merge done from the dlang/phobos repository.
diff --git a/libphobos/src/std/complex.d b/libphobos/src/std/complex.d
index 6b1541a2532..5a155387adc 100644
--- a/libphobos/src/std/complex.d
+++ b/libphobos/src/std/complex.d
@@ -1695,9 +1695,9 @@ Complex!T log(T)(Complex!T x) @safe pure nothrow @nogc
  */
 Complex!T log10(T)(Complex!T x) @safe pure nothrow @nogc
 {
-    static import std.math;
+    import std.math.constants : LN10;
 
-    return log(x) / Complex!T(std.math.log(10.0));
+    return log(x) / Complex!T(LN10);
 }
 
 ///
diff --git a/libphobos/src/std/file.d b/libphobos/src/std/file.d
index 05fad6724e0..b8b4a8ce6c6 100644
--- a/libphobos/src/std/file.d
+++ b/libphobos/src/std/file.d
@@ -176,9 +176,9 @@ class FileException : Exception
     private this(scope const(char)[] name, scope const(char)[] msg, string file, size_t line, uint errno) @safe pure
     {
         if (msg.empty)
-            super(name.idup, file, line);
+            super(name is null ? "(null)" : name.idup, file, line);
         else
-            super(text(name, ": ", msg), file, line);
+            super(text(name is null ? "(null)" : name, ": ", msg), file, line);
 
         this.errno = errno;
     }
@@ -1067,11 +1067,38 @@ private void removeImpl(scope const(char)[] name, scope const(FSChar)* namez) @t
         if (!name)
         {
             import core.stdc.string : strlen;
-            auto len = strlen(namez);
+
+            auto len = namez ? strlen(namez) : 0;
             name = namez[0 .. len];
         }
         cenforce(core.stdc.stdio.remove(namez) == 0,
-            "Failed to remove file " ~ name);
+            "Failed to remove file " ~ (name is null ? "(null)" : name));
+    }
+}
+
+@safe unittest
+{
+    import std.exception : collectExceptionMsg, assertThrown;
+
+    string filename = null; // e.g. as returned by File.tmpfile.name
+
+    version (linux)
+    {
+        // exact exception message is OS-dependent
+        auto msg = filename.remove.collectExceptionMsg!FileException;
+        assert("Failed to remove file (null): Bad address" == msg, msg);
+    }
+    else version (Windows)
+    {
+        import std.algorithm.searching : startsWith;
+
+        // don't test exact message on windows, it's language dependent
+        auto msg = filename.remove.collectExceptionMsg!FileException;
+        assert(msg.startsWith("(null):"), msg);
+    }
+    else
+    {
+        assertThrown!FileException(filename.remove);
     }
 }
 
diff --git a/libphobos/src/std/math/exponential.d b/libphobos/src/std/math/exponential.d
index daf2cecbebe..e32330fd932 100644
--- a/libphobos/src/std/math/exponential.d
+++ b/libphobos/src/std/math/exponential.d
@@ -2862,14 +2862,16 @@ float ldexp(float n, int exp)   @safe pure nothrow @nogc { return core.math.ldex
 
 private
 {
-    import std.math : floatTraits, RealFormat;
-
-    version (INLINE_YL2X) {} else
+    // Coefficients shared across log(), log2(), log10().
+    template LogCoeffs(T)
     {
-        static if (floatTraits!real.realFormat == RealFormat.ieeeQuadruple)
+        import std.math : floatTraits, RealFormat;
+
+        static if (floatTraits!T.realFormat == RealFormat.ieeeQuadruple)
         {
-            // Coefficients for log(1 + x) = x - x**2/2 + x**3 P(x)/Q(x)
-            static immutable real[13] logCoeffsP = [
+            // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
+            // Theoretical peak relative error = 5.3e-37
+            static immutable real[13] logP = [
                 1.313572404063446165910279910527789794488E4L,
                 7.771154681358524243729929227226708890930E4L,
                 2.014652742082537582487669938141683759923E5L,
@@ -2884,7 +2886,7 @@ private
                 4.998469661968096229986658302195402690910E-1L,
                 1.538612243596254322971797716843006400388E-6L
             ];
-            static immutable real[13] logCoeffsQ = [
+            static immutable real[13] logQ = [
                 3.940717212190338497730839731583397586124E4L,
                 2.626900195321832660448791748036714883242E5L,
                 7.777690340007566932935753241556479363645E5L,
@@ -2900,9 +2902,18 @@ private
                 1.0
             ];
 
-            // Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2)
+            // log2 uses the same coefficients as log.
+            alias log2P = logP;
+            alias log2Q = logQ;
+
+            // log10 uses the same coefficients as log.
+            alias log10P = logP;
+            alias log10Q = logQ;
+
+            // Coefficients for log(x) = z + z^^3 P(z^^2)/Q(z^^2)
             // where z = 2(x-1)/(x+1)
-            static immutable real[6] logCoeffsR = [
+            // Theoretical peak relative error = 1.1e-35
+            static immutable real[6] logR = [
                 1.418134209872192732479751274970992665513E5L,
                 -8.977257995689735303686582344659576526998E4L,
                 2.048819892795278657810231591630928516206E4L,
@@ -2910,7 +2921,7 @@ private
                 8.057002716646055371965756206836056074715E1L,
                 -8.828896441624934385266096344596648080902E-1L
             ];
-            static immutable real[7] logCoeffsS = [
+            static immutable real[7] logS = [
                 1.701761051846631278975701529965589676574E6L,
                 -1.332535117259762928288745111081235577029E6L,
                 4.001557694070773974936904547424676279307E5L,
@@ -2922,8 +2933,9 @@ private
         }
         else
         {
-            // Coefficients for log(1 + x) = x - x**2/2 + x**3 P(x)/Q(x)
-            static immutable real[7] logCoeffsP = [
+            // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
+            // Theoretical peak relative error = 2.32e-20
+            static immutable real[7] logP = [
                 2.0039553499201281259648E1L,
                 5.7112963590585538103336E1L,
                 6.0949667980987787057556E1L,
@@ -2932,7 +2944,7 @@ private
                 4.9854102823193375972212E-1L,
                 4.5270000862445199635215E-5L,
             ];
-            static immutable real[7] logCoeffsQ = [
+            static immutable real[7] logQ = [
                 6.0118660497603843919306E1L,
                 2.1642788614495947685003E2L,
                 3.0909872225312059774938E2L,
@@ -2942,15 +2954,42 @@ private
                 1.0000000000000000000000E0L,
             ];
 
-            // Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2)
+            // Coefficients for log(1 + x) = x - x^^2/2 + x^^3 P(x)/Q(x)
+            // Theoretical peak relative error = 6.2e-22
+            static immutable real[7] log2P = [
+                1.0747524399916215149070E2L,
+                3.4258224542413922935104E2L,
+                4.2401812743503691187826E2L,
+                2.5620629828144409632571E2L,
+                7.7671073698359539859595E1L,
+                1.0767376367209449010438E1L,
+                4.9962495940332550844739E-1L,
+            ];
+            static immutable real[8] log2Q = [
+                3.2242573199748645407652E2L,
+                1.2695660352705325274404E3L,
+                2.0307734695595183428202E3L,
+                1.6911722418503949084863E3L,
+                7.7952888181207260646090E2L,
+                1.9444210022760132894510E2L,
+                2.3479774160285863271658E1L,
+                1.0000000000000000000000E0,
+            ];
+
+            // log10 uses the same coefficients as log2.
+            alias log10P = log2P;
+            alias log10Q = log2Q;
+
+            // Coefficients for log(x) = z + z^^3 P(z^^2)/Q(z^^2)
             // where z = 2(x-1)/(x+1)
-            static immutable real[4] logCoeffsR = [
+            // Theoretical peak relative error = 6.16e-22
+            static immutable real[4] logR = [
                -3.5717684488096787370998E1L,
                 1.0777257190312272158094E1L,
                -7.1990767473014147232598E-1L,
                 1.9757429581415468984296E-3L,
             ];
-            static immutable real[4] logCoeffsS = [
+            static immutable real[4] logS = [
                -4.2861221385716144629696E2L,
                 1.9361891836232102174846E2L,
                -2.6201045551331104417768E1L,
@@ -2972,92 +3011,100 @@ private
  */
 real log(real x) @safe pure nothrow @nogc
 {
-    import std.math.constants : LN2, LOG2, SQRT1_2;
-    import std.math.traits : isInfinity, isNaN, signbit;
-    import std.math.algebraic : poly;
-
     version (INLINE_YL2X)
+    {
+        import std.math.constants : LN2;
         return core.math.yl2x(x, LN2);
+    }
     else
-    {
-        // C1 + C2 = LN2.
-        enum real C1 = 6.93145751953125E-1L;
-        enum real C2 = 1.428606820309417232121458176568075500134E-6L;
+        return logImpl(x);
+}
 
-        // Special cases.
-        if (isNaN(x))
-            return x;
-        if (isInfinity(x) && !signbit(x))
-            return x;
-        if (x == 0.0)
-            return -real.infinity;
-        if (x < 0.0)
-            return real.nan;
+///
+@safe pure nothrow @nogc unittest
+{
+    import std.math.operations : feqrel;
+    import std.math.constants : E;
 
-        // Separate mantissa from exponent.
-        // Note, frexp is used so that denormal numbers will be handled properly.
-        real y, z;
-        int exp;
+    assert(feqrel(log(E), 1) >= real.mant_dig - 1);
+}
 
-        x = frexp(x, exp);
+private T logImpl(T)(T x) @safe pure nothrow @nogc
+{
+    import std.math.constants : SQRT1_2;
+    import std.math.algebraic : poly;
+    import std.math.traits : isInfinity, isNaN, signbit;
 
-        // Logarithm using log(x) = z + z^^3 R(z) / S(z),
-        // where z = 2(x - 1)/(x + 1)
-        if ((exp > 2) || (exp < -2))
-        {
-            if (x < SQRT1_2)
-            {   // 2(2x - 1)/(2x + 1)
-                exp -= 1;
-                z = x - 0.5;
-                y = 0.5 * z + 0.5;
-            }
-            else
-            {   // 2(x - 1)/(x + 1)
-                z = x - 0.5;
-                z -= 0.5;
-                y = 0.5 * x  + 0.5;
-            }
-            x = z / y;
-            z = x * x;
-            z = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS));
-            z += exp * C2;
-            z += x;
-            z += exp * C1;
+    alias coeffs = LogCoeffs!T;
 
-            return z;
-        }
+    // C1 + C2 = LN2.
+    enum T C1 = 6.93145751953125E-1L;
+    enum T C2 = 1.428606820309417232121458176568075500134E-6L;
 
-        // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
+    // Special cases.
+    if (isNaN(x))
+        return x;
+    if (isInfinity(x) && !signbit(x))
+        return x;
+    if (x == 0.0)
+        return -T.infinity;
+    if (x < 0.0)
+        return T.nan;
+
+    // Separate mantissa from exponent.
+    // Note, frexp is used so that denormal numbers will be handled properly.
+    T y, z;
+    int exp;
+
+    x = frexp(x, exp);
+
+    // Logarithm using log(x) = z + z^^3 R(z) / S(z),
+    // where z = 2(x - 1)/(x + 1)
+    if ((exp > 2) || (exp < -2))
+    {
         if (x < SQRT1_2)
-        {
+        {   // 2(2x - 1)/(2x + 1)
             exp -= 1;
-            x = 2.0 * x - 1.0;
+            z = x - 0.5;
+            y = 0.5 * z + 0.5;
         }
         else
-        {
-            x = x - 1.0;
+        {   // 2(x - 1)/(x + 1)
+            z = x - 0.5;
+            z -= 0.5;
+            y = 0.5 * x  + 0.5;
         }
+        x = z / y;
         z = x * x;
-        y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ));
-        y += exp * C2;
-        z = y - 0.5 * z;
-
-        // Note, the sum of above terms does not exceed x/4,
-        // so it contributes at most about 1/4 lsb to the error.
+        z = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS));
+        z += exp * C2;
         z += x;
         z += exp * C1;
 
         return z;
     }
-}
 
-///
-@safe pure nothrow @nogc unittest
-{
-    import std.math.operations : feqrel;
-    import std.math.constants : E;
+    // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
+    if (x < SQRT1_2)
+    {
+        exp -= 1;
+        x = 2.0 * x - 1.0;
+    }
+    else
+    {
+        x = x - 1.0;
+    }
+    z = x * x;
+    y = x * (z * poly(x, coeffs.logP) / poly(x, coeffs.logQ));
+    y += exp * C2;
+    z = y - 0.5 * z;
 
-    assert(feqrel(log(E), 1) >= real.mant_dig - 1);
+    // Note, the sum of above terms does not exceed x/4,
+    // so it contributes at most about 1/4 lsb to the error.
+    z += x;
+    z += exp * C1;
+
+    return z;
 }
 
 /**************************************
@@ -3072,95 +3119,103 @@ real log(real x) @safe pure nothrow @nogc
  */
 real log10(real x) @safe pure nothrow @nogc
 {
-    import std.math.constants : LOG2, LN2, SQRT1_2;
-    import std.math.algebraic : poly;
-    import std.math.traits : isNaN, isInfinity, signbit;
-
     version (INLINE_YL2X)
+    {
+        import std.math.constants : LOG2;
         return core.math.yl2x(x, LOG2);
+    }
     else
-    {
-        // log10(2) split into two parts.
-        enum real L102A =  0.3125L;
-        enum real L102B = -1.14700043360188047862611052755069732318101185E-2L;
+        return log10Impl(x);
+}
 
-        // log10(e) split into two parts.
-        enum real L10EA =  0.5L;
-        enum real L10EB = -6.570551809674817234887108108339491770560299E-2L;
+///
+@safe pure nothrow @nogc unittest
+{
+    import std.math.algebraic : fabs;
 
-        // Special cases are the same as for log.
-        if (isNaN(x))
-            return x;
-        if (isInfinity(x) && !signbit(x))
-            return x;
-        if (x == 0.0)
-            return -real.infinity;
-        if (x < 0.0)
-            return real.nan;
+    assert(fabs(log10(1000) - 3) < .000001);
+}
 
-        // Separate mantissa from exponent.
-        // Note, frexp is used so that denormal numbers will be handled properly.
-        real y, z;
-        int exp;
+private T log10Impl(T)(T x) @safe pure nothrow @nogc
+{
+    import std.math.constants : SQRT1_2;
+    import std.math.algebraic : poly;
+    import std.math.traits : isNaN, isInfinity, signbit;
 
-        x = frexp(x, exp);
+    alias coeffs = LogCoeffs!T;
 
-        // Logarithm using log(x) = z + z^^3 R(z) / S(z),
-        // where z = 2(x - 1)/(x + 1)
-        if ((exp > 2) || (exp < -2))
-        {
-            if (x < SQRT1_2)
-            {   // 2(2x - 1)/(2x + 1)
-                exp -= 1;
-                z = x - 0.5;
-                y = 0.5 * z + 0.5;
-            }
-            else
-            {   // 2(x - 1)/(x + 1)
-                z = x - 0.5;
-                z -= 0.5;
-                y = 0.5 * x  + 0.5;
-            }
-            x = z / y;
-            z = x * x;
-            y = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS));
-            goto Ldone;
-        }
+    // log10(2) split into two parts.
+    enum T L102A =  0.3125L;
+    enum T L102B = -1.14700043360188047862611052755069732318101185E-2L;
+
+    // log10(e) split into two parts.
+    enum T L10EA =  0.5L;
+    enum T L10EB = -6.570551809674817234887108108339491770560299E-2L;
+
+    // Special cases are the same as for log.
+    if (isNaN(x))
+        return x;
+    if (isInfinity(x) && !signbit(x))
+        return x;
+    if (x == 0.0)
+        return -T.infinity;
+    if (x < 0.0)
+        return T.nan;
+
+    // Separate mantissa from exponent.
+    // Note, frexp is used so that denormal numbers will be handled properly.
+    T y, z;
+    int exp;
+
+    x = frexp(x, exp);
 
-        // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
+    // Logarithm using log(x) = z + z^^3 R(z) / S(z),
+    // where z = 2(x - 1)/(x + 1)
+    if ((exp > 2) || (exp < -2))
+    {
         if (x < SQRT1_2)
-        {
+        {   // 2(2x - 1)/(2x + 1)
             exp -= 1;
-            x = 2.0 * x - 1.0;
+            z = x - 0.5;
+            y = 0.5 * z + 0.5;
         }
         else
-            x = x - 1.0;
-
+        {   // 2(x - 1)/(x + 1)
+            z = x - 0.5;
+            z -= 0.5;
+            y = 0.5 * x  + 0.5;
+        }
+        x = z / y;
         z = x * x;
-        y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ));
-        y = y - 0.5 * z;
-
-        // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
-        // This sequence of operations is critical and it may be horribly
-        // defeated by some compiler optimizers.
-    Ldone:
-        z = y * L10EB;
-        z += x * L10EB;
-        z += exp * L102B;
-        z += y * L10EA;
-        z += x * L10EA;
-        z += exp * L102A;
+        y = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS));
+        goto Ldone;
+    }
 
-        return z;
+    // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
+    if (x < SQRT1_2)
+    {
+        exp -= 1;
+        x = 2.0 * x - 1.0;
     }
-}
+    else
+        x = x - 1.0;
 
-///
-@safe pure nothrow @nogc unittest
-{
-    import std.math.algebraic : fabs;
+    z = x * x;
+    y = x * (z * poly(x, coeffs.log10P) / poly(x, coeffs.log10Q));
+    y = y - 0.5 * z;
 
-    assert(fabs(log10(1000) - 3) < .000001);
+    // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
+    // This sequence of operations is critical and it may be horribly
+    // defeated by some compiler optimizers.
+Ldone:
+    z = y * L10EB;
+    z += x * L10EB;
+    z += exp * L102B;
+    z += y * L10EA;
+    z += x * L10EA;
+    z += exp * L102A;
+
+    return z;
 }
 
 /**
@@ -3179,29 +3234,15 @@ real log10(real x) @safe pure nothrow @nogc
  */
 real log1p(real x) @safe pure nothrow @nogc
 {
-    import std.math.traits : isNaN, isInfinity, signbit;
-    import std.math.constants : LN2;
-
     version (INLINE_YL2X)
     {
         // On x87, yl2xp1 is valid if and only if -0.5 <= lg(x) <= 0.5,
         //    ie if -0.29 <= x <= 0.414
+        import std.math.constants : LN2;
         return (core.math.fabs(x) <= 0.25)  ? core.math.yl2xp1(x, LN2) : core.math.yl2x(x+1, LN2);
     }
     else
-    {
-        // Special cases.
-        if (isNaN(x) || x == 0.0)
-            return x;
-        if (isInfinity(x) && !signbit(x))
-            return x;
-        if (x == -1.0)
-            return -real.infinity;
-        if (x < -1.0)
-            return real.nan;
-
-        return log(x + 1.0);
-    }
+        return log1pImpl(x);
 }
 
 ///
@@ -3220,6 +3261,23 @@ real log1p(real x) @safe pure nothrow @nogc
     assert(log1p(real.infinity) == real.infinity);
 }
 
+private T log1pImpl(T)(T x) @safe pure nothrow @nogc
+{
+    import std.math.traits : isNaN, isInfinity, signbit;
+
+    // Special cases.
+    if (isNaN(x) || x == 0.0)
+        return x;
+    if (isInfinity(x) && !signbit(x))
+        return x;
+    if (x == -1.0)
+        return -T.infinity;
+    if (x < -1.0)
+        return T.nan;
+
+    return logImpl(x + 1.0);
+}
+
 /***************************************
  * Calculates the base-2 logarithm of x:
  * $(SUB log, 2)x
@@ -3233,78 +3291,10 @@ real log1p(real x) @safe pure nothrow @nogc
  */
 real log2(real x) @safe pure nothrow @nogc
 {
-    import std.math.traits : isNaN, isInfinity, signbit;
-    import std.math.constants : SQRT1_2, LOG2E;
-    import std.math.algebraic : poly;
-
     version (INLINE_YL2X)
         return core.math.yl2x(x, 1.0L);
     else
-    {
-        // Special cases are the same as for log.
-        if (isNaN(x))
-            return x;
-        if (isInfinity(x) && !signbit(x))
-            return x;
-        if (x == 0.0)
-            return -real.infinity;
-        if (x < 0.0)
-            return real.nan;
-
-        // Separate mantissa from exponent.
-        // Note, frexp is used so that denormal numbers will be handled properly.
-        real y, z;
-        int exp;
-
-        x = frexp(x, exp);
-
-        // Logarithm using log(x) = z + z^^3 R(z) / S(z),
-        // where z = 2(x - 1)/(x + 1)
-        if ((exp > 2) || (exp < -2))
-        {
-            if (x < SQRT1_2)
-            {   // 2(2x - 1)/(2x + 1)
-                exp -= 1;
-                z = x - 0.5;
-                y = 0.5 * z + 0.5;
-            }
-            else
-            {   // 2(x - 1)/(x + 1)
-                z = x - 0.5;
-                z -= 0.5;
-                y = 0.5 * x  + 0.5;
-            }
-            x = z / y;
-            z = x * x;
-            y = x * (z * poly(z, logCoeffsR) / poly(z, logCoeffsS));
-            goto Ldone;
-        }
-
-        // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
-        if (x < SQRT1_2)
-        {
-            exp -= 1;
-            x = 2.0 * x - 1.0;
-        }
-        else
-            x = x - 1.0;
-
-        z = x * x;
-        y = x * (z * poly(x, logCoeffsP) / poly(x, logCoeffsQ));
-        y = y - 0.5 * z;
-
-        // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
-        // This sequence of operations is critical and it may be horribly
-        // defeated by some compiler optimizers.
-    Ldone:
-        z = y * (LOG2E - 1.0);
-        z += x * (LOG2E - 1.0);
-        z += y;
-        z += x;
-        z += exp;
-
-        return z;
-    }
+        return log2Impl(x);
 }
 
 ///
@@ -3323,6 +3313,79 @@ real log2(real x) @safe pure nothrow @nogc
     assert(isClose(log2(1024.0L), 10, 1e-18));
 }
 
+private T log2Impl(T)(T x) @safe pure nothrow @nogc
+{
+    import std.math.traits : isNaN, isInfinity, signbit;
+    import std.math.constants : SQRT1_2, LOG2E;
+    import std.math.algebraic : poly;
+
+    alias coeffs = LogCoeffs!T;
+
+    // Special cases are the same as for log.
+    if (isNaN(x))
+        return x;
+    if (isInfinity(x) && !signbit(x))
+        return x;
+    if (x == 0.0)
+        return -T.infinity;
+    if (x < 0.0)
+        return T.nan;
+
+    // Separate mantissa from exponent.
+    // Note, frexp is used so that denormal numbers will be handled properly.
+    T y, z;
+    int exp;
+
+    x = frexp(x, exp);
+
+    // Logarithm using log(x) = z + z^^3 R(z) / S(z),
+    // where z = 2(x - 1)/(x + 1)
+    if ((exp > 2) || (exp < -2))
+    {
+        if (x < SQRT1_2)
+        {   // 2(2x - 1)/(2x + 1)
+            exp -= 1;
+            z = x - 0.5;
+            y = 0.5 * z + 0.5;
+        }
+        else
+        {   // 2(x - 1)/(x + 1)
+            z = x - 0.5;
+            z -= 0.5;
+            y = 0.5 * x  + 0.5;
+        }
+        x = z / y;
+        z = x * x;
+        y = x * (z * poly(z, coeffs.logR) / poly(z, coeffs.logS));
+        goto Ldone;
+    }
+
+    // Logarithm using log(1 + x) = x - .5x^^2 + x^^3 P(x) / Q(x)
+    if (x < SQRT1_2)
+    {
+        exp -= 1;
+        x = 2.0 * x - 1.0;
+    }
+    else
+        x = x - 1.0;
+
+    z = x * x;
+    y = x * (z * poly(x, coeffs.log2P) / poly(x, coeffs.log2Q));
+    y = y - 0.5 * z;
+
+    // Multiply log of fraction by log10(e) and base 2 exponent by log10(2).
+    // This sequence of operations is critical and it may be horribly
+    // defeated by some compiler optimizers.
+Ldone:
+    z = y * (LOG2E - 1.0);
+    z += x * (LOG2E - 1.0);
+    z += y;
+    z += x;
+    z += exp;
+
+    return z;
+}
+
 /*****************************************
  * Extracts the exponent of x as a signed integral value.
  *
@@ -3337,35 +3400,23 @@ real log2(real x) @safe pure nothrow @nogc
  *      $(TR $(TD $(PLUSMN)0.0)      $(TD -$(INFIN)) $(TD yes) )
  *      )
  */
-real logb(real x) @trusted nothrow @nogc
+pragma(inline, true)
+real logb(real x) @trusted pure nothrow @nogc
 {
     version (InlineAsm_X87_MSVC)
-    {
-        version (X86_64)
-        {
-            asm pure nothrow @nogc
-            {
-                naked                       ;
-                fld     real ptr [RCX]      ;
-                fxtract                     ;
-                fstp    ST(0)               ;
-                ret                         ;
-            }
-        }
-        else
-        {
-            asm pure nothrow @nogc
-            {
-                fld     x                   ;
-                fxtract                     ;
-                fstp    ST(0)               ;
-            }
-        }
-    }
+        return logbAsm(x);
     else
-        return core.stdc.math.logbl(x);
+        return logbImpl(x);
 }
 
+/// ditto
+pragma(inline, true)
+double logb(double x) @trusted pure nothrow @nogc { return logbImpl(x); }
+
+/// ditto
+pragma(inline, true)
+float logb(float x) @trusted pure nothrow @nogc { return logbImpl(x); }
+
 ///
 @safe @nogc nothrow unittest
 {
@@ -3377,6 +3428,83 @@ real logb(real x) @trusted nothrow @nogc
     assert(logb(-real.infinity) == real.infinity);
 }
 
+@safe @nogc nothrow unittest
+{
+    import std.meta : AliasSeq;
+    import std.typecons : Tuple;
+    import std.math.traits : isNaN;
+    static foreach (F; AliasSeq!(float, double, real))
+    {{
+        alias T = Tuple!(F, F);
+        T[17] vals =   // x, logb(x)
+        [
+            T(1.0          , 0          ),
+            T(100.0        , 6          ),
+            T(0.0          , -F.infinity),
+            T(-0.0         , -F.infinity),
+            T(1024         , 10         ),
+            T(-2000        , 10         ),
+            T(0x0.1p-127   , -131       ),
+            T(0x0.01p-127  , -135       ),
+            T(0x0.011p-127 , -135       ),
+            T(F.nan        , F.nan      ),
+            T(-F.nan       , F.nan      ),
+            T(F.infinity   , F.infinity ),
+            T(-F.infinity  , F.infinity ),
+            T(F.min_normal , F.min_exp-1),
+            T(-F.min_normal, F.min_exp-1),
+            T(F.max        , F.max_exp-1),
+            T(-F.max       , F.max_exp-1),
+        ];
+
+        foreach (elem; vals)
+        {
+            if (isNaN(elem[1]))
+                assert(isNaN(logb(elem[1])));
+            else
+                assert(logb(elem[0]) == elem[1]);
+        }
+    }}
+}
+
+version (InlineAsm_X87_MSVC)
+private T logbAsm(T)(T x) @trusted pure nothrow @nogc
+{
+    version (X86_64)
+    {
+        asm pure nothrow @nogc
+        {
+            naked                       ;
+            fld     real ptr [RCX]      ;
+            fxtract                     ;
+            fstp    ST(0)               ;
+            ret                         ;
+        }
+    }
+    else
+    {
+        asm pure nothrow @nogc
+        {
+            fld     x                   ;
+            fxtract                     ;
+            fstp    ST(0)               ;
+        }
+    }
+}
+
+private T logbImpl(T)(T x) @trusted pure nothrow @nogc
+{
+    import std.math.traits : isFinite;
+
+    // Handle special cases.
+    if (!isFinite(x))
+        return x * x;
+    if (x == 0)
+        return -1 / (x * x);
+
+    return ilogb(x);
+}
+
 /*************************************
  * Efficiently calculates x * 2$(SUPERSCRIPT n).
  *