Problem 352: Blood Tests

Detailed mathematical approach

Problem Summary

A flock of \(s\) sheep must be screened for a virus that every animal carries independently with probability \(p\). A test on a pooled sample of any set of animals is negative exactly when every animal in it is virus-free, and every sheep must end up with a verdict. One restriction shapes the problem: whenever a mixed sample is tested, all the animals contributing to it must be fully screened before any other animal is examined. Then \(T(s,p)\) is the average number of tests of the cheapest strategy. For \(25\) sheep at \(p=0.02\) the statement computes \(7.40198008\) tests for five groups of five against only \(4.155452\) for the optimal strategy, and records \(T(25,0.10)=12.702124\). We must find the sum of \(T(10000,p)\) over \(p=0.01,0.02,\dots,0.50\), to six decimal places. ^[1]

The whole solution in one sentence: the restriction forces every strategy to be a nest of blocks, so the cheapest cost depends only on a block's size and on whether we already know that it contains an infected animal, and those two quantities satisfy a pair of recurrences that a table fills in \(O(s^2)\) steps per probability.

Mathematical Approach

1. The warm-up scheme

Write \(q=1-p\) for the probability that one animal is virus-free. Since the animals are independent, a pool of \(k\) of them is negative with probability \(q^k\). In the statement's flock \(q=0.98\), so a group of five is negative with probability \(0.98^5=0.9039207968\); otherwise five individual tests follow, and the group costs

$$1+5\,(1-0.98^5)=1.480396016$$

tests on average, so five groups cost \(7.40198008\), a saving of \(70.39\%\) over \(25\) separate tests. This two-stage plan is Dorfman's procedure of 1943. ^[2] The statement hints at three improvements to combine: pool all \(25\) samples first, since \(0.98^{25}=0.6034647297\ldots\) finishes \(60.35\%\) of flocks with one test; leave the last animal of a positive group untested when the others are clear; and choose the group sizes freely.

2. What the restriction forbids

A pool need only hold animals whose status is still open: a known virus-free animal leaves the outcome unchanged, a known infected one makes it certainly positive. While a tested pool still holds unresolved animals, every further test must stay inside it; a negative pool is settled at once, and attention moves outside it. So the work always happens inside one block of \(n\) interchangeable animals that has to be finished first, and a test keeps that block in one of two states: a negative sub-pool hands the block's guarantee on to the rest of it, while a positive sub-pool satisfies that guarantee itself and leaves the rest with its original independent distribution. About a block we therefore know either nothing, or one fact, that at least one of its animals is infected.

3. Two sequences

Let \(A_n\) be the cheapest average number of tests that finishes a block of \(n\) animals about which nothing is known, and \(B_n\) the same for a block known to contain at least one infected animal. The flock is a block of the first kind, so \(T(s,p)=A_s\). The base cases are forced: \(A_0=0\); \(A_1=1\), since a lone animal with no information must be tested; and \(B_1=0\), since a lone animal known to be infected needs no test — the statement's remark about the fifth animal of a positive group. The code also stores \(B_0=0\), never read, since conditioning an empty block on an infection means nothing.

4. The recurrence without prior information

In a block of \(n\) animals about which nothing is known, let the first pool hold \(k\le n\) of them. With probability \(q^k\) it is negative, those \(k\) are virus-free, and the block shrinks to the \(n-k\) others with no information. With probability \(1-q^k\) it is positive; the restriction forces us to finish those \(k\) first, knowing only that one of them is infected, which costs \(B_k\). The outcome concerns only the pool, so the animals outside stay unconditioned. By the law of total expectation ^[5] the cost is \(1+(1-q^k)B_k+A_{n-k}\), where \(A_{n-k}\) has weight \(1\) because it is paid in both branches. Hence

$$A_n=\min_{1\le k\le n}\bigl(1+(1-q^k)B_k+A_{n-k}\bigr).$$

The case \(k=n\) is the choice "pool the whole block". Naming the cost of finishing a block of \(k\) animals in one piece,

$$\begin{aligned}C_k&=1+(1-q^k)B_k,\\ A_n&=\min_{1\le k\le n}\bigl(C_k+A_{n-k}\bigr),\end{aligned}$$

makes \(C_1=1\) plain individual testing, and shows what \(A_n\) is: the cheapest total over all ways of cutting \(n\) animals into blocks, with cost \(C_k\) for a block of size \(k\).

5. The recurrence with one infected animal guaranteed

Now the block is known to contain an infected animal. Pooling all \(n\) would be positive for certain and waste a test, so the pool holds \(k\lt n\) animals. By the definition of conditional probability ^[4] it is negative with probability

