Simulating Auctions in C

An auction simulator is four hundred lines and produces a number. The number is almost always wrong the first time, and it is wrong in a way that looks right: a single run ranks two mechanisms, prints a decimal, and says nothing about whether the ranking would survive a different seed. This note builds the simulator — value draws from six distributions using nothing but libc, a bidder as a function pointer, one driver per auction format — and then spends most of its length on the part that actually decides whether the output is a result or a coin flip: fixed seeds, many replications, reported confidence intervals, and common random numbers. The centrepiece is a measured demonstration that 41.3 % of 1,000-trial seeds rank two reserve prices backwards, that one specific seed does so at t = −3.10 with a 95 % confidence interval that excludes the true answer entirely, and that the same seed with common random numbers reports the correct sign and honestly says it cannot tell. Along the way: rand() is not portable, and the POSIX specification itself says so, in a section titled “Generating the Same Sequence on Different Machines.”

This is the engineering companion to Ad Auctions and GSP, First-Price and Second-Price Auctions, Revenue Equivalence and Reserve Prices and Optimal Auctions — those notes carry the theory; this one carries the code and the statistics. It is also the auction-shaped sibling of Running a Strategy Tournament in C, which covers the function-pointer strategy interface, history representation and cross-compiler determinism for iterated games. Those sections are not repeated here. What is genuinely new in the auction setting is the sampling — you need draws from continuous distributions, which a Prisoner’s Dilemma tournament never does — and the variance structure, because auction formats with identical expected revenue have wildly different revenue variance, which changes how many trials you need by more than an order of magnitude.

Mental Model — Three Layers, Each Independently Replaceable

flowchart TB
    subgraph L1["Layer 1 · Randomness"]
        R["rng (PCG32)<br/>rng_unit → [0,1)"]
        D["dist → value draws<br/>inverse CDF · polar method"]
    end
    subgraph L2["Layer 2 · Behaviour"]
        E["env { n, reserve, dist* }"]
        B["bid_fn(value, env*, rng*)<br/>one function pointer per strategy"]
    end
    subgraph L3["Layer 3 · Rules"]
        F["fmt_fn(bids[], values[], env*)<br/>→ { revenue, winner_value, sold }"]
    end
    subgraph L4["Layer 4 · Statistics"]
        W["Welford accumulator<br/>mean · sd · sem"]
        CI["95 % CI · paired t · CRN"]
    end
    R --> D --> B
    E --> B
    B --> F --> W --> CI

The simulator’s four layers. The insight: the value draw and the bid function must be separate objects, because the single most useful experimental technique in this note — common random numbers — works by holding layer 1 fixed while varying layers 2 and 3. A simulator that generates values inside the format driver cannot do that, and will need roughly seventeen times as many trials for the same answer.

The env struct is the piece people leave out and regret. A bidder’s equilibrium strategy depends on how many rivals there are and what distribution their values come from — the first-price equilibrium bid for n bidders with values uniform on [0,1] is b(v) = (n−1)/n · v, and both n and “uniform” are in that formula. Passing the environment to the bid function rather than baking it into a global is what lets a single build sweep n from 2 to 20:

typedef struct {
    int    n;          /* number of bidders in THIS auction */
    double reserve;    /* reserve price, 0 if none */
    dist  *vd;         /* the value distribution, common knowledge */
} env;
 
typedef double (*bid_fn)(double value, const env *e, rng *r);
typedef struct { const char *name; bid_fn bid; } bidder;

The signature differs from the tournament’s in exactly one way that matters: a bidder sees only its own value, never anyone else’s. In Running a Strategy Tournament in C the strategy receives a const view * onto the shared history precisely so TIT FOR TAT can copy its opponent; here, passing anything about the other bidders would be handing a sealed-bid bidder the sealed bids. The absence of that parameter is the model.

The Random Source, and Why rand() Is Disqualified

The claim “rand() is not reproducible across platforms” is usually asserted. It can be proved from the specification, and the specification is unusually direct about it. Here is the POSIX Issue 8 page for rand (IEEE Std 1003.1-2024, Open Group, fetched 2026-08-28). The DESCRIPTION guarantees only:

“The rand() function shall compute a sequence of pseudo-random integers in the range [0,{RAND_MAX}] with a period of at least 2³². …If srand() is then called with the same seed value, the sequence of pseudo-random numbers shall be repeated.”

Read that carefully. The repeatability guarantee is within one implementation, for one seed. Nothing constrains what the sequence is. And the specification then includes an EXAMPLES section headed, verbatim:

“Generating the Same Sequence on Different Machines” “The following code defines a pair of functions that could be incorporated into applications wishing to ensure that the same sequence of numbers is generated across different machines.”

static unsigned long next = 1;
int myrand(void)   /* RAND_MAX assumed to be 32767. */
{ next = next * 1103515245 + 12345;
  return((unsigned)(next/65536) % 32768); }
void mysrand(unsigned seed) { next = seed; }

The standard’s own remedy for cross-machine reproducibility is “stop calling rand() and vendor a generator.” That is the whole argument, and it comes from the specification rather than from folklore.

How different are they in practice? Compiling that exact myrand alongside glibc’s rand(), seeding both with 1:

seed 1: first 8 values
  glibc rand() : 1804289383 846930886 1681692777 1714636915 1957747793 424238335 719885386 1649760492
  C11 example  :      16838       5758      10113      17515      31051       5627      23010       7419
  RAND_MAX glibc = 2147483647, C11 example range = 0..32767
  agreement over 100000 draws (mod 32768): 7 (0.0070%)

They agree 7 times in 100,000 — which is 1/32768 × 100000 ≈ 3, i.e. chance. glibc’s rand() is not that LCG at all; it is a large additive-feedback generator with RAND_MAX = 2³¹−1, and the man 3 rand page notes that rand() “is not reentrant, since it uses hidden state that is modified on each call” (man7.org). POSIX’s own APPLICATION USAGE says these functions “should be avoided whenever non-trivial requirements (including safety) have to be fulfilled.”

The consequence for an auction experiment is concrete: a paper that says “seed 12345, 10,000 trials, first-price revenue 0.6672” is unreproducible by anyone not running your libc, and irreproducible by you after a distribution upgrade.

The second charge: rand() also fails the tests

Non-reproducibility is the argument from the specification. There is an independent argument from measurement, and it is the one that matters if you only ever run on one machine.

How generator quality is actually assessed. The standard instrument is TestU01 — L’Ecuyer & Simard, “TestU01: A C library for empirical testing of random number generators”, ACM Transactions on Mathematical Software 33(4), Article 22 (2007). A test in TestU01 performs a statistically well-understood task using the candidate generator and checks whether the outcome is plausible. L’Ecuyer’s own worked example (WSC 2015 §2) is the collision test: partition the unit hypercube (0,1)ˢ into k = 2^{ds} equal subcubes, generate n points, and count C, the number of points landing in a subcube that was already occupied. Under the null hypothesis that the generator emits independent U(0,1) values, C is approximately Poisson with mean λ = n²/(2k), so you can compute both tails:

p⁺ = P[X ≥ c]   very small ⇒ too many collisions  ⇒ lack of uniformity
p⁻ = P[X ≤ c]   very small ⇒ too few collisions   ⇒ lack of INDEPENDENCE (too uniform)

