Problem 536: Modulo Power Identity

Project Euler Problem 536 Solution

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

Problem Summary Define $S(N)$ as the sum of all integers $m\le N$ such that $m$ is squarefree and every prime divisor $p\mid m$ satisfies $$p-1\mid m+3.$$ The target value is $N=10^{12}$, so testing every integer up to $N$ is infeasible. Useful checkpoints are $$S(100)=32,\qquad S(10^6)=22868117.$$ The implementations therefore convert the prime-divisor conditions into modular constraints and search only inside branches that remain arithmetically possible. Mathematical Approach The search is organized around the factorization of $m$. Instead of scanning all integers, we explicitly choose some large prime factors, deduce a congruence for the missing part, and only enumerate that remainder when the resulting arithmetic progression is sparse enough. Step 1: Immediate consequences of the condition The definition already requires $$m\ \text{to be squarefree}.$$ Two further observations are crucial. If $m=q$ is prime, then the only nontrivial condition is $$q-1\mid q+3,$$ so $q-1\mid 4$. Hence the prime solutions are exactly $$q\in\{2,3,5\}.$$ The value $1$ is also valid vacuously because it has no prime divisors. Next, a composite solution cannot be even. If $2\mid m$ and some odd prime $p\mid m$, then $m+3$ is odd, while $p-1$ is an even number greater than $1$, so $p-1\nmid m+3$....

Detailed mathematical approach

Problem Summary

Define $S(N)$ as the sum of all integers $m\le N$ such that $m$ is squarefree and every prime divisor $p\mid m$ satisfies

$$p-1\mid m+3.$$

The target value is $N=10^{12}$, so testing every integer up to $N$ is infeasible. Useful checkpoints are

$$S(100)=32,\qquad S(10^6)=22868117.$$

The implementations therefore convert the prime-divisor conditions into modular constraints and search only inside branches that remain arithmetically possible.

Mathematical Approach

The search is organized around the factorization of $m$. Instead of scanning all integers, we explicitly choose some large prime factors, deduce a congruence for the missing part, and only enumerate that remainder when the resulting arithmetic progression is sparse enough.

Step 1: Immediate consequences of the condition

The definition already requires

$$m\ \text{to be squarefree}.$$

Two further observations are crucial.

If $m=q$ is prime, then the only nontrivial condition is

$$q-1\mid q+3,$$

so $q-1\mid 4$. Hence the prime solutions are exactly

$$q\in\{2,3,5\}.$$

The value $1$ is also valid vacuously because it has no prime divisors.

Next, a composite solution cannot be even. If $2\mid m$ and some odd prime $p\mid m$, then $m+3$ is odd, while $p-1$ is an even number greater than $1$, so $p-1\nmid m+3$. Therefore every composite solution is an odd squarefree product of distinct odd primes.

Step 2: Separate selected large prime factors from the rest

Choose a descending set of prime factors

$$\mathcal{L}=\{\ell_1\gt \ell_2\gt \cdots \gt \ell_k\ge 7\}$$

that we want to treat explicitly, and write

$$Q=\prod_{\ell\in\mathcal{L}}\ell,\qquad m=Q\,u.$$

The remaining factor $u$ must then be odd, squarefree, and have all of its prime factors strictly smaller than the smallest selected prime $\ell_k$.

For every $\ell\in\mathcal{L}$, the original condition gives

$$\ell-1\mid Q\,u+3.$$

Since $\ell\equiv 1\pmod{\ell-1}$, the selected factor $\ell$ disappears modulo $\ell-1$:

$$Q\equiv \frac{Q}{\ell}\pmod{\ell-1}.$$

So each selected prime yields one linear congruence for the unknown remainder factor:

$$u\cdot \frac{Q}{\ell}\equiv -3\pmod{\ell-1}.$$

Once the large prime factors are fixed, the rest of the problem has become a congruence problem in one variable.

Step 3: Reduce each congruence and merge them with CRT

A congruence of the form

$$a\,u\equiv -3\pmod n$$

has a solution only when

$$\gcd(a,n)\mid 3.$$

If that divisibility fails, the entire search branch is impossible and can be discarded immediately.

