Problem 932: $2025$

Project Euler Problem 932 Solution

EulerSolve provides an optimized solution for Project Euler Problem 932, $2025$, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We seek square numbers $n$ below $10^{16}$ whose decimal expansion can be split into a nonempty left part and a nonempty right part so that the square root is the sum of those two parts. Writing $n=s^2$, the split has the form $$n=a10^d+b,\qquad s=a+b,$$ where $b$ has exactly $d$ digits. The title example is $2025=20\mid 25$, because $20+25=45=\sqrt{2025}$. The task is to add all such squares, not all successful splits, so duplicates must be counted only once. This is closely related to base-10 Kaprekar numbers, but the efficient way to solve the Project Euler problem is to search over the root $s$ and describe the allowed values of $s$ by modular arithmetic. Mathematical Approach Fix a split length $d$ and write $$M_d=10^d-1.$$ For that fixed $d$, the whole problem is to characterize the roots $s$ for which the split exists. Eliminating the two decimal pieces From $n=s^2=a10^d+b$ and $s=a+b$, substitute $b=s-a$: $$s^2=a10^d+(s-a)=s+a(10^d-1).$$ Therefore $$a=\frac{s(s-1)}{M_d},\qquad b=s-a.$$ This already gives the key divisibility condition $$M_d\mid s(s-1).$$ Once $a$ is an integer and the digit conditions hold, the reconstruction is automatic, because $$a10^d+b=a(10^d-1)+s=s(s-1)+s=s^2.$$ So the search is not really over $(a,b)$; it is over roots $s$ satisfying one congruence and a few digit inequalities....

Detailed mathematical approach

Problem Summary

We seek square numbers $n$ below $10^{16}$ whose decimal expansion can be split into a nonempty left part and a nonempty right part so that the square root is the sum of those two parts. Writing $n=s^2$, the split has the form

$$n=a10^d+b,\qquad s=a+b,$$

where $b$ has exactly $d$ digits. The title example is $2025=20\mid 25$, because $20+25=45=\sqrt{2025}$. The task is to add all such squares, not all successful splits, so duplicates must be counted only once.

This is closely related to base-10 Kaprekar numbers, but the efficient way to solve the Project Euler problem is to search over the root $s$ and describe the allowed values of $s$ by modular arithmetic.

Mathematical Approach

Fix a split length $d$ and write

$$M_d=10^d-1.$$

For that fixed $d$, the whole problem is to characterize the roots $s$ for which the split exists.

Eliminating the two decimal pieces

From $n=s^2=a10^d+b$ and $s=a+b$, substitute $b=s-a$:

$$s^2=a10^d+(s-a)=s+a(10^d-1).$$

Therefore

$$a=\frac{s(s-1)}{M_d},\qquad b=s-a.$$

This already gives the key divisibility condition

$$M_d\mid s(s-1).$$

Once $a$ is an integer and the digit conditions hold, the reconstruction is automatic, because

$$a10^d+b=a(10^d-1)+s=s(s-1)+s=s^2.$$

So the search is not really over $(a,b)$; it is over roots $s$ satisfying one congruence and a few digit inequalities.

Why every prime power allows only residue $0$ or $1$

Factor $M_d$ into pairwise coprime prime powers:

$$M_d=\prod_{i=1}^{t} p_i^{e_i}.$$

Consecutive integers are coprime, so $\gcd(s,s-1)=1$. If $p_i^{e_i}$ divides the product $s(s-1)$, then the full prime power must divide exactly one of the two consecutive factors. Hence every valid root satisfies

$$s\equiv 0,1 \pmod{p_i^{e_i}},\qquad 1\le i\le t.$$

That is the crucial structural simplification. Instead of examining all roots below $10^8$, we only have to choose, for every distinct prime-power factor of $10^d-1$, whether $s$ is congruent to $0$ or to $1$.

Combining the local choices with the Chinese remainder theorem

Each prime power contributes two local options, so a fixed $d$ leads to a collection of global residue classes. The Chinese remainder theorem merges those local choices into unique congruences

$$s\equiv r \pmod{M_d}.$$

Every admissible root for that split length lies in one arithmetic progression

$$s=r+kM_d,\qquad k\ge 0.$$

The implementations therefore build all CRT residue classes first and only then test the corresponding candidate roots.

The digit bounds that finish the reduction

The modular condition is necessary, but not sufficient. We still need the split to be a genuine decimal split with both parts nonzero and with no leading zeros in the right block.

Because $b=s-a\gt 0$, we need $a\lt s$. Using the formula for $a$,

