# 问题

计算对于给定的整数$k$，求$f(n, k)=\sum_{i=1}^n i^k$关于$n$的多项式表达式。例如

\begin{align} f(n, 0)&=n\\ f(n, 1)&=n(n+1)/2\\ f(n, 2)&=n(n+1)(2n+1)/6\\ f(n, 3)&=n^2(n+1)^2/4 \end{align}

# 差分的等价描述

基于差分，我们可以将$f(n, k)$转换为等价形式，初值条件和递推关系

\begin{align} f(1, k)&=1\\ \Delta f(n, k) &= f(n, k)-f(n-1, k)=n^k \end{align}

我们也可以通过递推关系定义初值$f(0, k)=0$来代替$f(1, k)=1$

# 拟设(Ansatz)

我们拟设$f(n, k)$有形如$\sum_{i=0}^{k+1} c_i^k n^i$的待定系数解。

代入初值条件$f(0, k)$，我们就可以得到 $$c_0^k=0$$

代入递推关系，可得 \begin{align} n^k&=\Delta f(n, k)\\ &=\sum_{i=0}^{k+1} c_i^k [n^i-(n-1)^i]\\ &=\sum_{i=1}^{k+1} c_i^k \sum_{j=0}^{i-1} C_i^j(-1)^{i-j}n^j\\ &=\sum_{i=1}^{k+1} c_i^k \sum_{j=1}^{i} C_i^{j-1}(-1)^{i-j+1}n^{j-1}\\ &= \delta_j^{k+1} n^{j-1} \end{align}

因此 $\delta_j^{k+1}=L_{ji}c_i^k=\big[\sum_{i=1}^{k+1} \sum_{j=1}^{i} C_i^{j-1}(-1)^{i-j+1}\big] c_i^k$

简写求和符号，令$i,j$求和都在$1,\ldots,k+1$范围，但是规定$i\geq j$，得到上三角矩阵：

$$L_{ji}=(-1)^{i-j+1}C_i^{j-1}$$

因此，只需要求解线性方程组$L_{ji}c_i^k=\delta_j^{k+1}$即可得到$c_i^k=L^{-1}_{ij}\delta_j^{k+1}$

# 具体计算

利用$C_{i}^j=C_{i-1}^{j}+C_{i-1}^{j-1}$可以进行递推运算

In [13]:

from fractions import Fraction
from linear_solver import solve_triangular, fractize

def LMatrix(k):
    N = arange(k+1)
    L = diag(N+1)
    L[0] = 1
    for i in range(2, k+1):
        L[1:i, i]=L[0:i-1, i-1]+L[1:i, i-1]
    L *= (-1)**(N[:, np.newaxis]+N[np.newaxis, :])
    return fractize(L)

def delta(k):
    d = zeros([k+1, 1])
    d[-1]=1
    return fractize(d)

def coef(k):
    return solve_triangular(LMatrix(k), delta(k))[1].flatten()

def poly(c):
    def _poly(n):
        return dot(c, cumprod(ones_like(c, dtype='int')*n))
    return _poly

In [73]:

for i in range(6):
    print(coef(i))

[1]
[1/2 1/2]
[1/6 1/2 1/3]
[0 1/4 1/2 1/4]
[-1/30 0 1/3 1/2 1/5]
[0 -1/12 0 5/12 1/2 1/6]

整数幂的求和

# 问题

# 差分的等价描述

# 拟设(Ansatz)

# 具体计算

Comments