FFT（多项式）专题水分小记（FFT注意事项）

最新推荐文章于 2023-07-17 18:47:23 发布

jokerwyt

最新推荐文章于 2023-07-17 18:47:23 发布

阅读量440

点赞数

分类专栏：题解数学

本文链接：https://blog.csdn.net/jokerwyt/article/details/80199336

版权

题解同时被 2 个专栏收录

186 篇文章 0 订阅

订阅专栏

数学

7 篇文章 0 订阅

订阅专栏

传送

fft的更深理解。
我们都知道，fft是做一个循环卷积。
而多项式快速幂是nlog n^2的（一开始次数界就是n）
考虑一个方法：假设次方后的次数界是m，我们直接做次数界为m的dft，然后对每个点取幂。
当m很大的时候，显然这是很愚蠢的。
但是，假若本身就是求循环卷积，这个方法就可以应用上了
求dft后直接对每个点取幂，再idft就是C^k了。
当然这题不需要再idft了。

CF960G Bandit Blues

即为求 第一类斯特林数 S(n,m)（n个互异元素放m个环中方案）
显然地，S(n)的生成函数为x(x+1)(x+2)..(x+n-1)
分治fft求上式，时间复杂度n log n^2

序列统计

keyword: 乘法+离散对数=加法

https://blog.csdn.net/jokerwyt/article/details/80456871

城市规划

keyword: 带系数卷积拆分，多项式求逆的种种细节

推一下式子化出一个 666
将C拆成阶乘形式，发现是关于k或n-k的。这样可以考虑将他们合并。

需要用一个多项式求逆。推导方法：
(错漏较多)https://blog.csdn.net/xianhaoming/article/details/78761623
(标程太繁琐)http://blog.miskcoo.com/2015/05/polynomial-inverse#mjx-eqn-eqn1
一起食用更佳。
思想是类似快速幂的倍增，已知A(x)B(x)=1(mod x^(n/2))，求C(x)使得A(x)C(x)=1(mod x^n)。
这个n/2是可以下取整或上取整的。
x^n开到需要的次数界（含）以上就行。通常就取相同的，这样最后是1.
注意时间是n log n而不是n log2 n的。
因为 $T(n)=O(n log n)+T(n/2)$ ，将log固定为logn，等比数列求下和。
这里大概有个2的常数？

注意到B的次数是m/2，tmpA的次数是m。
B与tmpA做运算时次数界应该开到m*2.
因为是mo x^m下的逆元，所以m之后的系数可以全部清掉。
B必须全部清掉，不然B的次数就变成m+m/2，下一次做fft的时候可能会循环卷积爆到前面就GG。
tmpA不清也行，因为有memcpy赋值了。

void inverse(ll *A) {
    B[0]=ksm(A[0],mo-2);
    for (ll m=2; m<=M; m*=2) {
        ll csj=m*2;
        init(csj);
        memcpy(tmpA,A,m*sizeof A[0]);
        dft(B,csj,1);
        dft(tmpA,csj,1);
        for (ll i=0; i<csj; i++) B[i]=(2-B[i]*tmpA[i]%mo)*B[i]%mo;
        dft(B,csj,-1);
        for (ll i=m; i<csj; i++) tmpA[i]=B[i]=0;
    }
    memcpy(A,B,sizeof B);
}

分零食

keyword: 多次fft记得清空数组

构造多项式 $F=0+f(1)x+f(2)x^2...$
一个比较显然的事实， $F^n$ 的m次项系数就是n个人给m个，每个人都非空的贡献。
有点生成函数的意思
那么答案就是求 $p(n)=\sum_{i=1}^{n} F^i[m]$
显然 $F^i$ 是一个可以类似快速幂来做的东西
由 $F^{a+b}=F^a * F^b$ ，多项式乘法满足分配律可知

p(n)也可以类似快速幂的做，即
$p(n)=p(n/2)+F^{n/2}\sum F^{i}$
$~~~~~~~~=p(n/2)+F^{n/2}p(n/2)+F^n*(nmod2==1)$
其实是个显然的东西（等比数列分治求和），套个多项式乘法就不认识了。
注意的地方：这题打dft可能会有精度问题，需要ntt.
可以发现只要每轮求完都%p，保证整个过程不会爆ntt模数。
一个事实是：预处理ntt主根并无法带来多大优化。

多次fft记得将0..M-1先都清空。

梦醒

keyword:二维卷积压成一维，循环累加trick，负数ntt-trick.

枚举每一个右下角
将式子拆开，发现预处理不出来的只有 $\sum \sum Cij * Bij$
考虑对于A中一个位置(x,y)与B中一个位置(i,j)重合时，对右下角累加贡献

该累加到的位置是：(x+(h-i),y+(w-j))

显然是一个二维卷积
将B[i][j]放到Bp[h-i][w-j]的位置。就相当于A与B做二维卷积。
不知道二维卷积表述准不准确，我脑补的二维卷积就是这样的..
$Fij=\sum Axy * B(i-x)(j-y)$

