Problem 190: Maximising a Weighted Product

Project Euler Problem 190 Solution

EulerSolve provides an optimized solution for Project Euler Problem 190, Maximising a Weighted Product, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer $m$ from 2 through 15, define $$P_m=\max \left\{x_1^1x_2^2\cdots x_m^m : x_1+x_2+\cdots+x_m=m,\ x_i \gt 0\right\}.$$ The exponents $1,2,\dots,m$ make the variables asymmetric: increasing a large-index variable matters much more than increasing a small-index one. The task is to determine the maximizing configuration for every $m$, take $\lfloor P_m\rfloor$, and add those integer parts. Mathematical Approach This problem looks like a many-variable nonlinear search, but the structure is very rigid. The maximizing point is forced to be proportional to the weight vector $(1,2,\dots,m)$, so the entire problem collapses to a closed formula. Take logarithms and move to a concave objective Because the logarithm is strictly increasing, maximizing $P_m$ is equivalent to maximizing $$F(x_1,\dots,x_m)=\log P_m=\sum_{i=1}^{m} i\log x_i$$ under the same linear constraint $x_1+\cdots+x_m=m$. This reformulation is useful for two reasons. First, it turns a product into a sum, so differentiation becomes straightforward. Second, the Hessian is diagonal: $$\frac{\partial^2 F}{\partial x_i^2}=-\frac{i}{x_i^2},\qquad \frac{\partial^2 F}{\partial x_i\partial x_j}=0 \quad (i\ne j).$$ Every diagonal entry is negative on the feasible region $x_i \gt 0$, so $F$ is strictly concave....

Detailed mathematical approach

Problem Summary

For a whole number $m\ge2$ we consider $m$ positive real numbers $x_1,x_2,\dots,x_m$ with the sum $x_1+x_2+\dots+x_m=m$ and the weighted product $x_1\cdot x_2^2\cdots x_m^m$, in which the $i$-th number carries the exponent $i$. The $m$-tuple $S_m=(x_1,\dots,x_m)$ that makes this product as large as possible gives the maximal value $P_m$. With $\lfloor\,\rfloor$ the integer part, the statement gives $\lfloor P_{10}\rfloor=4112$, and we must find the sum of the integer parts $\lfloor P_m\rfloor$ for $m=2,3,\dots,15$. ^[1]

The whole solution in one sentence: the product is largest when every $x_i$ is proportional to its exponent, $x_i=2i/(m+1)$, as the inequality of arithmetic and geometric means proves, so $P_m$ is an explicit product that the programs evaluate through its logarithm for $m=2,\dots,15$ before they add the integer parts.

Mathematical Approach

We solve the smallest case by hand, guess the general maximiser with a balancing argument, prove it with the inequality of arithmetic and geometric means, and finally check that floating-point numbers are accurate enough for the integer parts.

1. The case $m=2$ by hand

For $m=2$ the numbers satisfy $x_1+x_2=2$. Write $x_2=t$ and $x_1=2-t$ with $0\lt t\lt2$; the product is $f(t)=(2-t)\,t^2$. A few values: $f(1)=1$, $f(1.2)=1.152$ and $f(1.5)=1.125$. The derivative $f'(t)=4t-3t^2=t(4-3t)$ is positive for $t\lt4/3$ and negative for $t\gt4/3$, so $f$ rises up to $t=4/3$ and falls after it. Hence $S_2=(2/3,\,4/3)$ and

$$P_2=\frac23\cdot\Bigl(\frac43\Bigr)^2=\frac{32}{27}\approx1.185,$$

so $\lfloor P_2\rfloor=1$. Notice that $x_2=2x_1$: the number with the exponent $2$ receives twice as much as the number with the exponent $1$.

2. Guessing the maximiser by balancing

Logarithms turn the product into a sum: $\ln(x_1x_2^2\cdots x_m^m)$ equals $\sum_i i\ln x_i$, and since $\ln$ is increasing, making the sum largest is the same as making the product largest. Now move a small amount $\delta$ from $x_j$ to $x_i$; the total stays $m$. Because the derivative of $\ln x$ is $1/x$, the term $i\ln x_i$ grows by about $\delta\,i/x_i$, and $j\ln x_j$ shrinks by about $\delta\,j/x_j$. So the move pays off when $i/x_i\gt j/x_j$: the ratio $i/x_i$ is the gain per unit of the $i$-th number. At the point $(1,1)$ for $m=2$ the gains are $1$ and $2$, so some of $x_1$ should go to $x_2$; at $(2/3,\,4/3)$ both gains equal $3/2$, and no small move helps.

