Problem 141: Investigating Progressive Numbers Which Are Also Square

Project Euler Problem 141 Solution

EulerSolve provides an optimized solution for Project Euler Problem 141, Investigating Progressive Numbers Which Are Also Square, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A progressive number is a positive integer $n$ for which there exists a divisor $d$, a quotient $q$, and a remainder $r$ such that $n=dq+r$, $0<r<d$, and the three integers $d$, $q$, and $r$ are consecutive terms of a geometric progression in some order. Problem 141 asks for the sum of all such $n\le 10^{12}$ that are also perfect squares. The difficulty is not checking one fixed $n$; it is generating only the values that can possibly qualify. The implementations therefore avoid scanning all squares up to $10^{12}$. Instead they parameterize every admissible geometric progression, turn the division identity into two explicit formulas for $n$, and test only those candidates. Mathematical Approach The entire solution rests on a rigid description of three consecutive integer terms in a geometric progression. Once that description is written down, the condition $n=dq+r$ produces only two algebraically distinct families. Parameterizing the geometric progression Let the common ratio be $a/b>1$ in lowest terms, with $\gcd(a,b)=1$ and $a>b$. If the smallest term is $x$, then the three consecutive terms are $$x,\quad x\frac{a}{b},\quad x\frac{a^2}{b^2}.$$ For all three to be integers, $x$ must be a multiple of $b^2$, so we may write $x=tb^2$ for some integer $t\ge 1$....

Detailed mathematical approach

Problem Summary

Dividing a positive integer $n$ by a positive integer $d$ gives a quotient $q$ and a remainder $r$, so that $n=dq+r$ with $0\le r\lt d$. The number $n$ is called progressive if, for some $d$, the three numbers $d$, $q$ and $r$ are consecutive positive integer terms of a geometric sequence, though not necessarily in that order. For example, $58$ divided by $6$ gives the quotient $9$ and the remainder $4$, and $4,6,9$ is a geometric sequence with common ratio $3/2$. Some progressive numbers are perfect squares, such as $9$ and $10404=102^2$, and the progressive perfect squares below $10^5$ add up to $124657$. We must find the sum of all progressive perfect squares below $10^{12}$. ^[1]

The whole solution in one sentence: three positive integers in geometric progression have the form $tb^2$, $tab$, $ta^2$, and the remainder, being smaller than the divisor, is either the middle term, which makes $n$ a product of two consecutive integers and never a square, or the smallest term, which makes $n=t^2a^3b+tb^2$, so we run through all triples $(a,b,t)$ that keep this value within $10^{12}$ and add up the squares.

Mathematical Approach

We look at small cases, describe all geometric progressions of three positive integers and place the remainder in them; this gives two formulas, and only one of them can produce squares.

1. Small progressive numbers

Trying every $d$ from $1$ to $n$ for every $n\lt100$ finds $14$ progressive numbers: $6$, $9$, $12$, $20$, $28$, $30$, $34$, $42$, $56$, $58$, $65$, $72$, $75$ and $90$. They come in two kinds:

$6=4\cdot1+2$, $12=9\cdot1+3$, $20=16\cdot1+4$, …: the remainder is the middle term of $1,2,4$, $1,3,9$, $1,4,16$, …, and these eight numbers, $6$ to $90$, are the products $2\cdot3$, $3\cdot4$, …, $9\cdot10$;
$9=2\cdot4+1$, $28=3\cdot9+1$, $34=4\cdot8+2$, $58=6\cdot9+4$, $65=4\cdot16+1$ and $75=6\cdot12+3$: the remainder is the smallest term, of $1,2,4$, $1,3,9$, $2,4,8$, $4,6,9$, $1,4,16$ and $3,6,12$.

A number can have several such divisions: $9=4\cdot2+1$ as well, with divisor and quotient exchanged, and $42$ has four: $42=36\cdot1+6$${}=18\cdot2+6$${}=12\cdot3+6$${}=9\cdot4+6$. The only square among the fourteen numbers is $9$. The progressive numbers form the sequence A335065 of the OEIS, which splits them into exactly these two kinds. ^[2]

2. Three positive integers in geometric progression

Three positive numbers $x\lt y\lt z$ are consecutive terms of a geometric sequence exactly when $y/x=z/y$, that is, $y^2=xz$: the square of the middle term is the product of its two neighbours. ^[3] Our three terms are all different, because two equal terms would force the ratio $1$ and make all three equal, while $r\lt d$; so we list them in increasing order, with ratio $y/x\gt1$.

Claim. Positive integers $x\lt y\lt z$ satisfy $y^2=xz$ exactly when

