Problem 132: Large Repunit Factors

Project Euler Problem 132 Solution

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

Problem Summary Let $R_n$ denote the base-10 repunit with $n$ digits, so $$R_n = 11\ldots1 = 1 + 10 + 10^2 + \cdots + 10^{n-1} = \frac{10^n - 1}{9}.$$ The task is to find the sum of the first 40 prime numbers dividing $R_{10^9}$. Since $R_{10^9}$ has a billion decimal digits, direct construction or factorization is completely out of the question. The solution works by testing primes through modular arithmetic instead of touching the gigantic integer itself. Mathematical Approach The whole problem becomes manageable once divisibility of repunits is rewritten in terms of powers of 10 modulo a small number. After that, the only remaining job is to scan primes in increasing order and keep the first 40 that satisfy the derived condition. Repunits as a geometric series Because $R_n$ is a finite geometric sum, multiplying by 9 gives $$10^n - 1 = 9R_n.$$ If $p$ is a prime different from $2$ and $5$, then $10$ is invertible modulo $9p$, and the divisibility test becomes $$p \mid R_n \iff 9p \mid (10^n - 1) \iff 10^n \equiv 1 \pmod{9p}.$$ This is exactly the condition used by the Python and Java implementations. The primes $2$ and $5$ can be discarded immediately, because every repunit ends in the digit 1 and therefore is divisible by neither 2 nor 5....

Detailed mathematical approach

Problem Summary

A repunit is a number written with the digit $1$ only, and $R(k)$ is the repunit with $k$ digits. The statement gives the example $R(10)=1111111111$${}=11\cdot41\cdot271\cdot9091$, whose prime factors add up to $9414$. We must find the sum of the first forty prime factors of $R(10^9)$, that is, of the forty smallest different primes that divide this number of one billion digits. ^[1]

The whole solution in one sentence: a prime $p$ divides $R(n)$ exactly when $10^n$ leaves the remainder $1$ after division by $9p$, this remainder takes only about thirty squarings to compute even for $n=10^9$, and so we test the primes $3,7,11,13,\dots$ in increasing order and add up the first forty that pass.

Mathematical Approach

1. Small repunits

A repunit is a sum of powers of ten, a geometric series: ^[2]

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

Factoring the first repunits shows a pattern:

$R(2)=11$, $R(3)=3\cdot37$, $R(4)=11\cdot101$, $R(5)=41\cdot271$;
$R(6)=3\cdot7\cdot11\cdot13\cdot37$, $R(8)=11\cdot73\cdot101\cdot137$;
$R(10)=11\cdot41\cdot271\cdot9091$, $R(12)=3\cdot7\cdot11\cdot13$${}\cdot37\cdot101\cdot9901$.

In this list $11$ divides exactly the repunits of even length, $41$ and $271$ appear at the lengths $5$ and $10$, and $3$ and $37$ at the multiples of $3$. Each prime seems to come with a length of its own and to divide exactly the repunits whose length is a multiple of it; Sections 2 to 4 prove this.

2. A test with remainders

We write $a\equiv b\pmod m$ when $m$ divides $a-b$, that is, when $a$ and $b$ leave the same remainder after division by $m$; such congruences can be added, multiplied and raised to powers like equations. ^[3] Since $9R(n)=10^n-1$, every prime $p$ satisfies

$$p\mid R(n)\iff 9p\mid 10^n-1\iff 10^n\equiv1\pmod{9p}.$$

Indeed, $R(n)=p\,t$ gives $10^n-1=9p\,t$, and conversely. For example $10^6-1=999999=63\cdot15873$, so $10^6\equiv1\pmod{63}$, and indeed $R(6)=111111=7\cdot15873$. The primes $2$ and $5$ never divide a repunit, because a repunit ends in the digit $1$.

3. The prime 3

Since $10\equiv1\pmod9$, every power of ten leaves the remainder $1$ after division by $9$, and the $n$ terms of the sum in Section 1 give $R(n)\equiv n\pmod9$. So $3$ divides $R(n)$ exactly when $3$ divides $n$, as for $R(3)$, $R(6)$ and $R(12)$ above. ^[2] Because $10^9=3\cdot333333333+1$, the prime $3$ does not divide $R(10^9)$. The test of Section 2 agrees: $10^3=1000=37\cdot27+1$, so $10^3\equiv1\pmod{27}$ and $10^{10^9}=(10^3)^{333333333}\cdot10$${}\equiv10\pmod{27}$.