At a maximum no move helps, so all gains $i/x_i$ are equal, say to $\lambda$, and $x_i=i/\lambda$. The condition $x_1+\dots+x_m=m$ becomes $(1+2+\dots+m)/\lambda=m$, and since $1+2+\dots+m=m(m+1)/2$, we get $\lambda=(m+1)/2$ and

$$x_i=\frac{2i}{m+1},\qquad i=1,2,\dots,m.$$

The method of Lagrange multipliers reaches the same point: it sets the partial derivatives $i/x_i$ of $\sum_k k\ln x_k$ equal to $\lambda$ times those of $\sum_k x_k$, which are all $1$. ^[2] Like that method, the balancing argument only says where a maximum must lie if there is one; it does not show that the product has a maximum at all. The next section proves both at once.

3. Proof with the inequality of arithmetic and geometric means

The inequality of arithmetic and geometric means, AM–GM for short, says that for nonnegative numbers $a_1,\dots,a_n$

$$\frac{a_1+a_2+\dots+a_n}{n}\ge\sqrt[n]{a_1a_2\cdots a_n},$$

with equality only when all $a_k$ are equal. ^[3] For $m=2$, split $x_2$ into two halves: the three numbers $x_1$, $x_2/2$ and $x_2/2$ have the sum $2$, so AM–GM gives $(2/3)^3\ge x_1(x_2/2)^2$, that is $x_1x_2^2\le4\cdot(2/3)^3=32/27$, with equality only for $x_1=x_2/2$, at $(2/3,\,4/3)$. This is the result of Section 1 without any derivative.

In general, split every $x_i$ into $i$ equal parts $x_i/i$. There are $T=1+2+\dots+m=m(m+1)/2$ parts, their sum is $x_1+\dots+x_m=m$, and so their mean is $m/T=2/(m+1)$. AM–GM, raised to the power $T$, gives

$$\prod_{i=1}^{m}\Bigl(\frac{x_i}{i}\Bigr)^{i}\le\Bigl(\frac{2}{m+1}\Bigr)^{T},$$

and multiplying by $1^1\cdot2^2\cdots m^m$ shows that every admissible $m$-tuple satisfies

$$x_1x_2^2\cdots x_m^m\le\prod_{i=1}^{m}\Bigl(\frac{2i}{m+1}\Bigr)^{i}.$$

Equality holds exactly when all parts are equal, $x_i/i=2/(m+1)$. So the bound is reached, and at one point only: the maximiser is the point of Section 2,

$$S_m=\Bigl(\frac{2}{m+1},\frac{4}{m+1},\dots,\frac{2m}{m+1}\Bigr),$$

and $P_m$ equals the right-hand side of the bound. In the language of the weighted AM–GM inequality, the numbers $x_i/i$ carry the weights $i$. ^[3]

4. Worked examples

$m=3$: $S_3=(1/2,\,1,\,3/2)$ and $P_3=\tfrac12\cdot1^2\cdot(3/2)^3$${}=27/16=1.6875$, so $\lfloor P_3\rfloor=1$. The nearby admissible point $(0.6,\,1,\,1.4)$ gives only $0.6\cdot1.4^3=1.6464$.
$m=4$: $S_4=(2/5,\,4/5,\,6/5,\,8/5)$ and $P_4=28311552/9765625$${}=2.8991029248$, so $\lfloor P_4\rfloor=2$.
$m=5$: $S_5=(1/3,\,2/3,\,1,\,4/3,\,5/3)$ and $P_5=3200000/531441\approx6.0214$, so $\lfloor P_5\rfloor=6$.
$m=10$, the checkpoint: $S_{10}=(2/11,\,4/11,\,\dots,\,20/11)$ and $P_{10}\approx4112.0850$, so $\lfloor P_{10}\rfloor=4112$, as in the statement.

The point $(1,1,\dots,1)$ is admissible and gives the product $1$. For $m\ge2$ it is not the maximiser, so $P_m\gt1$ and every integer part is at least $1$.

5. An exact fraction

The $T$ factors $2$ and the $T$ denominators $m+1$ can be collected. With the hyperfactorial $H(m)=1^1\cdot2^2\cdots m^m$, whose values for $m=1,\dots,5$ are $1$, $4$, $108$, $27648$ and $86400000$, ^[4]

$$P_m=\frac{2^{T}\,H(m)}{(m+1)^{T}},\qquad T=\frac{m(m+1)}{2}.$$

