Problem 542: Geometric Progression with Maximum Sum

View on Project Euler

Project Euler Problem 542 Solution

EulerSolve provides an optimized solution for Project Euler Problem 542, Geometric Progression with Maximum Sum, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer \(k \ge 4\), let \(M(k)\) be the maximum possible sum of a geometric progression of positive integers with at least three terms, common ratio greater than \(1\), and largest term at most \(k\). The final target is the alternating total $$A(N)=\sum_{k=4}^{N}(-1)^k M(k).$$ The real input is extremely large, so the solution cannot scan all geometric progressions and certainly cannot evaluate \(M(k)\) one value of \(k\) at a time. The key is to compress both the family of relevant progressions and the places where the maximum can actually change. Mathematical Approach The three implementations are based on the same reduction: first write every admissible integer geometric progression in a canonical form, then keep only the candidates that can ever be optimal, and finally evaluate the alternating sum block by block. Step 1: Parametrize Every Integer Geometric Progression Take a geometric progression with \(p+1\) terms, where \(p \ge 2\), common ratio \(a/b\) in lowest terms, and \(a > b \ge 1\). Because all terms must be integers, the progression can be written as $$q b^p,\ q a b^{p-1},\ q a^2 b^{p-2},\ \dots,\ q a^p,$$ with some integer scale factor \(q \ge 1\)....

Detailed mathematical approach

Problem Summary

For each integer \(k \ge 4\), let \(M(k)\) be the maximum possible sum of a geometric progression of positive integers with at least three terms, common ratio greater than \(1\), and largest term at most \(k\). The final target is the alternating total

$$A(N)=\sum_{k=4}^{N}(-1)^k M(k).$$

The real input is extremely large, so the solution cannot scan all geometric progressions and certainly cannot evaluate \(M(k)\) one value of \(k\) at a time. The key is to compress both the family of relevant progressions and the places where the maximum can actually change.

Mathematical Approach

The three implementations are based on the same reduction: first write every admissible integer geometric progression in a canonical form, then keep only the candidates that can ever be optimal, and finally evaluate the alternating sum block by block.

Step 1: Parametrize Every Integer Geometric Progression

Take a geometric progression with \(p+1\) terms, where \(p \ge 2\), common ratio \(a/b\) in lowest terms, and \(a > b \ge 1\). Because all terms must be integers, the progression can be written as

$$q b^p,\ q a b^{p-1},\ q a^2 b^{p-2},\ \dots,\ q a^p,$$

with some integer scale factor \(q \ge 1\). Its largest term is \(q a^p\), and its sum is

$$q\sum_{i=0}^{p} a^i b^{p-i}=q\,\frac{a^{p+1}-b^{p+1}}{a-b}.$$

So for fixed \(a\), \(b\), and \(p\), the best scale factor allowed by the bound \(k\) is

$$q=\left\lfloor \frac{k}{a^p} \right\rfloor.$$

Step 2: The Best Ratio Always Has Consecutive Numerator and Denominator

For fixed \(a\), \(p\), and \(q\), the sum

$$q\sum_{i=0}^{p} a^i b^{p-i}$$

is strictly increasing in \(b\), because every term with exponent \(p-i > 0\) grows when \(b\) grows. Since \(b\) must satisfy \(1 \le b < a\), the best admissible choice is always

$$b=a-1.$$

Therefore only progressions of the form

$$q(a-1)^p,\ q a (a-1)^{p-1},\ \dots,\ q a^p$$

need to be considered. Their largest term and one-step contribution are

$$P(a,p)=a^p,\qquad \Delta(a,p)=a^{p+1}-(a-1)^{p+1}.$$

For a fixed \(k\), this candidate contributes

$$V_{a,p}(k)=\left\lfloor \frac{k}{P(a,p)} \right\rfloor \Delta(a,p).$$

Hence

$$M(k)=\max_{a \ge 2,\ p \ge 2} V_{a,p}(k).$$

Step 3: Discard Dominated Candidates and Keep Only Records

If two candidates satisfy

$$P_1 \le P_2,\qquad \Delta_1 \ge \Delta_2,$$

then for every \(k\),

$$\left\lfloor \frac{k}{P_1} \right\rfloor \Delta_1 \ge \left\lfloor \frac{k}{P_2} \right\rfloor \Delta_2,$$

