Problem 323: Bitwise-OR Operations on Random Integers

Detailed mathematical approach

Problem Summary

Let \(y_0,y_1,y_2,\dots\) be independent random unsigned 32-bit integers, each of the values \(0\le y_i\lt 2^{32}\) equally likely, and put \(x_0=0\) and \(x_i=x_{i-1}\mathbin{|}y_{i-1}\) for \(i\gt 0\), with \(\mathbin{|}\) the bitwise OR. A bit once set is never cleared, and with probability \(1\) there is a first index \(N\) with \(x_N=2^{32}-1\), the pattern of all ones; from then on \(x_i=2^{32}-1\) for every \(i\ge N\). We must find the expected value of \(N\) to ten decimal places. ^[1] The statement gives no checkpoint values, so the small cases below are computed from scratch.

The whole solution in one sentence: \(N\) is the largest of \(32\) independent waiting times, one per bit, each of them the number of coin tosses until the first head, so \(N\le n\) has probability \((1-2^{-n})^{32}\), and adding the complementary probabilities over all \(n\ge 0\) gives the expectation as a short series of positive terms.

Mathematical Approach

1. What the recursion does

Write \(b\) for the word length: \(b=32\) in the problem, \(b=1,2,3,4\) in the small cases below. The OR never turns a one back into a zero, so bit \(j\) of \(x_i\) is one exactly when at least one of the draws \(y_0,\dots,y_{i-1}\) had bit \(j\) set. Hence \(x_i\) is all ones as soon as every bit has been seen once, and stays there; the \(N\) of the statement is the first index at which the pattern is complete,

$$N=\min\{\,i:\ x_i=2^{b}-1\,\}.$$

Here is a complete run with \(b=4\). The draws \(y_0=9\), \(y_1=4\), \(y_2=2\) have the bit patterns \(1001\), \(0100\), \(0010\), and the ORs accumulate as \(x_1=9=(1001)_2\), then \(x_2=13=(1101)_2\), then \(x_3=15=(1111)_2\). For this history \(N=3\). Note that \(N\ge 1\) always, because \(x_0=0\).

2. One bit at a time

Fix a bit position \(j\) and let \(G_j\) be the index of the first draw that sets it: the smallest \(i\ge 1\) with bit \(j\) set in \(y_{i-1}\). A uniformly random value in \(\{0,\dots,2^{b}-1\}\) is the same thing as \(b\) independent fair coin tosses, one per bit, since the uniform distribution on a product set is the product of its factors' uniform distributions. So \(G_j\) counts tosses until the first head: the geometric distribution on \(\{1,2,3,\dots\}\) with \(p=\tfrac12\), of mean \(1/p=2\) and distribution function \(1-(1-p)^{n}\). ^[2] In particular

$$\Pr(G_j\gt n)=2^{-n},\qquad n\ge 0.$$

By Section 1 the pattern is complete exactly when every bit has appeared, so \(N=\max_j G_j\). The case \(b=1\) is already finished: \(N=G_1\), of expectation \(2\).

3. Why the bits may be multiplied

The same argument applied to each draw makes the \(bn\) bits of \(y_0,\dots,y_{n-1}\) a table of independent fair coin tosses. Whether bit \(j\) is still zero in \(x_n\) depends only on the \(n\) tosses in column \(j\), and different columns use disjoint tosses, so the \(b\) events \(\{G_j\gt n\}\) are independent, and so are their complements.

4. The distribution of \(N\)

Independence turns the maximum into a product:

$$\Pr(N\le n)=\prod_{j=1}^{b}\Pr(G_j\le n)=\bigl(1-2^{-n}\bigr)^{b}.$$

At \(n=0\) this gives \(0\), as it must, since \(x_0=0\). A check by hand with \(b=2\), \(n=2\): the pair \((y_0,y_1)\) ranges over the \(16\) elements of \(\{0,1,2,3\}^2\), and \(x_2=y_0\mathbin{|}y_1\) equals \(3\) for \(16-4-4+1=9\) of them — drop the \(4\) pairs in which neither draw has bit \(0\) and the \(4\) in which neither has bit \(1\), then add back \((0,0)\), counted twice. So \(\Pr(N\le 2)=9/16\), which is \((1-2^{-2})^2\).