The factor $9$ in the modulus matters here: modulo $3$ alone every power of ten leaves the remainder $1$, so the shorter test $10^n\equiv1\pmod p$ would wrongly accept $3$ for every $n$.

4. The order of 10

For a prime $p$ other than $2$, $3$ and $5$ the factor $9$ can be dropped: $9$ always divides $10^n-1$ by Section 3, and $9$ and $p$ have no common factor, so $9p$ divides $10^n-1$ exactly when $p$ does. For $p=7$ the powers $10^1,\dots,10^6$ leave the remainders $3,2,6,4,5,1$, each the previous one times $10$ reduced modulo $7$, so $1$ first appears at the exponent $6$. Such an exponent always exists, because $10^{p-1}\equiv1\pmod p$ by Fermat's little theorem. ^[4] The smallest exponent $k\ge1$ with $10^k\equiv1\pmod p$ is called the multiplicative order $\operatorname{ord}_p(10)$. ^[5]

Claim. For such a prime, $10^n\equiv1\pmod p$ holds exactly when $\operatorname{ord}_p(10)$ divides $n$.

Proof. Let $k=\operatorname{ord}_p(10)$ and divide $n$ by $k$ with remainder: $n=qk+r$ with $0\le r\lt k$. Then $10^n=(10^k)^q\cdot10^r\equiv10^r\pmod p$. If $r=0$, this is $1$. If $0\lt r\lt k$, it is not $1$, because $k$ is the smallest exponent that gives $1$. $\square$

So $p$ divides $R(n)$ exactly when $\operatorname{ord}_p(10)$ divides $n$, which is the pattern of Section 1: the orders of $7,11,13,37,41,73,101,137$ are $6,2,6,3,5,8,4,8$. The claim with $n=p-1$ shows that $\operatorname{ord}_p(10)$ divides $p-1$. The order is also the period of the decimal expansion of $1/p$; for example $1/7=0.142857\,142857\ldots$ repeats after $6$ digits. ^[6] Finally, if $m$ divides $n$, then $R(m)$ divides $R(n)$, ^[2] because the $n$ ones split into $n/m$ blocks of $m$ ones:

$$R(n)=R(m)\,\bigl(1+10^m+10^{2m}+\dots+10^{n-m}\bigr).$$

Since $10$ divides $10^9$, the four primes of the statement's example all divide $R(10^9)$.

5. Which primes divide $R(10^9)$

The divisors of $10^9=2^9\cdot5^9$ are the numbers $2^a5^b$ with $0\le a,b\le9$. A prime $p$ other than $2$, $3$ and $5$ therefore divides $R(10^9)$ exactly when

$$\operatorname{ord}_p(10)=2^a\,5^b,\qquad 0\le a,b\le9.$$

$7$ and $13$ (order $6$), $19$ ($18$), $31$ ($15$), $37$ ($3$) and $61$ ($60$) fail, because their orders contain the factor $3$.
The ten smallest prime factors of $R(10^9)$ are $11$, $17$, $41$, $73$, $101$, $137$, $251$, $257$, $271$ and $353$, with the orders $2,16,5,8,4,8,50,256,5,32$.
$65537$ has the order $65536=2^{16}$, a power of $2$, but $2^{16}$ does not divide $10^9$, so $65537$ fails: the bounds on $a$ and $b$ matter.

There are more than forty such primes: by Zsigmondy's theorem, ^[7] for every $d\ge2$ some prime divides $10^d-1$ but no $10^k-1$ with $1\le k\lt d$. It has order exactly $d$ and is not $3$, which divides $10^1-1$, so it divides $R(d)$, and hence $R(10^9)$ when $d$ divides $10^9$. With $10$ choices for each of $a$ and $b$, $10^9$ has $100$ divisors, $99$ of them with $d\ge2$, and since a prime has only one order, they give $99$ different prime factors. The fortieth smallest turns out to be $160001$, with order $1250$.

6. Powers by repeated squaring

Exponentiation by squaring finds the remainder of $10^{10^9}$ modulo $9p$ without forming the huge power. ^[8] Write the exponent in binary and start with the value $10$, which is $10^e$ for the leading digit, $e=1$. For each further digit, from left to right, square the value, which doubles $e$, and if the digit is $1$ also multiply by $10$, which adds $1$ to $e$; after every step only the remainder modulo $m=9p$ is kept. For $6=110_2$ and $p=7$, so $m=63$, the value $10$ becomes $10^2\cdot10=1000\equiv55$, with $e=3$, and then $55^2=3025=48\cdot63+1\equiv1$, with $e=6$. The exponent $10^9$ has $30$ binary digits, $13$ of them ones, so one test needs $29$ squarings and $12$ multiplications by $10$, all with numbers below $9p$.