$$\alpha_k=\frac{q^k\,(1-q^{n-k})}{1-q^n},$$

because a negative pool together with the known fact forces an infected animal among the other \(n-k\); a positive pool already implies one, so its conditional probability is \(1-\alpha_k=(1-q^k)/(1-q^n)\). In the negative branch the tested animals are clear and the remaining \(n-k\) inherit the guarantee, costing \(B_{n-k}\); in the positive branch the pool is finished first, costing \(B_k\), and the other \(n-k\) are unconditioned, costing \(A_{n-k}\). Therefore

$$\begin{aligned}B_n=\min_{1\le k\lt n}\Bigl(&1+\alpha_k B_{n-k}\\&+(1-\alpha_k)\bigl(B_k+A_{n-k}\bigr)\Bigr).\end{aligned}$$

Each expression is a lower bound for every legal strategy whose first pool has that size, and a minimising pool followed by optimal continuations attains it, so induction on \(n\) proves that \(A_n\) and \(B_n\) are the true optima. Both recurrences use indices below \(n\), except that \(A_n\) uses \(B_n\), so the table is filled for \(n=2,\dots,s\) with \(B_n\) before \(A_n\): bottom-up dynamic programming, each subproblem solved once. ^[6] These are the recurrences of Sobel and Groll's 1959 procedure \(R_3\), the variant that forbids recombining sets, with \(H_3=A\), \(G_3=B\). ^[3]

6. Two and three sheep by hand

For two animals known to contain an infected one, the only choice is to test one of them. It is clear with probability \(q(1-q)/(1-q^2)=q/(1+q)\), and then its partner is infected and needs no test; otherwise the tested animal is infected and the partner is still open, costing \(A_1=1\). Hence

$$\begin{aligned}B_2&=1+\frac1{1+q},\\ A_2&=\min\bigl(3-q-q^2,\,2\bigr),\end{aligned}$$

the first entry of the minimum being \(C_2\) and the second two individual tests. At \(p=0.1\): \(B_2=29/19=1.526315789\ldots\) and \(A_2=1.29\) exactly.

For three animals with the guarantee, testing one first gives \(\alpha_1=q(1+q)/(1+q+q^2)\) and the cost \(1+\alpha_1B_2+(1-\alpha_1)A_2\), while testing two first gives \(1+(3+2q)/(1+q+q^2)\); at \(q=0.9\) these are \(661/271\) and \(751/271\), so \(B_3=661/271=2.439114391\ldots\), and \(A_3\) then compares \(C_1+A_2=2.29\), \(C_2+A_1=2.29\) and \(C_3\), where the whole-block option wins:

$$A_3=1+0.271\cdot\tfrac{661}{271}=1.661.$$

7. Five and twenty-five sheep at \(p=0.02\)

Now \(q=0.98\). For five animals the table says: pool all five, finishing the \(90.39\%\) of flocks whose test is negative. Otherwise we are in the guaranteed state with \(n=5\), whose best first pool holds two animals, \(\alpha_2=0.5878400457\ldots\):

$$B_5=1+\alpha_2B_3+(1-\alpha_2)\bigl(B_2+A_3\bigr).$$

With \(B_3=2.353557339\ldots\), \(B_2=1.505050505\ldots\) and \(A_3=1.138408\) this is \(B_5=3.473042990\ldots\), so

$$\begin{aligned}A_5&=1+0.0960792032\cdot B_5\\&=1.333687203\ldots,\end{aligned}$$

against \(1.480396016\) for "pool the five, then test all five". For \(25\) animals the table again pools the whole flock first; the best second pool holds \(8\) of the \(25\), and \(B_{25}=7.957557479\ldots\) gives

$$\begin{aligned}A_{25}&=1+(1-0.98^{25})B_{25}\\&=4.155452205\ldots,\end{aligned}$$

which rounds to the checkpoint \(4.155452\) and saves \(83.38\%\) of \(25\) separate tests.

8. Blocks, and the cost per animal

Since \(A_n\) is the cheapest sum of block costs, the ratio \(C_k/k\) is the number of tests per animal that a block of size \(k\) charges, and \(A_n\ge n\cdot\min_k(C_k/k)\). At \(p=0.10\) the cheapest ratio belongs to \(k=8\), where \(C_8=4.05204279\), and the second checkpoint falls out as three blocks, \(25=8+8+9\):

$$2\,C_8+C_9=12.702124091$$

with \(C_9=4.598038511\), while \(8+8+8+1\) would cost \(13.15612837\), because a leftover single animal pays a whole test: the best ratio alone does not decide the answer, the best composition does. As \(p\) falls the ideal block grows: for \(s=10000\) its size is \(2\) at \(p=0.30\), \(8\) at \(p=0.10\) and \(89\) at \(p=0.01\), where \(T(10000,0.01)\approx964.186678\), about \(0.0964\) tests per animal.

