Problem 541: Divisibility of Harmonic Number Denominators

Project Euler Problem 541 Solution

EulerSolve provides an optimized solution for Project Euler Problem 541, Divisibility of Harmonic Number Denominators, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $$H_n=\sum_{m=1}^{n}\frac{1}{m}$$ be the nth harmonic number. For an odd prime $p$, the important exceptional indices are those for which $H_n$ is divisible by $p$ in the $p$-adic sense. Equivalently, after reducing $H_n$ to lowest terms, its denominator is coprime to $p$ and its numerator is divisible by $p$. Denote this set by $J_p$. Once $\max J_p$ is known, the target quantity is $$M(p)=p\,\max J_p+(p-1).$$ The implementations first verify the checkpoint $M(7)=719102$, then use the same method for $p=137$. Mathematical Approach The search is performed inside the rings $\mathbb{Z}/p^r\mathbb{Z}$ rather than with huge rational numbers. That makes the divisibility by $p$ explicit and lets the algorithm lift surviving indices level by level. Step 1: Define the Exceptional Set $J_p$ Write the reduced harmonic number as $$H_n=\frac{A_n}{B_n}, \qquad \gcd(A_n,B_n)=1.$$ If $p\nmid B_n$, then $H_n$ has a well-defined residue modulo $p$. The relevant indices are exactly those for which $$H_n\in p\mathbb{Z}_p,$$ or equivalently $$p\mid A_n \qquad \text{and} \qquad p\nmid B_n.$$ These are the members of $J_p$. The denominator problem is therefore converted into a $p$-adic divisibility problem for harmonic numbers....

Detailed mathematical approach

Problem Summary

Let

$$H_n=\sum_{m=1}^{n}\frac{1}{m}$$

be the nth harmonic number. For an odd prime $p$, the important exceptional indices are those for which $H_n$ is divisible by $p$ in the $p$-adic sense. Equivalently, after reducing $H_n$ to lowest terms, its denominator is coprime to $p$ and its numerator is divisible by $p$. Denote this set by $J_p$. Once $\max J_p$ is known, the target quantity is

$$M(p)=p\,\max J_p+(p-1).$$

The implementations first verify the checkpoint $M(7)=719102$, then use the same method for $p=137$.

Mathematical Approach

The search is performed inside the rings $\mathbb{Z}/p^r\mathbb{Z}$ rather than with huge rational numbers. That makes the divisibility by $p$ explicit and lets the algorithm lift surviving indices level by level.

Step 1: Define the Exceptional Set $J_p$

Write the reduced harmonic number as

$$H_n=\frac{A_n}{B_n}, \qquad \gcd(A_n,B_n)=1.$$

If $p\nmid B_n$, then $H_n$ has a well-defined residue modulo $p$. The relevant indices are exactly those for which

$$H_n\in p\mathbb{Z}_p,$$

or equivalently

$$p\mid A_n \qquad \text{and} \qquad p\nmid B_n.$$

These are the members of $J_p$. The denominator problem is therefore converted into a $p$-adic divisibility problem for harmonic numbers.

Step 2: Decompose $H_{pn+k}$ into a Core and a Tail

For $n\ge 1$ and $0\le k<p$, split the sum up to $pn+k$ into multiples of $p$, non-multiples up to $pn$, and the final tail:

$$H_{pn+k}=\sum_{a=1}^{n}\frac{1}{ap}+\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}+\sum_{t=1}^{k}\frac{1}{pn+t}.$$

The first sum is $H_n/p$, so the exact recurrence is

$$H_{pn+k}=\frac{H_n}{p}+B_p(n)+T_{p,k}(n),$$

where

$$B_p(n)=\sum_{\substack{1\le m\le pn \\ p\nmid m}}\frac{1}{m}, \qquad T_{p,k}(n)=\sum_{t=1}^{k}\frac{1}{pn+t}.$$

