Problem 157: Solving the Diophantine Equation 1/a+1/b= p/10^n

Project Euler Problem 157 Solution

EulerSolve provides an optimized solution for Project Euler Problem 157, Solving the Diophantine Equation 1/a+1/b= p/10^n, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $f(n)$ be the number of positive-integer solutions of $$\frac{1}{a}+\frac{1}{b}=\frac{p}{10^n}.$$ The Project Euler task is to compute $\sum_{n=1}^{9} f(n)$. The variable $p$ is part of the solution rather than an input, so the real difficulty is to understand which pairs $(a,b)$ can occur and how many values of $p$ they generate. A direct search is not realistic because $a$ and $b$ are not bounded in advance. The key simplification is to separate the common factor of $a$ and $b$, reduce the problem to a coprime core, and then count all lifts of that core through the divisor function. Mathematical Approach The entire solution turns on one observation: once the common factor of $a$ and $b$ is removed, the remaining coprime pieces can only involve the primes 2 and 5, because the right-hand side has denominator $10^n=2^n5^n$. Strip off the gcd Write $$a=dx,\qquad b=dy,\qquad \gcd(x,y)=1.$$ Substituting into the equation gives $$\frac{1}{dx}+\frac{1}{dy}=\frac{x+y}{dxy}=\frac{p}{10^n},$$ so after clearing denominators we get $$10^n(x+y)=pdxy.$$ For a fixed coprime pair $(x,y)$, all remaining freedom is carried by the common factor $d$ and the numerator $p$....

Detailed mathematical approach

Problem Summary

We count the solutions of

$$\frac1a+\frac1b=\frac{p}{10^n}$$

in positive integers $a$, $b$, $p$ and $n$ with $a\le b$. The numerator $p$ is not given in advance, and the fraction on the right need not be in lowest terms: $\frac12+\frac12=\frac{10}{10}$ is one of the solutions. The exponent $n$ belongs to a solution as well, so $\frac11+\frac11=\frac{20}{10}$ and $\frac11+\frac11=\frac{200}{100}$ are two different solutions. The statement lists the $20$ solutions with $n=1$ and asks for the number of solutions with $1\le n\le9$. ^[1] If $f(n)$ is the number of solutions for one fixed $n$, we need $f(1)+f(2)+\dots+f(9)$, and $f(1)=20$.

The whole solution in one sentence: after dividing $a$ and $b$ by their greatest common divisor $g$, the two coprime numbers that remain are built from the primes $2$ and $5$ alone, so there are only $2n^2+2n+1$ possible pairs, and for each pair the possible values of $g$ are exactly the divisors of one number, so $f(n)$ is a short sum of divisor counts.

Mathematical Approach

Throughout, $N=10^n=2^n5^n$. We split every solution into a shape, the ratio $b/a$ in lowest terms, and a size, the common factor of $a$ and $b$. Only a few shapes are possible, and the sizes that belong to one shape are counted by the number of divisors of a single integer.

1. Why a direct search is too slow

From $a\le b$ we get $p/N=1/a+1/b\le2/a$, so $a\le2N/p\le2N$. Solving for $b$ gives

$$b=\frac{Na}{pa-N},$$

so $pa\gt N$, because $b$ is positive, and $pa\le2N$, because $b\ge a$. The denominator is at least $1$, so $b\le Na$, and $b$ can really be that large: $a=N+1$ and $p=1$ give $b=N(N+1)$, as in $\frac1{11}+\frac1{110}=\frac1{10}$. For $n=9$ this $b$ exceeds $10^{18}$, so a loop over all pairs $(a,b)$ is hopeless. A direct search over $a$ and $p$ is finite: it takes every $a\le2N$ and every integer $p$ with $N/a\lt p\le2N/a$, and keep the candidate when $pa-N$ divides $Na$, so that $b$ is an integer. Writing $\lfloor t\rfloor$ for the largest integer not exceeding $t$, this search tests

$$\sum_{a=1}^{2N}\Bigl(\Bigl\lfloor\frac{2N}{a}\Bigr\rfloor-\Bigl\lfloor\frac{N}{a}\Bigr\rfloor\Bigr)$$

candidates, $39$ for $n=1$ but $22\,263\,991\,779$ for $n=9$. We use this search only as an independent check and count the solutions instead.

2. The twenty solutions of the statement, sorted by shape

Sorting the $20$ solutions for $n=1$ by the ratio $b/a$ shows the structure:

