[51nod 1258] [伯努利数] [多项式求逆] [任意模数NTT] 序列求和 V4
来源:互联网 发布:vscode 如何打开sln 编辑:程序博客网 时间:2024/06/06 09:44
接http://blog.csdn.net/coldef/article/details/57908865
上次做一套模拟赛的时候,其中需要求自然数k次幂和,然后我只会n^2的…我记得n^2有20分,nlogn求可以爆到90分……
——鏼鏼鏼2015年国家集训队论文
大概就这样多项式求个逆,求出生成函数就可以了
这样是
至于为什么要大于
不过直接上CRT要爆longlong,可以用一个小技巧处理
——图自AntiLeaf的博客
这样每次NTT分成三个素数的NTT,再合并,就可以了
#include <cstdio>#include <iostream>#include <algorithm>#include <string>#include <cstring>using namespace std;typedef long long ll;const int N=300010,P=1e9+7;const int p1=998244353,p2=1004535809,p3=469762049;inline char nc(){ static char buf[100000],*p1=buf,*p2=buf; return p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++;}inline void rea(int &x){ char c=nc(); x=0; for(;c>'9'||c<'0';c=nc());for(;c>='0'&&c<='9';x=x*10+c-'0',c=nc());}inline void rea(ll &x){ char c=nc(); x=0; for(;c>'9'||c<'0';c=nc());for(;c>='0'&&c<='9';x=x*10+c-'0',c=nc());}inline ll Pow(ll x,ll y,ll p){ x%=p; ll ret=1; for(;y;y>>=1,x=x*x%p) if(y&1) ret=ret*x%p; return ret;}const ll M_=1LL*p1*p2,inv1=Pow(p1,p2-2,p2),inv2=Pow(p2,p1-2,p1),inv3=Pow(M_%p3,p3-2,p3);int rev[N],fac[N],inv[N],a[N],b[N];struct NTT{ int P,num; int w[2][N]; void Pre(int p,int n){ P=p; num=n; int ng=Pow(3,(P-1)/num,P),ig=Pow(ng,P-2,P); w[1][0]=w[0][0]=1; for(int i=1;i<num;i++) w[1][i]=1LL*w[1][i-1]*ng%P,w[0][i]=1LL*w[0][i-1]*ig%P; } void FFT(int *a,int n,int r){ for(int i=1;i<n;i++) if(rev[i]>i) swap(a[i],a[rev[i]]); for(int i=1;i<n;i<<=1) for(int j=0;j<n;j+=(i<<1)) for(int k=0;k<i;k++){ int x=a[j+k],y=1LL*a[j+k+i]*w[r][num/(i<<1)*k]%P; a[j+k]=(x+y)%P; a[j+k+i]=(x+P-y)%P; } if(!r) for(int i=0,inv=Pow(n,P-2,P);i<n;i++) a[i]=1LL*a[i]*inv%P; }}ntt[3];inline void Pre(int n){ fac[0]=1; for(int i=1;i<=n;i++) fac[i]=1LL*fac[i-1]*i%P; inv[1]=1; for(int i=2;i<=n;i++) inv[i]=1LL*(P-P/i)*inv[P%i]%P; inv[0]=1; for(int i=1;i<=n;i++) inv[i]=1LL*inv[i]*inv[i-1]%P;}inline ll mul(ll a,ll b,ll p){ a%=p; b%=p; return ((a*b-(ll)((ll)((long double)a/p*b+1e-3)*p))%p+p)%p;}inline int CRT(int a1,int a2,int a3){ ll A=(mul(1LL*a1*p2%M_,inv2,M_)+mul(1LL*a2*p1%M_,inv1,M_))%M_; ll K=(1LL*a3+p3-A%p3)*inv3%p3; return (K*(M_%P)+A)%P;}int _b[3][N],__b[N];void Inv(int *a,int *b,int n){ static int tmp[N]; if(n==1) return void(b[0]=CRT(Pow(a[0],p1-2,p1),Pow(a[0],p2-2,p2),Pow(a[0],p3-2,p3))); Inv(a,b,n>>1); int L=0; while(!(n>>L&1)) L++; for(int i=1;i<(n<<1);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L); for(int k=0;k<3;k++) for(int i=0;i<n;i++) _b[k][i]=b[i],_b[k][i+n]=0; for(int k=0;k<3;k++){ for(int i=0;i<n;i++) tmp[i]=a[i],tmp[i+n]=0; ntt[k].FFT(tmp,n<<1,1); ntt[k].FFT(_b[k],n<<1,1); for(int i=0;i<(n<<1);i++) _b[k][i]=1LL*_b[k][i]*tmp[i]%ntt[k].P; ntt[k].FFT(_b[k],n<<1,0); } for(int i=0;i<n;i++) _b[0][i]=_b[1][i]=_b[2][i]=(P-CRT(_b[0][i],_b[1][i],_b[2][i]))%P; _b[0][0]=_b[1][0]=_b[2][0]=(_b[0][0]+2)%P; for(int k=0;k<3;k++){ for(int i=0;i<n;i++) tmp[i]=b[i],tmp[i+n]=0; ntt[k].FFT(tmp,n<<1,1); ntt[k].FFT(_b[k],n<<1,1); for(int i=0;i<(n<<1);i++) _b[k][i]=1LL*_b[k][i]*tmp[i]%ntt[k].P; ntt[k].FFT(_b[k],n<<1,0); } for(int i=0;i<n;i++) b[i]=CRT(_b[0][i],_b[1][i],_b[2][i]);}int t,k;ll n;inline int C(int x,int y){ return 1LL*fac[x]*inv[y]%P*inv[x-y]%P;}int main(){ int M=1; for(;M<=50000;M<<=1); ntt[0].Pre(p1,M<<1); ntt[1].Pre(p2,M<<1); ntt[2].Pre(p3,M<<1); Pre(M); for(int i=0;i<M;i++) a[i]=inv[i+1]; Inv(a,b,M); for(int i=0;i<M;i++) b[i]=1LL*b[i]*fac[i]%P; b[1]=P-b[1]; rea(t); while(t--){ rea(n); rea(k); n%=P; int ans=0; for(int i=k,prod=n;~i;prod=1LL*prod*n%P,i--) ans=(1LL*prod*b[i]%P*C(k+1,i)%P+ans)%P; ans=1LL*ans*Pow(k+1,P-2,P)%P; printf("%d\n",ans); } return 0;}
阅读全文
0 0
- [51nod 1258] [伯努利数] [多项式求逆] [任意模数NTT] 序列求和 V4
- [伯努利数] 51Nod 1258 序列求和 V4
- 51nod 1258 序列求和 V4 拉格朗日插值法求自然数幂和
- 51nod 1258 FFT 多项式求逆 伯努利数
- bzoj 3456: 城市规划 NTT+多项式求逆
- bzoj3625(NTT+多项式求逆+多项式开根)
- NesCafe41 异化多肽 NTT实现多项式求逆
- bzoj 3456: 城市规划 (NTT+多项式求逆)
- [BZOJ2259]异化多肽(生成函数+NTT+多项式求逆)
- [生成函数][NTT][多项式求逆]BZOJ 3456: 城市规划
- 51nod 1228 序列求和
- 51Nod-1228-序列求和
- 【BZOJ3625】【CF438E】小朋友和二叉树 NTT 生成函数 多项式开根 多项式求逆
- 51NOD 序列求和 V5题解
- 51nod 1172 Partial Sums V2 任意模数FFT
- 任意十个数求和
- [BZOJ3456]城市规划(组合数学+容斥原理+NTT+多项式求逆)
- bzoj 4555:[Tjoi2016&Heoi2016]求和 多项式求逆
- Mysql主从配置
- Android Studio下使用ndkBuild
- maven链接私服方法
- 微信公众号三方平台开发【component_access_token篇】
- 【JavaSE学习笔记】算术运算符
- [51nod 1258] [伯努利数] [多项式求逆] [任意模数NTT] 序列求和 V4
- YII2 advanced 高级版本项目搭建-添加API应用以及多应用(一)
- 用verilog实现AES密码算法1---一些理论准备
- Jack server 错误
- nrf51822-配对绑定实现过程
- C++11右值 &&引用
- Python 程序员需要知道的 30 个技巧
- centos彻底删除文件夹、文件命令
- 7