This is the backbone of the lifting tree. First compute the value at $pn$, then generate $pn+1,\dots,pn+p-1$ by adding one inverse at a time.

Step 3: Represent the Block Term by an Even Polynomial

The only expensive part is $B_p(n)$. Rewriting the non-multiples of $p$ in complete blocks gives

$$B_p(n)=\sum_{a=0}^{n-1}\sum_{t=1}^{p-1}\frac{1}{ap+t}.$$

For fixed precision $p^S$, each reciprocal admits a $p$-adic expansion

$$\frac{1}{ap+t}=\frac{1}{t}\cdot\frac{1}{1+ap/t}=\frac{1}{t}\sum_{j\ge 0}(-1)^j\left(\frac{ap}{t}\right)^j \pmod{p^S}.$$

After summing over the full residue system $t=1,2,\dots,p-1$, the odd contributions cancel, so modulo $p^S$ the whole block term lies in the span of even powers of $n$. Thus there exist coefficients $c_1,\dots,c_{(S-1)/2}$ such that

$$B_p(n)\equiv \sum_{j=1}^{(S-1)/2} c_j n^{2j} \pmod{p^S}.$$

The implementations do not derive these coefficients symbolically. Instead, they evaluate $B_p(n)$ for enough small sample values and solve a short linear system modulo $p^S$.

Step 4: Lift One Level and Lose One Power of $p$

Suppose a surviving node $n$ is known together with $H_n\bmod p^r$, and suppose $H_n\equiv 0\pmod p$. Then $H_n/p$ is meaningful modulo $p^{r-1}$. Using the polynomial block correction, the first child satisfies

$$H_{pn}\equiv \frac{H_n}{p}+B_p(n) \pmod{p^{r-1}}.$$

The remaining children are then generated by cumulative updates

$$H_{pn+k}=H_{pn+k-1}+\frac{1}{pn+k} \pmod{p^{r-1}}, \qquad 1\le k<p.$$

Only children with

$$H_{pn+k}\equiv 0\pmod p$$

survive to the next layer. Each lift consumes one power of $p$, so a starting precision $S$ permits at most $S-1$ meaningful layers.

Step 5: Recover the Quantity $M(p)$

If $n\in J_p$, then $H_n/p$ is still $p$-integral, and both $B_p(n)$ and every short tail $T_{p,k}(n)$ are sums of $p$-adic units. Hence the whole block

$$H_{pn},H_{pn+1},\dots,H_{pn+p-1}$$

remains $p$-integral. Therefore the largest index attached to a surviving parent $n$ is $pn+(p-1)$. Once $\max J_p$ is known, the final answer is exactly

$$M(p)=p\,\max J_p+(p-1).$$

Worked Example: $p=7$

For $p=7$, the seed search only needs $1\le n<7$. We have

$$H_6=1+\frac12+\frac13+\frac14+\frac15+\frac16=\frac{49}{20},$$

so $H_6\equiv 0\pmod 7$. Therefore

$$J_7^{(0)}=\{6\}.$$

After one lift the surviving children are

$$J_7^{(1)}=\{42,48\}.$$

The next layers found by the implementation are

$$J_7^{(2)}=\{295,299,337,341\},$$

$$J_7^{(3)}=\{2096,2390\},$$

$$J_7^{(4)}=\{14675,16731,16735\},$$

$$J_7^{(5)}=\{102728\}.$$

No further child survives, so

$$\max J_7=102728.$$

Hence

$$M(7)=7\cdot 102728+6=719102,$$

which matches the validation used by the implementations.

How the Code Works

The C++, Python, and Java implementations follow the same sequence. First they fix the precision at $S=8$ and precompute $p,p^2,\dots,p^S$. Then they sample the block sum $B_p(n)$ at a few small values of $n$ and solve a modular linear system to recover the coefficients of the even polynomial.

