Story of an interesting bug, found (and solved) by Claude Code

Most of the-stats-duck — our open-source DuckDB extension for statistics — is written with Claude Code doing the typing. This post is the story of the most interesting bug of the v0.6.0 release: the kind that returns plausible-looking output while being statistically wrong, the kind that slips past every code review and could sit in production for months.

It turned out to be two bugs, stacked. One was ours. One was the C++ standard library's. Claude Code wrote the first one, then found and fixed both — and the diagnostic trail is a nice little detective story, so we wrote it down.

The setup

v0.6.0 added the random-sampling family: rnorm, rt, rchisq, rf, rgamma, rbeta, rexp, rweibull, rlnorm, rpois. The design is classic inverse-CDF sampling: draw a uniform $u \sim \mathrm{Uniform}(0,1)$, push it through the distribution's quantile function, done. The quantile functions already existed and were cross-checked against R — the r* family just had to compose them with a uniform generator.

The first draft of rexp(rate) looked like every other distribution scalar in the codebase:

static void RExp1Exec(DataChunk &args, ExpressionState &, Vector &result) {
    UnaryExecutor::ExecuteWithNulls<double, double>(
        args.data[0], result, args.size(),
        [](double rate, ValidityMask &mask, idx_t idx) {
            return SafeSample(
                [&](double u) { return stats_duck::ExponentialQuantile(u, rate); },
                mask, idx);
        });
}

If you have ever written a DuckDB scalar function, this is exactly what you would write. UnaryExecutor is the officially sanctioned helper that deals with all the vector plumbing for you. The function was marked VOLATILE, the flag that tells DuckDB "this is not a pure function, don't optimise calls away." Every box ticked.

It compiled. It ran. It returned random-looking numbers.

A smoke test with a smell

Then Claude Code ran the smoke test it runs on every sampler — check the moments against theory:

SELECT round(avg(rexp(1.0)), 2) AS mean,
       round(max(rexp(1.0)), 1) AS mx
FROM range(500000);

An exponential with rate λ = 1 has mean 1.0, and the expected maximum of 500,000 independent draws is a textbook extreme-value fact: about $\ln(500{,}000) + \gamma \approx 13.7$.

We got a mean drifting between 0.92 and 1.1 across runs, and a max of about 8.

Not absurd numbers. Numbers a tired human would wave through — random is random, right? But the mean of 500k draws should sit within a hair of 1.00, and a max of 8 is a right tail that has been amputated. Something was wrong, and one more observation made it interesting: rnorm() — the zero-argument form — had perfect moments, while rgamma(2, 3) was slightly off. The error correlated with which executor pattern the function used, not with which distribution it sampled.

500,000 rows, 245 random numbers

The killer diagnostic was one query:

SELECT count(*) AS total, count(DISTINCT v) AS distinct_cnt
FROM (SELECT rexp(1.0) AS v FROM range(500000));

total: 500000
distinct_cnt: 245

Half a million rows. 245 unique values. And the ratio gives the whole game away: 500,000 / 245 ≈ 2041 — which is DuckDB's vector size of 2048. Each chunk of 2048 rows was getting exactly one random draw, copied 2048 times.

Here is why. DuckDB executes queries in chunks of up to 2048 rows, and UnaryExecutor has a fast path for constant inputs: when the argument is the literal 1.0, the executor reads the constant once, invokes your lambda once, and broadcasts the single result across the chunk. For a pure function, that is a perfectly sound optimisation. For a sampling function, it is catastrophic — and it also explains the clue above: rnorm() takes no arguments, so it never went through UnaryExecutor and had a hand-written per-row loop all along.

But we marked it VOLATILE!

We did. And this is the part worth writing on a wall:

FunctionStability::VOLATILE instructs the planner — don't constant-fold, don't common-subexpression-eliminate. It says nothing to the executor, whose constant-input fast path is a runtime optimisation that fires whether the function is volatile or not.

So VOLATILE correctly stops rnorm() + rnorm() from being folded into 2 * rnorm() at plan time, and then the executor cheerfully replicates one draw across each vector at run time anyway. Two different layers, two different notions of "don't cache this."

The fix: drop the helper and write the per-row loop explicitly via UnifiedVectorFormat, which still handles constant / flat / dictionary input encodings — but in a loop we control:

template <typename Quantile>
static inline void RunPerRow1(DataChunk &args, Vector &result, Quantile &&q) {
    idx_t count = args.size();
    UnifiedVectorFormat a0;
    args.data[0].ToUnifiedFormat(count, a0);
    auto v0 = UnifiedVectorFormat::GetData<double>(a0);
    result.SetVectorType(VectorType::FLAT_VECTOR);
    auto out = FlatVector::GetData<double>(result);
    auto &mask = FlatVector::Validity(result);
    for (idx_t i = 0; i < count; i++) {
        auto i0 = a0.sel->get_index(i);
        if (!a0.validity.RowIsValid(i0)) {
            mask.SetInvalid(i);
            (void)NextUniform(); // keep RNG advance independent of NULLs
            out[i] = 0.0;
            continue;
        }
        out[i] = SafeSampleOne(
            [&](double u) { return q(u, v0[i0]); }, mask, i);
    }
}

