Problem 307: Chip Defects

View on Project Euler

Project Euler Problem 307 Solution

EulerSolve provides an optimized solution for Project Euler Problem 307, Chip Defects, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We throw \(k\) defects independently and uniformly onto \(n\) chips. Each defect is placed separately, so there are exactly \(n^k\) equally likely outcomes. The goal is to find the probability that at least one chip receives three or more defects. Mathematical Approach 1) Work with the complement Let \(E\) be the event that some chip receives at least three defects. It is much easier to count the complementary event $$E^c=\{\text{every chip receives at most two defects}\}.$$ Then $$\Pr(E)=1-\Pr(E^c).$$ So the real task is to count all placements in which occupancies are only \(0\), \(1\), or \(2\). 2) Classify by the number of double-occupied chips Suppose exactly \(i\) chips receive two defects. Then exactly \(k-2i\) chips receive one defect, and the remaining \(n-(k-i)\) chips receive none. The parameter \(i\) can range only over $$0 \le i \le \left\lfloor \frac{k}{2}\right\rfloor.$$ 3) Counting a fixed \(i\)-layer We count placements with exactly \(i\) double chips in two steps. Step A: choose how the \(k\) labeled defects are grouped into \(i\) unordered pairs and \(k-2i\) singles. The number of such groupings is $$\frac{k!}{2^i\,i!\,(k-2i)!}.$$ The factor \(2^i\) removes the order inside each pair, and \(i!\) removes the order among the pairs themselves. Step B: send these \(k-i\) groups to distinct chips....

Detailed mathematical approach

Problem Summary

A factory makes \(n\) integrated-circuit chips and \(k\) defects fall on them at random: each lands on one chip, every chip is equally likely, the defects are independent, and a chip may carry any number. Let \(p(k,n)\) be the probability that some chip has at least three defects. The statement gives \(p(3,7)\approx0.0204081633\) and asks for \(p(20\,000,1\,000\,000)\) to ten decimals, as 0.abcdefghij. [1]

Independence and uniformity make all \(n^k\) ways of sending the \(k\) labelled defects to chips equally likely: the occupancy model of \(k\) balls in \(n\) bins. [3] The checkpoint confirms that reading, since three defects reach a load of three only when all three share a chip, which happens in \(7\) of the \(7^3=343\) assignments, and \(7/343=1/49\) matches the quoted value. Drawing uniformly from the \(\binom93=84\) load vectors would give \(7/84\approx0.0833\) instead, but that is a different model; forgetting the labels changes no probability.

The whole solution in one sentence: we group the arrangements that keep every chip under three defects by how many chips carry a pair, add those groups' probabilities with logarithms so that the arithmetic stays in range, and subtract the total from one.

Mathematical Approach

1. Why the complement is the easier event

The event "some chip has at least three defects" is a union of overlapping events, and one chip may hold anything from three defects to all \(k\), so inclusion–exclusion over the chips would be unwieldy. Its complement is one easily counted condition. Write an arrangement as a map from the \(k\) defects to the \(n\) chips, and let the load \(\ell_j\) of chip \(j\) be the number of defects sent to it, so \(\ell_1+\cdots+\ell_n=k\). Then

$$p(k,n)=1-\Pr\Bigl(\max_{j}\ell_j\le2\Bigr).$$

2. The shape of such an arrangement

