Problem 381: $(\text{prime}-k)$ Factorial

Project Euler Problem 381 Solution

EulerSolve provides an optimized solution for Project Euler Problem 381, $(\text{prime}-k)$ Factorial, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For every prime $p \ge 5$, define $$S(p)=\sum_{k=1}^{5}(p-k)!\pmod{p}=(p-1)!+(p-2)!+(p-3)!+(p-4)!+(p-5)!\pmod{p}.$$ The problem asks for the prime sum $$\sum_{5 \le p \lt 10^8} S(p).$$ A naive implementation would recompute factorial residues for every prime, but the local C++, Python, and Java solutions all exploit the same fact: after Wilson's theorem, the five factorial terms collapse to one constant-size modular expression. Mathematical Approach Step 1: Start from Wilson's theorem For every prime $p$, Wilson's theorem gives $$(p-1)! \equiv -1 \pmod{p}.$$ Because the problem only uses primes $p \ge 5$, the numbers $2$, $3$, $4$, $8$, and $24$ are all invertible modulo $p$. That lets us divide congruences safely. Step 2: Rewrite the descending factorials one by one Using $p-1 \equiv -1$, $p-2 \equiv -2$, $p-3 \equiv -3$, and $p-4 \equiv -4 \pmod{p}$, we obtain $$(p-2)! \equiv \frac{(p-1)!}{p-1} \equiv \frac{-1}{-1} \equiv 1 \pmod{p},$$ $$(p-3)! \equiv \frac{(p-2)!}{p-2} \equiv \frac{1}{-2} \equiv -\frac{1}{2} \pmod{p},$$ $$(p-4)! \equiv \frac{(p-3)!}{p-3} \equiv \frac{-1/2}{-3} \equiv \frac{1}{6} \pmod{p},$$ $$(p-5)! \equiv \frac{(p-4)!}{p-4} \equiv \frac{1/6}{-4} \equiv -\frac{1}{24} \pmod{p}.$$ So the five required terms are known without ever computing a factorial explicitly....

Detailed mathematical approach

Problem Summary

For a prime $p$ the problem defines

$$S(p)=\Bigl(\sum_{k=1}^{5}(p-k)!\Bigr)\bmod p,$$

the remainder modulo $p$ of the five descending factorials $(p-1)!,\ldots,(p-5)!$; it needs $p\ge5$, so that $(p-5)!$ exists. The statement works out $p=7$, where

$$\begin{aligned}&6!+5!+4!+3!+2!\\&=720+120+24+6+2=872,\end{aligned}$$

and $872\bmod7=4$, so $S(7)=4$. It also records $\sum S(p)=480$ over the primes with $5\le p\lt100$. We must evaluate the same sum over all primes below $10^8$: ^[1]

$$\sum_{5\le p\lt 10^8}S(p).$$

The whole solution in one sentence: Wilson's theorem turns each of the five factorials into a tiny fraction modulo $p$, those fractions add up to the single constant $-\tfrac38$ for every prime $p\ge5$, so the task shrinks to sieving the primes below $10^8$ and adding one closed-form value per prime.

Mathematical Approach

1. What the definition gives on small primes

It pays to evaluate the definition literally:

$p=5$: with $0!=1$, $24+6+2+1+1=34$ and $34\bmod5=4$, so $S(5)=4$;
$p=7$: $872$ and $872\bmod7=4$, so $S(7)=4$;
$p=11$: $4037760$, so $S(11)=1$;
$p=13$: $522950400$, so $S(13)=11$;
$p=17$: $S(17)=6$; $p=19$: $S(19)=2$.

The factorials grow fast — $18!$ already has sixteen digits — so forming them exactly for every prime in the target range is hopeless, and even multiplying modulo $p$ from the start costs $O(p)$ per prime. The results $4,4,1,11,6,2$ also look random, which we must explain away.

2. Wilson's theorem

Throughout, $a\equiv b\pmod p$ means that $a$ and $b$ leave the same remainder on division by $p$, and dividing modulo $p$ means multiplying by an inverse: modulo $7$ the inverse of $2$ is $4$, since $2\cdot4\equiv1$, so $-1/2\equiv-4\equiv3\pmod 7$. Every fraction below carries that meaning. The one classical fact we need concerns the largest term. Wilson's theorem: for every prime $p$,