$b=a$: $(1,1)$, $(2,2)$, $(4,4)$, $(5,5)$, $(10,10)$, $(20,20)$;
$b=2a$: $(1,2)$, $(3,6)$, $(5,10)$, $(15,30)$;
$b=5a$: $(1,5)$, $(2,10)$, $(3,15)$, $(4,20)$, $(6,30)$, $(12,60)$;
$b=10a$: $(1,10)$, $(11,110)$;
$2b=5a$: $(2,5)$, $(14,35)$.

Each line consists of the multiples $g\cdot(x,y)$ of its first pair $(x,y)$, and the factors $g$ of a line are exactly the divisors of one number: $1,2,4,5,10,20$ are the divisors of $20$, and the next lines run through the divisors of $15$, $12$, $11$ and $7$. The numerators follow the same numbers, since in every line $pg$ is that number; in the first line $\frac11+\frac11=\frac{20}{10}$ has $p=20$ and $\frac1{20}+\frac1{20}=\frac1{10}$ has $p=1$. So $20=6+4+6+2+2$, and only the five shapes $(1,1)$, $(1,2)$, $(1,5)$, $(1,10)$ and $(2,5)$ occur. Sections 3 to 6 explain why.

3. Removing the common factor

Let $g=\gcd(a,b)$ and write $a=gx$ and $b=gy$. Then $x\le y$, and $x$ and $y$ are coprime, $\gcd(x,y)=1$: they have no prime factor in common. ^[2] We call $(x,y)$ the reduced pair of the solution; $y/x$ is the ratio $b/a$ in lowest terms. Since

$$\frac1{gx}+\frac1{gy}=\frac{x+y}{gxy},$$

multiplying the equation by $gxyN$ turns it into

$$N(x+y)=p\,g\,xy.\qquad(\ast)$$

For the shape $(1,1)$ and $n=1$ this reads $20=pg$, the rule of the first line of Section 2.

4. A sum and a product of coprime numbers are coprime

Claim. If $\gcd(x,y)=1$, then $\gcd(x+y,xy)=1$.

Proof. Suppose a prime $q$ divides both $xy$ and $x+y$. By Euclid's lemma, a prime that divides a product divides one of its factors, ^[3] say $q\mid x$. Then $q$ also divides $(x+y)-x=y$, so $q$ is a common prime factor of $x$ and $y$, which is impossible. Hence $x+y$ and $xy$ have no prime factor in common.

For example, $x=2$ and $y=5$ give $x+y=7$ and $xy=10$. Without coprimality the claim fails: $x=y=2$ give $4$ and $4$.

5. Only the primes 2 and 5 can occur

By $(\ast)$, $xy$ divides $N(x+y)$, and by Section 4 it is coprime to $x+y$. The general form of Euclid's lemma says that if $m$ divides $uv$ and is coprime to $u$, then $m$ divides $v$; ^[3] hence

$$xy\mid N=2^n5^n.$$

By the uniqueness of prime factorisation, ^[4] the divisors of $2^n5^n$ are the numbers $2^i5^j$ with $0\le i,j\le n$. So $x=2^i5^j$ and $y=2^k5^l$ with $i+k\le n$ and $j+l\le n$, and since $x$ and $y$ are coprime, each of the primes $2$ and $5$ divides at most one of them. Two cases remain:

One of the two numbers is $1$. Then the reduced pair is $(1,2^\alpha5^\beta)$ with $0\le\alpha,\beta\le n$, which gives $(n+1)^2$ pairs.
Neither number is $1$. Then each contains at least one of the two primes, and they cannot share one, so one number is $2^\alpha$ and the other is $5^\beta$ with $1\le\alpha,\beta\le n$; the smaller of the two is $x$. This gives $n^2$ pairs.

Altogether there are $(n+1)^2+n^2=2n^2+2n+1$ reduced pairs: $5$ for $n=1$, the five shapes of Section 2, $13$ for $n=2$, $181$ for $n=9$, and $669$ for all $n\le9$ together.

6. Each divisor gives exactly one solution

Fix a pair $(x,y)$ from Section 5 and put

$$K=\frac{N(x+y)}{xy}=N\Bigl(\frac1x+\frac1y\Bigr),$$

