多项式算法学习笔记

之前零零散散学了一些多项式的算法，准备总结一下。

先是大佬博客https://blog.csdn.net/WADuan2/article/details/79529900 (原理)

无敌的总结和题目https://www.cnblogs.com/zhouzhendong/p/Fast-Fourier-Transform.html

这篇博客没有什么原创的东西，几乎就是看大佬博客学习总结一下，方便自己回忆，也可以方便新手入门学习。

FFT：原理建议看上面大佬博客，感觉说得比较通俗易懂了。FFT就是利用复数的性质来快速计算多项式乘法的算法，两个多项式相乘暴力的时间复杂度是O(n^2)，用FFT来计算的时间复杂度是O(nlogn)。

模板题LOJ#108 多项式乘法：计算两个多项式乘法输出结果

 1 #include<iostream>
 2 #include<cstring>
 3 #include<cstdio>
 4 #include<cmath>
 5 using namespace std;
 6 const int N=400005;  //开四倍所有数组(A数组两倍*B数组两倍)
 7 const double pi=acos(-1.0);
 8 struct complex {
 9     double x,y;
10     complex(double X=0,double Y=0) { x=X,y=Y; }
11 }A[N],B[N];  //两个多项式 
12 complex operator + (complex a,complex b) { return complex(a.x+b.x,a.y+b.y); }
13 complex operator - (complex a,complex b) { return complex(a.x-b.x,a.y-b.y); }
14 complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); }
15 int S,T,n,m,L,R[N];
16 long long F[N];//存储多项式相乘结果的系数从0-n
17 double f[N],g[N];  //输入两个多项式 
18 
19 void FFT(complex A[N],int opt) {  //opt=1求值，opt=2插值 
20     for (int i=0; i<n; ++i)
21         if (i<R[i]) swap(A[i],A[R[i]]);
22     for (int k=1; k<n; k<<=1) {
23         complex wn=complex(cos(pi/k),opt*sin(pi/k));
24         for (int i=0; i<n; i+=(k<<1)) {
25             complex w=complex(1,0);
26             for (int j=0; j<k; ++j,w=w*wn) {
27                 complex x=A[i+j],y=w*A[i+j+k];
28                 A[i+j]=x+y,A[i+j+k]=x-y;
29             }
30         }
31     }
32 }
33 
34 void calc(int opt) {
35     FFT(A,1); FFT(B,1);  //分别求点 
36     for (int i=0; i<=n; ++i) A[i]=A[i]*B[i];  //点相乘 
37     FFT(A,-1);  //插值 
38     for (int i=0; i<=n; ++i) { //上限为n，在运用时灵活修改。
39         F[i]=(long long)(A[i].x/n+0.5)*opt;
40         //cout<<F[i]<<endl;
41     }
42 }
43 
44 int main() {
45     scanf("%d%d",&S,&T);
46     for(int i=0; i<=S; i++) scanf("%lf",&f[i]);  //第一个多项式 
47     for(int i=0; i<=T; i++) scanf("%lf",&g[i]);  //第二个多项式 
48     
49     m=S+T;  //结果长度lenS+lenT-1,例如：111*11=4位 
50     for (n=1; n<=m; n<<=1) ++L;
51     for (int i=0; i<n; ++i)
52         R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));  ///R[i]代表反转后的位置
53     for (int i=0; i<=n; ++i)
54         A[i]=complex(1.0*f[i],0.0),B[i]=complex(1.0*g[i],0.0);
55     calc(1);
56     for(int i=0; i<=S+T; i++) printf("%d ",F[i]);
57     return 0;
58 }

View Code

NTT：FFT利用复数进行运算在数较大的时候会产生比较大的精度误差，所以我们有了NTT。NTT的原理和FFT差不太多，但是NTT是利用原根来代替复数进行运算（当然也因为原根能有复数的一些性质才能用原根代替复数），原根是整数所以不会产生精度误差且能够取模。

