Problem 271: Modular Cubes, Part 1

Project Euler Problem 271 Solution

EulerSolve provides an optimized solution for Project Euler Problem 271, Modular Cubes, Part 1, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We want all residues $x$ satisfying $$x^3\equiv1\pmod n,$$ and then we sum all nontrivial solutions, meaning all residues with $1<x<n$. The default value used by the code is $$n=13082761331670030=2\cdot3\cdot5\cdot7\cdot11\cdot13\cdot17\cdot19\cdot23\cdot29\cdot31\cdot37\cdot41\cdot43,$$ which is squarefree. Mathematical Approach 1. Split the Problem Prime by Prime Because $n$ is squarefree, we may write $$n=\prod_{i=1}^{k}p_i$$ with distinct primes $p_i$. Then $$x^3\equiv1\pmod n$$ is equivalent to the simultaneous system $$x^3\equiv1\pmod{p_i}\qquad(i=1,\dots,k).$$ So the whole problem becomes: 1. find all cube roots of unity modulo each prime $p_i$, 2. combine one local choice from each prime modulus using the Chinese Remainder Theorem. 2. How Many Cube Roots Exist Modulo a Prime? For a prime $p$, the nonzero residues modulo $p$ form a cyclic multiplicative group of size $p-1$. In a cyclic group of order $p-1$, the equation $$u^3=1$$ has exactly $$\gcd(3,p-1)$$ solutions. Therefore: 1. if $p\equiv2\pmod3$, then $\gcd(3,p-1)=1$, so the only root is $1$, 2. if $p\equiv1\pmod3$, then $\gcd(3,p-1)=3$, so there are three roots. The code does not use a special generator argument; since every prime factor of the target $n$ is small, it simply brute-forces all residues $1\le x<p$ and keeps those with $x^3\equiv1\pmod p$. 3....

Detailed mathematical approach

Problem Summary

For a positive integer $n$, $S(n)$ is the sum of all integers $x$ with $1\lt x\lt n$ and $x^3\equiv1\pmod{n}$, that is, of all $x$ in this range for which $x^3-1$ is divisible by $n$. For $n=91$ there are eight such numbers, $9,16,22,29,53,74,79,81$, and $S(91)=363$. We must find $S(n)$ for $n=13082761331670030$. ^[1]

The whole solution in one sentence: this $n$ is the product of the fourteen primes from $2$ to $43$, so choosing a cube root of $1$ modulo $n$ is the same as choosing one cube root of $1$ modulo each of these primes, there are only $3^6=729$ such choices, and the Chinese remainder theorem turns each choice into one number below $n$; $S(n)$ is the sum of these numbers except $1$.

Mathematical Approach

We look at small moduli first, then split the congruence modulo $n$ into one congruence for each prime factor, find all cube roots of $1$ modulo a prime, and glue the roots for the different primes back together.

1. Small moduli

Modulo $7$, the cubes of $1,2,\dots,6$ leave the remainders $1,1,6,1,6,6$; for instance $4^3=64=9\cdot7+1$. So $x^3\equiv1\pmod7$ has the three solutions $1,2,4$, and $S(7)=2+4=6$. Modulo $13$ the solutions are $1,3,9$; for instance $3^3=27=2\cdot13+1$ and $9^3=729=56\cdot13+1$. Modulo $11$ the cubes of $1,\dots,10$ are $1,8,5,9,4,7,2,6,3,10$: every nonzero remainder appears once, so $1$ is the only cube root of $1$. ^[5]

The solutions modulo $91=7\cdot13$ fit these lists: $9$ leaves the remainders $2$ and $9$ on division by $7$ and $13$, $16$ leaves $2$ and $3$, and so on; together with $x=1$ there are $9=3\cdot3$ solutions, one for each pair of roots. Sections 3 and 4 explain why.

2. The modulus

Trial division shows that

