Problem 236: Luxury Hampers

Project Euler Problem 236 Solution

EulerSolve provides an optimized solution for Project Euler Problem 236, Luxury Hampers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The problem gives five product totals for two suppliers' luxury hampers. To keep the algebra neutral, write those totals as $$P=(5248,1312,2624,5760,3936),\qquad Q=(640,1888,3776,3776,5664).$$ Their grand totals are $$T_P=18880,\qquad T_Q=15744.$$ We seek the largest rational number $m>1$ for which one can choose positive integers $p_i$ and $q_i$ for all five products so that the same multiplicative bias $m$ holds both product-by-product and after the five products are combined. The direct search over all five pairs $(p_i,q_i)$ is far too large, so the solution looks for structure in the fixed supplier totals. Mathematical Approach Write the desired multiplier as $m=\frac{u}{v}$ in lowest terms. The implementations turn the problem into an exact arithmetic feasibility test: first generate all possible rational candidates for $m$, then decide which of them admit positive integer counts satisfying the product constraints and the global consistency equation. Ratio Equations for the Five Products For each product $i$, the implementation works with positive integers $p_i$ and $q_i$ constrained by $$vP_iq_i=uQ_ip_i.$$ The combined condition over all five products is $$vT_Q\sum_{i=1}^{5}p_i=uT_P\sum_{i=1}^{5}q_i.$$ So a candidate $m$ is valid if and only if these six integer equations can be satisfied simultaneously with every $p_i,q_i\ge 1$....

Detailed mathematical approach

Problem Summary

Two suppliers, A and B, deliver five products to the luxury hamper market: A delivers $5248$ items of Beluga Caviar, $1312$ of Christmas Cake, $2624$ of Gammon Joint, $5760$ of Vintage Port and $3936$ of Champagne Truffles, $18880$ in all, and B delivers $640$, $1888$, $3776$, $3776$ and $5664$ of them, $15744$ in all. Some items spoil. A spoilage rate is the number of spoiled items divided by the number delivered, for one product or, for the overall rate, for all five products of a supplier together. For each of the five products the rate of B is $m$ times the rate of A, with the same factor $m\gt1$ for all five, and yet the overall rate of A is $m$ times that of B. The statement says that $35$ values of $m\gt1$ make this possible, the smallest being $1476/1475$, and asks for the largest possible $m$ as a fraction $u/v$ in lowest terms. ^[1]

The numbers of spoiled items are whole numbers. In each product the rate of B is higher than that of A, so it is positive, and since it is $m$ times the rate of A, that rate is positive too. So each supplier has at least one spoiled item of every product, as the published programs assume. The situation is an extreme case of Simpson's paradox, in which a comparison that holds inside every group is reversed when the groups are put together. ^[2]

The whole solution in one sentence: the caviar leaves only about a million possible fractions for $m$; for each of them the per-product conditions fix all spoiled counts up to one whole-number multiplier per product, and the overall condition becomes one linear equation in three bounded whole numbers, which the extended Euclidean algorithm settles in a few steps.

Mathematical Approach

We first see on a tiny table how the paradox can happen and why it limits $m$, and then reduce the six conditions to one linear equation in three bounded whole numbers.

1. A paradox with fifteen items

Take two products. Supplier A delivers $5$ items of each; supplier B delivers $1$ item of the first product and $4$ of the second. Suppose A has $4$ spoiled items of the first product and $1$ of the second, and B has $1$ spoiled item of each. The rates of A are $4/5$ and $1/5$, those of B are $1$ and $1/4$, and in both products the rate of B is $5/4$ times that of A: $\tfrac54\cdot\tfrac45=1$ and $\tfrac54\cdot\tfrac15=\tfrac14$. Overall, A has $5$ spoiled items out of $10$, the rate $1/2$, and B has $2$ out of $5$, the rate $2/5$. Now the rate of A is $5/4$ times that of B, because $\tfrac54\cdot\tfrac25=\tfrac12$. The reason is the mix: half of the items of A belong to the badly spoiling first product, but only a fifth of the items of B. An overall rate is a weighted average of the per-product rates, and different weights can reverse a comparison.

2. The unknowns and the six conditions

Number the products $i=1,\dots,5$ in the order of the statement, so that product $1$ is the caviar and product $4$ the port. Let $A_i$ and $B_i$ be the numbers delivered by A and B, let $S_A=18880$ and $S_B=15744$ be their totals, and let $x_i$ and $y_i$ be the numbers of spoiled items, whole numbers with $1\le x_i\le A_i$ and $1\le y_i\le B_i$. We call $m\gt1$ valid if such numbers satisfy

$$\frac{y_i}{B_i}=m\,\frac{x_i}{A_i}\qquad(i=1,\dots,5),$$

$$\frac{x_1+\dots+x_5}{S_A}=m\,\frac{y_1+\dots+y_5}{S_B}.$$

3. How large m can be

Write $r_i=x_i/A_i$ for the rate of A in product $i$. The overall rate of A is $\sum_i(A_i/S_A)\,r_i$, and that of B is $\sum_i(B_i/S_B)\,m\,r_i$. So the overall condition says that $\sum_i(A_i/S_A)\,r_i$ equals $m^2\sum_i(B_i/S_B)\,r_i$. With the positive weights $w_i=(B_i/S_B)\,r_i$ and the numbers $\rho_i=(A_i/S_A)/(B_i/S_B)$, which compare the shares of product $i$ in the two deliveries, the left-hand side is $\sum_i\rho_i\,w_i$, and so

$$m^2=\frac{\sum_i\rho_i\,w_i}{\sum_i w_i}.$$

Thus $m^2$ is a weighted average of the $\rho_i$, a weighted mediant, ^[3] and it lies strictly between the smallest and the largest $\rho_i$: multiply $\rho_{\min}\le\rho_i\le\rho_{\max}$ by $w_i$ and add; equality would need all $\rho_i$ to be equal. For the hamper table, $\rho_i=10086/1475\approx6.838$ for the caviar, $10086/17405\approx0.579$ for cake, gammon and truffles, and $4428/3481\approx1.272$ for the port. Hence every valid $m$ satisfies