$$x=tb^2,\qquad y=tab,\qquad z=ta^2$$

for integers $a\gt b\ge1$ with $\gcd(a,b)=1$ and $t\ge1$, and these $a$, $b$, $t$ are unique.

Proof. Write the ratio $y/x$ in lowest terms as $a/b$, so that $\gcd(a,b)=1$ and $a\gt b$. Then $z=y\cdot a/b=x\,a^2/b^2$, so $b^2z=a^2x$. A prime $p$ of $b^2$ divides $b$, not $a$, so by unique prime factorisation it occurs in $a^2x$ exactly as often as in $x$; ^[4] hence its power in $b^2$, which divides $a^2x$, divides $x$, and so $b^2$ divides $x$. Writing $x=tb^2$ gives $y=xa/b=tab$ and $z=ta^2$. Conversely, $(tab)^2=tb^2\cdot ta^2$. The terms determine the reduced fraction $a/b$, and then $t=x/b^2$. For example, $4,6,9$ has the ratio $3/2$ and $(a,b,t)=(3,2,1)$, $2,4,8$ has $(2,1,2)$, and $3,6,12$ has $(2,1,3)$.

3. Where the remainder can stand

The remainder is smaller than the divisor, so it is not the largest of the three terms. This leaves two cases:

the remainder is the middle term; then the divisor, being larger, is the largest term and the quotient the smallest;
the remainder is the smallest term; then the divisor and the quotient are the other two terms, in either order.

Conversely, every such arrangement is a genuine division: if $n=dq+r$ with $0\lt r\lt d$, then $q$ and $r$ are the quotient and the remainder of $n$ divided by $d$, because division with remainder has only one result. ^[5] So the two cases together describe all progressive numbers.

4. Remainder in the middle: never a square

If $q\lt r\lt d$, the square of the middle term is the product of the outer terms, $qd=r^2$, and therefore

$$n=dq+r=r^2+r=r(r+1).$$

For $r\ge1$ we have $r^2\lt r(r+1)\lt(r+1)^2$, so $n$ lies strictly between two consecutive squares and is never a perfect square. Conversely, every product $r(r+1)$ with $r\ge2$ is progressive: divided by $d=r^2$ it leaves the quotient $1$ and the remainder $r\lt r^2$, and $1,r,r^2$ is geometric. These are the pronic numbers from $6$ on, the first kind of Section 1. ^[6]

5. Remainder smallest: the formula

If the remainder is the smallest term, Section 2 gives $r=tb^2$, and the divisor and the quotient are $tab$ and $ta^2$ in some order. Their product is $t^2a^3b$ either way, so

$$n=t^2a^3b+tb^2=tb\,(ta^3+b)$$

with $a\gt b\ge1$, $\gcd(a,b)=1$ and $t\ge1$; the condition $r\lt d$ holds automatically, because $r$ is the smallest term. The example of the statement is $(a,b,t)=(3,2,1)$: the terms are $4,6,9$ and $n=27\cdot2+4=58$. Together with Sections 3 and 4 this proves the key fact: a perfect square is progressive exactly when it equals $t^2a^3b+tb^2$ for such a triple $(a,b,t)$. With the middle term $m=tab$, the largest term is $m^2/r$, so $n=m^3/r+r$; conversely, every divisor $r\lt m$ of $m^2$ gives the progression $r,m,m^2/r$. These numbers form the OEIS sequence A127629. ^[2]

6. Worked example: the squares below $10^5$

Exactly $792$ triples $(a,b,t)$ give a value $t^2a^3b+tb^2$ below $10^5$. They produce only $785$ different numbers, because seven numbers arise twice; for instance $4002=20\cdot200+2$ comes from $(a,b,t)=(10,1,2)$, and $4002=33\cdot121+9$ from $(11,3,1)$. Only four of the values are squares:

$(2,1,1)$: terms $1,2,4$ and $2\cdot4+1=9=3^2$;
$(2,1,36)$: terms $36,72,144$ and $72\cdot144+36=10404=102^2$;
$(3,1,25)$: terms $25,75,225$ and $75\cdot225+25=16900=130^2$;
$(23,2,2)$: terms $8,92,1058$ and $92\cdot1058+8=97344=312^2$.

Their sum $9+10404+16900+97344$ is $124657$, the value given in the statement. In the last one, $97344$ divided by $92$ leaves the quotient $1058$ and the remainder $8$, and the ratio is $23/2$.

7. The case $b=1$ and Pell's equation