So $P_m$ is a fraction of two whole numbers, and $\lfloor P_m\rfloor$ is the quotient of an integer division. For $m=4$: $T=10$, $2^{10}\cdot27648=28311552$ and $5^{10}=9765625$, with the quotient $2$ of Section 4. With big integers this gives every $\lfloor P_m\rfloor$ exactly, as in the verification below; the published programs use floating point instead, which the next sections justify.

6. How fast $P_m$ grows

For $m=2,3,\dots,10$ the integer parts are $1, 1, 2, 6, 15, 46, 169, 759$ and $4112$. In $\ln P_m=\sum_i i\ln\bigl(2i/(m+1)\bigr)$ a term is negative for $2i\lt m+1$, zero for $2i=m+1$ and positive for $2i\gt m+1$; the total is positive, since $P_m\gt1$ by Section 4. The hyperfactorial grows like

$$H(n)\approx A\,n^{n^2/2+n/2+1/12}\,e^{-n^2/4},$$

where $A\approx1.28243$ is the Glaisher–Kinkelin constant; the relative error is of order $1/n^2$. ^[4] Here $O(1/m)$ is a term at most a constant times $1/m$ in size. From $\ln(1+1/m)$${}=1/m-1/(2m^2)+O(1/m^3)$ we get $T\ln(m+1)$${}=T\ln m+m/2+1/4+O(1/m)$, and the formula of Section 5, $\ln P_m=T\ln2+\ln H(m)$${}-T\ln(m+1)$, becomes

$$\begin{aligned}\ln P_m&=T\Bigl(\ln2-\frac12\Bigr)-\frac m4+\frac{\ln m}{12}\\&\quad+\ln A-\frac14+O\Bigl(\frac1m\Bigr).\end{aligned}$$

Since $\ln2-\tfrac12\approx0.193$, the logarithm grows like $0.0966\,m^2$: the number of digits of $P_m$ grows quadratically, and $P_{20}$ already has $16$ digits before the decimal point. For $m=10$, the formula without its error term, exponentiated, gives $4079.5$ instead of $P_{10}\approx4112.09$; the omitted term in the logarithm is $0.0080$ here, close to $1/(12m)$.

7. Computing through logarithms

The programs evaluate

$$L_m=\sum_{i=1}^{m}i\ln\frac{2i}{m+1},\qquad P_m=e^{L_m}.$$

Logarithms are the usual safeguard against overflow and underflow in long products. Here the values stay moderate either way: for $m=15$ the smallest partial product of the powers is about $1.7\cdot10^{-7}$, and a plain product would be about as accurate. The price of the logarithm is that an absolute error $\varepsilon$ in $L_m$ becomes a relative error of about $\varepsilon$ in $P_m$, because $e^{L+\varepsilon}=e^L(1+\varepsilon+\dots)$. A double, the IEEE-754 binary64 format, has $53$ significant bits, about $16$ decimal digits, and rounding a number of normal size to the nearest double changes it by a relative error of at most $2^{-53}\approx1.1\cdot10^{-16}$. ^[5] The logarithms and the exponential add errors of their own, which depend on the compiler and the maths library. On the machine used for this article, the computed $L_m$ of all three programs differ from the exact values by less than $5\cdot10^{-15}$ for $2\le m\le15$, and the computed $P_m$ by less than $1.6\cdot10^{-6}$, the largest errors occurring in Java.

8. When the integer part is safe

The integer part of a computed value is correct if the exact value lies farther from the nearest integer than the error of the computation. The exact fractions of Section 5 give these distances: for $2\le m\le15$ the smallest fractional part of $P_m$ is about $0.0214$, for $m=5$, and the largest about $0.9952$, for $m=12$. The closest approach to an integer is therefore about $0.0048$, at $m=12$, more than three thousand times the largest error of Section 7, so every computed integer part is exact.

The programs add $10^{-12}$ before they cut off the fractional part. Such a guard repairs a value that ought to be an integer $k$ but comes out as $k-10^{-15}$, say, but it could also push a value lying just below an integer over it. For $2\le m\le15$ neither can happen: no $P_m$ is an integer, and every one lies at least $0.0047$ from the nearest integer, far more than the guard and the rounding errors together. In double precision the guard does not even change values from $2^{14}=16384$ on, where it is less than half the gap between neighbouring doubles.