$$m\lt\sqrt{10086/1475}\approx2.615.$$

A large $m$ needs most of the weight on the caviar, $27.8\%$ of the delivery of A but only $4.07\%$ of that of B. In Section 1, $\rho_1=\tfrac12\big/\tfrac15=\tfrac52$ and $\rho_2=\tfrac12\big/\tfrac45=\tfrac58$, both weights are $4/25$, and $m^2=25/16$ is the average of $5/2$ and $5/8$.

4. The table in lowest terms

The entries of the table share large factors: $5248=128\cdot41$ and $640=128\cdot5$; $1312=32\cdot41$ and $1888=32\cdot59$; $2624=64\cdot41$ and $3776=64\cdot59$; $5760=64\cdot90$ and $3776=64\cdot59$; $3936=96\cdot41$ and $5664=96\cdot59$. The totals are $S_A=64\cdot295$ and $S_B=64\cdot246$. So the ratio $B_i/A_i$ takes only three values: $5/41$ for the caviar, $59/41$ for cake, gammon and truffles, and $59/90$ for the port. We call the three sets of products the groups $G_0=\{1\}$, $G_1=\{2,3,5\}$ and $G_2=\{4\}$, and we write the ratio of group $g$ as $\beta_g/\alpha_g$.

5. One multiplier per product

Write $m=u/v$ in lowest terms. For a product $i$ of group $g$, the per-product condition says $y_i/x_i=m\,B_i/A_i$${}=u\beta_g/(v\alpha_g)$. Let $n_g/d_g$ be this fraction in lowest terms; it depends only on $m$ and on the group. From $d_g\,y_i=n_g\,x_i$ and $\gcd(n_g,d_g)=1$, Euclid's lemma shows that $d_g$ divides $x_i$, ^[4] so

$$x_i=d_g\,t_i,\qquad y_i=n_g\,t_i$$

with a positive whole number $t_i$, the multiplier of product $i$; conversely, every such pair satisfies the per-product condition. The limits $x_i\le A_i$ and $y_i\le B_i$ become $1\le t_i\le U_i$ with

$$U_i=\min\bigl(\lfloor A_i/d_g\rfloor,\lfloor B_i/n_g\rfloor\bigr).$$

If $U_i=0$ for some product, $m$ is impossible. For example, for $m=369/295$ the caviar has $\frac{n_0}{d_0}=\frac{369\cdot5}{295\cdot41}$${}=\frac9{59}$, so $x_1=59\,t_1$ and $y_1=9\,t_1$, with $U_1=\min(88,71)=71$.

6. Finitely many candidates

The caviar gives $m=(y_1/B_1)/(x_1/A_1)$${}=41y_1/(5x_1)$, since $A_1/B_1=41/5$. So every valid $m$ is one of the reduced fractions

$$m=\frac{41b}{5a},\qquad1\le a\le5248,\quad1\le b\le640,$$

where $a=x_1$ and $b=y_1$. Of the $3\,358\,720$ pairs $(a,b)$, $1\,681\,600$ give $m\gt1$, and they produce $1\,021\,210$ different fractions, whose numerators and denominators are at most $41\cdot640=5\cdot5248=26240$. Only $7040$ of them give $U_i\ge1$ for all five products.

7. Three sums that fill whole intervals

For a fixed $m$ with $U_i\ge1$ for every product, the products of one group share $n_g$ and $d_g$, so only the sum of their multipliers matters:

$$s_0=t_1,\qquad s_1=t_2+t_3+t_5,\qquad s_2=t_4.$$

The smallest possible value of $s_g$ is the number of products in the group, $L_g=1$, $3$ and $1$, and the largest is the sum $H_g$ of the bounds $U_i$ of its products. Every whole number between $L_g$ and $H_g$ occurs as well: start with all $t_i=1$ and raise the multipliers by $1$ at a time, each one up to its bound. Every step raises the sum by exactly $1$, and the last step reaches $H_g$. For $m=369/295$ the bounds of cake, gammon and truffles are $209$, $419$ and $629$, so $s_1$ can be any number from $3$ to $1257$.

8. The overall condition as one linear equation

Without denominators, the overall condition reads $v\,S_B\sum_ix_i=u\,S_A\sum_iy_i$. With $x_i=d_g\,t_i$ and $y_i=n_g\,t_i$ the sums of Section 7 appear, and the condition becomes

$$c_0s_0+c_1s_1+c_2s_2=0,\qquad c_g=v\,S_B\,d_g-u\,S_A\,n_g.$$

Hence $m$ is valid exactly when this equation has a solution with $L_g\le s_g\le H_g$ for $g=0,1,2$. Every valid $m$ gives one; conversely, splitting each $s_g$ into multipliers within their bounds as in Section 7 and putting $x_i=d_g\,t_i$ and $y_i=n_g\,t_i$ gives counts that satisfy the per-product conditions by Section 5, the limits since $t_i\le U_i$, and the overall condition by the equation.

The signs of the coefficients repeat Section 3. Since $n_g/d_g=m\,\beta_g/\alpha_g$, we have $c_g\gt0$ exactly when $m^2\lt(S_B/S_A)(\alpha_g/\beta_g)$, and the right-hand side is the value $\rho_i$ of the products of group $g$. For $m\gt1$ the coefficient $c_1$ is always negative, $c_2$ is positive exactly when $m\lt\sqrt{4428/3481}\approx1.128$, and $c_0$ is positive exactly when $m\lt\sqrt{10086/1475}$, the bound of Section 3. Above that bound all three coefficients are negative and there is no solution. No $c_g$ is ever $0$, because no $\rho_i$ is the square of a fraction: in $10086/1475=2\cdot3\cdot41^2/(5^2\cdot59)$, $10086/17405=2\cdot3\cdot41^2/(5\cdot59^2)$ and $4428/3481=2^2\cdot3^3\cdot41/59^2$ some prime has an odd exponent. Of the $7040$ candidates of Section 6, $4430$ lie below the bound.