$$\begin{aligned}n&=13082761331670030\\&=2\cdot3\cdot5\cdot7\cdot11\cdot13\cdot17\\&\quad\cdot19\cdot23\cdot29\cdot31\cdot37\cdot41\cdot43,\end{aligned}$$

the product of the first fourteen primes, also called the primorial $43\#$. ^[6] Every prime factor occurs once, so $n$ is squarefree. Six of the primes leave the remainder $1$ on division by $3$, namely $7,13,19,31,37,43$; the other eight are $3$ and the seven primes $2,5,11,17,23,29,41$ of the form $3k+2$.

3. One congruence for each prime

The Chinese remainder theorem ^[2] says: if the moduli $m_1,\dots,m_k$ are pairwise coprime with product $M$, then for every choice of remainders $a_i$ modulo $m_i$ there is exactly one $x$ with $0\le x\lt M$ and $x\equiv a_i\pmod{m_i}$ for all $i$. In particular, a number is divisible by $M$ exactly when it is divisible by every $m_i$, since $0$ is the only residue whose remainders are all $0$.

Claim. Let $n=p_1p_2\cdots p_k$ be a product of different primes. Then $x^3\equiv1\pmod n$ holds exactly when $x^3\equiv1\pmod{p_i}$ for every $i$, and the solutions modulo $n$ correspond one-to-one to the lists $(r_1,\dots,r_k)$ in which each $r_i$ is a cube root of $1$ modulo $p_i$.

Proof. Different primes are coprime, so the first statement is the divisibility fact applied to $x^3-1$. The remainder of $x^3$ modulo $p_i$ depends only on the remainder $r_i$ of $x$, so $x$ is a solution exactly when every $r_i$ is a cube root of $1$ modulo $p_i$, and by the theorem every list of remainders comes from exactly one residue modulo $n$. $\square$

So the number of solutions modulo $n$ is the product of the numbers of local roots. The primes must be different: modulo $9$ the cube roots of $1$ are $1,4,7$, but modulo $3$ there is only $1$.

4. Cube roots of 1 modulo a prime

Claim. Let $p$ be a prime. If $p\equiv1\pmod3$, the congruence $x^3\equiv1\pmod p$ has exactly three solutions modulo $p$; otherwise its only solution is $x\equiv1$.

Proof. A solution $x$ is not divisible by $p$, so Fermat's little theorem gives $x^{p-1}\equiv1\pmod p$. ^[3] If $p=3$, the form $x^p\equiv x$ of the same theorem gives $x\equiv x^3\equiv1$. If $p\equiv2\pmod3$, which includes $p=2$, then $3$ divides $2p-1$, and

$$\begin{aligned}x&\equiv x\,(x^{p-1})^2=x^{2p-1}\\&=(x^3)^{(2p-1)/3}\equiv1\pmod p.\end{aligned}$$

For $p=11$ this reads $x\equiv x^{21}=(x^3)^7\equiv1$, which explains the list of Section 1. If $p\equiv1\pmod3$, write $p-1=3k$. Then

$$\begin{aligned}x^{3k}-1&=(x^3-1)\,Q(x),\\Q(x)&=\sum_{j=0}^{k-1}x^{3j}.\end{aligned}$$

By Fermat's theorem $p$ divides the left side at all $p-1$ nonzero residues, and since $p$ is prime, at each of them it divides $x^3-1$ or $Q(x)$. Lagrange's theorem ^[4] says that a polynomial of degree $d$ whose coefficients are not all divisible by $p$ has at most $d$ roots modulo $p$. All coefficients of $Q$ are $1$ and its degree is $3k-3=p-4$, so $p$ divides $Q(x)$ for at most $p-4$ of the $p-1$ nonzero residues. For each of the at least $3$ others, $p$ must divide $x^3-1$; by Lagrange's theorem this happens for at most $3$ residues, hence for exactly $3$. $\square$

