[关闭]
@ivorysi 2018-01-02T08:57:28.000000Z 字数 11440 阅读 729

FFT学习笔记

笔记


具体证明还是在纸上写吧,在这写太麻烦了

BZOJ 3527: [Zjoi2014]力

Description

给出n个数qi,给出Fj的定义如下:

令Ei=Fi/qi,求Ei.

Input

第一行一个整数n。
接下来n行每行输入一个数,第i行表示qi。
n≤100000,0

Output

n行,第i行输出Ei。与标准答案误差不超过1e-2即可。

Sample Input

5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

Sample Output

-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

题解




我们发现实际上就是

对应位置的系数
减去

相反位置的系数
用FFT优化卷积

有一点没想通就是为什么他们都把数组长度做n*2然后再扩展次数界
才反应过来……乘法后的数字长度不是最多往后扩展两倍吗……多项式乘法最大不就最高次项指数翻一番吗……所以次数界最小不就要比n*2还要大吗

代码

  1. #include <iostream>
  2. #include <cstdio>
  3. #include <algorithm>
  4. #include <cmath>
  5. #define siji(i,x,y) for(int i=(x);i<=(y);++i)
  6. #define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
  7. #define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
  8. #define MAXN 400005
  9. #define ivorysi
  10. using namespace std;
  11. typedef long long ll;
  12. struct complex {
  13. double r,i;
  14. complex(double real=0.0,double image=0.0) {
  15. r=real;i=image;
  16. }
  17. complex operator + (const complex &rhs) {
  18. return complex(r+rhs.r,i+rhs.i);
  19. }
  20. complex operator - (const complex &rhs) {
  21. return complex(r-rhs.r,i-rhs.i);
  22. }
  23. complex operator * (const complex &rhs) {
  24. return complex(r*rhs.r-i*rhs.i,r*rhs.i+i*rhs.r);
  25. }
  26. };
  27. const double PI=acos(-1.0);
  28. void brc(complex *p,int l) {
  29. int j=l/2;
  30. for(int i=1;i<l-1;++i) {
  31. if(i<j) swap(p[i],p[j]);
  32. int k=l/2;
  33. while(j>=k) {
  34. j-=k;
  35. k>>=1;
  36. }
  37. j+=k;
  38. }
  39. }
  40. void FFT(complex *p,int l,int on) {
  41. brc(p,l);
  42. complex u,t;
  43. for(int h=2;h<=l;h<<=1) {
  44. complex wn=complex(cos(on*2*PI/h),sin(on*2*PI/h));
  45. for(int k=0;k<l;k+=h) {
  46. complex w=complex(1.0,0);
  47. for(int j=k;j<k+h/2;++j) {
  48. u=p[j];
  49. t=w*p[j+h/2];
  50. p[j]=u+t;
  51. p[j+h/2]=u-t;
  52. w=w*wn;
  53. }
  54. }
  55. }
  56. }
  57. int n;
  58. complex f[MAXN],f_rev[MAXN],g[MAXN],ans1[MAXN],ans2[MAXN];
  59. int main() {
  60. #ifdef ivorysi
  61. freopen("f1.in","r",stdin);
  62. #endif
  63. scanf("%d",&n);--n;
  64. siji(i,0,n) {
  65. scanf("%lf",&f[i].r);
  66. f_rev[n-i]=f[i];
  67. }
  68. siji(i,1,n) {
  69. g[i].r=1.0/i/i;
  70. }
  71. int l=1;while(l-1<n*2) {l<<=1;}
  72. FFT(f,l,1);FFT(f_rev,l,1);FFT(g,l,1);
  73. xiaosiji(i,0,l) {
  74. ans1[i]=f[i]*g[i];
  75. ans2[i]=f_rev[i]*g[i];
  76. }
  77. FFT(ans1,l,-1);FFT(ans2,l,-1);
  78. xiaosiji(i,0,l) ans1[i].r/=l,ans2[i].r/=l;
  79. siji(i,0,n) {
  80. printf("%.3lf\n",ans1[i].r-ans2[n-i].r);
  81. }
  82. return 0;
  83. }

