Problem 302: Strong Achilles Numbers

Project Euler Problem 302 Solution

EulerSolve provides an optimized solution for Project Euler Problem 302, Strong Achilles Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A strong Achilles number is an integer $n$ such that both $n$ and $\varphi(n)$ are Achilles numbers. An integer is Achilles when it is powerful, but not a perfect power. The problem asks for the number of such $n$ below $10^{18}$. Mathematical Approach Write $$n=\prod_{i=1}^r p_i^{a_i},\qquad a_i\ge 1.$$ 1) Achilles condition in exponent language The number $n$ is powerful exactly when every exponent is at least $2$: $$a_i\ge 2\qquad \text{for all }i.$$ It is a perfect $k$-th power exactly when all exponents are divisible by the same $k \ge 2$, equivalently when $$\gcd(a_1,a_2,\dots,a_r)\ge 2.$$ Therefore $$n\text{ is Achilles}\iff a_i\ge 2\ \forall i\quad\text{and}\quad \gcd(a_1,\dots,a_r)=1.$$ This turns the first half of the problem into a statement about exponent vectors rather than about the decimal value of $n$. 2) What the totient contributes Euler's product formula gives $$\varphi(n)=\prod_{i=1}^r p_i^{a_i-1}(p_i-1).$$ Fix a prime $q$. Its exponent inside $\varphi(n)$ is $$v_q(\varphi(n))=\sum_{i=1}^r v_q(p_i-1)+\sum_{\substack{1\le i\le r\\p_i=q}}(a_i-1).$$ So every exponent of $\varphi(n)$ comes from two sources: 1. factors already forced by the numbers $p_i-1$; 2. the self-contribution $a_i-1$ when the prime $q$ itself appears in $n$....

Detailed mathematical approach

Problem Summary

A strong Achilles number is an integer $n$ such that both $n$ and $\varphi(n)$ are Achilles numbers. An integer is Achilles when it is powerful, but not a perfect power. The problem asks for the number of such $n$ below $10^{18}$.

Mathematical Approach

Write

$$n=\prod_{i=1}^r p_i^{a_i},\qquad a_i\ge 1.$$

1) Achilles condition in exponent language

The number $n$ is powerful exactly when every exponent is at least $2$:

$$a_i\ge 2\qquad \text{for all }i.$$

It is a perfect $k$-th power exactly when all exponents are divisible by the same $k \ge 2$, equivalently when

$$\gcd(a_1,a_2,\dots,a_r)\ge 2.$$

Therefore

$$n\text{ is Achilles}\iff a_i\ge 2\ \forall i\quad\text{and}\quad \gcd(a_1,\dots,a_r)=1.$$

This turns the first half of the problem into a statement about exponent vectors rather than about the decimal value of $n$.

2) What the totient contributes

Euler's product formula gives

$$\varphi(n)=\prod_{i=1}^r p_i^{a_i-1}(p_i-1).$$

Fix a prime $q$. Its exponent inside $\varphi(n)$ is

$$v_q(\varphi(n))=\sum_{i=1}^r v_q(p_i-1)+\sum_{\substack{1\le i\le r\\p_i=q}}(a_i-1).$$

So every exponent of $\varphi(n)$ comes from two sources:

1. factors already forced by the numbers $p_i-1$;

2. the self-contribution $a_i-1$ when the prime $q$ itself appears in $n$.

Thus $\varphi(n)$ is Achilles exactly when all final values $v_q(\varphi(n))$ are at least $2$ and their gcd is $1$.

3) Worked example: $n=500$

Take

$$500=2^2\cdot 5^3.$$

The exponent vector of $n$ is $(2,3)$, so $500$ is powerful and

$$\gcd(2,3)=1.$$

Hence $500$ is Achilles.

Now compute the totient:

$$\varphi(500)=500\left(1-\frac12\right)\left(1-\frac15\right)=200=2^3\cdot 5^2.$$

The exponent vector of $\varphi(500)$ is $(3,2)$, again all exponents are at least $2$ and

$$\gcd(3,2)=1.$$

So $500$ is a strong Achilles number. This is the smallest example, and it is a good mental model for what the recursion is building.

4) Why the search runs from large primes to small primes

The factorization of $p-1$ uses only primes strictly smaller than $p$. Therefore, if the recursion processes prime bases in descending order, then once we move below a prime $q$, no future decision can increase the exponent of $q$ inside $\varphi(n)$.