5. From the distribution to the mean

For a random variable \(N\) with values in \(\{1,2,3,\dots\}\) the expectation is the sum of its tail probabilities, ^[3]

$$\mathbb E[N]=\sum_{n\ge 0}\Pr(N\gt n).$$

The proof is one exchange of summation order: substituting \(\Pr(N\gt n)=\sum_{m\gt n}\Pr(N=m)\) and collecting the pairs \((n,m)\) with \(0\le n\lt m\) by \(m\) turns the double sum into \(\sum_{m\ge 1}m\Pr(N=m)\), the definition of the expectation; all terms are non-negative, so the rearrangement is legitimate. With Section 4 this gives the formula the programs evaluate,

$$\mathbb E[N]=\sum_{n\ge 0}\Bigl[1-\bigl(1-2^{-n}\bigr)^{b}\Bigr].$$

The first terms for \(b=32\) are \(1\), then \(0.999999999767\), \(0.999899547574\), \(0.986060162962\), \(0.873211213623\) and \(0.637944710744\). They start close to \(1\) — after three draws the pattern is still incomplete in \(98.6\%\) of runs — and then fall away geometrically.

6. Worked example: two bits, end to end

For \(b=2\) the distribution function \((1-2^{-n})^2\) is \(1/4\), \(9/16\), \(49/64\), \(225/256\), \(961/1024\) at \(n=1,\dots,5\); differencing gives the point probabilities \(1/4\), \(5/16\), \(13/64\), \(29/256\), \(61/1024\). Summing \(n\Pr(N=n)\) over \(n\le 40\) gives \(2.666666666590\), creeping up on \(8/3\). The tail series settles it exactly, since \(1-(1-2^{-n})^2=2\cdot 2^{-n}-4^{-n}\) and two geometric series give

$$\begin{aligned}\mathbb E[N]&=2\cdot\frac{1}{1-\frac12}-\frac{1}{1-\frac14}\\&=4-\frac43=\frac83.\end{aligned}$$

So two bits cost \(8/3\approx 2.666666666667\) draws on average. The same calculation for \(b=3\) and \(b=4\) gives \(22/7\) and \(368/105\).

7. An exact closed form by inclusion–exclusion

The tail probability can also be expanded. Let \(A_j\) be the event that bit \(j\) is still zero after \(n\) draws, so \(\Pr(N\gt n)=\Pr(A_1\cup\dots\cup A_b)\). A given set of \(k\) bit positions stays all-zero exactly when all \(kn\) of its tosses are tails, with probability \(2^{-kn}\), depending on \(k\) alone. The inclusion–exclusion principle ^[4] therefore gives

$$\Pr(N\gt n)=\sum_{k=1}^{b}(-1)^{k+1}\binom{b}{k}2^{-kn},$$

the same expansion the binomial theorem gives for \(1-(1-2^{-n})^b\). Summing over \(n\ge 0\) term by term, legitimate because the inner sum is finite, and using \(\sum_{n\ge 0}2^{-kn}=2^{k}/(2^{k}-1)\), gives a finite exact formula:

$$\mathbb E[N]=\sum_{k=1}^{b}(-1)^{k+1}\binom{b}{k}\frac{2^{k}}{2^{k}-1}.$$

For \(b=4\) it reads \(8-8+\tfrac{32}{7}-\tfrac{16}{15}=\tfrac{368}{105}\), matching Section 6. As exact fractions the first values are \(2\), \(\tfrac83\), \(\tfrac{22}{7}\), \(\tfrac{368}{105}\), \(\tfrac{2470}{651}\), \(\tfrac{7880}{1953}\), \(\tfrac{150266}{35433}\) and \(\tfrac{13315424}{3011805}\) for \(b=1,\dots,8\), the last two being \(4.240848926142\) and \(4.421077725816\).

