| ... | @@ -0,0 +1,147 @@ |
| 1 | //! drand48 functions are based off a 48-bit lcg prng: https://pubs.opengroup.org/onlinepubs/9799919799/functions/drand48.html |
| 2 | |
| 3 | const std = @import("std"); |
| 4 | const common = @import("../common.zig"); |
| 5 | const builtin = @import("builtin"); |
| 6 | const Lcg = std.Random.lcg.Wrapping(u48); |
| 7 | |
| 8 | comptime { |
| 9 | if (builtin.target.isMuslLibC() or builtin.target.isWasiLibC()) { |
| 10 | @export(&erand48, .{ .name = "erand48", .linkage = common.linkage, .visibility = common.visibility }); |
| 11 | @export(&jrand48, .{ .name = "jrand48", .linkage = common.linkage, .visibility = common.visibility }); |
| 12 | @export(&nrand48, .{ .name = "nrand48", .linkage = common.linkage, .visibility = common.visibility }); |
| 13 | @export(&drand48, .{ .name = "drand48", .linkage = common.linkage, .visibility = common.visibility }); |
| 14 | @export(&lrand48, .{ .name = "lrand48", .linkage = common.linkage, .visibility = common.visibility }); |
| 15 | @export(&mrand48, .{ .name = "mrand48", .linkage = common.linkage, .visibility = common.visibility }); |
| 16 | @export(&lcong48, .{ .name = "lcong48", .linkage = common.linkage, .visibility = common.visibility }); |
| 17 | @export(&seed48, .{ .name = "seed48", .linkage = common.linkage, .visibility = common.visibility }); |
| 18 | @export(&srand48, .{ .name = "srand48", .linkage = common.linkage, .visibility = common.visibility }); |
| 19 | } |
| 20 | } |
| 21 | |
| 22 | // NOTE: all "magic" numbers and tests are extracted and adapted from the source above |
| 23 | |
| 24 | const default_multiplier = 0x5DEECE66D; |
| 25 | const default_addend = 0xB; |
| 26 | |
| 27 | var lcg: Lcg = .init(0, default_multiplier, default_addend); |
| 28 | var seed48_xi: [3]c_ushort = undefined; |
| 29 | |
| 30 | fn erand48(xsubi: *[3]c_ushort) callconv(.c) f64 { |
| 31 | 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); |
| 32 | |
| 33 | var separate_lcg: Lcg = .init(xi, lcg.a, lcg.c); |
| 34 | const next_xi = separate_lcg.next(); |
| 35 | |
| 36 | xsubi.* = .{ @truncate(next_xi & 0xFFFF), @truncate((next_xi >> 16) & 0xFFFF), @truncate((next_xi >> 32) & 0xFFFF) }; |
| 37 | return @as(f64, next_xi) / @as(f64, std.math.maxInt(u48)); |
| 38 | } |
| 39 | |
| 40 | fn jrand48(xsubi: *[3]c_ushort) callconv(.c) c_long { |
| 41 | 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); |
| 42 | |
| 43 | var separate_lcg: Lcg = .init(xi, lcg.a, lcg.c); |
| 44 | const next_xi = separate_lcg.next(); |
| 45 | |
| 46 | xsubi.* = .{ @truncate(next_xi & 0xFFFF), @truncate((next_xi >> 16) & 0xFFFF), @truncate((next_xi >> 32) & 0xFFFF) }; |
| 47 | return @as(i32, @bitCast(@as(u32, @truncate(next_xi >> 16)))); |
| 48 | } |
| 49 | |
| 50 | fn nrand48(xsubi: *[3]c_ushort) callconv(.c) c_long { |
| 51 | 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); |
| 52 | |
| 53 | var separate_lcg: Lcg = .init(xi, lcg.a, lcg.c); |
| 54 | const next_xi = separate_lcg.next(); |
| 55 | |
| 56 | xsubi.* = .{ @truncate(next_xi & 0xFFFF), @truncate((next_xi >> 16) & 0xFFFF), @truncate((next_xi >> 32) & 0xFFFF) }; |
| 57 | return @intCast(next_xi >> 17); // a c_long is always at least 32-bits, this is never UB |
| 58 | } |
| 59 | |
| 60 | fn drand48() callconv(.c) f64 { |
| 61 | return 2e-48 * @as(f64, lcg.next()); |
| 62 | } |
| 63 | |
| 64 | fn lrand48() callconv(.c) c_long { |
| 65 | return @intCast(lcg.next() >> 17); |
| 66 | } |
| 67 | |
| 68 | fn mrand48() callconv(.c) c_long { |
| 69 | return @as(i32, @bitCast(@as(u32, @truncate(lcg.next() >> 16)))); |
| 70 | } |
| 71 | |
| 72 | // 0..3 is `Xi`, 3..6 is `a`, 6 is `c` |
| 73 | // first low 16-bits, then mid, then high. |
| 74 | fn lcong48(param: *[7]c_ushort) callconv(.c) void { |
| 75 | 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); |
| 76 | 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); |
| 77 | lcg.c = @as(u16, @truncate(param[6])); |
| 78 | } |
| 79 | |
| 80 | fn seed48(seed16v: *[3]c_ushort) callconv(.c) *[3]c_ushort { |
| 81 | seed48_xi = .{ @truncate(lcg.xi & 0xFFFF), @truncate((lcg.xi >> 16) & 0xFFFF), @truncate((lcg.xi >> 32) & 0xFFFF) }; |
| 82 | 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); |
| 83 | lcg = .init(xi, default_multiplier, default_addend); |
| 84 | return &seed48_xi; |
| 85 | } |
| 86 | |
| 87 | fn srand48(seedval: c_long) callconv(.c) void { |
| 88 | const xi = (@as(u32, @truncate(@as(c_ulong, @bitCast(seedval)))) << 16) | 0x330E; |
| 89 | lcg = .init(xi, default_multiplier, default_addend); |
| 90 | } |
| 91 | |
| 92 | test erand48 { |
| 93 | var xsubi: [3]c_ushort = .{ 37174, 64810, 11603 }; |
| 94 | |
| 95 | try std.testing.expectApproxEqAbs(0.8965, erand48(&xsubi), 0.0005); |
| 96 | try std.testing.expectEqualSlices(c_ushort, &.{ 22537, 47966, 58735 }, &xsubi); |
| 97 | |
| 98 | try std.testing.expectApproxEqAbs(0.3375, erand48(&xsubi), 0.0005); |
| 99 | try std.testing.expectEqualSlices(c_ushort, &.{ 37344, 32911, 22119 }, &xsubi); |
| 100 | |
| 101 | try std.testing.expectApproxEqAbs(0.6475, erand48(&xsubi), 0.0005); |
| 102 | try std.testing.expectEqualSlices(c_ushort, &.{ 23659, 29872, 42445 }, &xsubi); |
| 103 | |
| 104 | try std.testing.expectApproxEqAbs(0.5005, erand48(&xsubi), 0.0005); |
| 105 | try std.testing.expectEqualSlices(c_ushort, &.{ 31642, 7875, 32802 }, &xsubi); |
| 106 | |
| 107 | try std.testing.expectApproxEqAbs(0.5065, erand48(&xsubi), 0.0005); |
| 108 | try std.testing.expectEqualSlices(c_ushort, &.{ 64669, 14399, 33170 }, &xsubi); |
| 109 | } |
| 110 | |
| 111 | test jrand48 { |
| 112 | var xsubi: [3]c_ushort = .{ 25175, 11052, 45015 }; |
| 113 | |
| 114 | try std.testing.expectEqual(1699503220, jrand48(&xsubi)); |
| 115 | try std.testing.expectEqualSlices(c_ushort, &.{ 2326, 23668, 25932 }, &xsubi); |
| 116 | |
| 117 | try std.testing.expectEqual(-992276007, jrand48(&xsubi)); |
| 118 | try std.testing.expectEqualSlices(c_ushort, &.{ 41577, 4569, 50395 }, &xsubi); |
| 119 | |
| 120 | try std.testing.expectEqual(-19535776, jrand48(&xsubi)); |
| 121 | try std.testing.expectEqualSlices(c_ushort, &.{ 31936, 59488, 65237 }, &xsubi); |
| 122 | |
| 123 | try std.testing.expectEqual(79438377, jrand48(&xsubi)); |
| 124 | try std.testing.expectEqualSlices(c_ushort, &.{ 40395, 8745, 1212 }, &xsubi); |
| 125 | |
| 126 | try std.testing.expectEqual(-1258917728, jrand48(&xsubi)); |
| 127 | try std.testing.expectEqualSlices(c_ushort, &.{ 37242, 28832, 46326 }, &xsubi); |
| 128 | } |
| 129 | |
| 130 | test nrand48 { |
| 131 | var xsubi: [3]c_ushort = .{ 546, 33817, 23389 }; |
| 132 | |
| 133 | try std.testing.expectEqual(914920692, nrand48(&xsubi)); |
| 134 | try std.testing.expectEqualSlices(c_ushort, &.{ 29829, 10728, 27921 }, &xsubi); |
| 135 | |
| 136 | try std.testing.expectEqual(754104482, nrand48(&xsubi)); |
| 137 | try std.testing.expectEqualSlices(c_ushort, &.{ 6828, 28997, 23013 }, &xsubi); |
| 138 | |
| 139 | try std.testing.expectEqual(609453945, nrand48(&xsubi)); |
| 140 | try std.testing.expectEqualSlices(c_ushort, &.{ 58183, 3826, 18599 }, &xsubi); |
| 141 | |
| 142 | try std.testing.expectEqual(1878644360, nrand48(&xsubi)); |
| 143 | try std.testing.expectEqualSlices(c_ushort, &.{ 36678, 44304, 57331 }, &xsubi); |
| 144 | |
| 145 | try std.testing.expectEqual(2114923686, nrand48(&xsubi)); |
| 146 | try std.testing.expectEqualSlices(c_ushort, &.{ 58585, 22861, 64542 }, &xsubi); |
| 147 | } |