Problem 517: A Real Recursion

View on Project Euler

Project Euler Problem 517 Solution

EulerSolve provides an optimized solution for Project Euler Problem 517, A Real Recursion, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a positive integer \(n\), let \(a=\sqrt{n}\) and define a recursion by $$g_a(x)=\begin{cases} 1, & x<a,\\ g_a(x-1)+g_a(x-a), & x\ge a. \end{cases}$$ The quantity of interest is \(G(n)=g_a(n)\) with \(a=\sqrt n\). Problem 517 asks for $$\sum_{\substack{10^7<p<10^7+10^4\\ p\text{ prime}}} G(p)\pmod{10^9+7}.$$ A naive recursive evaluation would branch explosively, so the solution turns the recursion into a direct combinatorial formula and then evaluates that formula for every prime in the interval. Mathematical Approach The key observation is that every recursive branch corresponds to a sequence of jumps of size \(1\) and \(\sqrt n\). Counting leaves of the recursion tree is therefore the same as counting valid ordered jump patterns. Step 1: Interpret the Recursion as Ordered Jumps Start with remainder \(x=n\). Every recursive call replaces \(x\) by either \(x-1\) or \(x-\sqrt n\), and the branch stops as soon as the new remainder is \(<\sqrt n\). Thus each leaf is an ordered sequence of two move types: $$1\text{-step}: x\mapsto x-1,\qquad \sqrt n\text{-step}: x\mapsto x-\sqrt n.$$ If a branch uses \(k\) steps of size \(\sqrt n\), then automatically $$k\sqrt n\le n,$$ so only $$0\le k\le \left\lfloor\sqrt n\right\rfloor$$ can contribute. Step 2: Record Where the \(\sqrt n\)-Steps Occur Fix a branch using exactly \(k\) steps of size \(\sqrt n\)....

Detailed mathematical approach

Problem Summary

For a positive integer \(n\), let \(a=\sqrt{n}\) and define a recursion by

$$g_a(x)=\begin{cases} 1, & x<a,\\ g_a(x-1)+g_a(x-a), & x\ge a. \end{cases}$$

The quantity of interest is \(G(n)=g_a(n)\) with \(a=\sqrt n\). Problem 517 asks for

$$\sum_{\substack{10^7<p<10^7+10^4\\ p\text{ prime}}} G(p)\pmod{10^9+7}.$$

A naive recursive evaluation would branch explosively, so the solution turns the recursion into a direct combinatorial formula and then evaluates that formula for every prime in the interval.

Mathematical Approach

The key observation is that every recursive branch corresponds to a sequence of jumps of size \(1\) and \(\sqrt n\). Counting leaves of the recursion tree is therefore the same as counting valid ordered jump patterns.

Step 1: Interpret the Recursion as Ordered Jumps

Start with remainder \(x=n\). Every recursive call replaces \(x\) by either \(x-1\) or \(x-\sqrt n\), and the branch stops as soon as the new remainder is \(<\sqrt n\).

Thus each leaf is an ordered sequence of two move types:

$$1\text{-step}: x\mapsto x-1,\qquad \sqrt n\text{-step}: x\mapsto x-\sqrt n.$$

If a branch uses \(k\) steps of size \(\sqrt n\), then automatically

$$k\sqrt n\le n,$$

so only

$$0\le k\le \left\lfloor\sqrt n\right\rfloor$$

can contribute.

Step 2: Record Where the \(\sqrt n\)-Steps Occur

Fix a branch using exactly \(k\) steps of size \(\sqrt n\). For \(r=1,\dots,k\), let \(s_r\) be the number of \(1\)-steps that have already occurred before the \(r\)-th \(\sqrt n\)-step is taken.

These numbers satisfy

$$0\le s_1\le s_2\le \cdots \le s_k,$$

because the total number of \(1\)-steps seen so far can only increase.

Before taking the \(r\)-th \(\sqrt n\)-step, the current remainder is