将坐标表示为iR+j的形式，压成一维，做FFT即可。

有个小trick：虽然有一些点与B中一些点重合时，右下角不在矩阵内，这样就会导致循环累加。解决方法可以是开大次数界容纳这些多余数据，但是注意到这题可以被循环累加到的位置实质上是不存在的。因此不需要开大次数界。 (infleaking大佬的方法)
（开大次数界带来的是巨大且直观的常数变化，时间直接*2*3*4…）

还有就是虽然有负数，但是绝对值不会超过mo/2，因此mo意义下可以判出哪些是负数，可以用NTT.

#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
const int N=700,mo=998244353,MXM=4*(N+N+1)*N,g=3;
ll n,m,h,w,K;
ll a[N][N],b[N][N],bp[N][N];
ll sum2[N][N],sum[N][N],bs,bs2,tot;
ll yA[MXM],yB[MXM],M,X,Y;
ll QA,H[MXM];

pair<ll,pair<ll,ll> > px[N*N];

ll ksm(ll x,ll y) {
    if (y==0) return 1; if (y==1) return x%mo;
    ll t=ksm(x,y>>1); return t*t%mo*ksm(x,y&1)%mo;
}

void dft(ll *a,ll sig) {
    for (int i=0; i<M; i++) if (H[i]<i) swap(a[i],a[H[i]]);
    for (int m=2; m<=M; m*=2) {
        ll hf=m/2;
        ll zg=ksm(g,(mo-1)/m),w=1;
        if (sig==-1) zg=ksm(zg,mo-2);
        for (int i=0; i<hf; i++,w=w*zg%mo)
        for (int j=i; j<M; j+=m) {
            ll u=a[j],v=a[j+hf]*w%mo;
            a[j]=(u+v)%mo,a[j+hf]=((u-v)%mo+mo)%mo;
        }
    }
}

void fft(ll *u,ll *v) {
//  M=4;
//  for (int i=0; i<=3; i++) u[i]=i;
    for (M=1; M<=QA*(n+1); M<<=1);
    for (int i=0; i<M; i++) H[i]=(H[i>>1]>>1)+(i&1)*(M>>1);
    dft(u,1);
    dft(v,1);
    for (int i=0; i<M; i++) u[i]=u[i]*v[i]%mo;
    dft(u,-1);
    ll cs=ksm(M,mo-2);
    for (int i=0; i<M; i++) u[i]=(u[i]*cs%mo+mo)%mo;
}
ll query(ll sum[N][N],ll lx,ll ly,ll rx,ll ry) {
    return sum[rx][ry] - sum[lx-1][ry] - sum[rx][ly-1] + sum[lx-1][ly-1];
}

int main() {
    freopen("3343.in","r",stdin);
    freopen("3343.out","w",stdout);
    cin>>n>>m; 
    for (int i=1; i<=n; i++) for (int j=1; j<=m; j++) {
        scanf("%lld",&a[i][j]);
        sum[i][j]=sum[i-1][j]+sum[i][j-1]-sum[i-1][j-1]
            +a[i][j];
        sum2[i][j]=sum2[i-1][j]+sum2[i][j-1]-sum2[i-1][j-1]
            +a[i][j]*a[i][j];
    }

    cin>>h>>w; QA=m+1;
    for (int i=1; i<=n; i++) for (int j=1; j<=m; j++) {
        yA[i*QA + j] = a[i][j];
    }

    for (int i=1; i<=h; i++)for (int j=1; j<=w; j++) {
        scanf("%lld",&b[i][j]);
        bp[h-i][w-j]=b[i][j];
        bs+=b[i][j],bs2+=b[i][j] * b[i][j];
    }
    cin>>X>>Y>>K;
    for (int i=0; i<h; i++) for (int j=0; j<w; j++) {
        yB[i*QA + j]=bp[i][j];
    }
    fft(yA,yB);
    for (int i=0; i<M; i++) yA[i]=yA[i]>mo/2 ? yA[i]-mo : yA[i];

    for (int i=1; i<=n; i++) for (int j=1; j<=m; j++) {
        if ((i+h-1)>n || (j+w-1)>m) break;
        ll xy=a[i+X-1][j+Y-1];
        ll ans=xy*xy*h*w+query(sum2,i,j,i+h-1,j+w-1)+bs2;
        ans-=2*xy*query(sum,i,j,i+h-1,j+w-1);
        ans+=2*xy*bs;
        ans-=2*yA[(i+h-1)*QA+j+w-1];
        px[++tot]=make_pair(ans,make_pair(i,j));
    }
    sort(px+1,px+tot+1);
    for (int i=1; i<=K; i++) {
        printf("%lld %lld %lld\n",px[i].second.first,px[i].second.second,px[i].first);
    }
}