各类求自然数幂和方法

高斯消元

我们知道：

[sum_{i=1}^{n}i=frac{n(n+1)}{2} ]

以及：

[sum_{i=1}^{n}i^2=frac{n(n+1)(2n+1)}{6} ]

以及：

[sum_{i=1}^{n}i^3=(sum_{i=1}^{n}i)^2=(frac{n(n+1)}{2})^2 ]

那我们可以猜想，自然数的(k)次幂和对应的公式是一个次数为(k+1)的没有常数项的多项式（实际上也是的）。
证明吗，暂时不会。。。

However，我们可以拿这个猜想做题。

设这个(k+1)次的多项式(f(x)=sum_{i=1}^{k+1}a_ix^i)
利用待定系数法，我们只需要知道(k+2)对((x,f(x)))，列出方程组就能解出所有的(a_i)，从而就能代入更大的(x)求出(f(x))。

由于解方程组需要用到高斯消元算法，时间复杂度是(O(k^3))，在(kleq 100)的范围内还是能无压力解决的。

总结

时间复杂度：(O(k^3))
空间复杂度：(O(k^2))

由于高斯消元时要在模意义下做除法，对于模数不是质数的情况无法适应，而且时间复杂度难以接受，不是一种较常用的方法。

第二类斯特林数

分析

定义(S(n,m))表示(n)个有差别的球放入m个无差别的盒子中的方案数，要求盒子不能为空。

容易得到下面的递推式：

[S(n,m)=S(n-1,m-1)+mS(n-1,m) ]

考虑新加入的球，要么放在新的盒子里，要么放在之前的盒子里。因为球是有差别的，所以放在任意一个盒子里的方案都是不一样的，因此(S(n-1,m))要乘上一个(m)。

要用它解决自然数幂和问题，还是要用到第二类斯特林数的一个性质：

[a^k=sum_{i=0}^{k}S(k,i)i!C_a^i ]

这个性质还是很好解释的，我们可以把(a^k)当做(k)个有差别的球，放入(a)个有差别的盒子的方案数，盒子可以为空。
那么我们就枚举(i)个盒子被放满了，(S(k,i))只保证了球有差别，乘以(i!)相当于给盒子编号，令盒子也有差别，最后乘上一个(C_a^{i})表示在(a)个盒子中选(i)个的方案数。

那么就可以开始化自然数幂求和的式子：
(sum_{a=1}^{n}a^k)
(=sum_{a=1}^{n}sum_{i=0}^{k}S(k,i)i!C_a^i)
两个sigma没有关联，我们可以交换枚举顺序：
(=sum_{i=0}^{k}S(k,i)i!sum_{a=1}^{n}C_a^i)
由于(a<i)时(C_a^{i}=0)，又可以化成：
(=sum_{i=0}^{k}S(k,i)i!sum_{a=i}^{n}C_a^i)

继续化简需要用到一个性质：
(sum_{a=i}^{n}C_a^i=C_{n+1}^{i+1})
证明考虑运用组合数递推公式即：
(C_i^j=C_{i-1}^j+C_{i-1}^{j-1})
(C_{n+1}^{i+1})
(=C_n^i+C_n^{i+1})
(=C_n^i+C_{n-1}^i+C_{n-1}^{i+1})
(=C_n^i+C_{n-1}^i+C_{n-2}^i+C_{n-2}^{i+1})
继续化下去就会得到：
(=sum_{a=i}^{n}C_a^i)

性质就得证了，上面的式子就化简为：
(=sum_{i=0}^{k}S(k,i)i!C_{n+1}^{i+1})
组合数有点麻烦，我们展开为阶乘形式：
(=sum_{i=0}^{k}S(k,i)i!frac{(n+1)!}{(i+1)!(n-i)!})
拆开((i+1)!=i!*(i+1))：
(=sum_{i=0}^{k}S(k,i)frac{(n+1)!}{(i+1)(n-i)!})
发现(frac{(n+1)!}{(n-i)!})其实是(i+1)个连续整数相乘，其中必有一个是(i+1)的倍数，因此式子一定取整数，就不用考虑模数的问题了。

那么直接枚举这(i+1)个连续的整数，得到了时间复杂度为(O(k^2))的算法。计算斯特林数和求自然数幂的复杂度都是(O(k^2))，总复杂度就是(O(k^2))。

Code

附带分解乘法黑科技

#include <cstdio>
#include <cstring>

typedef long long ll;
const int N = 2007;

ll k, n, p, s[N][N];

ll multi(ll a, ll b)
{
	ll x1 = a / 1000000, x2 = a % 1000000, y1 = b / 1000000, y2 = b % 1000000;
	return (x1 * 1000000 % p * y1 % p * 1000000 % p + x1 * 1000000 % p * y2 % p + y1 * 1000000 % p * x2 % p + x2 * y2 % p) % p;
}

ll solve(ll n)
{
	if (n == 0) return 0;
	ll ret = 0;
	for (int i = 1; i <= k && i <= n; i++)
	{
		ll sum = s[k][i];
		for (ll j = n - i + 1; j <= n + 1; j++)
			if (j % (i + 1) == 0) sum = multi(sum, j / (i + 1));
			else sum = multi(sum, j);
		ret = (ret + sum) % p;
	}
	return ret;
}

int main()
{
	scanf("%lld%lld%lld%lld", &k, &n, &p);
	s[0][0] = 1;
	for (int i = 1; i <= k; i++)
		for (int j = 1; j <= i; j++)
			s[i][j] = (s[i - 1][j - 1] + multi(j, s[i - 1][j])) % p;
	printf("%lld
", solve(n) % p);
	return 0;
}

总结

时间复杂度：(O(k^2))
空间复杂度：(O(k^2))

这种做法由于不用除法而适用于模数为任意数的情况，但是求斯特林数复杂度是(O(k^2))的，当(k)较大时不再适用。

拉格朗日插值法

分析

考虑上面的高斯消元做法，我们用待定系数法求出多项式每一项的系数，然后再代入(x)求值。

但其实不必求出多项式每一项的系数也能代入(x)求值。

先给出拉格朗日插值法的公式：

设((x_i,y_i) iin [1,n+1])为(n)次多项式(f(x))经过的不同的(n+1)个点。
有：

[f(x)=sum_{i=1}^{n+1}y_i prod_{j=1,j≠i}^{n+1}frac{x-x_j}{x_i-x_j} ]

这个公式不难理解，我们考虑要求的(f(x))满足什么条件：

• (forall iin [1,n+1],f(x_i)=y_i)。（经过指定的(n+1)个点）
• (f(x))的次数为(n)。（是(n)次数多项式）

这样的多项式是唯一的，而上面公式所得的(f(x))明显是满足以上两个条件的，因此公式算出来的(f(x))就是我们要求的。

那么，就可以用拉格朗日插值法代替上面的高斯消元过程。但是求的(k+2)个点比较特殊，必须是((1,f(1)),(2,f(2)),...,(k+2,f(k+2)))，这样后面计算的时候就可以用阶乘优化到(O(klogk))，当然随便求(k+2)个点也是可以的，但复杂度就是(O(k^2+klogk))了。

总结

时间复杂度：(O(klogk))（快速幂）
空间复杂度：(O(k))

由于除法，模数必须是质数，适用范围相对较广。