A draw as a uniform: the top 24 bits, in [0, 1)¶

A draw as a uniform in [0, 1) (the host's seeded error's uniform): the top 24 bits of the 32-bit draw over 2^24, exact in binary32, so a kernel and the host's seeded error agree on it to the bit.

Applied to lusine.error.draw(seed, step, key, index), it is the host's seeded error's uniform(seed, step, key, index). Checked against a reference in numpy, which has its own tests.

Source

tools/lusine/error.py, uniform.

Proof

reference, limit, curve.

uniform.trait

namespace lusine.error;

float uniform(bits drawn)
{
  return integer(drawn >> 8) * (1.0 / 16777216.0);
}

(drawn >> 8) / 2²⁴, the top 24 bits of a 32-bit draw, exact in binary32, which is the host's seeded error's uniform. Of lusine.error.draw(seed, step, key, index) it is the host's uniform(seed, step, key, index) to the bit.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import error, trait_run
M32 = 0xFFFFFFFF
EDGES = np.array([0, 1, 2, 0x7FFFFFFF, 0x80000000, 0xFFFFFFFE, M32], dtype=np.uint32)
rng = np.random.default_rng(364)
def random_bits(n):
    return np.concatenate([EDGES, rng.integers(0, 2 ** 32, n, dtype=np.uint64).astype(np.uint32)])
def lanes(law, **inputs):
    """A kernel's law over many lanes, two thousand items a program (the host's queue holds a step's stores)."""
    n = max(np.size(v) for v in inputs.values())
    parts = [law(**{k: (v[i:i + 2000] if np.ndim(v) else v) for k, v in inputs.items()}) for i in range(0, n, 2000)]
    return np.concatenate(parts)
def planted(text, n, inputs, outputs):
    """A planted program's outputs: `text` compiled by trait and run on n items, each input loaded, not a host field."""
    run = trait_run.Program(['error/finalize.trait'], text, {'probe.items': [n]})
    for field, values in inputs.items():
        run.store('probe.items', field, values)
    return [o[0] for o in run.steps('probe.main', 1, [('probe.items', f) for f in outputs])]

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np import matplotlib.pyplot as plt from lusine import error, trait_run M32 = 0xFFFFFFFF EDGES = np.array([0, 1, 2, 0x7FFFFFFF, 0x80000000, 0xFFFFFFFE, M32], dtype=np.uint32) rng = np.random.default_rng(364) def random_bits(n): return np.concatenate([EDGES, rng.integers(0, 2 ** 32, n, dtype=np.uint64).astype(np.uint32)]) def lanes(law, **inputs): """A kernel's law over many lanes, two thousand items a program (the host's queue holds a step's stores).""" n = max(np.size(v) for v in inputs.values()) parts = [law(**{k: (v[i:i + 2000] if np.ndim(v) else v) for k, v in inputs.items()}) for i in range(0, n, 2000)] return np.concatenate(parts) def planted(text, n, inputs, outputs): """A planted program's outputs: `text` compiled by trait and run on n items, each input loaded, not a host field.""" run = trait_run.Program(['error/finalize.trait'], text, {'probe.items': [n]}) for field, values in inputs.items(): run.store('probe.items', field, values) return [o[0] for o in run.steps('probe.main', 1, [('probe.items', f) for f in outputs])]

The kernel¶

The law as the library states it.

In [2]:

print(trait_run.source('error/uniform.trait'))

print(trait_run.source('error/uniform.trait'))

// lusine.error.uniform: a draw as a uniform in [0, 1) (the host's seeded error's uniform): the top 24 bits of the
// 32-bit draw over 2^24, exact in binary32, so a kernel and the host's seeded error agree on it to the bit. Applied to
// lusine.error.draw(seed, step, key, index), it is the host's seeded error's uniform(seed, step, key, index).
// Checked against a reference in numpy, which has its own tests.

namespace lusine.error;

float uniform(bits drawn)
{
  return integer(drawn >> 8) * (1.0 / 16777216.0);
}

The program¶

Each section runs the kernel over a lattice of items, one a draw.

In [3]:

print(trait_run.law_program('error/uniform.trait', 'uniform'))

def uniform(drawn):
    return np.asarray(lanes(lambda **a: trait_run.law('error/uniform.trait', 'uniform', **a), drawn=drawn), dtype=np.float64)

