ZJT's Blog

大変に気分がいい~

有限交换群的分解算法

2018-9-13 13:142018-9-13 23:50
群论抽象代数

众所周知(?),有限交换群(即有限阿贝尔群)可被分解为若干个质数幂阶循环群的直积。然而沙雕zjt之前一直不知道怎么求出这个分解,最近看了LCA的论文和几篇paper,稍微学习了一个,这里就总结一下吧。

前置知识

有限阿贝尔群基本定理:所有有限阿贝尔群都可以被表示成若干质数幂阶循环子群的直积。(证明naidesu)

Sylow p-子群:现有任意有限群$G$及质数$p$,令$p^k$为能整除$|G|$的最高次幂,则$G$存在阶为$p^k$的子群,称为$G$的Sylow p-子群。在有限交换群中,$G$可以被分解为各个Sylow p-子群的直积。

记$\langle x\rangle$为$x$生成的循环子群,$|\langle x\rangle|$称为$x$的阶,记作$|x|$。以下内容都在某个有限交换群$G$中进行。记$n=|G|$。

线性无关:称$\{b_1,b_2,\cdots,b_k\}$线性无关,当且仅当对于任意$i$,有$\langle b_i\rangle\cap\prod_{j\ne i}\langle b_j\rangle=\{e\}$。若不满足则称为线性相关。元素$x$与集合$A$线性无关的定义类似,即$\{x\}\cup A$线性无关。

基(basis):称$\{b_1,b_2,\cdots,b_k\}$为$G$的一组基,当且仅当$\prod \langle b_i\rangle=G$,且$b_1,b_2,\cdots,b_k$线性无关。

我们要实现的事情是求出$G$的直积分解,并给出对应的基。由于在非交互题中给出整个群就已经是$O(n^2)$的了,这里先介绍一个$O(n^2)$的算法。事实上还存在着复杂度更低甚至线性的算法,这里先咕着(

算法介绍

我们先考虑只有一个Sylow p-子群的情况,即$|G|=p^k$。首先可以$O(n^2)$求出每个元素的阶(这显然不是复杂度下界,不过更低也没啥意义)。算法的基本思路是按照阶从大到小求出每个基,在求的过程中维护与当前基线性相关的元素集合,然后在剩下的集合里选出新的基。

记当前基的集合为$B$,$L$为$B$生成的子群,$M$为所有跟当前基线性相关的元素集合。我们有引理:对于$x\in M,x\notin L$,且$|x|\leq\min_{b\in B}|b|$,一定可以找到$x'\notin M$,使得$\prod_{b\in B\cup{\{x\}}}\langle b\rangle=\prod_{b\in B\cup{\{x'\}}}\langle b\rangle$,且$|x'|\big||x|$。

这个引理说明了,我们贪心在$G\setminus M$中选取阶最大的元素,不会产生$L\ne G$却找不到新的基的情况,即证明了该算法的正确性。

当$G$含有多个Sylow p-子群的情况呢?事实上这个算法依然是对的,因为每个Sylow p-子群内部发生的事情实际上可以看作是独立的,我们每次找出阶最大的元素其实等价于在每个Sylow p-子群内都找出一个阶最大的元素。我们在找出最大阶元素的同时对阶进行一下质因数分解,求出每个Sylow p-子群对应的基即可。

整理一下算法的具体流程:

  1. 求出每个元素的阶
  2. 维护$B,L,M$,初始时$B=\varnothing,L=M=\{e\}$
  3. 在$G\setminus M$中找到阶最大的元素$g$,对$|g|$进行质因数分解,求出每个Sylow p-子群对应的基$g_{p_1},g_{p_2},\cdots$,并加入$B$中
  4. 每当$B$中插入新的基,枚举原有的$L$进行运算,扩展出新的$L$(设新插入的基为$b$,枚举$x\in L,1\leq k<|b|$,将$xb^k$加入$L$中)。每当$L$中有新元素$x$加入,枚举所有$y$使得$y\notin M$且$x$可表示为$y$的幂,把$y$加入$M$中。
  5. 若$L=G$,则找到了所有基,算法结束;否则回到第$3$步