an integer because $xy\mid N$. Equation $(\ast)$ says $pg=K$: the common factor $g$ must be a divisor of $K$, and then $p=K/g$. Conversely, for every divisor $g$ of $K$ the numbers $a=gx$, $b=gy$ and $p=K/g$ form a solution with $a\le b$, because $(\ast)$ with $pg=K$ is just the definition of $K$. Moreover $\gcd(gx,gy)=g\gcd(x,y)=g$, ^[2] so Section 3 leads back to the same $g$ and the same pair. Different divisors give different $a$, and different reduced pairs give different ratios $b/a=y/x$, so every solution is counted exactly once:

$$f(n)=\sum_{(x,y)}\tau\Bigl(\frac{N(x+y)}{xy}\Bigr),$$

where the sum runs over the $2n^2+2n+1$ reduced pairs and $\tau(m)$ is the number of positive divisors of $m$. With the exponents of Section 5, $K=2^{n-i-k}\,5^{n-j-l}\,(x+y)$; for $(2,5)$ and $n=1$ this is $K=7$. Since $1/x+1/y\le2$, also $K\le2N$, with equality only for the pair $(1,1)$.

7. Counting divisors

If $m=q_1^{e_1}q_2^{e_2}\cdots q_r^{e_r}$ with distinct primes $q_t$, the divisors of $m$ are the numbers $q_1^{c_1}\cdots q_r^{c_r}$ with $0\le c_t\le e_t$. Each exponent $c_t$ can be chosen in $e_t+1$ ways, independently of the others, and different choices give different divisors by the uniqueness of factorisation. Hence ^[5]

$$\tau(m)=(e_1+1)(e_2+1)\cdots(e_r+1).$$

For example $20=2^2\cdot5$ has $\tau(20)=3\cdot2=6$ divisors, $1,2,4,5,10,20$, the six factors $g$ of the first line of Section 2. The formula needs the factorisation of $K$, and the factor $x+y$ has no predictable one: for $n=2$, the pair $(1,100)$ gives $x+y=101$, a prime, while $(1,25)$ gives $26=2\cdot13$. This is why the program factorises every $K$ separately.

8. Worked example: $n=1$ and $n=2$

For $n=1$ the five reduced pairs give the following values, where the first factor of $K$ is $N/(xy)$ and the second is $x+y$:

$(1,1)$: $K=10\cdot2=20=2^2\cdot5$, so $\tau(K)=6$;
$(1,2)$: $K=5\cdot3=15=3\cdot5$, so $\tau(K)=4$;
$(1,5)$: $K=2\cdot6=12=2^2\cdot3$, so $\tau(K)=6$;
$(1,10)$: $K=1\cdot11=11$, a prime, so $\tau(K)=2$;
$(2,5)$: $K=1\cdot7=7$, a prime, so $\tau(K)=2$.

So $f(1)=6+4+6+2+2=20$, as in the statement, and the divisors of the five values of $K$ are the five lists of factors $g$ in Section 2. For $n=2$ the nine pairs $(1,y)$ with $y=1$, $2$, $4$, $5$, $10$, $20$, $25$, $50$, $100$ give $K=200$, $150$, $125$, $120$, $110$, $105$, $104$, $102$, $101$ with $12$, $12$, $4$, $16$, $8$, $8$, $8$, $8$, $2$ divisors, $78$ in all. The four pairs $(2,5)$, $(2,25)$, $(4,5)$ and $(4,25)$ give $K=70$, $54$, $45$, $29$ with $8$, $8$, $6$, $2$ divisors, $24$ in all. Hence $f(2)=78+24=102$, and the direct search of Section 1 also finds $102$ solutions for $n=2$.

9. A check with a classical count

By Section 6, $p=K/g$ equals $1$ exactly when $g=K$, so every reduced pair gives exactly one solution with $p=1$, and $1/a+1/b=1/N$ has $2n^2+2n+1$ solutions with $a\le b$. A classical argument confirms this. Multiplying $1/a+1/b=1/N$ by $abN$ gives $ab-Na-Nb=0$, and adding $N^2$ to both sides gives

$$(a-N)(b-N)=N^2,$$

where both factors are positive because $1/a$ and $1/b$ are smaller than $1/N$. So the solutions with $a\le b$ correspond to the divisors $d\le N$ of $N^2$, through $a=N+d$ and $b=N+N^2/d$. The divisors of the square $N^2$ pair up as $d$ and $N^2/d$, except for $d=N$, so there are $(\tau(N^2)+1)/2$ of them. ^[6] Here $\tau(10^{2n})=(2n+1)^2$, and indeed $((2n+1)^2+1)/2=2n^2+2n+1$. For $n=1$ these are the five solutions with $p=1$ in the statement, from $\frac1{11}+\frac1{110}$ to $\frac1{20}+\frac1{20}$. In the same way, the choice $g=1$ shows that each reduced pair gives exactly one solution with $\gcd(a,b)=1$.

