Problem 133: Repunit Nonfactors

Project Euler Problem 133 Solution

EulerSolve provides an optimized solution for Project Euler Problem 133, Repunit Nonfactors, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $R(k)=\frac{10^k-1}{9}=11\ldots1$ be the repunit with $k$ digits. The problem asks for the sum of all primes $p<10^5$ such that there is no exponent $n\ge 1$ with $p\mid R(10^n)$. The key point is that one never has to build the enormous numbers $R(10^n)$ themselves. For each prime, the question can be decided entirely from the multiplicative behavior of $10$ modulo that prime. Mathematical Approach The implementations reduce the whole problem to a yes-or-no test on the multiplicative order of $10$ modulo $p$. Once that order is known, the classification is immediate. From repunits to a congruence For every integer $k\ge 1$, $$R(k)=1+10+10^2+\cdots+10^{k-1}=\frac{10^k-1}{9}.$$ If $p\notin\{2,3,5\}$, then $9$ is invertible modulo $p$, so multiplying by $9$ does not change divisibility by $p$. Therefore $$p\mid R(k)\iff 10^k\equiv 1 \pmod p.$$ In this problem the length is not an arbitrary $k$, but specifically $k=10^n$. So for primes other than $2$, $3$, and $5$, we are asking whether $$10^{10^n}\equiv 1 \pmod p$$ holds for at least one $n\ge 1$. Orders built only from $2$ and $5$ Let $$t=\operatorname{ord}_p(10),$$ the smallest positive integer such that $10^t\equiv 1\pmod p$. By definition, $10^m\equiv 1\pmod p$ happens exactly when $t\mid m$....

Detailed mathematical approach

Problem Summary

Let $R(k)=\frac{10^k-1}{9}=11\ldots1$ be the repunit with $k$ digits. The problem asks for the sum of all primes $p<10^5$ such that there is no exponent $n\ge 1$ with $p\mid R(10^n)$.

The key point is that one never has to build the enormous numbers $R(10^n)$ themselves. For each prime, the question can be decided entirely from the multiplicative behavior of $10$ modulo that prime.

Mathematical Approach

The implementations reduce the whole problem to a yes-or-no test on the multiplicative order of $10$ modulo $p$. Once that order is known, the classification is immediate.

From repunits to a congruence

For every integer $k\ge 1$,

$$R(k)=1+10+10^2+\cdots+10^{k-1}=\frac{10^k-1}{9}.$$

If $p\notin\{2,3,5\}$, then $9$ is invertible modulo $p$, so multiplying by $9$ does not change divisibility by $p$. Therefore

$$p\mid R(k)\iff 10^k\equiv 1 \pmod p.$$

In this problem the length is not an arbitrary $k$, but specifically $k=10^n$. So for primes other than $2$, $3$, and $5$, we are asking whether

$$10^{10^n}\equiv 1 \pmod p$$

holds for at least one $n\ge 1$.

Orders built only from $2$ and $5$

Let

$$t=\operatorname{ord}_p(10),$$

the smallest positive integer such that $10^t\equiv 1\pmod p$. By definition, $10^m\equiv 1\pmod p$ happens exactly when $t\mid m$. Hence

$$p\mid R(10^n)\iff t\mid 10^n.$$

Now $10^n=2^n5^n$, so its only prime divisors are $2$ and $5$. Therefore $t\mid 10^n$ for some $n$ if and only if every prime divisor of $t$ is either $2$ or $5$. Equivalently, after removing all factors $2$ and $5$ from $t$, nothing remains.

So for every prime $p\notin\{2,3,5\}$,

$$p\text{ divides some }R(10^n)\iff \operatorname{ord}_p(10)=2^a5^b\text{ for some }a,b\ge 0.$$

The converse is just as important as the forward implication: if $t=2^a5^b$, then choosing $n\ge \max(a,b)$ gives $t\mid 10^n$, so such a prime really does divide one of the required repunits.

Special primes and worked examples

The primes $2$ and $5$ are immediate. Every repunit ends in the digit $1$, so it is divisible by neither $2$ nor $5$. Both must therefore be included in the final sum.

The prime $3$ is special for a different reason: $9$ is not invertible modulo $3$, so the previous equivalence must not be used. Instead,

$$R(k)=1+10+\cdots+10^{k-1}\equiv 1+1+\cdots+1\equiv k \pmod 3.$$

Since $10^n\equiv 1\pmod 3$, we get $R(10^n)\equiv 1\pmod 3$. Thus $3$ never divides any $R(10^n)$, even though it does divide other repunits such as $R(3)$.

For a typical nontrivial example, $p=17$ has $\operatorname{ord}_{17}(10)=16=2^4$. Because $16\mid 10^4$, the prime $17$ divides $R(10^4)$, so it is excluded from the sum.

By contrast, $p=19$ has $\operatorname{ord}_{19}(10)=18=2\cdot 3^2$. The extra factor $3$ means $18\nmid 10^n$ for every $n$, so $19$ never divides any $R(10^n)$ and must be included.

