Problem 487: Sums of Power Sums

Project Euler Problem 487 Solution

EulerSolve provides an optimized solution for Project Euler Problem 487, Sums of Power Sums, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For an integer $k$, define the ordinary power sum and its cumulative sum by $$F_k(n)=\sum_{j=1}^{n} j^k,\qquad S_k(n)=\sum_{i=1}^{n}F_k(i).$$ The concrete task is to evaluate $$\sum_{p\in \mathcal{P}} \left(S_{10000}(10^{12}) \bmod p\right),\qquad \mathcal{P}=\{p\text{ prime}: 2\cdot 10^9 \le p \le 2\cdot 10^9+2000\}.$$ A direct computation is hopeless: $n=10^{12}$ is far too large for nested summation, and even evaluating $F_k(n)$ term by term would be impossible. The solution therefore transforms the double sum into polynomial evaluation modulo each prime. Mathematical Approach The implementations solve the problem prime by prime. For each prime $p$ in the interval, the goal is to obtain $S_k(n)\bmod p$ without summing up to $n$. Step 1: Collapse the Double Sum Start from the definition $$S_k(n)=\sum_{i=1}^{n}\sum_{j=1}^{i} j^k.$$ Swap the order of summation. A fixed value $j$ appears in every inner sum with $i\ge j$, so it is counted exactly $n+1-j$ times. Hence $$S_k(n)=\sum_{j=1}^{n}(n+1-j)j^k.$$ Now separate the two terms: $$S_k(n)=(n+1)\sum_{j=1}^{n}j^k-\sum_{j=1}^{n}j^{k+1}=(n+1)F_k(n)-F_{k+1}(n).$$ So the whole problem reduces to evaluating two ordinary power sums modulo $p$: $F_k(n)$ and $F_{k+1}(n)$....

Detailed mathematical approach

Problem Summary

For an integer $k$, define the ordinary power sum and its cumulative sum by

$$F_k(n)=\sum_{j=1}^{n} j^k,\qquad S_k(n)=\sum_{i=1}^{n}F_k(i).$$

The concrete task is to evaluate

$$\sum_{p\in \mathcal{P}} \left(S_{10000}(10^{12}) \bmod p\right),\qquad \mathcal{P}=\{p\text{ prime}: 2\cdot 10^9 \le p \le 2\cdot 10^9+2000\}.$$

A direct computation is hopeless: $n=10^{12}$ is far too large for nested summation, and even evaluating $F_k(n)$ term by term would be impossible. The solution therefore transforms the double sum into polynomial evaluation modulo each prime.

Mathematical Approach

The implementations solve the problem prime by prime. For each prime $p$ in the interval, the goal is to obtain $S_k(n)\bmod p$ without summing up to $n$.

Step 1: Collapse the Double Sum

Start from the definition

$$S_k(n)=\sum_{i=1}^{n}\sum_{j=1}^{i} j^k.$$

Swap the order of summation. A fixed value $j$ appears in every inner sum with $i\ge j$, so it is counted exactly $n+1-j$ times. Hence

$$S_k(n)=\sum_{j=1}^{n}(n+1-j)j^k.$$

Now separate the two terms:

$$S_k(n)=(n+1)\sum_{j=1}^{n}j^k-\sum_{j=1}^{n}j^{k+1}=(n+1)F_k(n)-F_{k+1}(n).$$

So the whole problem reduces to evaluating two ordinary power sums modulo $p$: $F_k(n)$ and $F_{k+1}(n)$.

Step 2: View $F_e(x)$ as a Polynomial

For a fixed exponent $e$, the function

$$F_e(x)=\sum_{t=1}^{x} t^e$$

is a polynomial in $x$ of degree $e+1$. One way to see this is through the forward difference

$$F_e(x)-F_e(x-1)=x^e,$$

whose right-hand side has degree $e$, so $F_e$ must have degree $e+1$.

That means $F_e(x)$ is uniquely determined by its values at any $e+2$ distinct points. Over the field $\mathbb{F}_p$, the implementations use the points

$$x=0,1,2,\dots,e+1.$$