9. Fixing the caviar sum: Bézout's identity

The program tries $s_0=1,2,\dots,H_0$ in turn, at most $640$ values, because $y_1=n_0\,s_0\le640$. For a fixed $s_0$ it needs integers $s_1$ and $s_2$ with

$$c_1s_1+c_2s_2=r,\qquad r=-c_0s_0.$$

By Bézout's identity, the numbers $c_1x+c_2y$ with integers $x$ and $y$ are exactly the multiples of $g=\gcd(c_1,c_2)$, and the extended Euclidean algorithm finds $x_0$ and $y_0$ with $c_1x_0+c_2y_0=g$. ^[5] ^[6] So there is no solution unless $g$ divides $r$. If it does, $s_1=x_0\,r/g$ and $s_2=y_0\,r/g$ form one solution, and all solutions are

$$\begin{aligned}s_1&=x_0\,\frac rg+\frac{c_2}g\,k,\\s_2&=y_0\,\frac rg-\frac{c_1}g\,k,\end{aligned}$$

with any integer $k$. Indeed, if $(s_1,s_2)$ and $(s_1',s_2')$ are two solutions, then $\frac{c_1}g(s_1-s_1')$${}=-\frac{c_2}g(s_2-s_2')$. The numbers $c_1/g$ and $c_2/g$ have no common factor, so by Euclid's lemma $s_1-s_1'$ is a multiple $k\,c_2/g$, and then $s_2-s_2'=-k\,c_1/g$.

10. Which k stay within the bounds

Both $s_1$ and $s_2$ have the form $p+qk$ with $q\ne0$, and each must lie in its interval $[L,H]$. Dividing $L\le p+qk\le H$ by $q$, and turning the inequalities around when $q$ is negative, gives

$$\Bigl\lceil\frac{L-p}q\Bigr\rceil\le k\le\Bigl\lfloor\frac{H-p}q\Bigr\rfloor\quad(q\gt0),$$

$$\Bigl\lceil\frac{H-p}q\Bigr\rceil\le k\le\Bigl\lfloor\frac{L-p}q\Bigr\rfloor\quad(q\lt0).$$

Here $\lfloor z\rfloor$ is the largest integer not above $z$, and $\lceil z\rceil=-\lfloor-z\rfloor$ is the smallest integer not below $z$; for example, $\lfloor-7/3\rfloor=-3$ and $\lceil-7/3\rceil=-2$. ^[7] The value $s_0$ works exactly when the two intervals of $k$ overlap, and $m$ is valid exactly when this happens for some $s_0$ between $1$ and $H_0$.

11. Worked example: m = 369/295

The fraction $369/295=41\cdot9/(5\cdot59)$ comes from the caviar with $a=59$ and $b=9$. The three groups have $n_g/d_g=9/59$, $9/5$ and $41/50$; for instance $\frac{369\cdot59}{295\cdot41}=\frac95$, because $369=9\cdot41$ and $295=5\cdot59$. The bounds are $\min(88,71)=71$ for the caviar, $\min(262,209)=209$, $\min(524,419)=419$ and $\min(787,629)=629$ for cake, gammon and truffles, and $\min(115,92)=92$ for the port, so

$$1\le s_0\le71,\qquad3\le s_1\le1257,\qquad1\le s_2\le92.$$

The coefficients are $c_g=295\cdot15744\,d_g$${}-369\cdot18880\,n_g$, that is, $c_g=2322240\,(2d_g-3n_g)$, which gives $(c_0,c_1,c_2)=2322240\cdot(91,-17,-23)$. The same equation follows directly: the overall condition demands $\sum x_i/\sum y_i=m\,S_A/S_B=3/2$, and with $\sum x_i=59s_0+5s_1+50s_2$ and $\sum y_i=9s_0+9s_1+41s_2$ this is

$$91\,s_0=17\,s_1+23\,s_2.$$

Try $s_0=1$. The extended Euclidean algorithm divides with remainder, $23=1\cdot17+6$, $17=2\cdot6+5$ and $6=1\cdot5+1$, and then substitutes back: $1=6-5=3\cdot6-17$${}=3\cdot23-4\cdot17$. Multiplying by $91$ gives the solution $s_1=-364$, $s_2=273$, and all solutions are $s_1=-364+23k$, $s_2=273-17k$. The bounds $3\le s_1\le1257$ allow $16\le k\le70$, the bounds $1\le s_2\le92$ allow $11\le k\le16$, and the only common value $k=16$ gives $s_1=4$ and $s_2=1$. So $m=369/295$ is valid. Splitting $s_1=4$ as $2+1+1$ gives a complete table of spoiled items:

Beluga Caviar: A $59$ of $5248$ ($1.124\%$), B $9$ of $640$ ($1.406\%$);
Christmas Cake: A $10$ of $1312$ ($0.762\%$), B $18$ of $1888$ ($0.953\%$);
Gammon Joint: A $5$ of $2624$ ($0.191\%$), B $9$ of $3776$ ($0.238\%$);
Vintage Port: A $50$ of $5760$ ($0.868\%$), B $41$ of $3776$ ($1.086\%$);
Champagne Truffles: A $5$ of $3936$ ($0.127\%$), B $9$ of $5664$ ($0.159\%$).

In every line the rate of B is exactly $369/295$ times that of A, for instance $\frac{9/640}{59/5248}=\frac{47232}{37760}$${}=\frac{369}{295}$. Overall, A has $129$ spoiled items, the rate $129/18880\approx0.683\%$, and B has $86$, the rate $86/15744\approx0.546\%$, so the comparison is reversed, and $\frac{129/18880}{86/15744}=\frac32\cdot\frac{246}{295}$${}=\frac{369}{295}$.

12. A candidate that fails: m = 2