Next they compute $H_1,H_2,\dots,H_{p-1}$ modulo $p^S$, keeping only the values that are congruent to $0\pmod p$. This forms the initial frontier. For each surviving node, the implementation evaluates the block polynomial, computes the child $pn$, scans the remaining $p-1$ children by successive inverse additions, and keeps only those whose harmonic values remain divisible by $p$. Throughout the search it updates the largest surviving index.

The C++ implementation additionally parallelizes frontier expansion when the frontier is large, but the arithmetic and survivor test are the same in all three languages.

Complexity Analysis

Let $L_r$ be the number of surviving nodes at lifting level $r$, and let $d=(S-1)/2$ be the number of polynomial coefficients. Recovering the block polynomial costs $O(d^3)$ modular operations after the sample sums are computed. Since $S=8$ is fixed, this precomputation is constant-size.

The main cost comes from expanding the survivor tree. Each surviving node requires one polynomial evaluation of cost $O(d)$ and then a scan of $p$ children, each using one modular inverse and one modular addition. Therefore the total search cost is

$$O\left(\sum_r L_r(d+p)\right),$$

which for fixed $S$ is effectively

$$O\left(p\sum_r L_r\right).$$

The memory usage is linear in the largest frontier:

$$O\left(\max_r L_r\right).$$

Footnotes and References

Problem page: https://projecteuler.net/problem=541
Harmonic numbers: Wikipedia - Harmonic number
p-adic numbers: Wikipedia - p-adic number
Hensel lifting background: Wikipedia - Hensel's lemma
Wolstenholme-type harmonic congruences: Wikipedia - Wolstenholme's theorem

Problem 541 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <stdexcept>
#include <thread>
#include <vector>

using namespace std;

using u128 = unsigned __int128;
using i128 = __int128_t;

static inline uint64_t mul_mod_u64(uint64_t a, uint64_t b, uint64_t mod) {
    return static_cast<uint64_t>((u128)a * (u128)b % (u128)mod);
}

static inline uint64_t sub_mod_u64(uint64_t a, uint64_t b, uint64_t mod) {
    return (a >= b) ? (a - b) : static_cast<uint64_t>((u128)mod + a - b);
}

static uint64_t modinv_u64(uint64_t a, uint64_t mod) {
    a %= mod;
    if (a == 0) throw runtime_error("modinv: a=0");

    i128 t = 0, newt = 1;
    i128 r = static_cast<i128>(mod), newr = static_cast<i128>(a);

    while (newr != 0) {
        i128 q = r / newr;
        i128 tmp_t = t - q * newt;
        t = newt;
        newt = tmp_t;
        i128 tmp_r = r - q * newr;
        r = newr;
        newr = tmp_r;
    }

    if (r != 1) throw runtime_error("modinv: not invertible");
    if (t < 0) t += static_cast<i128>(mod);
    return static_cast<uint64_t>(t);
}

static uint64_t pow_mod_u64(uint64_t a, uint64_t e, uint64_t mod) {
    uint64_t res = 1 % mod;
    a %= mod;
    while (e) {
        if (e & 1) res = mul_mod_u64(res, a, mod);
        a = mul_mod_u64(a, a, mod);
        e >>= 1;
    }
    return res;
}

static bool is_prime_trial(uint64_t p) {
    if (p < 2) return false;
    if (p % 2 == 0) return p == 2;
    for (uint64_t d = 3; d * d <= p; d += 2) {
        if (p % d == 0) return false;
    }
    return true;
}

static vector<uint64_t> pow_table(uint64_t p, int S) {
    vector<uint64_t> pp(S + 1);
    pp[0] = 1;
    for (int i = 1; i <= S; i++) pp[i] = pp[i - 1] * p;
    return pp;
}