Note the two-sidedness — a generator that spreads its points too evenly is as broken as one that clumps, and only the two-tailed test catches it. L’Ecuyer’s threshold for declaring failure is a p-value below 10⁻¹⁰ or 10⁻¹⁵; a middling value around 10⁻³10⁻⁸ means “rerun on a longer, disjoint segment of the cycle until the answer is clear.”

Those individual tests are bundled into the batteries everyone quotes: SmallCrush, Crush, BigCrush (for (0,1) reals) and Rabbit (for bit sequences). A generator passing all three Crush batteries is called Crush-resistant. Vigna, whose xoshiro/xoroshiro family is the other generator this project uses, states his protocol precisely (prng.di.unimi.it): run BigCrush from 100 equispaced points of the generator’s state space, count a failure as any test whose p-value statistic falls outside [0.001, 0.999], and call a failure systematic if it occurs at all 100 points. He also runs the suite on the bit-reversed generator, because TestU01 is a 32-bit suite with a known bias toward high bits — a generator with weak low bits will sail through unless you reverse it.

Where rand() lands. Figure 1(a) of O’Neill’s PCG technical report (HMC-CS-2014-0905) plots SmallCrush failures for the generators in widest use, and the bars include 32-bit LCG, Knuth '81 (LCG), Unix drand48 (LCG), Minstd (LCG), Kreutzer86, Unix random, XorShift 32 and XorShift* 32every one of them registers failures, and Unix random (the additive-feedback generator that glibc’s rand() is built on) is among the worst in the group. O’Neill’s summary is blunt and worth quoting because it sets the bar: “many of the generators in widest use cannot withstand even minimal statistical scrutiny, sometimes failing statistical tests after generating only a few thousand numbers”, and “the combined tests require less than fifteen seconds to run.” Her caption adds the standard everyone should apply: “Even one failure indicates a problem.”

How much of a generator’s period you may actually spend. This is the part of the RNG literature that is directly a Monte Carlo budgeting question and is almost never stated in simulation write-ups. O’Neill derives it from the generalized birthday problem: a generator with b bits of state emitting r-bit outputs runs out of legal repeats as you approach its period, a condition she names being overtaxed. Applying the calculation test-by-test to TestU01 gives clean thresholds — even a mathematically ideal generator “can reasonably fail SmallCrush when b < 32, fail Crush when b < 35, and fail BigCrush when b < 36”. For linear congruential generators specifically she quotes L’Ecuyer & Simard’s rule that an LCG is free of birthday-test issues only while

n < 16 · p^(1/3)        n = numbers consumed, p = period

and draws the consequence that even a 128-bit LCG can show statistical flaws after fewer than 2⁴⁷ draws — “an algorithm could plausibly use one number per nanosecond and 2⁴⁷ nanoseconds is less than two days.”

Put that cube root against this note’s own experiments. The reversal experiment below consumes 1,000 seeds × 100,000 trials × 2 configurations × ~7 draws ≈ 1.4 × 10⁹ ≈ 2³⁰·⁴ numbers. Against PCG32’s 2⁶⁴ period the rule gives a budget of 16 · 2^{64/3} ≈ 2²⁵which this experiment exceeds by five orders of magnitude. That is not a reason to panic (PCG is not a raw LCG; the output permutation is the entire point of the design, and the rule is stated for bare LCGs), but it is exactly the calculation that should be done before quoting a result, and it is the reason the note runs the same experiment through an unrelated generator family as a cross-check.

Uncertain

Verify: whether the n < 16·p^{1/3} budget rule applies to PCG32 as it does to a raw LCG. Reason: the rule is stated by O’Neill (quoting L’Ecuyer & Simard) for LCGs, and PCG32 is an LCG with an output permutation that specifically targets the flaws the rule describes. O’Neill’s overtaxing analysis is framed in terms of state bits b and output bits r and would give PCG32 (b = 64, r = 32) a much larger safe budget, but I did not find an explicit statement of the safe draw count for PCG32 in the report. The primary source that would settle the general form is L’Ecuyer & Simard’s 2007 TOMS paper, which I could not retrieve: iro.umontreal.ca serves an Anubis anti-bot interstitial (HTTP 200, title “Making sure you’re not a bot!”) in place of the PDF, so both this claim and the collision-test description above rest on O’Neill’s and L’Ecuyer’s later restatements rather than on the 2007 original. To resolve: read the TOMS paper, or run SmallCrush against PCG32 locally with 1.4 × 10⁹ draws.

PCG32, vendored

This simulator uses PCG32 — a 64-bit linear congruential state with an output permutation, by Melissa O’Neill. It is chosen here partly because it is small, well documented, and public, and partly because it is deliberately different from the xoshiro256** used by Running a Strategy Tournament in C: running the same experiment through two unrelated generator families is a cheap and real check that a result is not an artefact of one stream. The core is nine lines, transcribed from the reference pcg_basic.c (fetched 2026-08-28) and verified character-for-character against it:

typedef struct { uint64_t state, inc; } rng;
 
static inline uint32_t rng_u32(rng *r)
{
    uint64_t old = r->state;
    r->state = old * 6364136223846793005ULL + r->inc;   /* the LCG step        */
    uint32_t xorshifted = (uint32_t)(((old >> 18u) ^ old) >> 27u);  /* permute  */
    uint32_t rot = (uint32_t)(old >> 59u);              /* ...by a random rot  */
    return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
}
static void rng_seed(rng *r, uint64_t seed, uint64_t seq)
{
    r->state = 0u; r->inc = (seq << 1u) | 1u;
    rng_u32(r); r->state += seed; rng_u32(r);
}

Line by line: the state advances by a 64-bit LCG with the standard multiplier; a plain LCG has famously weak low bits, so the output is a permutation of the state rather than the state itself — XOR-fold the high bits down, then rotate by an amount taken from the top five bits of the same word, which is what makes the permutation data-dependent and defeats the lattice structure of the raw LCG. rng_seed takes two arguments, and the second is the reason PCG is a good fit here: seq selects a stream. inc is forced odd ((seq << 1) | 1) so the LCG has full period, and two generators seeded with the same seed but different seq produce different, independent-looking sequences. That gives an honest way to say “run the first-price format on stream 1 and the second-price format on stream 2 with the same seed” — which is exactly the wrong thing to do statistically, as measured below, but you need to be able to express it in order to measure how wrong it is.

Converting to a double is the usual trap. rng_u32 gives 32 bits; a double has 53 bits of mantissa. Two draws, shifted and concatenated:

static inline double rng_unit(rng *r)          /* uniform on [0,1) */
{
    uint64_t hi = (uint64_t)rng_u32(r) >> 5;   /* 27 bits */
    uint64_t lo = (uint64_t)rng_u32(r) >> 6;   /* 26 bits */
    return (double)((hi << 26) | lo) * 0x1.0p-53;
}
static inline double rng_unit_pos(rng *r)      /* uniform on (0,1] — never 0 */
{
    double u; do { u = rng_unit(r); } while (u == 0.0); return u;
}

The multiplier is written as the hexadecimal float literal 0x1.0p-53 — exactly 2⁻⁵³, exactly representable, no decimal-to-binary rounding anywhere. Every value the function can return is exactly representable as a double, so the conversion introduces no bias.