$$\frac{s(s-1)}{M_d}\lt s,$$

which simplifies to

$$s\lt 10^d.$$

So every valid root must satisfy

$$2\le s\le \min\!\bigl(\lfloor\sqrt{10^{16}-1}\rfloor,\;10^d-1\bigr).$$

There is also an exact range for the right part:

$$1\le b\lt 10 \quad(d=1),\qquad 10^{d-1}\le b\lt 10^d \quad(d\ge 2).$$

The second condition is what rules out right parts such as $01$ or $0017$: they are numerically positive, but they do not occupy exactly $d$ digits.

Worked example: the title number $2025$

Take $d=2$. Then

$$M_2=10^2-1=99=9\cdot 11.$$

For the prime powers $9$ and $11$, the only local residues are $0$ and $1$. CRT therefore produces four global residue classes modulo $99$:

$$s\equiv 0,\ 1,\ 45,\ 55 \pmod{99}.$$

Since $s\lt 100$, only the representatives $45$, $55$, and $99$ can matter.

For $s=45$,

$$a=\frac{45\cdot 44}{99}=20,\qquad b=45-20=25,$$

so $45^2=2025=20\mid 25$.

For $s=55$,

$$a=\frac{55\cdot 54}{99}=30,\qquad b=55-30=25,$$

so $55^2=3025=30\mid 25$.

For $s=99$,

$$a=\frac{99\cdot 98}{99}=98,\qquad b=1,$$

but $b=1$ is not a two-digit number, so the split $98\mid 01$ is invalid. This example shows exactly why the algorithm needs both pieces: CRT generates the modular candidates, and the digit checks remove the false positives.

Why the final sum must be de-duplicated

The Project Euler question asks for the sum of the numbers $n$, not the sum over all successful pairs $(n,d)$. A square can satisfy the defining property for more than one split length, so the accepted values $n=s^2$ must be stored as distinct integers before the final addition is performed.

How the Code Works

Precomputation

The C++, Python, and Java implementations precompute powers of $10$, the root bound $\lfloor\sqrt{10^{16}-1}\rfloor$, and a prime table large enough to factor every repunit-like modulus $10^d-1$ with $1\le d\le 15$.

Per-split residue construction

For each split length $d$, the implementation factors $10^d-1$ into prime powers. It starts from the trivial congruence modulo $1$ and repeatedly merges the two local possibilities $s\equiv0$ and $s\equiv1$ modulo each prime power. After all merges, the result is the complete CRT list of admissible residues $r$ modulo $10^d-1$.

Candidate testing and unique accumulation

Each residue class produces a candidate root $s$ in the allowed range. From that root, the implementation computes $a=\frac{s(s-1)}{10^d-1}$, then $b=s-a$, checks positivity, checks the exact digit range of $b$, reconstructs $n=s^2$, and inserts $n$ into a hash set. The C++ implementation also includes small direct-search self-checks for shorter digit limits before producing the final $10^{16}$ answer.

Complexity Analysis

Let $\omega(M_d)$ be the number of distinct prime divisors of $M_d=10^d-1$. For a fixed $d$, the CRT stage creates exactly $2^{\omega(M_d)}$ residue classes, because each prime power contributes the binary choice $0$ or $1$.

For this particular problem size, the bound $s\lt 10^d$ means that each residue class contributes at most one relevant root candidate. As a result, after factorization the remaining work is essentially one constant-time arithmetic check per CRT residue. Memory usage is small: a prime list, a temporary residue list, and the set of accepted squares. This is far smaller than a naive search over all roots up to $10^8$ and all possible split positions.

Footnotes and References

Problem page: Project Euler 932
Chinese remainder theorem: Wikipedia - Chinese remainder theorem
Kaprekar number: Wikipedia - Kaprekar number
Modular arithmetic: Wikipedia - Modular arithmetic
Repunit: Wikipedia - Repunit

Problem 932 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <unordered_set>
#include <utility>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;
using i128 = __int128_t;

u64 pow10_int(int e) {
    u64 v = 1;
    for (int i = 0; i < e; ++i) {
        v *= 10;
    }
    return v;
}

std::vector<int> sieve_primes(int limit) {
    std::vector<bool> composite(limit + 1, false);
    std::vector<int> primes;
    primes.reserve(limit / 10);

    for (int i = 2; i <= limit; ++i) {
        if (!composite[i]) {
            primes.push_back(i);
            if (static_cast<u64>(i) * static_cast<u64>(i) <= static_cast<u64>(limit)) {
                for (int j = i * i; j <= limit; j += i) {
                    composite[j] = true;
                }
            }
        }
    }

    return primes;
}

