Problem 274: Divisibility Multipliers

Project Euler Problem 274 Solution

EulerSolve provides an optimized solution for Project Euler Problem 274, Divisibility Multipliers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For every prime $p\neq2,5$, define the decimal divisibility multiplier $m(p)$ as the number that lets us replace a divisibility test on $$n=10q+r$$ by a test on the shorter number $q+m(p)r$. The problem asks for the sum of these multipliers over all primes below a large limit. The final numeric answer is omitted here; what matters is the mathematical reason this multiplier exists and why the code computes it so quickly. Mathematical Approach 1) Why primes $2$ and $5$ are excluded Any decimal integer can be split into quotient and last digit as $$n=10q+r,\qquad 0\le r\le 9.$$ If $p=2$ or $p=5$, then $10\equiv0\pmod p$, so $10$ has no multiplicative inverse modulo $p$. That means the whole “remove the last digit and adjust” trick cannot be expressed by multiplying with $10^{-1}$. For every other prime $p$, we have $\gcd(10,p)=1$, so a unique inverse of $10$ exists modulo $p$. 2) Deriving the multiplier formula Start from the divisibility condition $$p\mid n \iff p\mid (10q+r).$$ Let $t$ be the modular inverse of $10$, so $$10t\equiv1\pmod p.$$ Multiply the congruence $10q+r\equiv0\pmod p$ by $t$: $$q+tr\equiv0\pmod p.$$ Therefore $$p\mid n \iff p\mid (q+tr).$$ So the divisibility multiplier is exactly $$\boxed{m(p)\equiv10^{-1}\pmod p.}$$ The code returns the least nonnegative representative of this inverse....

Detailed mathematical approach

Problem Summary

For an integer $p\gt1$ with no common factor with $10$, a positive integer $m\lt p$ is called a divisibility multiplier for $p$ if, for every positive integer $n$, the number $f(n)$ obtained by deleting the last digit of $n$ and adding $m$ times that digit to what is left is divisible by $p$ exactly when $n$ is. For $113$ the multiplier is $34$: $f(76275)=7627+5\cdot34=7797$, and both numbers are divisible by $113$, while $f(12345)=1234+5\cdot34=1404$, and neither is. For the primes coprime to $10$, that is, all primes except $2$ and $5$, the multipliers below $1000$ add up to $39517$, and we must find their sum below $10^7$. ^[1]

The whole solution in one sentence: the multiplier of $p$ is the inverse of $10$ modulo $p$, the only $m$ between $1$ and $p-1$ for which $10m-1$ is a multiple of $p$, so we sieve the primes below $10^7$, skip $2$ and $5$, and add up these inverses, each found in a few steps of the extended Euclidean algorithm.

Mathematical Approach

We find the multipliers of small primes by trial, prove the pattern they show, and then explain how the programs find the inverses and the primes.

1. Cutting off the last digit

Write $n=10q+r$ with $0\le r\le9$. Then $r$ is the last digit of $n$, and $q=\lfloor n/10\rfloor$, the largest integer not exceeding $n/10$, is the number formed by the other digits, or $0$ if $n$ has one digit. So

$$f(n)=q+m\,r,$$

and for $76275$ we have $q=7627$ and $r=5$. We write $a\equiv b\pmod p$ when $a-b$ is a multiple of $p$, that is, when $a$ and $b$ leave the same remainder after division by $p$. ^[2]

2. Small cases by trial

For $p=7$ we can try all six candidates. With $m=1$ the multiple $14$ of $7$ goes to $1+4=5$, which is not a multiple; with $m=2$ the non-multiple $13$ goes to the multiple $1+3\cdot2=7$, which breaks the other direction of "exactly when". Likewise $m=3$, $4$ and $6$ fail at $n=12$, $14$ and $11$, and only $m=5$ passes for every $n\le700$. For $p=3,7,11,13,17,19,23$ the same trial gives $m=1,5,10,4,12,2,7$, and in every case $10m-1$ is a multiple of $p$: $9=3\cdot3$, $49=7\cdot7$, $99=9\cdot11$, $39=3\cdot13$, $119=7\cdot17$, $19=1\cdot19$ and $69=3\cdot23$.

So $10m\equiv1\pmod p$: for remainders after division by $p$, multiplying by $m$ undoes multiplying by $10$, and $m$ is called an inverse of $10$ modulo $p$. ^[2]

