Problem 427: $n$-sequences

Problem 427 source code

C++

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

long long MOD = 1000000009;

long long power(long long base, long long exp) {
    long long res = 1;
    base %= MOD;
    while (exp > 0) {
        if (exp % 2 == 1) res = (res * base) % MOD;
        base = (base * base) % MOD;
        exp /= 2;
    }
    return res;
}

long long modInverse(long long n) {
    return power(n, MOD - 2);
}

vector<long long> fact;
vector<long long> invFact;
vector<long long> powN;
vector<long long> powN_1; 

void precompute(int n) {
    fact.resize(n + 1);
    invFact.resize(n + 1);
    powN.resize(n + 1);
    powN_1.resize(n + 1);

    fact[0] = 1;
    powN[0] = 1;
    powN_1[0] = 1;
    long long n_mod = n % MOD;
    long long n_1_mod = (MOD - (n - 1) % MOD) % MOD; 

    for (int i = 1; i <= n; i++) {
        fact[i] = (fact[i - 1] * i) % MOD;
        powN[i] = (powN[i - 1] * n_mod) % MOD;
        powN_1[i] = (powN_1[i - 1] * n_1_mod) % MOD;
    }

    invFact[n] = modInverse(fact[n]);
    for (int i = n - 1; i >= 0; i--) {
        invFact[i] = (invFact[i + 1] * (i + 1)) % MOD;
    }
}

long long nCr(int n, int r) {
    if (r < 0 || r > n) return 0;
    return fact[n] * invFact[r] % MOD * invFact[n - r] % MOD;
}

long long get_C(int M, int k, int n_val) {
    if (M < 0) return 0;
    long long sum = 0;
    int max_p = M / (k + 1);
    
    for (int p = 0; p <= max_p; ++p) {
        long long term = nCr(M - p * k, p);
        term = (term * powN[M - p * (k + 1)]) % MOD;
        term = (term * powN_1[p]) % MOD;
        sum = (sum + term) % MOD;
    }
    return sum;
}

long long solve_f(int n) {
    precompute(n);
    
    long long total_sequences = power(n, n + 1); 
    long long sum_Nk = 0;

    #pragma omp parallel for reduction(+:sum_Nk)
    for (int k = 1; k < n; ++k) {
        long long c_n_1 = get_C(n - 1, k, n);
        long long c_n_1_k = get_C(n - 1 - k, k, n);
        
        long long W_n = (n * c_n_1) % MOD;
        long long sub = (n * c_n_1_k) % MOD;
        
        W_n = (W_n - sub + MOD) % MOD;
        
        sum_Nk = (sum_Nk + W_n) % MOD;
    }

    long long ans = (total_sequences - sum_Nk + MOD) % MOD;
    return ans;
}

int main() {
    struct TestCase {
        int n;
        long long expected;
    };
    
    vector<TestCase> tests = {
        {3, 45},
        {7, 1403689},
        {11, 496890792} 
    };

    bool all_passed = true;
    for (const auto& t : tests) {
        long long result = solve_f(t.n);
        if (result != t.expected) {
            cout << "Test failed for n=" << t.n << ". Expected " << t.expected << ", got " << result << endl;
            all_passed = false;
        } else {
            cout << "Test passed for n=" << t.n << endl;
        }
    }

    if (all_passed) {
        int target_n = 7500000;
        cout << "Computing f(" << target_n << ")..." << endl;
        long long final_ans = solve_f(target_n);
        cout << "f(" << target_n << ") mod " << MOD << " = " << final_ans << endl;
    }

    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.stream.IntStream;

public class Euler427 {
    static final long MOD = 1000000009L;

    static long power(long base, long exp) {
        long res = 1;
        base %= MOD;
        if (base < 0)
            base += MOD;
        while (exp > 0) {
            if (exp % 2 == 1)
                res = (res * base) % MOD;
            base = (base * base) % MOD;
            exp /= 2;
        }
        return res;
    }

    static long modInverse(long n) {
        return power(n, MOD - 2);
    }

    static long[] fact;
    static long[] invFact;
    static long[] powN;
    static long[] powN1;

    static void precompute(int n) {
        fact = new long[n + 1];
        invFact = new long[n + 1];
        powN = new long[n + 1];
        powN1 = new long[n + 1];

        fact[0] = 1;
        powN[0] = 1;
        powN1[0] = 1;

        long nMod = n % MOD;
        long n1Mod = (MOD - (n - 1) % MOD) % MOD;

        for (int i = 1; i <= n; i++) {
            fact[i] = (fact[i - 1] * i) % MOD;
            powN[i] = (powN[i - 1] * nMod) % MOD;
            powN1[i] = (powN1[i - 1] * n1Mod) % MOD;
        }

        invFact[n] = modInverse(fact[n]);
        for (int i = n - 1; i >= 0; i--) {
            invFact[i] = (invFact[i + 1] * (i + 1)) % MOD;
        }
    }

    static long nCr(int n, int r) {
        if (r < 0 || r > n)
            return 0;
        return fact[n] * invFact[r] % MOD * invFact[n - r] % MOD;
    }

    static long getC(int M, int k, int n) {
        if (M < 0)
            return 0;
        long sum = 0;
        int maxP = M / (k + 1);

        for (int p = 0; p <= maxP; p++) {
            long term = nCr(M - p * k, p);
            term = (term * powN[M - p * (k + 1)]) % MOD;
            term = (term * powN1[p]) % MOD;
            sum = (sum + term) % MOD;
        }
        return sum;
    }

    static long solveF(int n) {
        precompute(n);

        long totalSequences = power(n, n + 1);

        long sumNk = IntStream.range(1, n)
                .parallel()
                .mapToLong(k -> {
                    long cN1 = getC(n - 1, k, n);
                    long cN1k = getC(n - 1 - k, k, n);

                    long Wn = ((long) n * cN1) % MOD;
                    long sub = ((long) n * cN1k) % MOD;

                    Wn = (Wn - sub + MOD) % MOD;
                    return Wn;
                })
                .reduce(0L, (a, b) -> (a + b) % MOD);

        long ans = (totalSequences - sumNk + MOD) % MOD;
        return ans;
    }

    public static String solve() {
        return Long.toString(solveF(7500000));
    }

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