Problem 479: Roots on the Rise

Project Euler Problem 479 Solution

EulerSolve provides an optimized solution for Project Euler Problem 479, Roots on the Rise, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The goal is to evaluate $$S(n)=\sum_{k=1}^{n}\sum_{p=1}^{n}(1-k^2)^p \pmod{M},\qquad M=10^9+7.$$ For the actual input $n=10^6$, a direct double loop would require about $10^{12}$ term updates, so the computation must be reorganized. The key observation is that for each fixed $k$, the inner sum is a finite geometric series. Once that series is collapsed to a closed form, the problem becomes a single loop over $k$, with every contribution evaluated by modular exponentiation. A useful checkpoint is $$S(4)=51160,$$ which the implementation uses to confirm that the closed form matches direct evaluation on a small case. Mathematical Approach The entire method comes from simplifying the inner sum for one fixed value of $k$. Step 1: Isolate the Inner Sum for One $k$ Define $$u_k=1-k^2.$$ Then the contribution of this $k$ to the outer sum is $$T_k=\sum_{p=1}^{n}u_k^p.$$ Therefore $$S(n)=\sum_{k=1}^{n}T_k.$$ This separates the problem into $n$ independent geometric sums. Step 2: Collapse the Geometric Series For $u_k\neq 1$, the usual finite geometric-series identity gives $$T_k=\frac{u_k^{n+1}-u_k}{u_k-1}=u_k\cdot\frac{u_k^n-1}{u_k-1}.$$ This is exactly the form used by the implementation: one fast power computes $u_k^n$, and the remaining work is a small number of modular multiplications....

Detailed mathematical approach

Problem Summary

The goal is to evaluate

$$S(n)=\sum_{k=1}^{n}\sum_{p=1}^{n}(1-k^2)^p \pmod{M},\qquad M=10^9+7.$$

For the actual input $n=10^6$, a direct double loop would require about $10^{12}$ term updates, so the computation must be reorganized. The key observation is that for each fixed $k$, the inner sum is a finite geometric series. Once that series is collapsed to a closed form, the problem becomes a single loop over $k$, with every contribution evaluated by modular exponentiation.

A useful checkpoint is

$$S(4)=51160,$$

which the implementation uses to confirm that the closed form matches direct evaluation on a small case.

Mathematical Approach

The entire method comes from simplifying the inner sum for one fixed value of $k$.

Step 1: Isolate the Inner Sum for One $k$

Define

$$u_k=1-k^2.$$

Then the contribution of this $k$ to the outer sum is

$$T_k=\sum_{p=1}^{n}u_k^p.$$

Therefore

$$S(n)=\sum_{k=1}^{n}T_k.$$

This separates the problem into $n$ independent geometric sums.

Step 2: Collapse the Geometric Series

For $u_k\neq 1$, the usual finite geometric-series identity gives

$$T_k=\frac{u_k^{n+1}-u_k}{u_k-1}=u_k\cdot\frac{u_k^n-1}{u_k-1}.$$

This is exactly the form used by the implementation: one fast power computes $u_k^n$, and the remaining work is a small number of modular multiplications.

Step 3: Simplify the Denominator

Because $u_k=1-k^2$, we have

$$u_k-1=-k^2.$$

So the same quantity can also be written as

$$T_k=\frac{u_k-u_k^{n+1}}{k^2}.$$

This makes the arithmetic easier to interpret: the denominator is just $k^2$. Under the Euler modulus $M=10^9+7$, every $k$ in the range $1\le k\le 10^6$ is nonzero modulo $M$, so the denominator is invertible for the actual problem input.

Step 4: Replace Division by a Modular Inverse

All calculations are performed modulo the prime $M$. For any nonzero residue $a$, Fermat's little theorem gives

$$a^{-1}\equiv a^{M-2}\pmod{M}.$$

Hence, for $u_k\not\equiv 1\pmod{M}$,

$$T_k\equiv u_k\,(u_k^n-1)\,(u_k-1)^{-1}\pmod{M}.$$

This turns the whole computation into repeated modular exponentiation and multiplication, which is far cheaper than summing $n$ powers for each $k$.

Step 5: Handle the Degenerate Case Safely

If $u_k\equiv 1\pmod{M}$, then every term in the inner sum equals $1$, so

$$T_k\equiv n\pmod{M}.$$

For the actual Project Euler parameters this case does not occur, because $u_k\equiv 1$ would imply $k^2\equiv 0\pmod{M}$, impossible for $1\le k<M$. Still, the implementation keeps this branch so the formula never attempts to divide by zero modulo $M$.

Worked Example: $n=4$

Now compute the checkpoint directly.

For $k=1$, $u_1=1-1^2=0$, so

$$T_1=0+0+0+0=0.$$

For $k=2$, $u_2=-3$, hence

$$T_2=(-3)+(-3)^2+(-3)^3+(-3)^4=-3+9-27+81=60.$$

For $k=3$, $u_3=-8$, so

$$T_3=-8+64-512+4096=3640.$$

For $k=4$, $u_4=-15$, therefore

$$T_4=-15+225-3375+50625=47460.$$

Adding the four contributions gives

$$S(4)=0+60+3640+47460=51160,$$

which matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations all follow the same mathematical structure. They iterate once over $k=1,2,\dots,n$, reduce $1-k^2$ into the range $0$ to $M-1$, and then evaluate the closed form for $T_k$ modulo $M$.

For each $k$, the implementation computes the power $u_k^n \bmod M$ with binary exponentiation. If the denominator is nonzero, it computes the modular inverse by raising that denominator to the power $M-2$, again with binary exponentiation. The resulting contribution is added to the running total modulo $M$.