std::vector<u64> factor_prime_powers(u64 x, const std::vector<int>& primes) {
    std::vector<u64> parts;

    for (int p : primes) {
        const u64 pp = static_cast<u64>(p);
        if (static_cast<u128>(pp) * pp > x) {
            break;
        }
        if (x % pp != 0) {
            continue;
        }
        u64 pk = 1;
        while (x % pp == 0) {
            x /= pp;
            pk *= pp;
        }
        parts.push_back(pk);
    }

    if (x > 1) {
        parts.push_back(x);
    }

    return parts;
}

i128 egcd(i128 a, i128 b, i128& x, i128& y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    i128 x1 = 0;
    i128 y1 = 0;
    i128 g = egcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

u64 mod_inv(u64 a, u64 mod) {
    i128 x = 0;
    i128 y = 0;
    i128 g = egcd(static_cast<i128>(a), static_cast<i128>(mod), x, y);
    assert(g == 1);
    i128 r = x % static_cast<i128>(mod);
    if (r < 0) {
        r += static_cast<i128>(mod);
    }
    return static_cast<u64>(r);
}

std::pair<u64, u64> crt_merge(u64 r1, u64 m1, u64 r2, u64 m2) {
    if (m1 == 1) {
        return {r2 % m2, m2};
    }
    const u64 inv = mod_inv(m1 % m2, m2);

    i128 diff = static_cast<i128>(r2) - static_cast<i128>(r1);
    diff %= static_cast<i128>(m2);
    if (diff < 0) {
        diff += static_cast<i128>(m2);
    }

    const u64 t = static_cast<u64>((diff * static_cast<i128>(inv)) % static_cast<i128>(m2));
    const u128 mod = static_cast<u128>(m1) * static_cast<u128>(m2);
    const u128 r = static_cast<u128>(r1) + static_cast<u128>(m1) * static_cast<u128>(t);
    return {static_cast<u64>(r % mod), static_cast<u64>(mod)};
}

u64 digit_count(u64 x) {
    u64 d = 0;
    do {
        ++d;
        x /= 10;
    } while (x > 0);
    return d;
}

bool is_2025_number(u64 n, const std::vector<u64>& pow10) {
    const int digits = static_cast<int>(digit_count(n));

    for (int d = 1; d < digits; ++d) {
        const u64 p10 = pow10[d];
        const u64 a = n / p10;
        const u64 b = n % p10;
        if (a == 0 || b == 0) {
            continue;
        }
        if (d > 1 && b < pow10[d - 1]) {
            continue;
        }
        const u64 s = a + b;
        if (static_cast<u128>(s) * s == n) {
            return true;
        }
    }

    return false;
}

u64 brute_sum_digits(int nd) {
    const u64 max_n = pow10_int(nd) - 1;
    u64 root = static_cast<u64>(std::sqrt(static_cast<long double>(max_n)));
    while (static_cast<u128>(root + 1) * (root + 1) <= max_n) {
        ++root;
    }
    while (static_cast<u128>(root) * root > max_n) {
        --root;
    }

    std::vector<u64> pow10(nd + 1, 1);
    for (int i = 1; i <= nd; ++i) {
        pow10[i] = pow10[i - 1] * 10;
    }

    u128 total = 0;
    for (u64 s = 2; s <= root; ++s) {
        const u64 n = static_cast<u64>(static_cast<u128>(s) * s);
        if (is_2025_number(n, pow10)) {
            total += n;
        }
    }

    return static_cast<u64>(total);
}

u64 solve_digits(int nd) {
    const u64 max_n = pow10_int(nd) - 1;
    u64 s_limit = static_cast<u64>(std::sqrt(static_cast<long double>(max_n)));
    while (static_cast<u128>(s_limit + 1) * (s_limit + 1) <= max_n) {
        ++s_limit;
    }
    while (static_cast<u128>(s_limit) * s_limit > max_n) {
        --s_limit;
    }

    std::vector<u64> pow10(nd + 1, 1);
    for (int i = 1; i <= nd; ++i) {
        pow10[i] = pow10[i - 1] * 10;
    }

    const int sieve_limit = static_cast<int>(std::sqrt(static_cast<long double>(pow10[nd - 1] - 1))) + 10;
    const std::vector<int> primes = sieve_primes(sieve_limit);

    std::unordered_set<u64> seen;
    seen.reserve(256);

    for (int d = 1; d < nd; ++d) {
        const u64 mod = pow10[d] - 1;
        const std::vector<u64> prime_powers = factor_prime_powers(mod, primes);

        std::vector<std::pair<u64, u64>> residues = {{0, 1}};
        for (u64 pk : prime_powers) {
            std::vector<std::pair<u64, u64>> next;
            next.reserve(residues.size() * 2);
            for (const auto& [r, m] : residues) {
                next.push_back(crt_merge(r, m, 0, pk));
                next.push_back(crt_merge(r, m, 1 % pk, pk));
            }
            residues.swap(next);
        }

        const u64 s_max = std::min(s_limit, pow10[d] - 1);
        const u64 b_low = (d == 1) ? 1 : pow10[d - 1];
        const u64 p10d = pow10[d];

        for (const auto& [r0, step] : residues) {
            u64 s = r0;
            if (s <= 1) {
                const u64 need = 2 - s;
                const u64 k = (need + step - 1) / step;
                if (k > (s_max - s) / step) {
                    continue;
                }
                s += k * step;
            }

            for (; s <= s_max; s += step) {
                const u128 ss = static_cast<u128>(s) * static_cast<u128>(s - 1);
                if (ss % mod != 0) {
                    continue;
                }

                const u64 a = static_cast<u64>(ss / mod);
                if (a == 0 || a >= s) {
                    continue;
                }
                const u64 b = s - a;

                if (b < b_low || b >= p10d) {
                    continue;
                }

                const u64 n = static_cast<u64>(static_cast<u128>(s) * s);
                if (n > max_n) {
                    continue;
                }
                if (static_cast<u128>(a) * p10d + b != n) {
                    continue;
                }

                seen.insert(n);
            }
        }
    }

    u128 total = 0;
    for (u64 v : seen) {
        total += v;
    }
    return static_cast<u64>(total);
}

void run_validations() {
    assert(brute_sum_digits(4) == 5131);
    assert(solve_digits(4) == 5131);

    for (int nd : {5, 6, 7}) {
        assert(solve_digits(nd) == brute_sum_digits(nd));
    }
}

}  // namespace

