authorgravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2018-04-15 15:20:54-04:00
committergravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2018-04-15 15:22:27-04:00
log859b10d8bfcca3c4a30798b4522fd88ec6c66de6
treea62caedfaf8bf24e5eac24cc3ef630f01b536941
parenta8d794215e5031d60c5bd761a2478f382972bde6

std.math.ln and std.math.exp use float strict mode

closes #920

4 files changed, 68 insertions(+), 0 deletions(-)

std/math/exp.zig+5
...@@ -6,6 +6,7 @@...@@ -6,6 +6,7 @@
6const std = @import("../index.zig");6const std = @import("../index.zig");
7const math = std.math;7const math = std.math;
8const assert = std.debug.assert;8const assert = std.debug.assert;
9const builtin = @import("builtin");
910
10pub fn exp(x: var) @typeOf(x) {11pub fn exp(x: var) @typeOf(x) {
11 const T = @typeOf(x);12 const T = @typeOf(x);
...@@ -17,6 +18,8 @@ pub fn exp(x: var) @typeOf(x) {...@@ -17,6 +18,8 @@ pub fn exp(x: var) @typeOf(x) {
17}18}
1819
19fn exp32(x_: f32) f32 {20fn exp32(x_: f32) f32 {
21 @setFloatMode(this, builtin.FloatMode.Strict);
22
20 const half = []f32 { 0.5, -0.5 };23 const half = []f32 { 0.5, -0.5 };
21 const ln2hi = 6.9314575195e-1;24 const ln2hi = 6.9314575195e-1;
22 const ln2lo = 1.4286067653e-6;25 const ln2lo = 1.4286067653e-6;
...@@ -94,6 +97,8 @@ fn exp32(x_: f32) f32 {...@@ -94,6 +97,8 @@ fn exp32(x_: f32) f32 {
94}97}
9598
96fn exp64(x_: f64) f64 {99fn exp64(x_: f64) f64 {
100 @setFloatMode(this, builtin.FloatMode.Strict);
101
97 const half = []const f64 { 0.5, -0.5 };102 const half = []const f64 { 0.5, -0.5 };
98 const ln2hi: f64 = 6.93147180369123816490e-01;103 const ln2hi: f64 = 6.93147180369123816490e-01;
99 const ln2lo: f64 = 1.90821492927058770002e-10;104 const ln2lo: f64 = 1.90821492927058770002e-10;
std/math/ln.zig+2
...@@ -89,6 +89,8 @@ pub fn ln_32(x_: f32) f32 {...@@ -89,6 +89,8 @@ pub fn ln_32(x_: f32) f32 {
89}89}
9090
91pub fn ln_64(x_: f64) f64 {91pub fn ln_64(x_: f64) f64 {
92 @setFloatMode(this, @import("builtin").FloatMode.Strict);
93
92 const ln2_hi: f64 = 6.93147180369123816490e-01;94 const ln2_hi: f64 = 6.93147180369123816490e-01;
93 const ln2_lo: f64 = 1.90821492927058770002e-10;95 const ln2_lo: f64 = 1.90821492927058770002e-10;
94 const Lg1: f64 = 6.666666666666735130e-01;96 const Lg1: f64 = 6.666666666666735130e-01;
test/behavior.zig+1
...@@ -12,6 +12,7 @@ comptime {...@@ -12,6 +12,7 @@ comptime {
12 _ = @import("cases/bugs/655.zig");12 _ = @import("cases/bugs/655.zig");
13 _ = @import("cases/bugs/656.zig");13 _ = @import("cases/bugs/656.zig");
14 _ = @import("cases/bugs/828.zig");14 _ = @import("cases/bugs/828.zig");
15 _ = @import("cases/bugs/920.zig");
15 _ = @import("cases/cast.zig");16 _ = @import("cases/cast.zig");
16 _ = @import("cases/const_slice_child.zig");17 _ = @import("cases/const_slice_child.zig");
17 _ = @import("cases/coroutines.zig");18 _ = @import("cases/coroutines.zig");
test/cases/bugs/920.zig created+60
...@@ -0,0 +1,60 @@
1const std = @import("std");
2const math = std.math;
3const Random = std.rand.Random;
4
5const ZigTable = struct {
6 r: f64,
7 x: [257]f64,
8 f: [257]f64,
9
10 pdf: fn(f64) f64,
11 is_symmetric: bool,
12 zero_case: fn(&Random, f64) f64,
13};
14
15fn ZigTableGen(comptime is_symmetric: bool, comptime r: f64, comptime v: f64, comptime f: fn(f64) f64,
16 comptime f_inv: fn(f64) f64, comptime zero_case: fn(&Random, f64) f64) ZigTable {
17 var tables: ZigTable = undefined;
18
19 tables.is_symmetric = is_symmetric;
20 tables.r = r;
21 tables.pdf = f;
22 tables.zero_case = zero_case;
23
24 tables.x[0] = v / f(r);
25 tables.x[1] = r;
26
27 for (tables.x[2..256]) |*entry, i| {
28 const last = tables.x[2 + i - 1];
29 *entry = f_inv(v / last + f(last));
30 }
31 tables.x[256] = 0;
32
33 for (tables.f[0..]) |*entry, i| {
34 *entry = f(tables.x[i]);
35 }
36
37 return tables;
38}
39
40const norm_r = 3.6541528853610088;
41const norm_v = 0.00492867323399;
42
43fn norm_f(x: f64) f64 { return math.exp(-x * x / 2.0); }
44fn norm_f_inv(y: f64) f64 { return math.sqrt(-2.0 * math.ln(y)); }
45fn norm_zero_case(random: &Random, u: f64) f64 { return 0.0; }
46
47const NormalDist = blk: {
48 @setEvalBranchQuota(30000);
49 break :blk ZigTableGen(true, norm_r, norm_v, norm_f, norm_f_inv, norm_zero_case);
50};
51
52test "bug 920 fixed" {
53 const NormalDist1 = blk: {
54 break :blk ZigTableGen(true, norm_r, norm_v, norm_f, norm_f_inv, norm_zero_case);
55 };
56
57 for (NormalDist1.f) |_, i| {
58 std.debug.assert(NormalDist1.f[i] == NormalDist.f[i]);
59 }
60}