简単でしょう?(本気)

复杂度分析

每一轮扩展$L$的复杂度为$O(|L||b|)$,可以看作是$O(n|b|)$,而所有基的阶之和显然为$O(n)$,所以这一步的复杂度为$O(n^2)$。而每个元素只会被加入$L$中至多一次,所以后面那一步的复杂度也为$O(n^2)$。因此算法的总时间复杂度为$O(n^2)$。

算法实装

ないです。やめたらこの仕事!!

代码的话,等后面更例题的时候再贴吧,应当是会有的(心虚)

查看详细内容

[LOJ 2476] [2018 集训队互测 Day 3] 蒜头的奖杯

2018-7-3 3:152018-7-3 8:9
集训队互测LOJ莫比乌斯反演

题目链接

这似乎是SDOI那几道题的强化版?

考虑枚举$i$,并对$(i,j),(i,k)$进行反演(本文用$(x,y)$表示$x,y$的$\gcd$): $$ \begin{aligned} &\sum_{i,j,k}A_iB_jC_kD_{(i,j)}E_{(i,k)}F_{(j,k)}\\ =&\sum_{i,j,k} A_iB_jC_kF_{(j,k)}\sum_{p|(i,j)}(\mu D)_p\sum_{q|(i,k)}(\mu E)_q\\ =&\sum_{p,q}(\mu D)_p(\mu E)_q\sum_{\mathbb{lcm}(p,q)|i}A_i\sum_{p|j}B_j\sum_{q|k}C_kF_{(j,k)}\\ &(记d=(p,q),x=p/d,y=q/d,n'=\lfloor n/d\rfloor,\\ &A'_d=\sum_{d|x,x\leq n}A_x,\\ &D'=\mu D,E'=\mu E)\\ =&\sum_d\sum_{xy\leq n',(x,y)=1}D'_{dx}E'_{dy}A'_{dxy}\sum_{dx|j}\sum_{dy|k}B_jC_kF_{(j,k)} \end{aligned} $$

假设$d,x,y$已经确定,且$x\leq y$,考虑计算后面的求和: $$ \begin{aligned} &\sum_{dx|j}\sum_{dy|k}B_jC_kF_{(j,k)}\\ =&\sum_{j=1}^{n'}\sum_{k=1}^{n'}[x|j,y|k]B_{dj}C_{dk}F_{d(j,k)}\\ =&\sum_{j=1}^{n'}\sum_{k=1}^{n'}[x|j,y|k]B'_jC'_kF'_{(j,k)}\\ =&\sum_{y|k}C'_k\sum_{d|k}(\mu F')_d\sum_{d|j}[x|j]B'_j \end{aligned} $$

我们暴力枚举所有的$d,x$。上面这坨东西只涉及到对因子求和以及对倍数求和,通过枚举质因子,我们可以每次用$O(n/d\log \log n/d)$的时间计算这坨东西对于每个$y$的值。这样复杂度就是$O(\sum_dn/d * \sqrt{n/d}\log \log n/d)=O(n^{1.5}\log \log n)$。中间那个根号是因为每次要在$\sqrt{n/d}$的范围内枚举$x$的取值。$x>y$的情况同理,求和式子里把$j$提到前面,再每次暴力枚举$y$就行了。

#include <bits/stdc++.h>
#define MAXN 100010
#define ULL unsigned long long

ULL read(){
    char c;
    while((c=getchar())<'0' || c>'9');
    ULL res=c-'0';
    while((c=getchar())>='0' && c<='9') res=res*10+c-'0';
    return res;
}

int n,n1;
int p[MAXN],np;
ULL A[MAXN],B[MAXN],C[MAXN],D[MAXN],E[MAXN],F[MAXN];
ULL B1[MAXN],C1[MAXN],F1[MAXN],S[MAXN];

void init(){
    static bool flag[MAXN];
    for(int i=2;i<=n;i++){
        if(!flag[i]) p[++np]=i;
        for(int j=1;j<=np && i*p[j]<=n;j++){
            flag[i*p[j]]=1;
            if(i%p[j]==0) break;
        }
    }
}

//s_x = \sum_{d|x} f_d
void s1(ULL *a,int n){
    for(int i=1;i<=np && p[i]<=n;i++)
        for(int j=p[i],k=1;j<=n;j+=p[i],k++)
            a[j]+=a[k];
}

//s_x = \sum_{x|d} f_d
void s2(ULL *a,int n){
    for(int i=1;i<=np && p[i]<=n;i++)
        for(int j=n/p[i],k=j*p[i];j>=1;j--,k-=p[i])
            a[j]+=a[k];
}

//s_x = \sum_{d|x} \mu(x/d)*f_d
void s3(ULL *a,int n){
    for(int i=1;i<=np && p[i]<=n;i++)
        for(int j=n/p[i],k=j*p[i];j>=1;j--,k-=p[i])
            a[k]-=a[j];
}

int gcd(int x,int y){
    while(y){ int t=x%y; x=y; y=t; }
    return x;
}

void gao_x(int x,ULL *B1,ULL *C1){
    memset(S,0,sizeof(S[0])*(n1+1));
    for(int i=x;i<=n1;i+=x) S[i]=B1[i];
    s2(S,n1);
    for(int i=1;i<=n1;i++) S[i]*=F1[i];
    s1(S,n1);
    for(int i=1;i<=n1;i++) S[i]*=C1[i];
    s2(S,n1);
}

ULL solve_d(int d){
    n1=n/d;
    for(int i=1;i<=n1;i++) B1[i]=B[i*d],C1[i]=C[i*d],F1[i]=F[i*d];
    s3(F1,n1);
    ULL res=0;
    for(int x=1;x*x<=n1;x++){
        gao_x(x,B1,C1);
        for(int y=x;y*x<=n1;y++)
            if(gcd(x,y)==1)
                res+=S[y]*A[d*x*y]*D[d*x]*E[d*y];
    }
    for(int y=1;y*y<=n1;y++){
        gao_x(y,C1,B1);
        for(int x=y+1;y*x<=n1;x++)
            if(gcd(x,y)==1)
                res+=S[x]*A[d*x*y]*D[d*x]*E[d*y];
    }
    return res;
}

int main(){
#ifdef DEBUG
    freopen("2476.in","r",stdin);
#endif
    n=read();
    for(int i=1;i<=n;i++) A[i]=read();
    for(int i=1;i<=n;i++) B[i]=read();
    for(int i=1;i<=n;i++) C[i]=read();
    for(int i=1;i<=n;i++) D[i]=read();
    for(int i=1;i<=n;i++) E[i]=read();
    for(int i=1;i<=n;i++) F[i]=read();
    init();
    s2(A,n);
    s3(D,n); s3(E,n);
    ULL ans=0;
    for(int d=1;d<=n;d++)
        ans+=solve_d(d);
    printf("%llu\n",ans);
    return 0;
}

查看详细内容

[UOJ #356] [JOI2017春季合宿] Port Facility

2018-6-5 7:132018-6-5 7:16
UOJJOI并查集线段树

题目链接

本题的限制等价于把一堆区间分成两个集合,每个集合内的区间两两之间要么不相交,要么完全包含。

我们把限制建成一个无向图,每两个冲突的区间之间都连一条边。显然,如果建出来不是二分图,答案必定为$0$(一个奇环上无法黑白染色使得两两相邻不同色)。我们把孤立的点去掉(去掉每个点会使答案乘$2$),考虑左边的每个点$x$,记右边与它相邻的点为$S_x$,则$S_x$中的点肯定是同色的,且$x$的颜色由$S_x$的颜色唯一确定。所以我们只需要把每个$S_x$用并查集并起来,最后统计右边的连通块数就行了。

然而直接建图搞的复杂度是$O(n^2)$的,我们考虑优化。注意到二分图的两边都可以看成一个合法的括号序列,我们建出对应的树结构。设右边所有区间中包含$x$的左端点及右端点的最小区间分别为$y_1,y_2$,则$S_x$即为$y_1$到$y_2$路径上的所有点减去两个点的LCA。这个东西在树上用并查集随便搞搞就行了。

关于求初始的二分图,我的做法是直接广搜,搜完之后检查一下合法性。对区间左端点建一个线段树,线段树每个区间按右端点的顺序维护未访问的所有区间。这里可以做到$O(n\log n)$的复杂度,加上前面的求LCA,总复杂度为$O(n\log n)$。

#include <bits/stdc++.h>
#define MAXN 2000010
using namespace std;

const int P=1000000007;

int read(){
    char c;
    while((c=getchar())<'0' || c>'9');
    int res=c-'0';
    while((c=getchar())>='0' && c<='9') res=res*10+c-'0';
    return res;
}

int n;
pair<int,int> p[MAXN];
int p1[MAXN],pt1[MAXN],pt2[MAXN];
int *s[MAXN*3],len[MAXN*3],cur[MAXN*3];
int c[MAXN],pre[MAXN][21],dep[MAXN],top[MAXN];
int f[MAXN];
queue<int> Q;

int getf(int x){
    if(f[x]==x) return x;
    return f[x]=getf(f[x]);
}

void build(int x,int l,int r){
    if(l==r){
        if(p1[l]){
            s[x]=new int[len[x]=1];
            s[x][0]=p1[l];
        }
        return;
    }
    int mid=(l+r)>>1;
    build(x<<1,l,mid);
    build(x<<1|1,mid+1,r);
    int x0=0,x1=0,x2=0,l1=len[x<<1],l2=len[x<<1|1];
    s[x]=new int[len[x]=l1+l2];
    while(x1<l1 || x2<l2){
        if(x1==l1) s[x][x0++]=s[x<<1|1][x2++];
        else if(x2==l2) s[x][x0++]=s[x<<1][x1++];
        else if(p[s[x<<1][x1]].second>p[s[x<<1|1][x2]].second) s[x][x0++]=s[x<<1][x1++];
        else s[x][x0++]=s[x<<1|1][x2++];
    }
}

void extend(int x,int cl,int cr,int l,int r,int nc){
    if(l<=cl && cr<=r){
        while(cur[x]<len[x] && p[s[x][cur[x]]].second>r){
            int t=s[x][cur[x]++];
            if(c[t]==-1){
                c[t]=nc;
                Q.push(t);
            }
        }
        return;
    }
    int mid=(cl+cr)>>1;
    if(l<=mid) extend(x<<1,cl,mid,l,r,nc);
    if(r>mid) extend(x<<1|1,mid+1,cr,l,r,nc);
}

void bfs(int x){
    Q.push(x);
    c[x]=0;
    while(!Q.empty()){
        int x=Q.front();
        Q.pop();
        extend(1,1,n<<1,p[x].first,p[x].second,c[x]^1);
        extend(1,1,n<<1,1,p[x].first,c[x]^1);
    }
}

bool check(){
    static stack<int> S;
    for(int _i=1;_i<=n<<1;_i++){
        while(!S.empty() && p[S.top()].second<_i) S.pop();
        if(p1[_i] && !c[p1[_i]]){
            int i=p1[_i];
            if(!S.empty() && p[S.top()].second<p[i].second) return 0;
            if(!S.empty()) pre[i][0]=S.top();
            S.push(i);
            dep[i]=dep[pre[i][0]]+1;
            top[i]=pre[i][0]?top[pre[i][0]]:i;
        }
        if(!S.empty()) pt1[_i]=S.top();
    }
    while(!S.empty()) S.pop();
    for(int _i=1;_i<=n<<1;_i++){
        while(!S.empty() && p[S.top()].second<_i) S.pop();
        if(p1[_i] && c[p1[_i]]){
            int i=p1[_i];
            if(!S.empty() && p[S.top()].second<p[i].second) return 0;
            if(!S.empty()) pre[i][0]=S.top();
            S.push(i);
            dep[i]=dep[pre[i][0]]+1;
            top[i]=pre[i][0]?top[pre[i][0]]:i;
        }
        if(!S.empty()) pt2[_i]=S.top();
    }
    return 1;
}

void pre_gao(){
    for(int i=1;i<=20;i++)
        for(int j=1;j<=n;j++)
            pre[j][i]=pre[pre[j][i-1]][i-1];
}

int getLCA(int x,int y){
    if(top[x]^top[y]) return 0;
    if(dep[x]>dep[y]) swap(x,y);
    int dd=dep[y]-dep[x];
    for(int i=0,l=1;l<=dd;i++,l<<=1)
        if(dd&l) y=pre[y][i];
    if(x==y) return x;
    for(int i=20;i>=0;i--)
        if(pre[x][i]^pre[y][i])
            x=pre[x][i],y=pre[y][i];
    return pre[x][0];
}

void merge(int x,int y){
    int y0=y;
    y=getf(y);
    while((x=getf(x))^y){
        if(pre[x][0]^y0) f[x]=pre[x][0];
        x=pre[x][0];
    }
}

void query1(int x){
    printf("%d %d\n",pt2[p[x].first],pt2[p[x].second]);
}

void query2(int x){
    printf("%d %d\n",pt1[p[x].first],pt1[p[x].second]);
}

int gao(){
    int ans=1;
    static vector<pair<int,int> > qs;
    for(int i=1;i<=n;i++)
        if(!c[i]){
            int t1=pt2[p[i].first],t2=pt2[p[i].second];
            if((!t1 && !t2) || t1==t2) ans=ans*2%P;
            else if(!t1 || !t2){
                merge(t1+t2,0);
            }else{
                int lca=getLCA(t1,t2);
                if(t1==lca) merge(t2,t1);
                else if(t2==lca) merge(t1,t2);
                else{
                    merge(t1,lca); merge(t2,lca);
                    qs.push_back(make_pair(t1,t2));
                }
            }
        }
    for(auto i:qs)
        if(getf(i.first)^getf(i.second))
            f[f[i.first]]=f[i.second];
    for(int i=1;i<=n;i++)
        if(c[i] && getf(i)==i)
            ans=ans*2%P;
    return ans;
}

int main(){
#ifdef DEBUG
    freopen("356.in","r",stdin);
#endif
    n=read();
    for(int i=1;i<=n;i++) f[i]=i;
    for(int i=1;i<=n;i++) p[i].first=read(),p[i].second=read();
    sort(p+1,p+n+1);
    for(int i=1;i<=n;i++) p1[p[i].first]=i;
    build(1,1,n<<1);
    memset(c,-1,sizeof c);
    for(int i=1;i<=n;i++)
        if(c[i]==-1) 
            bfs(i);
    if(!check()){
        puts("0");
        return 0;
    }
    pre_gao();
    printf("%d\n",gao());
}

查看详细内容

自动脸部替换工具的构想与实现

2018-5-19 2:312018-5-20 13:13
python图像处理

最近比较有空,这几天又一直在P图,感觉换脸的过程其实就是一个枯燥的重复过程,于是就产生了做一个自动换脸工具的想法。


先来看一下最终的效果:

效果图

やりますねぇ!(赞赏)


准备工作

要实现自动P脸,肯定会用到人脸的识别。这里我直接使用了github上的一个轮子:face_recognition,这里面用到了dlib,在安装它之前要先装好cmake,接着在命令行里用pip安装即可(当然还需要用到PIL和numpy,就假装已经装好了吧)

这个库可以支持在一张图片中识别人脸,并返回每张人脸的所有特征点。后面做特征匹配的时候用到的特征点就是用它找出来的。

人脸定位

现在我们有两张人脸,每张人脸都有许多不同的特征点...

feature

我们要对其中一张人脸进行变换,使得变换后的所有特征点与另一张人脸的特征点尽可能接近。这里的变换只涉及到旋转,平移,缩放,平行扭曲(仿射变换),可以用$6$个参数$a,b,c,d,e,f$来表示:

$$x_1=ax_0+by_0+c,y_1=dx_0+ey_0+f$$

我们定义代价函数为$S=\frac 12\sum_i{({x_1}_i-{x_2}_i)^2+({y_1}_i-{y_2}_i)^2}$,其中${x_1}_i,{y_1}_i$是人脸$1$经过变换后的第$i$个特征点的坐标,${x_2}_i,{y_2}_i$是人脸$2$中第$i$个特征点的坐标。这里乘上$\frac 12$是为了后面求导方便。

我们代入变换的公式:

$$ S=\frac 12\sum(ax_0+by_0+c-x_2)^2+(dx_0+ey_0+f-y_2)^2 $$

这里省略了下标$i$。我们对六个参数求偏导:

$$ \frac{\partial S}{\partial a}=\sum(ax_0+by_0+c-x_2)x_0\\ \frac{\partial S}{\partial b}=\sum(ax_0+by_0+c-x_2)y_0\\ \frac{\partial S}{\partial c}=\sum ax_0+by_0+c-x_2\\ \frac{\partial S}{\partial d}=\sum(dx_0+ey_0+f-y_2)x_0\\ \frac{\partial S}{\partial e}=\sum(dx_0+ey_0+f-y_2)y_0\\ \frac{\partial S}{\partial f}=\sum dx_0+ey_0+f-y_2\\ $$

本来想直接梯度下降最优化这个东西的,后来发现令所有偏导数都为$0$,可以直接得到关于参数的线性方程组,直接解方程组就能得到最优的仿射变换参数了。具体方程组为:

$$ \begin{cases} (\sum x_0^2)a&+(\sum y_0x_0)b&+(\sum x_0)c&=\sum x_2x_0\\ (\sum x_0y_0)a&+(\sum y_0^2)b&+(\sum y_0)c&=\sum x_2y_0\\ (\sum x_0)a&+(\sum y_0)b&+nc&=\sum x_2\\ \end{cases} $$

($d,e,f$类似)

这样我们就能在$O(n)$的时间内求出对应的仿射变换了(虽然这里对性能并没有什么要求。。因为特征点并不多)。找到仿射变换后,可以用PIL.Image.transform把变换作用到图片上。

为了解决面部中心特征点过于密集的问题,我们还可以让特征点带上权值,降低较密集的特征点的权值。具体可以看代码。

我们把所有经过对应变换的人脸叠加在一张新图中,后面需要把这张图和原图进行合成。

图片合成

现在我们有两张图片:底图(记为图1)和一张由新的人脸组成的图(记为图2)。

一个比较直接的方法是直接把图2按alpha通道混合在图1中(这意味着什么都没有处理)。这样的话新图的颜色可能会跟原图对不上,人脸的边缘也会比较粗糙。

为了解决这两个问题,我们做两个额外的步骤:颜色调整和边缘羽化。

为了让颜色更加自然,我们需要调整图2中的颜色,让其接近图1的同时尽可能地保留细节。解决方案是把图1对图2遮罩得到图3(即削除了图1中不含人脸的部分),把图2和图3分别高斯模糊得到图4和图5,把图2对应像素除以图4再乘以图5,就能得到这个像素调整后的颜色。这么做是因为高斯模糊消除了图片细节,所以这样调整后即不会表现出覆盖掉的内容中的细节,也能体现其原有的大体颜色。

处理完颜色后,我们还需要对边缘进行羽化,让过渡更加自然。我没想到什么特别好的方法,直接对图2的alpha通道做高斯模糊后,再令新的alpha值$\alpha'$为$\min(0,\max(1,(\alpha-0.5)\times 2))$。这样做的原理是把让边缘的$alpha$值取$0$,且往内部扩展时$alpha$值变大,在边缘平滑的时候效果比较好,因为平滑边缘处的高斯模糊值可近似看作$0.5$,但在边缘鬼畜的时候可能会有一些奇妙问题。。不过一般情况下效果还是挺不错的。


具体实现

https://github.com/Alif-01/AutoHeadReplace

求一波star QwQ

神秘链接

查看详细内容

FFT相关的一些优化技巧

2018-4-17 21:132018-4-17 21:13
FFT卡常数

大概就是myy论文里的那些东西吧,闲着没事整理一下。

记FFT变换的矩阵为$F$(即单位根的范德蒙德矩阵),将FFT视作$F$乘上待变换的行向量:$FFT(a)=Fa$。我们记$\overline F$为$F$的复共轭矩阵(矩阵的每个元素都共轭),$a'$为$a$除了第$0$项外都翻转后的结果(本文的所有矩阵、向量的下标均从$0$开始),即第$a_i$与$a_{n-i}$交换($n$为$a$的长度)。显然$a'$可以写成$a$乘上一个翻转矩阵$T$,即$a'=aT$。

  • 性质 1:$\overline Fa=(Fa)'$

这个根据单位根的性质$\omega^{-k}=\omega^{n-k}$可以轻易得到。

  • 性质 2:$F^{-1}=\overline F/n$

这其实就是IFFT最后一步干的事情。。

技巧 1:快速求两个实多项式的点值(以及插值)

现在有两个实向量$a,b$,我们想计算$Fa,Fb$。按照一般的做法,直接跑两遍FFT就行了。其实这个过程可以做到只用一次FFT。

我们构造复向量$a+bi$,考虑把这个向量FFT: $$F(a+bi)=Fa+F(bi)=Fa+iFb$$

我们把这整个东西取复共轭,看看会发生什么: $$ \begin{aligned} \overline{F(a+bi)}&=\overline {Fa}+\overline {iFb}\\ &=\overline Fa-i\overline Fb\\ &=(Fa)'-i(Fb)'\\ \end{aligned} $$

两边同时乘上$T$,再联立之前FFT的结果: $$ Fa-iFb=(\overline{F(a+bi)})'\\ Fa+iFb=F(a+bi) $$

知道这两个东西之后,可以直接解出$Fa,Fb$。也就是说,我们直接算出$a+bi$的FFT,把它复共轭之后翻转系数(从第$1$项开始翻转),之后就可以解出$Fa$和$Fb$了。

同理,我们也可以同时做两个实多项式点值的插值。把这整个过程反过来做一遍即可。

技巧 2:4次FFT的任意模数多项式乘法

这其实是上面那个东西的一个应用。现在有两个整系数多项式$A(x),B(x)$,我们要算它们循环卷积的结果(如果要做普通卷积,直接把次数扩大一倍即可),所有系数都对一个$10^9$级别的整数$P$取模(这个整数可以什么性质都没有)。我们考虑选取整数$M\approx\sqrt P$,我们把$A,B$的系数对$M$作带余除法,得到$A(x)=A_0(x)+MA_1(x),B(x)=B_0(x)+MB_1(x)$。这里$A_0,A_1,B_0,B_1$的系数都是$M$级别的,他们之间乘起来之后系数是$O(M^2n)=O(Pn)$级别的,$n$不太大的话(大概$5*10^5$以内)可以认为是在double的精度范围内,可以直接用double做乘法后取模。

我们用上面提到的那个技巧,用$2$次FFT算出$A_0,A_1,B_0,B_1$的FFT点值$a_0,a_1,b_0,b_1$,之后根据$$\begin{aligned}A(x)B(x)&=A_0(x)B_0(x)\\&+M(A_0(x)B_1(x)+A_1(x)B_0(x))\\&+M^2A_1(x)B_1(x)\end{aligned}$$我们可以再次使用上面那个技巧,用$2$次FFT算出$a_0\otimes b_0,a_0\otimes b_1,a_1\otimes b_0,a_1\otimes b_1$的IFFT插值,就可以直接进行计算了。

总共$4$次FFT。

技巧 3:2次FFT的实多项式乘法

这个用技巧$1$也可以做到,但较技巧$1$更容易实现。
有两个实多项式$A(x),B(x)$,我们要计算它们的循环卷积。

我们可以用两次FFT计算$(A(x)+iB(x))^2=A^2(x)-B^2(x)+2iA(x)B(x)$。

直接取虚部并除以2即可。

查看详细内容