lib/accy/src/profiling/wos/render.zig
daab053ee43316e1809a84551d573ddd1e5bf3d2
1 const std = @import("std");
2 const gpu = @import("gpu");
3 const choir_abi = @import("choir_abi");
4 const accy = @import("accy");
5 const bench = @import("bench");
6 const sys = @import("sys");
7
8 const options_mod = @import("options.zig");
9
10 const Allocator = std.mem.Allocator;
11 const kernel = accy.kernel;
12 const pretty_json = bench.pretty.json;
13 const random = accy.kernel.library.random;
14
15 const tau: f32 = 2.0 * std.math.pi;
16
17 pub const Options = struct {
18 sdf: []const u8 = "",
19 dir: []const u8 = "",
20 res: u32 = 128,
21 image: u32 = 360,
22 frames: u32 = 140,
23 walks: u32 = 4,
24 slice: f32 = 0.02,
25 epsilon: f32 = 0.004,
26 cap: u32 = 96,
27 };
28
29 pub const oracle_sample_pixels: usize = 16;
30 pub const oracle_host_walks: usize = 4096;
31 pub const oracle_sigma: f64 = 6.0;
32 pub const oracle_slack: f64 = 1e-3;
33
34 const WalkState = struct {
35 x: kernel.Value,
36 y: kernel.Value,
37 z: kernel.Value,
38 steps: kernel.Value,
39 };
40
41 pub fn buildBunnyGraph(allocator: Allocator, options: Options) !kernel.Graph {
42 var builder = try kernel.Builder.init(allocator, kernel.Builder.Limits.standard, "accy_wos_bunny", &.{
43 kernel.dynamicBuffer(.f32),
44 kernel.dynamicBuffer(.f32),
45 kernel.dynamicBuffer(.f32),
46 kernel.dynamicBuffer(.f32),
47 kernel.scalar(.i32),
48 kernel.scalar(.i32),
49 });
50 errdefer builder.deinit();
51
52 const xs = builder.argument(0);
53 const ys = builder.argument(1);
54 const sdf = builder.argument(2);
55 const dst = builder.argument(3);
56 const seed_lo = builder.argument(4);
57 const seed_hi = builder.argument(5);
58
59 const gid = try builder.globalId(.x);
60 const walk_id = try builder.cast(gid, .i32);
61 const x0 = try builder.load(xs, gid);
62 const y0 = try builder.load(ys, gid);
63 const z0 = try builder.constantFloat(.f32, options.slice);
64 const one_counter = try builder.constantInt(.i32, 1);
65 const key_words = try random.block.philoxWords(
66 &builder,
67 random.philox_default_rounds,
68 walk_id,
69 one_counter,
70 seed_lo,
71 seed_hi,
72 );
73
74 const zero_i = try builder.constantInt(.i32, 0);
75 const one_i = try builder.constantInt(.i32, 1);
76 const cap_value = try builder.constantInt(.i32, @intCast(options.cap));
77 const epsilon_value = try builder.constantFloat(.f32, options.epsilon);
78
79 const Context = struct {
80 sdf: kernel.Value,
81 key_lo: kernel.Value,
82 key_hi: kernel.Value,
83 one_i: kernel.Value,
84 cap: kernel.Value,
85 epsilon: kernel.Value,
86 res: u32,
87 };
88
89 const walked = try builder.whileLoop(WalkState{
90 .x = x0,
91 .y = y0,
92 .z = z0,
93 .steps = zero_i,
94 }, Context{
95 .sdf = sdf,
96 .key_lo = key_words[0],
97 .key_hi = key_words[1],
98 .one_i = one_i,
99 .cap = cap_value,
100 .epsilon = epsilon_value,
101 .res = options.res,
102 }, struct {
103 fn keepGoing(inner: *kernel.Builder, walk: WalkState, ctx: Context) !kernel.Value {
104 const radius = try emitRadius(inner, ctx.sdf, walk.x, walk.y, walk.z, ctx.res);
105 const can_move = try inner.compare(.gt, radius, ctx.epsilon);
106 const below_cap = try inner.compare(.lt, walk.steps, ctx.cap);
107 return inner.and_(can_move, below_cap);
108 }
109 }.keepGoing, struct {
110 fn step(inner: *kernel.Builder, walk: WalkState, ctx: Context) !WalkState {
111 const radius = try emitRadius(inner, ctx.sdf, walk.x, walk.y, walk.z, ctx.res);
112 const words = try random.block.philoxWords(
113 inner,
114 random.philox_default_rounds,
115 walk.steps,
116 try inner.constantInt(.i32, 2),
117 ctx.key_lo,
118 ctx.key_hi,
119 );
120 const u1_value = try random.block.outputWord(inner, .f32, words[0]);
121 const u2_value = try random.block.outputWord(inner, .f32, words[1]);
122 const one = try inner.constantFloat(.f32, 1.0);
123 const zdir = try inner.sub(one, try inner.mul(u1_value, try inner.constantFloat(.f32, 2.0)));
124 const rho2 = try inner.sub(one, try inner.mul(zdir, zdir));
125 const rho = try inner.sqrt(try inner.max(rho2, try inner.constantFloat(.f32, 0.0)));
126 const angle = try inner.mul(u2_value, try inner.constantFloat(.f32, tau));
127 const dx = try inner.mul(rho, try inner.cos(angle));
128 const dy = try inner.mul(rho, try inner.sin(angle));
129 return .{
130 .x = try inner.add(walk.x, try inner.mul(radius, dx)),
131 .y = try inner.add(walk.y, try inner.mul(radius, dy)),
132 .z = try inner.add(walk.z, try inner.mul(radius, zdir)),
133 .steps = try inner.add(walk.steps, ctx.one_i),
134 };
135 }
136 }.step);
137
138 const boundary = try emitBoundaryValue(&builder, walked.x, walked.y);
139 try builder.store(boundary, dst, gid);
140 try builder.return_();
141
142 return builder.finish();
143 }
144
145 fn emitRadius(
146 builder: *kernel.Builder,
147 sdf: kernel.Value,
148 x: kernel.Value,
149 y: kernel.Value,
150 z: kernel.Value,
151 res: u32,
152 ) !kernel.Value {
153 const value = try emitTrilinear(builder, sdf, x, y, z, res);
154 const zero = try builder.constantFloat(.f32, 0.0);
155 return builder.max(try builder.sub(zero, value), zero);
156 }
157
158 fn emitTrilinear(
159 builder: *kernel.Builder,
160 sdf: kernel.Value,
161 x: kernel.Value,
162 y: kernel.Value,
163 z: kernel.Value,
164 res: u32,
165 ) !kernel.Value {
166 const resf: f32 = @floatFromInt(res - 1);
167 const scale = try builder.constantFloat(.f32, 0.5 * resf);
168 const limit = try builder.constantFloat(.f32, resf - 1.001);
169 const zero = try builder.constantFloat(.f32, 0.0);
170 const one = try builder.constantFloat(.f32, 1.0);
171
172 var cell: [3]kernel.Value = undefined;
173 var frac: [3]kernel.Value = undefined;
174 const coords = [3]kernel.Value{ x, y, z };
175 for (coords, 0..) |coord, axis| {
176 const g = try builder.min(try builder.max(try builder.mul(try builder.add(coord, one), scale), zero), limit);
177 const floored = try builder.floor(g);
178 cell[axis] = try builder.cast(floored, .i32);
179 frac[axis] = try builder.sub(g, floored);
180 }
181
182 const res_i = try builder.constantInt(.i32, @intCast(res));
183 var corners: [8]kernel.Value = undefined;
184 inline for (0..8) |corner| {
185 const ox: i64 = corner & 1;
186 const oy: i64 = (corner >> 1) & 1;
187 const oz: i64 = (corner >> 2) & 1;
188 const ix = try builder.add(cell[0], try builder.constantInt(.i32, ox));
189 const iy = try builder.add(cell[1], try builder.constantInt(.i32, oy));
190 const iz = try builder.add(cell[2], try builder.constantInt(.i32, oz));
191 const row = try builder.add(try builder.mul(iz, res_i), iy);
192 const index = try builder.add(try builder.mul(row, res_i), ix);
193 corners[corner] = try builder.load(sdf, try builder.castIndex(index));
194 }
195
196 const c00 = try lerp(builder, corners[0], corners[1], frac[0]);
197 const c10 = try lerp(builder, corners[2], corners[3], frac[0]);
198 const c01 = try lerp(builder, corners[4], corners[5], frac[0]);
199 const c11 = try lerp(builder, corners[6], corners[7], frac[0]);
200 const c0 = try lerp(builder, c00, c10, frac[1]);
201 const c1 = try lerp(builder, c01, c11, frac[1]);
202 return lerp(builder, c0, c1, frac[2]);
203 }
204
205 fn lerp(builder: *kernel.Builder, a: kernel.Value, b: kernel.Value, t: kernel.Value) !kernel.Value {
206 return builder.add(a, try builder.mul(try builder.sub(b, a), t));
207 }
208
209 fn emitBoundaryValue(builder: *kernel.Builder, x: kernel.Value, y: kernel.Value) !kernel.Value {
210 const wave = try builder.add(
211 try builder.mul(y, try builder.constantFloat(.f32, 7.0)),
212 try builder.mul(x, try builder.constantFloat(.f32, 2.5)),
213 );
214 const half = try builder.constantFloat(.f32, 0.5);
215 return builder.add(half, try builder.mul(half, try builder.sin(wave)));
216 }
217
218 pub fn hostWalk(
219 sdf: []const f32,
220 res: u32,
221 x0: f32,
222 y0: f32,
223 z0: f32,
224 epsilon: f32,
225 cap: u32,
226 rng: std.Random,
227 ) f32 {
228 var x = x0;
229 var y = y0;
230 var z = z0;
231 var steps: u32 = 0;
232 while (steps < cap) : (steps += 1) {
233 const radius = @max(-hostTrilinear(sdf, res, x, y, z), 0.0);
234 if (radius <= epsilon) break;
235 const zdir = 1.0 - 2.0 * rng.float(f32);
236 const rho = @sqrt(@max(1.0 - zdir * zdir, 0.0));
237 const angle = rng.float(f32) * tau;
238 x += radius * rho * @cos(angle);
239 y += radius * rho * @sin(angle);
240 z += radius * zdir;
241 }
242 return 0.5 + 0.5 * @sin(y * 7.0 + x * 2.5);
243 }
244
245 pub fn hostTrilinear(sdf: []const f32, res: u32, x: f32, y: f32, z: f32) f32 {
246 const resf: f32 = @floatFromInt(res - 1);
247 const coords = [3]f32{ x, y, z };
248 var cell: [3]u32 = undefined;
249 var frac: [3]f32 = undefined;
250 for (coords, 0..) |coord, axis| {
251 const g = std.math.clamp((coord + 1.0) * 0.5 * resf, 0.0, resf - 1.001);
252 const floored = @floor(g);
253 cell[axis] = @intFromFloat(floored);
254 frac[axis] = g - floored;
255 }
256 var acc: f32 = 0.0;
257 inline for (0..8) |corner| {
258 const ox: u32 = corner & 1;
259 const oy: u32 = (corner >> 1) & 1;
260 const oz: u32 = (corner >> 2) & 1;
261 const weight =
262 (if (ox == 1) frac[0] else 1.0 - frac[0]) *
263 (if (oy == 1) frac[1] else 1.0 - frac[1]) *
264 (if (oz == 1) frac[2] else 1.0 - frac[2]);
265 const index = ((cell[2] + oz) * res + (cell[1] + oy)) * res + (cell[0] + ox);
266 acc += weight * sdf[index];
267 }
268 return acc;
269 }
270
271 pub fn run(arena: Allocator, backing: Allocator, out: *std.Io.Writer, options: Options) !u8 {
272 const sdf_bytes = try sys.fs.readFileAlloc(arena, options.sdf, 64 * 1024 * 1024);
273 const expected = @as(usize, options.res) * options.res * options.res * @sizeOf(f32);
274 if (sdf_bytes.len != expected) return error.InvalidArguments;
275 const sdf = try arena.alloc(f32, expected / @sizeOf(f32));
276 @memcpy(std.mem.sliceAsBytes(sdf), sdf_bytes);
277
278 var cuda_state = gpu.cuda.State.initDevice(backing, 0) catch |err| switch (err) {
279 error.RuntimeUnavailable => {
280 try out.print("render: CUDA driver or device unavailable\n", .{});
281 try out.flush();
282 return 3;
283 },
284 else => return err,
285 };
286 defer cuda_state.deinit();
287 const handle = cuda_state.handle();
288
289 const image = options.image;
290 const pixel_count = @as(usize, image) * image;
291 var origin_x = try arena.alloc(f32, pixel_count);
292 var origin_y = try arena.alloc(f32, pixel_count);
293 var pixel_of_walk = try arena.alloc(u32, pixel_count);
294 var walk_count: usize = 0;
295 for (0..image) |iy| {
296 for (0..image) |ix| {
297 const wx = -1.0 + 2.0 * (@as(f32, @floatFromInt(ix)) + 0.5) / @as(f32, @floatFromInt(image));
298 const wy = 1.0 - 2.0 * (@as(f32, @floatFromInt(iy)) + 0.5) / @as(f32, @floatFromInt(image));
299 if (hostTrilinear(sdf, options.res, wx, wy, options.slice) < -options.epsilon) {
300 origin_x[walk_count] = wx;
301 origin_y[walk_count] = wy;
302 pixel_of_walk[walk_count] = @intCast(iy * image + ix);
303 walk_count += 1;
304 }
305 }
306 }
307 if (walk_count == 0) return error.InvalidArguments;
308 const block = options_mod.threads_per_block;
309 const padded = std.mem.alignForward(usize, walk_count, block);
310 for (walk_count..padded) |i| {
311 origin_x[i] = origin_x[walk_count - 1];
312 origin_y[i] = origin_y[walk_count - 1];
313 }
314
315 var graph = try buildBunnyGraph(backing, options);
316 defer graph.deinit();
317 var artifact = try kernel.createKernelArtifact(backing, handle, &graph, .{
318 .authored_kernel_diagnostic_id = "profiling/wos-bunny-render",
319 });
320 defer artifact.deinit();
321 const loaded = try handle.loadArtifact(&artifact);
322 defer handle.destroyObject(loaded.id);
323
324 const xs_binding = try writeDeviceBuffer(handle, .f32, std.mem.sliceAsBytes(origin_x[0..padded]), padded, .read_only);
325 defer handle.destroyObject(xs_binding.handle.id);
326 const ys_binding = try writeDeviceBuffer(handle, .f32, std.mem.sliceAsBytes(origin_y[0..padded]), padded, .read_only);
327 defer handle.destroyObject(ys_binding.handle.id);
328 const sdf_binding = try writeDeviceBuffer(handle, .f32, std.mem.sliceAsBytes(sdf), sdf.len, .read_only);
329 defer handle.destroyObject(sdf_binding.handle.id);
330 const dst_binding = try allocateDeviceBuffer(handle, .f32, padded * @sizeOf(f32), padded);
331 defer handle.destroyObject(dst_binding.handle.id);
332
333 const host_dst = try arena.alloc(f32, padded);
334 const sums = try arena.alloc(f64, pixel_count);
335 @memset(sums, 0.0);
336 const sumsqs = try arena.alloc(f64, pixel_count);
337 @memset(sumsqs, 0.0);
338 const img = try arena.alloc(f32, pixel_count);
339
340 try sys.fs.createDirPath(options.dir);
341
342 const geometry = choir_abi.LaunchGeometry{
343 .grid = .{ @intCast(padded / block), 1, 1 },
344 .threadgroup = .{ block, 1, 1 },
345 };
346
347 var launches: u32 = 0;
348 var frame: u32 = 0;
349 while (frame < options.frames) : (frame += 1) {
350 var walk: u32 = 0;
351 while (walk < options.walks) : (walk += 1) {
352 const seeds = [2]choir_abi.ScalarArgument{
353 .{ .i32 = @bitCast(launches *% 2654435761 +% srange) },
354 .{ .i32 = @bitCast(@as(u32, 0x9E3779B9)) },
355 };
356 try handle.launch(.{
357 .artifact = &artifact,
358 .loaded_artifact = loaded,
359 .buffers = &.{ xs_binding, ys_binding, sdf_binding, dst_binding },
360 .scalar_arguments = seeds[0..],
361 .geometry = geometry,
362 });
363 try handle.synchronize(.{ .scope = .device });
364 try handle.readBuffer(.{ .handle = dst_binding.handle, .bytes = std.mem.sliceAsBytes(host_dst) });
365 for (0..walk_count) |i| {
366 const value = host_dst[i];
367 if (!(value >= 0.0 and value <= 1.0)) return error.OracleMismatch;
368 sums[pixel_of_walk[i]] += value;
369 sumsqs[pixel_of_walk[i]] += @as(f64, value) * value;
370 }
371 launches += 1;
372 }
373
374 const scale = 1.0 / @as(f64, @floatFromInt(launches));
375 @memset(img, std.math.nan(f32));
376 for (0..walk_count) |i| {
377 img[pixel_of_walk[i]] = @floatCast(sums[pixel_of_walk[i]] * scale);
378 }
379 const path = try std.fmt.allocPrint(arena, "{s}/wos_bunny-f{d:0>4}.f32", .{ options.dir, frame });
380 defer arena.free(path);
381 try sys.fs.writeFile(path, std.mem.sliceAsBytes(img));
382 }
383
384 try verifySampledPixels(sdf, options, origin_x, origin_y, pixel_of_walk, walk_count, sums, sumsqs, launches);
385
386 const meta_path = try std.fmt.allocPrint(arena, "{s}/wos_bunny.json", .{options.dir});
387 defer arena.free(meta_path);
388 var meta_buffer: [512]u8 = undefined;
389 var meta = std.Io.Writer.fixed(&meta_buffer);
390 var stream = pretty_json.Writer.init(&meta, .minified);
391 const object = try stream.object();
392 try object.field("workload", "wos_bunny");
393 try object.field("kind", "wos");
394 try object.field("m", image);
395 try object.field("n", image);
396 try object.field("frames", options.frames);
397 try object.field("walks_per_frame", options.walks);
398 try object.print("slice", "{e:.3}", .{options.slice});
399 try object.print("epsilon", "{e:.3}", .{options.epsilon});
400 try object.field("cap", options.cap);
401 try object.field("inside_pixels", walk_count);
402 try object.field("layout", "row-major, nan outside");
403 try object.endLine();
404 try sys.fs.writeFile(meta_path, meta.buffered());
405
406 try out.print("render: wos_bunny frames={d} walks/frame={d} inside={d} dir={s}\n", .{
407 options.frames,
408 options.walks,
409 walk_count,
410 options.dir,
411 });
412 try out.flush();
413 return 0;
414 }
415
416 fn verifySampledPixels(
417 sdf: []const f32,
418 options: Options,
419 origin_x: []const f32,
420 origin_y: []const f32,
421 pixel_of_walk: []const u32,
422 walk_count: usize,
423 sums: []const f64,
424 sumsqs: []const f64,
425 launches: u32,
426 ) !void {
427 const samples: usize = @min(oracle_sample_pixels, walk_count);
428 const device_walks = @as(f64, @floatFromInt(launches));
429 const host_walks = @as(f64, @floatFromInt(oracle_host_walks));
430 var prng = std.Random.DefaultPrng.init(0x5eed0f5eed);
431 const rng = prng.random();
432
433 var sample: usize = 0;
434 while (sample < samples) : (sample += 1) {
435 const walk = (sample * walk_count) / samples + walk_count / (2 * samples);
436 const pixel = pixel_of_walk[walk];
437 const device_mean = sums[pixel] / device_walks;
438 const device_variance = @max(sumsqs[pixel] / device_walks - device_mean * device_mean, 0.0);
439
440 var host_sum: f64 = 0.0;
441 var host_sumsq: f64 = 0.0;
442 for (0..oracle_host_walks) |_| {
443 const value: f64 = hostWalk(
444 sdf,
445 options.res,
446 origin_x[walk],
447 origin_y[walk],
448 options.slice,
449 options.epsilon,
450 options.cap,
451 rng,
452 );
453 host_sum += value;
454 host_sumsq += value * value;
455 }
456 const host_mean = host_sum / host_walks;
457 const host_variance = @max(host_sumsq / host_walks - host_mean * host_mean, 0.0);
458
459 const bound = oracle_sigma * @sqrt(device_variance / device_walks + host_variance / host_walks) + oracle_slack;
460 if (@abs(device_mean - host_mean) > bound) return error.OracleMismatch;
461 }
462 }
463
464 const srange: u32 = 1013904223;
465
466 fn allocateDeviceBuffer(
467 handle: gpu.BackendHandle,
468 dtype: choir_abi.DType,
469 byte_size: usize,
470 element_count: usize,
471 ) !gpu.BufferBinding {
472 const device_buffer = try handle.allocateBuffer(.{
473 .byte_size = byte_size,
474 .alignment = 256,
475 .dtype = dtype,
476 .element_count = element_count,
477 });
478 return .{
479 .handle = device_buffer,
480 .access = .read_write,
481 .ownership = device_buffer.ownership,
482 .byte_size = device_buffer.byte_size,
483 };
484 }
485
486 fn writeDeviceBuffer(
487 handle: gpu.BackendHandle,
488 dtype: choir_abi.DType,
489 bytes: []const u8,
490 element_count: usize,
491 access: gpu.BufferAccess,
492 ) !gpu.BufferBinding {
493 var binding = try allocateDeviceBuffer(handle, dtype, bytes.len, element_count);
494 binding.access = access;
495 try handle.writeBuffer(.{ .handle = binding.handle, .bytes = bytes });
496 return binding;
497 }
498
499 test "hostTrilinear interpolates a linear ramp exactly" {
500 const res: u32 = 4;
501 var sdf: [64]f32 = undefined;
502 for (0..4) |z| for (0..4) |y| for (0..4) |x| {
503 sdf[(z * 4 + y) * 4 + x] = @as(f32, @floatFromInt(x));
504 };
505 const mid = hostTrilinear(sdf[0..], res, 0.0, -0.3, 0.4);
506 try std.testing.expectApproxEqAbs(@as(f32, 1.5), mid, 0.001);
507 }
508
509 test "buildBunnyGraph verifies" {
510 var graph = try buildBunnyGraph(std.testing.allocator, .{});
511 defer graph.deinit();
512 try graph.verify();
513 }
514
515 test "verifySampledPixels accepts self-consistent estimates" {
516 const res: u32 = 4;
517 const cells = @as(usize, res) * res * res;
518 const sdf = try std.testing.allocator.alloc(f32, cells);
519 defer std.testing.allocator.free(sdf);
520 @memset(sdf, 1.0);
521
522 const walk_count: usize = 30787;
523 const origin_x = try std.testing.allocator.alloc(f32, walk_count);
524 defer std.testing.allocator.free(origin_x);
525 @memset(origin_x, 0.0);
526 const origin_y = try std.testing.allocator.alloc(f32, walk_count);
527 defer std.testing.allocator.free(origin_y);
528 @memset(origin_y, 0.0);
529 const pixel_of_walk = try std.testing.allocator.alloc(u32, walk_count);
530 defer std.testing.allocator.free(pixel_of_walk);
531 @memset(pixel_of_walk, 0);
532
533 const launches: u32 = 2;
534 const boundary_at_origin: f64 = 0.5;
535 var sums = [_]f64{boundary_at_origin * @as(f64, launches)};
536 var sumsqs = [_]f64{boundary_at_origin * boundary_at_origin * @as(f64, launches)};
537
538 try verifySampledPixels(
539 sdf,
540 .{ .res = res },
541 origin_x,
542 origin_y,
543 pixel_of_walk,
544 walk_count,
545 sums[0..],
546 sumsqs[0..],
547 launches,
548 );
549 }
550
551 test "hostWalk from a sphere center estimates the spherical boundary mean" {
552 const res: u32 = 33;
553 const cells = @as(usize, res) * res * res;
554 const sdf = try std.testing.allocator.alloc(f32, cells);
555 defer std.testing.allocator.free(sdf);
556 const sphere_radius: f32 = 0.8;
557 for (0..res) |iz| for (0..res) |iy| for (0..res) |ix| {
558 const wx = -1.0 + 2.0 * @as(f32, @floatFromInt(ix)) / @as(f32, @floatFromInt(res - 1));
559 const wy = -1.0 + 2.0 * @as(f32, @floatFromInt(iy)) / @as(f32, @floatFromInt(res - 1));
560 const wz = -1.0 + 2.0 * @as(f32, @floatFromInt(iz)) / @as(f32, @floatFromInt(res - 1));
561 sdf[(iz * res + iy) * res + ix] = @sqrt(wx * wx + wy * wy + wz * wz) - sphere_radius;
562 };
563
564 var prng = std.Random.DefaultPrng.init(0xb0a7);
565 const rng = prng.random();
566 const walks = 20000;
567 var total: f64 = 0.0;
568 for (0..walks) |_| {
569 const value = hostWalk(sdf, res, 0.0, 0.0, 0.0, 0.02, 96, rng);
570 try std.testing.expect(value >= 0.0 and value <= 1.0);
571 total += value;
572 }
573 const mean = total / @as(f64, walks);
574 try std.testing.expectApproxEqAbs(@as(f64, 0.5), mean, 0.02);
575 }