$$(p-1)!\equiv-1\pmod p.$$

For $p=2$ this is $1!\equiv-1\pmod 2$. For odd primes it is a pairing argument. Every residue $x$ with $1\le x\le p-1$ has exactly one inverse $x^{-1}$ in that range, since $\gcd(x,p)=1$, ^[3] and is its own inverse exactly when $x^2\equiv1$, that is when $p$ divides $(x-1)(x+1)$, which for prime $p$ forces $x\equiv\pm1$. So the factors of $(p-1)!$ other than $1$ and $p-1$ split into disjoint pairs $\{x,x^{-1}\}$ with product $1$, and everything collapses to $1\cdot(p-1)\equiv-1$. For $p=11$ the pairs $\{2,6\}$, $\{3,4\}$, $\{5,9\}$, $\{7,8\}$ have products $12,12,45,56$, each one more than a multiple of $11$, leaving $10!\equiv-1$. The result goes back to Ibn al-Haytham; Lagrange proved it in 1771. ^[2]

3. Walking down from $(p-1)!$

The other four terms follow from the first by division. Since $(p-1)!=(p-1)\cdot(p-2)!$ and $p-1\equiv-1$, we may divide both sides by $-1$ and get $(p-2)!\equiv1$. Repeating with $p-2\equiv-2$, $p-3\equiv-3$, $p-4\equiv-4$ continues the chain, and the pattern is

$$(p-k)!\equiv\frac{(-1)^k}{(k-1)!}\pmod p.$$

This holds whenever $1\le k\le p$: then $(p-k)!$ exists and $(k-1)!$ is a product of factors below $p$, hence invertible.

Proof by induction on $k$. For $k=1$ it is Wilson's theorem, $(p-1)!\equiv-1/0!$. Assume it for some $k\lt p$; from $(p-k)!=(p-k)\cdot(p-k-1)!$ and $p-k\equiv-k$, dividing by $-k$ yields

$$(p-k-1)!\equiv\frac{(-1)^k}{-k\,(k-1)!}=\frac{(-1)^{k+1}}{k!},$$

which is the claim for $k+1$. $\square$

4. The five terms collapse to one constant

Putting $k=1,\ldots,5$ into Section 3 and adding gives, for every prime $p\ge5$,

$$S(p)\equiv-1+1-\frac12+\frac16-\frac1{24}\pmod p.$$

The integers cancel and the fractions combine over the denominator $24$:

$$-\frac{12}{24}+\frac{4}{24}-\frac{1}{24}=-\frac{9}{24}=-\frac38.$$

Hence

$$S(p)\equiv-\frac38\equiv-3\cdot8^{-1}\pmod p$$

for every prime $p\ge5$; indeed $p\ge5$ is exactly what makes $24$, hence $8$, invertible. The apparent randomness of Section 1 is gone: all the variation in $S(p)$ comes from how the fixed rational $-\tfrac38$ is represented modulo different primes.

5. Why $-\tfrac38$, and what the $5$ is doing

The constant is not an accident of the number $5$. A derangement is a permutation that moves every object: of the six orders of three objects only $(2,3,1)$ and $(3,1,2)$ qualify, so $D_3=2$. Inclusion–exclusion counts them, correcting all $n!$ orders alternately by the $\binom nj(n-j)!=n!/j!$ that fix $j$ chosen positions:

$$D_n=n!\sum_{j=0}^{n}\frac{(-1)^j}{j!},$$

the number of derangements of $n$ objects. ^[4] Summing the identity of Section 3 over $k=1,\ldots,K$ and putting $j=k-1$ gives the negative of exactly this alternating sum with $n=K-1$. Hence, for every prime $p\ge K$,

$$\sum_{k=1}^{K}(p-k)!\equiv-\frac{D_{K-1}}{(K-1)!}\pmod p.$$