7. The repunit itself by doubling

The C++ program finds the remainder of $R(n)$ modulo $p$ directly. Writing $a$ ones followed by $b$ ones gives $R(a+b)=R(a)\cdot10^b+R(b)$, and so

$$R(2m)=R(m)\,(10^m+1),\qquad R(2m+1)=10\,R(2m)+1.$$

The program keeps the pair $\bigl(R(m)\bmod p,\ 10^m\bmod p\bigr)$, starting from $m=0$ with $(0,1)$, since $R(0)=0$ is the empty sum and $10^0=1$. Every binary digit of $n$, read from the left, doubles $m$ with the first formula, and a digit $1$ then adds one more one with the second; the power of ten is squared, and multiplied by $10$, along with it. For $n=10=1010_2$ and $p=41$ the lengths are $1,2,4,5,10$:

$m=1$: $(1,10)$; $m=2$: $1\cdot10+1=11$ and $10^2=100\equiv18$, so $(11,18)$;
$m=4$: $11\cdot18+11=209\equiv4$ and $18^2=324\equiv37$, so $(4,37)$;
$m=5$: $4\cdot10+1=41\equiv0$ and $37\cdot10=370\equiv1$, so $(0,1)$; $m=10$: $(0,1)$ again.

So $41$ divides $R(5)$ and $R(10)$, in agreement with $\operatorname{ord}_{41}(10)=5$. This form of the test works for every prime, $3$ included, without the factor $9$.

8. The search

The first forty prime factors are the forty smallest, so the primes are tested in increasing order, $3,7,11,13,17,\dots$, skipping $2$ and $5$, until forty have passed. The primes come from trial division: a number $q\ge2$ that is not prime has a prime factor at most $\sqrt q$, because if $q=ab$ with $2\le a\le b$, then $a^2\le q$, and every prime factor of $a$ divides $q$. The search ends at the fortieth success, $p=160001$, after $14682$ tests.

A shortcut, not used by the programs: since $\operatorname{ord}_p(10)$ also divides $p-1$, for $p\ne3$ it would suffice to test $10^g\equiv1\pmod p$ with the smaller exponent $g=\gcd(p-1,10^9)$.

How the Code Works

The C++ program ^[9] starts in main with the default Options: forty primes, the exponent $10^9$ and the checkpoints on. parse_arguments accepts --wanted=, --exponent= and --skip-checkpoints and rejects other arguments and counts below $1$ with the exit code $1$; its digit parsers have no overflow check, which matters only for absurdly large values.

mul_mod multiplies two remainders as unsigned __int128, a compiler extension, ^[10] so no product overflows for any 64-bit modulus. repunit_mod(n, mod) is Section 7: from repunit $=0$ and pow10 $=1$ it walks from the highest binary digit of $n$ down to the last, doubling at every digit and adding one more one at every digit $1$; divides_repunit checks that the result is $0$.

solve(wanted, exponent) runs the search of Section 8 over the odd numbers $p=3,5,7,9,\dots$: is_prime divides $p$ by the stored primes until their square exceeds $p$, which suffices because the vector primes holds every prime below $p$. New primes are stored, $5$ is skipped, and every other prime, $3$ included, is tested; a success increments found and adds $p$ to a 64-bit sum. The loop has no upper bound: for $10^9$ it stops at $p=160001$, far inside the range of the int p, but for a repunit with fewer prime factors than requested, such as $R(1)=1$, it never stops. Before the search, run_checkpoints checks solve(4, 10) $=9414$, the statement's example, solve(4, 6) $=3+7+11+13=34$ from $R(6)=3\cdot7\cdot11\cdot13\cdot37$, and that $3$ divides $R(3)$ but not $R(10)$, with the exit code $2$ on failure; the message of the second check names $R(10^6)$, although it checks $R(6)$. Finally main prints the sum.