print(trait_run.law_program('error/uniform.trait', 'uniform')) def uniform(drawn): return np.asarray(lanes(lambda **a: trait_run.law('error/uniform.trait', 'uniform', **a), drawn=drawn), dtype=np.float64)

namespace probe;

lusine.version version = { { number = 1 } };

struct item {
  host bits drawn;
  float result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.error.uniform(items.drawn); }

Reference: the host's uniform, bit for bit, and a planted fault caught¶

Twenty thousand draws by the host's draw_bits over random seeds, steps, keys and indices: the kernel's uniform of each is the host's uniform on the same inputs, exactly. The planted program shifts by 9, the top 23 bits over 2²⁴: it misses on nearly every draw, and never reaches ½.

In [4]:

names = ['red', 'jack', 'jack/normal', 'café']
seed, step, index = random_bits(20000), random_bits(20000), random_bits(20000)
name = [names[i % len(names)] for i in range(len(seed))]
drawn = np.array([error.draw_bits(int(a), int(b), n, int(c)) for a, b, n, c in zip(seed, step, name, index)], dtype=np.uint32)
want = np.array([error.uniform(int(a), int(b), n, int(c)) for a, b, n, c in zip(seed, step, name, index)])
assert np.array_equal(uniform(drawn), want)
WRONG = '''namespace probe;

lusine.version version = { { number = 1 } };

struct item { bits drawn; float u; };
lattice<item, 1> items(spacing = 1 m, boundary = wall);

float shifted(bits drawn) { return integer(drawn >> 9) * (1.0 / 16777216.0); }

void main() { items.u = shifted(items.drawn); }
'''
w, = planted(WRONG, len(drawn), {'drawn': drawn}, ['u'])
w = np.asarray(w, dtype=np.float64)
assert np.mean(w != want) > 0.99 and w.max() < 0.5
print(f'{len(drawn)} draws: the kernel is error.py\'s uniform on every one; the shift by 9 misses on {np.mean(w != want):.2%}')

names = ['red', 'jack', 'jack/normal', 'café'] seed, step, index = random_bits(20000), random_bits(20000), random_bits(20000) name = [names[i % len(names)] for i in range(len(seed))] drawn = np.array([error.draw_bits(int(a), int(b), n, int(c)) for a, b, n, c in zip(seed, step, name, index)], dtype=np.uint32) want = np.array([error.uniform(int(a), int(b), n, int(c)) for a, b, n, c in zip(seed, step, name, index)]) assert np.array_equal(uniform(drawn), want) WRONG = '''namespace probe; lusine.version version = { { number = 1 } }; struct item { bits drawn; float u; }; lattice<item, 1> items(spacing = 1 m, boundary = wall); float shifted(bits drawn) { return integer(drawn >> 9) * (1.0 / 16777216.0); } void main() { items.u = shifted(items.drawn); } ''' w, = planted(WRONG, len(drawn), {'drawn': drawn}, ['u']) w = np.asarray(w, dtype=np.float64) assert np.mean(w != want) > 0.99 and w.max() < 0.5 print(f'{len(drawn)} draws: the kernel is error.py\'s uniform on every one; the shift by 9 misses on {np.mean(w != want):.2%}')

20007 draws: the kernel is error.py's uniform on every one; the shift by 9 misses on 100.00%

Limit: the low bits dropped, never 1¶

The low 8 bits are dropped: 0 and 255 give 0, 256 gives 2⁻²⁴, and 2³² − 1 gives 1 − 2⁻²⁴, the largest, below 1. Every uniform is a binary32 exactly, so the device's float and the host's float64 hold the same number.

In [5]:

got = uniform(np.array([0, 255, 256, 0x80000000, M32], dtype=np.uint32))
assert got[0] == 0.0 and got[1] == 0.0 and got[2] == 2.0 ** -24 and got[3] == 0.5 and got[4] == 1.0 - 2.0 ** -24
x = uniform(random_bits(20000))
assert np.array_equal(x.astype(np.float32).astype(np.float64), x) and np.all((x >= 0) & (x < 1))
print(f'0, 255, 256, 2^31, 2^32 - 1 -> {got.tolist()}')