$$n-s_r-(r-1)\sqrt n.$$

That step is legal only if the recursion has not stopped yet, namely if this remainder is at least \(\sqrt n\). Therefore

$$n-s_r-(r-1)\sqrt n\ge \sqrt n,$$

so

$$s_r\le n-r\sqrt n.$$

Since \(s_r\) is an integer, this becomes

$$s_r\le \left\lfloor n-r\sqrt n\right\rfloor.$$

For prime \(n\), the number \(\sqrt n\) is irrational, hence \(r\sqrt n\notin \mathbb Z\) for every \(r\ge1\), and we may rewrite the bound as

$$s_r\le n-\left\lceil r\sqrt n\right\rceil.$$

Because these upper bounds decrease as \(r\) increases, the single strongest condition is the last one:

$$s_k\le T_k:=n-\left\lceil k\sqrt n\right\rceil.$$

Step 3: Turn the Branch into a Multiset Count

Once the nondecreasing \(k\)-tuple \((s_1,\dots,s_k)\) is fixed, the whole branch is determined:

take \(s_1\) steps of size \(1\), then one step of size \(\sqrt n\); then \(s_2-s_1\) more \(1\)-steps, then another \(\sqrt n\)-step; continue in the same way until the \(k\)-th \(\sqrt n\)-step; after that, only \(1\)-steps remain until the remainder drops below \(\sqrt n\).

Conversely, every leaf with exactly \(k\) steps of size \(\sqrt n\) produces exactly one such tuple. So the problem is now:

How many tuples satisfy

$$0\le s_1\le s_2\le \cdots \le s_k\le T_k\ ?$$

This is the classical count of multisets of size \(k\) chosen from \(T_k+1\) values. Hence

$$\#\{\text{branches with exactly }k\text{ }\sqrt n\text{-steps}\}=\binom{T_k+k}{k}.$$

Substituting \(T_k=n-\lceil k\sqrt n\rceil\) gives

$$\binom{n-\lceil k\sqrt n\rceil+k}{k}.$$

Step 4: Sum over All Possible Values of \(k\)

Adding the contributions of all admissible \(k\) yields

$$G(n)=\sum_{k=0}^{\lfloor\sqrt n\rfloor}\binom{n-\lceil k\sqrt n\rceil+k}{k}.$$

This is the closed form used by the implementations. It replaces an exponential recursion tree with only \(O(\sqrt n)\) binomial terms.

Notice that the term for \(k=0\) is

$$\binom{n}{0}=1,$$

which corresponds to the unique branch that uses only \(1\)-steps until the stopping condition is reached.

Step 5: Prepare the Formula for Modular Computation

The required modulus is

$$M=10^9+7,$$

which is prime. Therefore each binomial coefficient can be evaluated as

$$\binom{m}{k}\equiv m!\,(k!)^{-1}\,((m-k)!)^{-1}\pmod M.$$

After factorials and inverse factorials have been precomputed once, every summand of \(G(n)\) is obtained in constant time.

The overall Project Euler sum is then

$$\sum_{\substack{10^7<p<10^7+10^4\\ p\text{ prime}}}\ \sum_{k=0}^{\lfloor\sqrt p\rfloor}\binom{p-\lceil k\sqrt p\rceil+k}{k}\pmod M.$$

Worked Example: \(n=5\)

Here \(a=\sqrt5\approx 2.236\), so \(k\) can only be \(0\), \(1\), or \(2\).

For \(k=0\), the contribution is

$$\binom{5}{0}=1.$$

This is the single branch that repeatedly subtracts \(1\):

$$5\to 4\to 3\to 2,<\sqrt5.$$

For \(k=1\), we have

$$T_1=5-\lceil \sqrt5\rceil=5-3=2,$$

so \(s_1\) can be \(0\), \(1\), or \(2\). That gives the three branches

$$\sqrt5,1;\qquad 1,\sqrt5;\qquad 1,1,\sqrt5,$$

hence the contribution