For $b=1$ the formula reads $n=t\,(ta^3+1)$. A common divisor of $t$ and $ta^3+1$ also divides $1$, so the two factors are coprime. Every prime then occurs in only one of them, and in a square every prime occurs an even number of times, so the product is a square only if both factors are squares: $t=u^2$ and $ta^3+1=v^2$, that is,

$$v^2-a^3u^2=1.$$

This is Pell's equation, which has infinitely many solutions whenever $a^3$ is not a perfect square, that is, whenever $a$ is not a square. ^[7] For $a=2$ the smallest positive solution is $(v,u)=(3,1)$, and the step $(v,u)\mapsto(3v+8u,\,v+3u)$ turns a solution into a larger one, because expanding gives $(3v+8u)^2-8(v+3u)^2=v^2-8u^2$; it yields $(17,6)$, $(99,35)$, … The squares $n=(uv)^2$ are $3^2$, $102^2$ and $3465^2$, and for $a=3$ the smallest solution $(26,5)$ gives $130^2$. So the two squares of the statement belong to one Pell family, and there are infinitely many progressive squares. For $b\gt1$ the factors can share divisors; the programs use none of this and test every value of Section 5.

8. How far the search must go

Let $N=10^{12}$. Since $a\gt b$, we have $t^2a^3b\gt t^2a^2b^2$, so every value of Section 5 that does not exceed $N$ has $tab\lt\sqrt N=10^6$. The programs use this crude bound in the slightly weaker form $t^2a^2b^2\le N$, which never skips a valid triple. The exact condition $t^2a^3b+tb^2\le N$ is much stronger: since the left side is at least $a^3+1$, it forces $a^3\lt N$, so $a\le9999$, and only $7\,302\,731$ triples satisfy it.

9. Repeated values and the boundary

The condition $\gcd(a,b)=1$ makes different triples describe different progressions; without it, $(2a,2b,t)$ would repeat the progression of $(a,b,4t)$, so the test removes repetitions, not solutions. Different progressions can still give the same number, as $4002$ shows, so the squares are collected in a set before they are added, although below $10^{12}$ no square arises twice. The statement asks for squares below $10^{12}$, while the programs keep values up to $10^{12}$ inclusive; this makes no difference, because $10^{12}=(10^6)^2$ is not progressive, as the complete search confirms.

How the Code Works

The C++ program ^[8] starts in main. parse_arguments accepts --limit=N, by default $10^{12}$, and --skip-checkpoints; any other argument, or a limit of $0$, ends the program with exit code $1$. Unless the checkpoint is skipped, run_checkpoints requires solve(100000) to return $124657$, the sum of Section 6, and otherwise the program stops with exit code $2$. Then main prints solve(limit).

solve collects the progressive squares in a std::set. Its outer loop runs over $a$ while $a^2\le N$, its middle loop over $b\lt a$ while $a^2b^2\le N$, skipping the pairs for which gcd_int is not $1$, and its inner loop over the scale factor k, the $t$ of Section 2, while $k^2a^2b^2\le N$; with $N$ the limit, these are the crude bounds of Section 8. For each triple it forms candidate1, equal to $k^2a^2b^2+kab$, the pronic value of Section 4, and candidate2, equal to $k^2a^3b+kb^2$, the value of Section 5. A candidate that does not exceed the limit and passes is_square goes into the set, and at the end the set is summed. By Section 4, candidate1 never passes: at $N=10^{12}$ its $36\,429\,311$ square tests are wasted, although they do not change the result. Likewise the outer loop runs on to $a=10^6$, although no $a\gt9999$ contributes.

is_square truncates the long double square root to an integer root and accepts the value if root or root + 1 squared equals it. Every tested value is at most $10^{12}\lt2^{53}$, so it is converted exactly and its root is far too accurate to miss a whole number; the second test covers a root just below a whole number. The Python version runs the same loops with while statements, math.gcd and the exact integer square root math.isqrt, and keeps a set. The Java version keeps a TreeSet<Long> and tests squares like the C++ program, with the correctly rounded Math.sqrt. ^[9] Neither has options or a checkpoint; both always use $N=10^{12}$.

Integer sizes: the largest product formed at $N=10^{12}$ is $k^2a^3b=10^{18}$, at $a=10^6$ and $b=k=1$, below $2^{63}\approx9.22\cdot10^{18}$; so the unsigned 64-bit u64 of C++ and the signed long of Java hold every value exactly, and Python integers have no fixed size. The option --limit is not safe for every value: from limits of about $6.98\cdot10^{12}$ on, the outer loop reaches $a=2\,642\,246$, whose cube exceeds $2^{64}$, and unsigned arithmetic silently wraps around modulo $2^{64}$. ^[10] At the limit $2\cdot10^{13}$, for instance, $a=2^{20}$, $b=1$ and $k=4$ give $2^{64}+4$, which wraps around to $4$, so the false squares $4$ and $1$ enter the set. The parser does not check for overflow either: a limit given as $2^{64}+1$ is read as $1$.