Re-verified 2026-08-29. The rng_u32, rng_seed and rng_unit listings above were compiled verbatim (gcc -O2 -std=gnu11) and checked against the reference implementation’s own canonical demo seeding, pcg32_srandom_r(&rng, 42u, 54u):

PCG32 seed(42,54) first 6 outputs: 0xa15c02b7 0x7b47f409 0xba1d3330 0x83d2f293 0xbfa4784b 0xcbed606e
rng_unit over 2,000,000 draws: min 9.1955e-07  max 0.99999976366  mean 0.499943  (expect 0.5)

and the pcg_basic.c source was re-fetched and diffed line by line: state = old*6364136223846793005ULL + inc, xorshifted = ((old >> 18u) ^ old) >> 27u, rot = old >> 59u, (xorshifted >> rot) | (xorshifted << ((-rot) & 31)), and inc = (initseq << 1u) | 1u all match character for character.

rng_unit_pos exists for one reason and it is not paranoia. Inverse-CDF sampling of an exponential is −log(U)/λ, and log(0) is −∞:

-log(0.0)            = inf
sqrt(-2*log(0.0))    = inf
BM sample with u1=0  = inf
smallest nonzero 53-bit unit draw = 1.1102230246251565e-16 (prob of exact 0 = 2^-53 = 1.11e-16)

The probability of drawing exactly zero is 2⁻⁵³ ≈ 1.11 × 10⁻¹⁶ per call. At 10⁹ draws per experiment and a thousand experiments, that is one chance in nine thousand of producing an inf that silently poisons a mean — and a mean of inf does not look like a bug, it looks like a crash three functions later. The rejection loop costs one predictable branch and removes the failure mode permanently.

Value Distributions from libc Alone

Auction theory is written in terms of a value distribution F, and every interesting theorem is a statement about FRevenue Equivalence needs independent private values from a common F; Reserve Prices and Optimal Auctions needs F’s hazard rate; the winner’s curse needs correlation. A simulator that can only draw uniforms can test one corner of the theory.

Two techniques cover everything needed, and both are a few lines.

Inverse transform sampling

If U ~ Uniform(0,1) then X = F⁻¹(U) has distribution F. That is the whole method, and it works whenever you can invert the CDF in closed form:

DistributionCDF F(x)DrawLine of C
Uniform [a,b](x−a)/(b−a)a + (b−a)Ud->a + (d->b - d->a) * rng_unit(r)
Exponential, rate λ1 − e^{−λx}−ln(U)/λ-log(rng_unit_pos(r)) / d->a
Pareto, shape α, scale x_m1 − (x_m/x)^αx_m · U^{−1/α}d->b / pow(rng_unit_pos(r), 1.0/d->a)
Power, F(x)=(x/b)^α on [0,b](x/b)^αb · U^{1/α}d->b * pow(rng_unit(r), 1.0/d->a)
Lognormal(µ,σ)exp(µ + σZ)exp(d->a + d->b * normal(r))

Note the substitution of 1−U for U throughout — since U and 1−U have the same distribution, −ln(U) is as good as −ln(1−U) and saves a subtraction. This is standard and correct, but it is not free: it means U = 0 maps to the tail rather than to zero, which is precisely why rng_unit_pos is required for exponential and Pareto and merely optional for the power distribution.

The power distribution deserves a mention because it is the most useful non-uniform family for auction work. F(x) = x^α on [0,1] has a closed-form first-price equilibrium — b(v) = v·α(n−1) / (α(n−1)+1) — so it gives a second independent closed-form cross-check for the whole pipeline, with α = 1 recovering the uniform case. Pareto matters for the opposite reason: with shape α ≤ 1 the mean does not exist, and with α ≤ 2 the variance does not exist, so a sample mean of revenue converges slowly or not at all. Every confidence interval in this note assumes a finite variance; drawing values from a heavy-tailed distribution quietly invalidates that assumption, and the symptom is a CI that refuses to shrink like 1/√n.

Normal draws: Box–Muller and the polar method

The normal CDF has no elementary inverse, so inverse transform is out. The classical answer is Box–Muller: take two independent uniforms and map them to two independent standard normals through polar coordinates,

Z₁ = √(−2 ln U₁) · cos(2π U₂)
Z₂ = √(−2 ln U₁) · sin(2π U₂)

Symbol by symbol: √(−2 ln U₁) is the radius — an exponentially distributed squared-radius is exactly what a 2-D standard normal has — and 2π U₂ is a uniformly random angle. A uniformly random direction at a Rayleigh-distributed radius is a 2-D standard normal, whose two Cartesian components are independent standard normals. Two uniforms in, two normals out.

The Marsaglia polar method does the same thing while avoiding sin and cos. Rejection-sample a point (u,v) uniformly inside the unit disc; then s = u²+v² is uniform on (0,1) and (u,v)/√s is a uniformly random direction already in Cartesian form:

static double d_normal_std(dist *d, rng *r)
{
    if (d->has_spare) { d->has_spare = 0; return d->spare; }   /* return the twin */
    double u, v, s;
    do {
        u = 2.0 * rng_unit(r) - 1.0;
        v = 2.0 * rng_unit(r) - 1.0;
        s = u * u + v * v;
    } while (s >= 1.0 || s == 0.0);          /* reject outside the disc, and s=0 */
    double m = sqrt(-2.0 * log(s) / s);
    d->spare = v * m; d->has_spare = 1;      /* cache the second normal */
    return u * m;
}

The has_spare cache is not an optimization detail, it is a correctness-of-reproducibility detail. Both methods produce normals in pairs. Throwing the second away doubles the number of uniforms consumed, which is fine — until you compare two configurations that consume different numbers of uniforms, at which point their random streams diverge and common random numbers stop working. Deciding how many uniforms a draw consumes is part of the experimental design.

Timing both on this machine (gcc 16.1.1, -O2, 2 × 10⁷ draws, /proc/loadavg 8.80 10.50 7.83 on a 32-core box — the machine was heavily loaded, so treat these as ratios, not absolutes):

Drawns/drawNotes
rng_unit (two PCG32 words → double)1.69the floor
Exponential, inverse CDF (-log(U))4.22one log
Normal, Marsaglia polar7.05~1.27 uniforms/draw amortized, one log, one sqrt
Normal, Box–Muller trig form19.35one log, one sqrt, one cos2.7× slower

The polar method is 2.7× faster despite rejecting about 21.5 % of its candidate pairs (the disc occupies π/4 ≈ 78.5 % of the square). Transcendentals dominate: cos costs more than the rejection loop saves.

Verifying the draws against closed forms

Every distribution has a mean you can write down, so the first thing the simulator does is check itself. Two million draws per distribution, seed 1:

dist               mean         sd        min        max  closed-form mean
uniform        0.500146   0.288664   0.000001   1.000000    0.500000
exponential    0.999778   1.000418   0.000000  13.696726    1.000000
normal+        1.000005   0.249962   0.000296   2.204168    1.000000
lognormal      1.133162   0.604312   0.082618  11.115440    1.133148
pareto         1.499948   0.852602   1.000000  96.117868    1.500000
power          0.666770   0.235720   0.001061   1.000000    0.666667