A useful sanity check is $p=7$. Its order is $6=2\cdot 3$, so $7$ divides $R(6)$ but never divides $R(10^n)$. That is exactly the distinction this problem is testing.

Why reducing divisors of $p-1$ finds the exact order

For every prime $p\notin\{2,5\}$, Fermat's little theorem gives $10^{p-1}\equiv 1\pmod p$. Therefore $\operatorname{ord}_p(10)$ is a divisor of $p-1$. This is why the implementations start from the candidate value $p-1$ instead of searching upward from $1$.

Suppose the current candidate is $d$, and $q$ is a prime divisor of $d$. If

$$10^{d/q}\equiv 1 \pmod p,$$

then the true order already divides $d/q$, so the factor $q$ can be removed. Repeating that test as long as it succeeds preserves the invariant that the true order divides the current candidate.

When no prime factor can be removed any further, the remaining candidate is the smallest positive exponent producing $1$, hence it is exactly $\operatorname{ord}_p(10)$. After that, the final classification is obtained by stripping powers of $2$ and $5$ from this order and checking whether the remainder is $1$.

How the Code Works

The C++, Python, and Java implementations first generate all primes below the limit with a sieve. As each prime is processed, the primes $2$, $3$, and $5$ are handled immediately as guaranteed nonfactors of every $R(10^n)$.

For every other prime $p$, the implementation starts from $p-1$, factors that number into its distinct prime divisors, and repeatedly lowers the candidate order whenever a modular-exponentiation test shows that a smaller divisor still works. This produces the exact multiplicative order of $10$ modulo $p$.

Once the order has been found, the code removes factors of $2$ and $5$. If the reduced value is $1$, then the prime can divide some $R(10^n)$ and is skipped. Otherwise the prime is added to the running total. The C++ and Python versions perform modular exponentiation by repeated squaring; the Java version uses the standard big-integer modular power routine, but the mathematical test is the same in all three languages.

Complexity Analysis

If the search limit is $L$, the sieve phase costs $O(L\log\log L)$ time and $O(L)$ space. Here $L=10^5$, so this part is tiny.

For a prime $p$, the order computation factors $p-1$ by trial division up to $\sqrt{p-1}$, and each successful or failed reduction test uses modular exponentiation in $O(\log p)$ multiplications. Since there are only $9592$ primes below $10^5$, the total running time is comfortably small, and the memory usage is dominated by the sieve array.

Footnotes and References

Problem page: Project Euler 133
Repunit: Wikipedia - Repunit
Multiplicative order: Wikipedia - Multiplicative order
Fermat's little theorem: Wikipedia - Fermat's little theorem
Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes

Problem 133 source code

C++

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

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