8. Why the closed form needs extra precision

The finite formula is ideal for exact rational checks, but in ordinary floating point its alternating terms cancel badly. Its largest single term is the one with \(k=16\),

$$\binom{32}{16}\frac{2^{16}}{2^{16}-1}\approx 601089561.9,$$

and the absolute values of all \(32\) terms add up to about \(4295400046\), just over \(2^{32}-1=4294967295\). These nine-digit numbers must cancel down to the value Section 10 keeps below \(7\): a sensitivity of order \(10^{9}\), so about nine significant decimal digits are lost. ^[5] Summed in the natural order in IEEE double arithmetic the formula comes out \(8.8\cdot 10^{-8}\) too high — seven correct decimals, where eleven significant digits are wanted.

The tail series of Section 5 cancels nothing between terms, since all of them are positive. Inside a term care is still needed, since a rounding of the base \(1-2^{-n}\) is amplified by the power: factorising \(a^{b}-c^{b}\) for \(a,c\in[0,1]\) gives \(b\) summands of modulus at most \(1\), so \(|a^{b}-c^{b}|\le b\,|a-c|\). A term thus carries an absolute error of order \(b\) machine roundings, at worst \(1.8\cdot 10^{-15}\) over all \(n\) in double arithmetic, and absolute error is what a sum of positive terms accumulates: the whole series in double comes out \(3.8\cdot 10^{-15}\) low.

9. How many terms are needed

The union bound gives a clean estimate of the dropped tail: the probability that at least one of several events happens is at most the sum of their probabilities, since a history in several of them is counted several times on the right. Each of the \(b\) events here is that one named bit is still zero after \(n\) draws, of probability \(2^{-n}\), so \(\Pr(N\gt n)\le b\,2^{-n}\) and

$$\sum_{n\gt n_0}\Pr(N\gt n)\le b\,2^{-n_0}.$$

Every omitted term is positive, so the exact partial sum through \(n_0\), and that sum plus \(b\,2^{-n_0}\), bracket \(\mathbb E[N]\); the bound is sharp to three digits for \(n_0\) past \(40\). As \(32\cdot 2^{-40}\approx 2.9\cdot 10^{-11}\), \(41\) exact terms bracket the value that tightly — though a correctly rounded tenth decimal needs that bracket, widened by the arithmetic error of the partial sum, to lie inside one rounding interval. The programs stop more generously, when a term falls below \(10^{-20}\): for \(b=32\) that first happens at \(n=72\), where the term is \(6.78\cdot 10^{-21}\) against \(1.36\cdot 10^{-20}\) at \(n=71\), so \(73\) passes suffice; at \(b=128\), the largest the C++ accepts, \(75\) do.

10. How large should the answer be?

Two bounds pin the size of the result without computing it. For the upper one, use \(\Pr(N\gt n)\le\min(1,b\,2^{-n})\) and split the series at \(m=\lceil\log_2 b\rceil\):

$$\mathbb E[N]\le m+\sum_{n\ge m}b\,2^{-n}=m+2b\,2^{-m}\le m+2.$$

For the lower one, \(1-x\le e^{-x}\) gives \((1-2^{-n})^{b}\le e^{-b2^{-n}}\), so each term is at least \(1-e^{-b2^{-n}}\). Keeping only \(n\le\lfloor\log_2 b\rfloor\) and writing \(j=\lfloor\log_2 b\rfloor-n\ge 0\), each kept term is at least \(1-e^{-2^{j}}\), and \(\sum_{j\ge 0}e^{-2^{j}}=0.52186593\ldots\lt 0.52187\), whence

$$\mathbb E[N]\ge\lfloor\log_2 b\rfloor+1-0.52187.$$

Both bounds hold for every \(b\ge 1\) and keep \(\mathbb E[N]\) within a bounded distance of \(\log_2 b\): for \(b=32\), between \(5.478\) and \(7\). Each doubling of the word length therefore costs roughly one more draw, and \(\mathbb E[N]\) is indeed \(5.3773779018\) at \(b=16\) and \(7.3439903348\) at \(b=64\), \(1.966612433\) apart.

