1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119

# rare.py: the arithmetic of rare failures. Standard library only; checks itself and prints PASS.
import math
import sys


def p_seen(p, n):
    """Chance that n independent runs show a bug that fails each run with probability p."""
    return 1.0 - (1.0 - p) ** n


def runs_needed(p, conf):
    """Fewest runs that show a rate-p bug at least once with probability conf."""
    return math.ceil(math.log(1.0 - conf) / math.log(1.0 - p))


def upper_if_clean(n, conf):
    """Exact one-sided upper bound on the failure rate after n runs with no failure."""
    return 1.0 - (1.0 - conf) ** (1.0 / n)


def binom_pmf(k, n, p):
    if p <= 0.0:
        return 1.0 if k == 0 else 0.0
    if p >= 1.0:
        return 1.0 if k == n else 0.0
    log = (math.lgamma(n + 1) - math.lgamma(k + 1) - math.lgamma(n - k + 1)
           + k * math.log(p) + (n - k) * math.log(1.0 - p))
    return math.exp(log)


def binom_cdf(k, n, p):
    return sum(binom_pmf(i, n, p) for i in range(k + 1))


def solve(f, target, lo=0.0, hi=1.0):
    """Bisection for an increasing f on [lo, hi]."""
    for _ in range(100):
        mid = (lo + hi) / 2
        if f(mid) < target:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2


def exact_interval(k, n, conf=0.95):
    """Clopper-Pearson two-sided interval for k failures in n runs."""
    a = (1.0 - conf) / 2
    lo = 0.0 if k == 0 else solve(lambda p: 1.0 - binom_cdf(k - 1, n, p), a)
    hi = 1.0 if k == n else solve(lambda p: -binom_cdf(k, n, p), -a)
    return lo, hi


def fisher(k1, n1, k2, n2):
    """Two-sided Fisher exact test: is k1/n1 really different from k2/n2?"""
    total_k, total_n = k1 + k2, n1 + n2

    def hyper(x):
        return math.exp(math.lgamma(total_k + 1) - math.lgamma(x + 1) - math.lgamma(total_k - x + 1)
                        + math.lgamma(total_n - total_k + 1) - math.lgamma(n1 - x + 1)
                        - math.lgamma(total_n - total_k - n1 + x + 1)
                        - (math.lgamma(total_n + 1) - math.lgamma(n1 + 1) - math.lgamma(n2 + 1)))

    lo, hi = max(0, n1 - (total_n - total_k)), min(n1, total_k)
    seen = hyper(k1)
    return min(1.0, sum(hyper(x) for x in range(lo, hi + 1) if hyper(x) <= seen * (1 + 1e-9)))


def check(cond, what):
    if not cond:
        print("FAIL:", what)
        sys.exit(1)


# 1. How many runs does it take to see a bug of rate p at least once?
print("runs needed to see the bug at least once")
print("  rate      95%     99%")
for p in (0.10, 0.013, 0.01, 0.001):
    print(f"  {p:6.3%} {runs_needed(p, 0.95):6d} {runs_needed(p, 0.99):7d}")

# 2. After n clean runs, how high could the rate still be? (exact, and the rule of three)
print("clean runs and what they rule out (95%, one-sided)")
for n in (10, 29, 100, 300, 1000, 3000):
    print(f"  0 in {n:4d}: rate below {upper_if_clean(n, 0.95):6.2%}   rule of three: {3 / n:6.2%}")

# 3. dbglab's B4, as measured: 100 failures in 1,000 seeds; stressed: 177 in 200
lo, hi = exact_interval(100, 1000)
print(f"B4, 100 in 1000: 10.0%, exact 95% interval {lo:.2%} to {hi:.2%}")
lo2, hi2 = exact_interval(177, 200)
print(f"B4 stressed, 177 in 200: 88.5%, exact 95% interval {lo2:.1%} to {hi2:.1%}")
print(f"  clean runs to rule B4 out at 95%: {runs_needed(0.10, 0.95)} default, "
      f"{runs_needed(lo2, 0.95)} stressed (using the interval's low end, {lo2:.1%})")

# the FIFO writer bug of the previous lessons: 13 in 1,000 nightly seeds, 196 in 200 with bursty stalls
lo3, hi3 = exact_interval(196, 200)
print(f"FIFO bug, 13 in 1000 nightly: {runs_needed(0.013, 0.95)} clean runs to rule out; "
      f"196 in 200 stressed (low end {lo3:.1%}): {runs_needed(lo3, 0.95)}")

# 4. Is one rate really different from another?
for k1, n1, k2, n2, what in ((21, 200, 21, 200, "default vs +OPS=6000"),
                             (21, 200, 177, 200, "default vs biased zero fills"),
                             (100, 1000, 0, 1000, "before vs after a fix"),
                             (3, 300, 0, 100, "3 in 300 before, 0 in 100 after"),
                             (12, 400, 4, 400, "12 in 400 before, 4 in 400 after")):
    print(f"  {what}: {k1}/{n1} vs {k2}/{n2}, p = {fisher(k1, n1, k2, n2):.3g}")

# self-checks against values worked by hand or known from the literature
check(runs_needed(0.10, 0.95) == 29 and runs_needed(0.10, 0.99) == 44, "runs for a 10% bug")
check(runs_needed(0.001, 0.95) == 2995, "runs for a 0.1% bug")
check(abs(upper_if_clean(1000, 0.95) - 0.002991) < 1e-5, "upper bound after 1000 clean runs")
check(all(abs(upper_if_clean(n, 0.95) * n / 3 - 1) < 0.05 for n in (100, 1000, 10000)), "rule of three")
check(abs(p_seen(0.10, 29) - 0.953) < 1e-3, "chance of seeing a 10% bug in 29 runs")
check(abs(lo - 0.0821) < 5e-4 and abs(hi - 0.1203) < 5e-4, "exact interval for 100 in 1000")
check(abs(binom_cdf(1000, 1000, 0.3) - 1.0) < 1e-9, "binomial sums to one")
check(fisher(21, 200, 21, 200) > 0.99, "equal rates are not different")
check(fisher(100, 1000, 0, 1000) < 1e-20, "10% vs 0% over 1000 runs each")
check(fisher(3, 300, 0, 100) > 0.5, "3 in 300 vs 0 in 100 proves nothing")
check(runs_needed(lo3, 0.95) == 2, "two clean runs rule out a 95% bug")
print("PASS: rare-failure arithmetic checks out")