$$\binom{2+1}{1}=3.$$

For \(k=2\),

$$T_2=5-\lceil 2\sqrt5\rceil=5-5=0,$$

so only \((s_1,s_2)=(0,0)\) is possible, corresponding to the branch

$$\sqrt5,\sqrt5,$$

and the contribution is

$$\binom{0+2}{2}=1.$$

Therefore

$$G(5)=1+3+1=5.$$

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline.

First, they precompute factorials and inverse factorials modulo \(10^9+7\) up to the largest integer that can appear in a binomial coefficient. This makes each later combination lookup constant time.

Next, they generate all primes in the interval \(10^7<p<10^7+10^4\) with a segmented sieve. Because the window width is only \(10^4\), this is much cheaper than sieving the entire range from \(1\) up to \(10^7+10^4\).

For each prime \(p\), the implementation loops over

$$0\le k\le \lfloor\sqrt p\rfloor$$

and evaluates the summand

$$\binom{p-\lceil k\sqrt p\rceil+k}{k}\pmod{10^9+7}.$$

The only delicate quantity is \(\lceil k\sqrt p\rceil\). Rather than simulating the original real recursion, the implementations compute this ceiling numerically from \(\sqrt{k^2p}\), then form the corresponding binomial term and add it to \(G(p)\).

Finally, all prime contributions are accumulated modulo \(10^9+7\) to produce the required answer.

Complexity Analysis

Let \(U=10^7+10^4-1\), let \(W=10^4\) be the width of the prime interval, and let \(P\) be the number of primes in that interval.

Precomputing factorials and inverse factorials up to \(U\) costs \(O(U)\) time and \(O(U)\) memory. The segmented sieve needs only the base primes up to \(\sqrt U\), so prime generation over the target window is roughly \(O(\sqrt U\log\log U+W\log\log U)\) time and \(O(W)\) extra memory.

For each prime \(p\), evaluating \(G(p)\) requires \(\lfloor\sqrt p\rfloor+1\) summands, so that phase costs \(O(\sqrt p)\) per prime. Overall the running time is therefore

$$O\!\left(U+P\sqrt U\right),$$

with memory dominated by the factorial tables:

$$O(U).$$

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=517
  2. Binomial coefficient: Wikipedia — Binomial coefficient
  3. Stars and bars / multiset counting: Wikipedia — Stars and bars
  4. Segmented sieve: Wikipedia — Segmented sieve
  5. Modular inverse and Fermat's little theorem: Wikipedia — Fermat's little theorem

Problem 517 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <functional>
#include <iostream>
#include <unordered_map>
#include <vector>

namespace {

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

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

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

std::vector<int> small_primes_up_to(const int n) {
    std::vector<bool> is_prime(static_cast<std::size_t>(n + 1), true);
    if (n >= 0) {
        is_prime[0] = false;
    }
    if (n >= 1) {
        is_prime[1] = false;
    }
    for (int p = 2; p * p <= n; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        for (int m = p * p; m <= n; m += p) {
            is_prime[static_cast<std::size_t>(m)] = false;
        }
    }
    std::vector<int> out;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) {
            out.push_back(i);
        }
    }
    return out;
}

std::vector<int> segmented_primes(const int low, const int high_exclusive) {
    if (high_exclusive <= low) {
        return {};
    }

    const int root = static_cast<int>(std::sqrt(static_cast<long double>(high_exclusive - 1))) + 1;
    const std::vector<int> base_primes = small_primes_up_to(root);

    const int size = high_exclusive - low;
    std::vector<bool> is_prime(static_cast<std::size_t>(size), true);

    for (const int p : base_primes) {
        i64 start = static_cast<i64>((low + p - 1) / p) * p;
        const i64 pp = static_cast<i64>(p) * p;
        if (start < pp) {
            start = pp;
        }
        for (i64 x = start; x < high_exclusive; x += p) {
            is_prime[static_cast<std::size_t>(x - low)] = false;
        }
    }

    if (low == 0) {
        is_prime[0] = false;
        if (size > 1) {
            is_prime[1] = false;
        }
    } else if (low == 1) {
        is_prime[0] = false;
    }

    std::vector<int> primes;
    for (int i = 0; i < size; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) {
            primes.push_back(low + i);
        }
    }
    return primes;
}