This is valid because every prime in the target interval is much larger than $e+1$, so these nodes remain distinct modulo $p$.

Step 3: Build the Sample Values

Let $m=e+1$. Define

$$y_i=F_e(i)\bmod p\qquad (0\le i\le m).$$

These values are easy to generate iteratively:

$$y_0=0,\qquad y_i=y_{i-1}+i^e \pmod p.$$

So for each prime $p$, the implementation computes $i^e\bmod p$ for $i=1,2,\dots,m$, accumulates those values, and obtains the full interpolation table $(0,y_0),(1,y_1),\dots,(m,y_m)$.

No Bernoulli numbers or closed forms are needed. The only ingredients are modular exponentiation and cumulative addition.

Step 4: Use Lagrange Interpolation over $\mathbb{F}_p$

Given the samples $y_0,\dots,y_m$, Lagrange interpolation reconstructs the value at any $x\in\mathbb{F}_p$:

$$F_e(x)=\sum_{i=0}^{m} y_i \prod_{\substack{0\le j\le m\\j\ne i}}\frac{x-j}{i-j}\pmod p.$$

The denominator can be simplified explicitly:

$$\prod_{\substack{0\le j\le m\\j\ne i}}(i-j)=i!\,(-1)^{m-i}(m-i)!.$$

This identity is what makes the implementation fast. Once factorials and inverse factorials are available modulo $p$, each basis denominator is obtained in constant time.

Step 5: Evaluate All Numerators in Linear Time

A naive interpolation would recompute

$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)$$

from scratch for every $i$, costing $O(m^2)$. The implementations avoid that by storing prefix and suffix products:

$$P_r=\prod_{j=0}^{r-1}(x-j),\qquad Q_r=\prod_{j=r}^{m}(x-j).$$

Then the numerator for node $i$ is simply

$$\prod_{\substack{0\le j\le m\\j\ne i}}(x-j)=P_i\,Q_{i+1}.$$

After one forward pass and one backward pass, every interpolation term is available in $O(1)$, so evaluating $F_e(x)$ costs only $O(m)$.

Step 6: Assemble $S_k(n)\bmod p$ and Sum over Primes

For each prime $p$, the implementation evaluates the two power sums at

$$x=n\bmod p,$$

which is sufficient because the polynomial is being evaluated in $\mathbb{F}_p$. It obtains

$$F_k(n)\bmod p\qquad \text{and}\qquad F_{k+1}(n)\bmod p,$$

then applies

$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p.$$

Primes in the interval are detected with deterministic Miller-Rabin for 64-bit integers, preceded by a few small trial divisions. Each prime contributes one residue, and the final answer is the ordinary sum of those residues.

Worked Example: $k=2$, $n=10$, $p=17$

Here we want $S_2(10)\bmod 17$. First use the identity

$$S_2(10)\equiv 11\,F_2(10)-F_3(10)\pmod{17}.$$

For $F_2$, the degree is $3$, so the samples at $x=0,1,2,3$ are

$$0,\ 1,\ 1+2^2=5,\ 1+2^2+3^2=14.$$

Thus the interpolation data modulo $17$ is

$$y^{(2)}=(0,1,5,14).$$

For $F_3$, the degree is $4$, so the samples at $x=0,1,2,3,4$ are

$$0,\ 1,\ 9,\ 36,\ 100 \equiv 0,\ 1,\ 9,\ 2,\ 15 \pmod{17}.$$

Hence

$$y^{(3)}=(0,1,9,2,15).$$

Lagrange interpolation at $x=10$ gives

$$F_2(10)\equiv 11 \pmod{17},\qquad F_3(10)\equiv 16 \pmod{17}.$$

Therefore

$$S_2(10)\equiv 11\cdot 11-16\equiv 3 \pmod{17}.$$

A direct check confirms this: $S_2(10)=1210$, and $1210\equiv 3\pmod{17}$.

How the Code Works

The C++, Python, and Java implementations all follow the same structure. They scan the interval from $2\cdot 10^9$ to $2\cdot 10^9+2000$, keep only primes, and process each surviving modulus independently.

