Project Euler 330

发表于 2026-01-07 分类于 Project Euler 阅读次数：本文字数： 6k 阅读时长 ≈ 5 分钟

Project Euler 330

题目

Euler’s Number

An infinite sequence of real numbers \(a(n)\) is defined for all integers \(n\) as follows:

\[a(n) = \begin{cases} 1 &n \lt 0\\ \sum \limits_{i = 1}^{\infty}{\dfrac{a(n - i)}{i!}} &n \ge 0 \end{cases}\]

For example,

\[\begin{aligned} &a(0)=\frac{1}{1!}+\frac{1}{2!}+\frac{1}{3!}+\dots=e-1\\ &a(1)=\frac{1}{1!}+\frac{1}{2!}+\frac{1}{3!}+\dots=e-1\\ &a(2)=\frac{2e-3}{1!}+\frac{e-1}{2!}+\frac{1}{3!}+\dots=\frac{7}{2}e-6 \end{aligned}\]

with \(e = 2.7182818\dots\) being Euler’s constant.

It can be shown that \(a(n)\) is of the form \(\dfrac{A(n)e+B(n)}{n!}\) for integers \(A(n)\) and \(B(n)\).

For example \(a(10) = \dfrac{328161643 e - 652694486}{10!}\). Find \(A(10^9) + B(10^9)\) and give your answer mod \(77 777 777\).

解决方案

令\(N=10^9\)。当 \(n\ge 0\) 时，若 \(i>n\) 则 \(n-i<0\)，因此 \(a(n-i)=1\)。于是\(\displaystyle{a(n)=\sum_{i=1}^{n}\frac{a(n-i)}{i!}+\sum_{i=n+1}^{\infty}\frac{1}{i!}.}\)

而\(\displaystyle{\sum_{i=n+1}^{\infty}\frac{1}{i!}=\sum_{i=0}^{\infty}\frac{1}{i!}-\sum_{i=0}^{n}\frac{1}{i!}=e-\sum_{i=0}^{n}\frac{1}{i!}.}\)

所以对 \(n\ge 0\)，有：

\[ a(n)=\sum_{i=1}^{n}\frac{a(n-i)}{i!}+e-\sum_{i=0}^{n}\frac{1}{i!}. \]

定义\(b(n)=n!a(n).\)两边级数乘 \(n!\) 得\(\displaystyle{b(n)=\sum_{i=1}^{n}\frac{n!}{i!}a(n-i)+n!e-\sum_{i=0}^{n}\frac{n!}{i!}.}\)

注意 \(a(n-i)=\dfrac{b(n-i)}{(n-i)!}\)，因此\(\dfrac{n!}{i!}a(n-i)=\dfrac{n!}{i!(n-i)!}b(n-i)=\dbinom{n}{i}b(n-i).\)

于是得到关于\(b(n),n\ge 0\)的有限递推：

\[ b(n)=\sum_{i=1}^{n}\binom{n}{i}b(n-i)+n!e-\sum_{i=0}^{n}\frac{n!}{i!} . \]

并且 \(a(0)=e-1\Rightarrow b(0)=e-1\)。

题目保证\(a(n)=\dfrac{A(n)e+B(n)}{n!},(A(n),B(n)\in\mathbb Z),\)等价于\(b(n)=n!a(n)=A(n)e+B(n).\)将上面的 \(b(n)\) 递推按 \(e\) 的系数与常数项分别比较：

初值：\(b(0)=e-1\Rightarrow A(0)=1,\ B(0)=-1\)。
对 \(n\ge 1\)：

\[ A(n)=n!+\sum_{i=1}^{n}\binom{n}{i}A(n-i),B(n)=\sum_{i=1}^{n}\binom{n}{i}B(n-i)-\sum_{i=0}^{n}\frac{n!}{i!}. \]

定义\(C(n)=2A(n)+B(n).\)接下来可以证明\(C(n)=n!\)。

由上两式：

\[ \begin{aligned} C(n)&=2\left(n!+\sum_{i=1}^n\binom{n}{i}A(n-i)\right) +\left(\sum_{i=1}^n\binom{n}{i}B(n-i)-\sum_{i=0}^n\frac{n!}{i!}\right)\\ &=\left(2n!-\sum_{i=0}^n\frac{n!}{i!}\right) +\sum_{i=1}^n\binom{n}{i}\bigl(2A(n-i)+B(n-i)\bigr)\\ &=\left(2n!-n!-\sum_{i=1}^n\frac{n!}{i!}\right) +\sum_{i=1}^n\binom{n}{i}C(n-i)\\ &=n!-\sum_{i=1}^n\frac{n!}{i!}+\sum_{i=1}^n\binom{n}{i}C(n-i). \end{aligned} \]

若猜测 \(C(k)=k!\) 对所有 \(k<n\) 成立，则\(\displaystyle{\sum_{i=1}^n\binom{n}{i}C(n-i)=\sum_{i=1}^n\binom{n}{i}(n-i)! = \sum_{i=1}^n\frac{n!}{i!}}\)。

