From 859b10d8bfcca3c4a30798b4522fd88ec6c66de6 Mon Sep 17 00:00:00 2001
From: Andrew Kelley <superjoe30@gmail.com>
Date: Sun, 15 Apr 2018 15:20:54 -0400
Subject: [PATCH] std.math.ln and std.math.exp use float strict mode

closes #920
---
 std/math/exp.zig        |  5 ++++
 std/math/ln.zig         |  2 ++
 test/behavior.zig       |  1 +
 test/cases/bugs/920.zig | 60 +++++++++++++++++++++++++++++++++++++++++
 4 files changed, 68 insertions(+)
 create mode 100644 test/cases/bugs/920.zig

diff --git a/std/math/exp.zig b/std/math/exp.zig
index 4032930a43..21aa558c57 100644
--- a/std/math/exp.zig
+++ b/std/math/exp.zig
@@ -6,6 +6,7 @@
 const std = @import("../index.zig");
 const math = std.math;
 const assert = std.debug.assert;
+const builtin = @import("builtin");
 
 pub fn exp(x: var) @typeOf(x) {
     const T = @typeOf(x);
@@ -17,6 +18,8 @@ pub fn exp(x: var) @typeOf(x) {
 }
 
 fn exp32(x_: f32) f32 {
+    @setFloatMode(this, builtin.FloatMode.Strict);
+
     const half = []f32 { 0.5, -0.5 };
     const ln2hi = 6.9314575195e-1;
     const ln2lo = 1.4286067653e-6;
@@ -94,6 +97,8 @@ fn exp32(x_: f32) f32 {
 }
 
 fn exp64(x_: f64) f64 {
+    @setFloatMode(this, builtin.FloatMode.Strict);
+
     const half = []const f64 { 0.5, -0.5 };
     const ln2hi: f64 = 6.93147180369123816490e-01;
     const ln2lo: f64 = 1.90821492927058770002e-10;
diff --git a/std/math/ln.zig b/std/math/ln.zig
index c349ed7c6f..d09494b998 100644
--- a/std/math/ln.zig
+++ b/std/math/ln.zig
@@ -89,6 +89,8 @@ pub fn ln_32(x_: f32) f32 {
 }
 
 pub fn ln_64(x_: f64) f64 {
+    @setFloatMode(this, @import("builtin").FloatMode.Strict);
+
     const ln2_hi: f64 = 6.93147180369123816490e-01;
     const ln2_lo: f64 = 1.90821492927058770002e-10;
     const Lg1: f64 = 6.666666666666735130e-01;
diff --git a/test/behavior.zig b/test/behavior.zig
index de39b20dad..2c10c6d71b 100644
--- a/test/behavior.zig
+++ b/test/behavior.zig
@@ -12,6 +12,7 @@ comptime {
     _ = @import("cases/bugs/655.zig");
     _ = @import("cases/bugs/656.zig");
     _ = @import("cases/bugs/828.zig");
+    _ = @import("cases/bugs/920.zig");
     _ = @import("cases/cast.zig");
     _ = @import("cases/const_slice_child.zig");
     _ = @import("cases/coroutines.zig");
diff --git a/test/cases/bugs/920.zig b/test/cases/bugs/920.zig
new file mode 100644
index 0000000000..13c03a304f
--- /dev/null
+++ b/test/cases/bugs/920.zig
@@ -0,0 +1,60 @@
+const std = @import("std");
+const math = std.math;
+const Random = std.rand.Random;
+
+const ZigTable = struct {
+    r: f64,
+    x: [257]f64,
+    f: [257]f64,
+
+    pdf: fn(f64) f64,
+    is_symmetric: bool,
+    zero_case: fn(&Random, f64) f64,
+};
+
+fn ZigTableGen(comptime is_symmetric: bool, comptime r: f64, comptime v: f64, comptime f: fn(f64) f64,
+       comptime f_inv: fn(f64) f64, comptime zero_case: fn(&Random, f64) f64) ZigTable {
+    var tables: ZigTable = undefined;
+
+    tables.is_symmetric = is_symmetric;
+    tables.r = r;
+    tables.pdf = f;
+    tables.zero_case = zero_case;
+
+    tables.x[0] = v / f(r);
+    tables.x[1] = r;
+
+    for (tables.x[2..256]) |*entry, i| {
+        const last = tables.x[2 + i - 1];
+        *entry = f_inv(v / last + f(last));
+    }
+    tables.x[256] = 0;
+
+    for (tables.f[0..]) |*entry, i| {
+        *entry = f(tables.x[i]);
+    }
+
+    return tables;
+}
+
+const norm_r = 3.6541528853610088;
+const norm_v = 0.00492867323399;
+
+fn norm_f(x: f64) f64 { return math.exp(-x * x / 2.0); }
+fn norm_f_inv(y: f64) f64 { return math.sqrt(-2.0 * math.ln(y)); }
+fn norm_zero_case(random: &Random, u: f64) f64 { return 0.0; }
+
+const NormalDist = blk: {
+    @setEvalBranchQuota(30000);
+    break :blk ZigTableGen(true, norm_r, norm_v, norm_f, norm_f_inv, norm_zero_case);
+};
+
+test "bug 920 fixed" {
+    const NormalDist1 = blk: {
+        break :blk ZigTableGen(true, norm_r, norm_v, norm_f, norm_f_inv, norm_zero_case);
+    };
+
+    for (NormalDist1.f) |_, i| {
+        std.debug.assert(NormalDist1.f[i] == NormalDist.f[i]);
+    }
+}