Problem 489: Common Factors Between Two Sequences

Project Euler Problem 489 Solution

EulerSolve provides an optimized solution for Project Euler Problem 489, Common Factors Between Two Sequences, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For fixed positive integers $a$ and $b$, consider the two sequences $$u_n=n^3+b,\qquad v_n=(n+a)^3+b.$$ Define $$M(a,b)=\max_{n\ge 0}\gcd(u_n,v_n).$$ Then $G(a,b)$ is the smallest nonnegative integer $n$ for which $\gcd(u_n,v_n)=M(a,b)$. The final target is $$H(m,n)=\sum_{a=1}^{m}\sum_{b=1}^{n}G(a,b).$$ The C++, Python, and Java implementations do not search $n$ directly. Instead, they determine the largest prime-power contribution that can occur in the gcd, then reconstruct the smallest index that realizes all of those local maxima simultaneously. Mathematical Approach Step 1: A common divisor is a simultaneous congruence condition For any prime power $p^k$, the statement $p^k\mid \gcd(u_n,v_n)$ is equivalent to the pair of congruences $$n^3+b\equiv 0\pmod{p^k},\qquad (n+a)^3+b\equiv 0\pmod{p^k}.$$ So the problem is local at each prime: we need to know for which primes and exponents the two cubics have a common residue class. Step 2: The resultant restricts the relevant primes If two integer polynomials have a common root modulo a prime $p$, then $p$ must divide their resultant. For $$f(x)=x^3+b,\qquad g(x)=(x+a)^3+b,$$ a direct Sylvester-determinant computation gives $$\operatorname{Res}_x(f,g)=a^3(a^6+27b^2).$$ Therefore every prime that can appear in the maximal gcd must divide $a^3(a^6+27b^2)$....

Detailed mathematical approach

Problem Summary

For fixed positive integers $a$ and $b$, consider the two sequences

$$u_n=n^3+b,\qquad v_n=(n+a)^3+b.$$

Define

$$M(a,b)=\max_{n\ge 0}\gcd(u_n,v_n).$$

Then $G(a,b)$ is the smallest nonnegative integer $n$ for which $\gcd(u_n,v_n)=M(a,b)$. The final target is

$$H(m,n)=\sum_{a=1}^{m}\sum_{b=1}^{n}G(a,b).$$

The C++, Python, and Java implementations do not search $n$ directly. Instead, they determine the largest prime-power contribution that can occur in the gcd, then reconstruct the smallest index that realizes all of those local maxima simultaneously.

Mathematical Approach

Step 1: A common divisor is a simultaneous congruence condition

For any prime power $p^k$, the statement $p^k\mid \gcd(u_n,v_n)$ is equivalent to the pair of congruences

$$n^3+b\equiv 0\pmod{p^k},\qquad (n+a)^3+b\equiv 0\pmod{p^k}.$$

So the problem is local at each prime: we need to know for which primes and exponents the two cubics have a common residue class.

Step 2: The resultant restricts the relevant primes

If two integer polynomials have a common root modulo a prime $p$, then $p$ must divide their resultant. For

$$f(x)=x^3+b,\qquad g(x)=(x+a)^3+b,$$

a direct Sylvester-determinant computation gives

$$\operatorname{Res}_x(f,g)=a^3(a^6+27b^2).$$

Therefore every prime that can appear in the maximal gcd must divide $a^3(a^6+27b^2)$. This is why the implementation factors $a^6+27b^2$ for each pair $(a,b)$ and then adds three times the prime exponents coming from $a$.

Step 3: A quadratic condition modulo $p$

Subtract the two congruences:

$$0\equiv (n+a)^3-n^3=a(3n^2+3an+a^2)\pmod{p^k}.$$

When $p\nmid a$, this forces

$$3n^2+3an+a^2\equiv 0\pmod{p}.$$

For odd primes $p>5$, we may complete the square:

$$4(3n^2+3an+a^2)=3(2n+a)^2+a^2,$$

hence

$$ (2n+a)^2\equiv -\frac{a^2}{3}\pmod{p}. $$

This reduces the base search modulo $p$ to a modular square-root problem. The implementation computes the square root when possible, converts it into up to two candidate residues, and then verifies those candidates in the original cubic congruences. For small primes and degenerate cases with $p\mid a$, it safely enumerates all residues modulo $p$ instead of relying on the quadratic shortcut.