The value $m=2$ lies below the bound of Section 3 and comes from the caviar with $a=41$ and $b=10$, but it is not valid. Here $n_g/d_g=10/41$, $118/41$ and $59/45$, the bounds are $H_0=64$, $H_1=16+32+48=96$ and $H_2=64$, and the overall condition $\sum x_i/\sum y_i=2S_A/S_B=295/123$ becomes

$$2093\,s_0=29767\,s_1+11870\,s_2.$$

The left-hand side is at most $2093\cdot64=133952$, because B delivers only $640$ items of caviar, and the right-hand side is at least $101171$. Only $(s_1,s_2)=(3,1)$, $(3,2)$, $(3,3)$ and $(4,1)$ keep the right-hand side small enough, and none of the values $101171$, $113041$, $124911$ and $130938$ is a multiple of $2093$, so no integer $s_0$ fits.

13. The smallest value of the statement

For $m=1476/1475$ the groups have $n_g/d_g=36/295$, $36/25$ and $82/125$, with $H_0=17$, $H_1=313$ and $H_2=46$. Since $m^2\lt4428/3481$, the coefficient $c_2$ is positive, and the equation, divided by $\gcd(c_0,c_1,c_2)=4644480$, reads $1259\,s_0+133\,s_2=91\,s_1$. There is no admissible solution for $s_0=1,\dots,6$, and $s_0=7$ works, for example with $s_1=110$ and $s_2=9$: $1259\cdot7+133\cdot9=10010$${}=91\cdot110$. The five smallest valid values are $1476/1475$, $60/59$, $902/885$, $3321/3245$ and $41/40$.

How the Code Works

The C++ program ^[8] stores the table in the arrays A and B, the groups in GROUP_PRODUCTS with the zero-based indices $\{0\}$, $\{1,2,4\}$ and $\{3\}$, and the ratios $\beta_g/\alpha_g$ as GROUP_B_NUM $=(5,59,59)$ over GROUP_A_DEN $=(41,41,90)$. A Fraction holds u and v as 64-bit integers; reduce_fraction makes the denominator positive and divides by the gcd, and fraction_less compares by cross-multiplying in the 128-bit type __int128_t, a compiler extension. ^[9]

main first calls generate_candidates, which is Section 6: it reserves room for all $3\,358\,720$ pairs, reduces $41b/(5a)$ for $a\le5248$ and $b\le640$, keeps the fractions above $1$, sorts them with fraction_less and removes repeats with std::unique. candidate_works is Sections 5, 7 and 8. For each group, reduce_fraction applied to $u\beta_g$ and $v\alpha_g$ yields N $=n_g$ and D $=d_g$; every product of the group gets the bound upper $=\min(A_i/D,\,B_i/N)$ by integer division, and a bound $0$ rejects the candidate at once. The bounds add up to max_sum[g] $=H_g$, min_sum $=(1,3,1)$ holds the $L_g$, and coef[g] is $c_g$. Then $s_0$ runs from $1$ to $H_0$ until exists_bounded_linear finds $s_1$ and $s_2$.

exists_bounded_linear is Sections 9 and 10. Its branch for a zero coefficient never runs here, by Section 8. extended_gcd computes $g$ and Bézout coefficients recursively with the C++ operators / and %, whose quotient is truncated towards zero. The back-substitution only needs the identity (a / b) * b + a % b == a, which holds for truncation too, ^[10] so negative coefficients are handled correctly. If $g$ divides the right-hand side, the particular solution $(x_0\,r/g,\,y_0\,r/g)$ is formed in 128-bit arithmetic, with a guard that answers no if it does not fit into 64 bits; the guard never fires here. tighten_range turns one bound $L\le p+qk\le H$ into an interval of $k$ as in Section 10 and intersects the intervals. Because C++ division truncates, floor_div lowers the quotient by $1$ when the remainder is not $0$ and has the opposite sign to the divisor, and ceil_div(a, b) is -floor_div(-a, b). For $m=369/295$ and $s_0=1$, exists_bounded_linear gets $x_0=4$ and $y_0=-3$ from the coefficients of Section 11, and the single admissible $k=-16$ gives $s_1=4$ and $s_2=1$.

solve_all applies candidate_works to every candidate and returns the valid fractions in increasing order. Unless --skip-checkpoints is given, run_checkpoints checks that there are exactly $35$ of them, that the smallest is $1476/1475$, and that candidate_works accepts $1476/1475$ directly. main prints the largest valid fraction as u/v and exits with code $2$ if there is none and $3$ if a checkpoint fails.

The Python and Java versions follow the same steps, Java with its names in camelCase, but have no checkpoints: their solve tests every candidate with candidate_works or candidateWorks and returns the last valid fraction. Python's integers have no size limit and its // already rounds down, so floor_div is simply a // b. ^[11] Python removes repeats with a set and sorts by the floating-point value $u/v$. The quotients are rounded, but the order is still exact: distinct candidates differ by at least $1/26240^2\approx1.5\cdot10^{-9}$, while neighbouring doubles between $1$ and $5248$ are at most $2^{-40}\approx9.1\cdot10^{-13}$ apart. Java's division also truncates towards zero, so floorDiv and ceilDiv are written out as in C++; its Fraction reduces itself in the constructor and compares by cross-multiplying in long, and there is no 128-bit guard. Plain long arithmetic is safe: instrumented copies of the C++ and Java programs show that the coefficients stay below $1.9\cdot10^{12}$ in absolute value and the particular solutions below $8.7\cdot10^{13}$.

Complexity and Verification

The candidate list dominates the work: $3\,358\,720$ gcd computations on numbers up to $26240$ and a sort of the $1\,681\,600$ fractions above $1$, or of the $1\,021\,210$ distinct ones in Python, in $O(P\log P)$ time for $P$ fractions. Only $7040$ candidates pass the bounds of Section 5, the loop over $s_0$ runs $56\,614$ times in all, at most $640$ times per candidate, and each run costs one extended Euclidean algorithm of logarithmic length, so the whole run takes $O(Q\log Q)$ time for the $Q=3\,358\,720$ caviar pairs. The C++ program reserves about 54 MB for the list. On the arm64 computer used for this article, the C++ and Java programs take about $0.4$ seconds and the Python program about $2$ seconds.