Every mean matches to at least three decimals and the lognormal to five (exp(µ + σ²/2) = exp(0.125) = 1.1331485). Two further things in that table are worth reading rather than skimming. The uniform’s sd is 0.288664 against the closed form 1/√12 = 0.2886751 — a check of the second moment, not just the first. And the Pareto’s maximum is 96.1 on a distribution with mean 1.5: with shape α = 3 the tail is thin enough for the mean and variance to exist, and it still threw a sample 64× the mean in two million draws. That is what a heavy tail looks like before it becomes a problem.

Format Drivers

An auction format is a pure function from bids and values to an outcome. Keeping values in the signature alongside bids is what lets the simulator report efficiency — did the bidder with the highest value actually win? — which is a different question from revenue and often the more interesting one.

typedef struct { double revenue; double winner_value; int sold; } outcome;
typedef outcome (*fmt_fn)(const double *bids, const double *vals, const env *e);
 
static outcome f_first(const double *b, const double *v, const env *e)
{
    int i1, i2; top2(b, e->n, &i1, &i2); (void)i2;
    outcome o = {0,0,0};
    if (b[i1] >= e->reserve) { o.revenue = b[i1]; o.winner_value = v[i1]; o.sold = 1; }
    return o;
}
static outcome f_second(const double *b, const double *v, const env *e)
{
    int i1, i2; top2(b, e->n, &i1, &i2);
    outcome o = {0,0,0};
    if (b[i1] >= e->reserve) {
        double p = b[i2] > e->reserve ? b[i2] : e->reserve;   /* reserve is a floor */
        o.revenue = p; o.winner_value = v[i1]; o.sold = 1;
    }
    return o;
}
static outcome f_allpay(const double *b, const double *v, const env *e)
{
    int i1, i2; top2(b, e->n, &i1, &i2); (void)i2;
    outcome o = {0,0,0};
    for (int i = 0; i < e->n; i++) o.revenue += b[i];    /* everyone pays */
    o.winner_value = v[i1]; o.sold = 1;
    return o;
}

top2 is a single linear pass finding the highest and second-highest bid. Sorting is the obvious implementation and the wrong one: you need two order statistics, not n, and O(n) beats O(n log n) on the hot path of a hundred-million-auction sweep. It also sidesteps a subtle reproducibility hazard — qsort is not required to be stable, so ties between equal bids could be broken differently by different libc implementations, and equal bids are common when bid functions are deterministic and values are drawn from a discretized grid.

The reserve price appears in both drivers and it is easy to get subtly wrong. In a second-price auction with reserve r, the winner pays max(second-highest bid, r), not r when the second bid is below it and not the second bid when it is below r. Getting this backwards produces a revenue curve that is monotone in r, which looks plausible and is wrong.

Three formats, one shared value draw, 500,000 trials, n = 5, values uniform [0,1], every bidder playing its symmetric equilibrium strategy (b(v)=v for second price, (n−1)/n·v for first price, (n−1)/n·vⁿ for all-pay):

n=5 trials=500000  closed form (n-1)/(n+1) = 0.666667
  first-price   0.666666 +- 0.000313 (sd 0.112743)
  second-price  0.666722 +- 0.000494 (sd 0.178303)
  all-pay       0.666835 +- 0.001246 (sd 0.449387)

Revenue Equivalence confirmed to four decimals across three structurally different mechanisms — and the standard deviations differ by a factor of four. That is the observation the rest of this note is built on.

xychart-beta
    title "Same expected revenue, very different spread (n=5, uniform values)"
    x-axis ["first-price", "second-price", "all-pay"]
    y-axis "std. deviation of per-auction revenue" 0 --> 0.5
    bar [0.112743, 0.178303, 0.449387]

Per-auction revenue standard deviation for three revenue-equivalent formats. The insight: equal means do not imply equal measurement cost. To pin revenue to ±0.001 with 95 % confidence you need (1.96·sd/0.001)² trials — about 49,000 for first price, 122,000 for second price, and 776,000 for all-pay. Budgeting the same number of trials for every format under test gives you three answers of three different qualities and no warning.

Why the spread differs is worth understanding rather than memorizing. First-price revenue is (n−1)/n × max(v): a scaled maximum, and the maximum of five uniforms is concentrated near 1. Second-price revenue is the second order statistic, which sits lower and moves around more. All-pay revenue is the sum of five bids, but the bid function vⁿ is so convex that almost all of the mass comes from whichever bidder happens to draw near 1 — the sum is dominated by one wildly variable term.

Statistics Done Properly

Everything above is plumbing. This is where a simulator becomes an experiment.

The accumulator

Never accumulate sum and sum_of_squares and subtract. Σx² − (Σx)²/n subtracts two nearly equal large numbers, and how much of the answer survives depends on how far the data sit from zero. Welford’s online algorithm costs one extra multiply and is numerically stable:

typedef struct { long n; double mean, m2; } acc;
static void acc_add(acc *a, double x)
{
    a->n++;
    double d = x - a->mean;
    a->mean += d / (double)a->n;
    a->m2   += d * (x - a->mean);      /* note: the UPDATED mean, deliberately */
}
static double acc_var(const acc *a){ return a->n > 1 ? a->m2 / (double)(a->n - 1) : 0.0; }
static double acc_sd (const acc *a){ return sqrt(acc_var(a)); }
static double acc_sem(const acc *a){ return a->n ? acc_sd(a) / sqrt((double)a->n) : 0.0; }

The asymmetry in m2 += d * (x - a->mean) — the first factor uses the old mean, the second the new — is not a typo and is the entire trick.

Correction — an earlier revision of this note got the reason wrong

This section previously claimed that on revenue values around 0.67 the naive Σx² formula “loses most of its significant digits by ~10⁸ samples.” That is measurably false, and I measured it. Running both accumulators on the actual second-price revenue stream (n = 5, uniform values, mean ≈ 0.667, sd ≈ 0.178):

           n        Welford var          naive var  rel.err naive
        1000     0.033122095734     0.033122095734       3.81e-14
      100000     0.031728003917     0.031728003917       1.45e-13
    10000000     0.031777321382     0.031777321382       8.59e-13
   100000000     0.031747463455     0.031747463456       1.12e-11
  1000000000     0.031746730424     0.031746730424       1.17e-11

At 10⁹ samples the naive formula still agrees with Welford to eleven significant digits. Sample count is essentially not the driver. The driver is the coefficient of variation, mean / sd — the two subtracted quantities are close in proportion to (mean/sd)², which is the condition number of the subtraction. Holding n = 10⁷ fixed and shifting the same revenue stream by a constant offset isolates it exactly:

N = 10⁷, x = offset + second-price revenue      sd(x) ≈ 0.178 throughout
      offset    mean/sd        Welford var          naive var   naive rel.err
       0e+00        3.7     0.031777321382           0.031777        8.59e-13
       1e+02      564.7     0.031777321382           0.031777        2.77e-08
       1e+04    56101.0     0.031777321382           0.031783        1.75e-04
       1e+05   560976.0     0.031777321382           0.030331        4.55e-02
       1e+06  5609726.0     0.031777321384           0.145818        3.59e+00
       1e+08 560972224.1     0.031777321613       -7673.899366        2.41e+05