3. A multiplier for every $p$

Since $p$ is coprime to $10$, its last digit is $1$, $3$, $7$ or $9$. Multiplying $p$ by $k=9$, $3$, $7$ or $1$ respectively gives a number $kp$ that ends in $9$, because the last digit of a product depends only on the last digits of the factors, and $1\cdot9$, $3\cdot3$, $7\cdot7=49$ and $9\cdot1$ end in $9$. Hence

$$m=\frac{kp+1}{10}$$

is an integer with $10m-1=kp$, and $1\le m\lt p$ because $10\le kp+1\le9p+1\lt10p$. For $113$ we get $k=3$, $kp=339$ and $m=34$; for $19$, $k=1$ and $m=2$. This is the generalised divisibility rule for divisors ending in $1$, $3$, $7$ or $9$, ^[4] and the OEIS lists the inverses of $10$ modulo the primes from $7$ on with the same formula. ^[5]

4. Why the inverse works

Claim. If $1\le m\lt p$ and $10m-1=kp$ for an integer $k$, then $m$ is a divisibility multiplier for $p$.

Proof. Multiplying $f(n)$ by $10$ gives

$$\begin{aligned}10f(n)&=10q+10m\,r\\&=n+(10m-1)\,r\\&=n+k\,r\,p.\end{aligned}$$

So $10f(n)-n$ is a multiple of $p$, and $p$ divides $n$ exactly when it divides $10f(n)$. As $p$ has no common factor with $10$, it divides $10f(n)$ exactly when it divides $f(n)$, by the general form of Euclid's lemma. ^[3] $\square$

Both examples of the statement have $r=5$, and $k=3$ for $113$, so $10f(n)-n=15\cdot113=1695$: $77970-76275$${}=14040-12345$${}=1695$. The remainder is not kept: multiplying $10f(n)\equiv n$ by $m$ gives $f(n)\equiv m\,n\pmod p$, so $12345$, which leaves $28$ after division by $113$, goes to $1404$, which leaves $48$, the remainder of $34\cdot28=952$.

5. No other multiplier works

Claim. Every divisibility multiplier $m'$ for $p$ equals the $m$ of Section 3.

Proof. Multiplying $p$ by $j=1$, $7$, $3$ or $9$ when its last digit is $1$, $3$, $7$ or $9$ gives a multiple $jp=10q+1$ that ends in $1$, since $1\cdot1$, $3\cdot7$, $7\cdot3$ and $9\cdot9$ end in $1$. Because $jp$ is a multiple of $p$, both $q+m'$ and $q+m$ are multiples of $p$, the first by the definition and the second by Section 4. So $p$ divides $m'-m$; as both lie between $1$ and $p-1$, $|m'-m|\lt p$, and therefore $m'=m$. $\square$

For $113$ the test number is $7\cdot113=791$: $79+m'$ must be a multiple of $113$, and the only one from $80$ to $191$ is $113$, so $m'=34$. So every $p\gt1$ coprime to $10$ has exactly one multiplier, the inverse of $10$ modulo $p$, which we call $m(p)$, and a wrong candidate already fails at $n=jp\le9p$.

6. Worked example: the primes below $50$

The primes below $50$ other than $2$ and $5$, grouped by their last digit, with $m(p)=(kp+1)/10$:

last digit $1$, $k=9$: $m(11)=10$, $m(31)=28$, $m(41)=37$;
last digit $3$, $k=3$: $m(3)=1$, $m(13)=4$, $m(23)=7$, $m(43)=13$;
last digit $7$, $k=7$: $m(7)=5$, $m(17)=12$, $m(37)=26$, $m(47)=33$;
last digit $9$, $k=1$: $m(19)=2$, $m(29)=3$.

The groups add up to $75$, $25$, $76$ and $5$, so the multipliers of the primes below $50$ sum to $181$. In the same way the primes below $100$ give $501$, and those below $1000$ give $39517$, the value of the statement.

7. Repeating $f$

The statement notes that $f(n)\lt n$ when $n$ is large, so that repeating $f$ gives a divisibility test. Precisely, if $n\ge10m$, then $q\ge m$, and since $r\le9$,

$$\begin{aligned}n-f(n)&=9q-(m-1)\,r\\&\ge9q-9(m-1)\ge9.\end{aligned}$$