Step 4: Lift solutions to the highest possible prime power

Suppose $r$ is a simultaneous solution modulo $p^t$. Any lift to the next power must have the form

$$n=r+u\,p^t,\qquad u\in\{0,1,\dots,p-1\}.$$

The implementation tests all $p$ descendants and keeps only those that still satisfy both cubic congruences modulo $p^{t+1}$. Repeating this determines the largest exponent $k_p$ for which simultaneous roots still exist. If lifting fails at the next stage, then $p^{k_p}$ is the largest power of $p$ that can divide the gcd.

Step 5: The maximal gcd is the product of the local maxima

Let $p^{e_p}$ be the exponent of $p$ in the resultant, and let $k_p\le e_p$ be the highest exponent actually reached by the lifting process. Then every value $\gcd(u_n,v_n)$ has $p$-adic valuation at most $k_p$, so

$$M(a,b)\le \prod_p p^{k_p}.$$

Conversely, choose one lifted residue class for each surviving prime. Different prime powers are coprime, so the Chinese remainder theorem merges them into a unique class modulo

$$Q(a,b)=\prod_p p^{k_p}.$$

Every $n$ in that class makes both sequences divisible by $p^{k_p}$ for every surviving prime, so

$$\boxed{M(a,b)=Q(a,b).}$$

Thus $G(a,b)$ is simply the smallest nonnegative integer in the CRT-merged residue classes.

Worked Examples

For $(a,b)=(1,1)$, the resultant is

$$1^3(1^6+27\cdot 1^2)=28=2^2\cdot 7.$$

There is no simultaneous root modulo $2$, while modulo $7$ the unique solution is $n\equiv 5\pmod{7}$. Hence the maximal gcd is $7$, and the first index that attains it is

$$G(1,1)=5.$$

A lifting example is $(a,b)=(1,5)$. Here

$$1^3(1^6+27\cdot 5^2)=676=2^2\cdot 13^2.$$

The only surviving prime is $13$. The residue $n\equiv 5\pmod{13}$ lifts further to $n\equiv 161\pmod{169}$, so the maximal gcd is $13^2=169$ and

$$G(1,5)=161.$$

How the Code Works

The implementations first build a reusable prime table, precompute the values $a^6$, $27b^2$, and the prime factorization of each allowed $a$. For each pair $(a,b)$, they factor $a^6+27b^2$, add the prime exponents contributed by $a^3$, and process each candidate prime separately.

For every prime block, the implementation finds the simultaneous residues modulo $p$, lifts them as far as possible, and stores the highest modulus that still has solutions. The prime blocks are then ordered by the number of local residues, which keeps the CRT merge smaller in practice. Every merged residue achieves the maximal gcd, and the smallest merged residue is returned as $G(a,b)$.

Finally, the values are summed over the full $(a,b)$-grid. Since every pair is independent, the C++, Python, and Java implementations parallelize this outer accumulation. Their published checkpoints are

$$G(1,1)=5,\qquad H(5,5)=128878,\qquad H(10,10)=32936544.$$

Complexity Analysis

For one pair $(a,b)$, the first cost is factoring $a^6+27b^2$ with the precomputed prime list. Then each candidate prime needs a local search modulo $p$, followed by at most one lifting stage for each additional exponent. If the surviving residue counts are $r_1,r_2,\dots,r_s$, then the CRT phase examines $r_1r_2\cdots r_s$ combinations in the worst case.

So the runtime is dominated by three ingredients: factoring the resultant part, lifting prime-power roots, and merging local residue classes. Over the full rectangle $1\le a\le m$, $1\le b\le n$, the total runtime is roughly linear in the number of parameter pairs, multiplied by the average local arithmetic cost. Memory usage is modest and is dominated by the temporary residue lists carried through lifting and CRT.

Footnotes and References

Problem page: https://projecteuler.net/problem=489
Resultant of polynomials: Wikipedia — Resultant
Chinese remainder theorem: Wikipedia — Chinese remainder theorem
Modular square root: Wikipedia — Tonelli-Shanks algorithm
Greatest common divisor and valuations: Wikipedia — $p$-adic valuation

Problem 489 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <thread>
#include <utility>
#include <vector>

using u128 = unsigned __int128;

