求 \(\sum_{x=1}^{n} x^m\) 的公式， \(n \leq 20\) 。
先看如下柿子：
通过这三个柿子可以大概发现规律：
\(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)\) ，比较系数可以得到如下方程组：
化简得
转化为矩阵乘法的形式：
高斯消元求解即可。时间复杂度 \(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;
}