The method has a limit. The C++ program accepts every range with $2\le$ from $\le$ to, however large. In the tested builds with $53$-bit doubles its integer parts agree with the exact ones up to $m=19$ and are wrong for $m=20$, $21$ and $22$. At $m=20$ every integer near $P_{20}\approx3.58\cdot10^{15}$ is still a double, but the computed values are off by $18$ to $44$ units, relative errors of $5\cdot10^{-15}$ to $1.2\cdot10^{-14}$. An emulation with a $64$-bit significand, the precision of the 80-bit x87 format, agrees up to $m=21$ and fails at $m=22$, and from $m=23$ on $P_m$ exceeds $2^{63}-1$ and no longer fits into a $64$-bit integer.

How the Code Works

The C++ program ^[6] keeps its settings in the struct Options: the first and the last value of $m$, from $=2$ and to $=15$, and the flag run_checkpoints. parse_arguments accepts --skip-checkpoints, --from= and --to=. For the two numbers, parse_int_after_prefix requires a nonempty string of decimal digits and collects them in an int without an overflow check, so an oversized number causes signed overflow, undefined behaviour in C++, before the range is checked. An unknown argument is reported on the error stream and, like a range with from below $2$ or to below from, ends the program with exit code $1$. Unless skipped, run_checkpoints requires P_floor(10) to return $4112$, or the program ends with exit code $2$. Then main prints the result of solve, which adds P_floor(m) for every $m$ of the range in the 64-bit integer type i64.

P_floor is Sections 3 and 7 in a few lines. For $i=1,\dots,m$ it forms the coordinate xi $=2i/(m+1)$ of $S_m$ and adds $i\ln x_i$ to log_sum with std::log; then it takes std::exp(log_sum), adds the guard $10^{-12}$ of Section 8, rounds down with std::floor and converts the result with static_cast<i64>. It computes in long double, whose format the C++ standard leaves open: many x86 and x86-64 compilers use the 80-bit x87 format, while MSVC uses the same 64-bit format as double. ^[7] On the Apple arm64 machine used for this article long double is the 64-bit format as well, with std::numeric_limits<long double>::digits equal to $53$. A floating-point value that does not fit into the target integer type makes the conversion undefined behaviour, ^[8] which happens for $m\ge23$, beyond the limit of Section 8.

Python and Java contain the same loop with the range $2$ to $15$ built in, without options or checkpoint: Python in the function solve, called at the end of the file, with math.log and math.exp; Java in main, with Math.log, Math.exp and a long total. Both compute in IEEE double precision. Python's int and Java's cast (long) drop the fractional part by rounding towards zero, ^[8] which for positive numbers is the same as rounding down, so all three programs take the integer part in the same way.

Complexity and Verification

For each $m$ the programs evaluate $m$ logarithms and one exponential, so the range $2\le m\le15$ costs $2+3+\dots+15=119$ logarithms and $14$ exponentials, plus ten logarithms and one exponential for the C++ checkpoint. A range up to $M$ costs $O(M^2)$ operations and $O(1)$ memory. The C++ program needs a few milliseconds, Python and Java a few hundredths of a second, mostly start-up time. Unless --skip-checkpoints is given, the C++ program checks the value $\lfloor P_{10}\rfloor=4112$ of the statement before it calls solve. Independent checks for this article, not part of the published solutions:

Exact integer division as in Section 5 gives every $\lfloor P_m\rfloor$ for $2\le m\le15$; GCC and Clang builds of the C++ program, run with --from= and --to= set to the same $m$, agree for every $m$, and the sums printed by all three programs equal the exact sum.
The optimiser SLSQP of SciPy, started at $(1,\dots,1)$ with numerical derivatives and the tolerance ftol=1e-15, maximises $\sum_i i\ln x_i$ under the constraint; for every $m$ from $2$ to $15$ it stops successfully within $4\cdot10^{-8}$ of $S_m$ in every coordinate. Twenty thousand random admissible points per $m$, uniformly distributed, never beat $P_m$, and every transfer of $\delta=0.001$, $0.01$ or $0.1$ from one coordinate of $S_m$ to another lowers the product. These checks only support the result; Section 3 proves it.
For Sections 7 and 8, copies of the three loops, the C++ one built with GCC and Clang at the optimisation levels $0$ and $2$, printed their values of $L_m$ and $P_m$, which were compared with the exact fractions evaluated to $100$ significant digits. In Python, a plain product of the powers has relative errors below $2.2\cdot10^{-15}$ for $m\le15$, against $1.4\cdot10^{-15}$ for the logarithmic loop.
For $m=10$, $100$, $1000$ and $10000$, $m$ times the error of the asymptotic formula of Section 6 is $0.0795$, $0.0829$, $0.08329$ and $0.083329$, approaching $1/12$.

Footnotes and References