洛谷P5395 第二类斯特林数：模板题第n行第二类斯特林数，用NTT计算式子即可。

 1 #include<bits/stdc++.h>
 2 #define LL long long
 3 using namespace std;
 4 const int N=1e6+10;
 5 const LL P=167772161,yg=3;
 6 LL n,fac[N],inv[N],f[N],g[N],S2[N];
 7 LL bin[N];
 8 
 9 LL power(LL x,LL p) {
10     LL ret=1;
11     for (;p;p>>=1) {
12         if (p&1) ret=(ret*x)%P;
13         x=(x*x)%P;
14     }
15     return ret;
16 }
17 
18 void NTT(LL *a,LL n,LL op) {  //NTT:系数a数组，长度为n，op=1求值op=-1插值 
19     for(LL i=0;i<n;i++) bin[i]=(bin[i>>1]>>1)|((i&1)*(n>>1));
20     for(LL i=0;i<n;i++) if(i<bin[i]) swap(a[i],a[bin[i]]);
21     for(LL i=1;i<n;i<<=1) {
22         LL wn=power(yg,op==1?(P-1)/(2*i):(P-1)-(P-1)/(2*i)),w,t;
23         for(LL j=0;j<n;j+=i<<1) {
24             w=1;
25             for(LL k=0;k<i;k++) {
26                 t=a[i+j+k]*w%P;w=w*wn%P;
27                 a[i+j+k]=(a[j+k]-t+P)%P;a[j+k]=(a[j+k]+t)%P;
28             }
29         }
30     }
31     if(op==-1) {
32         LL Inv=power(n,P-2);
33         for(LL i=0;i<n;i++) a[i]=a[i]*Inv%P;
34     }
35 }
36 
37 int main()
38 {
39     cin>>n;
40     fac[0]=inv[0]=1;
41     for (int i=1;i<=n;i++) fac[i]=fac[i-1]*i%P,inv[i]=power(fac[i],P-2);
42     for (int i=0;i<=n;i++) f[i]=(power(-1,i)+P)%P*inv[i]%P;
43     for (int i=0;i<=n;i++) g[i]=power(i,n)*inv[i]%P;
44     
45     LL len=1;while(len<(n+1)<<1) len<<=1;
46     
47     NTT(f,len,1); NTT(g,len,1);
48     for (int i=0;i<len;i++) S2[i]=(f[i]*g[i])%P;  //求f.g的卷积为S2 
49     NTT(S2,len,-1);
50     
51     for (int i=0;i<=n;i++) {
52         //S2[i]=(S2[i-1]+S2[i])%P;
53         printf("%lld ",S2[i]);
54     }
55     return 0;
56 }

View Code

MTT：当运算的数更加巨大且要求对结果取模的话我们用朴素的NTT也无法满足条件了。所以任意模数NTT（即MTT）诞生了，一般来说有两种常见方式能实现，一种是拆系数FFT，一种是三模数NTT（原理建议看大佬博客），博主学习的是三模数NTT。

洛谷P4245 任意模数NTT：MTT模板题，计算两个多项式乘法取模后的结果。

#include<bits/stdc++.h>
#define LL long long 
const int MAXN = 3 * 1e6 + 10;
using namespace std;
const int P1 = 469762049, P2 = 998244353, P3 = 1004535809, g = 3; 
const LL PP = 1ll * P1 * P2;
int N, M, P, limit = 1, L;
int A[MAXN], B[MAXN], C[MAXN], D[MAXN], Ans[3][MAXN], r[MAXN],ans[MAXN];
LL fastmul(LL a, LL b, LL mod) {
    a %= mod, b %= mod;
    return ((a * b - (LL)((LL)((long double)a / mod * b + 1e-3) * mod)) % mod + mod) % mod;
}
int fastpow(int a, int p, int mod) {
    int base = 1;
    while(p) {
        if(p & 1) base = 1ll * a * base % mod;
        a = 1ll * a * a % mod; p >>= 1;
    }
    return base % mod;
} 
void NTT(int *A, const int n, const int type, const int mod) {
    for(int i = 0; i < n; i++)
        if(i < r[i]) swap(A[i], A[r[i]]);
    for(int mid = 1; mid < n; mid <<= 1) {
        int W = fastpow(type == 1 ? g : fastpow(g, mod - 2, mod) , (mod - 1) / (mid << 1), mod);
        for(int j = 0; j < n; j += (mid << 1)) {
            int w = 1;
            for(int k = 0; k <mid; k++, w = 1ll * w * W % mod) {
                int x = A[j + k], y = 1ll * w * A[j + k + mid] % mod;
                A[j + k] = (x + y) % mod,
                A[j + k + mid] = (x - y + mod) % mod;
            }
        }
    }
    if(type == -1) {
        int inv = fastpow(n, mod - 2, mod);
        for(int i = 0; i < n; i++) 
            A[i] = 1ll * A[i] * inv % mod;
    }
}