For one prime $p$, the implementation first builds factorials and inverse factorials up to $k+2$. Only one modular inverse is needed: after inverting the largest factorial, the remaining inverse factorials are recovered by a backward pass. Because every target prime is far larger than $k+2$, all factorial values involved are invertible modulo $p$.

Next, it evaluates the degree-$k+1$ and degree-$k+2$ power-sum polynomials. Each sample table is generated by modular exponentiation plus prefix summation. If $n\bmod p$ already lies among the interpolation nodes, the value is returned immediately from the table; otherwise the prefix/suffix Lagrange formula is used.

Finally, the implementation combines the two polynomial values with

$$S_k(n)\equiv (n+1)F_k(n)-F_{k+1}(n)\pmod p,$$

and adds that residue to the running total. No global modulus is applied to the final accumulation.

Complexity Analysis

Let $k$ be the exponent parameter and let $\pi$ be the number of primes in the target interval. For one prime $p$, building factorial and inverse-factorial tables costs $O(k)$ modular operations and $O(k)$ memory. Constructing the sample values for $F_k$ and $F_{k+1}$ requires $O(k)$ modular exponentiation calls; counting modular multiplications explicitly, this is $O(k\log k)$. The two Lagrange evaluations themselves are linear, so their additional cost is $O(k)$.

Therefore the overall running time is $O(\pi k\log k)$ modular multiplications, with $O(k)$ memory. Since the prime interval is short, the per-prime interpolation work dominates the total runtime.

Footnotes and References

Problem page: https://projecteuler.net/problem=487
Power sums and Faulhaber-type identities: Wikipedia — Faulhaber's formula
Polynomial interpolation: Wikipedia — Lagrange polynomial
Primality testing: Wikipedia — Miller-Rabin primality test

Problem 487 source code

C++

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

namespace {

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

struct Options {
    int k = 10'000;
    u64 n = 1'000'000'000'000ULL;
    int p_lo = 2'000'000'000;
    int p_hi = 2'000'002'000;
    bool run_checkpoints = true;
};

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

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& out) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        out = static_cast<u64>(std::stoull(tail));
    } catch (...) {
        return false;
    }
    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, "--k=", options.k)) {
            continue;
        }
        if (parse_u64_after_prefix(arg, "--n=", options.n)) {
            continue;
        }
        if (parse_int_after_prefix(arg, "--p-lo=", options.p_lo)) {
            continue;
        }
        if (parse_int_after_prefix(arg, "--p-hi=", options.p_hi)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    if (options.k < 0 || options.p_lo > options.p_hi || options.p_lo < 2) {
        std::cerr << "Invalid arguments.\n";
        return false;
    }
    return true;
}

u64 mod_pow(u64 base, u64 exp, const u64 mod) {
    u64 out = 1ULL;
    base %= mod;
    while (exp > 0ULL) {
        if (exp & 1ULL) {
            out = static_cast<u64>((static_cast<u128>(out) * base) % mod);
        }
        base = static_cast<u64>((static_cast<u128>(base) * base) % mod);
        exp >>= 1ULL;
    }
    return out;
}

u64 mod_inv(const u64 x, const u64 p) { return mod_pow(x, p - 2ULL, p); }

u64 add_mod(const u64 a, const u64 b, const u64 p) {
    const u64 c = a + b;
    return (c >= p || c < a) ? (c % p) : c;
}

u64 sub_mod(const u64 a, const u64 b, const u64 p) {
    return (a >= b) ? (a - b) : (a + p - b);
}

bool is_prime_64(const u64 n) {
    if (n < 2ULL) {
        return false;
    }
    for (const u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL,
                        31ULL, 37ULL}) {
        if (n == p) {
            return true;
        }
        if (n % p == 0ULL) {
            return false;
        }
    }

    u64 d = n - 1ULL;
    int s = 0;
    while ((d & 1ULL) == 0ULL) {
        d >>= 1ULL;
        ++s;
    }

    auto witness = [&](const u64 a) {
        if (a % n == 0ULL) {
            return false;
        }
        u64 x = mod_pow(a, d, n);
        if (x == 1ULL || x == n - 1ULL) {
            return false;
        }
        for (int r = 1; r < s; ++r) {
            x = static_cast<u64>((static_cast<u128>(x) * x) % n);
            if (x == n - 1ULL) {
                return false;
            }
        }
        return true;
    };

    // Deterministic for 64-bit with this base set.
    for (const u64 a : {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL}) {
        if (witness(a)) {
            return false;
        }
    }
    return true;
}