Welford is still the right choice, for a different reason than the note gave. Correct statement: on a mean-zero-ish quantity like [0,1] auction revenue the naive formula is fine to a dozen digits at any sample size you will ever run — but the moment the quantity is offset from zero it degrades quadratically in mean/sd, and past mean/sd ≈ 10⁶ it returns a negative variance, which sqrt() turns into NaN. Simulations routinely produce such quantities: revenue in cents rather than fractions, a bidder’s cumulative bankroll, a timestamp, latency in nanoseconds. Since Welford costs one multiply and removes the whole class of failure, use it and stop thinking about it. And treat a negative variance as the diagnostic signature of naive-formula cancellation — it cannot arise from Welford. Write out the update: with d = x − mean_old the new mean is mean_old + d/n, so x − mean_new = d(1 − 1/n) and the increment is d·d(1 − 1/n) = d²(1 − 1/n) ≥ 0 for every n ≥ 1. Every term added to m2 is non-negative by construction, so a negative m2 is impossible short of a different bug.

The three numbers, and what each is for. The standard deviation describes one auction: how much revenue varies from auction to auction. It does not shrink with more trials, ever. The standard error of the mean (sem = sd/√n) describes your estimate: how far the sample mean is likely to be from the true mean. It shrinks as 1/√n. The 95 % confidence interval mean ± 1.96·sem is the reportable object. Confusing sd with sem is the most common error in simulation write-ups and it always makes the result look more certain than it is.

One run can rank two mechanisms backwards

Here is the demonstration, built so that the correct answer is known exactly and a wrong answer is unambiguously wrong.

Second-price auction, five bidders, values i.i.d. uniform [0,1], everyone truthful. Compare reserve price r = 0.40 against r = 0.50. Expected revenue has a closed form:

E[R | r] = n·rⁿ·(1−r) + (n−1)(1−rⁿ) − [n(n−1)/(n+1)]·(1−r^{n+1})

Term by term: n·rⁿ·(1−r) is the case where exactly one bidder clears the reserve and pays exactly r; the remaining two terms are E[second order statistic · 1{at least two clear r}]. Setting r = 0 recovers (n−1)/(n+1), the standard result. Evaluating in exact rational arithmetic with Python’s fractions.Fraction:

rE[R] exactdecimal
0.402094/31250.67008000
0.5043/640.67187500

r = 0.50 is better, by exactly 359/200000 = 0.001795. (It should be: for values uniform on [0,1] the Myerson-optimal reserve is 1/2 regardless of n — see Reserve Prices and Optimal Auctions.)

Now run the simulator at 1,000 trials per configuration, giving each configuration its own random stream — the natural, obvious, wrong design. Seed 2601:

INDEPENDENT STREAMS  seed=2601 trials=1000 n=5
  reserve 0.40 : 0.683568  sd 0.165385  95%CI [0.673317, 0.693819]
  reserve 0.50 : 0.658345  sd 0.199072  95%CI [0.646006, 0.670683]
  difference  : -0.025223  sd 0.257200  sem 0.008133  t -3.101  95%CI [-0.041165, -0.009282]

Read what that says. The lower reserve appears to earn 2.5 percentage points more, at t = −3.10. The two confidence intervals barely overlap. The interval on the difference, [−0.0412, −0.0093], does not contain the true value of +0.0018 — it does not even contain zero. Every conventional test declares the result significant, and the result is the opposite of the truth. A paper reporting this run would state, with a p-value under 0.002, that lowering the reserve raises revenue.

This is not a rare seed. Scanning seeds 1 … 1000 at several trial budgets and asking, for each, whether the run ranks the two reserves backwards:

trials/configbackwards“significant” at abs(t) > 1.96most misleading t
10048.80 %1.70 %−2.63
1,00042.70 %1.40 %−2.93
10,00022.90 %0.40 %−2.56
100,0001.40 %0.00 %−0.76
xychart-beta
    title "Fraction of seeds that rank r=0.40 above r=0.50 (truth: 0.50 wins)"
    x-axis ["100", "1000", "10000", "100000"]
    y-axis "% of 1000 seeds ranking backwards" 0 --> 55
    line [48.80, 42.70, 22.90, 1.40]

Backwards-ranking rate against trials per configuration, independent streams, 1,000 seeds each. The insight: at 100 trials the simulator is a coin flip — 48.8 % is indistinguishable from 50 % — and it takes 100,000 trials before the answer is reliably right. The effect being measured is 0.0018 against a per-auction standard deviation of 0.18, a ratio of 1:100, and no amount of clever coding compensates for that; only replications do.

Common random numbers: the same experiment, seventeen times cheaper

The fix is not more trials, it is the same trials. Draw the value vector once and evaluate both reserve prices on it. Every source of randomness that is not the thing under test cancels:

flowchart LR
    subgraph BAD["Independent streams — the obvious design"]
        S1["seed 2601<br/>stream 1"] --> V1["values V"] --> C1["config A<br/>r = 0.40"]
        S2["seed 2601<br/>stream 2"] --> V2["values W"] --> C2["config B<br/>r = 0.50"]
        C1 --> DIFF1["difference<br/>sd 0.2572"]
        C2 --> DIFF1
    end
    subgraph GOOD["Common random numbers"]
        S3["seed 2601<br/>stream 1"] --> V3["values V"]
        V3 --> C3["config A<br/>r = 0.40"]
        V3 --> C4["config B<br/>r = 0.50"]
        C3 --> DIFF2["difference<br/><b>sd 0.0657</b>"]
        C4 --> DIFF2
    end

The one structural change that buys a 17× variance reduction. The insight: in the top design the difference inherits two independent sources of draw noise; in the bottom design the value vector is bit-identical in both arms, so every fluctuation that is not caused by the reserve price cancels exactly in the subtraction. The per-configuration standard deviations are unchanged (0.165 and 0.177 either way) — it is only the standard deviation of the difference that collapses, and the difference is the thing you are trying to measure.

for (long t = 0; t < trials; t++) {
    for (int i = 0; i < n; i++) v[i] = d_draw(&vd, &ra);   /* ONE draw */
    double x = sp_rev(v, n, r1);                          /* both configs see it */
    double y = sp_rev(v, n, r2);
    acc_add(&a1, x); acc_add(&a2, y); acc_add(&ad, y - x); /* accumulate the DIFF */
}

This is common random numbers (CRN) — the auction analogue of the paired comparison in Running a Strategy Tournament in C, but stronger, because there the two strategies still faced different stochastic opponents whereas here the two configurations face a bit-identical world. The same seed 2601:

COMMON RANDOM NUMBERS  seed=2601 trials=1000 n=5
  reserve 0.40 : 0.683568  sd 0.165385  95%CI [0.673317, 0.693819]
  reserve 0.50 : 0.684641  sd 0.176713  95%CI [0.673688, 0.695594]
  difference  : +0.001073  sd 0.065703  sem 0.002078  t +0.517  95%CI [-0.002999, +0.005146]

The sign is now correct, the interval contains the true 0.001795, and the t of 0.517 says honestly “this run cannot resolve the difference.” That is the right answer at this budget. Note the mechanism: the two per-configuration standard deviations barely changed (0.165 and 0.177); it is the standard deviation of the difference that collapsed, from 0.2572 to 0.0657 — a 3.9× reduction in sd, 15.3× in variance. At seed 1 the same comparison gives 0.2566 → 0.0617, 4.16× in sd and 17.3× in variance. Roughly seventeen times fewer trials for the same precision, for a change that deletes a line of code.

Pushing the same seed to 100,000 trials:

