Problem 739: Summation of Summations

Project Euler Problem 739 Solution

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

Problem Summary The target quantity for Problem 739 is evaluated modulo $M=10^9+7$ for $n=10^8$. After simplifying the repeated-summation construction, the required value can be written as a single weighted sum of Lucas numbers: $$f(n)=\sum_{t=1}^{r} T_t\,L_{t+1}\pmod{M},\qquad r=n-1.$$ The real difficulty is the size of $r$. A direct quadratic summation is hopeless, and even an $O(n)$ scan becomes too slow if it performs one modular exponentiation per term. The successful approach is to stream the coefficients backward, update Lucas values in lockstep, and batch the modular inverses. Mathematical Approach The solution combines Lucas-number identities, a ballot-number coefficient sequence, and blockwise modular inversion. Step 1: Start from Lucas Numbers The Lucas sequence satisfies $L_0=2$, $L_1=1$, and $$L_{k+1}=L_k+L_{k-1}.$$ The implementations do not build this sequence from the beginning. Instead, they first obtain the Fibonacci pair $(F_k,F_{k+1})$ by fast doubling and then use $$L_k=F_{k-1}+F_{k+1}=2F_{k+1}-F_k.$$ That gives the starting values $L_r$ and $L_{r+1}$ in $O(\log n)$ time, which is small compared with the main scan....

Detailed mathematical approach

Problem Summary

The target quantity for Problem 739 is evaluated modulo $M=10^9+7$ for $n=10^8$. After simplifying the repeated-summation construction, the required value can be written as a single weighted sum of Lucas numbers:

$$f(n)=\sum_{t=1}^{r} T_t\,L_{t+1}\pmod{M},\qquad r=n-1.$$

The real difficulty is the size of $r$. A direct quadratic summation is hopeless, and even an $O(n)$ scan becomes too slow if it performs one modular exponentiation per term. The successful approach is to stream the coefficients backward, update Lucas values in lockstep, and batch the modular inverses.

Mathematical Approach

The solution combines Lucas-number identities, a ballot-number coefficient sequence, and blockwise modular inversion.

Step 1: Start from Lucas Numbers

The Lucas sequence satisfies $L_0=2$, $L_1=1$, and

$$L_{k+1}=L_k+L_{k-1}.$$

The implementations do not build this sequence from the beginning. Instead, they first obtain the Fibonacci pair $(F_k,F_{k+1})$ by fast doubling and then use

$$L_k=F_{k-1}+F_{k+1}=2F_{k+1}-F_k.$$

That gives the starting values $L_r$ and $L_{r+1}$ in $O(\log n)$ time, which is small compared with the main scan.

Step 2: Identify the Weight Sequence

The required sum uses coefficients $T_t$ generated backward from the endpoint:

$$T_r=1,$$

$$T_{t-1}=T_t\cdot \frac{(t-1)(2r-t)}{t(r-t+1)}\qquad (2\le t\le r).$$

This is not an arbitrary recurrence. It is exactly the ratio of ballot-number, or Catalan-triangle, coefficients. A closed form is

$$T_t=\frac{t}{r}\binom{2r-t-1}{r-1}=\binom{2r-t-1}{r-1}-\binom{2r-t-1}{r}.$$

So the repeated-summation problem collapses to a Lucas sum with very structured combinatorial weights.

Step 3: Turn the Sum into a Backward Stream

Because both the coefficients and the Lucas sequence admit backward updates, the whole computation can be done in one descending pass $t=r,r-1,\dots,1$. Rearranging the Lucas recurrence gives

$$L_{t-1}=L_{t+1}-L_t\pmod{M}.$$

During the scan, one term contributes

$$T_tL_{t+1}\pmod{M},$$

then the coefficient is multiplied by the rational step factor above, and the Lucas pair is shifted backward by one index. No table of all coefficients or all Lucas numbers is needed.

Step 4: Worked Example for $n=8$

For $n=8$, we have $r=7$. Starting from $T_7=1$, the recurrence yields

$$T_7=1,\quad T_6=6,\quad T_5=20,\quad T_4=48,\quad T_3=90,\quad T_2=132,\quad T_1=132.$$

The required Lucas values are

$$L_2=3,\quad L_3=4,\quad L_4=7,\quad L_5=11,\quad L_6=18,\quad L_7=29,\quad L_8=47.$$

Therefore

