zoukankan      html  css  js  c++  java
  • 卷积FFT、NTT、FWT

    先简短几句话说说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("%.7f
    ", 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: %d
    ", 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: %d
    ", ca++, (int)se.size());
    91     }
    92     return 0;
    93 }
    View Code

    =================================分割线=================================== 

    NTT

    NTT资料

    NTT资料2

    NTT与FFT类似,FFT用复数形式会有精度损失,而NTT则是在整数域内取模意义下,无精度损失。

    如果 P = r2+1 是个素数,G是模P下的一个原根,那么在mod P 意义下,可以处理 2k 以内规模的数据 。

    G在模P下的阶为 P-1,即 r2

    那么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 = r2+1 的问题,那么如何解决模任意数呢?

    先前的 NTT资料 里有提到,

    即用多个小素数跑NTT,最后用中国剩余定理求出 n(m-1)内满足条件的唯一值,当 各个素数积 > n(m-1)2 时中国剩余定理后显然可取得唯一值。

    =================================分割线===================================

    FWT资料

    FWT资料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回去即可。

  • 相关阅读:
    十四、数据库公共字段处理
    十、前端tag、自定义tag、filter和simple_tag
    kafka生产、消费py脚本
    django模板filter及自定义filter
    django基础,前后端分离数据传参
    django基础——使用django form校验数据
    django自带的后台管理框架django-admin
    django基础——前后端分离,页面分页
    django基础——数据库的增删改查
    django基础——models数据库操作
  • 原文地址:https://www.cnblogs.com/dirge/p/5887294.html
Copyright © 2011-2022 走看看