The C++ checkpoints pass, and all three programs print the same fraction. We also ran checks that are not part of the published solutions:

A separate test examines the $7040$ candidates with a double loop over $s_0$ and $s_2$, solving for $s_1$ by one division instead of the extended Euclidean algorithm. It finds the same $35$ values and the same largest one, and all three coefficients negative for the $2610$ candidates above the bound of Section 3.
For each of the $35$ values, a table of spoiled items built as in Section 11 satisfies the six conditions of Section 2 in exact rational arithmetic.
An exhaustive search over every spoiled count of every product, which collects all reachable values of $v\,S_B\sum x_i-u\,S_A\sum y_i$ product by product, agrees with Sections 5 to 10 on all $20736$ two-product tables with entries up to $12$, of which $66$ admit a value of $m$. It also agrees on the hamper table divided by $32$, $16$ and $8$, which admit $1$, $1$ and $4$ values.

Footnotes and References

The references below supply the background facts used above. The weighted-average bound for $m$, the reduction to three bounded group sums and the single linear equation are derived explicitly in this article.

Project Euler 236 — Luxury Hampers. The official statement gives the five products supplied by A and B and says that every per-product spoilage rate is worse for B by the same factor m > 1, while the overall rate is worse for A by the factor m. There are thirty-five such m, the smallest being 1476/1475, and the task asks for the largest m as a fraction u/v in lowest terms.
Wikipedia — Simpson's paradox. A trend appears in several groups of data but disappears or reverses when the groups are combined.
Wikipedia — Mediant (mathematics). The mediant of two different fractions with positive denominators lies strictly between them; the notion extends to n fractions with positive weights, and the mediant is a weighted arithmetic mean.
Wikipedia — Euclid's lemma. If an integer n divides a product ab and is coprime with a, then n divides b.
Wikipedia — Bézout's identity. For integers a and b with greatest common divisor d, the integers of the form az + bt are exactly the multiples of d, and from one solution (x, y) of ax + by = d all others are obtained as (x − kb/d, y + ka/d).
Wikipedia — Extended Euclidean algorithm. Besides gcd(a, b), the algorithm computes coefficients x and y with ax + by = gcd(a, b).
Wikipedia — Floor and ceiling functions. Negating the argument switches floor and ceiling and changes the sign, and truncation differs from the floor for negative numbers.
C++ — Euler236.cpp. The linked immutable C++ revision contains generate_candidates, candidate_works, exists_bounded_linear, extended_gcd, tighten_range and run_checkpoints. It is the implementation record for this article.
GCC manual — 128-bit Integers. The types __int128 and unsigned __int128 exist on targets with a 128-bit integer mode.
cppreference — Arithmetic operators. For integer operands the quotient is truncated towards zero, and (a / b) * b + a % b equals a whenever a / b is representable.
Python — Binary arithmetic operations: floor division of integers applies the floor function to the quotient, and x == (x//y)*y + (x%y). The Java Language Specification, §15.17.2 and §15.17.3: integer division rounds towards 0, and (a/b)*b+(a%b) is equal to a.

Problem 236 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <string>
#include <thread>
#include <vector>

namespace {

using i64 = std::int64_t;
using i128 = __int128_t;

struct Fraction {
    i64 u = 0;
    i64 v = 1;
};

struct Options {
    int threads = 0;
    bool run_checkpoints = true;
};

constexpr int NUM_PRODUCTS = 5;
constexpr std::array<i64, NUM_PRODUCTS> A = {5248, 1312, 2624, 5760, 3936};
constexpr std::array<i64, NUM_PRODUCTS> B = {640, 1888, 3776, 3776, 5664};
constexpr int GROUPS = 3;
constexpr std::array<std::array<int, 3>, GROUPS> GROUP_PRODUCTS = {{{0, -1, -1}, {1, 2, 4}, {3, -1, -1}}};
constexpr std::array<int, GROUPS> GROUP_SIZE = {1, 3, 1};
constexpr std::array<i64, GROUPS> GROUP_B_NUM = {5, 59, 59};
constexpr std::array<i64, GROUPS> GROUP_A_DEN = {41, 41, 90};
constexpr i64 SA = 18880;
constexpr i64 SB = 15744;

i64 abs_i64(const i64 x) {
    return x < 0 ? -x : x;
}

Fraction reduce_fraction(i64 u, i64 v) {
    if (v < 0) {
        u = -u;
        v = -v;
    }
    const i64 g = std::gcd(abs_i64(u), abs_i64(v));
    return {u / g, v / g};
}

bool fraction_equal(const Fraction& a, const Fraction& b) {
    return a.u == b.u && a.v == b.v;
}

bool fraction_less(const Fraction& a, const Fraction& b) {
    return static_cast<i128>(a.u) * static_cast<i128>(b.v) < static_cast<i128>(b.u) * static_cast<i128>(a.v);
}

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, "--threads=", options.threads)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.threads >= 0;
}

i64 floor_div(i64 a, i64 b) {
    i64 q = a / b;
    i64 r = a % b;
    if (r != 0 && ((r > 0) != (b > 0))) {
        --q;
    }
    return q;
}

i64 ceil_div(i64 a, i64 b) {
    return -floor_div(-a, b);
}