How the Code Works

The C++ program ^[7] starts in main. parse_arguments accepts --max-n= followed by a number $M$, which replaces the default upper limit $9$, and --skip-checkpoints; any other argument, or $M=0$, ends the program with exit code $1$. Unless the check is skipped, run_checkpoints compares count_for_n(1) with the $20$ solutions of the statement and stops with exit code $2$ on a mismatch. solve then adds count_for_n(n) for $n=1,\dots,M$, which by default is the requested sum over $n\le9$, and the program prints the result.

count_for_n(n) walks through the reduced pairs of Section 5. Its loop variables a2, a5, b2 and b5 are the exponents $i$, $j$, $k$ and $l$: the first two belong to the number left, the last two to right. The loop conditions a2 + b2 <= n and a5 + b5 <= n express $xy\mid N$, and the tests a2 != 0 && b2 != 0 and a5 != 0 && b5 != 0, which skip a tuple, express coprimality. The outer loop runs only for a2 $=0$, so every power of $2$ goes into right and the first test never fires. If b2 $=0$, the bound (b2 ? n : 0) allows only a5 $=0$ as well, which yields the pairs $(1,5^\beta)$ with $0\le\beta\le n$. If b2 $\ge1$, either a5 $=0$ and b5 runs from $0$ to $n$, which yields $(1,2^\alpha5^\beta)$ with $\alpha\ge1$, or a5 $\ge1$, and then the second test leaves only b5 $=0$, which yields the pair of $5^\beta$ and $2^\alpha$. These are the two families of Section 5, each pair exactly once. In the last case left may be the larger number, as in $(5,4)$ for $n=2$; this does not matter, because $K$ is symmetric in $x$ and $y$.

For each tuple that is not skipped, ipow forms the powers by repeated multiplication, scale is $2^{n-i-k}5^{n-j-l}=N/(xy)$, and (left + right) * scale is the number $K$ of Section 6, whose divisor count is added to the total.

divisor_count implements Section 7 by trial division. ^[8] It runs a trial divisor $t=2,3,4,\dots$, the variable p of the code, as long as $t^2$ does not exceed the remaining number, divides $t$ out as often as possible and multiplies count by the exponent plus one; this is why the counter e starts at $1$, and a $t$ that does not divide multiplies count by $1$. A composite $t$ never divides, because its prime factors were removed earlier. When the loop stops, the remaining number has no prime factor below $t$ and is smaller than $t^2$, so it is either $1$ or a prime, and a prime contributes a factor $2$. For $K=20$, $t=2$ divides twice and gives the factor $3$, the loop stops at $t=3$ because $9\gt5$, and the prime $5$ that remains doubles the count to $6$.

The Python version has the same functions ipow, divisor_count, count_for_n and solve, but reads the two options and runs the checkpoint inside main, with the same exit codes; its conversion of the number $M$ also accepts spellings such as +1 that the C++ parser rejects. Its integers have unlimited precision. ^[9] The Java version has the same three helpers, ipow, divisorCount and countForN, with long arithmetic, but its main always adds the counts for $n=1,\dots,9$ and has neither options nor the checkpoint. By Section 6 every $K$ is at most $2\cdot10^9$ for $n\le9$, far inside 64 bits. The 64-bit types of C++ would still hold every $K\le2\cdot10^{18}$ for $n\le18$, but for $n=19$ the value scale $=10^{19}$ of the pair $(1,1)$ already overflows, so values $M\ge19$ are not supported; the parser also reads the digits of $M$ into an int without a range check.

Complexity and Verification