INDEPENDENT   diff +0.001370  sem 0.000807  t +1.698  95%CI [-0.000211, +0.002952]
CRN           diff +0.001578  sem 0.000214  t +7.385  95%CI [+0.001159, +0.001997]

The CRN interval [0.001159, 0.001997] brackets the exact truth 0.001795 and is 3.8× narrower. The independent-stream run, at the same cost, still cannot exclude zero.

Repeating the seed scan with CRN:

trials/configbackwards (independent)backwards (CRN)“significant” and wrong (CRN)
10048.80 %37.80 %0.20 %
1,00042.70 %21.70 %0.10 %
10,00022.90 %0.40 %0.00 %
100,0001.40 %0.00 %0.00 %

CRN reaches at 10,000 trials the reliability that independent streams do not reach at 100,000.

The technique has a literature, and it names the failure mode this simulator has

CRN is not a trick invented for this note; it is the oldest and most reliable variance-reduction technique in discrete-event simulation, and the practitioner’s account of it is L’Ecuyer’s WSC 2015 tutorial “Random Number Generation with Multiple Streams for Sequential and Parallel Computing”. Two things in it are worth importing wholesale.

The magnitude is not special to auctions. L’Ecuyer’s worked example is an (s, S) inventory policy comparison — nothing like an auction — over n = 500 replications of m = 2000 days, comparing policies (80, 198) and (80, 200):

independent random numbers : mean diff 0.266   sd 1.530   90% CI (0.230, 0.302)
common random numbers      : mean diff 0.315   sd 0.352   90% CI (0.307, 0.324)
                             -> empirical variance of the difference divided by 18.85

A variance reduction factor of 18.85, against the 15.3–17.3 measured above for the reserve-price comparison. That two unrelated models land in the same range is the useful fact: for “compare two nearly identical configurations” problems, expect roughly an order of magnitude, and treat a much smaller factor as evidence that your streams are not actually synchronized. L’Ecuyer adds that the factor “can be arbitrarily large in some settings.”

The design rule, stated properly. The rule this note arrived at empirically — “give each configuration its own stream via PCG’s seq” — is half of the standard construction. L’Ecuyer’s full version is:

Use a separate stream for each source of random numbers needed in the system, and one different substream for each replication, resetting all streams to the start of the same substream before running each configuration.

Two streams because the purposes must not mix; substreams because the replications must line up. His own example shows exactly what goes wrong without the first half: in the inventory model one random number decides whether an order arrives, and orders are made on different days under different policies, so “a random number used to decide if the order arrives in one case could end up being used to generate a demand in the other case. This would greatly diminish the power of the CRN technology.”

That is precisely the failure mode listed below as “different configurations consuming different numbers of random draws” — and the substream discipline is the fix this note was missing. The reserve-price experiment survives without it only because both arms draw exactly n values per auction regardless of the reserve; add a randomizing bidder, an entry decision, or a stochastic number of bidders and the alignment breaks silently. Concretely, in this simulator that means: draw values from seq = 0, draw any bidder-side randomness from seq = 1, and before each trial reseed both to a substream index derived from the trial number rather than letting them run on. The cost is one extra rng_seed per trial; the benefit is that CRN keeps working when the model grows.

L’Ecuyer and Nadeau-Chamard’s follow-up (WSC 2021) covers the modern alternative: counter-based and splittable generators, where the state transition is trivial and all the work happens in the output function, so a substream is just an index rather than a jump-ahead polynomial. If this simulator ever needs streams it cannot enumerate in advance — an adaptive-agent model that spawns randomness on demand — that is the family to move to. PCG’s seq gives 2⁶³ streams but no substream mechanism, which is its one real limitation here.

…and then distrust the paired result too

CRN is variance reduction, not immunity. Scanning for the worst CRN seed at 1,000 trials finds seed 154:

COMMON RANDOM NUMBERS  seed=154 trials=1000 n=5
  reserve 0.40 : 0.670809 ...
  reserve 0.50 : 0.663206 ...
  difference  : -0.007603  sd 0.086361  sem 0.002731  t -2.784  95%CI [-0.012955, -0.002250]

Wrong sign, t = −2.78, interval excluding the truth. It is rarer — 1 seed in 1,000 rather than 14 — but it exists, and it will exist for exactly the reason it is supposed to: a 5 %-level test on a true effect too small to detect produces confident errors at roughly its nominal rate. The discipline that follows is the same as the tournament note’s: quote the effect size beside the t, and decide in advance what difference would matter. Here the whole revenue range across reserves from 0 to 0.8 is 0.667 → 0.672 → 0.557, a span of 0.115. An effect of 0.0018 is 1.6 % of that span. If your decision would not change over 1.6 %, do not spend 100,000 trials resolving it.

Cross-Checking Against Closed Forms

Simulation validates against arithmetic, never against plausibility. Four independent closed-form checks, all computed with fractions.Fraction so that a mismatch cannot be blamed on rounding:

CheckClosed formExactSimulated (200k–500k trials)
SP/FP/AP revenue, uniform, n=5(n−1)/(n+1)2/3 = 0.6666670.666666 / 0.666722 / 0.666835
SP revenue, uniform, r=0.5, n=5see formula above43/64 = 0.6718750.671935 ± 0.000261
FP eq. revenue, power F(x)=x², n=5k/(k+1) · αn/(αn+1), k=α(n−1)80/99 = 0.8080810.807914 ± 0.000325
SP revenue, power F(x)=x², n=5αn(n−1)·[1/(α(n−1)+1) − 1/(αn+1)]80/99 = 0.8080810.807985 ± 0.000513

The power-distribution pair is the strongest test in the set, because it exercises a different inverse-CDF path, a different equilibrium bid function, and a different order statistic, and still lands on the same rational number 80/99. It also demonstrates revenue equivalence away from the uniform case, which is where most textbook exercises stop.

A fifth check catches a class of bug the others miss. Run first price with a deliberately wrong bid function — a naive fixed 90 % shade instead of the equilibrium α(n−1)/(α(n−1)+1) = 8/9 ≈ 88.9 %:

power a=2.00 n=5 trials=200000
  FP equilibrium bid  revenue 0.807914 +- 0.000325
  SP truthful         revenue 0.807985 +- 0.000513
  FP naive 0.9*v      revenue 0.818013 +- 0.000329  <- WRONG bid function

The naive bidders hand the seller 0.9 × E[max] = 0.9 × 10/11 = 9/11 = 0.818182, and the simulator reports 0.818013 ± 0.000329 — an interval containing 9/11. Revenue equivalence is a statement about equilibrium, and a simulator that reports “first price beats second price” has almost certainly got the bid function wrong rather than discovered something. This check — bidders off equilibrium raise revenue and lose surplus — is the fastest way to distinguish a mechanism finding from a coding error.

The reserve-price sweep is the sixth check, and it reproduces a theorem rather than a number:

reserveSP revenuesemefficiencyunsold
0.000.6667220.0002521.00000.00 %
0.200.6669540.0002510.99990.03 %
0.400.6701390.0002470.99591.03 %
0.500.6719350.0002610.98433.14 %
0.600.6667070.0003150.95337.79 %
0.700.6388000.0004190.882516.79 %
0.800.5574290.0005530.737932.77 %

