diff --git a/std/complex.d b/std/complex.d
index b0780512e..8e488db41 100644
--- a/std/complex.d
+++ b/std/complex.d
@@ -832,8 +832,13 @@ Complex!T sin(T)(Complex!T z)  @safe pure nothrow @nogc
 @safe pure nothrow unittest
 {
     static import std.math;
+    import std.math : feqrel;
     assert(sin(complex(0.0)) == 0.0);
-    assert(sin(complex(2.0L, 0)) == std.math.sin(2.0L));
+    assert(sin(complex(2.0, 0)) == std.math.sin(2.0));
+    auto c1 = sin(complex(2.0L, 0));
+    auto c2 = complex(std.math.sin(2.0L), 0);
+    assert(feqrel(c1.re, c2.re) >= real.mant_dig - 1 &&
+        feqrel(c1.im, c2.im) >= real.mant_dig - 1);
 }
 
 
@@ -849,17 +854,20 @@ Complex!T cos(T)(Complex!T z)  @safe pure nothrow @nogc
 ///
 @safe pure nothrow unittest
 {
-    import std.complex;
-    import std.math;
+    static import std.math;
+    import std.math : feqrel;
     assert(cos(complex(0.0)) == 1.0);
-    assert(cos(complex(1.3L)) == std.math.cos(1.3L));
+    assert(cos(complex(1.3)) == std.math.cos(1.3));
     auto c1 = cos(complex(0, 5.2L));
-    auto c2 = cosh(5.2L);
+    auto c2 = complex(std.math.cosh(5.2L), 0.0L);
     assert(feqrel(c1.re, c2.re) >= real.mant_dig - 1 &&
         feqrel(c1.im, c2.im) >= real.mant_dig - 1);
+    auto c3 = cos(complex(1.3L));
+    auto c4 = complex(std.math.cos(1.3L), 0.0L);
+    assert(feqrel(c3.re, c4.re) >= real.mant_dig - 1 &&
+        feqrel(c3.im, c4.im) >= real.mant_dig - 1);
 }
 