The references below supply the background facts used above; the maximiser, the exact fraction, the growth formula for $P_m$ and all numbers are derived or computed in this article.

Project Euler 190 — Maximising a Weighted Product. The official statement: S_m = (x₁, x₂, …, x_m) is the m-tuple of positive real numbers with x₁ + x₂ + ⋯ + x_m = m for which P_m = x₁ · x₂² · ⋯ · x_m^m is maximised; it gives ⌊P₁₀⌋ = 4112 and asks for the sum of ⌊P_m⌋ for m = 2 to 15.
Wikipedia — Lagrange multiplier. For continuously differentiable functions, at a constrained extremum where the gradients of the constraints are linearly independent, the gradient of the objective is a linear combination of them, with the multipliers as coefficients; this is a necessary condition for optimality, not a sufficient one. Here the single constraint has the nonzero gradient (1, …, 1).
Wikipedia — AM–GM inequality. The arithmetic mean of a list of nonnegative real numbers is at least their geometric mean, with equality if and only if all the numbers are equal; in the weighted version with nonnegative weights wₖ, equality holds if and only if all xₖ with wₖ > 0 are equal.
Wikipedia — Hyperfactorial. H(n) = 1¹ · 2² · ⋯ · nⁿ, with the values 1, 1, 4, 108, 27648, 86400000, … from n = 0 on, and H(n) ≈ A · n^((6n² + 6n + 1)/12) · e^(−n²/4), with the Glaisher–Kinkelin constant A ≈ 1.28243 and a correction factor 1 + 1/(720n²) + …. OEIS A002109 lists the same numbers.
Wikipedia — Double-precision floating-point format. The 53-bit significand gives 15 to 17 significant decimal digits; rounding a number of normal size to the nearest double has a relative error of at most 2⁻⁵³, while subnormal numbers trade precision for range; the integers up to 2⁵³ in absolute value are exactly representable.
C++ — Euler190.cpp. This immutable revision, with P_floor, solve, the options and the checkpoint for P₁₀, is the implementation record for this article.
cppreference — Fundamental types. long double is an extended precision type that need not map to an IEEE-754 type; the 80-bit x87 format is used by many x86 and x86-64 implementations, while MSVC implements long double in the same format as double.
cppreference — Implicit conversions: a floating-point value converted to an integer type loses its fractional part, and if the truncated value does not fit, the behaviour is undefined. Python — int: for floating-point numbers, int truncates towards zero. Java Language Specification, §5.1.3: a double is converted to long by rounding towards zero, and a value that is too large becomes the largest long.

Problem 190 source code

C++

#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>

namespace {

using i64 = std::int64_t;

struct Options {
    int from = 2;
    int to = 15;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--from=", options.from) ||
            parse_int_after_prefix(arg, "--to=", options.to)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.from >= 2 && options.to >= options.from;
}

i64 P_floor(const int m) {
    long double log_sum = 0.0L;
    for (int i = 1; i <= m; ++i) {
        const long double xi = 2.0L * static_cast<long double>(i) / static_cast<long double>(m + 1);
        log_sum += static_cast<long double>(i) * std::log(xi);
    }

    const long double value = std::exp(log_sum);
    return static_cast<i64>(std::floor(value + 1e-12L));
}

i64 solve(const int from, const int to) {
    i64 sum = 0;
    for (int m = from; m <= to; ++m) {
        sum += P_floor(m);
    }
    return sum;
}

bool run_checkpoints() {
    if (P_floor(10) != 4112LL) {
        std::cerr << "Checkpoint failed for P(10)" << '\n';
        return false;
    }
    return true;
}

}  // namespace

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

    std::cout << solve(options.from, options.to) << '\n';
    return 0;
}

Python

# Problem 190: Maximising a weighted product
# For each m, maximize Prod(x_i^i) subject to sum(x_i) = m. Optimal: x_i = 2i/(m+1).

import math

def solve():
    total = 0
    for m in range(2, 16):
        log_sum = 0.0
        for i in range(1, m+1):
            xi = 2.0 * i / (m + 1)
            log_sum += i * math.log(xi)
        total += int(math.exp(log_sum) + 1e-12)
    print(total)

solve()

Java

public class Euler190 {
    public static void main(String[] args) {
        long total = 0;
        for (int m = 2; m <= 15; m++) {
            double logSum = 0;
            for (int i = 1; i <= m; i++) {
                double xi = 2.0 * i / (m + 1);
                logSum += i * Math.log(xi);
            }
            total += (long) (Math.exp(logSum) + 1e-12);
        }
        System.out.println(total);
    }
}