Since $x^3-1=(x-1)(x^2+x+1)$, the two roots other than $1$ are the roots of $x^2+x+1$. If $\omega$ is one of them, the other is $\omega^2$: it is a cube root of $1$, and $\omega^2\equiv1$ or $\omega^2\equiv\omega$ would force $\omega\equiv1$. From $\omega^2+\omega+1\equiv0$ the two add up to $-1$ modulo $p$; both lie between $2$ and $p-1$, so their sum is exactly $p-1$.

5. The roots for the fourteen primes

For the eight primes $2,3,5,11,17,23,29,41$ the only cube root of $1$ is $1$. For the other six, the roots besides $1$ are

$2,4$ modulo $7$ and $3,9$ modulo $13$;
$7,11$ modulo $19$ and $5,25$ modulo $31$;
$10,26$ modulo $37$ and $6,36$ modulo $43$.

Each pair adds up to $p-1$, and its second number is the square of the first modulo $p$; for example $7^2=49=2\cdot19+11$. By Section 3 there are $3^6=729$ solutions modulo $n$, all congruent to $1$ modulo each of the eight primes with a single root. One of them is $x=1$, and $x=0$ is never a solution, so $S(n)$ is a sum of $728$ terms. This agrees with OEIS A060839, ^[7] which counts the solutions of $x^3\equiv1\pmod n$ by a multiplicative rule: a prime power factor $p^e$ contributes $3$ if $p\equiv1\pmod3$ and $1$ if $p\equiv2\pmod3$, and a factor $3^e$ contributes $1$ for $e=1$ but $3$ for $e\ge2$.

6. Gluing two congruences