so the second candidate can never be optimal. The implementations therefore build only the record frontier: whenever the jump amount \(\Delta\) increases, they keep the smallest possible period \(P\) that achieves this new record. Because \(a^p \le N\), the exponent \(p\) only ranges up to \(\lfloor \log_2 N \rfloor\), and for each exponent the next record can be found by binary search on \(a\).

Step 4: Compute an Activity Bound for Each Record

Even among record candidates, many are useful only for a finite interval of \(k\). Suppose candidate \(j\) has a better asymptotic slope than candidate \(i\), meaning

$$\frac{\Delta_j}{P_j} > \frac{\Delta_i}{P_i} \iff \Delta_j P_i > \Delta_i P_j.$$

Then

$$V_j(k)\ge \left(\frac{k}{P_j}-1\right)\Delta_j,\qquad V_i(k)\le \frac{k}{P_i}\Delta_i.$$

So candidate \(j\) is guaranteed to beat candidate \(i\) whenever

$$\left(\frac{k}{P_j}-1\right)\Delta_j > \frac{k}{P_i}\Delta_i,$$

which simplifies to

$$k > \frac{\Delta_j P_j P_i}{\Delta_j P_i-\Delta_i P_j}.$$

For each record candidate, the implementations take the minimum such bound over all stronger competitors. This produces a proven upper activity limit \(R\): once \(k > R\), that candidate can never again attain the maximum.

Step 5: \(M(k)\) Can Change Only at Event Points

For one fixed candidate, the quantity \(\left\lfloor k / P \right\rfloor\) changes only when \(k\) reaches a multiple of \(P\). Therefore \(V_{a,p}(k)\) is constant between consecutive multiples of \(P\). After activity pruning, \(M(k)\) can change only at the union of all multiples of retained periods \(P\) up to \(\min(N,R)\). These integers are the event points. Once they are collected and sorted, every interval between consecutive events has a constant value of \(M(k)\).

Step 6: Evaluate the Alternating Sum by Constant Blocks

If \(M(k)=C\) for every \(k \in [L,R]\), then that whole block contributes

$$C\sum_{k=L}^{R}(-1)^k.$$

The inner sign sum is elementary:

$$\sum_{k=L}^{R}(-1)^k= \begin{cases} 0, & R-L+1 \text{ is even},\\ 1, & R-L+1 \text{ is odd and } L \text{ is even},\\ -1, & R-L+1 \text{ is odd and } L \text{ is odd}. \end{cases}$$

So the total is accumulated blockwise. Exact recomputation of the maximum is needed only at event points, not at every integer up to \(N\).

Worked Example: Small Values up to \(12\)

The first useful candidates are

$$a=2,\ p=2:\quad (1,2,4),\qquad V_{2,2}(k)=7\left\lfloor \frac{k}{4} \right\rfloor,$$

$$a=2,\ p=3:\quad (1,2,4,8),\qquad V_{2,3}(k)=15\left\lfloor \frac{k}{8} \right\rfloor,$$

$$a=3,\ p=2:\quad (4,6,9),\qquad V_{3,2}(k)=19\left\lfloor \frac{k}{9} \right\rfloor.$$

Up to \(12\), the event points are \(4,8,9,12\). Therefore

$$\begin{aligned} 4 \le k \le 7&:&&M(k)=7,\\ 8 \le k \le 8&:&&M(k)=14,\\ 9 \le k \le 11&:&&M(k)=19,\\ 12 \le k \le 12&:&&M(k)=21. \end{aligned}$$

In particular, \(M(4)=7\), \(M(10)=19\), and \(M(12)=21\). The alternating total up to \(12\) is

$$A(12)=0+14-19+21=16,$$

because the block \([4,7]\) has even length and contributes \(0\). This is exactly the kind of block compression the implementations exploit at the full scale.

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. First they determine all admissible exponents \(p\) with \(2^p \le N\), and for each such exponent they use integer roots to bound the search range for \(a\). Then they repeatedly raise the current record jump amount and use binary search to find the smallest base \(a\) that beats it for each exponent. Among those contenders, the implementation keeps the one with the smallest period \(a^p\), which yields the next record candidate on the frontier.

Next, the implementations compare every pair of retained candidates and assign each one an activity limit using the inequality from Step 4. After that, they generate all event points by listing the multiples of every retained period up to the smaller of the global limit and that candidate's activity bound. The list is sorted and duplicates are removed.

