diff --git a/std/format/internal/floats.d b/std/format/internal/floats.d
index 16c20655a..169249225 100644
--- a/std/format/internal/floats.d
+++ b/std/format/internal/floats.d
@@ -48,6 +48,7 @@ if (is(T == float) || is(T == double)
             import std.math.rounding : floor;
             import std.math.traits : isInfinity, isNaN;
             import std.math.exponential : log2;
+            import std.math.operations : nextUp;
 
             auto val2 = val;
 
@@ -55,6 +56,8 @@ if (is(T == float) || is(T == double)
                 exp = 32767;
             else if (abs(val2) < real.min_normal)
                 exp = 0;
+            else if (abs(val2) >= nextUp(real.max / 2))
+                exp = 32766;
             else
                 exp = cast(int) (val2.abs.log2.floor() + 16383);
 
@@ -63,22 +66,31 @@ if (is(T == float) || is(T == double)
                 // NaN or infinity
                 mnt = isNaN(val2) ? ((1L << 63) - 1) : 0;
             }
-            else if (exp > 16382 + 64) // bias + bits of ulong
-            {
-                val2 /= 2.0L ^^ (exp - (16382 + 64));
-                mnt = (cast(ulong) abs(val2)) & ((1L << 63) - 1);
-            }
             else
             {
-                auto delta = 16382 + 64 - (exp == 0 ? 1 : exp); // -1 in case of subnormals
-                if (delta > 16383)
+                if (exp > 16382 + 64) // bias + bits of ulong
+                    val2 /= 2.0L ^^ (exp - (16382 + 64));
+                else
                 {
-                    // need two steps to avoid overflow
-                    val2 *= 2.0L ^^ 16383;
-                    delta -= 16383;
+                    auto delta = 16382 + 64 - (exp == 0 ? 1 : exp); // -1 in case of subnormals
+                    if (delta > 16383)
+                    {
+                        // need two steps to avoid overflow
+                        val2 *= 2.0L ^^ 16383;
+                        delta -= 16383;
+                    }
+                    val2 *= 2.0L ^^ delta;
                 }
-                val2 *= 2.0L ^^ delta;
-                mnt = (cast(ulong) abs(val2)) & ((1L << 63) - 1);
+
+                ulong tmp = cast(ulong) abs(val2);
+                if (exp != 32767 && exp > 0 && tmp <= ulong.max / 2)
+                {
+                    // correction, due to log2(val2) being rounded up:
+                    exp--;
+                    val2 *= 2.0L;
+                    tmp = cast(ulong) abs(val2);
+                }
+                mnt = tmp & ((1L << 63) - 1);
             }
 
             double d = cast(double) val2;
diff --git a/std/format/internal/write.d b/std/format/internal/write.d
index 78afc6fe7..91a1a654c 100644
--- a/std/format/internal/write.d
+++ b/std/format/internal/write.d
@@ -918,6 +918,21 @@ if (is(FloatingPointTypeOf!T) && !is(T == enum) && !hasToString!(T, Char))
     }
 }
 
+// https://issues.dlang.org/show_bug.cgi?id=21853
+@safe pure unittest
+{
+    import std.math.exponential : log2;
+
+    // log2 is broken for x87-reals on some computers in CTFE
+    // the following test excludes these computers from the test
+    // (issue 21757)
+    enum test = cast(int) log2(3.05e2312L);
+    static if (real.mant_dig == 64 && test == 7681) // 80 bit reals
+    {
+        static assert(format!"%e"(real.max) == "1.189731e+4932");
+    }
+}
+
 // https://issues.dlang.org/show_bug.cgi?id=20536
 @safe pure unittest
 {