The numbers from $D_0$ on are $1,0,1,2,9,44,265$. ^[5] Our case $K=5$ gives $-D_4/4!=-9/24=-\tfrac38$, and the identity was confirmed by computation up to $K=7$. Since $D_n/n!\to e^{-1}$, the constant $-0.375$ is the five-term truncation of the series for $-e^{-1}=-0.3678\ldots$ — though what the problem adds up are the residues of that one constant modulo the different primes.

6. From a residue to an integer

A congruence is not yet a number. Being a remainder, $S(p)$ is the unique integer in $\{0,\ldots,p-1\}$ congruent to $-\tfrac38$. We need the $x$ with $8x\equiv-3\pmod p$ and $0\le x\lt p$, that is $8x+3=mp$ for some $m$; from $3\le 8x+3\lt 8p+3$ we get $1\le m\le8$, and $mp\equiv3\pmod 8$ pins $m$ down: writing $p=2u+1$ gives $p^2-1=4u(u+1)$, which $8$ divides because $u$ or $u+1$ is even, so $p^2\equiv1\pmod8$. Thus $p$ is its own inverse modulo $8$ and $m\equiv3p^{-1}\equiv3p\pmod8$. Therefore

$$S(p)=\frac{mp-3}{8},\qquad m=(3p)\bmod 8.$$

The division is exact since $mp\equiv3p^2\equiv3\pmod8$. Only four values of $m$ occur:

$p\equiv1\pmod 8\Rightarrow m=3$, $S(p)=(3p-3)/8$;
$p\equiv3\pmod 8\Rightarrow m=1$, $S(p)=(p-3)/8$;
$p\equiv5\pmod 8\Rightarrow m=7$, $S(p)=(7p-3)/8$;
$p\equiv7\pmod 8\Rightarrow m=5$, $S(p)=(5p-3)/8$.

In particular $1\le S(p)\le p-1$: the value never vanishes, since $-\tfrac38\equiv0$ would need $p\mid3$. Dirichlet's theorem puts infinitely many primes in each class and the prime number theorem for arithmetic progressions gives each an asymptotic quarter, ^[6] so the ratios $S(p)/p$ cluster near $\tfrac18$, $\tfrac38$, $\tfrac58$, $\tfrac78$. The primes with $5\le p\lt10^8$ split by class as $1439970$, $1440543$, $1440534$, $1440406$ for $p\bmod8=1,3,5,7$ — each within $0.03\%$ of a quarter.

7. The inverse of $8$ without the Euclidean algorithm

The programs take an equivalent route, building a representative of $8^{-1}$ rather than of $-3\cdot8^{-1}$. Let $r=p\bmod 8$, $k=8-r$ and

$$\mathrm{inv8}=\frac{kp+1}{8}.$$

This is an integer: $r\in\{1,3,5,7\}$ is odd, so $r^2\equiv1\pmod 8$ and

$$kp+1\equiv(8-r)r+1\equiv1-r^2\equiv0\pmod 8.$$

It inverts $8$, because $8\cdot\mathrm{inv8}=kp+1\equiv1\pmod p$, and it is already reduced: $1\le k\le7$ gives $1\le\mathrm{inv8}\le(7p+1)/8\lt p$. This replaces the extended Euclidean algorithm and the exponentiation $8^{p-2}$ — the two standard ways to invert modulo a prime ^[3] — by a fixed handful of integer operations, which matters when it runs nearly six million times.

8. A worked example, and the official checkpoint

Take $p=7$: $r=7$, $k=1$, $\mathrm{inv8}=(7+1)/8=1$, and indeed $8\equiv1\pmod 7$. So $3\cdot\mathrm{inv8}=3$ and $S(7)=7-3=4$, as the statement says; Section 6 agrees, with $m=21\bmod8=5$ and $(5\cdot7-3)/8=4$. For $p=13$: $\mathrm{inv8}=(3\cdot13+1)/8=5$ and $8\cdot5\equiv1\pmod{13}$, so $S(13)=13-(15\bmod13)=11$, and $(7\cdot13-3)/8$ gives the same. Both agree with Section 1.

Over the $23$ primes with $5\le p\lt100$ the closed form gives

$$\begin{aligned}&4,4,1,11,6,2,14,25,19,32,\\&15,5,29,46,7,53,8,44,\\&27,49,10,33,36,\end{aligned}$$