Finally, they sweep from \(4\) to \(N\). Between two successive event points the maximum sum is constant, so the code adds either \(0\), the block value, or its negation according to the parity rule above. At an event point it recomputes the exact maximum over the still-relevant candidates. The C++ version uses wide integer arithmetic and multiprecision only where overflow could occur; Python relies on arbitrary-precision integers naturally; Java uses big integers where the crossover calculation needs them.

Complexity Analysis

Let \(C\) be the number of retained record candidates and \(E\) the number of distinct event points. Let \(P_{\max}=\lfloor \log_2 N \rfloor\). Building the record frontier requires scanning those exponents and performing binary searches on the admissible base range, so the preprocessing cost is driven by about \(C\) record extractions across \(P_{\max}\) exponent families. The pairwise activity-bound computation costs \(O(C^2)\). Event generation costs \(O(E)\) insertions plus \(O(E \log E)\) for sorting and deduplication. The final sweep is \(O(E \cdot C)\) in the worst case because an event-point recomputation may inspect every retained candidate, although the activity limits keep the practical active set much smaller. The memory usage is \(O(C+E)\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=542
  2. Geometric progression: Wikipedia — Geometric progression
  3. Geometric series: Wikipedia — Geometric series
  4. Floor function: Wikipedia — Floor and ceiling functions
  5. Upper envelope: Wikipedia — Upper envelope

Problem 542 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <vector>

#include <boost/multiprecision/cpp_int.hpp>

using boost::multiprecision::cpp_int;

using u64 = std::uint64_t;
using u128 = unsigned __int128;
using i128 = __int128_t;

namespace {

struct Candidate {
    u64 n = 0;          // Jump period: n = u^p.
    u64 d = 0;          // Jump amount: d = u^(p+1) - (u-1)^(p+1).
    u64 relevance = 0;  // Proven upper bound where candidate can still be optimal.
};

u64 pow_with_cap(const u64 base, const int exp, const u64 cap) {
    u128 value = 1;
    for (int i = 0; i < exp; ++i) {
        value *= static_cast<u128>(base);
        if (value > static_cast<u128>(cap)) {
            return cap + 1;
        }
    }
    return static_cast<u64>(value);
}

u64 int_root(const u64 n, const int p) {
    u64 lo = 1;
    u64 hi = 2;
    while (pow_with_cap(hi, p, n) <= n) {
        if (hi > n / 2) {
            hi = n;
            break;
        }
        hi *= 2;
    }
    while (lo + 1 < hi) {
        const u64 mid = lo + (hi - lo) / 2;
        if (pow_with_cap(mid, p, n) <= n) {
            lo = mid;
        } else {
            hi = mid;
        }
    }
    return lo;
}

u64 delta_value(const u64 u, const int p) {
    // Computes u^(p+1) - (u-1)^(p+1) with a stable recurrence:
    // D_m = u * D_{m-1} + (u-1)^(m-1), D_1 = 1.
    const u64 v = u - 1;
    u128 D = 1;
    u128 vpow = 1;
    for (int m = 2; m <= p + 1; ++m) {
        vpow *= static_cast<u128>(v);
        D = static_cast<u128>(u) * D + vpow;
    }
    return static_cast<u64>(D);
}

int max_power_for_limit(const u64 n) {
    int p = 0;
    u64 x = 1;
    while (x <= n / 2) {
        x *= 2;
        ++p;
    }
    return p;
}

std::vector<Candidate> generate_record_candidates(const u64 limit) {
    if (limit < 4) {
        return {};
    }

    const int p_max = max_power_for_limit(limit);
    std::vector<u64> u_max(static_cast<std::size_t>(p_max + 1), 0);
    for (int p = 2; p <= p_max; ++p) {
        u_max[static_cast<std::size_t>(p)] = int_root(limit, p);
    }

    std::vector<Candidate> records;
    u64 current_best_d = 0;

    while (true) {
        bool found = false;
        u64 best_n = 0;
        u64 best_d = 0;

        for (int p = 2; p <= p_max; ++p) {
            const u64 umax = u_max[static_cast<std::size_t>(p)];
            if (umax < 2) {
                continue;
            }
            if (delta_value(umax, p) <= current_best_d) {
                continue;
            }

            u64 lo = 2;
            u64 hi = umax;
            while (lo < hi) {
                const u64 mid = lo + (hi - lo) / 2;
                if (delta_value(mid, p) > current_best_d) {
                    hi = mid;
                } else {
                    lo = mid + 1;
                }
            }

            const u64 u = lo;
            const u64 n = pow_with_cap(u, p, limit);
            const u64 d = delta_value(u, p);
            if (!found || n < best_n || (n == best_n && d > best_d)) {
                found = true;
                best_n = n;
                best_d = d;
            }
        }

        if (!found) {
            break;
        }

        records.push_back({best_n, best_d, limit});
        current_best_d = best_d;
    }

    return records;
}

void compute_relevance_bounds(std::vector<Candidate>& candidates, const u64 limit) {
    const std::size_t n = candidates.size();
    for (std::size_t i = 0; i < n; ++i) {
        u64 best = limit;
        const u64 ni = candidates[i].n;
        const u64 di = candidates[i].d;

        for (std::size_t j = 0; j < n; ++j) {
            const u64 nj = candidates[j].n;
            const u64 dj = candidates[j].d;

            const u128 lhs = static_cast<u128>(dj) * static_cast<u128>(ni);
            const u128 rhs = static_cast<u128>(di) * static_cast<u128>(nj);
            if (lhs <= rhs) {
                continue;
            }

            // For k > num/den, candidate j is guaranteed > candidate i:
            // num = dj*nj*ni, den = dj*ni - di*nj.
            const cpp_int num = cpp_int(dj) * cpp_int(nj) * cpp_int(ni);
            const cpp_int den = cpp_int(dj) * cpp_int(ni) - cpp_int(di) * cpp_int(nj);
            const cpp_int r = num / den;  // floor(num / den)

            if (r < cpp_int(best)) {
                best = static_cast<u64>(r);
            }
        }

        candidates[i].relevance = best;
    }
}

std::vector<u64> collect_events(const std::vector<Candidate>& candidates, const u64 limit) {
    std::vector<u64> events;
    events.reserve(16384);

    for (const Candidate& c : candidates) {
        const u64 lim = std::min(limit, c.relevance);
        if (lim < c.n) {
            continue;
        }

        u64 t = c.n;
        while (t <= lim) {
            events.push_back(t);
            if (t > lim - c.n) {
                break;
            }
            t += c.n;
        }
    }

    std::sort(events.begin(), events.end());
    events.erase(std::unique(events.begin(), events.end()), events.end());
    return events;
}

i128 to_i128(const u64 x) {
    return static_cast<i128>(x);
}

std::string i128_to_string(i128 value) {
    if (value == 0) {
        return "0";
    }
    bool neg = value < 0;
    if (neg) {
        value = -value;
    }
    std::string out;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        out.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    if (neg) {
        out.push_back('-');
    }
    std::reverse(out.begin(), out.end());
    return out;
}

int alternating_sign_sum(const u64 l, const u64 r) {
    if (l > r) {
        return 0;
    }
    const u64 len = r - l + 1;
    if ((len & 1ULL) == 0ULL) {
        return 0;
    }
    return (l & 1ULL) ? -1 : 1;
}

class Solver542 {
public:
    explicit Solver542(const u64 limit) : limit_(limit) {
        candidates_ = generate_record_candidates(limit_);
        compute_relevance_bounds(candidates_, limit_);
        events_ = collect_events(candidates_, limit_);
    }

    u64 S(const u64 k) const {
        if (k < 4) {
            return 0;
        }
        u128 best = 0;
        for (const Candidate& c : candidates_) {
            if (c.n > k) {
                break;
            }
            if (k > c.relevance) {
                continue;
            }
            const u128 value = static_cast<u128>(k / c.n) * static_cast<u128>(c.d);
            if (value > best) {
                best = value;
            }
        }
        return static_cast<u64>(best);
    }

    i128 T(const u64 n) const {
        if (n < 4) {
            return 0;
        }
        if (n > limit_) {
            throw std::runtime_error("Requested n exceeds solver limit.");
        }

        i128 answer = 0;
        u64 current = 4;
        u64 s_value = S(4);
        std::size_t event_idx = 0;

        while (current <= n) {
            const u64 next_event =
                (event_idx < events_.size() && events_[event_idx] <= n) ? events_[event_idx] : (n + 1);

            if (next_event > current) {
                const u64 r = std::min(n, next_event - 1);
                const int sign = alternating_sign_sum(current, r);
                if (sign != 0) {
                    answer += static_cast<i128>(sign) * to_i128(s_value);
                }
                current = r + 1;
                if (current > n) {
                    break;
                }
            }

            // At event points, at least one candidate jumps; recompute exact S(current).
            s_value = S(current);
            ++event_idx;
        }

        return answer;
    }

private:
    u64 limit_;
    std::vector<Candidate> candidates_;
    std::vector<u64> events_;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0) != 0) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    u64 parsed = 0;
    for (char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(ch - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }
    value = parsed;
    return true;
}

void run_checkpoints() {
    Solver542 check_solver(1000);
    if (check_solver.S(4) != 7) {
        throw std::runtime_error("Checkpoint failed: S(4) != 7");
    }
    if (check_solver.S(10) != 19) {
        throw std::runtime_error("Checkpoint failed: S(10) != 19");
    }
    if (check_solver.S(12) != 21) {
        throw std::runtime_error("Checkpoint failed: S(12) != 21");
    }
    if (check_solver.S(1000) != 3439) {
        throw std::runtime_error("Checkpoint failed: S(1000) != 3439");
    }
    if (check_solver.T(1000) != 2268) {
        throw std::runtime_error("Checkpoint failed: T(1000) != 2268");
    }
}

}  // namespace