The Python program, which fixes $n=10^9$ and forty primes and runs no checks, walks through the primes from $2$, skips $2$ and $5$ and tests pow(10, exp, 9 * p) == 1, the condition of Section 2, with the built-in modular power on arbitrary-precision integers; ^[11] the next prime comes from trial division by every $d=2,3,4,\dots$ up to $\sqrt p$. The Java program, also without checks, tries every int p from $2$ upwards with isPrime, which divides by $2$ and the odd numbers up to $\sqrt p$, skips $2$ and $5$, and tests TEN.modPow(EXP, mod) ^[11] with the BigInteger modulus 9L * p, adding to a long sum. Thanks to the modulus $9p$ both reject $3$; the C++ version needs no factor $9$ because it reduces $R(n)$ itself.

Complexity and Verification

Let $P=160001$ be the fortieth qualifying prime. With exponentiation by squaring each test costs $O(\log n)$ multiplications of numbers below $9P$: the C++ function calls mul_mod twice for each of the $30$ binary digits of $10^9$ and twice more for each of its $13$ ones, $86$ times per prime, while Python and Java use library implementations of modular exponentiation. ^[11] Trial division of $q$ needs at most about $\sqrt q$ divisions, so the whole search takes $O\bigl(P\sqrt P+\pi(P)\log n\bigr)$ operations, where $\pi(P)=14684$ is the number of primes up to $P$, all stored in the C++ vector; Python and Java store no list. In the main search the C++ program makes $1262652$ calls of mul_mod for $14682$ tests and $1260440$ trial divisions, by primes only; Java makes $2746342$ and Python $5260558$ trial divisions, Python's including one more next-prime search, up to $160009$. On a desktop computer each program finishes in well under a second.

For this article all three programs were run; they print the same number. Separate Python checks, not part of the published solutions, confirm the mathematics. For every $n$ from $1$ to $60$ and every prime $p\lt20000$, $135720$ cases, the doubling of Section 7 returns exactly $R(n)\bmod p$, which is $0$ exactly when $10^n\bmod 9p=1$; the two tests also agree on $20000$ random pairs with $n\lt10^{12}$. For every prime from $7$ to $160001$ they agree with the order criterion of Section 5, using independently computed orders, and with the shortcut of Section 8, and exactly forty primes pass, while both reject $3$. Finally $10^{10^9}\not\equiv1\pmod{9p^2}$ for each of the forty primes, so none divides $R(10^9)$ twice, and the list is the same with or without repetition.

Footnotes and References

The references below supply the background facts used above. The divisibility test with the modulus $9p$, the order criterion, the doubling formulas and the operation counts are derived in this article.

Project Euler 132 — Large Repunit Factors. The official statement defines R(k) as the repunit of length k, gives R(10) = 11 × 41 × 271 × 9091 with the prime factor sum 9414, and asks for the sum of the first forty prime factors of R(10⁹).
Wikipedia — Repunit. The decimal repunit Rₙ = (10ⁿ − 1)/9 consists of n ones; Rₐ divides Rₙ when a divides n; and Rₙ leaves the same remainder modulo 3 as n.
Wikipedia — Modular arithmetic. Integers are congruent modulo m when their difference is a multiple of m; congruences can be added, subtracted, multiplied and raised to powers.
Wikipedia — Fermat's little theorem. If p is prime and a is not divisible by p, then a^(p − 1) ≡ 1 (mod p).
Wikipedia — Multiplicative order. For a coprime to n, the order of a modulo n, written ordₙ(a), is the smallest positive integer k with aᵏ ≡ 1 (mod n).
Wikipedia — Repeating decimal. For a prime p other than 2 and 5, the length of the repetend of 1/p equals the order of 10 modulo p; 1/7 = 0.142857… has 6 repeating digits.
Wikipedia — Zsigmondy's theorem. For coprime a > b > 0, aⁿ − bⁿ has a prime divisor that divides no aᵏ − bᵏ with k < n, except for n = 1 with a − b = 1, n = 2 with a + b a power of two, and n = 6 with a = 2, b = 1; none of these applies to a = 10, b = 1.
Wikipedia — Exponentiation by squaring. xⁿ is computed with about log₂ n squarings and one multiplication fewer than there are ones in the binary expansion of n.
C++ — Euler132.cpp. The linked immutable C++ revision contains repunit_mod, solve and run_checkpoints. It is the implementation record for this article.
GCC manual — 128-bit Integers. The type unsigned __int128 is supported on targets with an integer mode wide enough to hold 128 bits.
Python — pow with three arguments returns base to the power exp modulo mod, computed more efficiently than pow(base, exp) % mod; Java — BigInteger.modPow returns this BigInteger to the power exponent, modulo m.