$$f(8)=132\cdot 3+132\cdot 4+90\cdot 7+48\cdot 11+20\cdot 18+6\cdot 29+1\cdot 47=2663,$$

which matches the sample check used by the implementations.

Step 5: Remove Per-Term Modular Inversions

The denominator in the coefficient update is

$$d_t=t(r-t+1).$$

Since $r=10^8-1\lt M$, both factors are nonzero modulo $M$, so every $d_t$ is invertible. Computing each inverse separately would be far too expensive. Instead, the denominators are grouped into blocks. For one block $d_1,\dots,d_m$, a single inverse of the total product is enough:

$$d_i^{-1}=\left(\prod_{j=1}^{m} d_j\right)^{-1}\prod_{j\ne i} d_j \pmod{M}.$$

This replaces $m$ costly modular exponentiations by one exponentiation plus $O(m)$ multiplications.

How the Code Works

The C++, Python, and Java implementations all follow the same plan. First they compute the Lucas starting pair with fast doubling under modulo $M$. Then they process the coefficient recurrence from high $t$ down to low $t$ in fixed-size blocks, materializing the denominators $t(r-t+1)$ for the current block and batch-inverting them.

Inside each block, the implementation adds the current contribution $T_tL_{t+1}$, updates the coefficient using the precomputed inverse for that step, and moves the Lucas pair backward using $L_{t-1}=L_{t+1}-L_t$. This streaming design keeps the memory footprint small while still reproducing the validation checkpoints $f(8)=2663$ and $f(20)=742296999$.

Complexity Analysis

Fast doubling costs $O(\log n)$ time. The main scan visits each $t$ once, so the total work is $O(n)$ modular multiplications. Batch inversion keeps the number of modular exponentiations proportional to the number of blocks rather than the number of terms. With block size $B$, the memory usage is $O(B)$, and the total runtime is linear in $n$ apart from the tiny logarithmic startup cost.

Footnotes and References

Problem page: https://projecteuler.net/problem=739
Lucas numbers: Wikipedia — Lucas number
Fast doubling for Fibonacci numbers: cp-algorithms — Fibonacci numbers
Ballot-number background: Wikipedia — Bertrand's ballot theorem
Modular inverses: Wikipedia — Modular multiplicative inverse

Problem 739 source code

C++

#include <cassert>
#include <cstdint>
#include <iostream>
#include <utility>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u32 = std::uint32_t;

constexpr u64 kMod = 1'000'000'007ULL;

inline u64 mul_mod(const u64 a, const u64 b) {
    return (a * b) % kMod;
}

u64 mod_pow(u64 base, u64 exp) {
    u64 result = 1ULL;
    while (exp > 0ULL) {
        if ((exp & 1ULL) != 0ULL) {
            result = mul_mod(result, base);
        }
        base = mul_mod(base, base);
        exp >>= 1ULL;
    }
    return result;
}

std::pair<u64, u64> fib_pair(const u64 n) {
    if (n == 0ULL) {
        return {0ULL, 1ULL};
    }
    const auto [a, b] = fib_pair(n >> 1ULL);

    const u64 two_b = (2ULL * b) % kMod;
    const u64 two_b_minus_a = (two_b + kMod - a) % kMod;
    const u64 c = mul_mod(a, two_b_minus_a);
    const u64 d = (mul_mod(a, a) + mul_mod(b, b)) % kMod;

    if ((n & 1ULL) == 0ULL) {
        return {c, d};
    }
    return {d, (c + d) % kMod};
}

u64 lucas(const u64 n) {
    if (n == 0ULL) {
        return 2ULL;
    }
    const auto [fn, fn1] = fib_pair(n);
    u64 value = (2ULL * fn1) % kMod;
    value = (value + kMod - fn) % kMod;
    return value;
}

void batch_invert(const std::vector<u32>& values,
                  std::vector<u32>& inverses,
                  std::vector<u32>& prefix) {
    const std::size_t m = values.size();
    inverses.resize(m);
    if (m == 0U) {
        return;
    }

    prefix.resize(m);
    prefix[0] = values[0];
    for (std::size_t i = 1; i < m; ++i) {
        prefix[i] = static_cast<u32>(mul_mod(prefix[i - 1], values[i]));
    }

    u64 suffix_inv = mod_pow(prefix[m - 1], kMod - 2ULL);
    for (std::size_t i = m; i-- > 0U;) {
        const u64 left = (i == 0U) ? 1ULL : prefix[i - 1U];
        inverses[i] = static_cast<u32>(mul_mod(suffix_inv, left));
        suffix_inv = mul_mod(suffix_inv, values[i]);
    }
}