BZOJ 3771: Triple

Description

我们讲一个悲伤的故事。
从前有一个贫穷的樵夫在河边砍柴。
这时候河里出现了一个水神,夺过了他的斧头,说:
“这把斧头,是不是你的?”
樵夫一看:“是啊是啊!”
水神把斧头扔在一边,又拿起一个东西问:
“这把斧头,是不是你的?”
樵夫看不清楚,但又怕真的是自己的斧头,只好又答:“是啊是啊!”
水神又把手上的东西扔在一边,拿起第三个东西问:
“这把斧头,是不是你的?”
樵夫还是看不清楚,但是他觉得再这样下去他就没法砍柴了。
于是他又一次答:“是啊是啊!真的是!”
水神看着他,哈哈大笑道:
“你看看你现在的样子,真是丑陋!”
之后就消失了。
樵夫觉得很坑爹,他今天不仅没有砍到柴,还丢了一把斧头给那个水神。
于是他准备回家换一把斧头。
回家之后他才发现真正坑爹的事情才刚开始。
水神拿着的的确是他的斧头。
但是不一定是他拿出去的那把,还有可能是水神不知道怎么偷偷从他家里拿走的。
换句话说,水神可能拿走了他的一把,两把或者三把斧头。
樵夫觉得今天真是倒霉透了,但不管怎么样日子还得过。
他想统计他的损失。
樵夫的每一把斧头都有一个价值,不同斧头的价值不同。总损失就是丢掉的斧头价值和。
他想对于每个可能的总损失,计算有几种可能的方案。
注意:如果水神拿走了两把斧头a和b,(a,b)和(b,a)视为一种方案。拿走三把斧头时,(a,b,c),(b,c,a),(c,a,b),(c,b,a),(b,a,c),(a,c,b)视为一种方案。

Input

第一行是整数N,表示有N把斧头。
接下来n行升序输入N个数字Ai,表示每把斧头的价值。

Output

若干行,按升序对于所有可能的总损失输出一行x y,x为损失值,y为方案数。

Sample Input

4
4
5
6
7

Sample Output

4 1
5 1
6 1
7 1
9 1
10 1
11 2
12 1
13 1
15 1
16 1
17 1
18 1

样例解释

11有两种方案是4+7和5+6,其他损失值都有唯一方案,例如4=4,5=5,10=4+6,18=5+6+7.

HINT

所有数据满足:Ai<=40000

题解

我FFT写出了一个很蠢的错误!把h写成了l……
这道题用fft,就是复杂化加法,因为我们知道,两个数相乘,指数相加
然后用a表示一个数的多项式,b表示两个相同的数搭配起来的值的多项式,c表示3个相同的数搭配出来的多项式
两把斧子的搭配方式 (a*a-b)/2
三把斧子搭配方式 (a*a*a-a*b*3+c*2)/6
用的是容斥原理

【打倒资本主义河神】