struct Options {
    int limit = 100000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(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_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.limit >= 2;
}

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

u64 pow_mod(u64 base, u64 exp, u64 mod) {
    u64 result = 1 % mod;
    base %= mod;
    while (exp > 0) {
        if ((exp & 1ULL) != 0ULL) {
            result = mul_mod(result, base, mod);
        }
        base = mul_mod(base, base, mod);
        exp >>= 1ULL;
    }
    return result;
}

int multiplicative_order_10(const int p) {
    int order = p - 1;
    int x = order;

    std::vector<int> factors;
    for (int f = 2; static_cast<std::int64_t>(f) * f <= x; ++f) {
        if (x % f != 0) {
            continue;
        }
        factors.push_back(f);
        while (x % f == 0) {
            x /= f;
        }
    }
    if (x > 1) {
        factors.push_back(x);
    }

    for (int q : factors) {
        while ((order % q) == 0 && pow_mod(10ULL, static_cast<u64>(order / q), static_cast<u64>(p)) == 1ULL) {
            order /= q;
        }
    }

    return order;
}

std::vector<int> primes_below(const int limit) {
    std::vector<bool> sieve(static_cast<std::size_t>(limit), true);
    if (limit > 0) {
        sieve[0] = false;
    }
    if (limit > 1) {
        sieve[1] = false;
    }

    for (int i = 2; static_cast<std::int64_t>(i) * i < limit; ++i) {
        if (!sieve[static_cast<std::size_t>(i)]) {
            continue;
        }
        for (int j = i * i; j < limit; j += i) {
            sieve[static_cast<std::size_t>(j)] = false;
        }
    }

    std::vector<int> primes;
    for (int i = 2; i < limit; ++i) {
        if (sieve[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
        }
    }
    return primes;
}

bool can_divide_some_R10n(const int p) {
    if (p == 2 || p == 3 || p == 5) {
        return false;
    }

    int ord = multiplicative_order_10(p);
    while ((ord % 2) == 0) {
        ord /= 2;
    }
    while ((ord % 5) == 0) {
        ord /= 5;
    }

    return ord == 1;
}

std::int64_t solve(const int limit) {
    const std::vector<int> primes = primes_below(limit);

    std::int64_t sum = 0;
    for (int p : primes) {
        if (!can_divide_some_R10n(p)) {
            sum += p;
        }
    }

    return sum;
}

bool run_checkpoints() {
    if (solve(100) != 918) {
        std::cerr << "Checkpoint failed for limit=100" << '\n';
        return false;
    }
    if (!can_divide_some_R10n(17) || can_divide_some_R10n(19) || can_divide_some_R10n(3)) {
        std::cerr << "Checkpoint failed for prime behavior (3, 17, 19)" << '\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.limit) << '\n';
    return 0;
}

Python

import sys
from typing import List

def parse_int_after_prefix(arg: str, prefix: str) -> (bool, int):
    if not arg.startswith(prefix):
        return False, 0
    tail = arg[len(prefix):]
    if not tail:
        return False, 0
    
    try:
        parsed = int(tail)
        return True, parsed
    except ValueError:
        return False, 0

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

def pow_mod(base: int, exp: int, mod: int) -> int:
    result = 1 % mod
    base %= mod
    while exp > 0:
        if exp & 1:
            result = mul_mod(result, base, mod)
        base = mul_mod(base, base, mod)
        exp >>= 1
    return result

def multiplicative_order_10(p: int) -> int:
    order = p - 1
    x = order
    
    factors = []
    f = 2
    while f * f <= x:
        if x % f == 0:
            factors.append(f)
            while x % f == 0:
                x //= f
        f += 1
    if x > 1:
        factors.append(x)
    
    for q in factors:
        while order % q == 0 and pow_mod(10, order // q, p) == 1:
            order //= q
    
    return order

def primes_below(limit: int) -> List[int]:
    if limit <= 2:
        return []
    
    sieve = [True] * limit
    sieve[0] = False
    sieve[1] = False
    
    i = 2
    while i * i < limit:
        if sieve[i]:
            j = i * i
            while j < limit:
                sieve[j] = False
                j += i
        i += 1
    
    return [i for i in range(2, limit) if sieve[i]]

def can_divide_some_R10n(p: int) -> bool:
    if p in (2, 3, 5):
        return False
    
    ord_val = multiplicative_order_10(p)
    
    while ord_val % 2 == 0:
        ord_val //= 2
    while ord_val % 5 == 0:
        ord_val //= 5
    
    return ord_val == 1

def solve(limit: int) -> int:
    primes = primes_below(limit)
    
    total = 0
    for p in primes:
        if not can_divide_some_R10n(p):
            total += p
    
    return total

def run_checkpoints() -> bool:
    if solve(100) != 918:
        return False
    if not can_divide_some_R10n(17) or can_divide_some_R10n(19) or can_divide_some_R10n(3):
        return False
    return True

def main():
    limit = 100000
    run_checkpoints_flag = True
    
    args = sys.argv[1:]
    for arg in args:
        if arg == "--skip-checkpoints":
            run_checkpoints_flag = False
        elif arg.startswith("--limit="):
            success, value = parse_int_after_prefix(arg, "--limit=")
            if success:
                limit = value
        else:
            print(f"Unknown argument: {arg}", file=sys.stderr)
            sys.exit(1)
    
    if limit < 2:
        print("Limit must be at least 2", file=sys.stderr)
        sys.exit(1)
    
    if run_checkpoints_flag and not run_checkpoints():
        print("Checkpoint failed", file=sys.stderr)
        sys.exit(2)
    
    print(solve(limit))

if __name__ == "__main__":
    main()

Java

import java.math.BigInteger;
import java.util.*;

public class Euler133 {
    public static void main(String[] args) {
        int limit = 100000;
        boolean[] sieve = new boolean[limit];
        Arrays.fill(sieve, true);
        sieve[0] = sieve[1] = false;
        for (int i = 2; i * i < limit; i++)
            if (sieve[i])
                for (int j = i * i; j < limit; j += i)
                    sieve[j] = false;

        BigInteger TEN = BigInteger.TEN;
        long sum = 0;
        for (int p = 2; p < limit; p++) {
            if (!sieve[p])
                continue;
            if (p == 2 || p == 3 || p == 5) {
                sum += p;
                continue;
            } // these can't divide
            // Find multiplicative order of 10 mod p
            int order = p - 1;
            // Factor order
            List<Integer> factors = new ArrayList<>();
            int x = order;
            for (int f = 2; (long) f * f <= x; f++) {
                if (x % f == 0) {
                    factors.add(f);
                    while (x % f == 0)
                        x /= f;
                }
            }
            if (x > 1)
                factors.add(x);
            BigInteger mod = BigInteger.valueOf(p);
            for (int q : factors)
                while (order % q == 0 && TEN.modPow(BigInteger.valueOf(order / q), mod).equals(BigInteger.ONE))
                    order /= q;
            // Check if order is 2^a * 5^b
            int o = order;
            while (o % 2 == 0)
                o /= 2;
            while (o % 5 == 0)
                o /= 5;
            if (o != 1)
                sum += p;
        }
        System.out.println(sum);
    }
}