u64 f(const u64 n) {
    const u64 r = n - 1ULL;

    u64 t = r;
    u64 coeff = 1ULL;  // T_r = 1

    u64 l_t = lucas(r);
    u64 l_t_plus_1 = lucas(r + 1ULL);

    u64 answer = 0ULL;

    constexpr u64 kBlock = 1'000'000ULL;
    std::vector<u32> denoms;
    std::vector<u32> inv_denoms;
    std::vector<u32> prefix;
    denoms.reserve(kBlock);
    inv_denoms.reserve(kBlock);
    prefix.reserve(kBlock);

    while (t >= 1ULL) {
        const u64 hi = t;
        const u64 lo = (hi > kBlock) ? (hi - kBlock + 1ULL) : 1ULL;
        const std::size_t m = static_cast<std::size_t>(hi - lo + 1ULL);

        denoms.resize(m);
        for (std::size_t i = 0; i < m; ++i) {
            const u64 cur_t = hi - static_cast<u64>(i);
            if (cur_t == 1ULL) {
                denoms[i] = 1ULL;
            } else {
                const u64 u = r - cur_t + 1ULL;
                denoms[i] = static_cast<u32>(mul_mod(cur_t % kMod, u % kMod));
            }
        }
        batch_invert(denoms, inv_denoms, prefix);

        for (std::size_t i = 0; i < m; ++i) {
            const u64 cur_t = hi - static_cast<u64>(i);

            answer += mul_mod(coeff, l_t_plus_1);
            if (answer >= kMod) {
                answer -= kMod;
            }

            if (cur_t == 1ULL) {
                break;
            }

            const u64 num_a = cur_t - 1ULL;
            const u64 num_b = 2ULL * r - cur_t;
            const u64 num = mul_mod(num_a % kMod, num_b % kMod);

            coeff = mul_mod(coeff, num);
            coeff = mul_mod(coeff, inv_denoms[i]);

            const u64 next_l_t_plus_1 = l_t;
            const u64 next_l_t = (l_t_plus_1 + kMod - l_t) % kMod;
            l_t_plus_1 = next_l_t_plus_1;
            l_t = next_l_t;
        }

        if (lo == 1ULL) {
            break;
        }
        t = lo - 1ULL;
    }

    return answer;
}

}  // namespace

int main() {
    assert(f(8ULL) == 2663ULL);
    assert(f(20ULL) == 742'296'999ULL);

    std::cout << f(100'000'000ULL) << '\n';
    return 0;
}

Python

def solve():
    MOD = 1000000007
    n = 100000000
    r = n - 1

    def mul_mod(a, b): return a * b % MOD
    def mod_pow(base, exp):
        result = 1
        while exp > 0:
            if exp & 1: result = result * base % MOD
            base = base * base % MOD; exp >>= 1
        return result

    def fib_pair(n):
        if n == 0: return (0, 1)
        a, b = fib_pair(n >> 1)
        c = a * ((2*b - a) % MOD) % MOD
        d = (a*a + b*b) % MOD
        return (d, (c+d) % MOD) if n & 1 else (c, d)

    def lucas(n):
        if n == 0: return 2
        fn, fn1 = fib_pair(n)
        return (2*fn1 - fn) % MOD

    # Batch inverse
    def batch_inv(vals):
        m = len(vals)
        if m == 0: return []
        prefix = [0]*m; prefix[0] = vals[0]
        for i in range(1, m):
            prefix[i] = prefix[i-1] * vals[i] % MOD
        sinv = mod_pow(prefix[-1], MOD-2)
        inv = [0]*m
        for i in range(m-1, -1, -1):
            left = prefix[i-1] if i > 0 else 1
            inv[i] = sinv * left % MOD
            sinv = sinv * vals[i] % MOD
        return inv

    l_t = lucas(r); l_tp1 = lucas(r+1)
    coeff = 1; answer = 0; t = r
    BLOCK = 1000000

    while t >= 1:
        hi = t; lo = max(1, hi - BLOCK + 1)
        m = hi - lo + 1
        denoms = [0]*m
        for i in range(m):
            ct = hi - i
            if ct == 1: denoms[i] = 1
            else:
                u = r - ct + 1
                denoms[i] = ct % MOD * (u % MOD) % MOD
        inv_d = batch_inv(denoms)

        for i in range(m):
            ct = hi - i
            answer = (answer + coeff * l_tp1) % MOD
            if ct == 1: break
            num_a = (ct - 1) % MOD
            num_b = (2 * r - ct) % MOD
            num = num_a * num_b % MOD
            coeff = coeff * num % MOD * inv_d[i] % MOD
            nl_tp1 = l_t; nl_t = (l_tp1 - l_t) % MOD
            l_tp1 = nl_tp1; l_t = nl_t

        if lo == 1: break
        t = lo - 1

    return str(answer)

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