代码

  1. #include <iostream>
  2. #include <cstdio>
  3. #include <algorithm>
  4. #include <cmath>
  5. #define siji(i,x,y) for(int i=(x);i<=(y);++i)
  6. #define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
  7. #define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
  8. #define MAXN 40005
  9. //#define ivorysi
  10. using namespace std;
  11. typedef long long ll;
  12. struct complex {
  13. double r,i;
  14. complex(double real=0.0,double image=0.0) {
  15. r=real;i=image;
  16. }
  17. friend complex operator +(const complex &a,const complex &b) {
  18. return complex(a.r+b.r,a.i+b.i);
  19. }
  20. friend complex operator -(const complex &a,const complex &b) {
  21. return complex(a.r-b.r,a.i-b.i);
  22. }
  23. friend complex operator *(const complex &a,const complex &b) {
  24. return complex(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);
  25. }
  26. }a[MAXN*4],b[MAXN*4],c[MAXN*4],e2[MAXN*4],e3[MAXN*4];
  27. const double PI=acos(-1.0);
  28. ll ans[MAXN*4];
  29. int tot,n;
  30. void brc(complex *p,int l) {
  31. int j=l/2;
  32. xiaosiji(i,1,l-1) {
  33. if(i<j) swap(p[i],p[j]);
  34. int k=l/2;
  35. while(j>=k) {
  36. j-=k;
  37. k>>=1;
  38. }
  39. j+=k;
  40. }
  41. }
  42. void FFT(complex *p,int l,int on) {
  43. brc(p,l);
  44. complex u,t;
  45. for(int h=2;h<=l;h<<=1) {
  46. complex wn(cos(on*2*PI/h),sin(on*2*PI/h));
  47. for(int k=0;k<l;k+=h) {
  48. complex w(1.0,0);
  49. for(int j=k;j<k+h/2;++j) {
  50. u=p[j];
  51. t=w*p[j+h/2];
  52. p[j]=u+t;
  53. p[j+h/2]=u-t;
  54. w=w*wn;
  55. }
  56. }
  57. }
  58. if(on==-1) {
  59. xiaosiji(i,0,l) p[i].r/=l;
  60. }
  61. }
  62. int main() {
  63. #ifdef ivorysi
  64. freopen("f1.in","r",stdin);
  65. #endif
  66. scanf("%d",&n);
  67. int x;
  68. siji(i,1,n) {
  69. scanf("%d",&x);
  70. ans[x]++;
  71. tot=max(tot,x);
  72. a[x].r=1.0;b[2*x].r=1.0;c[3*x].r=1.0;
  73. }
  74. int l=1;while(l-1<3*tot) l<<=1;
  75. FFT(a,l,1);FFT(b,l,1);FFT(c,l,1);
  76. xiaosiji(i,0,l) {
  77. e2[i]=(a[i]*a[i]-b[i])*complex(1.0/2.0,0);
  78. e3[i]=(a[i]*a[i]*a[i]-a[i]*b[i]*complex(3,0)+c[i]*complex(2,0))
  79. *complex(1.0/6.0,0);
  80. }
  81. FFT(e2,l,-1);FFT(e3,l,-1);
  82. xiaosiji(i,0,l) {
  83. ans[i]+=(ll)(e2[i].r+0.5);
  84. ans[i]+=(ll)(e3[i].r+0.5);
  85. if(ans[i]!=0) {
  86. printf("%d %d\n",i,ans[i]);
  87. }
  88. }
  89. return 0;
  90. }

BZOJ 3509: [CodeChef] COUNTARI

Description