从而 \(C(n)=n!\)。并且有初值\(C(0)=2A(0)+B(0)=2\cdot 1+(-1)=1=0!.\)因此可由归纳得到：

\[ 2A(n)+B(n)=n!. \]

定义\(D(n)=A(n)+B(n).\)

由 \(2A(n)+B(n)=n!\) 可推出\(A(n)=n!-D(n), B(n)=2D(n)-n!.\)

所以题目所求即为\(D(N)=A(N)+B(N)\)。因此接下来只需计算 \(D(n)\)。

将 \(A,B\) 的递推相加：

\[ \begin{aligned} D(n) &=\sum_{i=1}^{n}\binom{n}{i}D(n-i)+n!-\sum_{i=0}^{n}\frac{n!}{i!}\\ &=\sum_{i=1}^{n}\binom{n}{i}D(n-i)-\sum_{i=1}^{n}\frac{n!}{i!}. \end{aligned} \]

把第二项换元：令 \(k=n-i\)，则\(\displaystyle{\sum_{i=1}^{n}\frac{n!}{i!}=\sum_{k=0}^{n-1}\frac{n!}{(n-k)!}=\sum_{k=0}^{n-1}\binom{n}{k}k!.}\)

同时\(\displaystyle{\sum_{i=1}^{n}\binom{n}{i}D(n-i)=\sum_{k=0}^{n-1}\binom{n}{k}D(k).}\)

于是得到：

\[ D(n)=\sum_{k=0}^{n-1}\binom{n}{k}\bigl(D(k)-k!\bigr), (n\ge 1). \]

其中\(D(0)=0\)。

接下来通过引入有序Bell数，化为一阶递推与显式求和

令\(\displaystyle{d_n=\frac{D(n)}{n!}, \mathcal D(x)=\sum_{n\ge 0}d_n x^n.}\)那么有\(\displaystyle{d_n=\sum_{k=0}^{n-1}\frac{D(k)-k!}{k!(n-k)!}=\sum_{k=0}^{n-1}\frac{d_k-1}{(n-k)!}=\sum_{m=1}^{n}\frac{d_{n-m}-1}{m!}(n\ge 1).}\)

将\(d_n\)对 \(n\ge 1\) 乘 \(x^n\) 求和，有

\[\begin{aligned} \sum_{n\ge 1} d_n x^n &= \mathcal{D}(x)-d_0\\ &=\sum_{n\ge 1}\left(\sum_{m=1}^{n}\frac{d_{n-m}}{m!}\right)x^n-\sum_{n\ge 1}\left(\sum_{m=1}^{n}\frac{1}{m!}\right)x^n\\ &=\sum_{m\ge 1}\sum_{k\ge 0}\frac{d_k}{m!}x^{k+m}-\sum_{n\ge 1}\left(\sum_{m=1}^{n}\frac{1}{m!}\right)x^n\\ &=\left(\sum_{k\ge 0} d_k x^k\right)\left(\sum_{m\ge 1}\frac{x^m}{m!}\right)-\sum_{n\ge 1}\left(\sum_{m=1}^{n}\frac{1}{m!}\right)x^n\\ &=\mathcal{D}(x)\left(\sum_{m\ge 1}\frac{x^m}{m!}\right)-\sum_{n\ge 1}\left(\sum_{m=1}^{n}\frac{1}{m!}\right)x^n\\ \end{aligned} \]

利用级数的两个结果\(\displaystyle{\sum_{m\ge 1}\frac{x^m}{m!}=e^x-1,\sum_{n\ge 0}\left(\sum_{m=1}^n\frac1{m!}\right)x^n=\frac{e^x-1}{1-x},}\)可得\(\mathcal{D}(x)-d_0=\mathcal{D}(x)(e^x-1)-\dfrac{e^x-1}{1-x}\).代入\(d_0=0\)，整理后得到：

\[ \mathcal D(x)=\frac{e^x-1}{(1-x)(e^x-2)}. \]

现在定义新序列\(X(n)=nD(n-1)-D(n)(n\ge 1).\)可见\(\dfrac{X(n)}{n!}=d_{n-1}-d_n,\)

因此其普通生成函数为

\[ \sum_{n\ge 1}\frac{X(n)}{n!}x^n=\sum_{n\ge 1}(d_{n-1}-d_n)x^n=x\mathcal D(x)-(\mathcal D(x)-d_0)=(x-1)\mathcal D(x)=\dfrac{1}{2-e^x}-1. \]

代入 \(\mathcal D(x)\)：\((x-1)\mathcal D(x)=\dfrac{e^x-1}{2-e^x}=\dfrac{1}{2-e^x}-1.\)

定义

\[ \frac1{2-e^x}=\sum_{n\ge 0}\frac{F(n)}{n!}x^n. \]

这就是著名的有序 Bell数。于是比较系数可得：\(X(n)=F(n)\ (n\ge 1),\)即\(D(n)=nD(n-1)-F(n)(n\ge 1), D(0)=0.\)继续展开得到显式求和：