Problem 132 source code

C++

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

namespace {

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

struct Options {
    int wanted = 40;
    u64 exponent = 1000000000ULL;
    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_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, "--wanted=", options.wanted) ||
            parse_u64_after_prefix(arg, "--exponent=", options.exponent)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.wanted >= 1;
}

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

u64 repunit_mod(u64 n, u64 mod) {
    if (mod == 1) {
        return 0;
    }

    u64 pow10 = 1 % mod;
    u64 repunit = 0;
    const u64 ten_mod = 10 % mod;

    int msb = 63;
    while (msb > 0 && ((n >> msb) & 1ULL) == 0ULL) {
        --msb;
    }

    for (int bit = msb; bit >= 0; --bit) {
        repunit = (mul_mod(repunit, pow10, mod) + repunit) % mod;
        pow10 = mul_mod(pow10, pow10, mod);

        if (((n >> bit) & 1ULL) != 0ULL) {
            repunit = (mul_mod(repunit, ten_mod, mod) + 1ULL) % mod;
            pow10 = mul_mod(pow10, ten_mod, mod);
        }
    }

    return repunit;
}

bool is_prime(int n, const std::vector<int>& primes) {
    if (n <= 1) {
        return false;
    }
    for (int p : primes) {
        if (static_cast<std::int64_t>(p) * p > n) {
            break;
        }
        if ((n % p) == 0) {
            return false;
        }
    }
    return true;
}

bool divides_repunit(int p, u64 exponent) {
    return repunit_mod(exponent, static_cast<u64>(p)) == 0ULL;
}

u64 solve(const int wanted, const u64 exponent) {
    int found = 0;
    u64 sum = 0;
    std::vector<int> primes;
    primes.reserve(4096);
    primes.push_back(2);

    if (wanted <= 0) {
        return 0;
    }

    for (int p = 3; found < wanted; p += 2) {
        if (!is_prime(p, primes)) {
            continue;
        }
        primes.push_back(p);

        if (p == 5) {
            continue;
        }
        if (divides_repunit(p, exponent)) {
            ++found;
            sum += static_cast<u64>(p);
        }
    }

    return sum;
}

bool run_checkpoints() {
    if (solve(4, 10ULL) != 9414ULL) {
        std::cerr << "Checkpoint failed for first 4 factors of R(10)" << '\n';
        return false;
    }
    if (solve(4, 6ULL) != 34ULL) {
        std::cerr << "Checkpoint failed for first 4 factors of R(10^6)" << '\n';
        return false;
    }
    if (divides_repunit(3, 10ULL) || !divides_repunit(3, 3ULL)) {
        std::cerr << "Checkpoint failed for prime 3 divisibility behavior" << '\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.wanted, options.exponent) << '\n';
    return 0;
}

Python

# Problem 132: Large repunit factors
# Sum of first 40 primes that are factors of R(10^9).

def solve():
    exp = 10**9
    found = 0
    total = 0
    p = 2
    while found < 40:
        if p not in (2, 5):
            # Check if R(exp) % p == 0, i.e. pow(10, exp, 9*p) == 1
            if pow(10, exp, 9 * p) == 1:
                found += 1
                total += p
        # next prime
        p += 1
        while True:
            is_p = p >= 2
            d = 2
            while d * d <= p:
                if p % d == 0: is_p = False; break
                d += 1
            if is_p: break
            p += 1
    print(total)

solve()

Java

import java.math.BigInteger;

public class Euler132 {
    static boolean isPrime(int n) {
        if (n < 2)
            return false;
        if (n % 2 == 0)
            return n == 2;
        for (int p = 3; (long) p * p <= n; p += 2)
            if (n % p == 0)
                return false;
        return true;
    }

    public static void main(String[] args) {
        long exp = 1000000000L;
        int found = 0;
        long sum = 0;
        BigInteger TEN = BigInteger.TEN, EXP = BigInteger.valueOf(exp);
        for (int p = 2; found < 40; p++) {
            if (!isPrime(p))
                continue;
            if (p == 2 || p == 5)
                continue;
            BigInteger mod = BigInteger.valueOf(9L * p);
            if (TEN.modPow(EXP, mod).equals(BigInteger.ONE)) {
                found++;
                sum += p;
            }
        }
        System.out.println(sum);
    }
}