So the chain $n,f(n),f(f(n)),\dots$ decreases until it falls below $10m=kp+1$. By Section 4 each of its numbers is divisible by $p$ exactly when $n$ is, and the last one is compared with $p,2p,\dots,kp$. For $113$, where $10m=340$:

$76275\to7797\to1017\to339$, and $339=3\cdot113$, so $113$ divides $76275$;
$12345\to1404\to276$, and $276$ is neither $113$, $226$ nor $339$, so $113$ does not divide $12345$.

Below $10m$ the chain may stop shrinking: $f(10m-1)$${}=(m-1)+9m=10m-1$, for example $f(339)=339$.

8. Negative multipliers and school rules

School rules often subtract: "subtract twice the last digit from the rest" for $7$ uses $q-2r$ with $-2=m-p$, and works because $q+(m-p)\,r=f(n)-p\,r$. Likewise $m-p=-1$ for $11$ and $-5$ for $17$, while $m=4$ for $13$ and $m=2$ for $19$ are rules that add. ^[4] For $3$, $m=1$ because $10\equiv1\pmod3$, so every power of $10$ is congruent to $1$ and $n$ is congruent to its digit sum.

9. The extended Euclidean algorithm

The C++ and Java programs do not use the formula of Section 3; they compute the same inverse with the extended Euclidean algorithm. ^[6] By hand for $113$: $113=11\cdot10+3$ and $10=3\cdot3+1$, so

$$\begin{aligned}1&=10-3\cdot3\\&=10-3\,(113-11\cdot10)\\&=34\cdot10-3\cdot113,\end{aligned}$$

and $34\cdot10\equiv1\pmod{113}$. The algorithm does this bookkeeping forwards. It keeps two numbers, starting with $a=10$ and $b=p$, and with each a coefficient $x$ such that the number is congruent to $10x$ modulo $p$: at first $x=1$ for $a$ and $x=0$ for $b$. A round divides $a$ by $b$ with quotient $t$ and replaces $a,b$ by $b,a-tb$, and the coefficients in the same way, so the congruences survive. The pair follows the ordinary Euclidean algorithm and ends with $b=0$ and $a=\gcd(10,p)=1$. Then the coefficient $x$ of $a$ satisfies $10x\equiv1\pmod p$, and its remainder after division by $p$ is $m(p)$.

Since one number is $10$, the loop is short. For $p\gt10$ the first round only swaps the numbers, the second leaves $10$ and the last digit $d$ of $p$, and from $(10,d)$ at most three rounds remain, the longest being $(10,7)$, $(7,3)$, $(3,1)$, $(1,0)$. So no prime needs more than five rounds; $113$ needs four.

10. Listing the primes: the linear sieve

The primes come from the linear sieve, ^[7] a variant of the sieve of Eratosthenes that crosses out each composite number exactly once. It runs through $i=2,3,\dots,n$. An $i$ that is not crossed out is prime and joins the list. Then the sieve crosses out $i\cdot s$ for the listed primes $s=2,3,5,\dots$ in increasing order, and stops after the smallest prime factor of $i$ or when $i\cdot s\gt n$: for $i=12$ it crosses out only $24$, for $i=25$ the numbers $50$, $75$ and $125$.

A composite $c\le n$ with smallest prime factor $s$ is crossed out at the earlier step $i=c/s$, as every prime factor of this $i$ is at least $s$; so an $i$ that is not crossed out is indeed prime. Any other factorisation $c=i'\cdot s'$ with a prime $s'\gt s$ has $s$ dividing $i'$, so the loop for $i'$ stops before $s'$; thus $36$ is crossed out at $i=18$, not at $i=12$. One crossing per composite, plus at most one failed bound check for each $i$, makes the time proportional to $n$.

How the Code Works

The C++ program ^[9] reads the options --prime-limit=, by default $10^7$, and --skip-checkpoints, and exits with code $1$ on an unknown argument or a limit below $3$. Its digit reader parse_int_after_prefix accumulates in a signed int without an overflow check, so an oversized limit is undefined behaviour; a test build read --prime-limit=4294968296 as $1000$. Unless skipped, run_checkpoints requires inverse_mod_10(113) to be $34$ and solve(1000) to be $39517$, or the program exits with code $2$; then main prints solve for the limit.