给定一个长度为N的数组A[],求有多少对i, j, k(1<=i

Input

第一行一个整数N(N<=10^5)。
接下来一行N个数A[i](A[i]<=30000)。

Output

一行一个整数。

Sample Input

10
3 5 3 6 3 4 10 4 5 2

Sample Output

9

题解

fft运算的数组开四倍!
这道题可以很容易想到一个用FFT暴力的思路,枚举中间点,两边用卷积快速算一下,复杂度是显然超时
然后很神奇的是我们可以用空间换时间,这种方法就是分块
我们把这n个数分成个为一块,一共有块,额外开数组记下每个块中数出现的次数,对于一个点,在如果i,k在块内都可以通过预处理的数值计算出来
例如当前块的编号是id,我们要记录块[1,id-1]之间每个数字出现的个数,块[id+1,tot]之间每个数的个数,然后块里排在它之前每个数出现的个数,当前数字编号为p,枚举块内除了p以外所有点进行转移
这样有什么好处呢,假如i和k都不在块内我们才需要卷积,这样可以降低卷积的次数,也就是说,卷积的次数减少到了块的个数
那么复杂度就变成了
这里最影响复杂度的算是了,大概可以支持S到2000左右,如果S是根号n的话会超时

代码

  1. #include <iostream>
  2. #include <cstdio>
  3. #include <algorithm>
  4. #include <cmath>
  5. #define siji(i,x,y) for(int i=(x);i<=(y);++i)
  6. #define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
  7. #define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
  8. #define MAXN 100005
  9. //#define ivorysi
  10. using namespace std;
  11. typedef long long ll;
  12. struct complex {
  13. double r,i;
  14. complex(double real=0.0,double image=0.0) {
  15. r=real;i=image;
  16. }
  17. friend complex operator +(const complex &a,const complex &b) {
  18. return complex(a.r+b.r,a.i+b.i);
  19. }
  20. friend complex operator -(const complex &a,const complex &b) {
  21. return complex(a.r-b.r,a.i-b.i);
  22. }
  23. friend complex operator *(const complex &a,const complex &b) {
  24. return complex(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);
  25. }
  26. }Lnum[200005],Rnum[200005],match[200005];
  27. const double PI=acos(-1.0);
  28. void brc(complex *p,int l) {
  29. int j=l/2;
  30. xiaosiji(i,1,l-1) {
  31. if(i<j) swap(p[i],p[j]);
  32. int k=l/2;
  33. while(j>=k) {
  34. j-=k;
  35. k>>=1;
  36. }
  37. j+=k;
  38. }
  39. }
  40. void FFT(complex *p,int l,int on) {
  41. brc(p,l);
  42. complex u,t;
  43. for(int h=2;h<=l;h<<=1) {
  44. complex wn(cos(on*2*PI/h),sin(on*2*PI/h));
  45. for(int k=0;k<l;k+=h) {
  46. complex w(1.0,0);
  47. for(int j=k;j<k+h/2;++j) {
  48. u=p[j];
  49. t=w*p[j+h/2];
  50. p[j]=u+t;
  51. p[j+h/2]=u-t;
  52. w=w*wn;
  53. }
  54. }
  55. }
  56. if(on==-1) {
  57. xiaosiji(i,0,l) p[i].r/=l;
  58. }
  59. }
  60. ll ans;
  61. int n,a[MAXN],lim,Lz;
  62. int R[200005],L[200005],cnt[60005][350],cnt_pre[60005];
  63. int bl[350],br[350],tot,S;
  64. void Block_Init() {
  65. S=1000;
  66. xiaosiji(i,0,n) {
  67. if(i%S==0) br[tot]=i-1,bl[++tot]=i;
  68. }
  69. br[tot]=n-1;
  70. siji(i,1,tot) {
  71. siji(j,bl[i],br[i]) {
  72. cnt[a[j]][i]++;
  73. }
  74. }
  75. }
  76. void calc(int p) {
  77. int id=p/S+1;
  78. siji(j,bl[id],p-1) {
  79. if(2*a[p]>=a[j]) ans+=R[2*a[p]-a[j]];
  80. }
  81. siji(j,p+1,br[id]) {
  82. if(2*a[p]>=a[j]) ans+=L[2*a[p]-a[j]]+cnt_pre[2*a[p]-a[j]];
  83. }
  84. if(p==0 || p==n-1) return;
  85. ans+=(ll)(match[2*a[p]].r+0.5);
  86. }
  87. int main() {
  88. #ifdef ivorysi
  89. freopen("f1.in","r",stdin);
  90. #endif
  91. scanf("%d",&n);
  92. xiaosiji(i,0,n) {
  93. scanf("%d",&a[i]);
  94. lim=max(lim,a[i]);
  95. R[a[i]]++;
  96. }
  97. Block_Init();
  98. Lz=1;while(Lz-1<2*lim) Lz<<=1;
  99. siji(i,1,tot) {
  100. siji(j,0,lim) {R[j]-=cnt[j][i];cnt_pre[j]=0;}
  101. if(i!=1 || i!=tot) {
  102. siji(j,0,Lz) {
  103. Lnum[j]=complex(L[j],0);
  104. Rnum[j]=complex(R[j],0);
  105. }
  106. FFT(Lnum,Lz,1);FFT(Rnum,Lz,1);
  107. xiaosiji(j,0,Lz) {
  108. match[j]=Lnum[j]*Rnum[j];
  109. }
  110. FFT(match,Lz,-1);
  111. }
  112. siji(j,bl[i],br[i]) {
  113. calc(j);
  114. ++cnt_pre[a[j]];
  115. }
  116. siji(j,0,lim) {L[j]+=cnt[j][i];}
  117. }
  118. printf("%lld\n",ans);
  119. return 0;
  120. }

BZOJ 3456: 城市规划

Description

刚刚解决完电力网络的问题, 阿狸又被领导的任务给难住了.
刚才说过, 阿狸的国家有n个城市, 现在国家需要在某些城市对之间建立一些贸易路线, 使得整个国家的任意两个城市都直接或间接的连通. 为了省钱, 每两个城市之间最多只能有一条直接的贸易路径. 对于两个建立路线的方案, 如果存在一个城市对, 在两个方案中是否建立路线不一样, 那么这两个方案就是不同的, 否则就是相同的. 现在你需要求出一共有多少不同的方案.
好了, 这就是困扰阿狸的问题. 换句话说, 你需要求出n个点的简单(无重边无自环)无向连通图数目.
由于这个数字可能非常大, 你只需要输出方案数mod 1004535809(479 * 2 ^ 21 + 1)即可.

Input

仅一行一个整数n(<=130000)

Output

仅一行一个整数, 为方案数 mod 1004535809.

Sample Input

3

Sample Output

4

题解

这道题的质数比较神奇是
然后2^n还正好比n要大

然后我们进行神奇的,多项式求逆元
首先求逆元的式子是

我们要求
假如我们求好了

然后去求

为什么模数换了还成立啊……因为当然是的约数啊
然后后两个式子相减

我们把两边平方

为什么平方还成立啊
因为平方以前里面式子的每个项前面的系数都是0
平方以后i和k-i必然有一个东西是小于
所以平方后还成立

然后式子两边乘上

然后用FNT加速一下这个东西就好了

然后再说FNT是个什么鬼
因为我们的FFT只能满足在复数领域随便的乱搞,万一有了取模呢……然后就会变得特别的烦人
对于这道题,模数是一个非常神奇的东西,也就是说是,它还是个质数
同时这个
假如2×n的次数界就是的话
然后对于质数我们求它的原根(这道题的原根是3),原根和单位复根比较像,转到很多个不同的值后,根据费马小定理,还能会到1
这样保证了,同时还使这一圈取值各不相同,和单位复根的性质是一样的
然后我们就把当成单位复根做FFT就好,实际上还是很简单的

把基础知识都说了一遍之后,我们再来推卷积
我们设是不考虑连通性的简单无向图的个数
然后设是考虑连通性后简单无向图的个数,也就是我们需要的答案


然后等式两边同时除以同时把代入

我们发现这个东西满足卷积

然后



然后
的逆元
然后

代码

  1. #include <iostream>
  2. #include <cstdio>
  3. #include <algorithm>
  4. #include <cstring>
  5. #include <vector>
  6. #include <cmath>
  7. #define siji(i,x,y) for(int i=(x);i<=(y);++i)
  8. #define gongzi(j,x,y) for(int j=(x);j>=(y);--j)
  9. #define xiaosiji(i,x,y) for(int i=(x);i<(y);++i)
  10. //#define ivorysi
  11. using namespace std;
  12. typedef long long ll;
  13. const int MAXN =1<<20;
  14. const ll MOD=(479<<21)+1;
  15. const int primitive_root=3;
  16. int eps_length;
  17. int n;
  18. struct node {
  19. ll x[MAXN];
  20. int deg;
  21. node(){
  22. deg=0;
  23. x[0]=0;
  24. }
  25. };
  26. ll inv_fac[MAXN],inv[MAXN];
  27. ll w[MAXN],inv_w[MAXN],X[MAXN],choose[MAXN],ans;
  28. ll fpow(ll x,ll c) {
  29. ll res=1,t=x;
  30. while(c) {
  31. if(c&1) res=res*t%MOD;
  32. t=t*t%MOD;
  33. c>>=1;
  34. }
  35. return res;
  36. }
  37. void init_w(int L) {
  38. ll base=fpow(primitive_root,(MOD-1)/L);
  39. ll inv_base=fpow(base,MOD-2);
  40. w[0]=1;inv_w[0]=1;
  41. xiaosiji(i,1,L) {
  42. w[i]=w[i-1]*base%MOD;
  43. inv_w[i]=inv_w[i-1]*inv_base%MOD;
  44. }
  45. eps_length=L;
  46. }
  47. void FNT(ll *p,int L,ll *eps,int on) {
  48. for(int i=1,j=L/2;i<L-1;++i) {
  49. if(i<j) swap(p[i],p[j]);
  50. int k=L/2;
  51. while(j>=k) j-=k,k>>=1;
  52. j+=k;
  53. }
  54. ll u,t;
  55. for(int h=2;h<=L;h<<=1) {
  56. for(int k=0;k<L;k+=h) {
  57. for(int j=k,q=0;j<k+h/2;++j,q+=eps_length/h) {
  58. u=p[j];
  59. t=eps[q]*p[j+h/2]%MOD;
  60. p[j]=(u+t)%MOD;
  61. p[j+h/2]=((u-t)%MOD+MOD)%MOD;
  62. }
  63. }
  64. }
  65. if(on==-1) {
  66. ll inv_L=fpow(L,MOD-2);
  67. xiaosiji(i,0,L) {
  68. p[i]=p[i]*inv_L%MOD;
  69. }
  70. }
  71. }
  72. void calc_Inverse(int n,node &A,node &B) {
  73. if(n==1) {
  74. B.x[0]=fpow(A.x[0],MOD-2);
  75. B.deg=1;
  76. return;
  77. }
  78. calc_Inverse((n+1)>>1,A,B);
  79. int L=1;while(L-1 < n*2) L<<=1;
  80. xiaosiji(i,0,n) X[i]=A.x[i];
  81. xiaosiji(i,n,L) X[i]=0;
  82. xiaosiji(i,B.deg,L) B.x[i]=0;
  83. FNT(X,L,w,1);FNT(B.x,L,w,1);
  84. xiaosiji(i,0,L) {
  85. B.x[i]=B.x[i]*(2-X[i]*B.x[i]%MOD+MOD)%MOD;
  86. }
  87. FNT(B.x,L,inv_w,-1);
  88. B.deg=n;
  89. }
  90. node A,B,C;
  91. int main() {
  92. #ifdef ivorysi
  93. freopen("f1.in","r",stdin);
  94. #endif
  95. scanf("%d",&n);
  96. int L=1;while(L-1 < (n+1)*2) L<<=1;
  97. inv[1]=1;inv_fac[1]=inv_fac[0]=1;
  98. siji(i,2,n) {
  99. inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD;
  100. inv_fac[i]=inv[i]*inv_fac[i-1]%MOD;
  101. }
  102. init_w(L);
  103. calc_Inverse(n,A,C);
  104. choose[0]=choose[1]=1;
  105. siji(i,2,n) {
  106. choose[i]=fpow(2,(ll)i*(i-1)/2%(MOD-1));
  107. }
  108. siji(i,0,n) {
  109. A.x[i]=choose[i]*inv_fac[i]%MOD;
  110. }
  111. siji(i,1,n) {
  112. B.x[i]=choose[i]*inv_fac[i-1]%MOD;
  113. }
  114. calc_Inverse(n+1,A,C);
  115. xiaosiji(i,C.deg,L) {
  116. C.x[i]=0;
  117. }
  118. FNT(C.x,L,w,1);
  119. FNT(B.x,L,w,1);
  120. xiaosiji(i,0,L) {
  121. C.x[i]=C.x[i]*B.x[i]%MOD;
  122. }
  123. FNT(C.x,L,inv_w,-1);
  124. ans=C.x[n]*fpow(inv_fac[n-1],MOD-2)%MOD;
  125. printf("%lld\n",ans);
  126. return 0;
  127. }
添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注