For one $n$ the program evaluates $2n^2+2n+1$ divisor counts, $669$ for all $n\le9$ together. Trial division of a number $K\le2\cdot10^n$ takes at most about $\sqrt{2\cdot10^n}$ steps. The loops also visit $(n+1)+n(n+1)(n+2)/2$ exponent tuples, $505$ for $n=9$, and build each power by repeated multiplication, which adds only $O(n^3)$ operations. Counting every operation on machine integers as one step, the total for all exponents up to $M$ is therefore $O(M^2\,10^{M/2})$: the bounds $10^{n/2}$ grow geometrically, so their sum is dominated by its last term. The functions keep a fixed number of integer variables, so the memory is $O(1)$. In practice the loop bound applies to the remaining number, which shrinks as soon as the powers of $2$ and $5$ and other small factors are divided out: over all $669$ values the loop of divisor_count runs $70\,508$ times, while the worst case for a single value would be about $\sqrt{2\cdot10^9}\approx44\,721$. The longest single run, $8\,163$ steps, is for $K=10^9+5=3\cdot5\cdot66\,666\,667$, from the pair $(1,2\cdot10^8)$ at $n=9$. All three programs finish in a fraction of a second.

Before the main sum, the C++ and Python programs check $f(1)=20$ against the statement. For this article we also checked:

an independent brute force, which loops over all candidates $(a,p)$ of Section 1 and knows nothing about reduced pairs, agrees with the formula of Section 6 for every $n\le9$; at $n=9$ it tests $22\,263\,991\,779$ candidates;
for every $n\le9$, the tuples of count_for_n, with left and right put in increasing order, are exactly the coprime pairs $x\le y$ with $xy\mid10^n$, each once, compared with a separate search over all pairs of divisors of $10^n$;
with $M=1,\dots,6$ the C++ and Python programs print $20$, $122$, $478$, $1436$, $3628$ and $8084$, the cumulative counts of the brute force;
for $n\le4$ the brute force finds $5$, $13$, $25$ and $41$ solutions with $p=1$, and as many with $\gcd(a,b)=1$, as Section 9 predicts;
the three programs print the same total.

Footnotes and References

The references below supply the background facts used above. The reduction to reduced pairs, the divisor-count formula for $f(n)$ and the operation counts are derived explicitly in this article.

Project Euler 157 — Base-10 Diophantine Reciprocal. The official statement considers 1/a + 1/b = p/10ⁿ with positive integers a, b, p, n and a ≤ b, lists the 20 solutions for n = 1 and asks how many solutions the equation has for 1 ≤ n ≤ 9.
Wikipedia — Greatest common divisor. Among the listed properties: if gcd(a, b) = d, then gcd(a/d, b/d) = 1, and gcd(m·a, m·b) = m·gcd(a, b) for every positive integer m. Sections 3 and 6 use these two facts.
Wikipedia — Euclid's lemma. If a prime divides the product of two integers, it divides at least one of them; more generally, if an integer divides a product of two integers and is coprime with one of them, it divides the other. Section 4 uses the first form and Section 5 the second.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a product of primes in only one way, apart from the order of the factors. This gives the form 2ⁱ5ʲ of the divisors of 10ⁿ in Section 5 and the divisor count of Section 7.
Wikipedia — Divisor function. The number of divisors σ₀(n), also written d(n) or τ(n), is the product of the numbers aᵢ + 1 over the prime factorisation n = ∏ pᵢ^(aᵢ), the formula of Section 7.
OEIS A018892 — Number of ways to write 1/n as a sum of exactly 2 unit fractions. The entry gives a(n) = (τ(n²) + 1)/2, the classical count that Section 9 compares with the number of reduced pairs.
C++ — Euler157.cpp. The linked immutable C++ revision contains count_for_n, divisor_count, solve and the checkpoint for n = 1. It is the implementation record for this article.
Wikipedia — Trial division. The simplest factoring programs try successive integers, prime and composite, as divisors, and it suffices to go up to the square root of the number, because a composite number has a factor no larger than its square root.
Python — numeric types. The Python documentation states that integers have unlimited precision, so the Python version cannot overflow.

Problem 157 source code

C++

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

namespace {

using i64 = long long;

struct Options {
    int max_n = 9;
    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, "--max-n=", options.max_n)) {
            continue;
        }

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

    return options.max_n >= 1;
}

i64 ipow(const int base, int exp) {
    i64 out = 1;
    while (exp > 0) {
        out *= base;
        --exp;
    }
    return out;
}

int divisor_count(i64 x) {
    int count = 1;
    for (i64 p = 2; p * p <= x; ++p) {
        int e = 1;
        while (x % p == 0) {
            x /= p;
            ++e;
        }
        count *= e;
    }
    if (x != 1) {
        count *= 2;
    }
    return count;
}