summing to $480$ — the statement's value. A wrong constant, a flipped sign or a careless reduction to $\{0,\ldots,p-1\}$ would all fail this test.

9. All that remains is to list the primes

Each prime now costs constant work, so the computation is dominated by producing the primes below $L=10^8$. The sieve of Eratosthenes does that: leave every integer in $[2,L)$ unmarked and, for each unmarked $i$ in increasing order, mark $i^2,i^2+i,i^2+2i,\ldots$ as composite. Starting at $i^2$ is safe: for $2\le j\lt i$ the product $ij$ has a prime factor $q\le j$, and $ij\ge iq\gt q^2$, so $ij$ was already marked while $q$ was being processed. And since every composite below $L$ has a prime factor at most its square root, marking may stop once $i^2\ge L$. ^[7] That is all we need: the rest is one closed-form value per prime.

How the Code Works

The C++ file is the fullest version. parse_arguments reads --limit= and --skip-checkpoints, with parse_int_after_prefix turning decimal digits into an int, and rejects any limit below $6$; it checks neither that conversion nor the sieve's increment for overflow, so only limits well below INT_MAX are safe — the default $10^8$ is. primes_below(limit) is the sieve of Section 9 over a std::vector<bool> named is_composite, collecting the survivors into a std::vector<int>; its guard i <= (limit - 1) / i is the overflow-safe form of $i^2\lt\mathrm{limit}$, never forming $i\cdot i$.

inv8_mod_prime(p) is Section 7 verbatim: r = p & 7 takes $p\bmod 8$ from the low three bits and returns (k * p + 1) / 8. S_of_prime(p) then returns (p + p - (3 * inv8) % p) % p in the unsigned type u64. With $t=(3\cdot\mathrm{inv8})\bmod p$ we have $1\le t\le p-1$, so $2p-t$ lies in $[p+1,2p-1]$ and reduces to $p-t$, the representative of $-3\cdot8^{-1}$. The doubled $p$ is redundant — $t\lt p$ already makes $p-t$ a safe unsigned subtraction in range — but harmless.

solve(limit) sieves once, skips $2$ and $3$ with if (p < 5) continue; and accumulates into a u64. brute_S_of_prime(p) is for testing only: it builds $(p-1)!\bmod p$ by multiplying up, then walks down the five terms, each divisor's inverse found by scanning residues $x$ until $(p-k)x\equiv1$ — $O(p)$ per call, being $p-1$ multiplications plus at most $5(p-1)$ search steps. run_checkpoints verifies $S(7)=4$, compares the closed form with it for $p\in\{5,7,11,13,17,19\}$ and checks solve(100) == 480.

Python keeps the flags in a bytearray, builds the same prime list, and defines inv8_mod_prime and S_of_prime inside solve. Java is the leanest: it marks composites in a boolean[] for unmarked i up to isqrt = (int) Math.sqrt(limit), keeps no prime list, and computes two long values inv8 and s inline for each prime it meets. Both fix the limit inside solve and carry no checkpoints. The integer types differ — unsigned std::uint64_t, signed long, Python's arbitrary precision — but none needs big-integer arithmetic, since no factorial is ever formed: every per-prime value is at most $kp+1\lt8\cdot10^8$, and the running total is bounded separately below.

Complexity and Verification

The sieve runs in $O(L\log\log L)$ time, the classical bound from $\sum_{i\le\sqrt L}L/i$ over primes $i$. ^[7] For $L=10^8$ its marking writes number $242570202$, close to $L\ln\ln L\approx2.9\cdot10^8$. There are $\pi(10^8)=5761455$ primes below the limit, ^[8] of which $5761453$ are at least $5$. C++ and Python then walk their stored prime lists, $O(\pi(L))$ work; Java instead scans every integer from $5$ to $L-1$, $O(L)$ work, doing the modular arithmetic only where a flag marks a prime. Either way the sieve dominates and the method stays $O(L\log\log L)$.