solve calls primes_up_to(prime_limit - 1), so exactly the primes below the limit are used, skips $2$ and $5$, and adds inverse_mod_10(p) to the unsigned 64-bit sum. primes_up_to is the sieve of Section 10 on a std::vector<std::uint8_t>, with $i\cdot s$ formed as a long long. inverse_mod_10 is the loop of Section 9, with the numbers a, b, their coefficients x0, x1 and the quotient q; for $113$, x0 runs through $1,0,1,-11,34$. C++ rounds quotients towards zero, so x0 % p can be negative, and then $p$ is added: for $p=7$ the loop ends with $-2$, which becomes $5$. All numbers in the loop stay below $10^7$ in absolute value, so int suffices. The $664\,579$ primes below $10^7$ ^[10] have multipliers below $10^7$, so the sum is below $10^{13}$, far below $2^{63}$.

Python and Java have the same three functions, in Java primesUpTo, inverseMod10 and solve, and no options or checkpoints; Python's solve takes the limit as a parameter with default $10^7$, Java's fixes it. Python sieves a bytearray, takes the inverse from pow(10, -1, p), available since Python 3.8, ^[8] and returns a string. Java copies the C++ loop, including the correction, since its % also takes the sign of the dividend; it stores the primes in an ArrayList<Integer> and adds in a long.

Complexity and Verification

For the limit $L=10^7$ the sieve runs up to $n=L-1$, crosses out each of the $9\,335\,419$ composites once, and its inner loop runs $16\,034\,518$ times, about $1.6$ times per number. The inverses take at most five rounds per prime. So the method takes $O(L)$ time and $O(L)$ memory: about 10 MB of flags, one byte per number, and about 2.7 MB for the primes as 32-bit integers in C++, 4.2 MB with the spare capacity of the vector, or about 25 MB for a Python list with its integer objects. On an Apple M4 Max a full run takes about $0.02$ s in C++ (clang -O2), $0.1$ s in Java and $0.8$ s in Python.

Only the C++ program checks itself, with the two values of the statement. For this article, a brute force over the definition tried every $m\lt p$ on every $n\le100p$, for each $p$ from $3$ to $59$ coprime to $10$, and kept exactly the inverse of $10$; $n\le10p$ already suffices, as Section 5 predicts. For every limit from $3$ to $3000$, the C++ program agreed with a separate sum that finds the primes by trial division. For all $664\,577$ primes below $10^7$ other than $2$ and $5$, the formula of Section 3, Python's pow and a copy of inverse_mod_10 agree, the copy never needing more than five rounds; at the real limit the three programs and a separate sieve of Eratosthenes with that formula print the same number. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts; the two claims about the multiplier, the bound of Section 7 and the operation counts are derived in this article.

Project Euler 274 — Divisibility Multipliers. The official statement defines the multiplier m < p of an integer p > 1 coprime to 10 through f(n) = (all but the last digit of n) + (the last digit of n) · m, gives m = 34 for 113 and the sum 39517 below 1000, and asks for the sum for the primes below 10⁷.
Wikipedia — Modular multiplicative inverse. Two integers are congruent modulo m if m divides their difference; an inverse of a modulo m exists exactly when gcd(a, m) = 1 and is then unique modulo m.
Wikipedia — Euclid's lemma. In its general form: if an integer n divides ab and is coprime to a, then n divides b.
Wikipedia — Divisibility rule. Last-digit rules for 7, 11, 13, 17 and 19, the digit sum for 3, and a generalised rule: multiply a divisor ending in 1, 3, 7 or 9 by 9, 3, 7 or 1, add 1 and divide by 10 to get the multiplier.
OEIS A357913 — Inverse of 10 modulo prime(n). The terms 5, 10, 4, 12, 2, … start at prime(4) = 7; a(n) = (A114013(n) + 1)/10, with A114013(n) the least multiple of prime(n) ending in 9.
Wikipedia — Extended Euclidean algorithm. Remainders and Bézout coefficients are updated with the same quotient until the remainder is 0; for coprime a and n the coefficient of a, reduced modulo n, is the inverse of a.
cp-algorithms — Linear Sieve. Crossing out i · p only for the primes p up to the smallest prime factor of i reaches every composite exactly once, so the sieve runs in linear time.
Python — built-in function pow. With mod present, a negative exp is allowed for a base coprime to mod and uses an inverse of base modulo mod (since Python 3.8).
C++ — Euler274.cpp. The immutable revision with primes_up_to, inverse_mod_10, solve and run_checkpoints, the implementation record for this article.
OEIS A006880 — Number of primes < 10^n. Its term for n = 7 is 664579.