u64 lagrange_from_0_to_m(const std::vector<u64>& y,
                         const int m,
                         const u64 x,
                         const std::vector<u64>& fac,
                         const std::vector<u64>& invfac,
                         const u64 p) {
    if (x <= static_cast<u64>(m)) {
        return y[static_cast<std::size_t>(x)];
    }

    std::vector<u64> pref(static_cast<std::size_t>(m + 2), 1ULL);
    std::vector<u64> suf(static_cast<std::size_t>(m + 2), 1ULL);

    for (int i = 0; i <= m; ++i) {
        pref[static_cast<std::size_t>(i + 1)] =
            static_cast<u64>((static_cast<u128>(pref[static_cast<std::size_t>(i)]) *
                              sub_mod(x, static_cast<u64>(i), p)) %
                             p);
    }
    for (int i = m; i >= 0; --i) {
        suf[static_cast<std::size_t>(i)] =
            static_cast<u64>((static_cast<u128>(suf[static_cast<std::size_t>(i + 1)]) *
                              sub_mod(x, static_cast<u64>(i), p)) %
                             p);
    }

    u64 out = 0ULL;
    for (int i = 0; i <= m; ++i) {
        u64 num =
            static_cast<u64>((static_cast<u128>(pref[static_cast<std::size_t>(i)]) *
                              suf[static_cast<std::size_t>(i + 1)]) %
                             p);
        u64 den =
            static_cast<u64>((static_cast<u128>(invfac[static_cast<std::size_t>(i)]) *
                              invfac[static_cast<std::size_t>(m - i)]) %
                             p);
        if (((m - i) & 1) != 0) {
            den = (den == 0ULL) ? 0ULL : (p - den);
        }
        u64 term =
            static_cast<u64>((static_cast<u128>(y[static_cast<std::size_t>(i)]) * num) % p);
        term = static_cast<u64>((static_cast<u128>(term) * den) % p);
        out = add_mod(out, term, p);
    }
    return out;
}

u64 power_sum_mod(const int exp, const u64 n, const u64 p, const std::vector<u64>& fac,
                  const std::vector<u64>& invfac) {
    const int m = exp + 1;
    std::vector<u64> y(static_cast<std::size_t>(m + 1), 0ULL);
    for (int i = 1; i <= m; ++i) {
        const u64 pw = mod_pow(static_cast<u64>(i), static_cast<u64>(exp), p);
        y[static_cast<std::size_t>(i)] = add_mod(y[static_cast<std::size_t>(i - 1)], pw, p);
    }
    return lagrange_from_0_to_m(y, m, n % p, fac, invfac, p);
}

u64 S_mod_prime(const int k, const u64 n, const u64 p) {
    const int max_m = k + 2;
    std::vector<u64> fac(static_cast<std::size_t>(max_m + 1), 1ULL);
    for (int i = 1; i <= max_m; ++i) {
        fac[static_cast<std::size_t>(i)] =
            static_cast<u64>((static_cast<u128>(fac[static_cast<std::size_t>(i - 1)]) * i) % p);
    }
    std::vector<u64> invfac(static_cast<std::size_t>(max_m + 1), 1ULL);
    invfac[static_cast<std::size_t>(max_m)] = mod_inv(fac[static_cast<std::size_t>(max_m)], p);
    for (int i = max_m; i >= 1; --i) {
        invfac[static_cast<std::size_t>(i - 1)] =
            static_cast<u64>((static_cast<u128>(invfac[static_cast<std::size_t>(i)] * i)) % p);
    }

    const u64 fk = power_sum_mod(k, n, p, fac, invfac);
    const u64 fk1 = power_sum_mod(k + 1, n, p, fac, invfac);
    const u64 n1 = (n + 1ULL) % p;
    const u64 part = static_cast<u64>((static_cast<u128>(n1) * fk) % p);
    return sub_mod(part, fk1, p);
}