int main() {
    run_validations();
    std::cout << solve_digits(16) << '\n';
    return 0;
}

Python

import math

def pow10_int(e):
    return 10 ** e

def sieve_primes(limit):
    composite = [False] * (limit + 1)
    primes = []
    for i in range(2, limit + 1):
        if not composite[i]:
            primes.append(i)
            if i * i <= limit:
                for j in range(i * i, limit + 1, i):
                    composite[j] = True
    return primes

def factor_prime_powers(x, primes):
    parts = []
    for p in primes:
        if p * p > x:
            break
        if x % p != 0:
            continue
        pk = 1
        while x % p == 0:
            x //= p
            pk *= p
        parts.append(pk)
    if x > 1:
        parts.append(x)
    return parts

def mod_inv(a, mod):
    return pow(a, -1, mod)

def crt_merge(r1, m1, r2, m2):
    if m1 == 1:
        return r2 % m2, m2
    inv = mod_inv(m1 % m2, m2)
    diff = (r2 - r1) % m2
    t = (diff * inv) % m2
    mod = m1 * m2
    r = (r1 + m1 * t) % mod
    return r, mod

def solve_digits(nd):
    max_n = pow10_int(nd) - 1
    s_limit = int(math.isqrt(max_n))
    
    pow10 = [1] * (nd + 1)
    for i in range(1, nd + 1):
        pow10[i] = pow10[i - 1] * 10
        
    sieve_limit = int(math.isqrt(pow10[nd - 1] - 1)) + 10
    primes = sieve_primes(sieve_limit)
    
    seen = set()
    
    for d in range(1, nd):
        mod = pow10[d] - 1
        prime_powers = factor_prime_powers(mod, primes)
        
        residues = [(0, 1)]
        for pk in prime_powers:
            nxt = []
            for r, m in residues:
                nxt.append(crt_merge(r, m, 0, pk))
                nxt.append(crt_merge(r, m, 1 % pk, pk))
            residues = nxt
            
        s_max = min(s_limit, pow10[d] - 1)
        b_low = 1 if d == 1 else pow10[d - 1]
        p10d = pow10[d]
        
        for r0, step in residues:
            s = r0
            if s <= 1:
                need = 2 - s
                k = (need + step - 1) // step
                if k > (s_max - s) // step:
                    continue
                s += k * step
                
            while s <= s_max:
                ss = s * (s - 1)
                if ss % mod != 0:
                    s += step
                    continue
                    
                a = ss // mod
                if a == 0 or a >= s:
                    s += step
                    continue
                b = s - a
                
                if b < b_low or b >= p10d:
                    s += step
                    continue
                    
                n = s * s
                if n > max_n:
                    s += step
                    continue
                if a * p10d + b != n:
                    s += step
                    continue
                    
                seen.add(n)
                s += step
                
    return sum(seen)