namespace {

struct PrimeSolutions {
    uint64_t mod = 1;
    std::vector<uint64_t> residues;
};

struct Precompute {
    std::vector<int> primes;
    std::vector<uint64_t> a_pow6;
    std::vector<uint64_t> b_term;
    std::vector<std::vector<std::pair<uint64_t, int>>> a_factors;
};

uint64_t mul_mod(uint64_t a, uint64_t b, uint64_t mod) {
    return static_cast<uint64_t>((u128)a * b % mod);
}

uint64_t mod_pow(uint64_t base, uint64_t exp, uint64_t mod) {
    uint64_t result = 1 % mod;
    uint64_t cur = base % mod;
    while (exp > 0) {
        if (exp & 1) result = mul_mod(result, cur, mod);
        cur = mul_mod(cur, cur, mod);
        exp >>= 1;
    }
    return result;
}

uint64_t mod_inv_prime(uint64_t a, uint64_t p) {
    return mod_pow(a, p - 2, p);
}

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

uint64_t mod_inv_coprime(uint64_t a, uint64_t mod) {
    int64_t x = 0;
    int64_t y = 0;
    egcd(static_cast<int64_t>(a), static_cast<int64_t>(mod), x, y);
    int64_t res = x % static_cast<int64_t>(mod);
    if (res < 0) res += static_cast<int64_t>(mod);
    return static_cast<uint64_t>(res);
}

bool mod_sqrt(uint64_t n, uint64_t p, uint64_t& out) {
    if (p == 2) {
        out = n % p;
        return true;
    }
    n %= p;
    if (n == 0) {
        out = 0;
        return true;
    }
    if (mod_pow(n, (p - 1) / 2, p) != 1) return false;
    if (p % 4 == 3) {
        out = mod_pow(n, (p + 1) / 4, p);
        return true;
    }

    uint64_t q = p - 1;
    int s = 0;
    while ((q & 1) == 0) {
        q >>= 1;
        ++s;
    }

    uint64_t z = 2;
    while (mod_pow(z, (p - 1) / 2, p) != p - 1) ++z;

    uint64_t c = mod_pow(z, q, p);
    uint64_t x = mod_pow(n, (q + 1) / 2, p);
    uint64_t t = mod_pow(n, q, p);
    int m = s;
    while (t != 1) {
        int i = 1;
        uint64_t t2 = mul_mod(t, t, p);
        while (t2 != 1) {
            t2 = mul_mod(t2, t2, p);
            ++i;
        }
        uint64_t b = mod_pow(c, 1ULL << (m - i - 1), p);
        x = mul_mod(x, b, p);
        uint64_t b2 = mul_mod(b, b, p);
        t = mul_mod(t, b2, p);
        c = b2;
        m = i;
    }
    out = x;
    return true;
}

bool is_root_mod(uint64_t n, int a, int b, uint64_t mod) {
    u128 nn = n;
    u128 na = nn + static_cast<u128>(a);
    u128 val1 = nn * nn * nn + static_cast<u128>(b);
    u128 val2 = na * na * na + static_cast<u128>(b);
    return (val1 % mod) == 0 && (val2 % mod) == 0;
}

std::vector<int> sieve_primes(int limit) {
    std::vector<bool> is_prime(limit + 1, true);
    std::vector<int> primes;
    for (int i = 2; i <= limit; ++i) {
        if (!is_prime[i]) continue;
        primes.push_back(i);
        if (static_cast<int64_t>(i) * i > limit) continue;
        for (int j = i * i; j <= limit; j += i) {
            is_prime[j] = false;
        }
    }
    return primes;
}

std::vector<std::pair<uint64_t, int>> factorize(uint64_t n, const std::vector<int>& primes) {
    std::vector<std::pair<uint64_t, int>> out;
    uint64_t tmp = n;
    for (int p : primes) {
        if (static_cast<uint64_t>(p) * p > tmp) break;
        if (tmp % p != 0) continue;
        int e = 0;
        while (tmp % p == 0) {
            tmp /= p;
            ++e;
        }
        out.push_back({static_cast<uint64_t>(p), e});
    }
    if (tmp > 1) out.push_back({tmp, 1});
    return out;
}

void add_factor(std::vector<std::pair<uint64_t, int>>& factors, uint64_t p, int e) {
    for (auto& kv : factors) {
        if (kv.first == p) {
            kv.second += e;
            return;
        }
    }
    factors.push_back({p, e});
}

std::vector<uint64_t> solutions_mod_p(uint64_t p, int a, int b) {
    std::vector<uint64_t> sols;
    if (p <= 5 || (a % static_cast<int>(p)) == 0) {
        for (uint64_t n = 0; n < p; ++n) {
            if (is_root_mod(n, a, b, p)) sols.push_back(n);
        }
        return sols;
    }

    if (p == 3) return sols;

    uint64_t inv3 = mod_inv_prime(3, p);
    uint64_t inv2 = mod_inv_prime(2, p);
    uint64_t a_mod = static_cast<uint64_t>(a) % p;
    uint64_t D = (p + p - mul_mod(a_mod, a_mod, p)) % p;
    D = mul_mod(D, inv3, p);

    uint64_t sqrtD = 0;
    if (!mod_sqrt(D, p, sqrtD)) return sols;

    uint64_t roots[2] = {sqrtD, (p - sqrtD) % p};
    for (int i = 0; i < 2; ++i) {
        uint64_t s = roots[i];
        uint64_t tmp = (p + s - a_mod) % p;
        uint64_t n = mul_mod(tmp, inv2, p);
        if (is_root_mod(n, a, b, p)) {
            if (std::find(sols.begin(), sols.end(), n) == sols.end()) sols.push_back(n);
        }
    }
    return sols;
}

PrimeSolutions solve_prime_power(uint64_t p, int max_exp, int a, int b) {
    PrimeSolutions res;
    res.mod = 1;
    res.residues = {0};

    std::vector<uint64_t> sols = solutions_mod_p(p, a, b);
    if (sols.empty()) return res;

    uint64_t mod = p;
    for (int exp = 1; exp < max_exp; ++exp) {
        uint64_t next_mod = mod * p;
        std::vector<uint64_t> next;
        next.reserve(sols.size() * p);
        for (uint64_t r : sols) {
            for (uint64_t t = 0; t < p; ++t) {
                uint64_t n = r + t * mod;
                if (is_root_mod(n, a, b, next_mod)) {
                    next.push_back(n);
                }
            }
        }
        if (next.empty()) break;
        sols.swap(next);
        mod = next_mod;
    }

    res.mod = mod;
    res.residues = std::move(sols);
    return res;
}

uint64_t crt(uint64_t a, uint64_t mod_a, uint64_t b, uint64_t mod_b) {
    uint64_t inv = mod_inv_coprime(mod_a % mod_b, mod_b);
    uint64_t t = (b + mod_b - (a % mod_b)) % mod_b;
    uint64_t k = mul_mod(t, inv, mod_b);
    return static_cast<uint64_t>(a + static_cast<u128>(mod_a) * k);
}

uint64_t compute_G(int a, int b, const Precompute& pre) {
    uint64_t R = pre.a_pow6[a] + pre.b_term[b];
    std::vector<std::pair<uint64_t, int>> factors = factorize(R, pre.primes);
    for (const auto& kv : pre.a_factors[a]) {
        add_factor(factors, kv.first, 3 * kv.second);
    }

    std::vector<PrimeSolutions> blocks;
    blocks.reserve(factors.size());
    for (const auto& kv : factors) {
        PrimeSolutions sol = solve_prime_power(kv.first, kv.second, a, b);
        if (sol.mod > 1) blocks.push_back(std::move(sol));
    }
    if (blocks.empty()) return 0;

    std::sort(blocks.begin(), blocks.end(),
              [](const PrimeSolutions& lhs, const PrimeSolutions& rhs) {
                  return lhs.residues.size() < rhs.residues.size();
              });

    uint64_t mod = 1;
    std::vector<uint64_t> residues = {0};
    for (const auto& block : blocks) {
        std::vector<uint64_t> next;
        next.reserve(residues.size() * block.residues.size());
        for (uint64_t r : residues) {
            for (uint64_t s : block.residues) {
                next.push_back(crt(r, mod, s, block.mod));
            }
        }
        mod *= block.mod;
        residues.swap(next);
    }

    return *std::min_element(residues.begin(), residues.end());
}

uint64_t compute_H(int m, int n, const Precompute& pre, unsigned threads) {
    const int total = m * n;
    if (threads == 0) threads = 1;
    std::atomic<int> index(0);
    std::vector<std::thread> pool;
    std::vector<uint64_t> partial(threads, 0);

    for (unsigned t = 0; t < threads; ++t) {
        pool.emplace_back([&, t]() {
            while (true) {
                int i = index.fetch_add(1);
                if (i >= total) break;
                int a = i / n + 1;
                int b = i % n + 1;
                partial[t] += compute_G(a, b, pre);
            }
        });
    }
    for (auto& th : pool) th.join();

    uint64_t sum = 0;
    for (uint64_t v : partial) sum += v;
    return sum;
}

Precompute build_precompute(int max_a, int max_b) {
    Precompute pre;
    pre.primes = sieve_primes(200000);
    pre.a_pow6.assign(max_a + 1, 0);
    pre.a_factors.assign(max_a + 1, {});
    pre.b_term.assign(max_b + 1, 0);

    for (int a = 1; a <= max_a; ++a) {
        u128 v = 1;
        for (int i = 0; i < 6; ++i) v *= static_cast<u128>(a);
        pre.a_pow6[a] = static_cast<uint64_t>(v);
        pre.a_factors[a] = factorize(static_cast<uint64_t>(a), pre.primes);
    }

    for (int b = 1; b <= max_b; ++b) {
        pre.b_term[b] = 27ULL * static_cast<uint64_t>(b) * static_cast<uint64_t>(b);
    }

    return pre;
}

} // namespace

