Problem 533: Minimum Values of the Carmichael Function

Project Euler Problem 533 Solution

EulerSolve provides an optimized solution for Project Euler Problem 533, Minimum Values of the Carmichael Function, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For Carmichael's function $\lambda(m)$, define $$\mathcal L(n)=\max\{m:\lambda(m) \lt n\}+1.$$ The goal is to find the last nine digits of $\mathcal L(2\cdot 10^7)$. A direct search over $m$ is hopelessly large, so the implementation works in the reverse direction: it searches over promising values of a target $t$ with $1\le t\lt n$, and for each such $t$ constructs the largest integer whose Carmichael value divides $t$. Mathematical Approach For a fixed integer $t\ge 1$, let $$M(t)=\max\{m:\lambda(m)\mid t\}.$$ Then the original problem becomes $$\mathcal L(n)-1=\max_{1\le t\lt n} M(t).$$ This reformulation is valid because every integer $m$ with $\lambda(m)\lt n$ appears when we choose $t=\lambda(m)$, and conversely every $m$ counted by $M(t)$ satisfies $\lambda(m)\mid t\lt n$, hence $\lambda(m)\lt n$. Step 1: Carmichael Values on Prime Powers The construction starts from the standard formulas for prime powers....

Detailed mathematical approach

Problem Summary

For Carmichael's function $\lambda(m)$, define

$$\mathcal L(n)=\max\{m:\lambda(m) \lt n\}+1.$$

The goal is to find the last nine digits of $\mathcal L(2\cdot 10^7)$. A direct search over $m$ is hopelessly large, so the implementation works in the reverse direction: it searches over promising values of a target $t$ with $1\le t\lt n$, and for each such $t$ constructs the largest integer whose Carmichael value divides $t$.

Mathematical Approach

For a fixed integer $t\ge 1$, let

$$M(t)=\max\{m:\lambda(m)\mid t\}.$$

Then the original problem becomes

$$\mathcal L(n)-1=\max_{1\le t\lt n} M(t).$$

This reformulation is valid because every integer $m$ with $\lambda(m)\lt n$ appears when we choose $t=\lambda(m)$, and conversely every $m$ counted by $M(t)$ satisfies $\lambda(m)\mid t\lt n$, hence $\lambda(m)\lt n$.

Step 1: Carmichael Values on Prime Powers

The construction starts from the standard formulas for prime powers. For an odd prime $p$,

$$\lambda(p^a)=p^{a-1}(p-1).$$

For powers of $2$,

$$\lambda(2)=1,\qquad \lambda(4)=2,\qquad \lambda(2^a)=2^{a-2}\quad (a\ge 3).$$

$$m=\prod_i p_i^{a_i},$$

then Carmichael's function combines these local contributions by least common multiple:

$$\lambda(m)=\operatorname{lcm}\bigl(\lambda(p_1^{a_1}),\lambda(p_2^{a_2}),\dots\bigr).$$

So the entire problem reduces to understanding which prime powers may occur when $\lambda(m)$ is forced to divide a chosen value $t$.

Step 2: Maximal Exponent of an Odd Prime for Fixed $t$

Assume $p$ is odd and $p^a\parallel m$. If $\lambda(m)\mid t$, then in particular

$$\lambda(p^a)=p^{a-1}(p-1)\mid t.$$

Because $p$ and $p-1$ are coprime, this forces the two independent divisibility conditions

$$p-1\mid t,\qquad p^{a-1}\mid t.$$

Hence the exponent $a$ cannot exceed

$$a\le v_p(t)+1,$$

where $v_p(t)$ is the exponent of $p$ in the factorization of $t$. Therefore an odd prime $p$ can appear at all only when $p-1\mid t$, and when it does appear its largest admissible power is

$$p^{v_p(t)+1}.$$

This already explains a key feature of the search: divisors of $t$ of the form $p-1$ generate possible prime factors $p$ of the maximizing integer.

Step 3: The Special Role of $2$