Complexity and Verification

Write $M=\sqrt N$, which is $10^6$ here. The inner loop runs $\lfloor M/(ab)\rfloor$ times for each coprime pair $(a,b)$ with $b\lt a$, so its total is at most

$$\sum_{a=1}^{M}\sum_{b=1}^{M}\frac{M}{ab}=M\,H_M^2,$$

where $H_M=1+\tfrac12+\dots+\tfrac1M$ is the $M$th harmonic number; $H_M\le1+\ln M$, because for $j\ge2$ the term $1/j$ is at most the integral of $1/x$ from $j-1$ to $j$. The middle loop visits at most $M\,H_M$ pairs, each with a $\gcd$ of $O(\log N)$ steps. Counting arithmetic and square tests as single steps, the loops take $O(\sqrt N\log^2N)$ time, plus $O(\log S)$ for each insertion into the ordered sets of C++ and Java, where $S$ is the size of the set. ^[11] For $N=10^{12}$ the outer loop makes $10^6$ passes, the middle loop visits $6\,984\,517$ pairs, of which $4\,592\,842$ are coprime, and the inner loop runs $36\,429\,395$ times, well below the bound $M\,H_M^2\approx2.07\cdot10^8$, with a few multiplications and at most two square tests per pass. The set ends with only $S=23$ values, so the insertions are negligible and the memory is $O(S)$ plus a few loop variables. On a desktop computer the C++ and Java programs finish in well under a second and the Python program in about $14$ seconds. With the exact bound of Section 8 and without candidate1, the loops would visit only the $7\,302\,731$ useful triples, a fifth of the inner passes.

The C++ program checks the value $124657$ of the statement before the main run. The following checks were made for this article and are not part of the published solutions:

a brute force from the definition, trying every $d$ for every $n\lt10^5$, finds $1090$ progressive numbers: the $314$ pronic numbers $r(r+1)$ with $r\ge2$ and the $785$ values of Section 6, nine of which, such as $110=10\cdot11=6\cdot18+2$, are of both kinds;
the same test of every $d$, applied to every square below $10^7$, finds the progressive squares $9$, $10404$, $16900$, $97344$, $576081$, $6230016$ and $7322436$, exactly those given by Section 5;
at $N=10^{12}$, a separate enumeration of $n=m^3/r+r$ from Section 5, over every middle term $m$ and every divisor $r\lt m$ of $m^2$, finds the same set of squares as the loops of the programs, none of them twice, and not $10^{12}$ itself;
none of the $36\,429\,311$ values of candidate1 up to $10^{12}$ is a square, as Section 4 predicts;
all three programs print the same sum as this enumeration, and a 128-bit version confirms the wrap-around at the limit $2\cdot10^{13}$.

Footnotes and References

The references below supply the background facts used above. The description of the progressions, the two formulas, the reduction of the case $b=1$ to Pell's equation and the operation counts are derived in this article.

Project Euler 141 — Square Progressive Numbers. The official statement: n divided by d has quotient q and remainder r, and d, q, r are consecutive positive integer terms of a geometric sequence in some order, as in 58 = 6·9 + 4 with 4, 6, 9. It names 9 and 10404 = 102² as progressive squares, gives 124657 as the sum of the progressive perfect squares below one hundred thousand and asks for the sum below one trillion.
OEIS A335065 — d, q, r in geometric progression. The entry, inspired by Problem 141, splits the terms into the cases r < q < d and r < d < q, which give the same numbers, the sequence A127629, and q < r < d, which gives the oblong numbers from 6 on.
Wikipedia — Geometric progression. Each term after the first is the previous term multiplied by a fixed number, the common ratio, and geometric sequences satisfy aₙ = aₙ₋₁²/aₙ₋₂: the square of a term is the product of its two neighbours.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a product of primes in exactly one way, up to the order of the factors.
Wikipedia — Euclidean division. For integers a and b with b ≠ 0 there are unique integers q and r with a = bq + r and 0 ≤ r < |b|; q is called the quotient and r the remainder.
Wikipedia — Pronic number. A pronic or oblong number is a product n(n + 1) of two consecutive integers; the article gives the chain n² ≤ n(n + 1) < (n + 1)² < (n + 1)(n + 2).
Wikipedia — Pell's equation. For a positive integer n that is not a perfect square, x² − ny² = 1 has infinitely many integer solutions; the fundamental solution (x₁, y₁) is the one in positive integers with the least x, and the powers xₖ + yₖ√n = (x₁ + y₁√n)ᵏ give further solutions. The table lists (3, 1) for n = 8 and (26, 5) for n = 27.
C++ — Euler141.cpp. The linked immutable revision, with solve, is_square and run_checkpoints, is the implementation record for this article.
Python — math.isqrt returns the floor of the exact square root of a nonnegative integer, and Java — Math.sqrt returns the correctly rounded square root, the double value closest to the true mathematical square root.
cppreference — Arithmetic operators. Under "Overflows": unsigned integer arithmetic is always performed modulo 2ⁿ, where n is the number of bits of the type.
cppreference — std::set: search, removal and insertion have logarithmic complexity. Java — TreeSet: guaranteed log(n) time for add, remove and contains.