void MTT() {
    while(limit <= N + M) limit <<= 1, L++;
    for(int i = 0; i <= limit; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
    
    copy(A, A + N + 1, C); copy(B, B + M + 1, D);
    NTT(C, limit, 1, P1); NTT(D, limit, 1, P1);
    for(int i = 0; i <= limit; i++) Ans[0][i] = 1ll * C[i] * D[i] % P1;
    
    memset(C, 0, sizeof(C)); memset(D, 0, sizeof(D));
    copy(A, A + N + 1, C); copy(B, B + M + 1, D);
    NTT(C, limit, 1, P2); NTT(D, limit, 1, P2);
    for(int i = 0; i <= limit; i++) Ans[1][i] = 1ll * C[i] * D[i] % P2;
    
    memset(C, 0, sizeof(C)); memset(D, 0, sizeof(D));
    copy(A, A + N + 1, C); copy(B, B + M + 1, D);
    NTT(C, limit, 1, P3); NTT(D, limit, 1, P3);
    for(int i = 0; i <= limit; i++) Ans[2][i] = 1ll * C[i] * D[i] % P3;
    
    NTT(Ans[0], limit, -1, P1);
    NTT(Ans[1], limit, -1, P2);
    NTT(Ans[2], limit, -1, P3);
    
    for(int i = 0; i <= N + M; i++) {
        LL A = (fastmul(1ll * Ans[0][i] * P2 % PP, fastpow(P2 % P1, P1 - 2, P1), PP) + 
                fastmul(1ll * Ans[1][i] * P1 % PP, fastpow(P1 % P2, P2 - 2, P2), PP) ) % PP;
        LL K = ((Ans[2][i] - A) % P3 + P3) % P3 * fastpow(PP % P3, P3 - 2, P3) % P3;
        ans[i]=(A % P + ((K % P) * (PP % P)) % P ) % P;    //结果     
    }
}

int main() 
{
    scanf("%d%d%d",&N,&M,&P);  //A项次数，B项次数，模数P 
    for(int i = 0; i <= N; i++) scanf("%d",&A[i]);  //A项 
    for(int i = 0; i <= M; i++) scanf("%d",&B[i]);  //B项 
    
    MTT();
    
    for (int i=0;i<=N+M;i++) printf("%d ",ans[i]);
    return 0;
}

View Code

分治FFT：朴素的FFT/NTT是用来计算卷积的。

https://blog.csdn.net/VictoryCzt/article/details/82939586

了解原理还没有什么用，关键是知道这些算法能做什么。只能通过做题来了解一些做题常见套路。

题目练习：

计蒜客2018宁夏网络赛H

题意：给定游戏的五种手势和十种克制关系。每种手势克制其他两种手势并被另外两种手势克制。给你两个字符串分别表示A,B的手势顺序，A长B短，你可以随意挪动B相对于A的位置，求B最多能赢多少次？

解法：匹配问题是FFT的经典应用。我们可以通过一次多项式乘法做到A串在B串的每个位置开始匹配的结果。且这种匹配不要求连续(不像KMP)，只要求对应位置相等，当然这不是真正的字符串匹配只是通过人为的赋值然后通过卷积计算得到从每个位置开始的答案，一次FFT就达到了多次匹配的目的（当然这种人为的赋值就得相当的有智慧，具体问题具体分析）。

那么回到此题：比较容易想到的办法就是分成5次计算，每次令B串为一个手势，令A串为被B克制的另两个手势，然后FFT算贡献。具体怎么做：假如J克制L,J也克制K。那么令A串位置为J的项值为1，其他为0，同时令B串位置为K/L的值为1其他为0。AB串做一次FFT就能得到的结果就是从每个位置开始匹配的贡献值。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
const int N=4e6+10;  //¿ªËÄ±¶ËùÓÐÊý×é(AÊý×éÁ½±¶*BÊý×éÁ½±¶)
const double pi=acos(-1.0);
struct complex {
    double x,y;
    complex(double X=0,double Y=0) { x=X,y=Y; }
}A[N],B[N];  //Á½¸ö¶àÏîÊ½ 
complex operator + (complex a,complex b) { return complex(a.x+b.x,a.y+b.y); }
complex operator - (complex a,complex b) { return complex(a.x-b.x,a.y-b.y); }
complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); }
int S,T,n,m,L,R[N];
long long F[N];//´æ´¢¶àÏîÊ½Ïà³Ë½á¹ûµÄÏµÊý´Ó0-n
double f[N],g[N];  //ÊäÈëÁ½¸ö¶àÏîÊ½ 

void FFT(complex A[N],int opt) {  //opt=1ÇóÖµ£¬opt=2²åÖµ 
    for (int i=0; i<n; ++i)
        if (i<R[i]) swap(A[i],A[R[i]]);
    for (int k=1; k<n; k<<=1) {
        complex wn=complex(cos(pi/k),opt*sin(pi/k));
        for (int i=0; i<n; i+=(k<<1)) {
            complex w=complex(1,0);
            for (int j=0; j<k; ++j,w=w*wn) {
                complex x=A[i+j],y=w*A[i+j+k];
                A[i+j]=x+y,A[i+j+k]=x-y;
            }
        }
    }
}

void calc(int opt) {
    FFT(A,1); FFT(B,1);  //·Ö±ðÇóµã 
    for (int i=0; i<=n; ++i) A[i]=A[i]*B[i];  //µãÏà³Ë 
    FFT(A,-1);  //²åÖµ 
    for (int i=0; i<=n; ++i) { //ÉÏÏÞÎªn£¬ÔÚÔËÓÃÊ±Áé»îÐÞ¸Ä¡£
        F[i]+=(long long)(A[i].x/n+0.5)*opt;
        //cout<<F[i]<<endl;
    }
}