Java

public class Euler739 {
    static final long kMod = 1000000007L;

    static long mulMod(long a, long b) {
        return (a * b) % kMod;
    }

    static long modPow(long base, long exp) {
        long result = 1L;
        while (exp > 0L) {
            if ((exp & 1L) != 0L) {
                result = mulMod(result, base);
            }
            base = mulMod(base, base);
            exp >>= 1L;
        }
        return result;
    }

    static class Pair {
        long a, b;

        Pair(long a, long b) {
            this.a = a;
            this.b = b;
        }
    }

    static Pair fibPair(long n) {
        if (n == 0L) {
            return new Pair(0L, 1L);
        }
        Pair p = fibPair(n >> 1L);
        long a = p.a;
        long b = p.b;

        long twoB = (2L * b) % kMod;
        long twoBMinusA = (twoB + kMod - a) % kMod;
        long c = mulMod(a, twoBMinusA);
        long d = (mulMod(a, a) + mulMod(b, b)) % kMod;

        if ((n & 1L) == 0L) {
            return new Pair(c, d);
        }
        return new Pair(d, (c + d) % kMod);
    }

    static long lucas(long n) {
        if (n == 0L) {
            return 2L;
        }
        Pair p = fibPair(n);
        long fn = p.a;
        long fn1 = p.b;
        long value = (2L * fn1) % kMod;
        value = (value + kMod - fn) % kMod;
        return value;
    }

    static void batchInvert(int[] values, int[] inverses, int[] prefix, int m) {
        if (m == 0)
            return;

        prefix[0] = values[0];
        for (int i = 1; i < m; ++i) {
            prefix[i] = (int) mulMod(prefix[i - 1], values[i]);
        }

        long suffixInv = modPow(prefix[m - 1], kMod - 2L);
        for (int i = m - 1; i >= 0; --i) {
            long left = (i == 0) ? 1L : prefix[i - 1];
            inverses[i] = (int) mulMod(suffixInv, left);
            suffixInv = mulMod(suffixInv, values[i]);
        }
    }

    public static String solve() {
        long n = 100000000L;
        long r = n - 1L;

        long t = r;
        long coeff = 1L;

        long lT = lucas(r);
        long lTPlus1 = lucas(r + 1L);

        long answer = 0L;

        int kBlock = 1000000;
        int[] denoms = new int[kBlock];
        int[] invDenoms = new int[kBlock];
        int[] prefix = new int[kBlock];

        while (t >= 1L) {
            long hi = t;
            long lo = (hi > kBlock) ? (hi - kBlock + 1L) : 1L;
            int m = (int) (hi - lo + 1L);

            for (int i = 0; i < m; ++i) {
                long curT = hi - i;
                if (curT == 1L) {
                    denoms[i] = 1;
                } else {
                    long u = r - curT + 1L;
                    denoms[i] = (int) mulMod(curT % kMod, u % kMod);
                }
            }

            batchInvert(denoms, invDenoms, prefix, m);

            for (int i = 0; i < m; ++i) {
                long curT = hi - i;

                answer += mulMod(coeff, lTPlus1);
                if (answer >= kMod) {
                    answer -= kMod;
                }

                if (curT == 1L) {
                    break;
                }

                long numA = curT - 1L;
                long numB = 2L * r - curT;
                long num = mulMod(numA % kMod, numB % kMod);

                coeff = mulMod(coeff, num);
                coeff = mulMod(coeff, invDenoms[i]);

                long nextLTPlus1 = lT;
                long nextLT = (lTPlus1 + kMod - lT) % kMod;
                lTPlus1 = nextLTPlus1;
                lT = nextLT;
            }

            if (lo == 1L) {
                break;
            }
            t = lo - 1L;
        }

        return Long.toString(answer);
    }

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