The prime $2$ is different because its Carmichael values do not follow the odd-prime formula. If $t$ is odd, then $\lambda(2)=1\mid t$ but $\lambda(4)=2\nmid t$, so the largest allowed power of $2$ is just $2^1$.

If $v_2(t)=r\ge 1$, then from $\lambda(2^a)=2^{a-2}$ for $a\ge 3$ we obtain

$$2^{a-2}\mid t\iff a\le r+2.$$

So the largest admissible exponent of $2$ is

$$a_2(t)= \begin{cases} 1,& v_2(t)=0,\\ v_2(t)+2,& v_2(t)\ge 1. \end{cases}$$

This is why the implementation always inserts a power of $2$ first and treats it separately from the odd primes.

Step 4: Closed Formula for the Largest Integer with $\lambda(m)\mid t$

Combining the odd-prime rule with the special $2$-power rule gives

$$M(t)=2^{a_2(t)}\prod_{\substack{p\text{ odd prime}\\p-1\mid t}} p^{v_p(t)+1}.$$

The code evaluates the same formula in a divisor-based form. Since the condition $p-1\mid t$ is equivalent to saying that $d=p-1$ is a divisor of $t$, we can rewrite it as

$$M(t)=2^{a_2(t)}\prod_{\substack{d\mid t\\d+1\text{ prime}\\d+1\gt 2}} (d+1)^{v_{d+1}(t)+1}.$$

This is exactly why divisor generation is central: each divisor $d$ gives one potential prime $d+1$, and if that number is prime it contributes a whole prime power to the product.

Step 5: Why the Search Focuses on Divisor-Rich $t$

The formula for $M(t)$ shows two ways to make it large.

First, $t$ should have many divisors, because every divisor $d$ is a chance that $d+1$ is prime. Second, $t$ should contain small prime factors with reasonably large exponents, because the term $v_p(t)+1$ can raise already-admissible primes to higher powers.

For that reason, the implementation does not scan every integer below $n$. Instead, it explores a structured family of candidates built from small primes with nonincreasing exponents. This is the same pattern that appears in highly composite style searches: it strongly favors numbers with rich divisor structure, which are exactly the values of $t$ most likely to maximize $M(t)$.

Step 6: Worked Example $\mathcal L(6)=241$

To illustrate the method, take $n=6$. Then we must maximize $M(t)$ over $1\le t\lt 6$.

For $t=1$, only the factor $2$ is possible, so $M(1)=2$.

For $t=2$, we have $a_2(2)=3$, and the divisor $2$ yields the prime $3$. Therefore

$$M(2)=2^3\cdot 3=24.$$

For $t=3$, there is again no odd prime with $p-1\mid 3$ except $p=2$, so $M(3)=2$.

For $t=4$, we have $a_2(4)=4$, and the divisors $2$ and $4$ yield the primes $3$ and $5$. Hence

$$M(4)=2^4\cdot 3\cdot 5=240.$$

For $t=5$, no new odd prime appears, so $M(5)=2$. The maximum is therefore $240$, which means

$$\mathcal L(6)=240+1=241.$$

This is the small checkpoint verified by the implementations before they tackle the full input.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. They first build a primality sieve up to $n$, because every tested number has the form $d+1$ with $d\lt t\lt n$.

Next they generate candidate values of $t$ by depth-first search over products of small primes with nonincreasing exponents and with $t\lt n$. During this search they keep the factorization of the current candidate, so every divisor of $t$ can later be reconstructed efficiently.

For one candidate $t$, the implementation enumerates all divisors of $t$, starts with the special $2$-power contribution, and then checks each number $d+1$. If $d+1$ is prime, the corresponding prime power $(d+1)^{v_{d+1}(t)+1}$ is multiplied into the current value of $M(t)$.

The true integers involved are enormous, so the implementations store two parallel views of each candidate: a logarithmic score for accurate comparison of magnitudes, and the value modulo $10^9$ for the requested last nine digits. After the best candidate has been found, they return

$$\bigl(\max M(t)+1\bigr)\bmod 10^9,$$