Had each draw delivered one uniformly random bit position instead of a random subset, this would be the coupon collector's problem, of expected waiting time \(b\,H_b\) with \(H_b=\sum_{j=1}^{b}1/j\) ^[7] — \(32H_{32}\approx 129.8718463\) for \(32\) coupons, growing like \(b\ln b\) since \(H_b=\ln b+\gamma+O(1/b)\), with \(\gamma\) Euler's constant. ^[6] Receiving about half the missing bits at each draw brings that down to logarithmic growth.

How the Code Works

The C++ program ^[8] reads --bits= (default \(32\)) and --skip-checkpoints into an Options struct, rejecting unknown arguments and bit counts outside \(1\dots128\). One limitation: parse_int_after_prefix accumulates the digits in an int with no overflow check before the range test, so an argument above \(2^{31}-1\) overflows a signed integer, which is undefined behaviour — the default run never reaches it.

solve is Section 5 written out, with the running total in a long double. For \(n=0,1,2,\dots\) it sets bit_zero_prob to std::ldexp(1.0L, -n) — scaling the exponent yields \(2^{-n}\) with no rounding at all, which a general power function would not guarantee — then takes all_ones_prob to be std::powl(1 - bit_zero_prob, bits), the distribution function of Section 4, adds term, which is one minus it, and breaks once term drops below 1e-20. The exponent reaches std::powl as a long double, so the accuracy of the power is whatever the C library provides; the for header caps the loop at \(10000\) passes, far beyond the stopping indices of Section 9.

run_checkpoints re-derives the two worked examples, requiring solve(1) to be \(2\) and solve(2) to be \(8/3\) within \(10^{-15}\); main then prints the result with std::fixed and setprecision(10) after a cast to double.

The arithmetic is where the three languages genuinely differ. long double is not a fixed format: the x87 80-bit type on many x86 platforms, a wider type on some others, the format of double elsewhere. ^[9] On the machine used here its significand has 53 bits, so 1 - bit_zero_prob rounds to exactly \(1\) at \(n=54\) (at \(n=65\) for an 80-bit long double); the term becomes exactly zero and the loop leaves after \(55\) passes rather than \(73\). The dropped remainder is at most \(32\cdot 2^{-53}\approx 3.6\cdot 10^{-15}\), far below what the tenth decimal can feel, but the comparison fires only because the computed term has rounded to zero while the exact term at \(n=54\) is still \(1.78\cdot 10^{-15}\).

Python runs the identical loop with Decimal at getcontext().prec = 50. Decimal arithmetic takes its inputs as exact and rounds only results, to a precision the caller chooses, ^[10] so \(1-2^{-n}\) keeps \(50\) significant digits over the range that matters and the loop really runs to \(n=72\), \(73\) passes. Java hands MathContext.DECIMAL128 — \(34\) digits, round-half-even ^[10] — to its BigDecimal division and power, while its two subtract calls carry no context and are exact; it takes \(73\) passes too, and String.format(Locale.US, …) fixes the decimal separator. Neither port carries the checkpoints, so the self-check lives only in the C++ file.

Complexity and Verification

With exact terms the loop makes \(L=O(\log(b/\varepsilon))\) passes for a threshold \(\varepsilon\): \(73\) at \(b=32\), and \(55\) for the C++ on the 53-bit platform described above. The C++ makes one ldexp and one powl call per pass; Python and Java instead recompute both powers every time, with exponents \(n\) and \(b\), which under a repeated-squaring model is \(O(\log n+\log b)\) multiplications, so \(O(L(\log L+\log b))\) in all. These count arithmetic operations, whose cost and storage follow the fixed working precision. No version keeps a table, only a few numbers, and all three finish in milliseconds.

Besides the C++ checkpoints, the following were verified for this article, in exact rational or 50-digit arithmetic unless noted.