static vector<uint64_t> compute_coeff_C(uint64_t p, int S) {
    const vector<uint64_t> pPow = pow_table(p, S);
    const uint64_t mod = pPow[S];
    const int N = (S - 1) / 2;
    if (N <= 0) return {};

    auto gcd_u64 = [](uint64_t a, uint64_t b) {
        while (b) {
            uint64_t t = a % b;
            a = b;
            b = t;
        }
        return a;
    };

    vector<uint64_t> sample_n;
    sample_n.reserve(N);
    for (uint64_t n = 1; sample_n.size() < static_cast<size_t>(N); ++n) {
        if (n % p != 0) sample_n.push_back(n);
    }

    vector<uint64_t> b(N);
    for (int i = 0; i < N; i++) {
        uint64_t n = sample_n[static_cast<size_t>(i)];
        uint64_t sum = 0;
        const uint64_t up = p * n;
        for (uint64_t m = 1; m <= up; m++) {
            if (m % p == 0) continue;
            sum += modinv_u64(m, mod);
            if (sum >= mod) sum %= mod;
        }
        b[n - 1] = sum % mod;
    }

    vector<vector<uint64_t>> M(N, vector<uint64_t>(N + 1, 0));
    for (int i = 0; i < N; i++) {
        uint64_t n = sample_n[static_cast<size_t>(i)];
        for (int k = 0; k < N; k++) {
            M[i][k] = pow_mod_u64(n, static_cast<uint64_t>(2 * (k + 1)), mod);
        }
        M[i][N] = b[i];
    }

    for (int col = 0; col < N; col++) {
        int pivot = -1;
        for (int r = col; r < N; r++) {
            if (gcd_u64(M[r][col], mod) == 1) {
                pivot = r;
                break;
            }
        }
        if (pivot == -1) {
            throw runtime_error("Coefficient solve failed: no invertible pivot.");
        }
        if (pivot != col) swap(M[pivot], M[col]);

        uint64_t invPivot = modinv_u64(M[col][col], mod);
        for (int c = col; c <= N; c++) {
            M[col][c] = mul_mod_u64(M[col][c], invPivot, mod);
        }

        for (int r = 0; r < N; r++) {
            if (r == col) continue;
            uint64_t factor = M[r][col];
            if (factor == 0) continue;
            for (int c = col; c <= N; c++) {
                uint64_t sub = mul_mod_u64(factor, M[col][c], mod);
                M[r][c] = sub_mod_u64(M[r][c], sub, mod);
            }
        }
    }

    vector<uint64_t> C(N);
    for (int i = 0; i < N; i++) C[i] = M[i][N] % mod;

    {
        uint64_t testn = sample_n.back() + 1;
        while (testn % p == 0) ++testn;
        uint64_t btest = 0;
        const uint64_t up = p * testn;
        for (uint64_t m = 1; m <= up; m++) {
            if (m % p == 0) continue;
            btest += modinv_u64(m, mod);
            if (btest >= mod) btest %= mod;
        }
        btest %= mod;

        uint64_t poly = 0;
        for (int k = 0; k < N; k++) {
            uint64_t term = mul_mod_u64(C[k],
                                        pow_mod_u64(testn, static_cast<uint64_t>(2 * (k + 1)), mod),
                                        mod);
            poly += term;
            if (poly >= mod) poly %= mod;
        }
        poly %= mod;

        if (btest != poly) {
            throw runtime_error("Coefficient validation failed: b_n mismatch.");
        }
    }

    return C;
}

struct Node {
    uint64_t n;
    uint64_t H;
};

