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

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$.

Why Product 1 Generates Every Candidate

From the first product we have $P_1=5248$ and $Q_1=640$, so

$$\frac{u}{v}=\frac{P_1q_1}{Q_1p_1}=\frac{5248\,q_1}{640\,p_1}=\frac{41q_1}{5p_1}.$$

Because $1\le p_1\le P_1$ and $1\le q_1\le Q_1$, every valid multiplier must appear among the reduced fractions

$$m=\frac{41b}{5a},\qquad 1\le a\le 5248,\quad 1\le b\le 640,\quad m>1.$$

This is the key finiteness argument: the search never needs to guess arbitrary rationals. It only needs to enumerate, reduce, sort, and deduplicate the fractions forced by product 1.

Five Products Collapse to Three Ratio Groups

The five supplier ratios are not all different:

$$\frac{Q_1}{P_1}=\frac{5}{41},\qquad \frac{Q_2}{P_2}=\frac{Q_3}{P_3}=\frac{Q_5}{P_5}=\frac{59}{41},\qquad \frac{Q_4}{P_4}=\frac{59}{90}.$$

So the products split naturally into three groups:

$$G_0=\{1\},\qquad G_1=\{2,3,5\},\qquad G_2=\{4\}.$$

For a fixed candidate $m=\frac{u}{v}$, let the base ratio of group $g$ be $\frac{\beta_g}{\alpha_g}$, so the three values are $\frac{5}{41}$, $\frac{59}{41}$, and $\frac{59}{90}$. Reduce

$$\frac{u\beta_g}{v\alpha_g}=\frac{n_g}{d_g}$$

to lowest terms. Then every product $i\in G_g$ must satisfy

$$\frac{q_i}{p_i}=\frac{uQ_i}{vP_i}=\frac{n_g}{d_g},$$

hence there is a positive integer $t_i$ such that

$$p_i=d_gt_i,\qquad q_i=n_gt_i.$$

This is the central simplification: once $m$ is fixed, each product contributes only one new integer parameter $t_i$, and products that share the same supplier ratio also share the same pair $(n_g,d_g)$.

Bounds and the Group-Sum Interval

The chosen counts cannot exceed the supplier totals, so for every product $i\in G_g$,

$$1\le t_i\le \min\left(\left\lfloor\frac{P_i}{d_g}\right\rfloor,\left\lfloor\frac{Q_i}{n_g}\right\rfloor\right).$$

Now define one sum per group:

$$s_g=\sum_{i\in G_g}t_i.$$

If a group contains $|G_g|$ products, then the smallest possible value is

$$L_g=|G_g|,$$

because each $t_i\ge 1$. The largest possible value is

$$H_g=\sum_{i\in G_g}\min\left(\left\lfloor\frac{P_i}{d_g}\right\rfloor,\left\lfloor\frac{Q_i}{n_g}\right\rfloor\right).$$

A useful observation is that every integer between $L_g$ and $H_g$ is achievable. There is no subset-sum complication here: all $t_i$ move in steps of 1, so starting from the minimum assignment $t_i=1$, we can increase variables one by one until any target group sum in the interval is reached.

Thus the five original parameters $t_1,\dots,t_5$ can be replaced by the three interval-constrained sums

$$s_0\in[L_0,H_0],\qquad s_1\in[L_1,H_1],\qquad s_2\in[L_2,H_2],$$

with $L_0=1$, $L_1=3$, and $L_2=1$.

The Final Feasibility Test Is a Bounded Linear Diophantine Equation

Substituting $p_i=d_gt_i$ and $q_i=n_gt_i$ into the global condition gives

$$vT_Q\sum_{g=0}^{2}d_gs_g=uT_P\sum_{g=0}^{2}n_gs_g.$$

Move everything to one side and define

$$c_g=vT_Qd_g-uT_Pn_g.$$

Then $m$ is feasible if and only if there exist integers $s_0,s_1,s_2$ in their intervals such that

$$c_0s_0+c_1s_1+c_2s_2=0.$$

The implementation fixes $s_0$ and rewrites the problem as

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

If $g=\gcd(c_1,c_2)$ does not divide $r$, there is no solution. Otherwise, if

$$c_1x_0+c_2y_0=g,$$

then every solution is of the form

$$s_1=x_0\frac{r}{g}+\frac{c_2}{g}t,\qquad s_2=y_0\frac{r}{g}-\frac{c_1}{g}t,\qquad t\in\mathbb{Z}.$$

The remaining task is purely interval arithmetic: intersect the range of $t$ forced by $L_1\le s_1\le H_1$ with the range forced by $L_2\le s_2\le H_2$. A non-empty intersection means the candidate $m$ works.