u64 ceil_k_sqrt_n(const int n, const int k, const long double root) {
    long double approx = std::floor(static_cast<long double>(k) * root);
    u64 x = static_cast<u64>(approx);
    const u128 target = static_cast<u128>(k) * static_cast<u128>(k) * static_cast<u128>(n);

    while (static_cast<u128>(x + 1ULL) * static_cast<u128>(x + 1ULL) <= target) {
        ++x;
    }
    while (static_cast<u128>(x) * static_cast<u128>(x) > target) {
        --x;
    }
    if (static_cast<u128>(x) * static_cast<u128>(x) == target) {
        return x;
    }
    return x + 1ULL;
}

u64 nCk_mod(const int n, const int k, const std::vector<u32>& fact, const std::vector<u32>& invfact) {
    if (k < 0 || k > n) {
        return 0ULL;
    }
    return static_cast<u64>(fact[static_cast<std::size_t>(n)]) *
           invfact[static_cast<std::size_t>(k)] % kMod *
           invfact[static_cast<std::size_t>(n - k)] % kMod;
}

u64 G_mod(const int n, const std::vector<u32>& fact, const std::vector<u32>& invfact) {
    const long double root = std::sqrt(static_cast<long double>(n));
    const int k_max = static_cast<int>(std::floor(root));

    u64 total = 0ULL;
    for (int k = 0; k <= k_max; ++k) {
        const u64 ceil_term = ceil_k_sqrt_n(n, k, root);
        const int m = static_cast<int>(static_cast<u64>(n) - ceil_term + static_cast<u64>(k));
        if (m < k) {
            continue;
        }
        total += nCk_mod(m, k, fact, invfact);
        total %= kMod;
    }
    return total;
}

u64 solve_sum_primes(const int low_exclusive, const int high_exclusive) {
    const int max_n = high_exclusive - 1;
    std::vector<u32> fact(static_cast<std::size_t>(max_n + 1), 1U);
    std::vector<u32> invfact(static_cast<std::size_t>(max_n + 1), 1U);

    for (int i = 1; i <= max_n; ++i) {
        fact[static_cast<std::size_t>(i)] =
            static_cast<u32>(static_cast<u64>(fact[static_cast<std::size_t>(i - 1)]) * i % kMod);
    }
    invfact[static_cast<std::size_t>(max_n)] =
        static_cast<u32>(mod_pow(fact[static_cast<std::size_t>(max_n)], kMod - 2ULL, kMod));
    for (int i = max_n; i >= 1; --i) {
        invfact[static_cast<std::size_t>(i - 1)] =
            static_cast<u32>(static_cast<u64>(invfact[static_cast<std::size_t>(i)]) * i % kMod);
    }

    const std::vector<int> primes = segmented_primes(low_exclusive + 1, high_exclusive);
    u64 total = 0ULL;
    for (const int p : primes) {
        total += G_mod(p, fact, invfact);
        total %= kMod;
    }
    return total;
}

u64 brute_G_small(const int n) {
    const long double a = std::sqrt(static_cast<long double>(n));
    std::unordered_map<u64, u64> memo;
    memo.reserve(2048);

    std::function<u64(int, int)> dfs = [&](const int i, const int j) -> u64 {
        const long double x =
            static_cast<long double>(n) - static_cast<long double>(i) - static_cast<long double>(j) * a;
        if (x < a - 1e-15L) {
            return 1ULL;
        }
        const u64 key = (static_cast<u64>(static_cast<u32>(i)) << 32) |
                        static_cast<u64>(static_cast<u32>(j));
        const auto it = memo.find(key);
        if (it != memo.end()) {
            return it->second;
        }
        const u64 value = dfs(i + 1, j) + dfs(i, j + 1);
        memo.emplace(key, value);
        return value;
    };

    return dfs(0, 0);
}