One implementation also includes two small validation steps: the known checkpoint $S(4)=51160$, and a comparison between the optimized formula and direct summation for a small input. Those checks confirm that the algebraic reduction is correct before the large target value is evaluated.

Complexity Analysis

There is one outer loop over $k=1$ to $n$. For each $k$, the implementation performs a constant amount of arithmetic plus modular exponentiation for $u_k^n$, and modular exponentiation for the inverse. Binary exponentiation costs $O(\log n)$ for the first power and $O(\log M)$ for the inverse, so the full running time is

$$O\bigl(n(\log n+\log M)\bigr).$$

With the fixed Euler modulus $M=10^9+7$, this is effectively $O(n\log n)$. The memory usage is

$$O(1),$$

since the program stores only a handful of current residues. This is a major improvement over the naive $O(n^2)$ direct summation.

Footnotes and References

Problem page: https://projecteuler.net/problem=479
Geometric series: Wikipedia — Geometric series
Fermat's little theorem: Wikipedia — Fermat's little theorem
Modular exponentiation: Wikipedia — Modular exponentiation

Problem 479 source code

C++

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

namespace {

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

struct Options {
    int n = 1'000'000;
    u64 mod = 1'000'000'007ULL;
    bool run_checkpoints = true;
};

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

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

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--n=", options.n) ||
            parse_u64_after_prefix(arg, "--mod=", options.mod)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1 && options.mod >= 2ULL;
}

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

u64 mod_inverse(const u64 a, const u64 mod) {
    return mod_pow(a % mod, mod - 2ULL, mod);
}

u64 solve(const int n, const u64 mod) {
    u64 ans = 0ULL;
    for (int k = 1; k <= n; ++k) {
        const i64 kk = static_cast<i64>(k);
        i64 t = 1LL - kk * kk;
        t %= static_cast<i64>(mod);
        if (t < 0) {
            t += static_cast<i64>(mod);
        }
        const u64 tm = static_cast<u64>(t);

        u64 contribution = 0ULL;
        if (tm == 1ULL) {
            contribution = static_cast<u64>(n) % mod;
        } else {
            const u64 pn = mod_pow(tm, static_cast<u64>(n), mod);
            const u64 numerator = static_cast<u64>((static_cast<u128>(tm) * ((pn + mod - 1ULL) % mod)) % mod);
            const u64 denom = (tm + mod - 1ULL) % mod;
            contribution = static_cast<u64>((static_cast<u128>(numerator) * mod_inverse(denom, mod)) % mod);
        }
        ans += contribution;
        if (ans >= mod) {
            ans -= mod;
        }
    }
    return ans;
}

u64 brute(const int n, const u64 mod) {
    u64 s = 0ULL;
    for (int k = 1; k <= n; ++k) {
        const i64 kk = static_cast<i64>(k);
        i64 u = 1LL - kk * kk;
        u %= static_cast<i64>(mod);
        if (u < 0) {
            u += static_cast<i64>(mod);
        }
        u64 term = 1ULL;
        for (int p = 1; p <= n; ++p) {
            term = static_cast<u64>((static_cast<u128>(term) * static_cast<u64>(u)) % mod);
            s += term;
            if (s >= mod) {
                s -= mod;
            }
        }
    }
    return s;
}

bool run_checkpoints() {
    const u64 mod = 1'000'000'007ULL;
    if (solve(4, mod) != 51'160ULL) {
        std::cerr << "Checkpoint failed: S(4)=51160" << '\n';
        return false;
    }
    if (solve(20, mod) != brute(20, mod)) {
        std::cerr << "Checkpoint failed: brute-force cross-check for n=20" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << solve(options.n, options.mod) << '\n';
    return 0;
}

Python

def solve():
    N = 1_000_000
    MOD = 1_000_000_007

    def mod_pow(base, exp, mod):
        result = 1
        base %= mod
        while exp > 0:
            if exp & 1: result = result * base % mod
            base = base * base % mod
            exp >>= 1
        return result

    def mod_inv(a, mod): return mod_pow(a % mod, mod - 2, mod)

    ans = 0
    for k in range(1, N + 1):
        t = (1 - k * k) % MOD
        if t == 1:
            ans = (ans + N % MOD) % MOD
        else:
            pn = mod_pow(t, N, MOD)
            num = t * ((pn - 1) % MOD) % MOD
            den = (t - 1) % MOD
            ans = (ans + num * mod_inv(den, MOD)) % MOD

    return str(ans % MOD)

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

Java

public class Euler479 {

    private static long modPow(long base, long exp, long mod) {
        long result = 1;
        base %= mod;
        while (exp > 0) {
            if ((exp & 1) == 1) {
                result = (result * base) % mod;
            }
            base = (base * base) % mod;
            exp >>= 1;
        }
        return result;
    }

    private static long modInverse(long a, long mod) {
        return modPow(a % mod, mod - 2, mod);
    }

    private static long solve(int n, long mod) {
        long ans = 0;
        for (long k = 1; k <= n; ++k) {
            long t = 1 - k * k;
            t %= mod;
            if (t < 0) {
                t += mod;
            }

            long contribution = 0;
            if (t == 1) {
                contribution = n % mod;
            } else {
                long pn = modPow(t, n, mod);
                long numerator = (t * ((pn + mod - 1) % mod)) % mod;
                long denom = (t + mod - 1) % mod;
                contribution = (numerator * modInverse(denom, mod)) % mod;
            }

            ans += contribution;
            if (ans >= mod) {
                ans -= mod;
            }
        }
        return ans;
    }

    public static void main(String[] args) {
        int n = 1000000;
        long mod = 1000000007L;
        System.out.println(solve(n, mod));
    }
}