Worked Example: The Maximal Candidate

The implementations ultimately find that the largest feasible multiplier is

$$m=\frac{123}{59}.$$

For this value, the three group reductions are

$$\frac{123\cdot 5}{59\cdot 41}=\frac{15}{59},\qquad \frac{123\cdot 59}{59\cdot 41}=3=\frac{3}{1},\qquad \frac{123\cdot 59}{59\cdot 90}=\frac{41}{30}.$$

So the group parameters are

$$ (d_0,n_0)=(59,15),\qquad (d_1,n_1)=(1,3),\qquad (d_2,n_2)=(30,41). $$

The corresponding upper bounds are

$$H_0=\min\!\left(\left\lfloor\frac{5248}{59}\right\rfloor,\left\lfloor\frac{640}{15}\right\rfloor\right)=42,$$

$$H_1=\min\!\left(\left\lfloor\frac{1312}{1}\right\rfloor,\left\lfloor\frac{1888}{3}\right\rfloor\right)+\min\!\left(\left\lfloor\frac{2624}{1}\right\rfloor,\left\lfloor\frac{3776}{3}\right\rfloor\right)+\min\!\left(\left\lfloor\frac{3936}{1}\right\rfloor,\left\lfloor\frac{5664}{3}\right\rfloor\right)=629+1258+1888=3775,$$

$$H_2=\min\!\left(\left\lfloor\frac{5760}{30}\right\rfloor,\left\lfloor\frac{3776}{41}\right\rfloor\right)=92.$$

The coefficients become

$$ (c_0,c_1,c_2)=(19971264,-6037824,-67344960)=464448\,(43,-13,-145). $$

So feasibility is equivalent to

$$43s_0-13s_1-145s_2=0.$$

One valid choice is

$$ (s_0,s_1,s_2)=(7,12,1), $$

because $43\cdot 7-13\cdot 12-145\cdot 1=0$. The middle-group sum $s_1=12$ is easy to realize, for example by taking the three product parameters there as $1,1,10$. This concrete construction shows why $\frac{123}{59}$ passes the feasibility test.

How the Code Works

The C++, Python, and Java implementations follow the derivation directly. They never enumerate arbitrary five-tuples of product counts. Instead, they search the rational candidate set forced by the first product and test each candidate with exact integer arithmetic.

Candidate Enumeration

The implementation loops over all possible positive counts for product 1, forms the reduced fraction $\frac{41b}{5a}$, discards candidates with $m\le 1$, sorts the remaining fractions, and removes duplicates. This step constructs the complete finite search space.

Feasibility Test for One Candidate

For a fixed candidate $m=\frac{u}{v}$, the implementation computes the three reduced group ratios $\frac{n_g}{d_g}$, then the upper bound for each product and the interval $[L_g,H_g]$ for each group sum $s_g$. If any product has no positive admissible value, the candidate is rejected immediately.

Otherwise it forms the coefficients $c_g=vT_Qd_g-uT_Pn_g$, loops over all admissible $s_0$, and for each one checks whether the remaining two-variable equation has a solution inside the intervals. That check uses the extended Euclidean algorithm together with interval tightening for the free parameter $t$.

Collecting the Answer

Each candidate can be tested independently. The C++ and Java implementations distribute those tests across multiple threads; the Python implementation performs the same arithmetic serially. After all feasible candidates are collected, the largest reduced fraction is returned.

Complexity Analysis

Candidate generation starts from the $5248\cdot 640=3{,}358{,}720$ raw fractions forced by product 1. Sorting and deduplicating them costs $O(M\log M)$ with $M=3{,}358{,}720$.

For each distinct candidate, the group preprocessing is constant-time, and the main loop scans only the possible values of $s_0$. Since the first group contains a single product, $s_0$ never ranges over more than a few hundred values. The two-variable Diophantine check is $O(1)$ once the gcd and parameter intervals are computed. So the practical running time is dominated by candidate enumeration plus an $O(C\cdot R)$ feasibility phase, where $C$ is the number of distinct candidates and $R$ is the width of the first group's interval. Memory usage is $O(C)$ for storing the candidate list, plus constant auxiliary state for each test. Parallelism improves wall-clock time but does not change the asymptotic bounds.

Footnotes and References

Problem page: https://projecteuler.net/problem=236
Rational number: Wikipedia - Rational number
Linear Diophantine equation: Wikipedia - Linear Diophantine equation
Bézout's identity: Wikipedia - Bézout's identity
Extended Euclidean algorithm: Wikipedia - Extended Euclidean algorithm
Greatest common divisor: Wikipedia - Greatest common divisor

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());
    }
}