From 70e5e942a0e03745edbea705b03f5d7522931cee Mon Sep 17 00:00:00 2001
From: Martin <noone@nowhere.com>
Date: Fri, 15 May 2015 21:42:24 +0200
Subject: [PATCH] std.math: proper support for 64-bit real in exp()

---
 std/math.d | 105 ++++++++++++++++++++++++++++++++++++++---------------
 1 file changed, 75 insertions(+), 30 deletions(-)

diff --git a/std/math.d b/std/math.d
index 4982ab1c0d1..5409589d5ae 100644
--- a/std/math.d
+++ b/std/math.d
@@ -1592,10 +1592,21 @@ real exp(real x) @trusted pure nothrow @nogc
         enum real C2 = 1.428606820309417232121458176568075500134E-6L;
 
         // Overflow and Underflow limits.
-        enum real OF =  1.1356523406294143949492E4L;
-        enum real UF = -1.1432769596155737933527E4L;
+        static if (real.mant_dig == 64) // 80-bit reals
+        {
+            enum real OF =  1.1356523406294143949492E4L; // ln((1-2^-64) * 2^16384) ≈ 0x1.62e42fefa39efp+13
+            enum real UF = -1.1398805384308300613366E4L; // ln(2^-16445) ≈ -0x1.6436716d5406ep+13
+        }
+        else static if (real.mant_dig == 53) // 64-bit reals
+        {
+            enum real OF =  0x1.62e42fefa39efp+9L; // ln((1-2^-53) * 2^1024) = 709.78271...
+            enum real UF = -0x1.74385446d71c3p+9L; // ln(2^-1074) = -744.44007...
+        }
+        else
+            static assert(0, "No over-/underflow limits for real type!");
 
         // Special cases.
+        // FIXME: set IEEE flags accordingly
         if (isNaN(x))
             return x;
         if (x > OF)
@@ -2140,39 +2151,73 @@ unittest
         ctrl.disableExceptions(FloatingPointControl.allExceptions);
     ctrl.rounding = FloatingPointControl.roundToNearest;
 
-    // @@BUG@@: Non-immutable array literals are ridiculous.
-    // Note that these are only valid for 80-bit reals: overflow will be different for 64-bit reals.
-    static const real [2][] exptestpoints =
-    [ // x,            exp(x)
-        [1.0L,           E                           ],
-        [0.5L,           0x1.A612_98E1_E069_BC97p+0L ],
-        [3.0L,           E*E*E                       ],
-        [0x1.1p13L,      0x1.29aeffefc8ec645p+12557L ], // near overflow
-        [-0x1.18p13L,    0x1.5e4bf54b4806db9p-12927L ], // near underflow
-        [-0x1.625p13L,   0x1.a6bd68a39d11f35cp-16358L],
-        [-0x1p30L,       0                           ], // underflow - subnormal
-        [-0x1.62DAFp13L, 0x1.96c53d30277021dp-16383L ],
-        [-0x1.643p13L,   0x1p-16444L                 ],
-        [-0x1.645p13L,   0                           ], // underflow to zero
-        [0x1p80L,        real.infinity               ], // far overflow
-        [real.infinity,  real.infinity               ],
-        [0x1.7p13L,      real.infinity               ]  // close overflow
-    ];
+    static if (real.mant_dig == 64) // 80-bit reals
+    {
+        static immutable real[2][] exptestpoints =
+        [ //  x               exp(x)
+            [ 1.0L,           E                            ],
+            [ 0.5L,           0x1.a61298e1e069bc97p+0L     ],
+            [ 3.0L,           E*E*E                        ],
+            [ 0x1.1p+13L,     0x1.29aeffefc8ec645p+12557L  ], // near overflow
+            [ 0x1.7p+13L,     real.infinity                ], // close overflow
+            [ 0x1p+80L,       real.infinity                ], // far overflow
+            [ real.infinity,  real.infinity                ],
+            [-0x1.18p+13L,    0x1.5e4bf54b4806db9p-12927L  ], // near underflow
+            [-0x1.625p+13L,   0x1.a6bd68a39d11f35cp-16358L ], // ditto
+            [-0x1.62dafp+13L, 0x1.96c53d30277021dp-16383L  ], // near underflow - subnormal
+            [-0x1.643p+13L,   0x1p-16444L                  ], // ditto
+            [-0x1.645p+13L,   0                            ], // close underflow
+            [-0x1p+30L,       0                            ], // far underflow
+        ];
+    }
+    else static if (real.mant_dig == 53) // 64-bit reals
+    {
+        static immutable real[2][] exptestpoints =
+        [ //  x,             exp(x)
+            [ 1.0L,          E                        ],
+            [ 0.5L,          0x1.a61298e1e069cp+0L    ],
+            [ 3.0L,          E*E*E                    ],
+            [ 0x1.6p+9L,     0x1.93bf4ec282efbp+1015L ], // near overflow
+            [ 0x1.7p+9L,     real.infinity            ], // close overflow
+            [ 0x1p+80L,      real.infinity            ], // far overflow
+            [ real.infinity, real.infinity            ],
+            [-0x1.6p+9L,     0x1.44a3824e5285fp-1016L ], // near underflow
+            [-0x1.64p+9L,    0x0.06f84920bb2d3p-1022L ], // near underflow - subnormal
+            [-0x1.743p+9L,   0x0.0000000000001p-1022L ], // ditto
+            [-0x1.8p+9L,     0                        ], // close underflow
+            [-0x1p30L,       0                        ], // far underflow
+        ];
+    }
+    else
+        static assert(0, "No exp() tests for real type!");
+
+    const minEqualMantissaBits = real.mant_dig - 2;
     real x;
     IeeeFlags f;
-    for (int i=0; i<exptestpoints.length;++i)
+    foreach (ref pair; exptestpoints)
     {
         resetIeeeFlags();
-        x = exp(exptestpoints[i][0]);
+        x = exp(pair[0]);
         f = ieeeFlags;
-        assert(x == exptestpoints[i][1]);
-        // Check the overflow bit
-        assert(f.overflow == (fabs(x) == real.infinity));
-        // Check the underflow bit
-        assert(f.underflow == (fabs(x) < real.min_normal));
-        // Invalid and div by zero shouldn't be affected.
-        assert(!f.invalid);
-        assert(!f.divByZero);
+        assert(feqrel(x, pair[1]) >= minEqualMantissaBits);
+
+        version (IeeeFlagsSupport)
+        {
+            // Check the overflow bit
+            if (x == real.infinity)
+            {
+                // don't care about the overflow bit if input was inf
+                // (e.g., the LLVM intrinsic doesn't set it on Linux x86_64)
+                assert(pair[0] == real.infinity || f.overflow);
+            }
+            else
+                assert(!f.overflow);
+            // Check the underflow bit
+            assert(f.underflow == (fabs(x) < real.min_normal));
+            // Invalid and div by zero shouldn't be affected.
+            assert(!f.invalid);
+            assert(!f.divByZero);
+        }
     }
     // Ideally, exp(0) would not set the inexact flag.
     // Unfortunately, fldl2e sets it!