char s1[N],s2[N];
int main() {
    //scanf("%d%d",&S,&T);
    //for(int i=0; i<=S; i++) scanf("%lf",&f[i]);  //µÚÒ»¸ö¶àÏîÊ½ 
    //for(int i=0; i<=T; i++) scanf("%lf",&g[i]);  //µÚ¶þ¸ö¶àÏîÊ½ 
    scanf("%s%s",s1,s2);
    S=strlen(s1); T=strlen(s2);
    reverse(s2,s2+T);
    
    
    m=S+T-2;  //m+1=½á¹û³¤¶È 
    for (n=1; n<=m; n<<=1) ++L;
    for (int i=0; i<n; ++i) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));  ///R[i]´ú±í·´×ªºóµÄÎ»ÖÃ
    
    //S>P / S>L
    memset(f,0,sizeof(f)); memset(g,0,sizeof(g));
    for (int i=0;i<S;i++) f[i]=(s1[i]=='P' || s1[i]=='L')?1:0;
    for (int i=0;i<T;i++) g[i]=(s2[i]=='S')?1:0;
    for (int i=0; i<=n; ++i) A[i]=complex(1.0*f[i],0.0),B[i]=complex(1.0*g[i],0.0);
    calc(1);
   
       //P>R&&P>K
       memset(f,0,sizeof(f)); memset(g,0,sizeof(g));
    for (int i=0;i<S;i++) f[i]=(s1[i]=='R' || s1[i]=='K')?1:0;
    for (int i=0;i<T;i++) g[i]=(s2[i]=='P')?1:0;
    for (int i=0; i<=n; ++i) A[i]=complex(1.0*f[i],0.0),B[i]=complex(1.0*g[i],0.0);
    calc(1);
    
       //R>S&&R>L
       memset(f,0,sizeof(f)); memset(g,0,sizeof(g));
    for (int i=0;i<S;i++) f[i]=(s1[i]=='S' || s1[i]=='L')?1:0;
    for (int i=0;i<T;i++) g[i]=(s2[i]=='R')?1:0;
    for (int i=0; i<=n; ++i) A[i]=complex(1.0*f[i],0.0),B[i]=complex(1.0*g[i],0.0);
    calc(1);
       
       //L>P&&L>K
       memset(f,0,sizeof(f)); memset(g,0,sizeof(g));
    for (int i=0;i<S;i++) f[i]=(s1[i]=='P' || s1[i]=='K')?1:0;
    for (int i=0;i<T;i++) g[i]=(s2[i]=='L')?1:0;
    for (int i=0; i<=n; ++i) A[i]=complex(1.0*f[i],0.0),B[i]=complex(1.0*g[i],0.0);
    calc(1);
       
       //K>R&&K>S
       memset(f,0,sizeof(f)); memset(g,0,sizeof(g));
    for (int i=0;i<S;i++) f[i]=(s1[i]=='R' || s1[i]=='S')?1:0;
    for (int i=0;i<T;i++) g[i]=(s2[i]=='K')?1:0;
    for (int i=0; i<=n; ++i) A[i]=complex(1.0*f[i],0.0),B[i]=complex(1.0*g[i],0.0);
    calc(1);
    
    long long ans=0;
    for (int i=0;i<=n;i++) ans=max(ans,F[i]);
    cout<<ans<<endl;
    return 0;
}

View Code

BZOJ4259: 残缺的字符串

题意：给你两个字符串原串和匹配串，串里只包含小写字母和 $*$