Otherwise, writing

$$d=\gcd(a,n),$$

we divide through by $d$ and obtain

$$\frac{a}{d}\,u\equiv -\frac{3}{d}\pmod{\frac{n}{d}}.$$

Now the coefficient is invertible modulo $n/d$, so this determines a unique residue class for $u$ modulo $n/d$.

Repeating the same reduction for every selected prime produces several residue classes, and the implementations merge them with the generalized Chinese Remainder Theorem. Either they are inconsistent, in which case the branch dies, or they collapse to a single progression

$$u\equiv r\pmod M.$$

Step 4: Enumerate the remainder only when the progression is sparse

After the congruence merge, all candidates have the form

$$u=r+tM,\qquad t\ge 0.$$

They must satisfy the obvious size bound

$$u\le \left\lfloor\frac{N}{Q}\right\rfloor.$$

There is also a sharper structural bound: because $u$ is squarefree and may use each odd prime below $\ell_k$ at most once,

$$u\le \prod_{\substack{q\lt \ell_k\\ q\ \text{odd prime}}} q.$$

Hence the real search interval is

$$u\le U=\min\left(\left\lfloor\frac{N}{Q}\right\rfloor,\ \prod_{\substack{q\lt \ell_k\\ q\ \text{odd prime}}} q\right).$$

If the arithmetic progression contributes only a modest number of terms up to $U$, the implementations enumerate those terms directly. Each candidate remainder is then checked to be odd and squarefree, factored only with primes below $\ell_k$, combined with the selected large primes, and finally tested against the original condition for every prime factor of the completed $m$.

If the progression is still too dense, the search does not enumerate yet. Instead it selects one more large prime factor, appends it to $\mathcal{L}$, and repeats the congruence-merging step. This branch-and-bound recursion shifts work from brute-force scanning into modular pruning.

Step 5: Worked example for $N=100$

The trivial valid values are

$$1,\ 2,\ 3,\ 5.$$

Now look at the branch where the only selected large prime is $7$. Then $Q=7$, and the congruence coming from $7-1=6$ is

$$u\cdot 1\equiv -3\equiv 3\pmod 6,$$

$$u\equiv 3\pmod 6.$$

Also

$$u\le \left\lfloor\frac{100}{7}\right\rfloor=14.$$

Because the remaining prime factors must be smaller than $7$, the only candidates in this progression are

$$u=3,\ 9.$$

The value $9$ is rejected because it is not squarefree. The value $3$ gives

$$m=7\cdot 3=21,$$

and indeed

$$21+3=24$$

is divisible by both $2$ and $6$, corresponding to the prime divisors $3$ and $7$.

For the next possible largest prime, $11$, the congruence becomes $u\equiv 7\pmod{10}$, so the only candidate below $100/11$ is $u=7$, which gives $m=77$. That value satisfies the congruence forced by $11$, but it fails the full condition for the smaller prime $7$, because $77+3=80$ is not divisible by $6$. No larger branch produces a valid composite below $100$.

Therefore the full list up to $100$ is

$$1,\ 2,\ 3,\ 5,\ 21,$$

and hence

$$S(100)=1+2+3+5+21=32.$$

How the Code Works

The C++, Python, and Java implementations all follow the same number-theoretic decomposition. They generate primes up to about $\sqrt{N}$, insert the trivial values $1,2,3,5$, and then recurse over descending sets of selected odd primes $\ge 7$.

For each recursive branch, the implementation turns the selected-prime constraints into linear congruences, applies the gcd solvability test, and merges the surviving residue classes with generalized CRT. If the combined modulus already exceeds the remaining numeric search interval, only the smallest nonnegative solution can still matter, so the search range is effectively clamped there.

Next it estimates how many terms of the progression remain below the current bound. A branch is enumerated directly only when that estimate is around $10^4$ terms or fewer; otherwise another large prime is appended and the recursion continues, which usually increases the modulus and decreases the candidate density.

Every enumerated remainder is factored with the branch-specific bound, checked for oddness and squarefreeness, combined with the selected primes, and then verified again against the original condition for every prime divisor. A set of accepted values removes duplicates that can arise from overlapping branches. The Python implementation is only a thin wrapper around the compiled search and parses the final numeric output.