That means the exponent of $q$ is now final. The code immediately folds that final exponent into a running gcd for $\varphi(n)$ and removes $q$ from the active state. This is the key structural idea that makes the recursion manageable.

5) The meaning of `active_factors`

The array active_factors is a sorted multiset of prime indices. If the index of a prime $q$ appears there exactly $t$ times, then $t$ is the part of the exponent of $q$ in $\varphi(n)$ that has already been created by factors of the form $p-1$ from larger chosen primes.

Those copies are not finalized yet because the recursion may still decide to include $q$ itself as a prime factor of $n$, which would add the extra term $a_q-1$ to $v_q(\varphi(n))$.

Alongside that multiset, the code stores two gcd summaries:

$$g_n=\gcd(\text{chosen exponents in }n),\qquad g_\varphi=\gcd(\text{already finalized exponents in }\varphi(n)).$$

6) Why a brand-new prime must start with exponent $3$

Suppose the recursion introduces a prime $p$ that does not yet occur in active_factors. Then larger primes have contributed no copy of $p$ to $\varphi(n)$. The only guaranteed contribution of $p$ is the self-term $a-1$ coming from $p^{a-1}$ in Euler's formula.

To make $\varphi(n)$ powerful we need

$$a-1\ge 2,$$

hence

$$a\ge 3.$$

That is exactly why the code tries deg = 3,4,5,\dots for a genuinely new prime. This also explains the cube-root pruning later on: every fresh prime costs at least $p^3$.

7) Why an already active prime may start with exponent $2$

Now suppose $p$ already appears in active_factors with multiplicity $t\ge 1$. Then larger primes have already forced a factor $p^t$ into $\varphi(n)$. If we now choose $p^a$ in $n$, the final exponent of $p$ in $\varphi(n)$ becomes

$$t+(a-1).$$

The powerful condition is only

$$t+(a-1)\ge 2.$$

Since $t\ge 1$, the smallest legal exponent is already $a=2$. This is why the code has a separate branch where an already active prime is tried with deg = 2,3,4,\dots.

The example $500=2^2\cdot 5^3$ shows this mechanism clearly: after choosing $5^3$, the factor $5-1=4=2^2$ places two copies of the prime $2$ into active_factors. Then choosing $2^2$ is enough, because the exponent of $2$ in $\varphi(500)$ becomes $2+(2-1)=3$.

8) When a recursive state is counted

A recursive state already represents a complete candidate $n$. It is counted exactly when:

$$g_n=1,$$

every remaining multiplicity inside active_factors is at least $2$, and

$$\gcd\bigl(g_\varphi,\text{all remaining multiplicities}\bigr)=1.$$

The meaning is direct:

1. $g_n=1$ means $n$ is not a perfect power.

2. Remaining multiplicities at least $2$ mean $\varphi(n)$ is powerful.

3. The final gcd equal to $1$ means $\varphi(n)$ is not a perfect power.

9) Why the cube-root bound is correct

If we still want to insert a completely new prime $p$, we must spend at least $p^3$. Therefore a fresh prime can only be used when

$$p^3\le \text{remaining limit}.$$

Equivalently,

$$p\le \left\lfloor\sqrt[3]{\text{remaining limit}}\right\rfloor.$$

This is why the solver computes a cube-root upper bound before opening the “new prime” branch. It removes almost all hopeless candidates very early.

10) Small checkpoints

The implementation contains two explicit tests:

$$\text{count}(10^4)=7,\qquad \text{count}(10^8)=656.$$

The first seven strong Achilles numbers are

$$500,\ 864,\ 1944,\ 2000,\ 2592,\ 3456,\ 5000.$$

These checkpoints are important because the full target $10^{18}$ is far too large for brute force.

How the Code Works

The sieve first generates primes up to a safe bound. Then, for every integer $m$ in the sieve range, it precomputes the prime factorization of $m$ as a multiset of prime indices. This lets the recursion merge in the factors of $p-1$ very cheaply.

The function count_recursive(...) performs three tasks at each state:

1. test whether the current state already forms a strong Achilles number;

2. try primes that already appear in active_factors;

3. try genuinely new primes, subject to the exponent-$3$ rule and the cube-root bound.

Because the search order is canonical and descending, each valid $n$ is generated exactly once.