i64 extended_gcd(i64 a, i64 b, i64& x, i64& y) {
    if (b == 0) {
        x = (a >= 0) ? 1 : -1;
        y = 0;
        return abs_i64(a);
    }
    i64 x1 = 0;
    i64 y1 = 0;
    const i64 g = extended_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

bool tighten_range(i64 p, i64 q, i64 low, i64 high, i64& t_low, i64& t_high) {
    if (low > high) {
        return false;
    }
    if (q == 0) {
        return (p >= low && p <= high);
    }

    i64 l = 0;
    i64 r = 0;
    if (q > 0) {
        l = ceil_div(low - p, q);
        r = floor_div(high - p, q);
    } else {
        l = ceil_div(high - p, q);
        r = floor_div(low - p, q);
    }
    if (l > r) {
        return false;
    }
    if (l > t_low) {
        t_low = l;
    }
    if (r < t_high) {
        t_high = r;
    }
    return t_low <= t_high;
}

bool exists_bounded_linear(i64 c1, i64 l1, i64 h1, i64 c2, i64 l2, i64 h2, i64 rhs) {
    if (l1 > h1 || l2 > h2) {
        return false;
    }
    if (c1 == 0 && c2 == 0) {
        return rhs == 0;
    }
    if (c1 == 0) {
        if (rhs % c2 != 0) {
            return false;
        }
        const i64 s2 = rhs / c2;
        return s2 >= l2 && s2 <= h2;
    }
    if (c2 == 0) {
        if (rhs % c1 != 0) {
            return false;
        }
        const i64 s1 = rhs / c1;
        return s1 >= l1 && s1 <= h1;
    }

    i64 x0 = 0;
    i64 y0 = 0;
    const i64 g = extended_gcd(c1, c2, x0, y0);
    if (rhs % g != 0) {
        return false;
    }

    const i64 scale = rhs / g;
    const i64 step1 = c2 / g;
    const i64 step2 = -c1 / g;

    const i128 base1_128 = static_cast<i128>(x0) * static_cast<i128>(scale);
    const i128 base2_128 = static_cast<i128>(y0) * static_cast<i128>(scale);
    if (base1_128 < static_cast<i128>(std::numeric_limits<i64>::min()) ||
        base1_128 > static_cast<i128>(std::numeric_limits<i64>::max()) ||
        base2_128 < static_cast<i128>(std::numeric_limits<i64>::min()) ||
        base2_128 > static_cast<i128>(std::numeric_limits<i64>::max())) {
        return false;
    }
    const i64 base1 = static_cast<i64>(base1_128);
    const i64 base2 = static_cast<i64>(base2_128);

    i64 t_low = std::numeric_limits<i64>::min() / 4;
    i64 t_high = std::numeric_limits<i64>::max() / 4;
    if (!tighten_range(base1, step1, l1, h1, t_low, t_high)) {
        return false;
    }
    if (!tighten_range(base2, step2, l2, h2, t_low, t_high)) {
        return false;
    }
    return t_low <= t_high;
}

std::vector<Fraction> generate_candidates() {
    std::vector<Fraction> out;
    out.reserve(static_cast<std::size_t>(A[0] * B[0]));

    for (i64 a = 1; a <= A[0]; ++a) {
        for (i64 b = 1; b <= B[0]; ++b) {
            Fraction m = reduce_fraction(41 * b, 5 * a);
            if (m.u > m.v) {
                out.push_back(m);
            }
        }
    }

    std::sort(out.begin(), out.end(), fraction_less);
    out.erase(std::unique(out.begin(), out.end(), fraction_equal), out.end());
    return out;
}

bool candidate_works(const Fraction& m) {
    if (m.u <= m.v) {
        return false;
    }

    std::array<i64, GROUPS> min_sum = {1, 3, 1};
    std::array<i64, GROUPS> max_sum = {0, 0, 0};
    std::array<i64, GROUPS> coef = {0, 0, 0};

    for (int g = 0; g < GROUPS; ++g) {
        const Fraction group_ratio = reduce_fraction(m.u * GROUP_B_NUM[g], m.v * GROUP_A_DEN[g]);
        const i64 N = group_ratio.u;
        const i64 D = group_ratio.v;

        i64 upper_sum = 0;
        for (int j = 0; j < GROUP_SIZE[g]; ++j) {
            const int idx = GROUP_PRODUCTS[g][j];
            const i64 upper = std::min(A[idx] / D, B[idx] / N);
            if (upper < 1) {
                return false;
            }
            upper_sum += upper;
        }
        max_sum[g] = upper_sum;
        coef[g] = m.v * SB * D - m.u * SA * N;
    }

    for (i64 s0 = min_sum[0]; s0 <= max_sum[0]; ++s0) {
        const i64 rhs = -(coef[0] * s0);
        if (exists_bounded_linear(coef[1], min_sum[1], max_sum[1], coef[2], min_sum[2], max_sum[2], rhs)) {
            return true;
        }
    }

    return false;
}

std::vector<Fraction> solve_all(const std::vector<Fraction>& candidates, int threads) {
    if (threads <= 0) {
        threads = static_cast<int>(std::thread::hardware_concurrency());
        if (threads <= 0) {
            threads = 1;
        }
    }

    const int n = static_cast<int>(candidates.size());
    threads = std::max(1, std::min(threads, n));

    std::vector<std::vector<Fraction>> local(static_cast<std::size_t>(threads));
    std::vector<std::thread> workers;
    workers.reserve(static_cast<std::size_t>(threads));

    for (int t = 0; t < threads; ++t) {
        workers.emplace_back([&, t]() {
            for (int i = t; i < n; i += threads) {
                if (candidate_works(candidates[static_cast<std::size_t>(i)])) {
                    local[static_cast<std::size_t>(t)].push_back(candidates[static_cast<std::size_t>(i)]);
                }
            }
        });
    }
    for (std::thread& worker : workers) {
        worker.join();
    }

    std::vector<Fraction> solutions;
    for (auto& chunk : local) {
        solutions.insert(solutions.end(), chunk.begin(), chunk.end());
    }
    std::sort(solutions.begin(), solutions.end(), fraction_less);
    solutions.erase(std::unique(solutions.begin(), solutions.end(), fraction_equal), solutions.end());
    return solutions;
}

bool run_checkpoints(const std::vector<Fraction>& solutions) {
    if (solutions.size() != 35U) {
        std::cerr << "Checkpoint failed for number of valid m values" << '\n';
        return false;
    }
    if (!fraction_equal(solutions.front(), Fraction{1476, 1475})) {
        std::cerr << "Checkpoint failed for smallest valid m value" << '\n';
        return false;
    }
    if (!candidate_works(Fraction{1476, 1475})) {
        std::cerr << "Checkpoint failed for explicit smallest-m verification" << '\n';
        return false;
    }
    return true;
}

}  // namespace

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

    const std::vector<Fraction> candidates = generate_candidates();
    const std::vector<Fraction> solutions = solve_all(candidates, options.threads);

    if (solutions.empty()) {
        std::cerr << "No solution found" << '\n';
        return 2;
    }
    if (options.run_checkpoints && !run_checkpoints(solutions)) {
        return 3;
    }

    const Fraction best = solutions.back();
    std::cout << best.u << '/' << best.v << '\n';
    return 0;
}