int main(int argc, char** argv) {
    u64 n = 100000000000000000ULL;
    bool run_checks = true;

    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            run_checks = false;
            continue;
        }
        u64 parsed = 0;
        if (parse_u64_after_prefix(arg, "--n=", parsed)) {
            n = parsed;
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return 1;
    }

    try {
        if (run_checks) {
            run_checkpoints();
        }
        Solver542 solver(n);
        const i128 answer = solver.T(n);
        std::cout << i128_to_string(answer) << '\n';
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

def pow_with_cap(base, exp, cap):
    val = 1
    for _ in range(exp):
        val *= base
        if val > cap: return cap + 1
    return val

def int_root(n, p):
    lo = 1
    hi = 2
    while pow_with_cap(hi, p, n) <= n:
        if hi > n // 2:
            hi = n
            break
        hi *= 2
    
    while lo + 1 < hi:
        mid = (lo + hi) // 2
        if pow_with_cap(mid, p, n) <= n:
            lo = mid
        else:
            hi = mid
    return lo

def delta_value(u, p):
    v = u - 1
    D = 1
    vpow = 1
    for m in range(2, p + 2):
        vpow *= v
        D = u * D + vpow
    return D

def max_power_for_limit(n):
    p = 0
    x = 1
    while x <= n // 2:
        x *= 2
        p += 1
    return p

class Candidate:
    def __init__(self, n, d, limit):
        self.n = n
        self.d = d
        self.relevance = limit

def generate_record_candidates(limit):
    if limit < 4: return []
    p_max = max_power_for_limit(limit)
    u_max = [0] * (p_max + 1)
    for p in range(2, p_max + 1):
        u_max[p] = int_root(limit, p)
        
    records = []
    current_best_d = 0
    
    while True:
        found = False
        best_n = 0
        best_d = 0
        
        for p in range(2, p_max + 1):
            umax = u_max[p]
            if umax < 2: continue
            if delta_value(umax, p) <= current_best_d: continue
            
            lo, hi = 2, umax
            while lo < hi:
                mid = (lo + hi) // 2
                if delta_value(mid, p) > current_best_d:
                    hi = mid
                else:
                    lo = mid + 1
                    
            u = lo
            n = pow_with_cap(u, p, limit)
            d = delta_value(u, p)
            if not found or n < best_n or (n == best_n and d > best_d):
                found = True
                best_n = n
                best_d = d
                
        if not found: break
        records.append(Candidate(best_n, best_d, limit))
        current_best_d = best_d
        
    return records

def compute_relevance_bounds(candidates, limit):
    for i, ci in enumerate(candidates):
        best = limit
        for j, cj in enumerate(candidates):
            if cj.d * ci.n <= ci.d * cj.n: continue
            num = cj.d * cj.n * ci.n
            den = cj.d * ci.n - ci.d * cj.n
            r = num // den
            if r < best:
                best = r
        ci.relevance = best

def collect_events(candidates, limit):
    events = set()
    for c in candidates:
        lim = min(limit, c.relevance)
        if lim < c.n: continue
        t = c.n
        while t <= lim:
            events.add(t)
            if t > lim - c.n: break
            t += c.n
    return sorted(list(events))

def alternating_sign_sum(l, r):
    if l > r: return 0
    length = r - l + 1
    if length % 2 == 0: return 0
    return -1 if l % 2 else 1

class Solver542:
    def __init__(self, limit):
        self.limit = limit
        self.candidates = generate_record_candidates(limit)
        compute_relevance_bounds(self.candidates, limit)
        self.events = collect_events(self.candidates, limit)
        
    def S(self, k):
        if k < 4: return 0
        best = 0
        for c in self.candidates:
            if c.n > k: break
            if k > c.relevance: continue
            val = (k // c.n) * c.d
            if val > best: best = val
        return best
        
    def T(self, n):
        if n < 4: return 0
        ans = 0
        current = 4
        s_val = self.S(4)
        event_idx = 0
        events_len = len(self.events)
        
        while current <= n:
            next_event = self.events[event_idx] if event_idx < events_len and self.events[event_idx] <= n else n + 1
            if next_event > current:
                r = min(n, next_event - 1)
                sign = alternating_sign_sum(current, r)
                if sign != 0:
                    ans += sign * s_val
                current = r + 1
                if current > n: break
            s_val = self.S(current)
            event_idx += 1
        return ans

def solve():
    s = Solver542(100000000000000000)
    return str(s.T(100000000000000000))

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

Java

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

public class Euler542 {

    static class Candidate {
        long n;
        long d;
        long relevance;

        Candidate(long n, long d, long relevance) {
            this.n = n;
            this.d = d;
            this.relevance = relevance;
        }
    }

    static long powWithCap(long base, int exp, long cap) {
        BigInteger val = BigInteger.ONE;
        BigInteger b = BigInteger.valueOf(base);
        BigInteger c = BigInteger.valueOf(cap);
        for (int i = 0; i < exp; i++) {
            val = val.multiply(b);
            if (val.compareTo(c) > 0) {
                return cap + 1;
            }
        }
        return val.longValue();
    }

    static long intRoot(long n, int p) {
        long lo = 1;
        long hi = 2;
        while (powWithCap(hi, p, n) <= n) {
            if (hi > n / 2) {
                hi = n;
                break;
            }
            hi *= 2;
        }
        while (lo + 1 < hi) {
            long mid = lo + (hi - lo) / 2;
            if (powWithCap(mid, p, n) <= n) {
                lo = mid;
            } else {
                hi = mid;
            }
        }
        return lo;
    }

    static long deltaValue(long u, int p) {
        long v = u - 1;
        BigInteger D = BigInteger.ONE;
        BigInteger vpow = BigInteger.ONE;
        BigInteger bigU = BigInteger.valueOf(u);
        BigInteger bigV = BigInteger.valueOf(v);
        for (int m = 2; m <= p + 1; m++) {
            vpow = vpow.multiply(bigV);
            D = bigU.multiply(D).add(vpow);
        }
        return D.longValue();
    }

    static int maxPowerForLimit(long n) {
        int p = 0;
        long x = 1;
        while (x <= n / 2) {
            x *= 2;
            p++;
        }
        return p;
    }

    static List<Candidate> generateRecordCandidates(long limit) {
        if (limit < 4)
            return new ArrayList<>();

        int pMax = maxPowerForLimit(limit);
        long[] uMax = new long[pMax + 1];
        for (int p = 2; p <= pMax; p++) {
            uMax[p] = intRoot(limit, p);
        }

        List<Candidate> records = new ArrayList<>();
        long currentBestD = 0;

        while (true) {
            boolean found = false;
            long bestN = 0;
            long bestD = 0;

            for (int p = 2; p <= pMax; p++) {
                long umax = uMax[p];
                if (umax < 2)
                    continue;
                if (deltaValue(umax, p) <= currentBestD)
                    continue;

                long lo = 2;
                long hi = umax;
                while (lo < hi) {
                    long mid = lo + (hi - lo) / 2;
                    if (deltaValue(mid, p) > currentBestD) {
                        hi = mid;
                    } else {
                        lo = Math.min(umax, mid + 1); // Avoid endless loops with mid+1 properly bounded
                    }
                }

                long u = lo;
                long n = powWithCap(u, p, limit);
                long d = deltaValue(u, p);
                if (!found || n < bestN || (n == bestN && d > bestD)) {
                    found = true;
                    bestN = n;
                    bestD = d;
                }
            }

            if (!found)
                break;

            records.add(new Candidate(bestN, bestD, limit));
            currentBestD = bestD;
        }

        return records;
    }

    static void computeRelevanceBounds(List<Candidate> candidates, long limit) {
        int n = candidates.size();
        for (int i = 0; i < n; i++) {
            long best = limit;
            long ni = candidates.get(i).n;
            long di = candidates.get(i).d;

            for (int j = 0; j < n; j++) {
                long nj = candidates.get(j).n;
                long dj = candidates.get(j).d;

                BigInteger lhs = BigInteger.valueOf(dj).multiply(BigInteger.valueOf(ni));
                BigInteger rhs = BigInteger.valueOf(di).multiply(BigInteger.valueOf(nj));
                if (lhs.compareTo(rhs) <= 0)
                    continue;

                BigInteger num = BigInteger.valueOf(dj).multiply(BigInteger.valueOf(nj))
                        .multiply(BigInteger.valueOf(ni));
                BigInteger den = BigInteger.valueOf(dj).multiply(BigInteger.valueOf(ni))
                        .subtract(BigInteger.valueOf(di).multiply(BigInteger.valueOf(nj)));

                BigInteger r = num.divide(den);
                if (r.compareTo(BigInteger.valueOf(best)) < 0) {
                    best = r.longValue();
                }
            }
            candidates.get(i).relevance = best;
        }
    }

    static List<Long> collectEvents(List<Candidate> candidates, long limit) {
        List<Long> events = new ArrayList<>();

        for (Candidate c : candidates) {
            long lim = Math.min(limit, c.relevance);
            if (lim < c.n)
                continue;

            long t = c.n;
            while (t <= lim) {
                events.add(t);
                if (t > lim - c.n)
                    break;
                t += c.n;
            }
        }

        Collections.sort(events);
        List<Long> uniqueEvents = new ArrayList<>();
        if (!events.isEmpty()) {
            uniqueEvents.add(events.get(0));
            for (int i = 1; i < events.size(); i++) {
                if (!events.get(i).equals(events.get(i - 1))) {
                    uniqueEvents.add(events.get(i));
                }
            }
        }
        return uniqueEvents;
    }

    static int alternatingSignSum(long l, long r) {
        if (l > r)
            return 0;
        long len = r - l + 1;
        if ((len & 1) == 0)
            return 0;
        return (l & 1) != 0 ? -1 : 1;
    }

    static class Solver {
        long limit;
        List<Candidate> candidates;
        List<Long> events;

        Solver(long limit) {
            this.limit = limit;
            this.candidates = generateRecordCandidates(limit);
            computeRelevanceBounds(this.candidates, limit);
            this.events = collectEvents(this.candidates, limit);
        }

        long S(long k) {
            if (k < 4)
                return 0;
            long best = 0;
            for (Candidate c : candidates) {
                if (c.n > k)
                    break;
                if (k > c.relevance)
                    continue;
                long val = (k / c.n) * c.d;
                if (val > best) {
                    best = val;
                }
            }
            return best;
        }

        BigInteger T(long n) {
            if (n < 4)
                return BigInteger.ZERO;

            BigInteger ans = BigInteger.ZERO;
            long current = 4;
            long sValue = S(4);
            int eventIdx = 0;

            while (current <= n) {
                long nextEvent = (eventIdx < events.size() && events.get(eventIdx) <= n) ? events.get(eventIdx)
                        : (n + 1);

                if (nextEvent > current) {
                    long r = Math.min(n, nextEvent - 1);
                    int sign = alternatingSignSum(current, r);
                    if (sign != 0) {
                        if (sign == 1) {
                            ans = ans.add(BigInteger.valueOf(sValue));
                        } else {
                            ans = ans.subtract(BigInteger.valueOf(sValue));
                        }
                    }
                    current = r + 1;
                    if (current > n)
                        break;
                }

                sValue = S(current);
                eventIdx++;
            }

            return ans;
        }
    }

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

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