Revenue peaks at exactly the theoretical optimum r = 0.5 (0.671935 against the exact 0.671875), while efficiency falls monotonically from 1.0000 to 0.7379 and the item goes unsold a third of the time. That is Myerson’s result made visible: the revenue-maximizing auction is not the efficient one, and the price of the extra 0.8 % revenue is destroying 1.6 % of the surplus. Note also that revenue is remarkably flat near the optimum — 0.4, 0.5 and 0.6 differ by less than 1 % — which is precisely why the reversal experiment above is hard, and precisely why real reserve-price decisions need field experiments rather than arguments. Reserve Prices and Optimal Auctions and the Ostrovsky–Schwarz experiment discussed in Ad Auctions and GSP are the deployed version of this table.

Cross-Checking Against Reality — and Why the Simulator Fails That Test

Everything above validates the simulator against arithmetic. Arithmetic only certifies that the code computes what the model says. The separate question is whether the model describes anything, and for auctions there is thirty years of laboratory evidence to check it against. The standing survey is Kagel & Levin, “Auctions: A Survey of Experimental Research, 1995–2008” (Ohio State), the sequel to the 1995 survey in the Handbook of Experimental Economics. Its findings are unkind to every number this note has produced.

Revenue equivalence does not hold with human bidders. Kagel & Levin open by recording that as of the 1995 survey “it was clear that both the revenue equivalence theorem as well as the strategic equivalence between each of the two pair of auction formats failed”, with “persistent reports of significant bidding above the risk neutral Nash equilibrium (RNNE) benchmark in first-price sealed bid auctions.” The closed-form table above shows first-price and second-price revenue agreeing on 80/99 to five decimals; in the laboratory first-price revenue is systematically higher, and thirty years of papers have argued about whether the cause is risk aversion, regret, or quantal-response noise without settling it.

Truthful bidding in a second-price auction — the assumption this whole simulator rests on — is the exception, not the rule. The survey’s numbers, counting any bid within five cents of the induced value as sincere:

Studysincereoverbidunderbidsetting
Kagel & Levin (1993)27.0 %67.2 %5.7 %student subjects
Garratt, Walker & Wooders (2004)21.2 %37.5 %41.3 %experienced eBay bidders
Andreoni, Che & Kim (2007)77.3 % (85.5 % in the last 10 periods)4 bidders, uniform values; highest sincere rate the survey is aware of

Even the best-case study leaves nearly a quarter of bids off-equilibrium, and the eBay-veteran study is worse than the student baseline on sincerity — its subjects are simply split, with habitual sellers underbidding (50.9 % of bids) and buyers overbidding (45.5 %). Kagel & Levin’s reading is that field experience transfers only when the field closely resembles the lab: “there is no particular reason to think that experienced professionals will perform much better than student subjects when placed in a laboratory type setting.”

The efficiency column in the reserve-price table is the claim that breaks hardest. This simulator reports efficiency 1.0000 at r = 0 — with truthful bidders the highest-value bidder wins by construction, always. Shogren, Parkhurst & McIntosh (2006), running 20 second-price auctions with 10 bidders and values i.i.d. uniform on [0, 20], report that the highest-value bidder won 42.5 % of standard second-price auctions (55.0 % under a tournament payoff structure), and that the highest or second-highest value bidder won 70.0 %. Their mean deviation bid − value was 6.28 with a standard deviation of 63.51 — on a [0, 20] support, an average overbid of nearly a third of the whole range, with a standard deviation three times the range itself.

flowchart LR
    M["Model<br/>truthful bidders<br/>i.i.d. values, known F, known n"]
    M --> C["Closed form<br/>(n−1)/(n+1), 43/64, 80/99"]
    M --> S["This simulator<br/>Monte Carlo + CRN"]
    C <-->|"agree to 4–5 decimals<br/><b>validates the CODE</b>"| S
    L["Laboratory<br/>Kagel and Levin survey<br/>21–77% sincere bidding<br/>42.5% efficiency"]
    S -.->|"efficiency 1.0000 vs 42.5%<br/>revenue equivalence holds vs fails<br/><b>invalidates the MODEL</b>"| L

What it shows: the two independent validation axes and what each can and cannot certify. The insight: passing every closed-form check in the previous section says nothing whatsoever about whether the simulator’s conclusions transfer. Agreement with arithmetic is a test of the implementation; disagreement with the laboratory is a verdict on the assumptions, and no amount of extra trials fixes it.

The practical discipline this imposes is not “stop simulating” — it is to state the behavioural assumption as a parameter rather than as a fact, and to sweep it. The simulator already has the hook: the bidder is a function pointer. The deliberately-wrong 0.9 × v bid function used above as a bug detector is, read differently, the cheapest possible model of an off-equilibrium bidder — and it moved revenue from 0.808 to 0.818, a 1.2 % effect from a 1.2 % misbehaviour. Sweeping a shade multiplier, or drawing per-bidder shading from a distribution fitted to the lab’s bid − value spread, converts “here is the equilibrium answer” into “here is how fast the answer degrades as bidders deviate” — which is the question a mechanism designer actually has.

Reproducibility Discipline

The tournament note establishes the standard — the same command must produce the same bytes under every compiler — and the argument is not repeated here. What is worth recording is that the discipline holds for this simulator too, including through the transcendental functions, which is not obvious a priori:

--- determinism, mode=paired 1 20000 5
auc        79a221b1d768defd06bebbc69b2b7a5b     gcc -O2
auc_O0     79a221b1d768defd06bebbc69b2b7a5b     gcc -O0
auc_fast   79a221b1d768defd06bebbc69b2b7a5b     gcc -O3 -march=native -ffast-math
auc_clang  79a221b1d768defd06bebbc69b2b7a5b     clang -O2
--- determinism, mode=dist   (exercises log, exp, pow, sqrt)
auc / auc_O0 / auc_fast / auc_clang   c7977349065b5ee01797eb404a19b944
--- determinism, mode=resv2  (the reversal experiment)
auc / auc_O0 / auc_fast / auc_clang   141d0fee40780d664ccf87567ae1e81d

Byte-identical across four builds (gcc 16.1.1 and clang 22.1.8), including -ffast-math -march=native, and including the mode that calls log, exp, pow and sqrt two million times each.

Uncertain — this is weaker evidence than it looks

Verify: whether the simulator is bit-reproducible across libm implementations, not just across compilers. Reason: all four builds link the same glibc libm, so the test above proves the compiler does not reassociate the arithmetic — it does not prove that pow, log and exp return identical bits on musl, on macOS, or on a different glibc version. Only sqrt is required by IEEE-754 to be correctly rounded; log, exp and pow are not, and glibc’s accuracy tables list them as within a few ulp rather than exact. A build with -m32 (which historically changes x87 vs SSE evaluation) could not be produced here — the 32-bit toolchain is not installed. To resolve: run the same binaries’ output against a musl-linked build, or replace pow(u, 1/α) with exp(log(u)/α) and check whether the digests move. Until then, “reproducible across compilers on this machine” is the claim, not “reproducible everywhere.” #uncertain

The practical mitigation, if bit-exactness across platforms matters: keep transcendentals out of anything that feeds a decision. In this simulator the only decisions are comparisons of bids (> on doubles produced by multiplication) and the reserve test; log/pow feed only the value draws. A drifting last-ulp in pow changes a value by 10⁻¹⁶ and changes an ordering only if two values were already within 10⁻¹⁶, which happens with probability around 10⁻¹⁶ per comparison.