with leading zeros preserved in the final output format when necessary.

Complexity Analysis

Let $n$ be the input bound, and let $S$ be the structured set of $t$-candidates explored by the search. Building the primality sieve up to $n$ costs $O(n\log\log n)$ time and $O(n)$ memory.

For one candidate $t$, if $\tau(t)$ denotes its number of divisors, then generating all divisors and scanning them costs $O(\tau(t))$ time and $O(\tau(t))$ temporary space. The remaining work per divisor is small: a primality lookup, a short factor-exponent query, a modular exponentiation with a tiny exponent, and one logarithmic accumulation.

Therefore the overall running time is

$$O\left(n\log\log n+\sum_{t\in S}\tau(t)\right),$$

and the memory usage is

$$O\left(n+\max_{t\in S}\tau(t)\right).$$

In practice the dominant cost is not modular arithmetic but the number of divisor-rich candidates explored by the depth-first search.

Footnotes and References

Problem page: https://projecteuler.net/problem=533
Carmichael function: Wikipedia — Carmichael function
Prime factorization and valuations: Wikipedia — $p$-adic valuation
Divisor function: Wikipedia — Divisor function
Highly composite number: Wikipedia — Highly composite number

Problem 533 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u64 kLast = 1'000'000'000ULL;  // last 9 digits

u64 pow_mod(u64 a, int e, u64 mod) {
    u64 r = 1 % mod;
    a %= mod;
    while (e > 0) {
        if (e & 1) {
            r = static_cast<u128>(r) * a % mod;
        }
        a = static_cast<u128>(a) * a % mod;
        e >>= 1;
    }
    return r;
}

std::vector<bool> sieve_is_prime(const int n) {
    std::vector<bool> is_prime(static_cast<std::size_t>(n + 1), true);
    if (n >= 0) {
        is_prime[0] = false;
    }
    if (n >= 1) {
        is_prime[1] = false;
    }
    for (int p = 2; (u64)p * p <= static_cast<u64>(n); ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        for (int m = p * p; m <= n; m += p) {
            is_prime[static_cast<std::size_t>(m)] = false;
        }
    }
    return is_prime;
}

struct BestState {
    long double log_value = -1e300L;
    u64 mod_value = 0;
    u64 L = 0;
};

u64 exponent_in_factorization(const std::vector<std::pair<u64, int>>& fac, const u64 p) {
    for (const auto& [q, e] : fac) {
        if (q == p) {
            return static_cast<u64>(e);
        }
    }
    return 0;
}

void evaluate_L(const u64 M, const std::vector<bool>& is_prime, const u64 L,
                const std::vector<std::pair<u64, int>>& fac, BestState& best) {
    // Build full divisor list of L.
    std::vector<u64> divisors;
    divisors.reserve(4096);
    divisors.push_back(1);
    for (const auto& [p, e] : fac) {
        const std::size_t cur = divisors.size();
        u64 pe = 1;
        for (int i = 1; i <= e; ++i) {
            pe *= p;
            for (std::size_t j = 0; j < cur; ++j) {
                divisors.push_back(divisors[j] * pe);
            }
        }
    }

    const u64 v2 = exponent_in_factorization(fac, 2);
    const int exp2 = (v2 == 0) ? 1 : static_cast<int>(v2 + 2);

    long double logN = static_cast<long double>(exp2) * std::log(2.0L);
    u64 modN = pow_mod(2, exp2, kLast);

    for (const u64 d : divisors) {
        const u64 p = d + 1;
        if (p == 2) {
            continue;
        }
        if (p > M + 1) {
            continue;
        }
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        const int exp = static_cast<int>(exponent_in_factorization(fac, p)) + 1;
        logN += static_cast<long double>(exp) * std::log(static_cast<long double>(p));
        modN = static_cast<u128>(modN) * pow_mod(p, exp, kLast) % kLast;
    }

    if (logN > best.log_value) {
        best.log_value = logN;
        best.mod_value = modN;
        best.L = L;
    }
}

