首页 > 代码库 > 【BZOJ3453】XLkxc [拉格朗日插值法]
【BZOJ3453】XLkxc [拉格朗日插值法]
XLkxc
Time Limit: 20 Sec Memory Limit: 128 MB[Submit][Status][Discuss]
Description
给定 k,a,n,d,p
f(i)=1^k+2^k+3^k+......+i^k
g(x)=f(1)+f(2)+f(3)+....+f(x)
求(g(a)+g(a+d)+g(a+2d)+......+g(a+nd))mod p
Input
第一行数据组数,(保证小于6)
以下每行四个整数 k,a,n,d
Output
每行一个结果。
Sample Input
5
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
Sample Output
5
5
5
5
5
5
5
5
5
HINT
1<=k<=123
0<=a,n,d<=123456789
p==1234567891
0<=a,n,d<=123456789
p==1234567891
Main idea
给定k,a,n,d,求
Source
我们可以令
然后推一波式子,再令
那么显然有
然后我们通过若干次差分,发现g在差分k+3次时全为0,那么g就是一个k+2次多项式;f在差分k+5次时全为0,那么f就是一个k+4次多项式。
我们通过拉格朗日插值法插g,得到k+5个f的值,然后再插值f就可以得到答案了。
Code
1 #include<iostream> 2 #include<string> 3 #include<algorithm> 4 #include<cstdio> 5 #include<cstring> 6 #include<cstdlib> 7 #include<cmath> 8 using namespace std; 9 typedef long long s64;10 const int ONE=1001;11 const s64 MOD=1234567891;12 13 int T;14 int k,a,n,d;15 int g[ONE],f[ONE];16 int inv[ONE],U[ONE],Jc[ONE];17 int pre[ONE],suc[ONE];18 19 int get() 20 {21 int res=1,Q=1; char c;22 while( (c=getchar())<48 || c>57)23 if(c==‘-‘)Q=-1;24 if(Q) res=c-48; 25 while((c=getchar())>=48 && c<=57) 26 res=res*10+c-48; 27 return res*Q; 28 }29 30 int Quickpow(int a,int b)31 {32 int res=1;33 while(b)34 {35 if(b&1) res=(s64)res*a%MOD;36 a=(s64)a*a%MOD;37 b>>=1;38 }39 return res;40 }41 42 int P(int k,int i)43 {44 if((k-i)&1) return -1+MOD;45 return 1;46 }47 48 namespace First49 {50 void Deal_jc(int k)51 {52 Jc[0]=1;53 for(int i=1;i<=k;i++) Jc[i]=(s64)Jc[i-1]*i%MOD;54 }55 56 void Deal_inv(int k)57 {58 inv[0]=1; inv[k]=Quickpow(Jc[k],MOD-2);59 for(int i=k-1;i>=1;i--) inv[i]=(s64)inv[i+1]*(i+1)%MOD;60 }61 }62 63 int Final(int f[],int n,int k)64 {65 pre[0]=1; for(int i=1;i<=k;i++) pre[i]=(s64)pre[i-1] * (n-i+MOD) % MOD;66 suc[0]=1; for(int i=1;i<=k;i++) suc[i]=(s64)suc[i-1] * (s64)(n-k+i-1+MOD) % MOD;67 68 s64 Ans=0;69 for(int i=1;i<=k;i++)70 {71 int Up= (s64) pre[i-1]*suc[k-i] % MOD * f[i] % MOD;72 int Down= (s64) inv[i-1]*inv[k-i] % MOD;73 74 Ans=(s64)(Ans + (s64) Up*Down % MOD * P (k,i) %MOD) % MOD;75 }76 77 return Ans;78 }79 80 int main()81 { 82 First::Deal_jc(150); First::Deal_inv(150);83 T=get();84 while(T--)85 {86 k=get(); a=get(); n=get(); d=get();87 88 for(int i=0;i<=k+3;i++) g[i]=Quickpow(i,k);89 for(int i=1;i<=k+3;i++) g[i]=((s64)g[i-1]+g[i])%MOD;90 for(int i=1;i<=k+3;i++) g[i]=((s64)g[i-1]+g[i])%MOD;91 for(int i=0;i<=k+5;i++)92 f[i]=((s64)f[i-1]+Final(g,(a+(s64)i*d)%MOD,k+3)) % MOD;93 94 printf("%d\n",Final(f,n,k+5)%MOD);95 }96 }
【BZOJ3453】XLkxc [拉格朗日插值法]
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。