Complexity Analysis

There is no clean closed formula for the running time, because the search tree depends on the arithmetic structure of the numbers $p-1$. In the worst case the recursion is exponential, but in practice it is kept under control by four strong filters: prime exponents must be large, gcd conditions kill many states, fresh primes satisfy the cube-root bound, and finalized totient exponents are removed from the active state immediately.

The memory usage is modest: recursion depth, a short active multiset, and the sieve/factor tables.

Problem 302 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>

namespace {

using u64 = std::uint64_t;

constexpr int kPrimeSieveN = 6 * 105000;
int prime_count = 0;
int primes[kPrimeSieveN / 4];
unsigned char erat[kPrimeSieveN + 1];

int divisors[kPrimeSieveN + 1][24];
unsigned char divisors_count[kPrimeSieveN + 1];

u64 answer = 0;

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

int cube_root_floor(const u64 n) {
    if (n <= 1ULL) {
        return static_cast<int>(n);
    }
    int r = static_cast<int>(std::pow(static_cast<long double>(n), 1.0L / 3.0L));
    while (static_cast<u64>(r + 1) * static_cast<u64>(r + 1) * static_cast<u64>(r + 1) <= n) ++r;
    while (static_cast<u64>(r) * static_cast<u64>(r) * static_cast<u64>(r) > n) --r;
    return r;
}

void count_recursive(const int max_prime_idx,
                     const int* active_factors,
                     const int active_count,
                     int gcd_exp_n,
                     int gcd_exp_phi,
                     const u64 lim) {
    if (lim == 0ULL) {
        return;
    }

    if (gcd_exp_n == 1) {
        int g = gcd_exp_phi;
        int i = 0;
        while (i < active_count) {
            int j = i + 1;
            while (j < active_count && active_factors[j] == active_factors[i]) ++j;
            const int cnt = j - i;
            if (cnt < 2) {
                g = 0;
                break;
            }
            g = gcd_small(g, cnt);
            i = j;
        }
        if (g == 1) {
            ++answer;
        }
    }

    int min_idx = 0;
    if (active_count >= 2) {
        if (active_factors[active_count - 1] != active_factors[active_count - 2]) {
            min_idx = active_factors[active_count - 1];
        } else {
            for (int k = active_count - 2; k >= 1; --k) {
                if (active_factors[k] != active_factors[k + 1] && active_factors[k] != active_factors[k - 1]) {
                    min_idx = active_factors[k];
                    break;
                }
            }
        }
    } else if (active_count == 1) {
        min_idx = active_factors[0];
    }

    int max_idx = cube_root_floor(lim);
    if (primes[max_prime_idx - 1] > max_idx) {
        int a = 0;
        int b = prime_count - 1;
        while (a < b - 1) {
            const int c = (a + b) / 2;
            if (primes[c] > max_idx) {
                b = c;
            } else {
                a = c;
            }
        }
        max_idx = a;
    }
    if (max_idx > max_prime_idx - 1) {
        max_idx = max_prime_idx - 1;
    }

    int merged[70];
    for (int j = active_count - 1; j >= 0;) {
        int overdeg = 0;
        for (int y = j; y >= 0 && active_factors[y] == active_factors[j]; --y) ++overdeg;

        if (active_factors[j] < max_prime_idx) {
            int k1 = 0;
            int k2 = 0;
            int merged_count = 0;
            const int idx = active_factors[j];
            const int p_minus_one = primes[idx] - 1;
            const int div_cnt = divisors_count[p_minus_one];

            while (k1 < active_count && k2 < div_cnt) {
                if (active_factors[k1] < divisors[p_minus_one][k2]) {
                    merged[merged_count++] = active_factors[k1++];
                } else {
                    merged[merged_count++] = divisors[p_minus_one][k2++];
                }
            }
            while (k2 < div_cnt) merged[merged_count++] = divisors[p_minus_one][k2++];
            while (k1 < active_count) merged[merged_count++] = active_factors[k1++];

            int g2 = gcd_exp_phi;
            while (merged_count > 0 && merged[merged_count - 1] > active_factors[j]) {
                int y = merged_count - 1;
                while (y >= 0 && merged[y] == merged[merged_count - 1]) --y;
                g2 = gcd_small(g2, merged_count - y - 1);
                merged_count = y + 1;
            }
            if (merged_count > 0 && merged[merged_count - 1] == active_factors[j]) {
                merged_count -= overdeg;
            }

            const u64 p = static_cast<u64>(primes[idx]);
            u64 pow = p;
            int deg = 2;
            while (pow <= lim / p) {
                const int new_g1 = gcd_small(gcd_exp_n, deg);
                const int new_g2 = gcd_small(g2, deg - 1 + overdeg);
                count_recursive(idx, merged, merged_count, new_g1, new_g2, lim / pow / p);

                if (pow > (std::numeric_limits<u64>::max() / p)) break;
                pow *= p;
                ++deg;
            }
        }

        j -= overdeg;
        if (overdeg == 1) break;
    }

    int now = 0;
    while (now < active_count && active_factors[now] < min_idx) ++now;

    for (int idx = min_idx; idx <= max_idx; ++idx) {
        if (now < active_count && idx == active_factors[now]) {
            while (now < active_count && active_factors[now] == idx) ++now;
            continue;
        }

        int k1 = 0;
        int k2 = 0;
        int merged_count = 0;
        const int p_minus_one = primes[idx] - 1;
        const int div_cnt = divisors_count[p_minus_one];
        while (k1 < active_count && k2 < div_cnt) {
            if (active_factors[k1] < divisors[p_minus_one][k2]) {
                merged[merged_count++] = active_factors[k1++];
            } else {
                merged[merged_count++] = divisors[p_minus_one][k2++];
            }
        }
        while (k2 < div_cnt) merged[merged_count++] = divisors[p_minus_one][k2++];
        while (k1 < active_count) merged[merged_count++] = active_factors[k1++];

        int g2 = gcd_exp_phi;
        while (merged_count > 0 && merged[merged_count - 1] > idx) {
            int y = merged_count - 1;
            while (y >= 0 && merged[y] == merged[merged_count - 1]) --y;
            g2 = gcd_small(g2, merged_count - y - 1);
            merged_count = y + 1;
        }

        const u64 p = static_cast<u64>(primes[idx]);
        u64 pow = p * p;
        int deg = 3;
        while (pow <= lim / p) {
            const int new_g1 = gcd_small(gcd_exp_n, deg);
            const int new_g2 = gcd_small(g2, deg - 1);
            count_recursive(idx, merged, merged_count, new_g1, new_g2, lim / pow / p);

            if (pow > (std::numeric_limits<u64>::max() / p)) break;
            pow *= p;
            ++deg;
        }
    }
}

void prepare_data() {
    primes[0] = 2;
    primes[1] = 3;
    prime_count = 2;
    erat[1] = 1;

    const int sq = static_cast<int>(std::sqrt(static_cast<long double>(kPrimeSieveN)));
    for (int i = 0; i <= sq; i += 6) {
        if (erat[i + 1] == 0) {
            const int p = i + 1;
            const int step = p * 6;
            for (int j = p * p; j <= kPrimeSieveN; j += step) erat[j] = 1;
            for (int j = p * (p + 4); j <= kPrimeSieveN; j += step) erat[j] = 1;
            primes[prime_count++] = p;
        }
        if (erat[i + 5] == 0) {
            const int p = i + 5;
            const int step = p * 6;
            for (int j = p * p; j <= kPrimeSieveN; j += step) erat[j] = 1;
            for (int j = p * (p + 2); j <= kPrimeSieveN; j += step) erat[j] = 1;
            primes[prime_count++] = p;
        }
    }

    const int start = sq - (sq % 6) + 6;
    for (int i = start; i < kPrimeSieveN; i += 6) {
        if (erat[i + 1] == 0) primes[prime_count++] = i + 1;
        if (erat[i + 5] == 0) primes[prime_count++] = i + 5;
    }

    for (int i = 0; i < prime_count; ++i) {
        for (u64 deg = static_cast<u64>(primes[i]); deg <= static_cast<u64>(kPrimeSieveN);) {
            for (int j = static_cast<int>(deg); j <= kPrimeSieveN; j += static_cast<int>(deg)) {
                divisors[j][divisors_count[j]++] = i;
            }
            if (deg <= static_cast<u64>(kPrimeSieveN / primes[i])) {
                deg *= static_cast<u64>(primes[i]);
            } else {
                break;
            }
        }
    }
}

u64 solve(const u64 limit) {
    answer = 0;
    int root_factors[70] = {0};
    count_recursive(prime_count, root_factors, 0, 0, 0, limit);
    return answer;
}

}  // namespace