$*$

 1 #include<iostream>
 2 #include<cstring>
 3 #include<cstdio>
 4 #include<cmath>
 5 #include<algorithm>
 6 using namespace std;
 7 const int N=4*(3e5+10);  //开四倍所有数组(A数组两倍*B数组两倍)
 8 const double pi=acos(-1.0);
 9 struct complex {
10     double x,y;
11     complex(double X=0,double Y=0) { x=X,y=Y; }
12 }A[N],B[N];  //两个多项式 
13 complex operator + (complex a,complex b) { return complex(a.x+b.x,a.y+b.y); }
14 complex operator - (complex a,complex b) { return complex(a.x-b.x,a.y-b.y); }
15 complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); }
16 int S,T,n,m,L,R[N];
17 long long ans[N],F[N];//存储多项式相乘结果的系数从0-n
18 double f[N],g[N];  //输入两个多项式 
19 
20 void FFT(complex A[N],int opt) {  //opt=1求值，opt=2插值 
21     for (int i=0; i<n; ++i)
22         if (i<R[i]) swap(A[i],A[R[i]]);
23     for (int k=1; k<n; k<<=1) {
24         complex wn=complex(cos(pi/k),opt*sin(pi/k));
25         for (int i=0; i<n; i+=(k<<1)) {
26             complex w=complex(1,0);
27             for (int j=0; j<k; ++j,w=w*wn) {
28                 complex x=A[i+j],y=w*A[i+j+k];
29                 A[i+j]=x+y,A[i+j+k]=x-y;
30             }
31         }
32     }
33 }
34 
35 void calc(int opt) {
36     FFT(A,1); FFT(B,1);  //分别求点 
37     for (int i=0; i<=n; ++i) A[i]=A[i]*B[i];  //点相乘 
38     FFT(A,-1);  //插值 
39     for (int i=0; i<=n; ++i) { //上限为n，在运用时灵活修改。
40         F[i]+=(long long)(A[i].x/n+0.5)*opt;
41         //cout<<F[i]<<endl;
42     }
43 }
44 
45 char s[N],t[N];
46 int main() {
47     scanf("%d%d",&T,&S);
48     scanf("%s%s",t,s);
49     //for(int i=0; i<=S; i++) scanf("%lf",&f[i]);  //第一个多项式 
50     //for(int i=0; i<=T; i++) scanf("%lf",&g[i]);  //第二个多项式 
51     S--; T--;
52     for (int i=0;i<=S;i++) if (s[i]=='*') s[i]='a'-1;
53     for (int j=0;j<=T;j++) if (t[j]=='*') t[j]='a'-1;
54     reverse(t,t+T+1);
55     
56     m=S+T;  //结果长度 
57     for (n=1; n<=m; n<<=1) ++L;
58     for (int i=0; i<n; ++i) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));  ///R[i]代表反转后的位置
59     
60     memset(f,0,sizeof(f)); memset(g,0,sizeof(g));
61     for (int i=0;i<=S;i++) f[i]=(s[i]-'a'+1)*(s[i]-'a'+1)*(s[i]-'a'+1);
62     for (int i=0;i<=T;i++) g[i]=t[i]-'a'+1;
63     for (int i=0; i<=n; ++i) A[i]=complex(1.0*f[i],0.0),B[i]=complex(1.0*g[i],0.0);
64     calc(1);
65     
66     memset(f,0,sizeof(f)); memset(g,0,sizeof(g));
67     for (int i=0;i<=S;i++) f[i]=s[i]-'a'+1;
68     for (int i=0;i<=T;i++) g[i]=(t[i]-'a'+1)*(t[i]-'a'+1)*(t[i]-'a'+1);
69     for (int i=0; i<=n; ++i) A[i]=complex(1.0*f[i],0.0),B[i]=complex(1.0*g[i],0.0);
70     calc(1);
71     
72     memset(f,0,sizeof(f)); memset(g,0,sizeof(g));
73     for (int i=0;i<=S;i++) f[i]=(s[i]-'a'+1)*(s[i]-'a'+1);
74     for (int i=0;i<=T;i++) g[i]=(t[i]-'a'+1)*(t[i]-'a'+1);
75     for (int i=0; i<=n; ++i) A[i]=complex(1.0*f[i],0.0),B[i]=complex(1.0*g[i],0.0);
76     calc(-2);
77     
78    
79        long long ans=0;
80        for (int i=T;i<=m;i++) if (!F[i]) ans++;
81        cout<<ans<<endl;
82        for (int i=T;i<=S;i++)
83            if (!F[i]) printf("%d ",i-T+1);
84     return 0;
85 }

View Code

Vjudge2016香港网络赛

题意：给你一个大小为2e5的数组，问有多少对a[i]+a[j]==a[k],i,j,k两两不同。

解法：这道题解法不是特别难想但是细节比较多容易想漏情况。解法是我们令a[i]为多项式的下标，a[i]的出现次数为多项式的系数，然后多项式自己和自己相乘。这样得到的结果是：F[i]的得到的指数即为所有可能加法的结果，而对应的系数即为加法结果出现次数。即F[i]系数就是从原数组a中选两个数(这两个数可以是同一个数)相加结果得到i的方案数为F[i]。然后我们遍历原数组看看原数组a[i]在F有几次出现就是有几个等式啦。

但是要注意数有负数，自己相加和，零相加等等情况，具体细节要考虑清楚。

 1 #include<algorithm>
 2 #include<iostream>
 3 #include<cstring>
 4 #include<cstdio>
 5 #include<cmath>
 6 using namespace std;
 7 #define maxn 800005
 8 typedef long long ll;
 9 const int LOW=50000;