Failure Modes and Gotchas

Using rand(). Covered above: the standard’s own remedy is to vendor a generator. The secondary hazard is subtler — rand() % k is biased whenever k does not divide RAND_MAX+1, and (double)rand()/RAND_MAX produces a distribution on [0,1] with only 31 bits of resolution and a nonzero probability of returning exactly 1.0, which breaks 1/(1−U) style transforms.

Generating values inside the format driver. This is the design mistake that costs a factor of seventeen. If f_first calls the RNG, you can never feed two formats the same world, and CRN is unavailable. Draw first, pass arrays.

Different configurations consuming different numbers of random draws. CRN works only if the streams stay aligned. A bid function that calls the RNG (a randomizing bidder, or a normal draw with a spare-value cache in a different state) will desynchronize the two arms after the first auction where they differ, and CRN silently degrades to independent sampling — with no error message, just a wider interval than you expected. If a configuration must consume randomness, give it its own stream via PCG’s seq parameter rather than sharing the stream you are trying to hold fixed — and, per L’Ecuyer’s rule, reset every stream to the same substream index at the start of each replication so that a length mismatch in trial k cannot leak into trial k+1. Without the reset, one desynchronizing auction desynchronizes the entire remainder of the run. The diagnostic: compute the variance reduction factor and compare it against the ~15–19× that a properly synchronized paired comparison delivers; a factor near 1 means the streams have come apart.

Reporting sd where sem belongs. mean ± sd is not a confidence interval and is roughly √n times too wide. Reporting mean ± sem without the 1.96 is a 68 % interval labelled as if it were 95 %.

Assuming a finite variance. Every CI in this note is a normal approximation resting on the central limit theorem. With Pareto-distributed values at shape α ≤ 2 the variance is infinite, the CLT does not apply, and the interval is meaningless — while looking completely normal on screen. The symptom is an sd that grows as you add trials instead of stabilizing. Check it explicitly: plot sd against n, not just the mean.

Comparing across seeds instead of within. Running format A at seed 1 and format B at seed 2 and comparing the means is the independent-stream design that ranks backwards 42.7 % of the time at 1,000 trials.

Ignoring the multiple-comparisons problem. A sweep over nine reserve prices is 36 pairwise comparisons. At a 5 % level you should expect roughly two “significant” results from a table where nothing differs at all.

Catastrophic cancellation in the accumulator. Naive Σx² bookkeeping. Use Welford.

log(0). Guard the uniform draw for every inverse transform involving a logarithm or a negative power.

M_PI is not in <math.h> under a strict feature-test macro. Compiling with -std=gnu11 and #define _POSIX_C_SOURCE 200809L suppresses _DEFAULT_SOURCE, and M_PI — which is an X/Open extension, not ISO C — disappears. This cost a build cycle here. Either define the constant yourself or do not narrow the feature-test macro. (The related trap that -std=c11 hides clock_gettime is documented in Running a Strategy Tournament in C.)

Sorting bids with qsort. Unnecessary (O(n log n) for two order statistics) and a tie-breaking reproducibility risk. Write the two-pass top2.

Reserve-price arithmetic in second price. The price is max(second bid, reserve) conditional on the top bid clearing the reserve. Two other plausible-looking formulas both give monotone-in-r revenue curves that hide the optimum.

Alternatives and When to Choose Them

ApproachGood forCost
Closed-form / exact rational (fractions.Fraction)small instances, verifying the simulator, Ad Auctions and GSP-style counterexamplesonly works where a formula exists
Numerical integration over order statisticsexpected revenue for i.i.d. values, any F, any nneeds a quadrature routine; no equilibrium dynamics
Monte Carlo (this note)anything: asymmetric bidders, budgets, correlated values, learning agentsstatistical error, and all the ways to get it wrong above
Antithetic variateshalving variance for monotone outcomesbreaks when the outcome is not monotone in U
Field experimentthe only way to learn what bidders actually dosee Ostrovsky–Schwarz: 460,000 keywords for a 4.9 % effect

Use exact arithmetic wherever it reaches. The whole GSP counterexample in Ad Auctions and GSP — an equilibrium verified by exhaustive best-response search — is twenty lines of Fraction and needs no simulator at all, and being exact is what makes it a proof rather than evidence. Reach for Monte Carlo when the state space is continuous or the agents adapt.

Antithetic variates deserve a mention as the obvious next variance-reduction technique and a warning about why they are not used here. The idea: for each draw U, also evaluate at 1−U; if the outcome is monotone in U, the two are negatively correlated and their average has lower variance. Auction revenue is monotone in each individual value, but with n bidders the antithetic vector (1−U₁, …, 1−Uₙ) reverses the ranking, and the second-highest order statistic is not a monotone function of the vector in the required sense. CRN is the technique that works cleanly here because it holds the world fixed rather than trying to mirror it.

A note on OpenSpiel, Gambit and nashpy. As recorded in Games and Strategic Systems in C MOC, none is installable on this machine (no pip, no numpy), and the substitute is not a library at all: write the twenty-line exact oracle in Python’s stdlib and make the C agree with it. Every closed-form column in this note came from that oracle.

Production Notes

Simulation is how auction changes actually get shipped, and it is not enough. The clearest documented case is the Ostrovsky–Schwarz reserve-price work discussed at length in Ad Auctions and GSP: the theory gave a per-keyword optimal reserve, simulation gave the expected revenue lift, and the field experiment then found that the effect was +4.88 % on frequently-searched keywords with high optimal reserves and −8.73 % on frequently-searched keywords with low ones — a sign flip that no simulation predicted, because its cause was that advertisers on low-volume keywords “simply do not spend much time optimizing bids.” The simulator assumed equilibrium; the market had not reached one.

Their statistical caution is worth copying verbatim into any simulation report. The naive headline was “+12.85 %”; the authors immediately note that “by excluding a single keyword from the control sample, this number can be reduced to around 8 %”, and elsewhere that removing two keywords moves a +10.3 % estimate to +0.19 %. Their fix was to change the outcome variable to one that outliers cannot dominate (revenue per search rather than revenue). The simulator equivalent: when a mean is driven by a heavy tail, report the median and a trimmed mean beside it, and state what happens when you drop the largest few observations. A result that moves by a third when you delete one data point is not a result.

Varian’s fit is the other production pattern: rather than simulating forward from assumed values, he took 2,425 real Google auctions and solved, per auction, for the minimal perturbation to ad quality that would make the observed bids satisfy the symmetric-Nash inequalities — finding a mean absolute deviation of 5.8 %. That is a simulator run backwards, and it is a much stronger test of a behavioural model than any forward simulation, because the data did not come from your own assumptions.

Throughput, for budgeting. On this machine (heavily loaded, /proc/loadavg 8.8), the reversal experiment at 100,000 trials × 1,000 seeds × 2 configurations is 2 × 10⁸ auctions of 5 bidders — a few minutes single-threaded. At 1.69 ns per uniform and roughly 7 draws plus a top2 per auction, a five-bidder auction costs on the order of tens of nanoseconds. The compute is free; the thinking about what to compute is not. Every wrong answer in this note came from a design choice, never from a shortage of cycles.

See Also