int main() {
    prepare_data();
    assert(solve(10'000ULL) == 7ULL);
    assert(solve(100'000'000ULL) == 656ULL);
    std::cout << solve(1'000'000'000'000'000'000ULL) << '\n';
    return 0;
}

Python

import sys

sys.setrecursionlimit(20000)

class Solver:
    def __init__(self):
        self.kPrimeSieveN = 6 * 105000
        self.primes = []
        
        erat = [0] * (self.kPrimeSieveN + 1)
        erat[1] = 1
        
        self.primes.append(2)
        self.primes.append(3)
        
        sq = int(self.kPrimeSieveN ** 0.5)
        for i in range(0, sq + 1, 6):
            if not erat[i + 1]:
                p = i + 1
                step = p * 6
                for j in range(p * p, self.kPrimeSieveN + 1, step): erat[j] = 1
                for j in range(p * (p + 4), self.kPrimeSieveN + 1, step): erat[j] = 1
                self.primes.append(p)
            if not erat[i + 5]:
                p = i + 5
                step = p * 6
                for j in range(p * p, self.kPrimeSieveN + 1, step): erat[j] = 1
                for j in range(p * (p + 2), self.kPrimeSieveN + 1, step): erat[j] = 1
                self.primes.append(p)
                
        start = sq - (sq % 6) + 6
        for i in range(start, self.kPrimeSieveN, 6):
            if not erat[i + 1]: self.primes.append(i + 1)
            if not erat[i + 5]: self.primes.append(i + 5)
            
        self.prime_count = len(self.primes)
        
        self.divisors = [[] for _ in range(self.kPrimeSieveN + 1)]
        for i, p in enumerate(self.primes):
            deg = p
            while deg <= self.kPrimeSieveN:
                for j in range(deg, self.kPrimeSieveN + 1, deg):
                    self.divisors[j].append(i)
                if deg <= self.kPrimeSieveN // p:
                    deg *= p
                else:
                    break
                    
        self.answer = 0

    def gcd_small(self, a, b):
        while b != 0:
            a, b = b, a % b
        return a

    def cube_root_floor(self, n):
        if n <= 1: return int(n)
        r = int(n ** (1.0 / 3.0))
        while (r + 1) * (r + 1) * (r + 1) <= n: r += 1
        while r * r * r > n: r -= 1
        return r

    def count_recursive(self, max_prime_idx, active_factors, active_count, gcd_exp_n, gcd_exp_phi, lim):
        if lim == 0: return

        if gcd_exp_n == 1:
            g = gcd_exp_phi
            i = 0
            while i < active_count:
                j = i + 1
                while j < active_count and active_factors[j] == active_factors[i]:
                    j += 1
                cnt = j - i
                if cnt < 2:
                    g = 0
                    break
                g = self.gcd_small(g, cnt)
                i = j
            if g == 1:
                self.answer += 1

        min_idx = 0
        if active_count >= 2:
            if active_factors[active_count - 1] != active_factors[active_count - 2]:
                min_idx = active_factors[active_count - 1]
            else:
                for k in range(active_count - 2, 0, -1):
                    if active_factors[k] != active_factors[k + 1] and active_factors[k] != active_factors[k - 1]:
                        min_idx = active_factors[k]
                        break
        elif active_count == 1:
            min_idx = active_factors[0]

        max_idx = self.cube_root_floor(lim)
        if self.primes[max_prime_idx - 1] > max_idx:
            a = 0
            b = self.prime_count - 1
            while a < b - 1:
                c = (a + b) // 2
                if self.primes[c] > max_idx:
                    b = c
                else:
                    a = c
            max_idx = a
        if max_idx > max_prime_idx - 1:
            max_idx = max_prime_idx - 1

        merged = [0] * 70
        j = active_count - 1
        while j >= 0:
            overdeg = 0
            y = j
            while y >= 0 and active_factors[y] == active_factors[j]:
                overdeg += 1
                y -= 1

            if active_factors[j] < max_prime_idx:
                k1 = 0
                k2 = 0
                merged_count = 0
                idx = active_factors[j]
                p_minus_one = self.primes[idx] - 1
                divs = self.divisors[p_minus_one]
                div_cnt = len(divs)

                while k1 < active_count and k2 < div_cnt:
                    if active_factors[k1] < divs[k2]:
                        merged[merged_count] = active_factors[k1]
                        k1 += 1
                    else:
                        merged[merged_count] = divs[k2]
                        k2 += 1
                    merged_count += 1
                while k2 < div_cnt:
                    merged[merged_count] = divs[k2]
                    k2 += 1
                    merged_count += 1
                while k1 < active_count:
                    merged[merged_count] = active_factors[k1]
                    k1 += 1
                    merged_count += 1

                g2 = gcd_exp_phi
                while merged_count > 0 and merged[merged_count - 1] > active_factors[j]:
                    y = merged_count - 1
                    while y >= 0 and merged[y] == merged[merged_count - 1]:
                        y -= 1
                    g2 = self.gcd_small(g2, merged_count - y - 1)
                    merged_count = y + 1
                if merged_count > 0 and merged[merged_count - 1] == active_factors[j]:
                    merged_count -= overdeg

                p = self.primes[idx]
                pow_val = p
                deg = 2
                while pow_val <= lim // p:
                    new_g1 = self.gcd_small(gcd_exp_n, deg)
                    new_g2 = self.gcd_small(g2, deg - 1 + overdeg)
                    self.count_recursive(idx, merged, merged_count, new_g1, new_g2, lim // pow_val // p)

                    pow_val *= p
                    deg += 1

            j -= overdeg
            if overdeg == 1: break

        now = 0
        while now < active_count and active_factors[now] < min_idx:
            now += 1

        for idx in range(min_idx, max_idx + 1):
            if now < active_count and idx == active_factors[now]:
                while now < active_count and active_factors[now] == idx:
                    now += 1
                continue

            k1 = 0
            k2 = 0
            merged_count = 0
            p_minus_one = self.primes[idx] - 1
            divs = self.divisors[p_minus_one]
            div_cnt = len(divs)

            while k1 < active_count and k2 < div_cnt:
                if active_factors[k1] < divs[k2]:
                    merged[merged_count] = active_factors[k1]
                    k1 += 1
                else:
                    merged[merged_count] = divs[k2]
                    k2 += 1
                merged_count += 1
            while k2 < div_cnt:
                merged[merged_count] = divs[k2]
                k2 += 1
                merged_count += 1
            while k1 < active_count:
                merged[merged_count] = active_factors[k1]
                k1 += 1
                merged_count += 1

            g2 = gcd_exp_phi
            while merged_count > 0 and merged[merged_count - 1] > idx:
                y = merged_count - 1
                while y >= 0 and merged[y] == merged[merged_count - 1]:
                    y -= 1
                g2 = self.gcd_small(g2, merged_count - y - 1)
                merged_count = y + 1

            p = self.primes[idx]
            pow_val = p * p
            deg = 3
            while pow_val <= lim // p:
                new_g1 = self.gcd_small(gcd_exp_n, deg)
                new_g2 = self.gcd_small(g2, deg - 1)
                self.count_recursive(idx, merged, merged_count, new_g1, new_g2, lim // pow_val // p)

                pow_val *= p
                deg += 1

    def solve_main(self, limit):
        self.answer = 0
        root_factors = [0] * 70
        self.count_recursive(self.prime_count, root_factors, 0, 0, 0, limit)
        return self.answer

def solve(limit=10**18):
    solver = Solver()
    return str(solver.solve_main(limit))

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

Java

public class Euler302 {
    static final int kPrimeSieveN = 6 * 105000;
    static int primeCount = 0;
    static int[] primes = new int[kPrimeSieveN / 4];
    static byte[] erat = new byte[kPrimeSieveN + 1];

    static int[][] divisors = new int[kPrimeSieveN + 1][24];
    static byte[] divisorsCount = new byte[kPrimeSieveN + 1];

    static long answer = 0;

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

    static int cubeRootFloor(long n) {
        if (n <= 1L) {
            return (int) n;
        }
        int r = (int) Math.pow(n, 1.0 / 3.0);
        while ((long) (r + 1) * (r + 1) * (r + 1) <= n)
            ++r;
        while ((long) r * r * r > n)
            --r;
        return r;
    }

    static void countRecursive(int maxPrimeIdx, int[] activeFactors, int activeCount, int gcdExpN, int gcdExpPhi,
            long lim) {
        if (lim == 0L)
            return;

        if (gcdExpN == 1) {
            int g = gcdExpPhi;
            int i = 0;
            while (i < activeCount) {
                int j = i + 1;
                while (j < activeCount && activeFactors[j] == activeFactors[i])
                    ++j;
                int cnt = j - i;
                if (cnt < 2) {
                    g = 0;
                    break;
                }
                g = gcdSmall(g, cnt);
                i = j;
            }
            if (g == 1) {
                ++answer;
            }
        }

        int minIdx = 0;
        if (activeCount >= 2) {
            if (activeFactors[activeCount - 1] != activeFactors[activeCount - 2]) {
                minIdx = activeFactors[activeCount - 1];
            } else {
                for (int k = activeCount - 2; k >= 1; --k) {
                    if (activeFactors[k] != activeFactors[k + 1] && activeFactors[k] != activeFactors[k - 1]) {
                        minIdx = activeFactors[k];
                        break;
                    }
                }
            }
        } else if (activeCount == 1) {
            minIdx = activeFactors[0];
        }

        int maxIdx = cubeRootFloor(lim);
        if (maxPrimeIdx > 0 && primes[maxPrimeIdx - 1] > maxIdx) {
            int a = 0;
            int b = maxPrimeIdx - 1;
            while (a < b - 1) {
                int c = (a + b) / 2;
                if (primes[c] > maxIdx) {
                    b = c;
                } else {
                    a = c;
                }
            }
            maxIdx = a;
        }
        if (maxIdx > maxPrimeIdx - 1) {
            maxIdx = maxPrimeIdx - 1;
        }

        int[] merged = new int[70];
        for (int j = activeCount - 1; j >= 0;) {
            int overdeg = 0;
            for (int y = j; y >= 0 && activeFactors[y] == activeFactors[j]; --y)
                ++overdeg;

            if (activeFactors[j] < maxPrimeIdx) {
                int k1 = 0, k2 = 0, mergedCount = 0;
                int idx = activeFactors[j];
                int pMinusOne = primes[idx] - 1;
                int divCnt = divisorsCount[pMinusOne];

                while (k1 < activeCount && k2 < divCnt) {
                    if (activeFactors[k1] < divisors[pMinusOne][k2]) {
                        merged[mergedCount++] = activeFactors[k1++];
                    } else {
                        merged[mergedCount++] = divisors[pMinusOne][k2++];
                    }
                }
                while (k2 < divCnt)
                    merged[mergedCount++] = divisors[pMinusOne][k2++];
                while (k1 < activeCount)
                    merged[mergedCount++] = activeFactors[k1++];

                int g2 = gcdExpPhi;
                while (mergedCount > 0 && merged[mergedCount - 1] > activeFactors[j]) {
                    int y = mergedCount - 1;
                    while (y >= 0 && merged[y] == merged[mergedCount - 1])
                        --y;
                    g2 = gcdSmall(g2, mergedCount - y - 1);
                    mergedCount = y + 1;
                }
                if (mergedCount > 0 && merged[mergedCount - 1] == activeFactors[j]) {
                    mergedCount -= overdeg;
                }

                long p = primes[idx];
                long pow = p;
                int deg = 2;
                while (pow <= lim / p) {
                    int newG1 = gcdSmall(gcdExpN, deg);
                    int newG2 = gcdSmall(g2, deg - 1 + overdeg);
                    countRecursive(idx, merged, mergedCount, newG1, newG2, lim / pow / p);

                    if (pow > (Long.MAX_VALUE / p))
                        break;
                    pow *= p;
                    ++deg;
                }
            }

            j -= overdeg;
            if (overdeg == 1)
                break;
        }

        int now = 0;
        while (now < activeCount && activeFactors[now] < minIdx)
            ++now;

        for (int idx = minIdx; idx <= maxIdx; ++idx) {
            if (now < activeCount && idx == activeFactors[now]) {
                while (now < activeCount && activeFactors[now] == idx)
                    ++now;
                continue;
            }

            int k1 = 0, k2 = 0, mergedCount = 0;
            int pMinusOne = primes[idx] - 1;
            int divCnt = divisorsCount[pMinusOne];

            while (k1 < activeCount && k2 < divCnt) {
                if (activeFactors[k1] < divisors[pMinusOne][k2]) {
                    merged[mergedCount++] = activeFactors[k1++];
                } else {
                    merged[mergedCount++] = divisors[pMinusOne][k2++];
                }
            }
            while (k2 < divCnt)
                merged[mergedCount++] = divisors[pMinusOne][k2++];
            while (k1 < activeCount)
                merged[mergedCount++] = activeFactors[k1++];

            int g2 = gcdExpPhi;
            while (mergedCount > 0 && merged[mergedCount - 1] > idx) {
                int y = mergedCount - 1;
                while (y >= 0 && merged[y] == merged[mergedCount - 1])
                    --y;
                g2 = gcdSmall(g2, mergedCount - y - 1);
                mergedCount = y + 1;
            }

            long p = primes[idx];
            long pow = p * p;
            int deg = 3;
            while (pow <= lim / p) {
                int newG1 = gcdSmall(gcdExpN, deg);
                int newG2 = gcdSmall(g2, deg - 1);
                countRecursive(idx, merged, mergedCount, newG1, newG2, lim / pow / p);

                if (pow > (Long.MAX_VALUE / p))
                    break;
                pow *= p;
                ++deg;
            }
        }
    }

    static void prepareData() {
        primes[0] = 2;
        primes[1] = 3;
        primeCount = 2;
        erat[1] = 1;

        int sq = (int) Math.sqrt(kPrimeSieveN);
        for (int i = 0; i <= sq; i += 6) {
            if (erat[i + 1] == 0) {
                int p = i + 1;
                int step = p * 6;
                for (int j = p * p; j <= kPrimeSieveN; j += step)
                    erat[j] = 1;
                for (int j = p * (p + 4); j <= kPrimeSieveN; j += step)
                    erat[j] = 1;
                primes[primeCount++] = p;
            }
            if (erat[i + 5] == 0) {
                int p = i + 5;
                int step = p * 6;
                for (int j = p * p; j <= kPrimeSieveN; j += step)
                    erat[j] = 1;
                for (int j = p * (p + 2); j <= kPrimeSieveN; j += step)
                    erat[j] = 1;
                primes[primeCount++] = p;
            }
        }

        int start = sq - (sq % 6) + 6;
        for (int i = start; i < kPrimeSieveN; i += 6) {
            if (erat[i + 1] == 0)
                primes[primeCount++] = i + 1;
            if (erat[i + 5] == 0)
                primes[primeCount++] = i + 5;
        }

        for (int i = 0; i < primeCount; ++i) {
            for (long deg = primes[i]; deg <= kPrimeSieveN;) {
                for (int j = (int) deg; j <= kPrimeSieveN; j += (int) deg) {
                    divisors[j][divisorsCount[j]++] = i;
                }
                if (deg <= kPrimeSieveN / primes[i]) {
                    deg *= primes[i];
                } else {
                    break;
                }
            }
        }
    }

    public static String solve() {
        prepareData();
        answer = 0;
        int[] rootFactors = new int[70];
        countRecursive(primeCount, rootFactors, 0, 0, 0, 1000000000000000000L);
        return String.valueOf(answer);
    }

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

Problem 302: Strong Achilles Numbers

Project Euler Problem 302 Solution

Detailed mathematical approach

Problem Summary

Mathematical Approach

1) Achilles condition in exponent language

2) What the totient contributes

3) Worked example: \(n=500\)

4) Why the search runs from large primes to small primes

5) The meaning of `active_factors`

6) Why a brand-new prime must start with exponent \(3\)

7) Why an already active prime may start with exponent \(2\)

8) When a recursive state is counted

9) Why the cube-root bound is correct

10) Small checkpoints

How the Code Works

Complexity Analysis

Further Reading

Problem 302 source code

C++

Python

Java

Problem 302: Strong Achilles Numbers

Project Euler Problem 302 Solution

Detailed mathematical approach

Problem Summary

Mathematical Approach

1) Achilles condition in exponent language

2) What the totient contributes

3) Worked example: \(n=500\)

4) Why the search runs from large primes to small primes

5) The meaning of active_factors

6) Why a brand-new prime must start with exponent \(3\)

7) Why an already active prime may start with exponent \(2\)

8) When a recursive state is counted

9) Why the cube-root bound is correct

10) Small checkpoints

How the Code Works

Complexity Analysis

Further Reading

Problem 302 source code

C++

Python

Java

5) The meaning of `active_factors`