10 const double pi=acos(-1.0);
11 struct complex
12 {
13     double x,y;
14     complex(double X=0,double Y=0)
15     {
16         x=X,y=Y;
17     }
18 }a[maxn],b[maxn];
19 complex operator + (complex a,complex b) {return complex(a.x+b.x,a.y+b.y);}
20 complex operator - (complex a,complex b) {return complex(a.x-b.x,a.y-b.y);}
21 complex operator * (complex a,complex b) {return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
22 int S,T,n,m,L,R[maxn],zero;
23 long long F[maxn],num[maxn],cnt[maxn];
24 void FFT(complex a[maxn],int opt)
25 {
26     for (int i=0;i<n;++i)
27         if (i<R[i]) swap(a[i],a[R[i]]);
28     for (int k=1;k<n;k<<=1)
29     {
30         complex wn=complex(cos(pi/k),opt*sin(pi/k));
31         for (int i=0;i<n;i+=(k<<1))
32         {
33             complex w=complex(1,0);
34             for (int j=0;j<k;++j,w=w*wn)
35             {
36                 complex x=a[i+j],y=w*a[i+j+k];
37                 a[i+j]=x+y,a[i+j+k]=x-y;
38             }
39         }
40     }
41 }
42 void calc(int opt)
43 {
44     FFT(a,1);FFT(b,1);
45     for (int i=0;i<=n;++i) a[i]=a[i]*b[i];
46     FFT(a,-1);
47     for (int i=0;i<n;++i)//对于统计两数组加法的都要完全统计0-n
48     {
49          F[i]+=(long long)(a[i].x/n+0.5)*opt;
50     }
51 }
52 int main()
53 {
54     scanf("%d",&S);
55     T=S;
56     for(int i=0;i<S;i++)
57     {
58         scanf("%lld",&num[i]);
59         if(num[i]==0) zero++;
60         cnt[num[i]+LOW]++;//离散化
61     }
62     for (n=1;n<=200000;n<<=1) ++L;
63     for (int i=0;i<n;++i)
64         R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
65     for (int i=0;i<=n;++i) a[i]=complex(1.0*cnt[i],0.0),b[i]=complex(1.0*cnt[i],0.0);
66     calc(1);//传递的参数为多项式乘法的常数可以为正负
67 
68     for(int i=0;i<T;i++)//去除本身与本身相加的结果
69     {
70         F[(num[i]+LOW)*2]--;
71     }
72     ll pp = 0;
73     for(int i=0;i<T;i++)
74     {
75         pp+=F[num[i]+LOW*2];
76         if(num[i]==0)
77             pp-=(zero-1)*2;//元素为0的情况
78         else
79             pp-=zero*2;//元素非0的情况
80     }
81     printf("%lld
",pp);
82     return 0;
83 }

View Code

HDU-4609

题意：给出n个数，问随机选三个数，这三个数能组成三角形的概率。

解法：感觉这道题还是蛮有意思的，解法是看kuangbin的。把数字当成下标(在卷积中数组下标即指数)，出现次数当成系数，自己跟自己做卷积。这样就得到选两个数组成的所有长度及其方案数了，当然还需要处理一下：除去一个数自己和自己相加的情况，方案数要除以2才是组合数。

得到这个方案数数组我们就可以开始计算答案了，我们考虑计算所有的不能组成三角形的方案数（具体就是枚举每一个a[i],然后计算以a[i]为最长边的不合法方案，即选两个数之和<=a[i]，就是sigma(F[a[i]])啦），然后总方案数C(n,3)-sum就是合法方案，那么这题就解决了。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
typedef long long LL;
const int N=4e5+10;  //开四倍所有数组(A数组两倍*B数组两倍)
const double pi=acos(-1.0);
struct complex {
    double x,y;
    complex(double X=0,double Y=0) { x=X,y=Y; }
}A[N],B[N];  //两个多项式 
complex operator + (complex a,complex b) { return complex(a.x+b.x,a.y+b.y); }
complex operator - (complex a,complex b) { return complex(a.x-b.x,a.y-b.y); }
complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); }
int S,T,n,m,L,nl,a[N],R[N];
long long F[N];//存储多项式相乘结果的系数从0-n
double f[N],g[N];  //输入两个多项式 

void FFT(complex A[N],int opt) {  //opt=1求值，opt=2插值 
    for (int i=0; i<n; ++i)
        if (i<R[i]) swap(A[i],A[R[i]]);
    for (int k=1; k<n; k<<=1) {
        complex wn=complex(cos(pi/k),opt*sin(pi/k));
        for (int i=0; i<n; i+=(k<<1)) {
            complex w=complex(1,0);
            for (int j=0; j<k; ++j,w=w*wn) {
                complex x=A[i+j],y=w*A[i+j+k];
                A[i+j]=x+y,A[i+j+k]=x-y;
            }
        }
    }
}

void calc(int opt) {
    FFT(A,1); FFT(B,1);  //分别求点 
    for (int i=0; i<=n; ++i) A[i]=A[i]*B[i];  //点相乘 
    FFT(A,-1);  //插值 
    for (int i=0; i<=n; ++i) { //上限为n，在运用时灵活修改。
        F[i]=(long long)(A[i].x/n+0.5)*opt;
        //cout<<F[i]<<endl;
    }
}