9. Where pooling stops paying

Pooling a pair beats two individual tests exactly when \(C_2\lt2\), that is \(3-q-q^2\lt2\), or \(q^2+q\gt1\). The positive root of \(q^2+q=1\) is \(1/\varphi=\varphi-1\) for the golden ratio \(\varphi=(1+\sqrt5)/2\), and the identity \(1=1/\varphi+1/\varphi^2\) turns the condition into ^[7]

$$\begin{aligned}p&\lt\varphi^{-2}=\frac{3-\sqrt5}{2}\\&=0.3819660112501\ldots\end{aligned}$$

Ungar proved in 1960 that this is the exact cutoff of the whole problem, unrestricted procedures included: at or above it individual testing is optimal — at the cutoff \(C_2=2\), so pooling a pair only ties — and below it pooling wins for every flock of at least two animals, while no general optimal procedure below the cutoff is known. ^[2] That is why this problem imposes its restriction. Since \(A_n\) is a sum of block costs, \(A_s=s\) holds exactly when \(C_k\ge k\) for every \(k\le s\), and the table confirms \(A_{10000}=10000\) for every \(p\ge0.39\): twelve of the fifty terms are plain individual testing. Just below the cutoff pooling pairs still pays, and \(p=0.38\) gives \(5000\,C_2=9978\) exactly.

How the Code Works

The C++ function T(const int s, const long double p) ^[9] is Sections 3 to 5 line by line. It returns \(0\) for \(s\le0\) and \(1\) for \(s=1\), builds qpow with qpow[i]\({}=q^i\) by repeated multiplication, and fills A, B of length \(s+1\) from A[1]\({}=1\), B[1]\({}=0\). For each \(n\) it computes pe\({}=1-q^n\) once, scans \(k=1,\dots,n-1\) with pa\({}=\alpha_k\), pb\({}=1-\alpha_k\) and keeps the smallest conditional cost in bestB, whose sentinel 1e100L is always replaced since \(n\ge2\); then it sets bestA to the \(k=n\) term \(1+(1-q^n)\)B[n] of Section 4 and compares it with the split terms \(1+\)A[n-k]\({}+(1-q^k)\)B[k]. solve adds the fifty values T(10000, i/100), each computed independently.

Only the C++ validates itself: run_checkpoints compares T(25,0.02) and T(25,0.10) with \(4.155452\) and \(12.702124\) through approx_equal_6dp, which demands a difference below \(5\cdot10^{-7}\); main takes only --skip-checkpoints, returns the error code \(2\) on a failed checkpoint, and otherwise prints the sum with setprecision(6). The Java port repeats the recurrence with double[] arrays but has no checkpoint; if collecting its results throws, it prints a stack trace and then formats the subtotal reached so far, which is not a finished computation. The Python file is no second implementation: it compiles the C++ at -O2 into a cached binary when that is missing or stale, runs it, and reads the last answer: line, else the last = line, else the last non-empty line — adding --skip-checkpoints whenever that string occurs in the C++ source, so the checkpoints are skipped on this path.

The two real implementations differ in their arithmetic type. The format of long double is left to the implementation ^[8]; on the arm64 Mac used here its significand has Java's \(53\) bits, on many x86-64 compilers \(64\). Equal precision is not identical rounding: a C++ compiler may fuse a multiplication and an addition into one rounded operation, and the C++ and Java values of \(T(25,0.02)\) differ in the sixteenth significant digit — compiling the C++ without that fusion reproduces the Java value exactly. In qpow the powers reach the subnormal range, where repeated multiplication may stall at a tiny positive value instead of reaching zero: at \(p=0.10\) qpow[10000] ends at about \(2.5\cdot10^{-323}\), at \(p=0.50\) at zero, both far below what six decimals see. The denominator stays safe, since \(p\ge0.01\) and \(n\ge2\) give pe\({}\ge1-0.99^2=0.0199\); only a \(p\) with \(1-p\) rounding to \(1\) would divide by zero.

Complexity and Verification

For one probability the two scans examine

$$\sum_{n=2}^{s}2\,(n-1)=s\,(s-1)=99\,990\,000$$

candidate pool sizes at \(s=10000\), so the time is \(\Theta(s^2)\) and the memory \(\Theta(s)\): three arrays of \(s+1\) numbers, about \(240\) kilobytes at eight bytes each. The fifty probabilities need \(4\,999\,500\,000\) candidates in all; measured here one takes about \(60\) milliseconds compiled with -O2, so the fifty are a few seconds of processor work, and Java is comparable.

