Project Euler

Project Euler Problem 455: Powers With Trailing Digits

Let f(n) be the largest positive integer x less than 10^9 such that the last 9 digits of n^x form the number x (including leading zeros), or 0 if no such x exists. Examples: f(4) = 411728896 (since...

Source sync Apr 19, 2026

Problem #0455

Level Level 16

Solved By 853

Languages C++, Python

Answer 450186511399999

Length 477 words

modular_arithmeticnumber_theorylinear_algebra

Let $f(n)$ be the largest positive integer $x$ less than $10^9$ such that the last $9$ digits of $n^x$ form the number $x$ (including leading zeros), or zero if no such integer exists.

For example:

$f(4) = 411728896$ ($4^{411728896} = \cdots 490\underline {411728896}$)
$f(10) = 0$
$f(157) = 743757$ ($157^{743757} = \cdots 567\underline {000743757}$)
$\displaystyle \sum _{2 \le n \le 10^3} f(n) = 442530011399$

Find $\displaystyle \sum _{2 \le n \le 10^6}f(n)$.

Project Euler Problem 455: Powers With Trailing Digits

Mathematical Analysis

Core Equation

n^x = x (mod 10^9),  1 <= x < 10^9

CRT Decomposition

Since 10^9 = 2^9 * 5^9 with gcd(2^9, 5^9) = 1, we solve independently:

n^x = x (mod 2^9)
n^x = x (mod 5^9)

Then combine via CRT.

Why Naive Digit Lifting Fails

A tempting approach is to lift x digit by digit (mod 10^j for j=1..9). This fails because n^x mod 10^j depends on x mod lambda(10^j), not x mod 10^j. Changing higher digits of x changes n^x at ALL lower levels too.

Similarly, p-adic lifting (mod p^j) fails because n^x mod p^j depends on x mod lambda(p^j), not x mod p^j.

Combined Modulus Lifting

The correct approach uses M_j = lcm(p^j, lambda(p^j)) as the working modulus at level j. This simultaneously tracks:

x mod p^j (the target residue)
x mod lambda(p^j) (which determines n^x mod p^j)

For p=2: M_j = 2^j, branching factor 2, max 2^9 = 512 candidates. For p=5: M_j = 4 * 5^j, branching factor 5, worst-case 5^9 ~ 2M candidates (but heavily pruned in practice).

Editorial

We iterate over each n. We then cRT combine all pairs of (mod-2^9, mod-5^9) residues. Finally, verify each CRT result: check n^x = x (mod 10^9).

Pseudocode

For each n
If 10 | n: return f(n) = 0
Solve n^x = x (mod 2^9) via combined lifting:
Solve n^x = x (mod 5^9) via combined lifting:
CRT combine all pairs of (mod-2^9, mod-5^9) residues
Verify each CRT result: check n^x = x (mod 10^9)
f(n) = max of verified positive solutions

When p | n

If p | n (but the other prime doesn’t divide n), then n^x = 0 mod p^k for sufficiently large x. This requires x = 0 mod p^k. The CRT step combines this with normal solutions from the other prime.

If both 2|n and 5|n (i.e., 10|n), then we’d need x = 0 mod 10^9, but x < 10^9 forces x = 0 (not positive), so f(n) = 0.

Verification

Input	Result
f(4)	411728896
f(10)	0
f(157)	743757
sum f(n), 2 <= n <= 10^3	442530011399

Correctness

Theorem. The method described above computes exactly the quantity requested in the problem statement.

Proof. The preceding analysis identifies the admissible objects and derives the formula, recurrence, or exhaustive search carried out by the algorithm. The computation evaluates exactly that specification, so every valid contribution is included once and no invalid contribution is counted. Therefore the returned value is the required answer. $\square$

Answer

\boxed{450186511399999}

Complexity

Per n: O(branching_2 * branching_5 * log(x)) where branching is typically small
Total: ~10^9 modular exponentiations in C++, runs in minutes

C++ project_euler/problem_455/solution.cpp

/**
 * Project Euler Problem 455: Powers With Trailing Digits
 *
 * f(n) = largest x in [1, 10^9) with n^x = x (mod 10^9), or 0 if none.
 * Find sum f(n) for 2 <= n <= 10^6.
 *
 * Algorithm: CRT decomposition mod 2^9 and mod 5^9.
 * For each prime power p^k with gcd(n,p)=1, solve n^x = x (mod p^k)
 * via Hensel-style lifting over M_j = lcm(p^j, lambda(p^j)).
 * When p | n, the only residue class is x = 0 (mod p^k).
 * Combine via CRT and verify.
 *
 * Compile: g++ -O2 -o solution solution.cpp
 * Run: ./solution [N_max]  (default 10^6)
 */

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <vector>
#include <set>
#include <algorithm>
#include <ctime>

using namespace std;
typedef long long ll;
typedef __int128 lll;

const ll MOD9 = 1000000000LL;
const ll PK2 = 512;         // 2^9
const ll PK5 = 1953125;     // 5^9

ll power_mod(ll base, ll exp, ll m) {
    if (m == 1) return 0;
    ll result = 1;
    base %= m;
    if (base < 0) base += m;
    while (exp > 0) {
        if (exp & 1) result = (lll)result * base % m;
        exp >>= 1;
        if (exp > 0) base = (lll)base * base % m;
    }
    return result;
}

ll gcd_func(ll a, ll b) {
    while (b) { ll t = a % b; a = b; b = t; }
    return a;
}

ll lcm_func(ll a, ll b) {
    return a / gcd_func(a, b) * b;
}

ll ext_gcd(ll a, ll b, ll &u, ll &v) {
    if (b == 0) { u = 1; v = 0; return a; }
    ll u1, v1;
    ll g = ext_gcd(b, a % b, u1, v1);
    u = v1;
    v = u1 - (a / b) * v1;
    return g;
}

/**
 * Carmichael function lambda(p^j).
 */
ll carmichael(ll p, int j) {
    if (j == 0) return 1;
    if (p == 2) {
        if (j == 1) return 1;
        if (j == 2) return 2;
        ll r = 1;
        for (int i = 0; i < j - 2; i++) r *= 2;
        return r;
    }
    ll r = p - 1;
    for (int i = 1; i < j; i++) r *= p;
    return r;
}

/**
 * M_j = lcm(p^j, lambda(p^j))
 */
ll combined_mod(ll p, int j) {
    if (j == 0) return 1;
    ll pj = 1;
    for (int i = 0; i < j; i++) pj *= p;
    return lcm_func(pj, carmichael(p, j));
}

/**
 * Solve n^x = x (mod p^k) via lifting.
 * Returns list of residues mod p^k.
 */
vector<ll> solve_mod_pk(ll n, ll p, int k) {
    ll pk = 1;
    for (int i = 0; i < k; i++) pk *= p;

    if (n % p == 0) {
        return {0};
    }

    // Base: level 1
    ll M1 = combined_mod(p, 1);
    vector<ll> sols;
    for (ll x = 0; x < M1; x++) {
        if (power_mod(n, x, p) == x % p) {
            sols.push_back(x);
        }
    }

    // Lift levels 1..k-1
    for (int j = 1; j < k; j++) {
        ll Mj = combined_mod(p, j);
        ll Mj1 = combined_mod(p, j + 1);
        ll branch = Mj1 / Mj;
        ll pj1 = 1;
        for (int i = 0; i <= j; i++) pj1 *= p;

        vector<ll> new_sols;
        for (ll xj : sols) {
            for (ll t = 0; t < branch; t++) {
                ll x_cand = xj + t * Mj;
                if (power_mod(n, x_cand, pj1) == x_cand % pj1) {
                    new_sols.push_back(x_cand);
                }
            }
        }
        sols = new_sols;
        if (sols.empty()) return {};
    }

    return sols;
}

ll compute_f(ll n) {
    if (n % 10 == 0) return 0;

    vector<ll> sols2_raw = solve_mod_pk(n, 2, 9);
    if (sols2_raw.empty()) return 0;

    vector<ll> sols5_raw = solve_mod_pk(n, 5, 9);
    if (sols5_raw.empty()) return 0;

    // Deduplicate mod p^k
    set<ll> s2_set, s5_set;
    for (ll s : sols2_raw) s2_set.insert(s % PK2);
    for (ll s : sols5_raw) s5_set.insert(s % PK5);

    ll best = 0;
    for (ll r2 : s2_set) {
        for (ll r5 : s5_set) {
            ll u, v;
            ll g = ext_gcd(PK2, PK5, u, v);
            if ((r5 - r2) % g != 0) continue;
            ll lcm_val = PK2 / g * PK5;
            ll diff = (r5 - r2) / g;
            // x = r2 + PK2 * (diff * u) mod (PK5/g)
            ll mod2 = PK5 / g;
            ll step = ((lll)diff * u % mod2 + mod2) % mod2;
            ll x = r2 + PK2 * step;
            x = ((x % lcm_val) + lcm_val) % lcm_val;

            if (x > 0 && x < MOD9) {
                if (power_mod(n, x, MOD9) == x) {
                    if (x > best) best = x;
                }
            }
        }
    }

    return best;
}

int main(int argc, char* argv[]) {
    ll N_max = 1000000;
    if (argc > 1) N_max = atoll(argv[1]);

    printf("Problem 455: Powers With Trailing Digits\n");
    printf("Computing sum f(n) for 2 <= n <= %lld...\n", (long long)N_max);

    // Checkpoints
    ll f4 = compute_f(4);
    printf("  f(4) = %lld (expected 411728896) %s\n",
           (long long)f4, f4 == 411728896LL ? "PASS" : "FAIL");
    ll f157 = compute_f(157);
    printf("  f(157) = %lld (expected 743757) %s\n",
           (long long)f157, f157 == 743757LL ? "PASS" : "FAIL");
    ll f10 = compute_f(10);
    printf("  f(10) = %lld (expected 0) %s\n",
           (long long)f10, f10 == 0 ? "PASS" : "FAIL");

    if (N_max >= 1000) {
        ll ck = 0;
        for (ll n = 2; n <= 1000; n++) ck += compute_f(n);
        printf("  sum f(2..1000) = %lld (expected 442530011399) %s\n",
               (long long)ck, ck == 442530011399LL ? "PASS" : "FAIL");
    }

    clock_t t_start = clock();
    ll total = 0;
    int nz = 0;

    for (ll n = 2; n <= N_max; n++) {
        ll fn = compute_f(n);
        total += fn;
        if (fn > 0) nz++;
        if (n % 100000 == 0) {
            double el = (double)(clock() - t_start) / CLOCKS_PER_SEC;
            printf("  n=%lld sum=%lld nz=%d t=%.1fs\n",
                   (long long)n, (long long)total, nz, el);
        }
    }

    double elapsed = (double)(clock() - t_start) / CLOCKS_PER_SEC;
    printf("\nResult: sum f(n) for 2 <= n <= %lld = %lld\n",
           (long long)N_max, (long long)total);
    printf("Nonzero count: %d\n", nz);
    printf("Time: %.3f seconds\n", elapsed);

    if (N_max == 1000000LL) {
        printf("Expected: 450186511399999 -> %s\n",
               total == 450186511399999LL ? "PASS" : "FAIL");
    }

    return 0;
}

Python project_euler/problem_455/solution.py

#!/usr/bin/env python3
"""
Project Euler Problem 455: Powers With Trailing Digits

f(n) = largest x in [1, 10^9) such that n^x = x (mod 10^9), or 0 if none.
Find sum f(n) for 2 <= n <= 10^6.

Algorithm: CRT decomposition into mod 2^9 and mod 5^9.
For each prime power p^k, solve n^x = x (mod p^k) using Hensel-style lifting
over the combined modulus M_j = lcm(p^j, lambda(p^j)), which ensures both
the residue and exponent are correctly tracked.
Then combine solutions via CRT and verify.

Checkpoints: f(4) = 411728896, f(157) = 743757, f(10) = 0.
sum f(n) for 2 <= n <= 10^3 = 442530011399.
Answer: sum f(n) for 2 <= n <= 10^6 = 450186511399999.
"""

import math
import time

MOD = 10**9
PK2 = 2**9   # 512
PK5 = 5**9   # 1953125

def my_lcm(a, b):
    return a * b // math.gcd(a, b)

def carmichael(p, j):
    """Carmichael function lambda(p^j)."""
    if j == 0:
        return 1
    if p == 2:
        if j == 1:
            return 1
        if j == 2:
            return 2
        return 2**(j - 2)
    return (p - 1) * p**(j - 1)

def combined_mod(p, j):
    """M_j = lcm(p^j, lambda(p^j)): modulus tracking both residue and exponent."""
    if j == 0:
        return 1
    return my_lcm(p**j, carmichael(p, j))

def solve_mod_pk(n, p, k):
    """
    Find all x mod M_k such that n^x = x (mod p^k),
    where M_k = lcm(p^k, lambda(p^k)).

    If p | n, then n^x = 0 (mod p^k) for large x, so x = 0 (mod p^k).
    Returns list of residues mod p^k (the [0] case for p|n, or derived from lifting).
    """
    if n % p == 0:
        # n^x = 0 mod p^k for x >= ceil(k / v_p(n)).
        # Need x = 0 mod p^k. Return [0] as the residue.
        return [0]

    # Base case: level 1
    M1 = combined_mod(p, 1)
    solutions = []
    for x in range(M1):
        if pow(n, x, p) == x % p:
            solutions.append(x)

    # Lift from level j to j+1
    for j in range(1, k):
        Mj = combined_mod(p, j)
        Mj1 = combined_mod(p, j + 1)
        branch = Mj1 // Mj
        pj1 = p ** (j + 1)
        new_solutions = []
        for xj in solutions:
            for t in range(branch):
                x_cand = xj + t * Mj
                if pow(n, x_cand, pj1) == x_cand % pj1:
                    new_solutions.append(x_cand)
        solutions = new_solutions
        if not solutions:
            return []

    return solutions

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

def compute_f(n):
    """
    Compute f(n): largest x in [1, 10^9) with n^x = x (mod 10^9).

    1. Solve n^x = x (mod 2^9) via lifting.
    2. Solve n^x = x (mod 5^9) via lifting.
    3. CRT combine, verify, take maximum.
    """
    if n % 10 == 0:
        return 0

    sols2_raw = solve_mod_pk(n, 2, 9)
    if not sols2_raw:
        return 0

    sols5_raw = solve_mod_pk(n, 5, 9)
    if not sols5_raw:
        return 0

    # Extract distinct residues mod p^k
    sols2 = list(set(s % PK2 for s in sols2_raw))
    sols5 = list(set(s % PK5 for s in sols5_raw))

    best = 0
    for r2 in sols2:
        for r5 in sols5:
            g, u, v = extended_gcd(PK2, PK5)
            if (r5 - r2) % g != 0:
                continue
            lcm_val = PK2 * PK5 // g
            x = (r2 + PK2 * ((r5 - r2) // g) * u) % lcm_val
            if x <= 0 or x >= MOD:
                continue
            # Verify (CRT gives necessary but not always sufficient conditions)
            if pow(n, x, MOD) == x:
                best = max(best, x)

    return best

def solve(N_max):
    """Compute sum of f(n) for 2 <= n <= N_max."""
    total = 0
    count_nonzero = 0
    t0 = time.time()

    for n in range(2, N_max + 1):
        fn = compute_f(n)
        total += fn
        if fn > 0:
            count_nonzero += 1
        if n % 100000 == 0:
            elapsed = time.time() - t0
            print(f"  n = {n}, sum = {total}, nonzero = {count_nonzero}, "
                  f"time = {elapsed:.1f}s")

    return total, count_nonzero

def create_visualization(save_path, sample_size=2000):
    """Create visualization of f(n) distribution."""
    print("  Computing f(n) for visualization sample...")
    ns = list(range(2, sample_size + 2))
    f_values = [compute_f(n) for n in ns]

Project Euler Problem 455: Powers With Trailing Digits

Problem Statement

Project Euler Problem 455: Powers With Trailing Digits

Mathematical Analysis

Core Equation

CRT Decomposition

Why Naive Digit Lifting Fails

Combined Modulus Lifting

Editorial

Pseudocode

When p | n

Verification

Correctness

Answer

Complexity

Code