Complexity Analysis

Let $N$ be the search limit. Building the prime table up to $\sqrt{N}$ costs $O(\sqrt{N}\log\log\sqrt{N})$ time and $O(\sqrt{N})$ memory. After that, there is no clean closed form for the total running time because the recursion is highly data-dependent: each branch either dies early due to incompatible congruences or produces an arithmetic progression whose size is roughly $U/M$ before final factoring and verification.

The implementation only enumerates a branch once that quantity is small, so most of the practical speed comes from modular pruning rather than from raw trial division. For one enumerated candidate, the remaining factor is checked by bounded trial division using primes below the current threshold, and the completed value is then validated against all of its distinct prime divisors. Memory usage is dominated by the prime list, the set of accepted values, and the recursion stack. In practice the admissible numbers are sparse enough that this strategy is fast for the target limit.

Footnotes and References

Problem page: https://projecteuler.net/problem=536
Chinese Remainder Theorem: Wikipedia — Chinese Remainder Theorem
Squarefree integer: Wikipedia — Squarefree integer
Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse
Euclidean algorithm: Wikipedia — Euclidean algorithm

Problem 536 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <unordered_set>
#include <utility>
#include <vector>
using namespace std;

using u64 = unsigned long long;
using u128 = unsigned __int128;

static u64 gcd_u64(u64 a, u64 b) {
    while (b) {
        u64 t = a % b;
        a = b;
        b = t;
    }
    return a;
}

static u64 mod_inv(u64 a, u64 mod) {
    // mod assumed > 1 and gcd(a,mod)=1.
    __int128 t = 0, newt = 1;
    __int128 r = (u64)mod, newr = (u64)a;
    while (newr != 0) {
        u64 q = (u64)(r / newr);
        __int128 tmp_t = t - (__int128)q * newt;
        t = newt;
        newt = tmp_t;
        __int128 tmp_r = r - (__int128)q * newr;
        r = newr;
        newr = tmp_r;
    }
    if (r != 1) return 0;
    t %= (__int128)mod;
    if (t < 0) t += (__int128)mod;
    return (u64)t;
}

struct CRTRes {
    u64 r;  // residue in [0, m)
    u64 m;  // modulus
    bool ok;
};

static CRTRes crt_merge(u64 r1, u64 m1, u64 r2, u64 m2, u64 maxN) {
    // Solve x ≡ r1 (mod m1), x ≡ r2 (mod m2)
    if (m1 == 1) return {r2 % m2, m2, true};
    if (m2 == 1) return {r1 % m1, m1, true};
    u64 g = gcd_u64(m1, m2);
    if ((r1 % g) != (r2 % g)) return {0, 0, false};

    u64 m1_g = m1 / g;
    u64 m2_g = m2 / g;

    // Compute t = ((r2-r1)/g) * inv(m1/g) mod (m2/g)
    long long diff = (long long)r2 - (long long)r1; // divisible by g
    long long diff_g = diff / (long long)g;

    u64 inv = mod_inv(m1_g % m2_g, m2_g);
    u64 t = (u64)(((__int128)diff_g % (long long)m2_g + (long long)m2_g) % (long long)m2_g);
    t = (u64)((u128)t * inv % m2_g);

    u128 newm128 = (u128)m1_g * (u128)m2; // m1/g * m2
    u128 newr128 = (u128)r1 + (u128)m1 * (u128)t;

    // If modulus exceeds the numeric search range, only the smallest nonnegative solution can matter.
    const u64 clamp_mod = maxN + 1;
    if (newm128 > (u128)clamp_mod) {
        u128 rr = newr128 % newm128;
        if (rr > (u128)maxN) return {0, 0, false};
        return {(u64)rr, clamp_mod, true};
    }
    u64 newm = (u64)newm128;
    u64 newr = (u64)(newr128 % newm128);
    return {newr, newm, true};
}