-
 /**
     Params: y = A real number.
     Returns: The value of cos(y) + i sin(y).
diff --git a/std/conv.d b/std/conv.d
index eaee62f24..743d203b2 100644
--- a/std/conv.d
+++ b/std/conv.d
@@ -1629,6 +1629,8 @@ private void testIntegralToFloating(Integral, Floating)()
 
 private void testFloatingToIntegral(Floating, Integral)()
 {
+    import std.math : floatTraits, RealFormat;
+
     bool convFails(Source, Target, E)(Source src)
     {
         try
@@ -1660,18 +1662,23 @@ private void testFloatingToIntegral(Floating, Integral)()
     {
         a = -a; // -Integral.min not representable as an Integral
         assert(convFails!(Floating, Integral, ConvOverflowException)(a)
-                || Floating.sizeof <= Integral.sizeof);
+                || Floating.sizeof <= Integral.sizeof
+                || floatTraits!Floating.realFormat == RealFormat.ieeeExtended53);
     }
     a = 0.0 + Integral.min;
     assert(to!Integral(a) == Integral.min);
     --a; // no more representable as an Integral
     assert(convFails!(Floating, Integral, ConvOverflowException)(a)
-            || Floating.sizeof <= Integral.sizeof);
+            || Floating.sizeof <= Integral.sizeof
+            || floatTraits!Floating.realFormat == RealFormat.ieeeExtended53);
     a = 0.0 + Integral.max;
-    assert(to!Integral(a) == Integral.max || Floating.sizeof <= Integral.sizeof);
+    assert(to!Integral(a) == Integral.max
+            || Floating.sizeof <= Integral.sizeof
+            || floatTraits!Floating.realFormat == RealFormat.ieeeExtended53);
     ++a; // no more representable as an Integral
     assert(convFails!(Floating, Integral, ConvOverflowException)(a)
-            || Floating.sizeof <= Integral.sizeof);
+            || Floating.sizeof <= Integral.sizeof
+            || floatTraits!Floating.realFormat == RealFormat.ieeeExtended53);
     // convert a value with a fractional part
     a = 3.14;
     assert(to!Integral(a) == 3);
@@ -3016,7 +3023,9 @@ if (isInputRange!Source && isSomeChar!(ElementType!Source) && !is(Source == enum
 @system unittest
 {
     // @system because strtod is not @safe.
-    static if (real.mant_dig == 53)
+    import std.math : floatTraits, RealFormat;
+
+    static if (floatTraits!real.realFormat == RealFormat.ieeeDouble)
     {
         import core.stdc.stdlib, std.exception, std.math;
 
@@ -3099,7 +3108,8 @@ if (isInputRange!Source && isSomeChar!(ElementType!Source) && !is(Source == enum
         {
             ushort[8] value;
         }
-        else static if (floatTraits!real.realFormat == RealFormat.ieeeExtended)
+        else static if (floatTraits!real.realFormat == RealFormat.ieeeExtended ||
+                        floatTraits!real.realFormat == RealFormat.ieeeExtended53)
         {
             ushort[5] value;
         }
@@ -3122,6 +3132,8 @@ if (isInputRange!Source && isSomeChar!(ElementType!Source) && !is(Source == enum
         enum s = "0x1.FFFFFFFFFFFFFFFEp-16382";
     else static if (floatTraits!real.realFormat == RealFormat.ieeeExtended)
         enum s = "0x1.FFFFFFFFFFFFFFFEp-16382";
+    else static if (floatTraits!real.realFormat == RealFormat.ieeeExtended53)
+        enum s = "0x1.FFFFFFFFFFFFFFFEp-16382";
     else static if (floatTraits!real.realFormat == RealFormat.ieeeDouble)
         enum s = "0x1.FFFFFFFFFFFFFFFEp-1000";
     else
@@ -3141,6 +3153,8 @@ if (isInputRange!Source && isSomeChar!(ElementType!Source) && !is(Source == enum
         else
             ld1 = strtold(s.ptr, null);
     }
+    else static if (floatTraits!real.realFormat == RealFormat.ieeeExtended53)
+        ld1 = 0x1.FFFFFFFFFFFFFFFEp-16382L; // strtold rounds to 53 bits.
     else
         ld1 = strtold(s.ptr, null);
 
diff --git a/std/internal/math/gammafunction.d b/std/internal/math/gammafunction.d
index dd206911b..c9677c724 100644
--- a/std/internal/math/gammafunction.d
+++ b/std/internal/math/gammafunction.d
@@ -253,6 +253,8 @@ static if (floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
     enum real MAXGAMMA = 1755.5483429L;
 else static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
     enum real MAXGAMMA = 1755.5483429L;
+else static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended53)
+    enum real MAXGAMMA = 1755.5483429L;
 else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
     enum real MAXGAMMA = 171.6243769L;
 else
@@ -603,6 +605,11 @@ else static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
     enum real MAXLOG = 0x1.62e42fefa39ef358p+13L;  // log(real.max)
     enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min_normal*real.epsilon) = log(smallest denormal)
 }
+else static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended53)
+{
+    enum real MAXLOG = 0x1.62e42fefa39ef358p+13L;  // log(real.max)
+    enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min_normal*real.epsilon) = log(smallest denormal)
+}
 else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
 {
     enum real MAXLOG = 0x1.62e42fefa39efp+9L;  // log(real.max)
diff --git a/std/math.d b/std/math.d
index 9ac0d740b..7ea65dd5a 100644
--- a/std/math.d
+++ b/std/math.d
@@ -478,7 +478,8 @@ T floorImpl(T)(const T x) @trusted pure nothrow @nogc
         else
             int pos = 3;
     }
-    else static if (F.realFormat == RealFormat.ieeeExtended)
+    else static if (F.realFormat == RealFormat.ieeeExtended ||
+                    F.realFormat == RealFormat.ieeeExtended53)
     {
         int exp = (y.vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff;
 
@@ -525,7 +526,10 @@ T floorImpl(T)(const T x) @trusted pure nothrow @nogc
     }
     else
     {
-        exp = (T.mant_dig - 1) - exp;
+        static if (F.realFormat == RealFormat.ieeeExtended53)
+            exp = (T.mant_dig + 11 - 1) - exp; // mant_dig is really 64
+        else
+            exp = (T.mant_dig - 1) - exp;
 
         // Zero 16 bits at a time.
         while (exp >= 16)
@@ -1062,13 +1066,13 @@ Lret: {}
         real t = tan(x);
 
         //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
-        if (!isIdentical(r, t)) assert(fabs(r-t) <= .0000001);
+        assert(approxEqual(r, t));
 
         x = -x;
         r = -r;
         t = tan(x);
         //printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
-        if (!isIdentical(r, t) && !(r != r && t != t)) assert(fabs(r-t) <= .0000001);
+        assert(approxEqual(r, t));
     }
     // overflow
     assert(isNaN(tan(real.infinity)));
@@ -1133,7 +1137,7 @@ float asin(float x) @safe pure nothrow @nogc  { return asin(cast(real) x); }
 
 @system unittest
 {
-    assert(equalsDigit(asin(0.5), PI / 6, useDigits));
+    assert(asin(0.5).approxEqual(PI / 6));
 }
 
 /***************
@@ -1362,7 +1366,7 @@ float atan2(float y, float x) @safe pure nothrow @nogc
 
 @system unittest
 {
-    assert(equalsDigit(atan2(1.0L, std.math.sqrt(3.0L)), PI / 6, useDigits));
+    assert(atan2(1.0, sqrt(3.0)).approxEqual(PI / 6));
 }
 
 /***********************************
@@ -1424,7 +1428,7 @@ float sinh(float x) @safe pure nothrow @nogc  { return sinh(cast(real) x); }
 
 @system unittest
 {
-    assert(equalsDigit(sinh(1.0), (E - 1.0 / E) / 2, useDigits));
+    assert(sinh(1.0).approxEqual((E - 1.0 / E) / 2));
 }
 
 /***********************************
@@ -1774,7 +1778,8 @@ real exp(real x) @trusted pure nothrow @nogc
             enum real OF =  7.09782712893383996732E2;  // ln((1-2^-53) * 2^1024)
             enum real UF = -7.451332191019412076235E2; // ln(2^-1075)
         }
-        else static if (F.realFormat == RealFormat.ieeeExtended)
+        else static if (F.realFormat == RealFormat.ieeeExtended ||
+                        F.realFormat == RealFormat.ieeeExtended53)
         {
             // Coefficients for exp(x)
             static immutable real[3] P = [
@@ -1865,7 +1870,7 @@ float exp(float x)  @safe pure nothrow @nogc   { return exp(cast(real) x); }
 
 @system unittest
 {
-    assert(equalsDigit(exp(3.0L), E * E * E, useDigits));
+    assert(exp(3.0).feqrel(E * E * E) > 16);
 }
 
 /**
@@ -2451,7 +2456,8 @@ private real exp2Impl(real x) @nogc @trusted pure nothrow
         ctrl.rounding = FloatingPointControl.roundToNearest;
     }
 
-    static if (real.mant_dig == 113)
+    enum realFormat = floatTraits!real.realFormat;
+    static if (realFormat == RealFormat.ieeeQuadruple)
     {
         static immutable real[2][] exptestpoints =
         [ //  x               exp(x)
@@ -2470,7 +2476,8 @@ private real exp2Impl(real x) @nogc @trusted pure nothrow
             [-0x1p+30L,       0                                        ], // far underflow
         ];
     }
-    else static if (real.mant_dig == 64) // 80-bit reals
+    else static if (realFormat == RealFormat.ieeeExtended ||
+                    realFormat == RealFormat.ieeeExtended53)
     {
         static immutable real[2][] exptestpoints =
         [ //  x               exp(x)
@@ -2489,7 +2496,7 @@ private real exp2Impl(real x) @nogc @trusted pure nothrow
             [-0x1p+30L,       0                            ], // far underflow
         ];
     }
-    else static if (real.mant_dig == 53) // 64-bit reals
+    else static if (realFormat == RealFormat.ieeeDouble)
     {
         static immutable real[2][] exptestpoints =
         [ //  x,             exp(x)
@@ -2510,14 +2517,14 @@ private real exp2Impl(real x) @nogc @trusted pure nothrow
     else
         static assert(0, "No exp() tests for real type!");
 
-    const minEqualDecimalDigits = real.dig - 3;
+    const minEqualMantissaBits = real.mant_dig - 13;
     real x;
     version (IeeeFlagsSupport) IeeeFlags f;
     foreach (ref pair; exptestpoints)
     {
         version (IeeeFlagsSupport) resetIeeeFlags();
         x = exp(pair[0]);
-        assert(equalsDigit(x, pair[1], minEqualDecimalDigits));
+        assert(feqrel(x, pair[1]) >= minEqualMantissaBits);
     }
 
     // Ideally, exp(0) would not set the inexact flag.
@@ -2633,7 +2640,8 @@ if (isFloatingPoint!T)
     alias F = floatTraits!T;
 
     ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
-    static if (F.realFormat == RealFormat.ieeeExtended)
+    static if (F.realFormat == RealFormat.ieeeExtended ||
+               F.realFormat == RealFormat.ieeeExtended53)
     {
         if (ex)
         {   // If exponent is non-zero
@@ -2921,7 +2929,8 @@ if (isFloatingPoint!T)
     y.rv = x;
 
     int ex = y.vu[F.EXPPOS_SHORT] & F.EXPMASK;
-    static if (F.realFormat == RealFormat.ieeeExtended)
+    static if (F.realFormat == RealFormat.ieeeExtended ||
+               F.realFormat == RealFormat.ieeeExtended53)
     {
         if (ex)
         {
@@ -3167,6 +3176,7 @@ float ldexp(float n, int exp) @safe pure nothrow @nogc { return ldexp(cast(real)
 @safe pure nothrow @nogc unittest
 {
     static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended ||
+               floatTraits!(real).realFormat == RealFormat.ieeeExtended53 ||
                floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
     {
         assert(ldexp(1.0L, -16384) == 0x1p-16384L);
@@ -4411,12 +4421,16 @@ long lrint(real x) @trusted pure nothrow @nogc
 
             return sign ? -result : result;
         }
-        else static if (F.realFormat == RealFormat.ieeeExtended)
+        else static if (F.realFormat == RealFormat.ieeeExtended ||
+                        F.realFormat == RealFormat.ieeeExtended53)
         {
             long result;
 
             // Rounding limit when casting from real(80-bit) to ulong.
-            enum real OF = 9.22337203685477580800E18L;
+            static if (F.realFormat == RealFormat.ieeeExtended)
+                enum real OF = 9.22337203685477580800E18L;
+            else
+                enum real OF = 4.50359962737049600000E15L;
 
             ushort* vu = cast(ushort*)(&x);
             uint* vi = cast(uint*)(&x);
@@ -5654,7 +5668,8 @@ bool isSubnormal(X)(X x) @trusted pure nothrow @nogc
         return (e == 0 &&
           ((ps[MANTISSA_LSB]|(ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF)) != 0));
     }
-    else static if (F.realFormat == RealFormat.ieeeExtended)
+    else static if (F.realFormat == RealFormat.ieeeExtended ||
+                    F.realFormat == RealFormat.ieeeExtended53)
     {
         ushort* pe = cast(ushort *)&x;
         long*   ps = cast(long *)&x;
@@ -5704,7 +5719,8 @@ if (isFloatingPoint!(X))
         return ((*cast(ulong *)&x) & 0x7FFF_FFFF_FFFF_FFFF)
             == 0x7FF0_0000_0000_0000;
     }
-    else static if (F.realFormat == RealFormat.ieeeExtended)
+    else static if (F.realFormat == RealFormat.ieeeExtended ||
+                    F.realFormat == RealFormat.ieeeExtended53)
     {
         const ushort e = cast(ushort)(F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]);
         const ulong ps = *cast(ulong *)&x;
@@ -5967,7 +5983,8 @@ F sgn(F)(F x) @safe pure nothrow @nogc
 real NaN(ulong payload) @trusted pure nothrow @nogc
 {
     alias F = floatTraits!(real);
-    static if (F.realFormat == RealFormat.ieeeExtended)
+    static if (F.realFormat == RealFormat.ieeeExtended ||
+               F.realFormat == RealFormat.ieeeExtended53)
     {
         // real80 (in x86 real format, the implied bit is actually
         // not implied but a real bit which is stored in the real)
@@ -6173,11 +6190,14 @@ real nextUp(real x) @trusted pure nothrow @nogc
         }
         return x;
     }
-    else static if (F.realFormat == RealFormat.ieeeExtended)
+    else static if (F.realFormat == RealFormat.ieeeExtended ||
+                    F.realFormat == RealFormat.ieeeExtended53)
     {
         // For 80-bit reals, the "implied bit" is a nuisance...
         ushort *pe = cast(ushort *)&x;
         ulong  *ps = cast(ulong  *)&x;
+        // EPSILON is 1 for 64-bit, and 2048 for 53-bit precision reals.
+        enum ulong EPSILON = 2UL ^^ (64 - real.mant_dig);
 
         if ((pe[F.EXPPOS_SHORT] & F.EXPMASK) == F.EXPMASK)
         {
@@ -6188,7 +6208,7 @@ real nextUp(real x) @trusted pure nothrow @nogc
         if (pe[F.EXPPOS_SHORT] & 0x8000)
         {
             // Negative number -- need to decrease the significand
-            --*ps;
+            *ps -= EPSILON;
             // Need to mask with 0x7FFF... so subnormals are treated correctly.
             if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_FFFF_FFFF_FFFF)
             {
@@ -6213,7 +6233,7 @@ real nextUp(real x) @trusted pure nothrow @nogc
         {
             // Positive number -- need to increase the significand.
             // Works automatically for positive zero.
-            ++*ps;
+            *ps += EPSILON;
             if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0)
             {
                 // change in exponent
@@ -6977,6 +6997,7 @@ if (isFloatingPoint!(X))
         static assert(F.realFormat == RealFormat.ieeeSingle
                     || F.realFormat == RealFormat.ieeeDouble
                     || F.realFormat == RealFormat.ieeeExtended
+                    || F.realFormat == RealFormat.ieeeExtended53
                     || F.realFormat == RealFormat.ieeeQuadruple);
 
         if (x == y)
@@ -7116,7 +7137,8 @@ body
 
     alias F = floatTraits!(T);
     T u;
-    static if (F.realFormat == RealFormat.ieeeExtended)
+    static if (F.realFormat == RealFormat.ieeeExtended ||
+               F.realFormat == RealFormat.ieeeExtended53)
     {
         // There's slight additional complexity because they are actually
         // 79-bit reals...