1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
|
const builtin = @import("builtin");
const std = @import("std");
const Random = std.rand.Random;
const zeroCaseFn = switch (builtin.zig_backend) {
.stage1 => fn (*Random, f64) f64,
else => *const fn (*Random, f64) f64,
};
const pdfFn = switch (builtin.zig_backend) {
.stage1 => fn (f64) f64,
else => *const fn (f64) f64,
};
const ZigTable = struct {
r: f64,
x: [257]f64,
f: [257]f64,
pdf: pdfFn,
is_symmetric: bool,
zero_case: zeroCaseFn,
};
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 @exp(-x * x / 2.0);
}
fn norm_f_inv(y: f64) f64 {
return @sqrt(-2.0 * @log(y));
}
fn norm_zero_case(random: *Random, u: f64) f64 {
_ = random;
_ = u;
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" {
if (builtin.zig_backend == .stage2_wasm) return error.SkipZigTest; // TODO
if (builtin.zig_backend == .stage2_c) return error.SkipZigTest; // TODO
if (builtin.zig_backend == .stage2_x86_64) return error.SkipZigTest; // TODO
if (builtin.zig_backend == .stage2_aarch64) return error.SkipZigTest; // TODO
if (builtin.zig_backend == .stage2_arm) return error.SkipZigTest; // TODO
const NormalDist1 = blk: {
break :blk ZigTableGen(true, norm_r, norm_v, norm_f, norm_f_inv, norm_zero_case);
};
for (NormalDist1.f) |_, i| {
// Here we use `expectApproxEqAbs` instead of `expectEqual` to account for the small
// differences in math functions of different libcs. For example, if the compiler
// links against glibc, but the target is musl libc, then these values might be
// slightly different.
// Arguably, this is a bug in the compiler because comptime should emulate the target,
// including rounding errors in libc math functions. However that behavior is not
// what this particular test is intended to cover.
try std.testing.expectApproxEqAbs(NormalDist1.f[i], NormalDist.f[i], @sqrt(std.math.floatEps(f64)));
}
}
|