\[ D(n)=-\sum_{k=1}^{n}\frac{n!}{k!}F(k). \]

由 \(\dfrac{1}{2-e^x}= \dfrac{1}{1-(e^x-1)}\)，将其视作几何级数可推出标准递推：\(F(0)=1,F(n)=\sum_{k=1}^{n}\dbinom{n}{k}F(n-k),(n\ge 1).\)

本题模数\(Q=77,777,777=7\cdot 11\cdot 73\cdot 101\cdot 137\)是互素素数乘积，因此可分别在每个素数模下计算 \(D(N)\)，最后用中国剩余定理合并。

在模素数 \(p\) 下，\(\dfrac{n!}{k!}=n(n-1)\cdots(k+1)\)是长度为 \(n-k\) 的连续整数乘积。若 \(n-k\ge p\)，则其中必包含 \(p\) 的倍数，从而\(\frac{n!}{k!}\equiv 0\pmod p.\)

因此当 \(n\ge p\) 时，和式只剩下 \(n-k\le p-1\) 的项。换元 \(t=n-k\)：\(\displaystyle{D(n)\equiv -\sum_{t=0}^{p-1} n^{\underline{t}}F(n-t)\pmod p,n^{\underline{t}}=n(n-1)\cdots(n-t+1).}\)

这里上界必须是 \(p-1\)（即允许 \(t=p-1\)），因为长度为 \(p-1\) 的连续区间未必含 \(p\) 的倍数，不能提前判零；只有 \(t\ge p\) 才能保证为 \(0\)。

利用 Stirling 数第二类 \(S(n,k)\) 的恒等式：

\[F(n)=\sum_{k=0}^{n}k!S(n,k),\]

\[k!S(n,k)=\sum_{j=0}^{k}(-1)^{k-j}\binom{k}{j}j^n.\]

在模素数 \(p\) 下，当 \(k\ge p\) 时 \(k!\equiv 0\pmod p\)，故

\[ F(n)\equiv \sum_{k=0}^{p-1}\sum_{j=0}^{k}(-1)^{k-j}\binom{k}{j}j^n\pmod p. \]

对任意 \(0\le j<p\)，当 \(n\ge 1\) 时有\(j^{n+p-1}\equiv j^n\pmod p.\)（\(j=0\) 时两边都是 \(0\)，\(j\neq 0\) 时由费马小定理 \(j^{p-1}\equiv 1\)）。因此，当\(n\ge 1\)时，有：

\[F(n+p-1)\equiv F(n)\pmod p.\]

所以只要预处理 \(F(0),F(1),\dots,F(p-1)\pmod p\)，即可回答所有大下标：

\[ F(n)\equiv \begin{cases} F(n) & n<p\\ F!\left(1+((n-1)\bmod (p-1))\right) & n\ge p,\ n\ge 1. \end{cases} \]

最终，对每个 \(p\in\{7,11,73,101,137\}\)：

用 \(O(p^2)\) 的递推预处理 \(F(0),\dots,F(p-1)\bmod p\)。
计算\(\displaystyle{D(N)\equiv -\sum_{t=0}^{p-1} (N)^{\underline{t}}F(N-t)\pmod p,}\)其中 falling factorial 逐项乘上 \((N-t)\) 即可 \(O(p)\) 完成。
得到 5 个同余，使用 CRT 逐步合并得到模 \(Q\) 的唯一解。

代码

from tools import CRT

N = 10**9

def fubini_mod_table(p: int):
    # Fubini numbers F(0..p-1) mod p
    F = [0] * p
    F[0] = 1 % p

    # binomial C(n,k) mod p for n,k <= p-1
    C = [[0] * p for _ in range(p)]
    for n in range(p):
        C[n][0] = 1
        C[n][n] = 1
        for k in range(1, n):
            C[n][k] = (C[n - 1][k - 1] + C[n - 1][k]) % p

    for n in range(1, p):
        s = 0
        for k in range(1, n + 1):
            s = (s + C[n][k] * F[n - k]) % p
        F[n] = s
    return F

def ordered_bell_mod_large(n: int, p: int, Ftab):
    # Exact for n < p; periodic with period (p-1) for n >= 1
    if n < p:
        return Ftab[n]
    r = (n - 1) % (p - 1) + 1  # in 1..p-1
    return Ftab[r]

def D_mod_prime(n: int, p: int) -> int:
    # D(n) = -sum_{k=1..n} (n!/k!) F(k)
    # mod p: only terms with t=n-k < p can survive, i.e. t=0..p-1
    Ftab = fubini_mod_table(p)

    max_t = min(p - 1, n - 1)
    falling = 1 % p             # n^{\underline{0}}
    s = 0

    for t in range(0, max_t + 1):
        k = n - t
        fk = ordered_bell_mod_large(k, p, Ftab)
        s = (s + falling * fk) % p
        falling = (falling * (k % p)) % p  # multiply by (n-t)

    return (-s) % p


primes = [7, 11, 73, 101, 137]
Q = 77777777
r = [D_mod_prime(N, p) for p in primes]
print(CRT(primes, r))