先简短几句话说说FFT....
多项式可用系数和点值表示,n个点可确定一个次数小于n的多项式。
多项式乘积为 f(x)*g(x),显然若已知f(x), g(x)的点值,O(n)可求得多项式乘积的点值。
我们所需要的就是O(nlogn)快速地将两个系数多项式表示成点值多项式,O(n)求得乘积的点值表示后O(nlogn)还原成系数多项式。
这里就需要套FFT板子了...
FFT中取n个单位根,需要n是2的幂。
又因为n个点可确定一个次数小于n的多项式,所以n > 乘积多项式的最高次数。
以上。
HDU4609 n个木棍任取三根能组成三角形的概率。
数组开小莫名T,WA.
1 #include <bits/stdc++.h>
2 #define ll long long
3 using namespace std;
4 const int N = 4e5+10;
5 struct comp{
6 double r,i;comp(double _r=0,double _i=0){r=_r;i=_i;}
7 comp operator+(const comp x){return comp(r+x.r,i+x.i);}
8 comp operator-(const comp x){return comp(r-x.r,i-x.i);}
9 comp operator*(const comp x){return comp(r*x.r-i*x.i,r*x.i+i*x.r);}
10 };
11 const double pi=acos(-1.0);
12 void FFT(comp a[],int n,int t){
13 for(int i=1,j=0;i<n-1;i++){
14 for(int s=n;j^=s>>=1,~j&s;);
15 if(i<j)swap(a[i],a[j]);
16 }
17 for(int d=0;(1<<d)<n;d++){
18 int m=1<<d,m2=m<<1;
19 double o=pi/m*t;comp _w(cos(o),sin(o));
20 for(int i=0;i<n;i+=m2){
21 comp w(1,0);
22 for(int j=0;j<m;j++){
23 comp &A=a[i+j+m],&B=a[i+j],t=w*A;
24 A=B-t;B=B+t;w=w*_w;
25 }
26 }
27 }
28 if(t==-1)for(int i=0;i<n;i++)a[i].r/=n;
29 }
30 comp x[N];
31 ll num[N], sum[N];
32 int main(){
33 int T; scanf("%d", &T);
34 while(T--){
35 memset(num, 0, sizeof num);
36 int n, u, maxnum = -1; scanf("%d", &n);
37 for(int i = 0; i < n; i++){
38 scanf("%d", &u);
39 maxnum = max(maxnum, u), num[u]++;
40 }
41 int len = 1;
42 while(len <= maxnum*2) len <<= 1;
43 for(int i = 0; i < len; i++)
44 x[i] = comp(num[i], 0);
45 FFT(x, len, 1);
46 for(int i = 0; i < len; i++)
47 x[i] = x[i]*x[i];
48 FFT(x, len , -1);
49 for(int i = 0; i < len; i++)
50 sum[i] = x[i].r+0.5;
51 for(int i = 0; i < len; i+=2)
52 sum[i] -= num[i>>1];//去掉两次取的木棍相同的
53 for(int i = 0; i < len; i ++)
54 sum[i] >>= 1;//算了2次
55 for(int i = 1; i < len; i++)
56 sum[i] += sum[i-1];
57 ll tot = (ll)n*(n-1)*(n-2)/6, ans = tot;
58 for(int i = 0; i <= maxnum; i++)
59 ans -= num[i]*sum[i];//去掉不能组成三角形的
60 printf("%.7fn", 1.0*ans/tot);
61 }
62 return 0;
63 }
View Code
LA4671 给出A串与B串,只含小写字母a、b。问:A串中有多少本质不同的子串满足 与B串长度相同 且 与B串相对应位置字符不同的数量小于k。
题解:a做1,b做0。将B倒着来一遍和A做卷积,可得有多少位置A串与B串都是a。a做0,b做1再来一遍即可。
hash!37做基1000173169做模会冲突!1e9+7做基1e9+9做模也会冲突!双哈希可以过,37做基2^64做模可过,37做基100000000173169LL做模也可过....
1 #include <bits/stdc++.h>
2 #define ull unsigned long long
3 using namespace std;
4 const int N = 4e5+10;
5 struct comp{
6 double r,i;comp(double _r=0,double _i=0){r=_r;i=_i;}
7 comp operator+(const comp x){return comp(r+x.r,i+x.i);}
8 comp operator-(const comp x){return comp(r-x.r,i-x.i);}
9 comp operator*(const comp x){return comp(r*x.r-i*x.i,r*x.i+i*x.r);}
10 };
11 const double pi=acos(-1.0);
12 void FFT(comp a[],int n,int t){
13 for(int i=1,j=0;i<n-1;i++){
14 for(int s=n;j^=s>>=1,~j&s;);
15 if(i<j)swap(a[i],a[j]);
16 }
17 for(int d=0;(1<<d)<n;d++){
18 int m=1<<d,m2=m<<1;
19 double o=pi/m*t;comp _w(cos(o),sin(o));
20 for(int i=0;i<n;i+=m2){
21 comp w(1,0);
22 for(int j=0;j<m;j++){
23 comp &A=a[i+j+m],&B=a[i+j],t=w*A;
24 A=B-t;B=B+t;w=w*_w;
25 }
26 }
27 }
28 if(t==-1)for(int i=0;i<n;i++)a[i].r/=n;
29 }
30 comp x[N], y[N];
31 char a[100010], b[100010];
32 int ans[N];
33 //hash
34 ull xp[122222] = {1}, H[122222];
35 void initHash(){
36 H[0] = a[0];
37 for(int i = 1; a[i]; i++)
38 H[i] = H[i-1]*31uLL+a[i];
39 }
40 ull askHash(int l, int r){
41 if(l == 0) return H[r];
42 return H[r]-H[l-1]*xp[r-l+1];
43 }
44
45 int main(){
46 for(int i = 1; i < 122222; i++) xp[i] = xp[i-1]*31uLL;
47
48 int ca = 1, K;
49 while(scanf("%d", &K), K != -1){
50 scanf(" %s %s", a, b);
51 int lenA = strlen(a), lenB = strlen(b), len = 1;
52 if(lenA < lenB){
53 printf("Case %d: %dn", ca++, 0);
54 continue ;
55 }
56
57 for(int i = 0; i < lenB-1-i; i++)
58 swap(b[i], b[lenB-1-i]);
59 while(len <= lenA+lenB) len <<= 1;
60
61 for(int i = 0; i < len; i++){
62 x[i] = comp(i < lenA? (a[i] == 'a'): 0, 0);
63 y[i] = comp(i < lenB? (b[i] == 'a'): 0, 0);
64 }
65 FFT(x, len, 1);
66 FFT(y, len, 1);
67 for(int i = 0; i < len; i++)
68 x[i] = x[i]*y[i];
69 FFT(x, len, -1);
70 for(int i = 0; i < len; i++)
71 ans[i] = x[i].r+0.5;
72
73 for(int i = 0; i < len; i++){
74 x[i] = comp(i < lenA? (a[i] == 'b') : 0, 0);
75 y[i] = comp(i < lenB? (b[i] == 'b') : 0, 0);
76 }
77 FFT(x, len, 1);
78 FFT(y, len, 1);
79 for(int i = 0; i < len; i++)
80 x[i] = x[i]*y[i];
81 FFT(x, len, -1);
82 for(int i = 0; i < len; i++)
83 ans[i] += x[i].r+0.5;
84
85 initHash();
86 set<ull> se;
87 for(int i = lenB-1; i < lenA; i++)
88 if(ans[i] >= lenB-K)
89 se.insert(askHash(i-lenB+1, i));
90 printf("Case %d: %dn", ca++, (int)se.size());
91 }
92 return 0;
93 }
View Code
=================================分割线===================================
NTT
NTT与FFT类似,FFT用复数形式会有精度损失,而NTT则是在整数域内取模意义下,无精度损失。
如果 P =r⋅2k+1是个素数,G是模P下的一个原根,那么在mod P意义下,可以处理2k以内规模的数据。
G在模P下的阶为 P-1,即r⋅2k
那么Gr在模P下的阶为2k,这里的 Gr即等价于FFT里的wn.
那么我们用模P下的卷积运算就不会产生精度损失。
P = 998244353 = 119*223+1, 能够处理223= 8e6+ 规模的数据,原根为3.
P = 1004535809=479*221+1, 能够处理221= 2e6+ 规模的数据,原根为3, 且1004535809 加起来不会爆 int.
NTT能解决模数 P =r⋅2k+1的问题,那么如何解决模任意数呢?
先前的 NTT资料 里有提到,
即用多个小素数跑NTT,最后用中国剩余定理求出 n(m-1)2内满足条件的唯一值,当 各个素数积 > n(m-1)2时中国剩余定理后显然可取得唯一值。
=================================分割线===================================
Ck =∑i⊕j=k (Ai*Bj),⊕是某种运算符号。当⊕是+时,即是傅里叶变换;当⊕是^, &, |等某种位运算时,即是FWT快速沃尔什变换。
FFT中,数组长度要大于结果的最高次幂,高位补0;FWT时,数组长度需要是2的整数次幂,不足补0。
原理不是怎么重要...
模板套用即可。
1 //快速沃尔什变换
2 void FWT(int*a,int n){
3 for(int d=1;d<n;d<<=1)for(int m=d<<1,i=0;i<n;i+=m)for(int j=0;j<d;j++){
4 int x=a[i+j],y=a[i+j+d];
5 //xor:a[i+j]=x+y,a[i+j+d]=x−y;
6 //and:a[i+j]=x+y;
7 //or:a[i+j+d]=x+y;
8 }
9 }
10 void UFWT(int*a,int n){
11 for(int d=1;d<n;d<<=1)for(int m=d<<1,i=0;i<n;i+=m)for(int j=0;j<d;j++){
12 int x=a[i+j],y=a[i+j+d];
13 //xor:a[i+j]=(x+y)/2,a[i+j+d]=(x−y)/2;
14 //and:a[i+j]=x−y;
15 //or:a[i+j+d]=y−x;
16 }
17 }
View Code
防溢出可mod 1e9+7大质数,则除以2的时候乘2的逆元。
FWT后,相乘,UFWT回去即可。
原文链接: https://www.cnblogs.com/dirge/p/5887294.html
欢迎关注
微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍
原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/240817
非原创文章文中已经注明原地址,如有侵权,联系删除
关注公众号【高性能架构探索】,第一时间获取最新文章
转载文章受原作者版权保护。转载请注明原作者出处!