Suppose every load is at most \(2\) and let \(i\) count the chips of load exactly \(2\). Counting defects gives \(2i+\#\{j:\ell_j=1\}=k\), so \(k-2i\) chips have load \(1\), the other \(n-k+i\) are clean, \(i\) runs over \(0\le i\le\lfloor k/2\rfloor\), and the \(k-i\) occupied chips cannot outnumber \(n\). This index cuts the arrangements into disjoint classes that exhaust them, so we may count one class at a time.

3. Counting one layer

Let \(N_i\) be the number of arrangements with exactly \(i\) chips of load \(2\) and \(k-2i\) of load \(1\). Write \((n)_m=n(n-1)\cdots(n-m+1)\) for the falling factorial, the number of injective maps from an \(m\)-element set into an \(n\)-element set. [4] Then

$$N_i=\frac{k!}{2^i\,i!\,(k-2i)!}\;(n)_{k-i}.$$

Proof. Such an arrangement is the same thing as a pair: a partition of the \(k\) defects into \(i\) blocks of size \(2\) and \(k-2i\) of size \(1\), plus an injective map from the blocks to the chips. From the arrangement, the blocks are the defect sets of the occupied chips and the map sends each to its chip, recovering the pair uniquely; conversely every such pair gives the required loads, and the map must be injective, since two blocks on one chip would either push its load past \(2\) or fuse two singletons into a pair the partition lacks.

The partitions number \(k!/(2^i\,i!\,(k-2i)!)\): list the defects in one of \(k!\) orders and read the first \(2i\) entries as \(i\) pairs, the rest as singletons. Two orders give the same partition exactly when they differ by swapping entries inside pairs (\(2^i\) ways), permuting the pairs (\(i!\)) or the singletons (\((k-2i)!\)), so each arises from \(2^i\,i!\,(k-2i)!\) orders; the blocks then admit \((n)_{k-i}\) placements. \(\square\)

No range condition is needed: if \(k-i\gt n\) then \((n)_{k-i}\) has the factor \(n-n=0\) and vanishes by itself. And \(k\ge2n+1\) forces \(k-i\ge n+1\) in every layer, so all \(N_i\) vanish and \(p(k,n)=1\) — the pigeonhole bound, reproduced automatically.

4. Worked example: four defects on three chips

Take \(k=4\) and \(n=3\), so there are \(3^4=81\) arrangements and \(i\) ranges over \(0,1,2\).

  • \(i=0\): four singletons would need four distinct chips, and \((3)_4=3\cdot2\cdot1\cdot0=0\), so \(N_0=0\).
  • \(i=1\): there are \(4!/(2\cdot1\cdot2)=6\) partitions, one per choice of the pair, and \((3)_3=6\) placements, so \(N_1=36\).
  • \(i=2\): there are \(4!/(4\cdot2\cdot1)=3\) partitions into two pairs and \((3)_2=6\) placements, so \(N_2=18\).

Hence \(N_0+N_1+N_2=54\) and \(p(4,3)=1-54/81=1/3\). A direct count of the bad arrangements confirms it: one chip holds all four in \(3\) ways, and one holds exactly three in \(3\cdot\binom43\cdot2=24\) ways — chip, three defects, another chip for the fourth; \(3+24=27=81/3\).

5. The checkpoint, and a closed form for three defects

The statement's checkpoint comes out as \(N_0=(7)_3=210\) and \(N_1=3\cdot(7)_2=126\), so \(336\) of the \(343\) arrangements are harmless and \(p(3,7)=7/343=1/49\). For general \(n\) the same two layers give \(N_0+N_1=(n)_3+3(n)_2=n^3-n\), so

$$p(3,n)=\frac{n^3-(n^3-n)}{n^3}=\frac1{n^2},$$

clear in hindsight: with three defects the only bad arrangements put all three on one chip, \(n\) of the \(n^3\). Being exact for every \(n\), it makes a good test case.

6. A cross-check by generating functions

Encode a chip that receives \(r\) defects by \(x^r/r!\), so a chip limited to \(r\le2\) carries the exponential generating function \(1+x+x^2/2\). The labelled product rule hands the \(k\) distinct labels out among the \(n\) chips and multiplies their generating functions — the loads are not independent, as they sum to \(k\) — giving [6]

$$\sum_{i}N_i=k!\,[x^k]\left(1+x+\frac{x^2}2\right)^{n},$$

with \([x^k]\) the coefficient of \(x^k\). Collecting the term with \(i\) factors \(x^2/2\) reproduces \(N_i\). We checked the identity for \((k,n)=(4,3)\), \((6,4)\), \((5,6)\), \((7,5)\): both sides give \(54\), \(1440\), \(6120\), \(25200\). The partition counts alone are \(k!\) times the coefficient of \(x^k\) in \(\exp(x+x^2/2)\), the involution numbers \(1,1,2,4,10,26,76,\ldots\) [5] — for \(k=4\), \(1+6+3=10\).

7. From counts to probabilities

Dividing by \(n^k\) turns the layer counts into probabilities \(P_i=N_i/n^k\). Pulling \(n^{k-i}\) out of the falling factorial, \((n)_{k-i}=n^{k-i}\prod_{j=0}^{k-i-1}(1-j/n)\), gives the numerically usable form

$$\begin{aligned}P_i&=\frac{k!}{2^i\,i!\,(k-2i)!\,n^i}\\&\quad{}\cdot\prod_{j=0}^{k-i-1}\left(1-\frac jn\right)\end{aligned}$$

and the answer is \(p(k,n)=1-\sum_{i=0}^{\lfloor k/2\rfloor}P_i\). The first layer is familiar: \(P_0=\prod_{j=0}^{k-1}(1-j/n)\) is the probability that all \(k\) defects land on different chips, the birthday-problem product with \(n\) in place of \(365\). [2] For \(k=20\,000\), \(n=10^6\) it is about \(3.6\cdot10^{-88}\): a collision is all but inevitable, and nearly all of the complement's probability lies in layers with many pairs.

8. Stepping from one term to the next

Consecutive layers share almost all their factors, so recomputing each \(P_i\) wastes work. Assume \(k\le n\), as the requested instance does; then \(k-i\le k\le n\), so every layer is positive and can be divided by. Dividing the formula of Section 7 by itself at \(i-1\), the partition factor contributes \((k-2i+2)(k-2i+1)/(2i)\), the power of \(n\) contributes \(1/n\), and the product loses its last factor \(1-(k-i)/n\). Hence, for \(1\le i\le\lfloor k/2\rfloor\),

$$\frac{P_i}{P_{i-1}}=\frac{(k-2i+2)(k-2i+1)}{2i(n-k+i)}.$$

The numerator forms a new pair from two of the \(k-2i+2\) still-single defects, and the \(2i\) removes the order inside it and among the \(i\) pairs; the code leaves the denominator split, as \(2ni\) times \(1-(k-i)/n\). For \(k\gt n\) the layer formula still holds but the recurrence cannot start, since \(P_0=0\). With \(k=4\), \(n=5\): \(P_0=24/125\), and the ratios \(3\) and \(1/6\) give \(P_1=72/125\), \(P_2=12/125\), whose total \(108/125\) leaves \(p(4,5)=17/125\). Once \(P_0\) is built in \(k\) multiplications each further term costs constant work, so the sum takes time proportional to \(k\), not \(k^2\).

9. Why the computation is done in logarithms

The individual pieces of \(P_i\) lie far outside binary64 arithmetic, which reaches about \(1.8\cdot10^{308}\) and \(4.9\cdot10^{-324}\). [8] At the requested parameters \(k!\) has \(77\,338\) decimal digits, \(n^k=10^{120\,000}\) has \(120\,001\) and \((n)_k\) has \(119\,913\); none can be formed, although their combination is an ordinary probability.

Logarithms turn every product into a sum and every quotient into a difference, all of moderate size. [8] Each factor \(1-j/n\) and each ratio of Section 8 is a plain positive double, so \(\log P_i\) accumulates safely and one \(\exp\) recovers the term. For large \(i\) the term is negligible anyway — the last, \(P_{10\,000}\), is near \(5\cdot10^{-21\,355}\) — and \(\exp\) of so negative a number returns \(0\).

10. Where the mass of the sum sits

The sum has \(\lfloor k/2\rfloor+1=10\,001\) terms at the requested parameters, very unevenly sized: those with \(i\le98\) contribute about \(6.6\cdot10^{-16}\) in total and those with \(i\ge313\) about \(6.3\cdot10^{-16}\), so only \(214\) matter at ten decimal places. The ratio of Section 8 locates the largest. Its numerator falls and its denominator \(2i(n-k+i)\) rises with \(i\), so the ratio decreases steadily, and it is \(1.00077\) at \(i=196\) but \(0.99549\) at \(i=197\); the terms therefore climb to \(P_{196}\approx0.0078093305\) and fall away. The unconditional mean number of chips of load exactly \(2\), \(n\binom k2 n^{-2}(1-1/n)^{k-2}\approx196.03\), sits in the same place.

11. A plausibility check

One chip's load is binomial with \(k\) trials and success probability \(1/n\), hence close to Poisson with mean \(k/n\). [7] With \(q=1-1/n\) the expected number of overloaded chips is exactly

$$\lambda=n\left(1-q^k-\frac kn q^{k-1}-\binom k2\frac{q^{k-2}}{n^2}\right),$$

so a given chip is overloaded with probability \(\lambda/n\). Were those \(n\) events independent the answer would be \(1-(1-\lambda/n)^n\), and replacing that by \(1-e^{-\lambda}\) approximates once more. At \(k=1000\), \(n=10\,000\): \(\lambda\approx1.5422348\), \(1-(1-\lambda/n)^n\approx0.78612291\), \(1-e^{-\lambda}\approx0.78609747\), against the exact \(p(1000,10\,000)\approx0.78831067\) — two decimals only, since chip loads are dependent, sharing a fixed total. The heuristic settles the order of magnitude, nothing more.

How the Code Works

The C++ program [9] keeps its parameters in an Options record defaulting to \(k=20\,000\), \(n=1\,000\,000\) with checkpoints enabled. parse_int_after_prefix accumulates a non-empty run of decimal digits after a given prefix into an int, and parse_arguments uses it for --defects= and --chips=, accepts --skip-checkpoints, rejects anything else and requires \(k\ge0\), \(n\ge1\). That accumulation is not range-checked, so an argument beyond the int range overflows, which is undefined behaviour; the defaults are safe.

solve implements Section 7 directly, returning \(0\) when \(k\le2\), since two defects cannot make a load of three. Otherwise answer starts at \(1\) and the outer loop runs over \(i=0,\dots,\lfloor k/2\rfloor\), rebuilding \(\log P_i\) in logp each time. The first inner loop adds the \(i\) logarithms of \((k-2j)(k-1-2j)/(2n(j+1))\) for \(j\lt i\), whose numerators give \(k!/(k-2i)!\) and denominators \((2n)^i\,i!\) — the partition factor over \(n^i\); the second adds the \(k-i\) logarithms of \(1-j/n\). One std::exp turns logp into \(P_i\), which is subtracted from answer.

run_checkpoints requires solve(3, 7) to land within \(10^{-10}\) of the quoted \(0.0204081633\); the true \(1/49=0.0204081632653\ldots\) is \(3.5\cdot10^{-11}\) away, so it passes, and main prints with std::fixed and std::setprecision(10), the requested form. The Python and Java solve functions also return \(0\) for \(k\le2\), then use the recurrence of Section 8: both initialise logp with the \(k\) logarithms making up \(\log P_0\), subtract its exponential from answer, and loop from \(i=1\) to \(\lfloor k/2\rfloor\), adding the logarithm of \((k-2i+2)(k-2i+1)\) over \(2ni\) and subtracting that of \(1-(k-i)/n\). Python calls those quantities num_term1, den_term1, term2_val, Java numTerm1, denTerm1, term2Val. Neither parses arguments or runs a checkpoint: Python's solve_problem formats the fixed call with "{:.10f}", Java's main prints through System.out.printf under Locale.US.

Otherwise the three differ little: parameters and loop counters are integers, every floating-point quantity is binary64, no big-integer type is needed, nothing recurses, and each program holds a few scalars. All three also share an unstated assumption, that the factors \(1-j/n\) stay positive, which needs \(k\le n\); outside that range none reports the pigeonhole value \(1\) of Section 3, and they fail differently. The C++ sum is still correct at \(k=n+1\), since \(\log0=-\infty\) makes the impossible layers vanish, and returns not-a-number only from \(k\ge n+2\), where the logarithm gets a negative argument; Java's recurrence combines infinities into not-a-number already at \(k=n+1\), and Python stops there with a ValueError from math.log(0.0). So \(p(4,3)=1/3\) and \(p(5,4)=53/128\) come out in C++ alone, \(p(6,4)\) in none of the three.

Complexity and Verification

The C++ outer loop takes \(i+(k-i)=k\) logarithms for every \(i\), hence \((\lfloor k/2\rfloor+1)k=200\,020\,000\) in all: time \(O(k^2)\), memory \(O(1)\). The recurrence needs \(k\) logarithms for \(P_0\) and two per step after, \(k+2\lfloor k/2\rfloor=40\,000\) altogether: time \(O(k)\), memory \(O(1)\). Measured on one machine: about \(0.66\) s for the C++ binary, \(0.04\) s for Python and \(0.06\) s for Java, so the asymptotic gap shows although all are fast.

The published code checks itself only through the C++ \(p(3,7)\) checkpoint. These independent checks were run for this article.

  • Enumerating all \(n^k\) arrangements for \(0\le k\le7\), \(1\le n\le6\) as exact fractions agrees with Section 3 in all \(48\) cases, including \(p(5,4)=53/128\) and the pigeonhole cases \(p(5,2)=p(7,3)=1\).
  • Exact rational evaluation gives \(p(10,20)=0.220963971250000\) exactly, \(p(100,1000)\approx0.140452258322459\) and \(p(1000,10\,000)\)\({}\approx0.788310668551485\), each within \(6\cdot10^{-15}\) of the published Python, and the identity of Section 6 matches the layer sum for its four pairs.
  • At the requested parameters a \(60\)-digit evaluation of the sum differs from the binary64 result by \(7.8\cdot10^{-14}\), about \(700\) units in the last place, while the exact value lies \(2.2\cdot10^{-11}\) from the nearest point where the tenth decimal would round differently — a margin \(280\) times the error, so the printed decimals are safe. The C++ direct sum and the recurrence differ by \(3.6\cdot10^{-14}\) there and round identically; all three programs print the same line.

That discrepancy collects the rounding of the factors, the \(40\,000\) logarithms, the \(10\,001\) exponentials and the subtractions; it stays small because the running variable, starting at \(1\), ends nowhere near \(0\), so nothing cancels catastrophically.

Footnotes and References

The references supply the background facts used above. The layer formula, its proof, the ratio between consecutive layers and every number quoted are derived or computed here.

  1. Project Euler 307 — Chip Defects. The official statement scatters k defects at random over n integrated-circuit chips, any number allowed per chip and each defect independent, defines p(k, n) as the probability of a chip with at least 3 defects, gives p(3, 7) ≈ 0.0204081633 and asks for p(20000, 1000000) to 10 decimal places.
  2. Wikipedia — Birthday problem. The probability that n people have pairwise different birthdays among d equally likely days is the falling factorial over dⁿ, a product of factors 1 − j/d: the layer P₀ of Section 7. Its three-way variant is answered only numerically.
  3. Wikipedia — Balls into bins problem. Each of m balls goes into one of n bins independently and uniformly, a bin's load is its number of balls, and the loads are multinomial; the article collects bounds on the largest load. That is this problem's model.
  4. Wikipedia — Falling and rising factorials. The falling factorial (x)ₙ is the product of x − j for j = 0, …, n − 1 and, for positive integer x, counts the injective functions from an n-element set into an x-element set — the count Section 3 needs.
  5. OEIS A000085 — involutions. The sequence 1, 1, 2, 4, 10, 26, 76, 232, 764, … counts the self-inverse permutations of n letters, satisfies a(n) = Σ n! / ((n − 2k)! · 2^k · k!) and has exponential generating function exp(x + x²/2). Its summands are the partition counts of Section 3, a partition into pairs and singletons being an involution.
  6. P. Flajolet and R. Sedgewick — Analytic Combinatorics. The booksite of the Cambridge University Press edition of 2009 offers the full text free. Its Chapter II, on labelled structures and exponential generating functions, gives the rule used in Section 6: the exponential generating function of an ordered tuple of components built on disjoint label sets is the product of theirs.
  7. Wikipedia — Poisson limit theorem. If n·pₙ → λ, then for each fixed j the binomial probability C(n, j)·pₙ^j·(1 − pₙ)^(n − j) tends to e^(−λ)·λ^j / j!, the Poisson approximation used in Section 11.
  8. Wikipedia — Log probability explains that keeping logarithms turns products into sums and helps when the probabilities are very small, as in Sections 8 and 9. Wikipedia — Double-precision floating-point format gives the binary64 limits of Section 9: largest finite value about 1.7976931348623157 × 10³⁰⁸, smallest positive subnormal about 4.9406564584124654 × 10⁻³²⁴, and 15 to 17 significant decimal digits.
  9. C++ — Euler307.cpp. The linked immutable revision defines parse_int_after_prefix, parse_arguments, solve and run_checkpoints with the defaults 20000 defects and 1000000 chips, and is the implementation record for this article.

Problem 307 source code

C++

#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>

namespace {

struct Options {
    int defects = 20000;
    int chips = 1000000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--defects=", options.defects) ||
            parse_int_after_prefix(arg, "--chips=", options.chips)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.defects >= 0 && options.chips >= 1;
}

double solve(const int defects, const int chips) {
    if (defects <= 2) {
        return 0.0;
    }

    double answer = 1.0;
    for (int i = 0; i <= defects / 2; ++i) {
        double logp = 0.0;

        for (int j = 0; j < i; ++j) {
            logp += std::log((static_cast<double>(defects - 2 * j) *
                              static_cast<double>(defects - 1 - 2 * j)) /
                             (2.0 * static_cast<double>(chips) * static_cast<double>(j + 1)));
        }
        for (int j = 0; j < defects - i; ++j) {
            logp += std::log(1.0 - static_cast<double>(j) / static_cast<double>(chips));
        }

        answer -= std::exp(logp);
    }

    return answer;
}

bool run_checkpoints() {
    const double sample = solve(3, 7);
    if (std::fabs(sample - 0.0204081633) > 1e-10) {
        std::cerr << "Checkpoint failed for p(3,7)" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << std::fixed << std::setprecision(10) << solve(options.defects, options.chips) << '\n';
    return 0;
}

Python

import math

def solve(defects, chips):
    if defects <= 2:
        return 0.0

    answer = 1.0
    
    logp = 0.0
    for j in range(defects):
        logp += math.log(1.0 - float(j) / float(chips))
        
    answer -= math.exp(logp)

    for i in range(1, defects // 2 + 1):
        num_term1 = float(defects - 2 * (i - 1)) * float(defects - 1 - 2 * (i - 1))
        den_term1 = 2.0 * float(chips) * float(i)
        
        term2_j = defects - i
        term2_val = 1.0 - float(term2_j) / float(chips)
        
        logp += math.log(num_term1 / den_term1) - math.log(term2_val)
        
        answer -= math.exp(logp)

    return answer

def solve_problem():
    ans = solve(20000, 1000000)
    return "{:.10f}".format(ans)

if __name__ == '__main__':
    print(solve_problem())

Java

import java.util.Locale;

public class Euler307 {

    static double solve(int defects, int chips) {
        if (defects <= 2) {
            return 0.0;
        }

        double answer = 1.0;

        double logp = 0.0;
        for (int j = 0; j < defects; j++) {
            logp += Math.log(1.0 - (double) j / (double) chips);
        }

        answer -= Math.exp(logp);

        for (int i = 1; i <= defects / 2; i++) {
            double numTerm1 = (double) (defects - 2 * (i - 1)) * (double) (defects - 1 - 2 * (i - 1));
            double denTerm1 = 2.0 * (double) chips * (double) i;

            int term2j = defects - i;
            double term2Val = 1.0 - (double) term2j / (double) chips;

            logp += Math.log(numTerm1 / denTerm1) - Math.log(term2Val);

            answer -= Math.exp(logp);
        }

        return answer;
    }

    public static void main(String[] args) {
        System.out.printf(Locale.US, "%.10f\n", solve(20000, 1000000));
    }
}