The checks behind this article are these.

• Both checkpoints come out, \(4.155452205\ldots\) and \(12.702124091\), and the C++ stops with the error code \(2\) if either is wrong. C++ and Java print the same six decimals; the Python wrapper reproduces them by running the C++ build, so it is not independent evidence.
• Every nested strategy for a block of at most five animals was enumerated in exact rational arithmetic at four probabilities, and its minimum equals the recurrence in all cases: \(A_5=416777251/312500000\) at \(p=0.02\), \(A_5=5\) at \(p=0.4\). A second computation over arbitrary subsets of up to nine animals, choosing freely which animals to pool, returns the same values: only a block's size matters, as Section 2 claims.
• A separate program drove the pool sizes the table selects on random flocks and counted the tests performed: three million flocks of \(25\) animals averaged \(4.155907\) at \(p=0.02\) and \(12.700457\) at \(p=0.10\), within two standard errors of the table at each of three seeds, one standard error of the mean being \(0.0024\) and \(0.0030\). A strategy attaining the recurrence therefore exists.
• Against exact rational arithmetic the double table is off by at most \(7.2\cdot10^{-15}\) at \(s=30\), and against \(40\)-digit arithmetic by \(4.3\cdot10^{-12}\) at \(s=1200\) and \(2.2\cdot10^{-11}\) at \(s=2400\). The error grows with \(s\), so these checks support the implementation at the sizes tested; certifying the sixth decimal would need a bound at \(s=10000\) and over all fifty terms.
• The block form of Section 4, evaluated separately for \(n\le600\) at four probabilities, agrees with the published order of the additions to within \(1.5\cdot10^{-14}\), and \(A_{10000}=10000\) holds exactly at the twelve probabilities above Ungar's cutoff, as it must.

Footnotes and References

The references below supply the background facts used above. The block decomposition, the two recurrences, the worked examples and the operation counts are derived explicitly in this article.

1. Project Euler 352 — Blood Tests. The official statement puts the five-groups-of-five scheme at 7.40198008 tests, requires every animal in a mixed sample to be fully screened before others are examined, gives T(25, 0.02) = 4.155452 and T(25, 0.10) = 12.702124, and asks for the sum of T(10000, p) over the fifty probabilities.
2. Wikipedia — Group testing. The article credits Robert Dorfman with introducing group testing in 1943 to screen army recruits, grouping the men and testing the members of positive groups individually. It also gives Peter Ungar's cutoff of 1960, p_u = (3 − √5)/2 ≈ 0.38, above which individual testing minimises the expected number of tests, and states that below it a general optimal procedure is still unknown.
3. M. Sobel and P. A. Groll, Group testing to eliminate efficiently all defectives in a binomial sample, Bell System Technical Journal 38 (1959), pages 1179–1252, in the Internet Archive. Each unit is an independent observation from a binomial population with a known probability of being good. Appendix B defines procedure R3, in which two sets once formed cannot be recombined, so sets are only broken into smaller ones, "yielding a nested set of partitions"; its equations (141) and (142) are the recurrences used here, H3(n) for a set with no known defective and G3(m) for one known to contain a defective, with H3(0) = G3(1) = 0.
4. Wikipedia — Conditional probability. The article defines P(A | B) = P(A ∩ B) / P(B) for P(B) greater than zero and calls A and B independent when the probability of their intersection is the product of their probabilities.
5. Wikipedia — Law of total expectation. For a partition into events Aᵢ the article gives E(X) = Σᵢ E(X | Aᵢ) P(Aᵢ), the form Sections 4 and 5 apply to the two branches of a pooled test.
6. Wikipedia — Dynamic programming. The article explains optimal substructure, where a solution is built from optimal solutions of subproblems, and overlapping subproblems, which dynamic programming solves once each, bottom-up in tabular form.
7. Wikipedia — Golden ratio. The article gives φ = (1 + √5)/2, states that φ² = φ + 1 and 1/φ = φ − 1 = 0.618033…, and records the identity 1 = 1/φ + 1/φ², which Section 9 uses to identify the cutoff with φ⁻².
8. cppreference — Fundamental types. The page calls long double an extended-precision type that need not map to an IEEE-754 format, and notes that many x86 and x86-64 implementations use the 80-bit x87 format, one notable implementation the binary64 format of double, and some ARM64 implementations binary128.
9. C++ — Euler352.cpp. The linked immutable C++ revision contains T, solve, approx_equal_6dp, run_checkpoints and main, and compares T(25, 0.02) and T(25, 0.10) with the two checkpoints. It is the implementation record for this article.