static uint64_t max_Jp(uint64_t p, int S, bool use_threads) {
    if (p < 3 || (p % 2 == 0)) throw runtime_error("Expected an odd prime p>=3.");
    if (!is_prime_trial(p)) throw runtime_error("p is not prime.");

    const vector<uint64_t> pPow = pow_table(p, S);
    const uint64_t modS = pPow[S];
    const vector<uint64_t> C = compute_coeff_C(p, S);

    vector<uint64_t> Hsmall(p);
    Hsmall[0] = 0;
    for (uint64_t n = 1; n < p; n++) {
        Hsmall[n] = (Hsmall[n - 1] + modinv_u64(n, modS)) % modS;
    }

    vector<Node> cur;
    cur.reserve(p);
    uint64_t maxn = 0;
    for (uint64_t n = 1; n < p; n++) {
        if (Hsmall[n] % p == 0) {
            cur.push_back({n, Hsmall[n]});
            maxn = max(maxn, n);
        }
    }

    auto expand_one = [&](const Node& parent, uint64_t mod_r, uint64_t mod_next) -> vector<Node> {
        const uint64_t n = parent.n;
        const uint64_t hn = parent.H % mod_r;

        const uint64_t hn_div_p = (hn / p) % mod_next;

        uint64_t poly = 0;
        uint64_t nmod = n % mod_next;
        uint64_t n2 = mul_mod_u64(nmod, nmod, mod_next);
        uint64_t pow_n2k = n2;
        for (size_t k = 0; k < C.size(); k++) {
            uint64_t ck = C[k] % mod_next;
            if (ck) {
                poly += mul_mod_u64(ck, pow_n2k, mod_next);
                if (poly >= mod_next) poly %= mod_next;
            }
            pow_n2k = mul_mod_u64(pow_n2k, n2, mod_next);
        }
        poly %= mod_next;

        uint64_t hp_n = (hn_div_p + poly) % mod_next;

        vector<Node> children;
        children.reserve(p);

        if (hp_n % p == 0) children.push_back({p * n, hp_n});

        uint64_t h = hp_n;
        for (uint64_t k = 1; k < p; k++) {
            uint64_t denom = p * n + k;
            uint64_t invd = modinv_u64(denom, mod_next);
            h += invd;
            if (h >= mod_next) h %= mod_next;
            if (h % p == 0) children.push_back({p * n + k, h});
        }

        return children;
    };

    int r = S;
    while (!cur.empty() && r >= 2) {
        const uint64_t mod_r = pPow[r];
        const uint64_t mod_next = pPow[r - 1];

        vector<Node> nxt;

        if (!use_threads || cur.size() < 4) {
            for (const auto& node : cur) {
                vector<Node> kids = expand_one(node, mod_r, mod_next);
                for (auto& ch : kids) {
                    maxn = max(maxn, ch.n);
                    nxt.push_back(ch);
                }
            }
        } else {
            unsigned hw = thread::hardware_concurrency();
            if (hw == 0) hw = 2;
            unsigned T = min<unsigned>(hw, static_cast<unsigned>(cur.size()));

            vector<vector<Node>> parts(T);
            vector<uint64_t> localMax(T, 0);
            vector<thread> threads;
            threads.reserve(T);

            for (unsigned t = 0; t < T; t++) {
                size_t L = static_cast<size_t>(t * cur.size() / T);
                size_t R = static_cast<size_t>((t + 1) * cur.size() / T);
                threads.emplace_back([&, t, L, R]() {
                    uint64_t mx = 0;
                    vector<Node> out;
                    for (size_t i = L; i < R; i++) {
                        vector<Node> kids = expand_one(cur[i], mod_r, mod_next);
                        for (auto& ch : kids) {
                            mx = max(mx, ch.n);
                            out.push_back(ch);
                        }
                    }
                    localMax[t] = mx;
                    parts[t].swap(out);
                });
            }
            for (auto& th : threads) th.join();

            for (unsigned t = 0; t < T; t++) {
                maxn = max(maxn, localMax[t]);
                nxt.insert(nxt.end(), parts[t].begin(), parts[t].end());
            }
        }

        cur.swap(nxt);
        r--;
    }

    if (!cur.empty() && r < 2) {
        throw runtime_error("Precision S too low: increase S.");
    }

    return maxn;
}

static uint64_t compute_M(uint64_t p, int S, bool use_threads) {
    uint64_t maxJ = max_Jp(p, S, use_threads);
    return p * maxJ + (p - 1);
}