Python

import math

NUM_PRODUCTS = 5
A = [5248, 1312, 2624, 5760, 3936]
B = [640, 1888, 3776, 3776, 5664]

GROUPS = 3
GROUP_PRODUCTS = [[0, -1, -1], [1, 2, 4], [3, -1, -1]]
GROUP_SIZE = [1, 3, 1]
GROUP_B_NUM = [5, 59, 59]
GROUP_A_DEN = [41, 41, 90]

SA = 18880
SB = 15744

def reduce_fraction(u, v):
    if v < 0:
        u = -u
        v = -v
    g = math.gcd(abs(u), abs(v))
    return (u // g, v // g)

def floor_div(a, b):
    return a // b

def ceil_div(a, b):
    return -((-a) // b)

def extended_gcd(a, b):
    if b == 0:
        return (1 if a >= 0 else -1, 0, abs(a))
    x1, y1, g = extended_gcd(b, a % b)
    x = y1
    y = x1 - (a // b) * y1
    return (x, y, g)

def tighten_range(p, q, low, high, t_low, t_high):
    if low > high:
        return False, t_low, t_high
    if q == 0:
        if p >= low and p <= high:
            return True, t_low, t_high
        return False, t_low, t_high

    if q > 0:
        l = ceil_div(low - p, q)
        r = floor_div(high - p, q)
    else:
        l = ceil_div(high - p, q)
        r = floor_div(low - p, q)

    if l > r:
        return False, t_low, t_high
        
    if l > t_low:
        t_low = l
    if r < t_high:
        t_high = r
        
    return t_low <= t_high, t_low, t_high

def exists_bounded_linear(c1, l1, h1, c2, l2, h2, rhs):
    if l1 > h1 or l2 > h2:
        return False
    if c1 == 0 and c2 == 0:
        return rhs == 0
    if c1 == 0:
        if rhs % c2 != 0:
            return False
        s2 = rhs // c2
        return l2 <= s2 <= h2
    if c2 == 0:
        if rhs % c1 != 0:
            return False
        s1 = rhs // c1
        return l1 <= s1 <= h1

    x0, y0, g = extended_gcd(c1, c2)
    if rhs % g != 0:
        return False

    scale = rhs // g
    step1 = c2 // g
    step2 = -c1 // g

    base1 = x0 * scale
    base2 = y0 * scale

    t_low = -10**20
    t_high = 10**20

    success, t_low, t_high = tighten_range(base1, step1, l1, h1, t_low, t_high)
    if not success:
        return False
    success, t_low, t_high = tighten_range(base2, step2, l2, h2, t_low, t_high)
    if not success:
        return False
        
    return t_low <= t_high

def generate_candidates():
    candidates = []
    for a in range(1, A[0] + 1):
        for b in range(1, B[0] + 1):
            m_u, m_v = reduce_fraction(41 * b, 5 * a)
            if m_u > m_v:
                candidates.append((m_u, m_v))
    
    # Sort correctly by fraction value
    def fraction_key(m):
        return m[0] / m[1]
        
    candidates = sorted(list(set(candidates)), key=fraction_key)
    return candidates

def candidate_works(m):
    m_u, m_v = m
    if m_u <= m_v:
        return False

    min_sum = [1, 3, 1]
    max_sum = [0, 0, 0]
    coef = [0, 0, 0]

    for g in range(GROUPS):
        N, D = reduce_fraction(m_u * GROUP_B_NUM[g], m_v * GROUP_A_DEN[g])

        upper_sum = 0
        for j in range(GROUP_SIZE[g]):
            idx = GROUP_PRODUCTS[g][j]
            upper = min(A[idx] // D, B[idx] // N)
            if upper < 1:
                return False
            upper_sum += upper
            
        max_sum[g] = upper_sum
        coef[g] = m_v * SB * D - m_u * SA * N

    for s0 in range(min_sum[0], max_sum[0] + 1):
        rhs = - (coef[0] * s0)
        if exists_bounded_linear(coef[1], min_sum[1], max_sum[1], coef[2], min_sum[2], max_sum[2], rhs):
            return True

    return False

def solve():
    candidates = generate_candidates()
    solutions = []
    
    for cand in candidates:
        if candidate_works(cand):
            solutions.append(cand)
            
    best = solutions[-1]
    return f"{best[0]}/{best[1]}"

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

Java

import java.util.*;

public class Euler236 {
    static final long[] A = { 5248, 1312, 2624, 5760, 3936 };
    static final long[] B = { 640, 1888, 3776, 3776, 5664 };

    static final int GROUPS = 3;
    static final int[][] GROUP_PRODUCTS = {
            { 0, -1, -1 },
            { 1, 2, 4 },
            { 3, -1, -1 }
    };
    static final int[] GROUP_SIZE = { 1, 3, 1 };
    static final long[] GROUP_B_NUM = { 5, 59, 59 };
    static final long[] GROUP_A_DEN = { 41, 41, 90 };

    static final long SA = 18880;
    static final long SB = 15744;

    static long gcd(long a, long b) {
        if (b == 0)
            return Math.abs(a);
        return gcd(b, a % b);
    }

    static class Fraction implements Comparable<Fraction> {
        long u, v;

        Fraction(long u, long v) {
            if (v < 0) {
                u = -u;
                v = -v;
            }
            long g = gcd(u, v);
            this.u = u / g;
            this.v = v / g;
        }

        @Override
        public int compareTo(Fraction o) {
            // u/v < o.u/o.v -> u*o.v < o.u*v
            // use BigInteger or double if overflow can happen, but here it fits in long?
            // Max u,v are around 5000. So 5000*5000 fits easily in long.
            long diff = this.u * o.v - o.u * this.v;
            return Long.compare(diff, 0);
        }

        @Override
        public boolean equals(Object obj) {
            if (!(obj instanceof Fraction))
                return false;
            Fraction o = (Fraction) obj;
            return u == o.u && v == o.v;
        }
    }

    static long floorDiv(long a, long b) {
        long q = a / b;
        long r = a % b;
        if (r != 0 && ((r > 0) != (b > 0))) {
            q--;
        }
        return q;
    }

    static long ceilDiv(long a, long b) {
        return -floorDiv(-a, b);
    }

    static long[] extendedGcd(long a, long b) {
        if (b == 0) {
            return new long[] { a >= 0 ? 1 : -1, 0, Math.abs(a) };
        }
        long[] res = extendedGcd(b, a % b);
        long x1 = res[0];
        long y1 = res[1];
        long g = res[2];
        return new long[] { y1, x1 - (a / b) * y1, g };
    }

    static long[] tightenRange(long p, long q, long low, long high, long tLow, long tHigh) {
        if (low > high)
            return new long[] { 0, tLow, tHigh };
        if (q == 0) {
            if (p >= low && p <= high)
                return new long[] { 1, tLow, tHigh };
            return new long[] { 0, tLow, tHigh };
        }

        long l = 0, r = 0;
        if (q > 0) {
            l = ceilDiv(low - p, q);
            r = floorDiv(high - p, q);
        } else {
            l = ceilDiv(high - p, q);
            r = floorDiv(low - p, q);
        }

        if (l > r)
            return new long[] { 0, tLow, tHigh };
        if (l > tLow)
            tLow = l;
        if (r < tHigh)
            tHigh = r;

        return new long[] { tLow <= tHigh ? 1 : 0, tLow, tHigh };
    }

    static boolean existsBoundedLinear(long c1, long l1, long h1, long c2, long l2, long h2, long rhs) {
        if (l1 > h1 || l2 > h2)
            return false;
        if (c1 == 0 && c2 == 0)
            return rhs == 0;
        if (c1 == 0) {
            if (rhs % c2 != 0)
                return false;
            long s2 = rhs / c2;
            return s2 >= l2 && s2 <= h2;
        }
        if (c2 == 0) {
            if (rhs % c1 != 0)
                return false;
            long s1 = rhs / c1;
            return s1 >= l1 && s1 <= h1;
        }

        long[] res = extendedGcd(c1, c2);
        long x0 = res[0];
        long y0 = res[1];
        long g = res[2];

        if (rhs % g != 0)
            return false;

        long scale = rhs / g;
        long step1 = c2 / g;
        long step2 = -c1 / g;

        // Ensure no overflow
        long base1 = x0 * scale;
        long base2 = y0 * scale;

        long tLow = Long.MIN_VALUE / 4;
        long tHigh = Long.MAX_VALUE / 4;

        long[] tr1 = tightenRange(base1, step1, l1, h1, tLow, tHigh);
        if (tr1[0] == 0)
            return false;
        tLow = tr1[1];
        tHigh = tr1[2];

        long[] tr2 = tightenRange(base2, step2, l2, h2, tLow, tHigh);
        if (tr2[0] == 0)
            return false;
        tLow = tr2[1];
        tHigh = tr2[2];

        return tLow <= tHigh;
    }

    static Fraction[] generateCandidates() {
        List<Fraction> out = new ArrayList<>();
        for (long a = 1; a <= A[0]; ++a) {
            for (long b = 1; b <= B[0]; ++b) {
                Fraction m = new Fraction(41 * b, 5 * a);
                if (m.u > m.v) {
                    out.add(m);
                }
            }
        }
        Collections.sort(out);
        List<Fraction> unique = new ArrayList<>();
        for (int i = 0; i < out.size(); ++i) {
            if (i == 0 || !out.get(i).equals(out.get(i - 1))) {
                unique.add(out.get(i));
            }
        }
        return unique.toArray(new Fraction[0]);
    }

    static boolean candidateWorks(Fraction m) {
        if (m.u <= m.v)
            return false;

        long[] minSum = { 1, 3, 1 };
        long[] maxSum = { 0, 0, 0 };
        long[] coef = { 0, 0, 0 };

        for (int g = 0; g < GROUPS; ++g) {
            Fraction groupRatio = new Fraction(m.u * GROUP_B_NUM[g], m.v * GROUP_A_DEN[g]);
            long N = groupRatio.u;
            long D = groupRatio.v;

            long upperSum = 0;
            for (int j = 0; j < GROUP_SIZE[g]; ++j) {
                int idx = GROUP_PRODUCTS[g][j];
                long upper = Math.min(A[idx] / D, B[idx] / N);
                if (upper < 1)
                    return false;
                upperSum += upper;
            }
            maxSum[g] = upperSum;
            coef[g] = m.v * SB * D - m.u * SA * N;
        }

        for (long s0 = minSum[0]; s0 <= maxSum[0]; ++s0) {
            long rhs = -(coef[0] * s0);
            if (existsBoundedLinear(coef[1], minSum[1], maxSum[1], coef[2], minSum[2], maxSum[2], rhs)) {
                return true;
            }
        }

        return false;
    }

    public static String solve() {
        Fraction[] candidates = generateCandidates();
        List<Fraction> solutions = new ArrayList<>();

        // Multi-threading could be used, but since we're porting directly and Python
        // runs it fairly fast
        for (Fraction cand : candidates) {
            if (candidateWorks(cand)) {
                solutions.add(cand);
            }
        }

        if (solutions.isEmpty()) {
            return "No solution";
        }

        Fraction best = solutions.get(solutions.size() - 1);
        return best.u + "/" + best.v;
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}