def solve():
    return str(solve_digits(16))

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

Java

import java.util.*;

public class Euler932 {

    static long pow10Int(int e) {
        long v = 1;
        for (int i = 0; i < e; ++i)
            v *= 10;
        return v;
    }

    static List<Integer> sievePrimes(int limit) {
        boolean[] composite = new boolean[limit + 1];
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= limit; ++i) {
            if (!composite[i]) {
                primes.add(i);
                if ((long) i * i <= limit) {
                    for (int j = i * i; j <= limit; j += i) {
                        composite[j] = true;
                    }
                }
            }
        }
        return primes;
    }

    static List<Long> factorPrimePowers(long x, List<Integer> primes) {
        List<Long> parts = new ArrayList<>();
        for (int p : primes) {
            if ((long) p * p > x)
                break;
            if (x % p != 0)
                continue;
            long pk = 1;
            while (x % p == 0) {
                x /= p;
                pk *= p;
            }
            parts.add(pk);
        }
        if (x > 1)
            parts.add(x);
        return parts;
    }

    static long[] egcd(long a, long b) {
        if (b == 0)
            return new long[] { a, 1, 0 };
        long[] g = egcd(b, a % b);
        long x1 = g[1];
        long y1 = g[2];
        return new long[] { g[0], y1, x1 - (a / b) * y1 };
    }

    static long modInv(long a, long mod) {
        long[] g = egcd(a, mod);
        long r = g[1] % mod;
        if (r < 0)
            r += mod;
        return r;
    }

    static class Pair {
        long r, m;

        Pair(long r, long m) {
            this.r = r;
            this.m = m;
        }
    }

    static Pair crtMerge(long r1, long m1, long r2, long m2) {
        if (m1 == 1)
            return new Pair(r2 % m2, m2);
        long inv = modInv(m1 % m2, m2);
        long diff = (r2 - r1) % m2;
        if (diff < 0)
            diff += m2;

        long t = (diff * inv) % m2;
        long mod = m1 * m2;
        long r = (r1 + m1 * t) % mod;
        return new Pair(r, mod);
    }

    static long solveDigits(int nd) {
        long maxN = pow10Int(nd) - 1;
        long sLimit = (long) Math.sqrt(maxN);

        long[] pow10 = new long[nd + 1];
        pow10[0] = 1;
        for (int i = 1; i <= nd; ++i)
            pow10[i] = pow10[i - 1] * 10;

        int sieveLimit = (int) Math.sqrt(pow10[nd - 1] - 1) + 10;
        List<Integer> primes = sievePrimes(sieveLimit);

        Set<Long> seen = new HashSet<>();

        for (int d = 1; d < nd; ++d) {
            long mod = pow10[d] - 1;
            List<Long> primePowers = factorPrimePowers(mod, primes);

            List<Pair> residues = new ArrayList<>();
            residues.add(new Pair(0, 1));

            for (long pk : primePowers) {
                List<Pair> next = new ArrayList<>();
                for (Pair rm : residues) {
                    next.add(crtMerge(rm.r, rm.m, 0, pk));
                    next.add(crtMerge(rm.r, rm.m, 1 % pk, pk));
                }
                residues = next;
            }

            long sMax = Math.min(sLimit, pow10[d] - 1);
            long bLow = d == 1 ? 1 : pow10[d - 1];
            long p10d = pow10[d];

            for (Pair p : residues) {
                long r0 = p.r;
                long step = p.m;
                long s = r0;
                if (s <= 1) {
                    long need = 2 - s;
                    long k = (need + step - 1) / step;
                    if (k > (sMax - s) / step)
                        continue;
                    s += k * step;
                }

                for (; s <= sMax; s += step) {
                    long ssMod = (s % mod) * ((s - 1) % mod) % mod;
                    if (ssMod != 0)
                        continue;

                    long ssDouble = s * (s - 1); // safe because max s limit is 10^8
                    long a = ssDouble / mod;
                    if (a == 0 || a >= s)
                        continue;
                    long b = s - a;

                    if (b < bLow || b >= p10d)
                        continue;

                    long n = s * s;
                    if (n > maxN)
                        continue;
                    if (a * p10d + b != n)
                        continue;

                    seen.add(n);
                }
            }
        }

        long total = 0;
        for (long v : seen) {
            total += v;
        }
        return total;
    }

    public static String solve() {
        return Long.toString(solveDigits(16));
    }

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