static vector<int> sieve_primes(int n) {
    vector<bool> is_prime(n + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (int i = 2; (long long)i * i <= n; ++i) {
        if (!is_prime[i]) continue;
        for (int j = i * i; j <= n; j += i) is_prime[j] = false;
    }
    vector<int> primes;
    for (int i = 2; i <= n; ++i) if (is_prime[i]) primes.push_back(i);
    return primes;
}

static bool factor_squarefree_bounded(u64 x, int bound_exclusive, const vector<int> &primes,
                                     vector<int> &factors_out) {
    // Factor x by trial division using primes list.
    // Return true iff x is squarefree, all prime factors are < bound_exclusive, and x is odd.
    factors_out.clear();
    if (x == 0) return false;
    if ((x & 1ULL) == 0) return false; // even not allowed (except overall m=2, handled separately)

    u64 n = x;
    for (int p : primes) {
        if ((u64)p * (u64)p > n) break;
        if (p >= bound_exclusive) break;
        if (n % (u64)p == 0) {
            factors_out.push_back(p);
            n /= (u64)p;
            if (n % (u64)p == 0) return false; // not squarefree
            while (n % (u64)p == 0) n /= (u64)p; // shouldn't happen, but keep safe
        }
    }
    if (n > 1) {
        if (n >= (u64)bound_exclusive) return false;
        // n must be prime at this point (or a product of >= bound primes, which would be caught by n>=bound if bound small)
        factors_out.push_back((int)n);
    }
    return true;
}

static bool check_solution(u64 m, const vector<int> &prime_factors) {
    // prime_factors are distinct and should be the full factorization of m.
    u64 mp3 = m + 3;
    for (int p : prime_factors) {
        int d = p - 1;
        if (d == 0) return false;
        if (mp3 % (u64)d != 0) return false;
    }
    return true;
}

struct Solver {
    u64 N;
    vector<int> primes;      // primes up to sqrt(N)+4
    vector<int> primes_small; // primes up to 1e6 (same list reused)
    unordered_set<u64> found;

    static constexpr u64 MAX_ENUM = 10000; // target max candidates per branch

    Solver(u64 N_) : N(N_) {}

    pair<u64,u64> congruence_for_remaining(const vector<int> &B, u64 P) {
        // Compute R ≡ r (mod M) implied by primes in B for remaining product R.
        u64 r = 0;
        u64 M = 1;
        for (int p : B) {
            u64 mod = (u64)(p - 1);
            if (mod == 0) return {0, 0};
            u64 a = (P / (u64)p) % mod;
            u64 g = gcd_u64(a, mod);
            if (3 % g != 0) return {0, 0};
            u64 mod2 = mod / g;
            if (mod2 == 1) {
                // no constraint
                continue;
            }
            u64 a2 = a / g;
            u64 rhs = (mod2 - ((u64)(3 / g) % mod2)) % mod2; // -3/g mod mod2
            u64 inv = mod_inv(a2 % mod2, mod2);
            if (inv == 0) return {0, 0};
            u64 rp = (u64)((u128)rhs * inv % mod2);

            CRTRes merged = crt_merge(r, M, rp, mod2, N);
            if (!merged.ok) return {0, 0};
            r = merged.r;
            M = merged.m;
        }
        return {r, M};
    }

    u64 max_remaining_product(int bound_exclusive, u64 cap) const {
        // Maximum squarefree odd product using primes < bound_exclusive, capped at cap.
        u128 prod = 1;
        for (int p : primes) {
            if (p == 2) continue;
            if (p >= bound_exclusive) break;
            if (prod > (u128)cap / (u64)p) return cap;
            prod *= (u64)p;
        }
        if (prod > (u128)cap) return cap;
        return (u64)prod;
    }

    void enumerate_and_check(const vector<int> &B, u64 P, int min_big, u64 r, u64 M) {
        u64 limit = N / P;
        limit = min(limit, max_remaining_product(min_big, limit));
        if (M == 0) return;
        // Enumerate R = r + t*M
        vector<int> factorsR;
        vector<int> factorsM;
        factorsM.reserve(B.size() + 16);

        for (u64 R = r; R <= limit; ) {
            if (R != 0) {
                if (factor_squarefree_bounded(R, min_big, primes_small, factorsR)) {
                    // Merge factors: B (largest primes) + factorsR
                    factorsM.clear();
                    for (int p : B) factorsM.push_back(p);
                    for (int p : factorsR) factorsM.push_back(p);

                    u64 m = P * R;
                    if (m <= N) {
                        // Ensure B primes are indeed >= min_big and remaining primes < min_big already checked.
                        if (check_solution(m, factorsM)) {
                            found.insert(m);
                        }
                    }
                }
            }
            // next
            if (R > limit - M) break;
            R += M;
        }
    }

    void dfs(int max_prime_idx, vector<int> &B, u64 P) {
        int min_big = B.back();
        // Compute congruence for remaining product R
        auto [r, M] = congruence_for_remaining(B, P);
        if (M == 0) return;

        u64 limit = N / P;
        limit = min(limit, max_remaining_product(min_big, limit));
        if (limit == 0) return;
        if (r > limit) {
            // There might be a larger representative r + t*M, but r is in [0,M). If r>limit and M>limit, no candidates.
            // If M <= limit, then r + t*M could still be <= limit. Handle by shifting to first >= r.
            // We'll compute using floor division.
        }

        // Quick upper bound on count: limit / M + 1
        u64 approx_cnt = limit / M + 1;
        if (approx_cnt <= MAX_ENUM) {
            enumerate_and_check(B, P, min_big, r, M);
            return;
        }

        // Need to include more large primes (next prime factor) to reduce enumeration.
        // Choose next prime q < min_big.
        for (int idx = max_prime_idx - 1; idx >= 0; --idx) {
            int q = primes[idx];
            if (q == 2) continue; // even factors not allowed in composite solutions
            u128 P2 = (u128)P * (u64)q;
            if (P2 > (u128)N) continue;
            B.push_back(q);
            dfs(idx, B, (u64)P2);
            B.pop_back();
        }
    }

    u64 solve() {
        // Handle trivial solutions.
        found.clear();
        found.insert(1);
        if (N >= 2) found.insert(2);
        if (N >= 3) found.insert(3);
        if (N >= 5) found.insert(5);

        int prime_limit = (int)(sqrt((long double)N) + 4);
        primes = sieve_primes(prime_limit);
        primes_small = primes; // same for trial division

        // Exclude 2 from being a largest prime factor in composite enumeration.
        // Iterate over odd primes >= 7 (since we already inserted 3 and 5).
        for (int i = (int)primes.size() - 1; i >= 0; --i) {
            int pmax = primes[i];
            if (pmax < 7) break;
            vector<int> B;
            B.push_back(pmax);
            dfs(i, B, (u64)pmax);
        }

        u128 sum = 0;
        for (u64 m : found) {
            if (m <= N) sum += (u128)m;
        }
        return (u64)sum;
    }
};

static u64 brute_sum(u64 N) {
    // Brute using SPF up to N (N should be small, <= 1e7 here).
    int n = (int)N;
    vector<int> spf(n + 1);
    for (int i = 0; i <= n; ++i) spf[i] = i;
    for (int i = 2; (long long)i * i <= n; ++i) {
        if (spf[i] != i) continue;
        for (long long j = 1LL * i * i; j <= n; j += i) {
            if (spf[(int)j] == (int)j) spf[(int)j] = i;
        }
    }

    u128 sum = 0;
    for (int m = 1; m <= n; ++m) {
        int x = m;
        bool ok = true;
        while (x > 1) {
            int p = spf[x];
            int cnt = 0;
            while (x % p == 0) {
                x /= p;
                if (++cnt > 1) {
                    ok = false;
                    break;
                }
            }
            if (!ok) break;
            if ((u64)(m + 3) % (u64)(p - 1) != 0) {
                ok = false;
                break;
            }
        }
        if (ok) sum += (u128)m;
    }
    return (u64)sum;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    // Validation checkpoints from the problem statement.
    {
        u64 s100 = brute_sum(100);
        if (s100 != 32) {
            cerr << "Checkpoint failed: S(100)=" << s100 << " (expected 32)\n";
            return 1;
        }
        u64 s1e6 = brute_sum(1000000ULL);
        if (s1e6 != 22868117ULL) {
            cerr << "Checkpoint failed: S(10^6)=" << s1e6 << " (expected 22868117)\n";
            return 1;
        }
    }

    const u64 N = 1000000000000ULL;
    Solver solver(N);
    u64 ans = solver.solve();
    cout << ans << "\n";
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.util.ArrayList;
import java.util.HashSet;
import java.util.List;
import java.util.Set;

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

    static long modInv(long a, long mod) {
        if (mod == 1)
            return 0;
        long m0 = mod;
        long y = 0, x = 1;
        a = a % mod;
        while (a > 1) {
            long q = a / mod;
            long t = mod;
            mod = a % mod;
            a = t;
            t = y;
            y = x - q * y;
            x = t;
        }
        if (x < 0)
            x += m0;
        return x;
    }

    static class CRTRes {
        long r, m;
        boolean ok;

        CRTRes(long r, long m, boolean ok) {
            this.r = r;
            this.m = m;
            this.ok = ok;
        }
    }

    static CRTRes crtMerge(long r1, long m1, long r2, long m2, long maxN) {
        if (m1 == 1)
            return new CRTRes(r2 % m2, m2, true);
        if (m2 == 1)
            return new CRTRes(r1 % m1, m1, true);
        long g = gcd(m1, m2);
        if ((r1 % g) != (r2 % g))
            return new CRTRes(0, 0, false);

        long m1_g = m1 / g;
        long m2_g = m2 / g;

        long diff = r2 - r1;
        long diff_g = diff / g;

        long inv = modInv(m1_g % m2_g, m2_g);
        if (inv == 0 && m2_g > 1)
            return new CRTRes(0, 0, false);

        long t = diff_g % m2_g;
        if (t < 0)
            t += m2_g;

        // Java handles 128-bit multiplication via BigInteger but here values fit in
        // long if we are careful.
        // Wait, m1_g * m2 can exceed long! But maxN = 10^12. If m exceeds 10^12, we
        // clamp.
        // Let's use BigInteger for CRT steps just to be safe.
        java.math.BigInteger big_t = java.math.BigInteger.valueOf(t).multiply(java.math.BigInteger.valueOf(inv))
                .mod(java.math.BigInteger.valueOf(m2_g));
        java.math.BigInteger newm = java.math.BigInteger.valueOf(m1_g).multiply(java.math.BigInteger.valueOf(m2));
        java.math.BigInteger newr = java.math.BigInteger.valueOf(r1)
                .add(java.math.BigInteger.valueOf(m1).multiply(big_t));

        long clamp_mod = maxN + 1;
        if (newm.compareTo(java.math.BigInteger.valueOf(clamp_mod)) > 0) {
            java.math.BigInteger rr = newr.remainder(newm);
            if (rr.compareTo(java.math.BigInteger.valueOf(maxN)) > 0)
                return new CRTRes(0, 0, false);
            return new CRTRes(rr.longValue(), clamp_mod, true);
        }
        return new CRTRes(newr.remainder(newm).longValue(), newm.longValue(), true);
    }

    static List<Integer> sievePrimes(int n) {
        boolean[] isPrime = new boolean[n + 1];
        java.util.Arrays.fill(isPrime, true);
        isPrime[0] = isPrime[1] = false;
        for (int i = 2; i * i <= n; i++) {
            if (isPrime[i]) {
                for (int j = i * i; j <= n; j += i)
                    isPrime[j] = false;
            }
        }
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; i++)
            if (isPrime[i])
                primes.add(i);
        return primes;
    }

    static class Solver {
        long N;
        List<Integer> primes;
        Set<Long> found;
        static final long MAX_ENUM = 10000;

        Solver(long n) {
            N = n;
            int limit = (int) Math.sqrt(N) + 4;
            primes = sievePrimes(limit);
            found = new HashSet<>();
        }

        long[] congruenceForRemaining(List<Integer> B, long P) {
            long r = 0;
            long M = 1;
            for (int p : B) {
                long mod = p - 1;
                if (mod == 0)
                    return null;
                long a = (P / p) % mod;
                long g = gcd(a, mod);
                if (3 % g != 0)
                    return null;
                long mod2 = mod / g;
                if (mod2 == 1)
                    continue;

                long a2 = a / g;
                long rhs = (mod2 - ((3 / g) % mod2)) % mod2;
                long inv = modInv(a2 % mod2, mod2);
                if (inv == 0)
                    return null;
                long rp = (rhs * inv) % mod2;

                CRTRes merged = crtMerge(r, M, rp, mod2, N);
                if (!merged.ok)
                    return null;
                r = merged.r;
                M = merged.m;
            }
            return new long[] { r, M };
        }

        long maxRemainingProduct(int bound_exclusive, long cap) {
            long prod = 1;
            for (int p : primes) {
                if (p == 2)
                    continue;
                if (p >= bound_exclusive)
                    break;
                if (prod > cap / p)
                    return cap;
                prod *= p;
            }
            if (prod > cap)
                return cap;
            return prod;
        }

        boolean factorSquarefreeBounded(long x, int bound_exclusive, List<Integer> factorsOut) {
            factorsOut.clear();
            if (x == 0 || (x & 1) == 0)
                return false;
            long n = x;
            for (int p : primes) {
                if ((long) p * p > n)
                    break;
                if (p >= bound_exclusive)
                    break;
                if (n % p == 0) {
                    factorsOut.add(p);
                    n /= p;
                    if (n % p == 0)
                        return false;
                }
            }
            if (n > 1) {
                if (n >= bound_exclusive)
                    return false;
                factorsOut.add((int) n);
            }
            return true;
        }

        boolean checkSolution(long m, List<Integer> factors) {
            long mp3 = m + 3;
            for (int p : factors) {
                long d = p - 1;
                if (d == 0 || mp3 % d != 0)
                    return false;
            }
            return true;
        }

        void enumerateAndCheck(List<Integer> B, long P, int min_big, long r, long M) {
            long limit = N / P;
            limit = Math.min(limit, maxRemainingProduct(min_big, limit));
            if (M == 0)
                return;

            List<Integer> factorsR = new ArrayList<>();
            List<Integer> factorsM = new ArrayList<>();

            for (long R = r; R <= limit;) {
                if (R != 0) {
                    if (factorSquarefreeBounded(R, min_big, factorsR)) {
                        factorsM.clear();
                        factorsM.addAll(B);
                        factorsM.addAll(factorsR);
                        long m = P * R;
                        if (m <= N && checkSolution(m, factorsM)) {
                            found.add(m);
                        }
                    }
                }
                if (R > limit - M)
                    break;
                R += M;
            }
        }

        void dfs(int max_prime_idx, List<Integer> B, long P) {
            int min_big = B.get(B.size() - 1);
            long[] rm = congruenceForRemaining(B, P);
            if (rm == null)
                return;
            long r = rm[0], M = rm[1];
            if (M == 0)
                return;

            long limit = N / P;
            limit = Math.min(limit, maxRemainingProduct(min_big, limit));
            if (limit == 0)
                return;

            long approx_cnt = limit / M + 1;
            if (approx_cnt <= MAX_ENUM) {
                enumerateAndCheck(B, P, min_big, r, M);
                return;
            }

            for (int idx = max_prime_idx - 1; idx >= 0; idx--) {
                int q = primes.get(idx);
                if (q == 2)
                    continue;
                if (P > N / q)
                    continue;
                long P2 = P * q;
                B.add(q);
                dfs(idx, B, P2);
                B.remove(B.size() - 1);
            }
        }

        long solve() {
            found.add(1L);
            if (N >= 2)
                found.add(2L);
            if (N >= 3)
                found.add(3L);
            if (N >= 5)
                found.add(5L);

            for (int i = primes.size() - 1; i >= 0; i--) {
                int pmax = primes.get(i);
                if (pmax < 7)
                    break;
                List<Integer> B = new ArrayList<>();
                B.add(pmax);
                dfs(i, B, pmax);
            }

            long sum = 0;
            for (long m : found) {
                if (m <= N)
                    sum += m;
            }
            return sum;
        }
    }

    public static String solve() {
        Solver solver = new Solver(1000000000000L);
        return Long.toString(solver.solve());
    }

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