求 \(\sum_{x=1}^{n} x^m\) 的公式， \(n \leq 20\) 。
先看如下柿子：
\[1 + 2 + 3 + \cdots + n = \frac{n \times (n+1)}{2} = \frac{n^2}{2} + \frac{n}{2} \]\[1^2 + 2^2 + 3^2 + \cdots + n^2 = \frac{n \times (n+1) \times (2n + 1)}{6} = \frac{n^3}{3} + \frac{n^2}{2} + \frac{n}{6} \]\[1^3 + 2^3 + 3^3 + \cdots + n^3 = \frac{[n(n+1)]^2}{4} = \frac{n^4}{4} + \frac{n^3}{2} + \frac{n^2}{4} \]通过这三个柿子可以大概发现规律：
\(F(n,m) = \sum\limits_{x=1}^{n} x^m\) 的最高次为 \(m+1\) ，且每一项都可以表示为 \(n^t \times p\) 的形式。
所以设 \(F(n,m) = a_1n^{m+1} + a_2n^{m} + \cdots + a_{m+1}n\) 。
向上式的两边同时加上 \((n+1)^m\) ，得
\(F(n+1,m) = a_1n^{m+1} + a_2n^{m} + \cdots + a_{m+1}n + (n+1)^m\) 。
用二项式定理展开末项，得
\(\color{CornFlowerBlue}F(n+1,m) = a_1n^{m+1} + (C_m^m + a_2)n^m + (C_m^{m-1}+a_3)n^{m-1} + \cdots + (C_m^1 + a_{m+1})n + C_m^0\)
而 \(F(n+1,m)\) 又等于 \(a_1(n+1)^{m+1} + a_2(n+1)^m + \cdots + a_{m+1}(n+1)\)
每一项都用二项式定理展开，得
\(\color{Salmon}F(n+1,m) = C_{m+1}^{m+1}a_1n^{m+1} + (C_{m+1}^{m}a_1 + C_{m}^{m}a_2)n^m + (C_{m+1}^{m-1}a_1 + C_m^{m-1}a_2 + C_{m-1}^{m-1}a_3)n^{m-1} + \cdots + (C_{m+1}^0a_1 + C_m^0a_2 + \cdots + C_1^0a_{m+1})\)
上面两个带颜色的柿子都等于 \(F(n+1,m)\) ，比较系数可以得到如下方程组：
\[ \begin{cases} a_1 = C_{m+1}^{m+1}a_1 \\ C_m^m + a_2 = C_{m+1}^ma_1 + C_m^ma_2 \\ C_m^{m-1} + a_3 = C_{m+1}^{m-1}a_1 + C_m^{m-1}a_2 + C_{m-1}^{m-1}a_3 \\ \dots \\ C_m^0 = C_{m+1}^0a_1 + C_m^0a_2 + \cdots + C_1^0a_{m+1} \end{cases} \]化简得
\[\begin{cases} C_{m+1}^ma_1 = C_m^m \\ C_{m+1}^{m-1}a_1 + C_m^{m-1}a_2 = C_m^{m-1} \\ C_{m+1}^{m-2}a_1 + C_{m}^{m-2}a_2 + C_{m-1}^{m-2}a_3 = C_m^{m-2} \\ \dots \\C_{m+1}^0a_1 + C_m^0a_2 + C_{m-1}^0a_3 + \cdots + C_1^0a_{m+1} = C_m^0\end{cases} \]转化为矩阵乘法的形式：
\[\begin{bmatrix} C_{m+1}^m \\ C_{m+1}^{m-1} & C_m^{m-1} \\ C_{m+1}^{m-2} & C_m^{m-2} & C_{m-1}^{m-2} \\ \vdots & & & \ddots \\ C_{m+1}^0 & C_m^0 & C_{m-1}^0 & \cdots & C_1^0 \end{bmatrix} \begin{bmatrix} a_1 \\ a_2 \\ a_3 \\ \vdots \\a_{m+1} \end{bmatrix} = \begin{bmatrix} C_m^m \\ C_m^{m-1} \\ C_m^{m-2} \\ \vdots \\ C_m^0 \end{bmatrix} \]高斯消元求解即可。时间复杂度 \(O(n^3)\) 。
#include<cmath>
#include<cstdio>
#include<algorithm>
using namespace std;
const int M=310;
typedef __int128 LL;
void write(__int128 x){
if(!x) return ;
if(x<0) putchar('-'),x=-x;
write(x/10);
putchar(x%10+'0');
}
int m;
LL fra[M];
void calc(){
fra[0]=1;
for(int i=1;i<=20;i++) fra[i]=fra[i-1]*i;
}
int C(int i,int j){
return fra[i]/fra[j]/fra[i-j];
}
LL gcd(LL A,LL B){
if(!B) return A;
return gcd(B,A%B);
}
struct Frac{
LL A,B; // A / B
Frac(){A=0,B=1;}
Frac(LL _A,LL _B):A(_A),B(_B){}
friend Frac operator+(Frac a,Frac b){
LL m=a.B*b.B,z=a.A*b.B+b.A*a.B;
LL qwq=gcd(m,z);
return Frac(z/qwq,m/qwq);
}
friend Frac operator-(Frac a,Frac b){
LL m=a.B*b.B,z=a.A*b.B-b.A*a.B;
LL qwq=gcd(m,z);
return Frac(z/qwq,m/qwq);
}
friend Frac operator*(Frac a,Frac b){
LL m=a.B*b.B,z=a.A*b.A;
LL qwq=gcd(m,z);
return Frac(z/qwq,m/qwq);
}
friend Frac operator/(Frac a,Frac b){
LL m=a.B*b.A,z=a.A*b.B;
LL qwq=gcd(m,z);
return Frac(z/qwq,m/qwq);
}
void print(){
bool negative=false;
if((A<0)^(B<0)) negative=true;
if(negative) printf("-");
write(A<0?-A:A);
printf("/");
write(B<0?-B:B);
}
double val(){
return 1.*A/B;
}
}a[M][M];
void Gauss(){
for(int i=1;i<=m+1;i++){
int cur=i;
for(int j=i+1;j<=m+1;j++)
if(fabs(a[cur][i].val())<fabs(a[j][i].val())) cur=j;
swap(a[i],a[cur]);
if(a[i][i].val()==0) return printf("No Way\n"),void();
for(int j=1;j<=m+1;j++){
if(i==j) continue;
Frac p=a[j][i]/a[i][i];
for(int k=1;k<=m+2;k++) a[j][k]=a[j][k]-(a[i][k]*p);
}
}
}
int main(){
calc();
scanf("%d",&m);
for(int i=1;i<=m+1;i++){
for(int j=1;j<=i;j++) a[i][j]=Frac(C((m+1)-j+1,m-i+1),1);
a[i][m+2]=Frac(C(m,m-i+1),1);
}
Gauss();
bool flag=false;
for(int i=m+1;i>=1;i--){
if(a[i][m+2].val()==0) continue;
if(flag) printf("+");
flag=true;
printf("(");
a[i][m+2]=a[i][m+2]/a[i][i];
a[i][m+2].print();
printf("*(n");
if((m+1)-i+1!=1) printf("^%d",(m+1)-i+1);
printf("))");
}
return 0;
}