[caioj1456][FFT][拆系数板子]累加
来源:互联网 发布:java弹出消息框 编辑:程序博客网 时间:2024/06/10 22:31
【题意】
给出一串数a[i],定义数组b经过累加变换a得到。
如a{1,2,3,4}
则第一个b是{1,3,6,10}
此后的b则为累加变换b得到
如b{1,3,6,10}
则第二个b{1,4,10,20}
我们想知道第k次变换后的b数组
【输入】
第1行,2个数N和K,中间用空格分隔,N表示数组的长度,K表示处理的次数(2 <= n <= 50000, 0 <= k <= 10^9, 0 <= a[i] <= 10^9)
【输出】
共N行,每行一个数,对应经过K次处理后的结果。输出结果 Mod 10^9 + 7。
【样例输入】
4 1
1
2
3
4
【样例输出】
1
3
6
10
【题解】
引入了新的内容,拆系数FFT(膜myy)
本题最关键的一点就是在于模数,因为模数过大,FFT可能会炸精度。所以需要拆掉系数
设模数为M,则每个数都可以表示为K*√M+B
设两个多项式分别为A(x)和B(x),其系数序列分别为ai和bi
则可以拆出来4个序列k[ai],k[bi],b[ai],b[bi]
分别对4个序列进行FFT操作然后合并起来就达到了要求啦
注意拆系数FFT的精度要求较高,需要long double!!!(坑点)
之后我们来看题
上图是第几次操作后ai对后面数的贡献
可以发现系数是一个斜着的杨辉三角,那么可以用组合数快速求出来
除法使用逆元来求
那么求出c数组,将给出的数列和c数组拆系数后做卷积运算即可
#include<cstdio>#include<cstring>#include<cstdlib>#include<algorithm>#include<cmath>using namespace std;#define double long double//***卡了我一天的地方,这精度得用long double!!!! typedef long long LL;const int N=50010;const LL mod=(1e9+7);const LL m=sqrt(mod)+5;const double PI=acos(-1.0);const int MAXN=(1<<17);struct Complex{ double r,i;//real imag Complex() {} Complex(double _r,double _i){r=_r;i=_i;} friend Complex operator + (const Complex &x,const Complex &y){return Complex(x.r+y.r,x.i+y.i);} friend Complex operator - (const Complex &x,const Complex &y){return Complex(x.r-y.r,x.i-y.i);} friend Complex operator * (const Complex &x,const Complex &y){return Complex(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r);} friend Complex operator * (const Complex &x,const double t){return Complex(x.r*t,x.i*t);}};int R[MAXN],L;Complex x1[MAXN],x2[MAXN],x3[MAXN],x4[MAXN];Complex x5[MAXN],x6[MAXN],x7[MAXN],x8[MAXN];void fft(Complex *y,int len,int on){ for(int i=0;i<len;i++)if(i<R[i])swap(y[i],y[R[i]]); for(int i=1;i<len;i<<=1) { Complex wn(cos(PI/i),sin(on*PI/i)); for(int j=0;j<len;j+=(i<<1)) { Complex w(1,0); for(int k=0;k<i;k++,w=w*wn) { Complex u=y[j+k]; Complex v=w*y[j+k+i]; y[j+k]=u+v; y[j+k+i]=u-v; } } } if(on==-1) { double t=1.0/len; for(int i=0;i<len;i++)y[i].r*=t;//保护精度 }}LL pow_mod(LL a,LL b)//a^b%mod{ LL ans=1;a%=mod; while(b) { if(b%2==1)ans=ans*a%mod; a=a*a%mod;b/=2; } return ans;}LL p[550000];int n,k,a[51000];int len;LL inv(int x){return pow_mod(x,mod-2);}void pre(){ p[0]=1; for(int i=1;i<len;i++)p[i]=p[i-1]*((i+k-1)%mod)%mod*inv(i)%mod; //p[i]=p[i-1]*(i+k-1)/i //inv(i)是i的逆元 乘i的逆元实质上就是除i }int main(){ scanf("%d%d",&n,&k); for(len=1;len<=n*2;len<<=1)L++; for(int i=0;i<len;i++)R[i]=(R[i>>1]>>1)|(i&1)<<(L-1); pre();//组合数 for(int i=1;i<=n;i++)scanf("%d",&a[i]); for(int i=0;i<=n;i++) { x1[i].r=a[i]/m; x2[i].r=a[i]%m; //a[i]=x1[i].r*m+x2.r x3[i].r=p[i]%mod/m; x4[i].r=p[i]%mod%m; //p[i]=x3[i].r*m+x4.r } //如上拆掉系数,分别做fft fft(x1,len,1);fft(x2,len,1);fft(x3,len,1);fft(x4,len,1); //卷积写为a*p=(x1*m+x2)*(x3*m+x4)=(x1*x3*m*m)+(x1*m*x4)+(x2*x3*m)+(x2*x4) //先把他们乘起来,m最后乘进去就行 for(int i=0;i<len;i++) { x5[i]=x1[i]*x3[i]; x6[i]=x1[i]*x4[i]; x7[i]=x2[i]*x3[i]; x8[i]=x2[i]*x4[i]; } fft(x5,len,-1);fft(x6,len,-1);fft(x7,len,-1);fft(x8,len,-1);//IFFT for(int i=1;i<=n;i++) { LL ans=0; ans+=(long long)(x5[i].r+0.5)*(m*m%mod)%mod; ans+=((long long)(x6[i].r+0.5)+(long long)(x7[i].r+0.5))%mod*m%mod; ans+=(long long)(x8[i].r+0.5); //套卷积式即可,为了防止爆精度我们多模几个 mod printf("%d\n",ans%mod); } return 0;}
阅读全文
0 0
- [caioj1456][FFT][拆系数板子]累加
- FFT黑科技(拆系数FFT)
- FFT 板子
- 快速傅里叶变换 FFT 板子
- [FFT] FFT的一些无聊板子题
- 任意模数FFT 板子
- 板子
- 板子
- 累加
- 累加
- 累加
- FFT
- "fft"
- FFT
- fft
- FFT
- fft
- FFT
- 【OI之路】02数论算法-2素数判断
- 重复值判断(堆排序的非递归使用) -- 算法小结
- 二叉树
- 详解C语言字节对齐
- java合并排序
- [caioj1456][FFT][拆系数板子]累加
- 轻松使用8266
- 存储方式
- 代码干货 | jdk1.6环境下struts2改spring boot方案
- Hadoop纹为分布式搭建及常见问题
- org.w3c.dom.Element调用问题
- 微信公众号开发(2)---ACCESS_TOKEN和用户信息,http调用工具类
- PhpStrom 对laravel的IDE实现
- 输入框输入数据格式合法性验证