Suppose we know $x\equiv a_1\pmod{m_1}$ with $0\le a_1\lt m_1$ and want also $x\equiv a_2\pmod{m_2}$, where $\gcd(m_1,m_2)=1$. The numbers $a_1+m_1t$ with $t=0,1,\dots,m_2-1$ all satisfy the first congruence, and no two of them are congruent modulo $m_2$, since $m_2$ would have to divide $m_1(t-t')$, hence $t-t'$. So exactly one of them satisfies the second congruence:

$$\begin{aligned}x&=a_1+m_1t,\\t&\equiv(a_2-a_1)\,m_1^{-1}\pmod{m_2},\end{aligned}$$

with $0\le t\lt m_2$, so that $0\le x\lt m_1m_2$. Here $m_1^{-1}$ is the inverse of $m_1$ modulo $m_2$, a number $u$ with $m_1u\equiv1\pmod{m_2}$; the extended Euclidean algorithm ^[8] finds integers $u,v$ with $m_1u+m_2v=1$. This is the two-modulus formula of the Chinese remainder theorem. ^[2] Starting from $x\equiv0\pmod1$ and gluing in one prime after another, we obtain a solution modulo the product of the primes used so far.

7. Worked example: $n=91$

The primes are $7$ and $13$, with the roots $1,2,4$ and $1,3,9$. Gluing $x\equiv0\pmod1$ with $x\equiv r_7\pmod7$ gives $t=r_7$ and $x=r_7$. The second step has $m_1=7$, $m_2=13$ and $m_1^{-1}=2$, because $7\cdot2=14=13+1$. For $r_7=2$ and $r_{13}=3$ we get $t\equiv(3-2)\cdot2=2$ and $x=2+7\cdot2=16$; indeed $16^3=4096=45\cdot91+1$. For $r_7=2$ and $r_{13}=1$ the difference is negative: $t\equiv(1-2)\cdot2\equiv11\pmod{13}$, so $x=2+7\cdot11=79$. All nine pairs $(r_7,r_{13})$ give

$(1,1)\to1$, $(1,3)\to29$, $(1,9)\to22$;
$(2,1)\to79$, $(2,3)\to16$, $(2,9)\to9$;
$(4,1)\to53$, $(4,3)\to81$, $(4,9)\to74$.

Leaving out $x=1$, the sum is $9+16+22+29+53$${}+74+79+81=363$, the value $S(91)$ of the statement.

8. The complete method and why it is correct

We take the primes $p_1\lt p_2\lt\dots\lt p_k$ of $n$ in turn and glue every root modulo the next prime onto every number obtained so far. Claim. After $j$ primes, the numbers obtained are exactly the solutions of $x^3\equiv1$ modulo $P_j=p_1\cdots p_j$ with $0\le x\lt P_j$, each once. Proof. For $j=0$ the only number is $0$, the single residue modulo $1$. If the claim holds for $j$, Section 6 turns every pair of an old number and a root modulo $p_{j+1}$ into the unique residue modulo $P_{j+1}$ with these two remainders, and by Section 3, applied to $P_{j+1}$, these residues are exactly the solutions modulo $P_{j+1}$. $\square$

After all fourteen primes we have the $729$ solutions with $0\le x\lt n$, and $S(n)$ is their sum without $x=1$; as a sum of $728$ numbers below $n$, it is less than $728\,n$. Testing all $x\lt n$ instead would take about $1.3\cdot10^{16}$ steps.

How the Code Works

The C++ program ^[10] starts in main. parse_arguments accepts --skip-checkpoints and --n= followed by decimal digits, read without an overflow check; an unknown argument or a value $n\le2$ gives exit code $1$. Unless skipped, run_checkpoints requires solve(91) to return $363$ and solve(7) to return $6$, or ends the program with exit code $2$. Then main prints solve(n).

solve follows Sections 2 to 8. distinct_prime_factors removes the factor $2$ and then tries odd $p$ while $p^2$ does not exceed the part of $n$ still unfactored, dividing each prime out completely; a leftover above $1$ is a prime. For our $n$ the trial divisors are the odd numbers from $3$ to $41$, and the leftover $43$ is added at the end. For every prime $p$, roots_by_prime receives the $x$ with $1\le x\lt p$ and mod_pow(x, 3, p) == 1: the program tests all $p-1$ candidates, and Section 4 predicts what it finds. mod_pow is square-and-multiply exponentiation ^[9] with 128-bit products. ^[11]

The recursion of Section 8 is the generic lambda dfs, which receives itself as self, the index idx of the next prime, the current residue and the product modulus of the primes glued in so far; the first call starts from $0$ modulo $1$. For every root of the next prime it calls crt_pair(residue, modulus, root, p) and recurses with modulus * p. At a leaf it adds residue to sum if it lies strictly between $1$ and $n$, which drops exactly $x=1$. crt_pair is the formula of Section 6: mod_inverse obtains $m_1^{-1}\bmod m_2$ from the recursive extended_gcd and would return $0$ for moduli that are not coprime, which cannot occur here; a2 + m2 - a1 % m2 keeps the difference nonnegative in unsigned arithmetic, and $t$ and $m_1t$ are computed in 128 bits. All calls at one depth have the same modulus and p, so the same inverse is computed again at each of them, which is wasteful but harmless at this size. The variable sum is an unsigned 64-bit integer, and $728\,n\approx9.52\cdot10^{18}$ is below $2^{64}\approx1.84\cdot10^{19}$, so it cannot overflow.

The Python program has the same structure, with solve_for_n and a nested dfs that updates ans through nonlocal; it uses the built-in pow(x, 3, p) and unbounded integers ^[11] and has no checkpoints or options. The Java program keeps sum, primes, rootsByPrime and N in static fields and forms modular products in multiplyMod with BigInteger, again without checkpoints. Its sum is a signed long, whose maximum $2^{63}-1\approx9.22\cdot10^{18}$ is smaller than $728\,n$, so the bound does not exclude an overflow. The exact total does fit: an overflow would have made the result negative, and Java prints the same number as Python.

All three programs need $n$ to be squarefree, because they solve the congruence only modulo each distinct prime. For a value of --n= that is not squarefree the C++ result can be wrong: for $n=9$ it prints $0$, although $4^3=64=7\cdot9+1$ and $7^3=343=38\cdot9+1$ give $S(9)=11$. The modulus of this problem is squarefree, so its result is not affected.

Complexity and Verification

For our $n$, factoring removes the factor $2$ and tests the $20$ odd candidates from $3$ to $41$, and the root search makes $\sum(p-1)=267$ calls of mod_pow. The recursion makes $1405$ calls of dfs, $729$ of them leaves, and $1404$ calls of crt_pair, each running Euclid's algorithm on numbers up to $43$ with at most $8$ nested calls of extended_gcd; the two checkpoints of the C++ program add $15$ more calls of crt_pair. The stored data are the $14$ primes and their root lists, and dfs is at most $15$ calls deep, plus the short stack of extended_gcd. On the test machine (GCC 16 with -O2, OpenJDK 17, Python 3.11) each program runs in at most a few hundredths of a second, mostly start-up time. In general, if $n$ has $k$ distinct prime factors and $R$ solutions, the recursion makes at most $kR$ calls of crt_pair, each computing an inverse in $O(\log p)$ steps. ^[8] Factoring needs up to about $\sqrt n/2$ trial divisions and the root search $\sum(p-1)$ cubings, so the method is fast only when all prime factors of $n$ are small, as here.

The C++ program checks $S(91)=363$ and $S(7)=6$ before the main computation. For this article we also ran checks that are not part of the published solutions:

a reimplementation of the published method agrees with a brute force over all $x$ for every squarefree $n$ with $3\le n\le5000$, and disagrees for $40$ of the non-squarefree $n\le200$, the smallest being $9$;
for our $n$ the $729$ leaves are distinct, and each satisfies $x^3\equiv1\pmod n$ in exact arithmetic;
a brute force over all $85276009$ residues modulo $7\cdot13\cdot19\cdot31\cdot37\cdot43$ finds exactly $729$ cube roots of $1$, and combining each with $x\equiv1$ modulo the product $153416670$ of the other eight primes, as Section 5 requires, gives the same sum;
the numbers of solutions for $n\le100$ agree with OEIS A060839, and the three programs print the same number.

Footnotes and References

The references below supply the background facts used above; the case analysis for the primes, the worked example and all counts are derived or computed in this article.

Project Euler 271 — Modular Cubes, Part 1. The official statement defines S(n) as the sum of the integers x with 1 < x < n and x³ ≡ 1 mod n, lists the eight values of x for n = 91 with S(91) = 363, and asks for S(13082761331670030).
Wikipedia — Chinese remainder theorem. For pairwise coprime moduli with product N, any given remainders belong to exactly one x with 0 ≤ x < N. For two moduli, Bézout coefficients with m₁n₁ + m₂n₂ = 1 give the solution x = a₁ + (a₂ − a₁)m₁n₁, and larger systems are solved two congruences at a time.
Wikipedia — Fermat's little theorem. For a prime p and any integer a, a^p ≡ a (mod p), and a^(p−1) ≡ 1 (mod p) when p does not divide a.
Wikipedia — Lagrange's theorem (number theory). For a prime p, a polynomial f with integer coefficients that are not all divisible by p has at most deg f roots modulo p.
Wikipedia — Root of unity modulo n. A k-th root of unity modulo n is an x with x^k ≡ 1 (mod n); the examples include the cube roots of unity 1, 2, 4 modulo 7 and the single cube root of unity 1 modulo 11.
OEIS A002110 — Primorial numbers. The products of the first n primes; the term with index 14 is 13082761331670030, the product of the primes from 2 to 43.
OEIS A060839 — Number of solutions to x^3 == 1 (mod n). The sequence is multiplicative with a(p^e) = 3 for primes p of the form 3k + 1, a(p^e) = 1 for primes of the form 3k + 2, a(3) = 1 and a(3^e) = 3 for e ≥ 2; its data give a(9) = 3 and a(91) = 9.
Wikipedia — Extended Euclidean algorithm. Besides gcd(a, b) it computes Bézout coefficients x, y with ax + by = gcd(a, b); for coprime a and b, x is the inverse of a modulo b. By Lamé's theorem, stated in Wikipedia — Euclidean algorithm, the number of division steps never exceeds five times the number of decimal digits of the smaller number.
Wikipedia — Modular exponentiation. The right-to-left binary method squares the base for every bit of the exponent and multiplies it into the result when the bit is 1, in time O(log e).
C++ — Euler271.cpp. The linked immutable revision contains distinct_prime_factors, crt_pair, solve with its recursive dfs and run_checkpoints. It is the implementation record for this article.
GCC manual — 128-bit Integers documents the __int128 extension; Java — Primitive Data Types gives the range −2⁶³ to 2⁶³ − 1 of long, and Java — BigInteger provides arbitrary-precision integers; Python — numeric types states that integers have unlimited precision.

Problem 271 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = long long;
using u64 = std::uint64_t;
using i128 = __int128_t;

struct Options {
    u64 n = 13082761331670030ULL;
    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 * 10ULL + 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, "--n=", options.n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n > 2;
}

u64 mod_pow(u64 base, u64 exp, u64 mod) {
    u64 result = 1 % mod;
    u64 cur = base % mod;
    while (exp > 0) {
        if ((exp & 1ULL) != 0ULL) {
            result = static_cast<u64>((static_cast<i128>(result) * cur) % mod);
        }
        cur = static_cast<u64>((static_cast<i128>(cur) * cur) % mod);
        exp >>= 1U;
    }
    return result;
}

std::vector<u64> distinct_prime_factors(u64 n) {
    std::vector<u64> factors;
    if ((n & 1ULL) == 0ULL) {
        factors.push_back(2ULL);
        while ((n & 1ULL) == 0ULL) {
            n >>= 1U;
        }
    }
    for (u64 p = 3; p * p <= n; p += 2ULL) {
        if (n % p != 0ULL) {
            continue;
        }
        factors.push_back(p);
        while (n % p == 0ULL) {
            n /= p;
        }
    }
    if (n > 1ULL) {
        factors.push_back(n);
    }
    return factors;
}

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

u64 mod_inverse(const u64 a, const u64 mod) {
    i64 x = 0;
    i64 y = 0;
    const i64 g = extended_gcd(static_cast<i64>(a), static_cast<i64>(mod), x, y);
    if (g != 1) {
        return 0;
    }
    i64 r = x % static_cast<i64>(mod);
    if (r < 0) {
        r += static_cast<i64>(mod);
    }
    return static_cast<u64>(r);
}

u64 crt_pair(const u64 a1, const u64 m1, const u64 a2, const u64 m2) {
    const u64 inv = mod_inverse(m1 % m2, m2);
    const u64 t = static_cast<u64>((static_cast<i128>((a2 + m2 - (a1 % m2)) % m2) * inv) % m2);
    return static_cast<u64>(a1 + static_cast<i128>(m1) * t);
}

u64 solve(const u64 n) {
    const std::vector<u64> primes = distinct_prime_factors(n);

    std::vector<std::vector<u64>> roots_by_prime;
    roots_by_prime.reserve(primes.size());
    for (u64 p : primes) {
        std::vector<u64> roots;
        for (u64 x = 1; x < p; ++x) {
            if (mod_pow(x, 3, p) == 1ULL) {
                roots.push_back(x);
            }
        }
        roots_by_prime.push_back(std::move(roots));
    }

    u64 sum = 0;
    const auto dfs = [&](auto&& self, std::size_t idx, u64 residue, u64 modulus) -> void {
        if (idx == primes.size()) {
            if (residue > 1ULL && residue < n) {
                sum += residue;
            }
            return;
        }
        const u64 p = primes[idx];
        for (u64 root : roots_by_prime[idx]) {
            const u64 next_residue = crt_pair(residue, modulus, root, p);
            self(self, idx + 1, next_residue, modulus * p);
        }
    };
    dfs(dfs, 0, 0ULL, 1ULL);
    return sum;
}

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

Python

import math

def distinct_prime_factors(n):
    factors = []
    if n % 2 == 0:
        factors.append(2)
        while n % 2 == 0:
            n //= 2
    p = 3
    while p * p <= n:
        if n % p == 0:
            factors.append(p)
            while n % p == 0:
                n //= p
        p += 2
    if n > 1:
        factors.append(n)
    return factors

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

def mod_inverse(a, mod):
    g, x, y = extended_gcd(a, mod)
    if g != 1: return 0
    return x % mod

def crt_pair(a1, m1, a2, m2):
    inv = mod_inverse(m1 % m2, m2)
    t = (((a2 - (a1 % m2)) % m2 + m2) % m2 * inv) % m2
    return a1 + m1 * t

def solve_for_n(n):
    primes = distinct_prime_factors(n)
    roots_by_prime = []
    for p in primes:
        roots = []
        for x in range(1, p):
            if pow(x, 3, p) == 1:
                roots.append(x)
        roots_by_prime.append(roots)

    ans = 0
    def dfs(idx, residue, modulus):
        nonlocal ans
        if idx == len(primes):
            if 1 < residue < n:
                ans += residue
            return
            
        p = primes[idx]
        for root in roots_by_prime[idx]:
            next_residue = crt_pair(residue, modulus, root, p)
            dfs(idx + 1, next_residue, modulus * p)

    dfs(0, 0, 1)
    return ans

def solve():
    n = 13082761331670030
    ans = solve_for_n(n)
    return str(ans)

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;

public class Euler271 {

    static List<Long> distinctPrimeFactors(long n) {
        List<Long> factors = new ArrayList<>();
        if ((n & 1) == 0) {
            factors.add(2L);
            while ((n & 1) == 0)
                n >>= 1;
        }
        for (long p = 3; p * p <= n; p += 2) {
            if (n % p == 0) {
                factors.add(p);
                while (n % p == 0)
                    n /= p;
            }
        }
        if (n > 1)
            factors.add(n);
        return factors;
    }

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

    static long modInverse(long a, long mod) {
        long[] res = extendedGcd(a, mod);
        if (res[0] != 1)
            return 0;
        long r = res[1] % mod;
        if (r < 0)
            r += mod;
        return r;
    }

    static long modPow(long base, long exp, long mod) {
        long res = 1 % mod;
        long cur = base % mod;
        while (exp > 0) {
            if ((exp & 1) == 1)
                res = multiplyMod(res, cur, mod);
            cur = multiplyMod(cur, cur, mod);
            exp >>= 1;
        }
        return res;
    }

    static long multiplyMod(long a, long b, long mod) {
        return BigInteger.valueOf(a).multiply(BigInteger.valueOf(b)).mod(BigInteger.valueOf(mod)).longValue();
    }

    static long crtPair(long a1, long m1, long a2, long m2) {
        long inv = modInverse(m1 % m2, m2);
        long diff = (a2 - (a1 % m2)) % m2;
        if (diff < 0)
            diff += m2;
        long t = multiplyMod(diff, inv, m2);
        return a1 + m1 * t;
    }

    static long sum = 0;
    static List<Long> primes;
    static List<List<Long>> rootsByPrime;
    static long N;

    static void dfs(int idx, long residue, long modulus) {
        if (idx == primes.size()) {
            if (residue > 1 && residue < N) {
                sum += residue;
            }
            return;
        }
        long p = primes.get(idx);
        for (long root : rootsByPrime.get(idx)) {
            long nextResidue = crtPair(residue, modulus, root, p);
            dfs(idx + 1, nextResidue, modulus * p);
        }
    }

    static long solve(long n) {
        N = n;
        primes = distinctPrimeFactors(n);
        rootsByPrime = new ArrayList<>();
        sum = 0;

        for (long p : primes) {
            List<Long> roots = new ArrayList<>();
            for (long x = 1; x < p; x++) {
                if (modPow(x, 3, p) == 1) {
                    roots.add(x);
                }
            }
            rootsByPrime.add(roots);
        }

        dfs(0, 0, 1);
        return sum;
    }

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