void dfs_generate(const std::vector<u64>& primes, const u64 M, const std::vector<bool>& is_prime, int idx, int max_exp,
                  u64 cur, std::vector<std::pair<u64, int>>& fac, BestState& best) {
    evaluate_L(M, is_prime, cur, fac, best);
    if (idx >= static_cast<int>(primes.size())) {
        return;
    }

    const u64 p = primes[static_cast<std::size_t>(idx)];
    u64 val = cur;
    for (int e = 1; e <= max_exp; ++e) {
        if (val > M / p) {
            break;
        }
        val *= p;
        fac.push_back({p, e});
        dfs_generate(primes, M, is_prime, idx + 1, e, val, fac, best);
        fac.pop_back();
    }
}

u64 solve_last9(const u64 bound_minus_1) {
    const u64 M = bound_minus_1;
    const auto is_prime = sieve_is_prime(static_cast<int>(M + 1));
    // Enough primes to generate all relevant "highly composite"-style candidates under 2e7.
    const std::vector<u64> primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43};

    BestState best;
    std::vector<std::pair<u64, int>> fac;
    dfs_generate(primes, M, is_prime, 0, 60, 1, fac, best);

    // L(n) = max{ m : lambda(m) < n } + 1, and last 9 digits are requested.
    return (best.mod_value + 1) % kLast;
}

bool run_checkpoints() {
    if (solve_last9(6 - 1) != 241ULL) {
        std::cerr << "Checkpoint failed: L(6)\n";
        return false;
    }
    if (solve_last9(100 - 1) != 174'525'281ULL) {
        std::cerr << "Checkpoint failed: L(100) last 9 digits\n";
        return false;
    }
    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }
    const u64 last9 = solve_last9(20'000'000ULL - 1);
    std::cout << std::setw(9) << std::setfill('0') << last9 << '\n';
    return 0;
}

Python

import math

kLast = 1000000000

def sieve_is_prime(n):
    is_prime = bytearray([1] * (n + 1))
    if n >= 0: is_prime[0] = 0
    if n >= 1: is_prime[1] = 0
    
    p = 2
    while p * p <= n:
        if is_prime[p]:
            for m in range(p * p, n + 1, p):
                is_prime[m] = 0
        p += 1
    return is_prime

best_log_value = -1e300
best_mod_value = 0
best_L = 0

def evaluate_L(M, is_prime, L, fac):
    global best_log_value, best_mod_value, best_L
    
    divisors = [1]
    for p, e in fac:
        cur_len = len(divisors)
        pe = 1
        for i in range(1, e + 1):
            pe *= p
            for j in range(cur_len):
                divisors.append(divisors[j] * pe)
                
    v2 = 0
    for q, e in fac:
        if q == 2:
            v2 = e
            break
            
    exp2 = 1 if v2 == 0 else v2 + 2
    
    logN = exp2 * math.log(2.0)
    modN = pow(2, exp2, kLast)
    
    for d in divisors:
        p = d + 1
        if p == 2: continue
        if p > M + 1: continue
        if not is_prime[p]: continue
        
        exp = 1
        for q, e in fac:
            if q == p:
                exp = e + 1
                break
                
        logN += exp * math.log(p)
        modN = (modN * pow(p, exp, kLast)) % kLast
        
    if logN > best_log_value:
        best_log_value = logN
        best_mod_value = modN
        best_L = L

def dfs_generate(primes, M, is_prime, idx, max_exp, cur, fac):
    evaluate_L(M, is_prime, cur, fac)
    if idx >= len(primes):
        return
        
    p = primes[idx]
    val = cur
    for e in range(1, max_exp + 1):
        if val > M // p:
            break
        val *= p
        fac.append((p, e))
        dfs_generate(primes, M, is_prime, idx + 1, e, val, fac)
        fac.pop()