Problem 141 source code

C++

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

namespace {

using u64 = std::uint64_t;

struct Options {
    u64 limit = 1000000000000ULL;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<u64>(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_u64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }

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

int gcd_int(int a, int b) {
    while (b != 0) {
        const int t = a % b;
        a = b;
        b = t;
    }
    return a;
}

bool is_square(const u64 value) {
    const u64 root = static_cast<u64>(std::sqrt(static_cast<long double>(value)));
    return root * root == value || (root + 1) * (root + 1) == value;
}

u64 solve(const u64 limit) {
    std::set<u64> progressive_squares;

    for (u64 a = 1; a * a <= limit; ++a) {
        for (u64 b = 1; b < a && a * a * b * b <= limit; ++b) {
            if (gcd_int(static_cast<int>(a), static_cast<int>(b)) != 1) {
                continue;
            }

            for (u64 k = 1; k * k * a * a * b * b <= limit; ++k) {
                const u64 candidate1 = k * a * a * k * b * b + k * a * b;
                if (candidate1 <= limit && is_square(candidate1)) {
                    progressive_squares.insert(candidate1);
                }

                const u64 candidate2 = k * a * a * k * a * b + k * b * b;
                if (candidate2 <= limit && is_square(candidate2)) {
                    progressive_squares.insert(candidate2);
                }
            }
        }
    }

    u64 sum = 0;
    for (u64 value : progressive_squares) {
        sum += value;
    }
    return sum;
}

bool run_checkpoints() {
    if (solve(100000ULL) != 124657ULL) {
        std::cerr << "Checkpoint failed for limit 100000" << '\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.limit) << '\n';
    return 0;
}

Python

import math

def is_square(value):
    root = math.isqrt(value)
    return root * root == value

def solve():
    limit = 10**12
    progressive_squares = set()
    a = 1
    while a * a <= limit:
        b = 1
        while b < a and a * a * b * b <= limit:
            if math.gcd(a, b) != 1:
                b += 1
                continue
            k = 1
            while k * k * a * a * b * b <= limit:
                candidate1 = k * a * a * k * b * b + k * a * b
                if candidate1 <= limit and is_square(candidate1):
                    progressive_squares.add(candidate1)
                candidate2 = k * a * a * k * a * b + k * b * b
                if candidate2 <= limit and is_square(candidate2):
                    progressive_squares.add(candidate2)
                k += 1
            b += 1
        a += 1
    return str(sum(progressive_squares))

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

Java

import java.util.*;

public class Euler141 {
    static boolean isSquare(long v) {
        long r = (long) Math.sqrt((double) v);
        return r * r == v || (r + 1) * (r + 1) == v;
    }

    static int gcd(int a, int b) {
        while (b != 0) {
            int t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    public static void main(String[] args) {
        long limit = 1000000000000L;
        Set<Long> set = new TreeSet<>();
        for (long a = 1; a * a <= limit; a++)
            for (long b = 1; b < a && a * a * b * b <= limit; b++) {
                if (gcd((int) a, (int) b) != 1)
                    continue;
                for (long k = 1; k * k * a * a * b * b <= limit; k++) {
                    long c1 = k * a * a * k * b * b + k * a * b;
                    if (c1 <= limit && isSquare(c1))
                        set.add(c1);
                    long c2 = k * a * a * k * a * b + k * b * b;
                    if (c2 <= limit && isSquare(c2))
                        set.add(c2);
                }
            }
        long sum = 0;
        for (long v : set)
            sum += v;
        System.out.println(sum);
    }
}