int main() {
    int cas; cin>>cas;
    while (cas--) {
        scanf("%d",&nl);
        int Max=0;
        for (int i=1;i<=nl;i++) scanf("%d",&a[i]),Max=max(Max,a[i]);
        S=T=Max; m=S+T;  //结果长度lenS+lenT-1,例如：111*11=4位 
        L=0; for (n=1; n<=m; n<<=1) ++L;
        for (int i=0; i<n; ++i)
            R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));  ///R[i]代表反转后的位置
          
          for (int i=0;i<=n;i++) f[i]=g[i]=0;
        for (int i=1;i<=nl;i++) f[a[i]]++,g[a[i]]++;
        
        for (int i=0; i<=n; ++i)
               A[i]=complex(1.0*f[i],0.0),B[i]=complex(1.0*g[i],0.0);
        calc(1);
        
        for (int i=1;i<=nl;i++) F[a[i]+a[i]]--;  //自己相加 
        for (int i=0;i<=n;i++) F[i]/=2;  //去掉顺序变组合 
        for (int i=0;i<=n;i++) F[i]+=F[i-1];  //前缀和 
        
        LL cant=0;
        for (int i=1;i<=nl;i++) cant+=F[a[i]];  //累加不能组成三角形情况 
        LL sum=(LL)nl*(nl-1)*(LL)(nl-2)/6;  //总方案数 
        printf("%.7lf
",(double)(sum-cant)/sum);
    }
    return 0;
}

View Code

洛谷P3338

题意就是有Fj=∑(i<j) qi*qj/(i−j)^2− ∑(i>j) qi*qj/(i−j)^2，求每一项Fj/qj 。
这题没什么好说的，推式子就完事了。要注意看卷积的式子张什么样子，一定要理解什么到底是卷积。

注意ai的这个i=j+(i-j) , 重点是是a的下标等于bc下标之和，像a3=a0*b3+a1*b2+a2*b1+a3*b0就是卷积，至于sigma的范围是否是0-i没太大所谓，因为像可以去掉a0*b3和a3*b0这一项sigma就变成1-i-1的范围。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
const int N=400005;  //开四倍所有数组(A数组两倍*B数组两倍)
const double pi=acos(-1.0);
struct complex {
    long double x,y;
    complex(double X=0,double Y=0) { x=X,y=Y; }
}A[N],B[N];  //两个多项式 
complex operator + (complex a,complex b) { return complex(a.x+b.x,a.y+b.y); }
complex operator - (complex a,complex b) { return complex(a.x-b.x,a.y-b.y); }
complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); }
int S,T,n,m,L,R[N];
long double f[N],g[N],tmp[N],F[N],G[N];  //输入两个多项式 

void FFT(complex A[N],int opt) {  //opt=1求值，opt=2插值 
    for (int i=0; i<n; ++i)
        if (i<R[i]) swap(A[i],A[R[i]]);
    for (int k=1; k<n; k<<=1) {
        complex wn=complex(cos(pi/k),opt*sin(pi/k));
        for (int i=0; i<n; i+=(k<<1)) {
            complex w=complex(1,0);
            for (int j=0; j<k; ++j,w=w*wn) {
                complex x=A[i+j],y=w*A[i+j+k];
                A[i+j]=x+y,A[i+j+k]=x-y;
            }
        }
    }
}

void calc(int opt) {
    FFT(A,1); FFT(B,1);  //分别求点 
    for (int i=0; i<=n; ++i) A[i]=A[i]*B[i];  //点相乘 
    FFT(A,-1);  //插值 
    if (opt==1)
    for (int i=0; i<=n; ++i) F[i]+=(A[i].x/n)*opt;
       else 
       for (int i=0; i<=n; ++i) G[i]+=(A[i].x/n)*opt;
}

int main() {
    scanf("%d",&n);
    for (int i=1;i<=n;i++) scanf("%Lf",&f[i]),g[i]=1.0/(long double)i/(long double)i,tmp[i]=f[i];
    S=n; T=n; m=S+T;
    L=0; for (n=1; n<=m; n<<=1) ++L;
    for (int i=0; i<n; ++i)
    R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    

    for (int i=0;i<=n;i++) A[i]=complex(1.0*f[i],0.0),B[i]=complex(g[i],0.0);
    calc(1);
    for (int i=1;i<=S;i++) f[i]=tmp[S-i+1];
    for (int i=0;i<=n;i++) A[i]=complex(1.0*f[i],0.0),B[i]=complex(g[i],0.0);
    calc(-1);
    
    for (int i=1;i<=S;i++) printf("%.5Lf
",F[i]+G[S-i+1]);
    return 0;
}

View Code

洛谷P4199

题意：给一个只有ab字母的字符串，求不是全连续的回文字串数量，例如aba的a a是满足条件，但是aba不满足。

解法：解法非常巧妙地运用卷积的特性。首先容易发现答案=所有可不连续回文串-回文字串。显然后面的可以用manache解决，难点在于怎么计算一个串的可不连续回文串数量呢？

我们观察卷积，例如f5=a0*b5+a1*b4+a2*b3+a3*b2+a4*b1+a5*b0，发现其实卷积是首位交叉相乘，这一点很对称很像回文。那么我们考虑从这个性质入手，构建一个数组A[i]=(st[i]=='a')。然后我们对这个A做自己和自己的卷积会是什么？例如f5=a1*b5+a1*a4+a2*a3+a3*a2+a4*a1+a5*a0=对称位都是a的对数的两倍！！！得到这个之后我们就可以计算可不连续回文字串数，如f5是回文对数，那么我们我们可以从f5/2对中选一些对出来构成可不连续回文串答案就是(1<<(f5/2))-1 (不能取空)。那么对字母b同理。

根据这个思路写这题就可解了，但是注意一些细节，例如当回文长度是奇数的时候对数其实可以+1。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
const int N=400005;  //开四倍所有数组(A数组两倍*B数组两倍)
const int MOD=1e9+7;
const double pi=acos(-1.0);
struct complex {
    double x,y;
    complex(double X=0,double Y=0) { x=X,y=Y; }
}A[N],B[N];  //两个多项式 
complex operator + (complex a,complex b) { return complex(a.x+b.x,a.y+b.y); }
complex operator - (complex a,complex b) { return complex(a.x-b.x,a.y-b.y); }
complex operator * (complex a,complex b) { return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); }
int S,T,n,m,L,R[N];
long long F[N];//存储多项式相乘结果的系数从0-n
double f[N],g[N];  //输入两个多项式 
char st[N],news[N];