Memory is $O(L)$. Bit-packed, the C++ flags hold about $12.5$ MB ^[9] and $5761455$ four-byte entries another $23$ MB, plus spare vector capacity; Python holds a $100$ MB bytearray and a list of prime objects; a runtime storing one byte per boolean, as Java's does, needs $100$ MB and no prime list. Measured on one machine, peak resident sizes were $71$, $146$ and $347$ MB and run times $0.39$ s, $0.35$ s and $5.8$ s, Python's per-object integers explaining both of its figures. A segmented sieve — base primes up to $\sqrt L$, blocks of about $\sqrt L$, each block's contributions added at once and the block then discarded, so that no full prime list is kept — brings memory to $O(\sqrt L)$ at unchanged time. ^[7]

Both accumulator widths are justified. Every term has $S(p)\ge1$, and already $\sum_{5\le p\lt10^6}S(p)=18773749932$ exceeds the signed 32-bit maximum $2147483647$, so $32$ bits cannot hold the target total; conversely $S(p)\le p-1\lt L$ bounds it by $L\cdot\pi(L)\approx5.8\cdot10^{14}$, well below $2^{63}-1\approx9.2\cdot10^{18}$. C++'s unsigned and Java's signed 64-bit accumulators are therefore both safe.

Beyond the packaged checkpoints, the following was verified for this article. The closed form agrees with the literal sum of the five exact factorials for all $301$ primes with $5\le p\lt2000$, and the three routes — modular inverse, $(mp-3)/8$ and inv8 — agree for all $17982$ primes with $5\le p\lt200000$. The descent identity of Section 3 holds for every prime $p\lt600$ and every $k\le p$, the derangement identity of Section 5 for $K\le7$ and every prime $7\le p\lt400$. The partial sums $480$ for $p\lt100$ and $38140$ for $p\lt1000$ came out the same from exact factorials and from the closed form. Finally the three programs, run to the full limit, returned identical results, reproduced by an independent Python sieve using $(mp-3)/8$.

Footnotes and References

The sources below supply the background facts; the descent identity, the collapse to $-\tfrac38$ and the closed form $(mp-3)/8$ are derived here.

Project Euler 381 — (prime−k) Factorial. Defines S(p) as the sum of (p−k)! for 1 ≤ k ≤ 5 reduced modulo p, works out S(7) = 4 from 6! + 5! + 4! + 3! + 2! = 872, gives 480 for the sum over the primes 5 ≤ p < 100 and asks for that sum over the primes 5 ≤ p < 10⁸.
Wikipedia — Wilson's theorem. States that n > 1 is prime if and only if (n − 1)! ≡ −1 (mod n), gives the pairing proof of Section 2 with p = 11, and credits the statement to Ibn al-Haytham and the proof to Lagrange in 1771.
Wikipedia — Modular multiplicative inverse. An inverse of a modulo m exists if and only if gcd(a, m) = 1 and is then unique in 0 … m − 1, which legitimises the divisions of Sections 2, 3 and 6. It names the two methods Section 7 avoids: the extended Euclidean algorithm, and raising a to the power m − 2.
Wikipedia — Derangement. Gives Dₙ = n! · Σ (−1)ʲ/j! for the permutations with no fixed point, derives it by inclusion–exclusion as in Section 5, and lists D₀ = 1 through D₆ = 265 with the limit Dₙ/n! → 1/e.
OEIS A000166 — Subfactorial or derangement numbers. Lists 1, 0, 1, 2, 9, 44, 265, 1854, … with offset 0, and states a(0) = 1 with a(n) = round(n!/e) for n > 0. Its D₄ = 9 gives D₄/4! = 9/24 = 3/8.
Wikipedia — Dirichlet's theorem on arithmetic progressions. States Dirichlet's theorem: for coprime a and d there are infinitely many primes a + nd. Its Distribution section adds the asymptotic proportion 1/φ(d) for each class coprime to d, pointing to the prime number theorem for arithmetic progressions. With d = 8 these are the two claims of Section 6.
Wikipedia — Sieve of Eratosthenes. Describes the algorithm of Section 9 and the start at i², gives the time O(n log log n) for constant-time array updates, calls the O(n) memory its real weakness, and gives the segmented variant with segments of size at least √n and space O(√n).
Wikipedia — Prime-counting function. Its table of π(x) at powers of ten gives π(10⁸) = 5,761,455, the count used above.
cppreference — std::vector<bool>. Documents it as a possibly space-efficient specialisation whose layout is implementation-defined, one potential optimisation being one bit per element instead of sizeof(bool) bytes. No packing is guaranteed; where it applies the flag payload is about 12.5 MB, not 100 MB.