bool run_checkpoints() {
    constexpr int checkpoint_n = 90;
    std::vector<u32> fact(static_cast<std::size_t>(checkpoint_n + 1), 1U);
    std::vector<u32> invfact(static_cast<std::size_t>(checkpoint_n + 1), 1U);
    for (int i = 1; i <= checkpoint_n; ++i) {
        fact[static_cast<std::size_t>(i)] =
            static_cast<u32>(static_cast<u64>(fact[static_cast<std::size_t>(i - 1)]) * i % kMod);
    }
    invfact[static_cast<std::size_t>(checkpoint_n)] =
        static_cast<u32>(mod_pow(fact[static_cast<std::size_t>(checkpoint_n)], kMod - 2ULL, kMod));
    for (int i = checkpoint_n; i >= 1; --i) {
        invfact[static_cast<std::size_t>(i - 1)] =
            static_cast<u32>(static_cast<u64>(invfact[static_cast<std::size_t>(i)]) * i % kMod);
    }

    if (G_mod(90, fact, invfact) != 7'564'511ULL) {
        std::cerr << "Checkpoint failed: G(90)\n";
        return false;
    }
    for (int n = 2; n <= 30; ++n) {
        if (G_mod(n, fact, invfact) != brute_G_small(n) % kMod) {
            std::cerr << "Checkpoint failed: formula/bruteforce mismatch at n=" << n << '\n';
            return false;
        }
    }
    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    constexpr int low = 10'000'000;
    constexpr int high = 10'010'000;
    std::cout << solve_sum_primes(low, high) << '\n';
    return 0;
}

Python

import math

kMod = 1000000007

def small_primes_up_to(n):
    is_prime = [True] * (n + 1)
    if n >= 0: is_prime[0] = False
    if n >= 1: is_prime[1] = False
    for p in range(2, int(math.isqrt(n)) + 1):
        if is_prime[p]:
            for m in range(p * p, n + 1, p):
                is_prime[m] = False
    return [i for i in range(2, n + 1) if is_prime[i]]

def segmented_primes(low, high_exclusive):
    if high_exclusive <= low:
        return []
    root = math.isqrt(high_exclusive - 1) + 1
    base_primes = small_primes_up_to(root)
    
    size = high_exclusive - low
    is_prime = [True] * size
    
    for p in base_primes:
        start = (low + p - 1) // p * p
        if start < p * p:
            start = p * p
        for x in range(start, high_exclusive, p):
            is_prime[x - low] = False
            
    if low == 0:
        is_prime[0] = False
        if size > 1: is_prime[1] = False
    elif low == 1:
        is_prime[0] = False
        
    primes = []
    for i in range(size):
        if is_prime[i]:
            primes.append(low + i)
    return primes

def ceil_k_sqrt_n(n, k):
    target = k * k * n
    x = math.isqrt(target)
    if x * x == target:
        return x
    return x + 1

def nCk_mod(n, k, fact, invfact):
    if k < 0 or k > n:
        return 0
    return (fact[n] * invfact[k] % kMod) * invfact[n - k] % kMod

def G_mod(n, fact, invfact):
    k_max = math.isqrt(n)
    total = 0
    for k in range(k_max + 1):
        ceil_term = ceil_k_sqrt_n(n, k)
        m = n - ceil_term + k
        if m < k:
            continue
        total = (total + nCk_mod(m, k, fact, invfact)) % kMod
    return total