• For every \(b=1,\dots,16\), exact rational arithmetic puts the closed form of Section 7 above the tail series truncated at \(n=300\) by at most \(b\,2^{-300}\), as Section 9 predicts.
• A Markov chain on the number of set bits, by first-step analysis in exact fractions, gives the same values for \(b=1,\dots,16\). At \(k\lt b\) ones the chain stays put with probability \(2^{-(b-k)}\), so that state's equation is divided by \(1-2^{-(b-k)}\).
• A linear system over all \(2^{b}\) bit patterns, assuming neither independence nor any symmetry, returns \(2\), \(\tfrac83\), \(\tfrac{22}{7}\) and \(\tfrac{368}{105}\) for \(b=1,\dots,4\).
• Exhaustive enumeration of histories confirms \(\Pr(N\le n)=(1-2^{-n})^b\) exactly, for \(b=2\) up to \(n=5\) and \(b=3\) up to \(n=4\).
• A simulation of the real 32-bit process, \(10^{6}\) runs from one fixed seed, reproduced the value of the series to within \(0.001\).
• All three programs print the same ten decimals as each other and as the exact closed form rounded to ten places. The C++ through --bits=, and the Python solve through its bits parameter, also match the exact values for \(b=1,2,3,4,8,16,64\); the Java port fixes the word length at \(32\).

Two margins are worth noting. The C++ checkpoints compare \(b=1\) and \(b=2\) with \(2\) and \(8/3\) at a tolerance of \(10^{-15}\); on the 53-bit long double used here they pass with errors \(0\) and \(4.4\cdot 10^{-16}\), only about a factor of two inside it, and say nothing about other word lengths. In decimal arithmetic the tail series lands within \(2.4\cdot 10^{-13}\) of the exact value at \(14\) working digits and \(5.4\cdot 10^{-11}\) at \(12\), while the closed form first reaches \(1.7\cdot 10^{-11}\) at \(20\).

Footnotes and References

The references supply the background facts used above; everything else is derived explicitly here.

1. Project Euler 323 — Bitwise-OR Operations on Random Integers. The official statement draws independent uniformly random unsigned 32-bit integers, sets x₀ = 0 and x_i = x_(i−1) OR y_(i−1), and asks for the expected first index at which x_i reaches the all-ones pattern 2³² − 1, to ten decimal places. No checkpoint values are given.
2. Wikipedia — Geometric distribution. For the version counting trials up to and including the first success, it gives the expected value 1/p and the distribution function 1 − (1 − p)^⌊x⌋ for x ≥ 1; Section 2 uses both with p = 1/2.
3. Wikipedia — Expected value. Its tail-sum subsection states that for a random variable with values in the non-negative integers the expected value equals the sum of Pr(X > k) over k ≥ 0.
4. Wikipedia — Inclusion–exclusion principle. It gives the alternating formula for the size of a union of n sets and states that it stays valid with finite probabilities in place of cardinalities, the version Section 7 applies.
5. Wikipedia — Loss of significance. It describes catastrophic cancellation: subtracting good approximations of two nearby numbers can give a very bad approximation of their difference, an effect inherent to subtraction when the inputs are approximations.
6. Wikipedia — Euler's constant. It defines γ as the limiting difference between the harmonic series and the natural logarithm and gives γ = 0.5772156649…, the constant Section 10 uses in H_b = ln b + γ + O(1/b).
7. Wikipedia — Coupon collector's problem. Collecting one uniformly random coupon per draw takes n·H_n draws in expectation, asymptotically n log n + γn + 1/2, a growth rate the article gives as Θ(n log n).
8. C++ — Euler323.cpp. The linked immutable revision contains solve, run_checkpoints and the --bits= option. It is the implementation record for this article.
9. cppreference — Fundamental types. It describes long double as an extended-precision type that need not map to an IEEE-754 type: the x87 80-bit format on many x86 implementations, binary128 on some others, the format of double on MSVC.
10. Python — the decimal module documents arithmetic with a user-alterable precision in which values count as exact and only results are rounded, the setting getcontext().prec = 50 chooses. Java — MathContext documents DECIMAL128 as 34 digits with rounding mode HALF_EVEN.