def solve_last9(bound_minus_1):
    global best_log_value, best_mod_value, best_L
    best_log_value = -1e300
    best_mod_value = 0
    best_L = 0
    
    M = bound_minus_1
    is_prime = sieve_is_prime(M + 1)
    primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43]
    
    fac = []
    dfs_generate(primes, M, is_prime, 0, 60, 1, fac)
    
    return (best_mod_value + 1) % kLast

def solve():
    last9 = solve_last9(20000000 - 1)
    return f"{last9:09d}"

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

Java

import java.util.ArrayList;
import java.util.List;

public class Euler533 {

    private static final long kLast = 1000000000L;

    private static long powMod(long a, int e, long mod) {
        long r = 1 % mod;
        a %= mod;
        while (e > 0) {
            if ((e & 1) != 0) {
                r = (r * a) % mod;
            }
            a = (a * a) % mod;
            e >>= 1;
        }
        return r;
    }

    private static byte[] sieveIsPrime(int n) {
        byte[] isPrime = new byte[n + 1];
        for (int i = 0; i <= n; i++)
            isPrime[i] = 1;
        if (n >= 0)
            isPrime[0] = 0;
        if (n >= 1)
            isPrime[1] = 0;
        for (int p = 2; (long) p * p <= n; ++p) {
            if (isPrime[p] == 1) {
                for (int m = p * p; m <= n; m += p) {
                    isPrime[m] = 0;
                }
            }
        }
        return isPrime;
    }

    static class BestState {
        double logValue = -1e300;
        long modValue = 0;
        long L = 0;
    }

    static class Pair {
        long p;
        int e;

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

    private static long exponentInFactorization(List<Pair> fac, long p) {
        for (Pair pair : fac) {
            if (pair.p == p)
                return pair.e;
        }
        return 0;
    }

    private static void evaluateL(long M, byte[] isPrime, long L, List<Pair> fac, BestState best) {
        long[] divisors = new long[8192];
        divisors[0] = 1;
        int size = 1;

        for (Pair pair : fac) {
            int cur = size;
            long pe = 1;
            for (int i = 1; i <= pair.e; ++i) {
                pe *= pair.p;
                for (int j = 0; j < cur; ++j) {
                    divisors[size++] = divisors[j] * pe;
                }
            }
        }

        long v2 = exponentInFactorization(fac, 2);
        int exp2 = (v2 == 0) ? 1 : (int) (v2 + 2);

        double logN = exp2 * Math.log(2.0);
        long modN = powMod(2, exp2, kLast);

        for (int i = 0; i < size; ++i) {
            long p = divisors[i] + 1;
            if (p == 2)
                continue;
            if (p > M + 1)
                continue;
            if (isPrime[(int) p] == 0)
                continue;

            int exp = (int) exponentInFactorization(fac, p) + 1;
            logN += exp * Math.log(p);
            modN = (modN * powMod(p, exp, kLast)) % kLast;
        }

        if (logN > best.logValue) {
            best.logValue = logN;
            best.modValue = modN;
            best.L = L;
        }
    }

    private static void dfsGenerate(long[] primes, long M, byte[] isPrime, int idx, int maxExp,
            long cur, List<Pair> fac, BestState best) {
        evaluateL(M, isPrime, cur, fac, best);
        if (idx >= primes.length)
            return;

        long p = primes[idx];
        long val = cur;
        for (int e = 1; e <= maxExp; ++e) {
            if (val > M / p)
                break;
            val *= p;
            fac.add(new Pair(p, e));
            dfsGenerate(primes, M, isPrime, idx + 1, e, val, fac, best);
            fac.remove(fac.size() - 1);
        }
    }

    private static long solveLast9(long boundMinus1) {
        long M = boundMinus1;
        byte[] isPrime = sieveIsPrime((int) (M + 1));
        long[] primes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43 };

        BestState best = new BestState();
        List<Pair> fac = new ArrayList<>();
        dfsGenerate(primes, M, isPrime, 0, 60, 1, fac, best);

        return (best.modValue + 1) % kLast;
    }

    public static void main(String[] args) {
        long last9 = solveLast9(20000000L - 1);
        System.out.printf("%09d\n", last9);
    }
}