i64 count_for_n(const int n) {
    i64 total = 0;
    for (int a2 = 0; a2 <= 0; ++a2) {
        for (int b2 = 0; a2 + b2 <= n; ++b2) {
            if (a2 != 0 && b2 != 0) {
                continue;
            }
            for (int a5 = 0; a5 <= (b2 ? n : 0); ++a5) {
                for (int b5 = 0; a5 + b5 <= n; ++b5) {
                    if (a5 != 0 && b5 != 0) {
                        continue;
                    }

                    const i64 left = ipow(2, a2) * ipow(5, a5);
                    const i64 right = ipow(2, b2) * ipow(5, b5);
                    const i64 scale = ipow(2, n - a2 - b2) * ipow(5, n - a5 - b5);
                    total += divisor_count((left + right) * scale);
                }
            }
        }
    }
    return total;
}

i64 solve(const int max_n) {
    i64 total = 0;
    for (int n = 1; n <= max_n; ++n) {
        total += count_for_n(n);
    }
    return total;
}

bool run_checkpoints() {
    if (count_for_n(1) != 20) {
        std::cerr << "Checkpoint failed for n=1" << '\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.max_n) << '\n';
    return 0;
}

Python

import sys

def ipow(base, exp):
    out = 1
    while exp > 0:
        out *= base
        exp -= 1
    return out

def divisor_count(x):
    count = 1
    p = 2
    while p * p <= x:
        e = 1
        while x % p == 0:
            x //= p
            e += 1
        count *= e
        p += 1
    if x != 1:
        count *= 2
    return count

def count_for_n(n):
    total = 0
    for a2 in range(0, 1):  # a2 <= 0 means only a2=0
        for b2 in range(0, n - a2 + 1):
            if a2 != 0 and b2 != 0:
                continue
            for a5 in range(0, (n if b2 else 0) + 1):
                for b5 in range(0, n - a5 + 1):
                    if a5 != 0 and b5 != 0:
                        continue

                    left = ipow(2, a2) * ipow(5, a5)
                    right = ipow(2, b2) * ipow(5, b5)
                    scale = ipow(2, n - a2 - b2) * ipow(5, n - a5 - b5)
                    total += divisor_count((left + right) * scale)
    return total

def solve(max_n):
    total = 0
    for n in range(1, max_n + 1):
        total += count_for_n(n)
    return total

def main():
    max_n = 9
    run_checkpoints = True

    args = sys.argv[1:]
    for arg in args:
        if arg == "--skip-checkpoints":
            run_checkpoints = False
        elif arg.startswith("--max-n="):
            try:
                max_n = int(arg[8:])
                if max_n < 1:
                    sys.exit(1)
            except ValueError:
                sys.stderr.write(f"Unknown argument: {arg}\n")
                sys.exit(1)
        else:
            sys.stderr.write(f"Unknown argument: {arg}\n")
            sys.exit(1)

    if run_checkpoints:
        if count_for_n(1) != 20:
            sys.stderr.write("Checkpoint failed for n=1\n")
            sys.exit(2)

    print(solve(max_n))

if __name__ == "__main__":
    main()

Java

public class Euler157 {
    static long ipow(int base, int exp) {
        long out = 1;
        while (exp-- > 0)
            out *= base;
        return out;
    }

    static int divisorCount(long x) {
        int count = 1;
        for (long p = 2; p * p <= x; p++) {
            int e = 1;
            while (x % p == 0) {
                x /= p;
                e++;
            }
            count *= e;
        }
        if (x != 1)
            count *= 2;
        return count;
    }

    static long countForN(int n) {
        long total = 0;
        for (int a2 = 0; a2 <= 0; a2++)
            for (int b2 = 0; a2 + b2 <= n; b2++) {
                if (a2 != 0 && b2 != 0)
                    continue;
                for (int a5 = 0; a5 <= (b2 > 0 ? n : 0); a5++)
                    for (int b5 = 0; a5 + b5 <= n; b5++) {
                        if (a5 != 0 && b5 != 0)
                            continue;
                        long left = ipow(2, a2) * ipow(5, a5);
                        long right = ipow(2, b2) * ipow(5, b5);
                        long scale = ipow(2, n - a2 - b2) * ipow(5, n - a5 - b5);
                        total += divisorCount((left + right) * scale);
                    }
            }
        return total;
    }

    public static void main(String[] args) {
        long total = 0;
        for (int n = 1; n <= 9; n++)
            total += countForN(n);
        System.out.println(total);
    }
}