u64 solve(const int k, const u64 n, const int p_lo, const int p_hi) {
    u64 total = 0ULL;
    for (int p = p_lo; p <= p_hi; ++p) {
        if (!is_prime_64(static_cast<u64>(p))) {
            continue;
        }
        total += S_mod_prime(k, n, static_cast<u64>(p));
    }
    return total;
}

u64 direct_S_small(const int k, const int n) {
    u128 out = 0;
    for (int i = 1; i <= n; ++i) {
        u128 fi = 0;
        for (int j = 1; j <= i; ++j) {
            u128 pw = 1;
            for (int t = 0; t < k; ++t) {
                pw *= static_cast<u128>(j);
            }
            fi += pw;
        }
        out += fi;
    }
    return static_cast<u64>(out);
}

bool run_checkpoints() {
    if (direct_S_small(4, 100) != 35'375'333'830ULL) {
        std::cerr << "Checkpoint failed: S_4(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 1;
    }

    std::cout << solve(options.k, options.n, options.p_lo, options.p_hi) << '\n';
    return 0;
}

Python

import math
import sys

def mod_pow(base, exp, mod):
    return pow(base, exp, mod)

def mod_inv(x, p):
    return mod_pow(x, p - 2, p)

def add_mod(a, b, p):
    c = a + b
    return c % p if (c >= p or c < a) else c

def sub_mod(a, b, p):
    return a - b if a >= b else a + p - b

def is_prime_64(n):
    if n < 2: return False
    for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]:
        if n == p: return True
        if n % p == 0: return False
    
    d = n - 1
    s = 0
    while (d & 1) == 0:
        d >>= 1
        s += 1
        
    def witness(a):
        if a % n == 0: return False
        x = mod_pow(a, d, n)
        if x == 1 or x == n - 1: return False
        for _ in range(1, s):
            x = (x * x) % n
            if x == n - 1: return False
        return True
        
    for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
        if witness(a): return False
    return True

def lagrange_from_0_to_m(y, m, x, fac, invfac, p):
    if x <= m:
        return y[x]
        
    pref = [1] * (m + 2)
    suf = [1] * (m + 2)
    
    for i in range(m + 1):
        pref[i + 1] = (pref[i] * sub_mod(x, i, p)) % p
    for i in range(m, -1, -1):
        suf[i] = (suf[i + 1] * sub_mod(x, i, p)) % p
        
    out = 0
    for i in range(m + 1):
        num = (pref[i] * suf[i + 1]) % p
        den = (invfac[i] * invfac[m - i]) % p
        if (m - i) & 1:
            den = 0 if den == 0 else p - den
        term = (y[i] * num) % p
        term = (term * den) % p
        out = add_mod(out, term, p)
    return out

def power_sum_mod(exp, n, p, fac, invfac):
    m = exp + 1
    y = [0] * (m + 1)
    for i in range(1, m + 1):
        pw = mod_pow(i, exp, p)
        y[i] = add_mod(y[i - 1], pw, p)
    return lagrange_from_0_to_m(y, m, n % p, fac, invfac, p)

def S_mod_prime(k, n, p):
    max_m = k + 2
    fac = [1] * (max_m + 1)
    for i in range(1, max_m + 1):
        fac[i] = (fac[i - 1] * i) % p
        
    invfac = [1] * (max_m + 1)
    invfac[max_m] = mod_inv(fac[max_m], p)
    for i in range(max_m, 0, -1):
        invfac[i - 1] = (invfac[i] * i) % p
        
    fk = power_sum_mod(k, n, p, fac, invfac)
    fk1 = power_sum_mod(k + 1, n, p, fac, invfac)
    n1 = (n + 1) % p
    part = (n1 * fk) % p
    return sub_mod(part, fk1, p)