got = uniform(np.array([0, 255, 256, 0x80000000, M32], dtype=np.uint32)) assert got[0] == 0.0 and got[1] == 0.0 and got[2] == 2.0 ** -24 and got[3] == 0.5 and got[4] == 1.0 - 2.0 ** -24 x = uniform(random_bits(20000)) assert np.array_equal(x.astype(np.float32).astype(np.float64), x) and np.all((x >= 0) & (x < 1)) print(f'0, 255, 256, 2^31, 2^32 - 1 -> {got.tolist()}')

0, 255, 256, 2^31, 2^32 - 1 -> [0.0, 0.0, 5.960464477539063e-08, 0.5, 0.9999999403953552]

Curve: flat over [0, 1)¶

Two hundred thousand uniforms of draws keyed "jack" over steps 0 to 199 999 from one seed: each fiftieth of [0, 1) holds a fiftieth of them within 0.0015, about five standard errors (√(p(1 − p) / n) ≈ 0.0003).

In [6]:

n = 200000
steps = np.arange(n, dtype=np.uint32)
drawn = np.array([error.draw_bits(20261005, int(t), 'jack', 0) for t in steps], dtype=np.uint32)
x = uniform(drawn)
share = np.histogram(x, bins=50, range=(0.0, 1.0))[0] / n
assert np.all(np.abs(share - 0.02) < 0.0015), (share.min(), share.max())
print(f'each fiftieth holds {share.min():.4f} to {share.max():.4f}')
fig, ax = plt.subplots(figsize=(7, 2.8))
ax.stairs(share, np.linspace(0, 1, 51))
ax.axhline(0.02, color='grey', lw=0.5)
ax.set_ylim(0, 0.03)
ax.set_xlabel('uniform')
ax.set_ylabel('share of draws')
plt.show()

n = 200000 steps = np.arange(n, dtype=np.uint32) drawn = np.array([error.draw_bits(20261005, int(t), 'jack', 0) for t in steps], dtype=np.uint32) x = uniform(drawn) share = np.histogram(x, bins=50, range=(0.0, 1.0))[0] / n assert np.all(np.abs(share - 0.02) < 0.0015), (share.min(), share.max()) print(f'each fiftieth holds {share.min():.4f} to {share.max():.4f}') fig, ax = plt.subplots(figsize=(7, 2.8)) ax.stairs(share, np.linspace(0, 1, 51)) ax.axhline(0.02, color='grey', lw=0.5) ax.set_ylim(0, 0.03) ax.set_xlabel('uniform') ax.set_ylabel('share of draws') plt.show()

each fiftieth holds 0.0194 to 0.0208

The normal's inputs¶

The host's normal is Box–Muller over the uniforms at indices 2i and 2i + 1 of the stream key + "/normal", through its float64 dlog and dcos. Trait's float is binary32 and its log and cos are the language's own fixed sequences, so the normal itself is not a library kernel: no binary32 law gives the host's float64 normal bit for bit. Its inputs are exact, though: a kernel draws them with lusine.error.draw and uniform, the key error.key("red/normal").

In [7]:

i = np.arange(500, dtype=np.uint32)
k = error.key('red/normal')
drawn = np.array([error.position_hash(error.position_hash(12, k, 7), int(j), 7) for j in np.concatenate([2 * i, 2 * i + 1])],
                 dtype=np.uint32)
u = uniform(drawn)
want = [error.uniform(7, 12, 'red/normal', int(j)) for j in np.concatenate([2 * i, 2 * i + 1])]
assert np.array_equal(u, want)
print(f'the normal\'s {len(u)} uniforms at seed 7, step 12, key "red", indices 0 to 499: error.py\'s, exactly')

i = np.arange(500, dtype=np.uint32) k = error.key('red/normal') drawn = np.array([error.position_hash(error.position_hash(12, k, 7), int(j), 7) for j in np.concatenate([2 * i, 2 * i + 1])], dtype=np.uint32) u = uniform(drawn) want = [error.uniform(7, 12, 'red/normal', int(j)) for j in np.concatenate([2 * i, 2 * i + 1])] assert np.array_equal(u, want) print(f'the normal\'s {len(u)} uniforms at seed 7, step 12, key "red", indices 0 to 499: error.py\'s, exactly')

the normal's 1000 uniforms at seed 7, step 12, key "red", indices 0 to 499: error.py's, exactly