int main() {
    const int max_a = 18;
    const int max_b = 1900;
    Precompute pre = build_precompute(max_a, max_b);

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) threads = 4;

    std::cout << "--- Validation Checkpoints ---\n";
    uint64_t g11 = compute_G(1, 1, pre);
    std::cout << "G(1, 1) = " << g11 << (g11 == 5 ? " [PASS]" : " [FAIL]") << "\n";

    uint64_t h55 = compute_H(5, 5, pre, threads);
    std::cout << "H(5, 5) = " << h55 << (h55 == 128878ULL ? " [PASS]" : " [FAIL]") << "\n";

    uint64_t h1010 = compute_H(10, 10, pre, threads);
    std::cout << "H(10, 10) = " << h1010 << (h1010 == 32936544ULL ? " [PASS]" : " [FAIL]") << "\n";

    std::cout << "\n--- Final Solution ---\n";
    uint64_t result = compute_H(18, 1900, pre, threads);
    std::cout << "H(18, 1900) = " << result << "\n";
    return 0;
}

Python

import multiprocessing

def mul_mod(a, b, mod):
    return (a * b) % mod

def mod_pow(base, exp, mod):
    return pow(base, exp, mod)

def mod_inv_prime(a, p):
    return pow(a, p - 2, p)

def egcd(a, b):
    if b == 0:
        return a, 1, 0
    g, x1, y1 = egcd(b, a % b)
    x = y1
    y = x1 - y1 * (a // b)
    return g, x, y

def mod_inv_coprime(a, mod):
    g, x, y = egcd(a, mod)
    return (x % mod + mod) % mod

def mod_sqrt(n, p):
    if p == 2:
        return True, n % p
    n %= p
    if n == 0:
        return True, 0
    if pow(n, (p - 1) // 2, p) != 1:
        return False, 0
    if p % 4 == 3:
        return True, pow(n, (p + 1) // 4, p)
        
    q = p - 1
    s = 0
    while (q & 1) == 0:
        q >>= 1
        s += 1
        
    z = 2
    while pow(z, (p - 1) // 2, p) != p - 1:
        z += 1
        
    c = pow(z, q, p)
    x = pow(n, (q + 1) // 2, p)
    t = pow(n, q, p)
    m = s
    while t != 1:
        i = 1
        t2 = mul_mod(t, t, p)
        while t2 != 1:
            t2 = mul_mod(t2, t2, p)
            i += 1
        b = pow(c, 1 << (m - i - 1), p)
        x = mul_mod(x, b, p)
        b2 = mul_mod(b, b, p)
        t = mul_mod(t, b2, p)
        c = b2
        m = i
    return True, x

def is_root_mod(n, a, b, mod):
    val1 = (n * n * n + b) % mod
    na = n + a
    val2 = (na * na * na + b) % mod
    return val1 == 0 and val2 == 0

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

def factorize(n, primes):
    out = []
    tmp = n
    for p in primes:
        if p * p > tmp: break
        if tmp % p == 0:
            e = 0
            while tmp % p == 0:
                tmp //= p
                e += 1
            out.append([p, e])
    if tmp > 1:
        out.append([tmp, 1])
    return out

def add_factor(factors, p, e):
    for kv in factors:
        if kv[0] == p:
            kv[1] += e
            return
    factors.append([p, e])

def solutions_mod_p(p, a, b):
    sols = []
    if p <= 5 or (a % p) == 0:
        for n in range(p):
            if is_root_mod(n, a, b, p):
                sols.append(n)
        return sols

    if p == 3: return sols

    inv3 = mod_inv_prime(3, p)
    inv2 = mod_inv_prime(2, p)
    a_mod = a % p
    D = (p + p - mul_mod(a_mod, a_mod, p)) % p
    D = mul_mod(D, inv3, p)

    ok, sqrtD = mod_sqrt(D, p)
    if not ok: return sols

    roots = [sqrtD, (p - sqrtD) % p]
    for s in roots:
        tmp = (p + s - a_mod) % p
        n = mul_mod(tmp, inv2, p)
        if is_root_mod(n, a, b, p):
            if n not in sols:
                sols.append(n)
    return sols

def solve_prime_power(p, max_exp, a, b):
    sols = solutions_mod_p(p, a, b)
    if not sols: return 1, [0]

    mod = p
    for exp in range(1, max_exp):
        next_mod = mod * p
        nxt = []
        for r in sols:
            for t in range(p):
                n = r + t * mod
                if is_root_mod(n, a, b, next_mod):
                    nxt.append(n)
        if not nxt: break
        sols = nxt
        mod = next_mod

    return mod, sols

def crt(a, mod_a, b, mod_b):
    inv = mod_inv_coprime(mod_a % mod_b, mod_b)
    t = (b + mod_b - (a % mod_b)) % mod_b
    k = mul_mod(t, inv, mod_b)
    return a + mod_a * k

def compute_G(a, b, primes, a_pow6, b_term, a_factors):
    R = a_pow6[a] + b_term[b]
    factors = factorize(R, primes)
    for p, e in a_factors[a]:
        add_factor(factors, p, 3 * e)

    blocks = []
    for p, e in factors:
        mod, sols = solve_prime_power(p, e, a, b)
        if mod > 1:
            blocks.append((mod, sols))
            
    if not blocks: return 0

    blocks.sort(key=lambda x: len(x[1]))

    mod = 1
    residues = [0]
    for b_mod, b_sols in blocks:
        nxt = []
        for r in residues:
            for s in b_sols:
                nxt.append(crt(r, mod, s, b_mod))
        mod *= b_mod
        residues = nxt

    return min(residues)

def worker_func(idx, primes, a_pow6, b_term, a_factors, n):
    max_b = n
    a = idx // max_b + 1
    b = idx % max_b + 1
    return compute_G(a, b, primes, a_pow6, b_term, a_factors)

def solve():
    max_a = 18
    max_b = 1900
    primes = sieve_primes(200000)
    a_pow6 = [0] * (max_a + 1)
    a_factors = [[] for _ in range(max_a + 1)]
    for a in range(1, max_a + 1):
        a_pow6[a] = a**6
        a_factors[a] = factorize(a, primes)
        
    b_term = [0] * (max_b + 1)
    for b in range(1, max_b + 1):
        b_term[b] = 27 * b * b
        
    total_tasks = max_a * max_b
    threads = multiprocessing.cpu_count()
    if threads == 0: threads = 4
    
    tasks = [(i, primes, a_pow6, b_term, a_factors, max_b) for i in range(total_tasks)]
    with multiprocessing.Pool(threads) as pool:
        results = pool.starmap(worker_func, tasks)
        
    total_sum = sum(results)
    return str(total_sum)

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.concurrent.Callable;
import java.util.concurrent.ExecutionException;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;

public class Euler489 {

    private static long mulMod(long a, long b, long mod) {
        BigInteger biA = BigInteger.valueOf(a);
        BigInteger biB = BigInteger.valueOf(b);
        BigInteger biMod = BigInteger.valueOf(mod);
        return biA.multiply(biB).mod(biMod).longValue();
    }

    private static long modPow(long base, long exp, long mod) {
        long result = 1 % mod;
        long cur = base % mod;
        while (exp > 0) {
            if ((exp & 1) != 0)
                result = mulMod(result, cur, mod);
            cur = mulMod(cur, cur, mod);
            exp >>= 1;
        }
        return result;
    }

    private static long modInvPrime(long a, long p) {
        return modPow(a, p - 2, p);
    }

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

    private static long modInvCoprime(long a, long mod) {
        long x = egcd(a, mod)[1];
        long res = x % mod;
        if (res < 0)
            res += mod;
        return res;
    }

    private static long[] modSqrt(long n, long p) {
        if (p == 2)
            return new long[] { 1, n % p };
        n %= p;
        if (n == 0)
            return new long[] { 1, 0 };
        if (modPow(n, (p - 1) / 2, p) != 1)
            return new long[] { 0, 0 };
        if (p % 4 == 3)
            return new long[] { 1, modPow(n, (p + 1) / 4, p) };

        long q = p - 1;
        int s = 0;
        while ((q & 1) == 0) {
            q >>= 1;
            ++s;
        }

        long z = 2;
        while (modPow(z, (p - 1) / 2, p) != p - 1)
            ++z;

        long c = modPow(z, q, p);
        long x = modPow(n, (q + 1) / 2, p);
        long t = modPow(n, q, p);
        int m = s;

        while (t != 1) {
            int i = 1;
            long t2 = mulMod(t, t, p);
            while (t2 != 1) {
                t2 = mulMod(t2, t2, p);
                ++i;
            }
            long b = modPow(c, 1L << (m - i - 1), p);
            x = mulMod(x, b, p);
            long b2 = mulMod(b, b, p);
            t = mulMod(t, b2, p);
            c = b2;
            m = i;
        }
        return new long[] { 1, x };
    }

    private static boolean isRootMod(long n, int a, int b, long mod) {
        BigInteger biN = BigInteger.valueOf(n);
        BigInteger biA = BigInteger.valueOf(a);
        BigInteger biB = BigInteger.valueOf(b);
        BigInteger biMod = BigInteger.valueOf(mod);

        BigInteger n3 = biN.multiply(biN).multiply(biN);
        BigInteger val1 = n3.add(biB).mod(biMod);
        if (!val1.equals(BigInteger.ZERO))
            return false;

        BigInteger na = biN.add(biA);
        BigInteger na3 = na.multiply(na).multiply(na);
        BigInteger val2 = na3.add(biB).mod(biMod);
        return val2.equals(BigInteger.ZERO);
    }

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

    static class Factor {
        long p;
        int e;

        Factor(long p, int e) {
            this.p = p;
            this.e = e;
        }
    }

    private static List<Factor> factorize(long n, List<Integer> primes) {
        List<Factor> out = new ArrayList<>();
        long tmp = n;
        for (int p : primes) {
            if ((long) p * p > tmp)
                break;
            if (tmp % p != 0)
                continue;
            int e = 0;
            while (tmp % p == 0) {
                tmp /= p;
                ++e;
            }
            out.add(new Factor(p, e));
        }
        if (tmp > 1)
            out.add(new Factor(tmp, 1));
        return out;
    }

    private static void addFactor(List<Factor> factors, long p, int e) {
        for (Factor f : factors) {
            if (f.p == p) {
                f.e += e;
                return;
            }
        }
        factors.add(new Factor(p, e));
    }

    private static List<Long> solutionsModP(long p, int a, int b) {
        List<Long> sols = new ArrayList<>();
        if (p <= 5 || (a % p) == 0) {
            for (long n = 0; n < p; ++n) {
                if (isRootMod(n, a, b, p))
                    sols.add(n);
            }
            return sols;
        }

        if (p == 3)
            return sols;

        long inv3 = modInvPrime(3, p);
        long inv2 = modInvPrime(2, p);
        long aMod = a % p;
        long D = (p + p - mulMod(aMod, aMod, p)) % p;
        D = mulMod(D, inv3, p);

        long[] sqrtRes = modSqrt(D, p);
        if (sqrtRes[0] == 0)
            return sols;

        long sqrtD = sqrtRes[1];
        long[] roots = { sqrtD, (p - sqrtD) % p };
        for (long s : roots) {
            long tmp = (p + s - aMod) % p;
            long n = mulMod(tmp, inv2, p);
            if (isRootMod(n, a, b, p)) {
                if (!sols.contains(n))
                    sols.add(n);
            }
        }
        return sols;
    }

    static class PrimeSolutions {
        long mod = 1;
        List<Long> residues = new ArrayList<>();
    }

    private static PrimeSolutions solvePrimePower(long p, int maxExp, int a, int b) {
        PrimeSolutions res = new PrimeSolutions();
        res.residues.add(0L);

        List<Long> sols = solutionsModP(p, a, b);
        if (sols.isEmpty())
            return res;

        long mod = p;
        for (int exp = 1; exp < maxExp; ++exp) {
            long nextMod = mod * p;
            List<Long> next = new ArrayList<>();
            for (long r : sols) {
                for (long t = 0; t < p; ++t) {
                    long n = r + t * mod;
                    if (isRootMod(n, a, b, nextMod)) {
                        next.add(n);
                    }
                }
            }
            if (next.isEmpty())
                break;
            sols = next;
            mod = nextMod;
        }

        res.mod = mod;
        res.residues = sols;
        return res;
    }

    private static long crt(long a, long modA, long b, long modB) {
        long inv = modInvCoprime(modA % modB, modB);
        long t = (b + modB - (a % modB)) % modB;
        long k = mulMod(t, inv, modB);
        BigInteger biA = BigInteger.valueOf(a);
        BigInteger biModA = BigInteger.valueOf(modA);
        BigInteger biK = BigInteger.valueOf(k);
        return biA.add(biModA.multiply(biK)).longValue();
    }

    static class Precompute {
        List<Integer> primes;
        long[] aPow6;
        long[] bTerm;
        List<Factor>[] aFactors;
    }

    private static Precompute buildPrecompute(int maxA, int maxB) {
        Precompute pre = new Precompute();
        pre.primes = sievePrimes(200000);
        pre.aPow6 = new long[maxA + 1];
        pre.aFactors = new ArrayList[maxA + 1];
        pre.bTerm = new long[maxB + 1];

        for (int a = 1; a <= maxA; ++a) {
            BigInteger v = BigInteger.ONE;
            for (int i = 0; i < 6; ++i)
                v = v.multiply(BigInteger.valueOf(a));
            pre.aPow6[a] = v.longValue();
            pre.aFactors[a] = factorize(a, pre.primes);
        }

        for (int b = 1; b <= maxB; ++b) {
            pre.bTerm[b] = 27L * b * b;
        }

        return pre;
    }

    private static long computeG(int a, int b, Precompute pre) {
        long R = pre.aPow6[a] + pre.bTerm[b];
        List<Factor> factors = factorize(R, pre.primes);
        for (Factor f : pre.aFactors[a]) {
            addFactor(factors, f.p, 3 * f.e);
        }

        List<PrimeSolutions> blocks = new ArrayList<>();
        for (Factor f : factors) {
            PrimeSolutions sol = solvePrimePower(f.p, f.e, a, b);
            if (sol.mod > 1)
                blocks.add(sol);
        }
        if (blocks.isEmpty())
            return 0;

        blocks.sort((x, y) -> Integer.compare(x.residues.size(), y.residues.size()));

        long mod = 1;
        List<Long> residues = new ArrayList<>();
        residues.add(0L);

        for (PrimeSolutions block : blocks) {
            List<Long> next = new ArrayList<>();
            for (long r : residues) {
                for (long s : block.residues) {
                    next.add(crt(r, mod, s, block.mod));
                }
            }
            mod *= block.mod;
            residues = next;
        }

        return Collections.min(residues);
    }

    private static long computeH(int maxA, int maxB, Precompute pre) throws Exception {
        int total = maxA * maxB;
        int threads = Runtime.getRuntime().availableProcessors();
        if (threads <= 0)
            threads = 4;

        ExecutorService executor = Executors.newFixedThreadPool(threads);
        List<Future<Long>> futures = new ArrayList<>();

        int chunk = (total + threads - 1) / threads;
        for (int t = 0; t < threads; t++) {
            final int start = t * chunk;
            final int end = Math.min(total, start + chunk);
            futures.add(executor.submit(() -> {
                long sum = 0;
                for (int i = start; i < end; i++) {
                    int a = i / maxB + 1;
                    int b = i % maxB + 1;
                    sum += computeG(a, b, pre);
                }
                return sum;
            }));
        }

        long totalSum = 0;
        for (Future<Long> f : futures) {
            totalSum += f.get();
        }
        executor.shutdown();
        return totalSum;
    }

    public static void main(String[] args) throws Exception {
        int maxA = 18;
        int maxB = 1900;
        Precompute pre = buildPrecompute(maxA, maxB);
        long result = computeH(maxA, maxB, pre);
        System.out.println(result);
    }
}