def solve_sum_primes(low_exclusive, high_exclusive):
    max_n = high_exclusive - 1
    fact = [1] * (max_n + 1)
    for i in range(1, max_n + 1):
        fact[i] = (fact[i - 1] * i) % kMod
        
    invfact = [1] * (max_n + 1)
    invfact[max_n] = pow(fact[max_n], kMod - 2, kMod)
    for i in range(max_n, 0, -1):
        invfact[i - 1] = (invfact[i] * i) % kMod
        
    primes = segmented_primes(low_exclusive + 1, high_exclusive)
    total = 0
    for p in primes:
        total = (total + G_mod(p, fact, invfact)) % kMod
    return total

def solve():
    low = 10000000
    high = 10010000
    ans = solve_sum_primes(low, high)
    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.List;

public class Euler517 {

    static final long kMod = 1000000007L;

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

    static List<Integer> smallPrimesUpTo(int n) {
        boolean[] isPrime = new boolean[n + 1];
        for (int i = 2; i <= n; i++)
            isPrime[i] = true;
        for (int p = 2; p * p <= n; p++) {
            if (isPrime[p]) {
                for (int m = p * p; m <= n; m += p) {
                    isPrime[m] = false;
                }
            }
        }
        List<Integer> out = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (isPrime[i])
                out.add(i);
        }
        return out;
    }

    static List<Integer> segmentedPrimes(int low, int highExclusive) {
        if (highExclusive <= low)
            return new ArrayList<>();
        int root = (int) Math.sqrt(highExclusive - 1) + 1;
        List<Integer> basePrimes = smallPrimesUpTo(root);

        int size = highExclusive - low;
        boolean[] isPrime = new boolean[size];
        for (int i = 0; i < size; i++)
            isPrime[i] = true;

        for (int p : basePrimes) {
            long start = ((long) low + p - 1) / p * p;
            long pp = (long) p * p;
            if (start < pp)
                start = pp;
            for (long x = start; x < highExclusive; x += p) {
                isPrime[(int) (x - low)] = false;
            }
        }

        if (low == 0) {
            isPrime[0] = false;
            if (size > 1)
                isPrime[1] = false;
        } else if (low == 1) {
            isPrime[0] = false;
        }

        List<Integer> primes = new ArrayList<>();
        for (int i = 0; i < size; i++) {
            if (isPrime[i])
                primes.add(low + i);
        }
        return primes;
    }

    static long ceilKSqrtN(int n, int k) {
        long target = (long) k * k * n;
        long x = (long) Math.sqrt(target);
        if (x * x == target)
            return x;
        return x + 1;
    }

    static long nCkMod(int n, int k, int[] fact, int[] invfact) {
        if (k < 0 || k > n)
            return 0;
        return ((long) fact[n] * invfact[k] % kMod) * invfact[n - k] % kMod;
    }

    static long gMod(int n, int[] fact, int[] invfact) {
        int kMax = (int) Math.sqrt(n);
        long total = 0;
        for (int k = 0; k <= kMax; k++) {
            long ceilTerm = ceilKSqrtN(n, k);
            int m = n - (int) ceilTerm + k;
            if (m < k)
                continue;
            total = (total + nCkMod(m, k, fact, invfact)) % kMod;
        }
        return total;
    }

    static long solveSumPrimes(int lowExclusive, int highExclusive) {
        int maxN = highExclusive - 1;
        int[] fact = new int[maxN + 1];
        int[] invfact = new int[maxN + 1];

        fact[0] = 1;
        for (int i = 1; i <= maxN; i++) {
            fact[i] = (int) (((long) fact[i - 1] * i) % kMod);
        }

        invfact[maxN] = (int) modPow(fact[maxN], kMod - 2, kMod);
        for (int i = maxN; i >= 1; i--) {
            invfact[i - 1] = (int) (((long) invfact[i] * i) % kMod);
        }

        List<Integer> primes = segmentedPrimes(lowExclusive + 1, highExclusive);
        long total = 0;
        for (int p : primes) {
            total = (total + gMod(p, fact, invfact)) % kMod;
        }
        return total;
    }

    public static void main(String[] args) {
        int low = 10000000;
        int high = 10010000;
        System.out.println(solveSumPrimes(low, high));
    }
}