int main(int argc, char** argv) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    uint64_t p = 137;
    if (argc >= 2) {
        p = stoull(argv[1]);
    }

    const int S = 8;

    {
        uint64_t M7 = compute_M(7, S, false);
        if (M7 != 719102ULL) {
            cerr << "Validation failed: M(7) expected 719102, got " << M7 << "\n";
            return 1;
        }
    }

    uint64_t ans = compute_M(p, S, true);
    cout << ans << "\n";
    return 0;
}

Python

import math

def modinv_u64(a, mod):
    a %= mod
    if a == 0: raise ValueError("modinv: a=0")
    try:
        return pow(a, -1, mod)
    except ValueError:
        raise ValueError("modinv: not invertible")

def pow_mod_u64(a, e, mod):
    return pow(a, e, mod)

def is_prime_trial(p):
    if p < 2: return False
    if p % 2 == 0: return p == 2
    for d in range(3, math.isqrt(p) + 1, 2):
        if p % d == 0: return False
    return True

def pow_table(p, S):
    pp = [1] * (S + 1)
    for i in range(1, S + 1):
        pp[i] = pp[i - 1] * p
    return pp

def compute_coeff_C(p, S):
    pPow = pow_table(p, S)
    mod = pPow[S]
    N = (S - 1) // 2
    if N <= 0: return []

    sample_n = []
    n = 1
    while len(sample_n) < N:
        if n % p != 0:
            sample_n.append(n)
        n += 1

    b = [0] * N
    for i in range(N):
        n = sample_n[i]
        sum_val = 0
        up = p * n
        for m in range(1, up + 1):
            if m % p == 0: continue
            sum_val = (sum_val + modinv_u64(m, mod)) % mod
        b[i] = sum_val

    M = [[0] * (N + 1) for _ in range(N)]
    for i in range(N):
        n = sample_n[i]
        for k in range(N):
            M[i][k] = pow_mod_u64(n, 2 * (k + 1), mod)
        M[i][N] = b[i]

    for col in range(N):
        pivot = -1
        for r in range(col, N):
            if math.gcd(M[r][col], mod) == 1:
                pivot = r
                break
        if pivot == -1:
            raise ValueError("Coefficient solve failed")
        if pivot != col:
            M[pivot], M[col] = M[col], M[pivot]

        invPivot = modinv_u64(M[col][col], mod)
        for c in range(col, N + 1):
            M[col][c] = (M[col][c] * invPivot) % mod

        for r in range(N):
            if r == col: continue
            factor = M[r][col]
            if factor == 0: continue
            for c in range(col, N + 1):
                sub = (factor * M[col][c]) % mod
                M[r][c] = (M[r][c] - sub + mod) % mod

    C = [0] * N
    for i in range(N):
        C[i] = M[i][N] % mod

    return C