Problem 274 source code

C++

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

namespace {

using u64 = std::uint64_t;

struct Options {
    int prime_limit = 10000000;
    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, "--prime-limit=", options.prime_limit)) {
            continue;
        }

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

    return options.prime_limit >= 3;
}

std::vector<int> primes_up_to(const int n) {
    std::vector<std::uint8_t> is_composite(static_cast<std::size_t>(n + 1), 0);
    std::vector<int> primes;

    for (int i = 2; i <= n; ++i) {
        if (!is_composite[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
        }
        for (int p : primes) {
            const long long v = 1LL * i * p;
            if (v > n) {
                break;
            }
            is_composite[static_cast<std::size_t>(v)] = 1;
            if (i % p == 0) {
                break;
            }
        }
    }

    return primes;
}

int inverse_mod_10(const int p) {
    int a = 10;
    int b = p;
    int x0 = 1;
    int x1 = 0;

    while (b != 0) {
        const int q = a / b;
        const int na = b;
        const int nb = a - q * b;
        a = na;
        b = nb;

        const int nx = x1;
        const int ny = x0 - q * x1;
        x0 = nx;
        x1 = ny;
    }

    int inv = x0 % p;
    if (inv < 0) {
        inv += p;
    }
    return inv;
}

u64 solve(const int prime_limit) {
    const std::vector<int> primes = primes_up_to(prime_limit - 1);

    u64 sum = 0;
    for (int p : primes) {
        if (p == 2 || p == 5) {
            continue;
        }
        sum += static_cast<u64>(inverse_mod_10(p));
    }
    return sum;
}

bool run_checkpoints() {
    if (inverse_mod_10(113) != 34) {
        std::cerr << "Checkpoint failed for p=113" << '\n';
        return false;
    }
    if (solve(1000) != 39517ULL) {
        std::cerr << "Checkpoint failed for prime limit 1000" << '\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.prime_limit) << '\n';
    return 0;
}

Python

def primes_up_to(n):
    is_composite = bytearray(n + 1)
    primes = []
    
    for i in range(2, n + 1):
        if not is_composite[i]:
            primes.append(i)
        
        for p in primes:
            v = i * p
            if v > n:
                break
            is_composite[v] = 1
            if i % p == 0:
                break
                
    return primes

def inverse_mod_10(p):
    return pow(10, -1, p)

def solve(prime_limit=10000000):
    primes = primes_up_to(prime_limit - 1)
    
    total_sum = 0
    for p in primes:
        if p == 2 or p == 5:
            continue
        total_sum += inverse_mod_10(p)
        
    return str(total_sum)

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

Java

import java.util.*;

public class Euler274 {
    static List<Integer> primesUpTo(int n) {
        byte[] isComposite = new byte[n + 1];
        List<Integer> primes = new ArrayList<>();

        for (int i = 2; i <= n; ++i) {
            if (isComposite[i] == 0) {
                primes.add(i);
            }
            for (int p : primes) {
                long v = (long) i * p;
                if (v > n) {
                    break;
                }
                isComposite[(int) v] = 1;
                if (i % p == 0) {
                    break;
                }
            }
        }

        return primes;
    }

    static int inverseMod10(int p) {
        int a = 10;
        int b = p;
        int x0 = 1;
        int x1 = 0;

        while (b != 0) {
            int q = a / b;
            int na = b;
            int nb = a - q * b;
            a = na;
            b = nb;

            int nx = x1;
            int ny = x0 - q * x1;
            x0 = nx;
            x1 = ny;
        }

        int inv = x0 % p;
        if (inv < 0) {
            inv += p;
        }
        return inv;
    }

    public static String solve() {
        int primeLimit = 10000000;
        List<Integer> primes = primesUpTo(primeLimit - 1);

        long sum = 0;
        for (int p : primes) {
            if (p == 2 || p == 5) {
                continue;
            }
            sum += inverseMod10(p);
        }

        return String.valueOf(sum);
    }

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