Problem 381 source code

C++

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

namespace {

using u64 = std::uint64_t;

struct Options {
    int limit = 100'000'000;
    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 ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(ch - '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 >= 6;
}

std::vector<int> primes_below(const int limit) {
    std::vector<bool> is_composite(static_cast<std::size_t>(limit), false);
    std::vector<int> primes;
    for (int i = 2; i < limit; ++i) {
        if (!is_composite[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
            if (i <= (limit - 1) / i) {
                for (int j = i * i; j < limit; j += i) {
                    is_composite[static_cast<std::size_t>(j)] = true;
                }
            }
        }
    }
    return primes;
}

u64 inv8_mod_prime(const u64 p) {
    const u64 r = p & 7ULL;
    const u64 k = 8ULL - r;
    return (k * p + 1ULL) / 8ULL;
}

u64 S_of_prime(const u64 p) {
    const u64 inv8 = inv8_mod_prime(p);
    return (p + p - ((3ULL * inv8) % p)) % p;
}

u64 brute_S_of_prime(const u64 p) {
    u64 fact = 1ULL;
    for (u64 i = 1; i < p; ++i) {
        fact = (fact * i) % p;
    }
    u64 term = fact;
    u64 sum = 0ULL;
    for (u64 k = 1; k <= 5; ++k) {
        sum = (sum + term) % p;
        const u64 denom = p - k;
        u64 inv = 1ULL;
        for (u64 x = 1; x < p; ++x) {
            if ((denom * x) % p == 1ULL) {
                inv = x;
                break;
            }
        }
        term = (term * inv) % p;
    }
    return sum;
}

u64 solve(const int limit) {
    const std::vector<int> primes = primes_below(limit);
    u64 sum = 0ULL;
    for (int p : primes) {
        if (p < 5) {
            continue;
        }
        sum += S_of_prime(static_cast<u64>(p));
    }
    return sum;
}

bool run_checkpoints() {
    if (S_of_prime(7ULL) != 4ULL) {
        std::cerr << "Checkpoint failed: S(7)" << '\n';
        return false;
    }
    for (u64 p : {5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL}) {
        if (S_of_prime(p) != brute_S_of_prime(p)) {
            std::cerr << "Checkpoint failed: direct factorial check for prime " << p << '\n';
            return false;
        }
    }
    if (solve(100) != 480ULL) {
        std::cerr << "Checkpoint failed: sum S(p), 5<=p<100" << '\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

def solve():
    limit = 100_000_000

    # Sieve primes
    is_composite = bytearray(limit)
    primes = []
    for i in range(2, limit):
        if not is_composite[i]:
            primes.append(i)
            if i <= (limit - 1) // i:
                for j in range(i*i, limit, i):
                    is_composite[j] = 1

    def inv8_mod_prime(p):
        r = p & 7
        k = 8 - r
        return (k * p + 1) // 8

    def S_of_prime(p):
        inv8 = inv8_mod_prime(p)
        return (p + p - (3 * inv8) % p) % p

    total = 0
    for p in primes:
        if p < 5:
            continue
        total += S_of_prime(p)

    return str(total)

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

Java

public class Euler381 {

    static String solve() {
        int limit = 100000000;
        boolean[] isComposite = new boolean[limit];

        int isqrt = (int) Math.sqrt(limit);
        for (int i = 2; i <= isqrt; i++) {
            if (!isComposite[i]) {
                for (int j = i * i; j < limit; j += i) {
                    isComposite[j] = true;
                }
            }
        }

        long total = 0;
        for (int p = 5; p < limit; p++) {
            if (!isComposite[p]) {
                long pl = p;
                long inv8 = ((8L - (pl & 7L)) * pl + 1L) / 8L;
                long s = (pl + pl - ((3L * inv8) % pl)) % pl;
                total += s;
            }
        }

        return Long.toString(total);
    }

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