1//! drand48 functions are based off a 48-bit lcg prng: https://pubs.opengroup.org/onlinepubs/9799919799/functions/drand48.html
2
3const builtin = @import("builtin");
4
5const std = @import("std");
6const Lcg = std.Random.lcg.Wrapping(u48);
7
8const symbol = @import("../../c.zig").symbol;
9
10comptime {
11 if (builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) {
12 symbol(&erand48, "erand48");
13 symbol(&jrand48, "jrand48");
14 symbol(&nrand48, "nrand48");
15 symbol(&drand48, "drand48");
16 symbol(&lrand48, "lrand48");
17 symbol(&mrand48, "mrand48");
18 symbol(&lcong48, "lcong48");
19 symbol(&seed48, "seed48");
20 symbol(&srand48, "srand48");
21 }
22}
23
24// NOTE: all "magic" numbers and tests are extracted and adapted from the source above
25
26const default_multiplier = 0x5DEECE66D;
27const default_addend = 0xB;
28
29var lcg: Lcg = .init(0, default_multiplier, default_addend);
30var seed48_xi: [3]c_ushort = undefined;
31
32fn erand48(xsubi: *[3]c_ushort) callconv(.c) f64 {
33 const xi = @as(u48, @as(u16, @truncate(xsubi[0])) | (@as(u48, @as(u16, @truncate(xsubi[1])))) << 16) | (@as(u48, @as(u16, @truncate(xsubi[2]))) << 32);
34
35 var separate_lcg: Lcg = .init(xi, lcg.a, lcg.c);
36 const next_xi = separate_lcg.next();
37
38 xsubi.* = .{ @truncate(next_xi & 0xFFFF), @truncate((next_xi >> 16) & 0xFFFF), @truncate((next_xi >> 32) & 0xFFFF) };
39 return @as(f64, @bitCast(0x3ff0000000000000 | (@as(u64, next_xi) << 4))) - 1.0;
40}
41
42fn jrand48(xsubi: *[3]c_ushort) callconv(.c) c_long {
43 const xi = @as(u48, @as(u16, @truncate(xsubi[0])) | (@as(u48, @as(u16, @truncate(xsubi[1])))) << 16) | (@as(u48, @as(u16, @truncate(xsubi[2]))) << 32);
44
45 var separate_lcg: Lcg = .init(xi, lcg.a, lcg.c);
46 const next_xi = separate_lcg.next();
47
48 xsubi.* = .{ @truncate(next_xi & 0xFFFF), @truncate((next_xi >> 16) & 0xFFFF), @truncate((next_xi >> 32) & 0xFFFF) };
49 return @as(i32, @bitCast(@as(u32, @truncate(next_xi >> 16))));
50}
51
52fn nrand48(xsubi: *[3]c_ushort) callconv(.c) c_long {
53 const xi = @as(u48, @as(u16, @truncate(xsubi[0])) | (@as(u48, @as(u16, @truncate(xsubi[1])))) << 16) | (@as(u48, @as(u16, @truncate(xsubi[2]))) << 32);
54
55 var separate_lcg: Lcg = .init(xi, lcg.a, lcg.c);
56 const next_xi = separate_lcg.next();
57
58 xsubi.* = .{ @truncate(next_xi & 0xFFFF), @truncate((next_xi >> 16) & 0xFFFF), @truncate((next_xi >> 32) & 0xFFFF) };
59 return @intCast(next_xi >> 17); // a c_long is always at least 32-bits, this is never UB
60}
61
62fn drand48() callconv(.c) f64 {
63 return @as(f64, @bitCast(0x3ff0000000000000 | (@as(u64, lcg.next()) << 4))) - 1.0;
64}
65
66fn lrand48() callconv(.c) c_long {
67 return @intCast(lcg.next() >> 17);
68}
69
70fn mrand48() callconv(.c) c_long {
71 return @as(i32, @bitCast(@as(u32, @truncate(lcg.next() >> 16))));
72}
73
74// 0..3 is `Xi`, 3..6 is `a`, 6 is `c`
75// first low 16-bits, then mid, then high.
76fn lcong48(param: *[7]c_ushort) callconv(.c) void {
77 lcg.xi = (@as(u48, @as(u16, @truncate(param[0]))) | (@as(u48, @as(u16, @truncate(param[1])))) << 16) | (@as(u48, @as(u16, @truncate(param[2]))) << 32);
78 lcg.a = (@as(u48, @as(u16, @truncate(param[3]))) | (@as(u48, @as(u16, @truncate(param[4])))) << 16) | (@as(u48, @as(u16, @truncate(param[5]))) << 32);
79 lcg.c = @as(u16, @truncate(param[6]));
80}
81
82fn seed48(seed16v: *[3]c_ushort) callconv(.c) *[3]c_ushort {
83 seed48_xi = .{ @truncate(lcg.xi & 0xFFFF), @truncate((lcg.xi >> 16) & 0xFFFF), @truncate((lcg.xi >> 32) & 0xFFFF) };
84 const xi = (@as(u48, @as(u16, @truncate(seed16v[0]))) | (@as(u48, @as(u16, @truncate(seed16v[1])))) << 16) | (@as(u48, @as(u16, @truncate(seed16v[2]))) << 32);
85 lcg = .init(xi, default_multiplier, default_addend);
86 return &seed48_xi;
87}
88
89fn srand48(seedval: c_long) callconv(.c) void {
90 const xi = (@as(u32, @truncate(@as(c_ulong, @bitCast(seedval)))) << 16) | 0x330E;
91 lcg = .init(xi, default_multiplier, default_addend);
92}