void FFT(complex A[N],int opt) {  //opt=1求值，opt=2插值 
    for (int i=0; i<n; ++i)
        if (i<R[i]) swap(A[i],A[R[i]]);
    for (int k=1; k<n; k<<=1) {
        complex wn=complex(cos(pi/k),opt*sin(pi/k));
        for (int i=0; i<n; i+=(k<<1)) {
            complex w=complex(1,0);
            for (int j=0; j<k; ++j,w=w*wn) {
                complex x=A[i+j],y=w*A[i+j+k];
                A[i+j]=x+y,A[i+j+k]=x-y;
            }
        }
    }
}

void calc(int opt) {
    FFT(A,1); FFT(B,1);  //分别求点 
    for (int i=0; i<=n; ++i) A[i]=A[i]*B[i];  //点相乘 
    FFT(A,-1);  //插值 
    for (int i=0; i<=n; ++i) { //上限为n，在运用时灵活修改。
        F[i]+=(long long)(A[i].x/n+0.5)/2*opt;
        //cout<<F[i]<<endl;
    }
}

int power(int x,int p) {
    int ret=1;
    for (;p;p>>=1) {
        if (p&1) ret=(long long)ret*x%MOD;
        x=(long long)x*x%MOD;
    }
    return ret;
}

int p[N];
int Init() {
    int len=strlen(st);
    news[0]='$';
    news[1]='#';
    int j=1;
    for (int i=0;i<len;i++) {
        news[++j]=st[i];
        news[++j]='#';
    }
    news[++j]='';  //最前最后不同字符 
    return j;
}
long long Manacher() {
    int len=Init(),id,mx=0;
    long long sum=0;
    int res=1;
    for (int i=1;i<len;i++) {  //news[0]='$'; news[len]='';
        if (i<mx) p[i]=min(p[2*id-i],mx-i); else p[i]=1;
        while (news[i-p[i]]==news[i+p[i]]) p[i]++;
        if (i+p[i]>mx) {
            id=i; mx=i+p[i];
        }
        sum=(sum+(p[i]>>1))%MOD;  //这样可以计算以字符i为结尾的回文串个数 
    }
    return sum;  
}

int main() {
    scanf("%s",st);
    S=T=strlen(st)-1;
    m=S+T;  //结果长度lenS+lenT-1,例如：111*11=4位 
    L=0; for (n=1; n<=m; n<<=1) ++L;
    for (int i=0; i<n; ++i)
        R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));  ///R[i]代表反转后的位置
    for (int i=0;i<=n;i++)
        if (st[i]=='a') A[i]=B[i]=complex(1.0,0.0); else A[i]=B[i]=complex(0.0,0.0);
    calc(1);
    for (int i=0;i<=n;i++)
        if (st[i]=='b') A[i]=B[i]=complex(1.0,0.0); else A[i]=B[i]=complex(0.0,0.0);
    calc(1);    
    
    long long ans=0;
    for (int i=0;i<=m;i++) {
        F[i]+=(i%2==0);
        if (F[i]>=0) ans=(ans+power(2,F[i]))%MOD,ans=(ans-1+MOD)%MOD;
    }
    ans-=Manacher();
    ans=(ans%MOD+MOD)%MOD;
    cout<<ans<<endl;
    return 0;
}

View Code

HDU-5829

题意：给定 $n$

$n$

$n$

$n$

解法：

相关阅读:
Vmware安装CentOs7+gitlab(一)
设计模式(一) 动态代理初尝试
 APP 技术支持
 App隐私条款
 Mac OS X中开启或关闭显示隐藏文件命令
 软件工程学习要点
 如何实现企业信息化
 记录一些flutter学习网址
 基于深度学习的语义分割
 对自动变速器的控制器建模
原文地址：https://www.cnblogs.com/clno1/p/11275049.html