feat: add version flag, fix --float parsing, and optimize distributions

src/dist.ts

@@ -4,16 +4,34 @@ import { type Options } from "./types";
 // Individual distribution samplers
 // ---------------------------------------------------------------------------
 
-/** Box-Muller transform. Each call consumes 2 uniform randoms. */
+// Cached second normal deviate from Box-Muller transform.
+let spareNormal: number | null = null;
+
+/** Box-Muller with pair caching: uses half the RNG calls of naive Box-Muller. */
 function normalRandom(mean: number, stddev: number): number {
+  if (spareNormal !== null) {
+    const v = spareNormal;
+    spareNormal = null;
+    return mean + stddev * v;
+  }
   let u1 = Math.random();
-  while (u1 === 0) u1 = Math.random(); // avoid log(0)
+  for (let attempt = 0; u1 === 0 && attempt < 100; attempt++) {
+    u1 = Math.random();
+  }
+  if (u1 === 0) u1 = Number.EPSILON;
   const u2 = Math.random();
-  return mean + stddev * Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
+  const r = Math.sqrt(-2 * Math.log(u1));
+  spareNormal = r * Math.sin(2 * Math.PI * u2);
+  return mean + stddev * r * Math.cos(2 * Math.PI * u2);
 }
 
-/** Sum of Bernoulli trials. */
+/** Bernoulli trials; Normal approximation when n>10_000 and np, n(1-p) both >5. */
 function binomialRandom(n: number, p: number): number {
+  if (n > 10_000 && n * p > 5 && n * (1 - p) > 5) {
+    const mean = n * p;
+    const stddev = Math.sqrt(n * p * (1 - p));
+    return Math.max(0, Math.min(n, Math.round(normalRandom(mean, stddev))));
+  }
   let s = 0;
   for (let i = 0; i < n; i++) {
     if (Math.random() < p) s++;
@@ -21,8 +39,11 @@ function binomialRandom(n: number, p: number): number {
   return s;
 }
 
-/** Knuth's algorithm. */
+/** Knuth's algorithm; Normal approximation for λ > 100 (avoids exp underflow). */
 function poissonRandom(lambda: number): number {
+  if (lambda > 100) {
+    return Math.max(0, Math.round(normalRandom(lambda, Math.sqrt(lambda))));
+  }
   const L = Math.exp(-lambda);
   let k = 0;
   let p = 1;
@@ -33,9 +54,9 @@ function poissonRandom(lambda: number): number {
   return k - 1;
 }
 
-/** Inverse CDF. */
+/** Inverse CDF. Caller must ensure λ > 0. */
 function exponentialRandom(lambda: number): number {
-  return -Math.log(Math.random()) / lambda;
+  return -Math.log(Math.random() || Number.EPSILON) / lambda;
 }
 
 /** Urn model — simulate drawing without replacement. */
@@ -58,31 +79,67 @@ function hypergeometricRandom(N: number, K: number, n: number): number {
 // Dispatcher
 // ---------------------------------------------------------------------------
 
-export function generate(opts: Options): number[] {
-  const results: number[] = [];
+export function validateOptions(opts: Options): void {
+  if (!Number.isInteger(opts.count) || opts.count < 1) {
+    throw new Error(`invalid count: ${opts.count}`);
+  }
+  if (!Number.isInteger(opts.decimals) || opts.decimals < 0 || opts.decimals > 100) {
+    throw new Error(`decimals must be 0–100, got ${opts.decimals}`);
+  }
+  switch (opts.dist) {
+    case "uniform":
+      if (opts.min > opts.max) throw new Error(`min (${opts.min}) > max (${opts.max})`);
+      break;
+    case "normal":
+      if (opts.stddev <= 0) throw new Error(`stddev must be > 0, got ${opts.stddev}`);
+      break;
+    case "binomial":
+      if (opts.trials < 0 || !Number.isInteger(opts.trials))
+        throw new Error(`trials must be a non-negative integer, got ${opts.trials}`);
+      if (opts.prob < 0 || opts.prob > 1)
+        throw new Error(`prob must be 0–1, got ${opts.prob}`);
+      break;
+    case "poisson":
+    case "exponential":
+      if (opts.lambda <= 0) throw new Error(`lambda must be > 0, got ${opts.lambda}`);
+      break;
+    case "hypergeometric":
+      if (opts.popSize < 0 || !Number.isInteger(opts.popSize))
+        throw new Error(`population size N must be a non-negative integer, got ${opts.popSize}`);
+      if (opts.successes < 0 || opts.successes > opts.popSize || !Number.isInteger(opts.successes))
+        throw new Error(`successes K must be 0–N, got ${opts.successes} (N=${opts.popSize})`);
+      if (opts.draws < 0 || opts.draws > opts.popSize || !Number.isInteger(opts.draws))
+        throw new Error(`draws n must be 0–N, got ${opts.draws} (N=${opts.popSize})`);
+      break;
+  }
+}
+
+export function* generate(opts: Options): Generator<number> {
+  validateOptions(opts);
   for (let i = 0; i < opts.count; i++) {
+    let v: number;
     switch (opts.dist) {
       case "uniform":
-        results.push(Math.random() * (opts.max - opts.min) + opts.min);
+        v = opts.decimals === 0
+          ? Math.floor(Math.random() * (opts.max - opts.min + 1)) + opts.min
+          : Math.random() * (opts.max - opts.min) + opts.min;
         break;
       case "normal":
-        results.push(normalRandom(opts.mean, opts.stddev));
+        v = normalRandom(opts.mean, opts.stddev);
         break;
       case "binomial":
-        results.push(binomialRandom(opts.trials, opts.prob));
+        v = binomialRandom(opts.trials, opts.prob);
         break;
       case "poisson":
-        results.push(poissonRandom(opts.lambda));
+        v = poissonRandom(opts.lambda);
         break;
       case "exponential":
-        results.push(exponentialRandom(opts.lambda));
+        v = exponentialRandom(opts.lambda);
         break;
       case "hypergeometric":
-        results.push(
-          hypergeometricRandom(opts.popSize, opts.successes, opts.draws),
-        );
+        v = hypergeometricRandom(opts.popSize, opts.successes, opts.draws);
         break;
     }
+    yield opts.decimals > 0 ? v : Math.round(v);
   }
-  return results;
 }