def max_Jp(p, S):
    if p < 3 or p % 2 == 0: raise ValueError("p must be odd prime")
    if not is_prime_trial(p): raise ValueError("p is not prime")

    pPow = pow_table(p, S)
    modS = pPow[S]
    C = compute_coeff_C(p, S)

    Hsmall = [0] * p
    for n in range(1, p):
        Hsmall[n] = (Hsmall[n - 1] + modinv_u64(n, modS)) % modS

    cur = []
    maxn = 0
    for n in range(1, p):
        if Hsmall[n] % p == 0:
            cur.append((n, Hsmall[n]))
            maxn = max(maxn, n)

    r = S
    while cur and r >= 2:
        mod_r = pPow[r]
        mod_next = pPow[r - 1]
        
        nxt = []
        for node in cur:
            n, hn = node
            hn_mod = hn % mod_r
            hn_div_p = (hn_mod // p) % mod_next

            poly = 0
            nmod = n % mod_next
            n2 = (nmod * nmod) % mod_next
            pow_n2k = n2
            
            for k in range(len(C)):
                ck = C[k] % mod_next
                if ck:
                    poly = (poly + ck * pow_n2k) % mod_next
                pow_n2k = (pow_n2k * n2) % mod_next

            hp_n = (hn_div_p + poly) % mod_next
            
            if hp_n % p == 0:
                nxt.append((p * n, hp_n))
                maxn = max(maxn, p * n)

            h = hp_n
            for k in range(1, p):
                denom = p * n + k
                invd = modinv_u64(denom, mod_next)
                h = (h + invd) % mod_next
                if h % p == 0:
                    nxt.append((p * n + k, h))
                    maxn = max(maxn, p * n + k)
                    
        cur = nxt
        r -= 1

    if cur and r < 2:
        raise ValueError("Precision S too low: increase S.")
    
    return maxn

def compute_M(p, S):
    maxJ = max_Jp(p, S)
    return p * maxJ + (p - 1)

def solve():
    return str(compute_M(137, 8))

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

public class Euler541 {

    static long mulMod(long a, long b, long mod) {
        if (a == 0 || b == 0)
            return 0;
        return BigInteger.valueOf(a).multiply(BigInteger.valueOf(b)).mod(BigInteger.valueOf(mod)).longValue();
    }

    static long subMod(long a, long b, long mod) {
        long res = a - b;
        if (res < 0)
            res += mod;
        return res;
    }

    static long modInv(long a, long mod) {
        long origMod = mod;
        a %= mod;
        long x = 0, y = 1, lastx = 1, lasty = 0, temp;
        long b = mod;
        while (b != 0) {
            long q = a / b;
            long r = a % b;
            a = b;
            b = r;
            temp = x;
            x = lastx - q * x;
            lastx = temp;
            temp = y;
            y = lasty - q * y;
            lasty = temp;
        }
        long res = lastx % origMod;
        if (res < 0)
            res += origMod;
        return res;
    }

    static long powMod(long a, long e, long mod) {
        long res = 1 % mod;
        a %= mod;
        while (e > 0) {
            if ((e & 1) != 0)
                res = mulMod(res, a, mod);
            a = mulMod(a, a, mod);
            e >>= 1;
        }
        return res;
    }

    static boolean isPrimeTrial(long p) {
        if (p < 2)
            return false;
        if (p % 2 == 0)
            return p == 2;
        for (long d = 3; d * d <= p; d += 2) {
            if (p % d == 0)
                return false;
        }
        return true;
    }

    static long[] powTable(long p, int S) {
        long[] pp = new long[S + 1];
        pp[0] = 1;
        for (int i = 1; i <= S; i++)
            pp[i] = pp[i - 1] * p;
        return pp;
    }

    static long gcd(long a, long b) {
        while (b != 0) {
            long t = a % b;
            a = b;
            b = t;
        }
        return a;
    }

    static long[] computeCoeffC(long p, int S) {
        long[] pPow = powTable(p, S);
        long mod = pPow[S];
        int N = (S - 1) / 2;
        if (N <= 0)
            return new long[0];

        List<Long> sampleN = new ArrayList<>();
        long n_val = 1;
        while (sampleN.size() < N) {
            if (n_val % p != 0)
                sampleN.add(n_val);
            n_val++;
        }

        long[] b = new long[N];
        for (int i = 0; i < N; i++) {
            long n = sampleN.get(i);
            long sum = 0;
            long up = p * n;
            for (long m = 1; m <= up; m++) {
                if (m % p == 0)
                    continue;
                sum += modInv(m, mod);
                if (sum >= mod)
                    sum %= mod;
            }
            b[i] = sum % mod;
        }

        long[][] M = new long[N][N + 1];
        for (int i = 0; i < N; i++) {
            long n = sampleN.get(i);
            for (int k = 0; k < N; k++) {
                M[i][k] = powMod(n, 2 * (k + 1), mod);
            }
            M[i][N] = b[i];
        }

        for (int col = 0; col < N; col++) {
            int pivot = -1;
            for (int r = col; r < N; r++) {
                if (gcd(M[r][col], mod) == 1) {
                    pivot = r;
                    break;
                }
            }
            if (pivot == -1)
                throw new RuntimeException("Coefficient solve failed");
            if (pivot != col) {
                long[] temp = M[pivot];
                M[pivot] = M[col];
                M[col] = temp;
            }

            long invPivot = modInv(M[col][col], mod);
            for (int c = col; c <= N; c++) {
                M[col][c] = mulMod(M[col][c], invPivot, mod);
            }

            for (int r = 0; r < N; r++) {
                if (r == col)
                    continue;
                long factor = M[r][col];
                if (factor == 0)
                    continue;
                for (int c = col; c <= N; c++) {
                    long sub = mulMod(factor, M[col][c], mod);
                    M[r][c] = subMod(M[r][c], sub, mod);
                }
            }
        }

        long[] C = new long[N];
        for (int i = 0; i < N; i++)
            C[i] = M[i][N] % mod;

        return C;
    }

    static class Node {
        long n;
        long H;

        Node(long n, long H) {
            this.n = n;
            this.H = H;
        }
    }

    static long maxJp(long p, int S) {
        if (p < 3 || p % 2 == 0)
            throw new RuntimeException("p must be odd prime");
        if (!isPrimeTrial(p))
            throw new RuntimeException("p is not prime");

        long[] pPow = powTable(p, S);
        long modS = pPow[S];
        long[] C = computeCoeffC(p, S);

        long[] Hsmall = new long[(int) p];
        Hsmall[0] = 0;
        for (int n = 1; n < p; n++) {
            Hsmall[n] = (Hsmall[n - 1] + modInv(n, modS)) % modS;
        }

        List<Node> cur = new ArrayList<>();
        long maxn = 0;
        for (int n = 1; n < p; n++) {
            if (Hsmall[n] % p == 0) {
                cur.add(new Node(n, Hsmall[n]));
                maxn = Math.max(maxn, n);
            }
        }

        int r = S;
        while (!cur.isEmpty() && r >= 2) {
            long mod_r = pPow[r];
            long mod_next = pPow[r - 1];

            List<Node> nxt = new ArrayList<>();
            for (Node node : cur) {
                long n = node.n;
                long hn = node.H % mod_r;
                long hn_div_p = (hn / p) % mod_next;

                long poly = 0;
                long nmod = n % mod_next;
                long n2 = mulMod(nmod, nmod, mod_next);
                long pow_n2k = n2;

                for (long ck_full : C) {
                    long ck = ck_full % mod_next;
                    if (ck != 0) {
                        poly = (poly + mulMod(ck, pow_n2k, mod_next)) % mod_next;
                    }
                    pow_n2k = mulMod(pow_n2k, n2, mod_next);
                }

                long hp_n = (hn_div_p + poly) % mod_next;

                if (hp_n % p == 0) {
                    nxt.add(new Node(p * n, hp_n));
                    maxn = Math.max(maxn, p * n);
                }

                long h = hp_n;
                for (long k = 1; k < p; k++) {
                    long denom = p * n + k;
                    long invd = modInv(denom, mod_next);
                    h = (h + invd) % mod_next;
                    if (h % p == 0) {
                        nxt.add(new Node(p * n + k, h));
                        maxn = Math.max(maxn, p * n + k);
                    }
                }
            }
            cur = nxt;
            r--;
        }

        if (!cur.isEmpty() && r < 2) {
            throw new RuntimeException("Precision S too low: increase S.");
        }

        return maxn;
    }

    static long computeM(long p, int S) {
        long maxJ = maxJp(p, S);
        return p * maxJ + (p - 1);
    }

    public static String solve() {
        return Long.toString(computeM(137, 8));
    }

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