(Note the small courtesy in the NULL branch: the RNG advances anyway, so the stream of draws doesn't depend on where the NULLs happen to sit.)

After the fix: 500,000 rows, 500,000 distinct values, mean ≈ 0.997. Bug one down.

The tail that wasn't there

Except the story wasn't over. The distribution's broad shape was now right, but the right tail of rexp(1.0) still came up light. Claude Code compared tail counts against what Exp(1) theory demands for 500k draws:

The largest uniform draw in half a million samples was about 0.999665. The probability of that happening by chance is $(0.999665)^{500{,}000} \approx e^{-167} \approx 10^{-72}$. That is not bad luck. Something was clipping u from above.

The uniform generator was this:

std::uniform_real_distribution<double> dist(
    std::nextafter(0.0, 1.0),
    std::nextafter(1.0, 0.0));
double u = dist(rng);

Textbook C++. Deliberately an open interval (0, 1), so the inverse-CDF transforms downstream — log(1-u), qnorm(u) — can never hit the endpoints. A code reviewer nods and moves on.

But on the MSVC STL our Windows builds ship against, that distribution does not deliver uniform draws across its declared support — the upper tail is squeezed in. Same story for std::generate_canonical<double, 53>. We didn't excavate the implementation to pin the exact mechanism; there are years-old reports against MSVC's <random> with similar symptoms, and we needed a working generator more than a culprit.

So the fix sidesteps the library entirely: take one 64-bit draw from mt19937_64, keep the top 53 bits, scale. Every representable double in [0, 1) at 53-bit precision gets equal probability mass:

static inline double NextUniform() {
    auto &rng = TLSRng();
    uint64_t bits = rng() >> 11;            // keep the top 53 bits
    double u = static_cast<double>(bits) * (1.0 / static_cast<double>(1ULL << 53));
    if (u <= 0.0) {
        u = std::numeric_limits<double>::min();
    }
    return u;
}

The upper bound is 1 − 2⁻⁵³ by construction, so no nudging needed there; the lower side is bumped off exact zero to keep log(0) and qnorm(0) out of reach. After this, the tail counts match Exp(1) theory to three significant figures out to the 10⁻⁵ quantile — and the smoke-test max lands right where extreme-value theory says it should, around 13.7.

(In hindsight, that very first max ≈ 8 was both bugs compressed into a single number: 245 effective draws would put the expected maximum near 6, and the clipped generator couldn't produce anything above 8.0 regardless.)

Both bugs pass code review

What makes this pair worth a blog post is that neither bug is visible in the code.

The first draft of rexp is idiomatic DuckDB — it uses the recommended executor helper and the documented volatility flag. The uniform generator is idiomatic C++ — std::uniform_real_distribution is the API every textbook tells you to use. You could put both snippets in front of ten strong engineers and collect ten approvals.

The bugs only exist statistically. And statistics is exactly the thing a review diff doesn't show. Worse: checking moments alone isn't enough either — a generator that clips its upper tail keeps the mean and standard deviation close to theory while silently ruining every Monte Carlo workload that cares about tails, which is roughly all of them.

The thing that caught both was a boring, mechanical verification loop: after implementing, check the output against theory — moments, then distinct counts, then tail counts — and refuse to shrug when 500k draws average 0.92 instead of 1.00. This is precisely the kind of discipline an agent is good at, because it never gets bored and never decides "eh, random is random." Claude Code wrote the naive first draft, like most of us would have. It also ran the checks that most of us, honestly, would have skipped — then chased the anomaly down through DuckDB's executor internals into the standard library, and wrote the engineering note so the next contributor — human or agent — starts from the pattern that works.

The checklist

Distilled from the wreckage, for anyone shipping sampling functions anywhere:

1. count(*) vs count(DISTINCT v). If a function is supposed to vary per row and these disagree, you have an executor caching problem. Full stop.
2. Tail counts against theoretical probabilities. count(v > threshold) vs n · P(V > threshold). Catches uniform-generator bias that moment checks can't see.
3. Inverse round-trip. pnorm(rnorm()) must be Uniform(0, 1): mean 0.5, stddev $1/\sqrt{12}$.
4. Run it more than once. Thread-local RNGs seed per thread, so symptoms drift between runs. One weird number might be sampling noise; a moving target across reruns is a strong signal.

Trust, but count(DISTINCT)

The fixes landed before v0.6.0 shipped, so every released build samples correctly. Don't take our word for it — this is the shipped build, running in your browser right now:

SELECT
  count(*) AS total
  , count(DISTINCT v) AS distinct_cnt
  , round(avg(v), 3) AS mean
  , round(max(v), 1) AS mx
FROM (SELECT rexp(1.0) AS v FROM range(500000));

Half a million rows, half a million distinct values, mean within a hair of 1, max where the theory puts it.

the-stats-duck is MIT-licensed and runs anywhere DuckDB runs — it's the statistics engine under Bedevere and KoliLang. The full engineering note, with the complete diagnostic trail, lives in the repo. And if there's a distribution or test you need that isn't there yet, get in touch — commissioned functions ship in the open-source extension, so everyone gets them. Carefully sampled.