def solve_impl(k, n, p_lo, p_hi):
    total = 0
    for p in range(p_lo, p_hi + 1):
        if not is_prime_64(p):
            continue
        total += S_mod_prime(k, n, p)
    return total

def solve():
    return str(solve_impl(10000, 1000000000000, 2000000000, 2000002000))

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

Java

import java.math.BigInteger;

public class Euler487 {
    public static String solve() {
        int k = 10000;
        long n = 1000000000000L;
        long pLo = 2000000000L, pHi = 2000002000L;
        long total = 0;
        for (long p = pLo; p <= pHi; p++) {
            if (!isPrime64(p))
                continue;
            total += sModPrime(k, n, p);
        }
        return String.valueOf(total);
    }

    static long sModPrime(int k, long n, long p) {
        int maxM = k + 2;
        long[] fac = new long[maxM + 1], ifac = new long[maxM + 1];
        fac[0] = 1;
        for (int i = 1; i <= maxM; i++)
            fac[i] = fac[i - 1] * i % p;
        ifac[maxM] = modPow(fac[maxM], p - 2, p);
        for (int i = maxM; i > 0; i--)
            ifac[i - 1] = ifac[i] * i % p;
        long fk = powerSumMod(k, n, p, fac, ifac);
        long fk1 = powerSumMod(k + 1, n, p, fac, ifac);
        long n1 = (n + 1) % p;
        long part = n1 * fk % p;
        return (part - fk1 + p) % p;
    }

    static long powerSumMod(int exp, long n, long p, long[] fac, long[] ifac) {
        int m = exp + 1;
        long[] y = new long[m + 1];
        for (int i = 1; i <= m; i++) {
            long pw = modPow(i, exp, p);
            y[i] = (y[i - 1] + pw) % p;
        }
        return lagrange(y, m, n % p, fac, ifac, p);
    }

    static long lagrange(long[] y, int m, long x, long[] fac, long[] ifac, long p) {
        if (x <= m)
            return y[(int) x];
        long[] pref = new long[m + 2], suf = new long[m + 2];
        pref[0] = 1;
        suf[m + 1] = 1;
        for (int i = 0; i <= m; i++)
            pref[i + 1] = pref[i] * ((x - i + p) % p) % p;
        for (int i = m; i >= 0; i--)
            suf[i] = suf[i + 1] * ((x - i + p) % p) % p;
        long out = 0;
        for (int i = 0; i <= m; i++) {
            long num = pref[i] * suf[i + 1] % p;
            long den = ifac[i] * ifac[m - i] % p;
            if ((m - i) % 2 != 0)
                den = (p - den) % p;
            out = (out + y[i] * num % p * den % p) % p;
        }
        return out;
    }

    static long modPow(long b, long e, long p) {
        b %= p;
        if (b < 0)
            b += p;
        long r = 1;
        while (e > 0) {
            if ((e & 1) != 0)
                r = BigInteger.valueOf(r).multiply(BigInteger.valueOf(b)).mod(BigInteger.valueOf(p)).longValue();
            b = BigInteger.valueOf(b).multiply(BigInteger.valueOf(b)).mod(BigInteger.valueOf(p)).longValue();
            e >>= 1;
        }
        return r;
    }

    static boolean isPrime64(long n) {
        if (n < 2)
            return false;
        for (long p : new long[] { 2, 3, 5, 7, 11, 13 }) {
            if (n % p == 0)
                return n == p;
        }
        long d = n - 1;
        int s = 0;
        while (d % 2 == 0) {
            d /= 2;
            s++;
        }
        for (long a : new long[] { 2, 325, 9375, 28178, 450775, 9780504, 1795265022L }) {
            if (a % n == 0)
                continue;
            long x = modPow(a, d, n);
            if (x == 1 || x == n - 1)
                continue;
            boolean witness = true;
            for (int i = 1; i < s; i++) {
                x = BigInteger.valueOf(x).multiply(BigInteger.valueOf(x)).mod(BigInteger.valueOf(n)).longValue();
                if (x == n - 1) {
                    witness = false;
                    break;
                }
            }
            if (witness)
                return false;
        }
        return true;
    }

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