ZJT's Blog

大変に気分がいい~

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

2018-7-2 23:152018-7-3 4: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 3:132018-6-5 3: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-18 22:312018-5-20 9: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 17:132018-4-17 17: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即可。

查看详细内容

动态图连通性问题的在线解法

2018-4-16 22:252018-4-16 22:25
动态图Euler-tour Tree毒瘤

最近居然在省选模拟赛里面做到了这个毒瘤问题。。

动态图连通性问题是指维护一个不断加边删边的无向图,要求查询两个点之间的连通性,以及支持对连通块的各种信息的维护。离线的话非常好做,直接处理出每条边的删除时间后,维护关于删除时间的最大生成树即可。强制在线的话,也有一个$O(\log^2n)$的做法,是通过将图分层并用ETT维护每一层的生成树来解决的。不过这个算法实现起来比较复杂,第一次写完+调试完大概花了6h,估计场上想写是不太可能的了。。

为了维护连通性的信息,我们可以考虑维护生成森林。我们给每条边钦点一个非负整数的等级(level),定义$G_i$为所有等级为$i$的边构成的子图,$T_i$为所有等级大于等于$i$的边构成的图中,关于等级的最大生成森林。最终询问的时候,直接在$T_0$上做询问就行了。

当一条新边加入时,我们先令它的等级为$0$。这时,我们会在$G_0$中加入这条边,如果$T_0$中加入这条边不会形成环,我们也把它加入$T_0$。

删除的情况比较复杂。如果删掉的是一条不在$T_0$中的边,那它也不可能在任何$T_i$中出现,我们直接在对应等级的$G_i$中把它删掉就行了。然而如果我们要删掉一条树边,则可能存在一条非树边,使得删掉这条树边形成两个连通块可以被这条非树边连接。这时候这条非树边就会取代删去的边,成为新的树边。

考虑如何在图中找到这条非树边。我们记删去的边为$e_0$,等级为$l_0$,替换它的非树边为$e_1$,等级为$l_1$。显然$l_1\leq l_0$,否则最大生成森林的性质就不满足了($e_1$可以替代$e_0$得到一个更大的生成森林)。

我们考虑从$l_0$到$0$枚举$l_1$,看看是否存在等级为$l_1$的非树边$e_1$能替换删去的$e_0$。一个比较显然的结论是$e_1$连接的两个端点一定在$T_{l_1}$上与$e_0$的端点连通,这个可以通过MST的性质得到。记$e_0$在$T_{l_1}$上连接的两个连通块为$X,Y$,且$|X|\leq|Y|$,考虑一个暴力做法,我们直接在$G_{l_1}$中枚举每条从$X$连出去的边,如果连出去的边位于$Y$中,则找到了$e_1$。如果枚举完每条边之后还是没找到,则说明不存在等级为$l_1$的满足条件的$e_1$。

这样做的复杂度显然是错的,我们考虑对它进行优化。我们在暴力枚举每条边的时候 ,如果这条边不能满足条件,就把它的等级$+1$(在$G_{l_1}$中删去它,在$G_{l_1+1}$中加入它,如果它是树边还要在$T_{l_1+1}$中也加入它)。这样每次访问到一条边时,要么找到了答案,要么这条边的等级增加,所以我们暴力枚举的复杂度就只跟每条边的最大等级有关了(等级增加次数肯定不会超过最大等级)。同时,我们还可以证明,$T_i$中的最大连通块大小是$O(n/2^i)$级别的。因为我们每次都是只处理从较小的连通块中连出的边($|X|\leq|Y|$),且在我们把所有边升级完之后,下一等级的生成树上形成的连通块大小也不会超过$|X|$(因为当前层的生成树一定包含下一层的生成树),所以每一层的连通块大小的上界都会减半。这也说明了边的最大等级是$O(\log n)$级别的,即每条边都最多只会被升级$O(\log n)$次。

我们用ETT维护$T_i$中有哪些点有$G_i$的连边,这样每次搜索的复杂度就是$O($搜到的边数$\times\log n)$的,在ETT中的link和cut也是$O(\log n)$的,所以总复杂度为$O(m\log^2 n)$。这样就成功解决了这个问题。

实现上可能有一些细节问题要注意,首先在搜索$X$的出边时,要先搜所有的树边。如果先搜了非树边,可能会因为非树边等级的提高而出现生成树心态发生改变的情况,这种情况是不好处理的。其次对边的搜索要一边搜索一边检查是否可行,如果先搜出了所有边再一个个判断,复杂度就是错的,在loj上会被卡(

这是模板题

附上代码实现一份(loj垫底qwq):

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

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,m,online;
int v[MAXN];

namespace Duliu{
    int ls[MAXN],numls;

    struct ETT{
        struct node{
            int lc,rc,p,key,size,pt;
            bool tag,stag;
        }a[MAXN<<1];

        int blist[MAXN],numb;
        int p[MAXN];
        map<int,int> e[MAXN];
        int tagp[MAXN];

        ETT(int n){
            for(int i=2*n;i>=1;i--) a[blist[++numb]=i].key=rand();
        }

        int newnode(int pt){
            int x=blist[numb--];
            a[x].p=a[x].lc=a[x].rc=a[x].tag=a[x].stag=0;
            a[x].size=1;
            a[x].pt=pt;
            return x;
        }

        void delnode(int x){
            blist[++numb]=x;
        }

        inline void pushup(int x){
            a[x].stag=a[a[x].lc].stag||a[a[x].rc].stag||a[x].tag;
            a[x].size=a[a[x].lc].size+a[a[x].rc].size+1;
        }

        int merge(int x,int y){
            if(!x || !y) return x+y;
            if(a[x].key<a[y].key){
                a[x].rc=merge(a[x].rc,y);
                a[a[x].rc].p=x;
                pushup(x);
                return x;
            }else{
                a[y].lc=merge(x,a[y].lc);
                a[a[y].lc].p=y;
                pushup(y);
                return y;
            }
        }

        pair<int,int> split(int x,int d){
            int s1=a[x].lc,s2=a[x].rc;
            if(d==-1){
                a[x].rc=0;
                pushup(x);
                s1=x;
                a[s2].p=0;
            }else if(d==1){
                a[x].lc=0;
                pushup(x);
                s2=x;
                a[s1].p=0;
            }else a[s1].p=a[s2].p=0;
            while(a[x].p){
                int p=a[x].p;
                if(a[p].lc==x){
                    a[p].lc=s2;
                    a[s2].p=p;
                    a[s1].p=0;
                    s2=p;
                }else{
                    a[p].rc=s1;
                    a[s1].p=p;
                    a[s2].p=0;
                    s1=p;
                }
                pushup(x=p);
            }
            return make_pair(s1,s2);
        }

        int getRoot(int x){
            while(a[x].p) x=a[x].p;
            return x;
        }

        int get_p(int x){
            if(e[x].empty()) return 0;
            return e[e[x].begin()->first][x];
        }

        void link(int x,int y){
            int p1=newnode(y),p2=newnode(x);
            if(tagp[x]==-1){
                tagp[x]=p2;
                a[p2].tag=a[p2].stag=1;
            }
            if(tagp[y]==-1){
                tagp[y]=p1;
                a[p1].tag=a[p1].stag=1;
            }
            int t1,t2,t3;
            if(e[x].empty()) t1=t3=0;
            else{
                pair<int,int> temp=split(e[x].begin()->second,1);
                t1=temp.first; t3=temp.second;
            }
            if(e[y].empty()) t2=0;
            else{
                pair<int,int> temp=split(e[y].begin()->second,1);
                t2=merge(temp.second,temp.first);
            }
            e[x][y]=p1; e[y][x]=p2;
            merge(merge(merge(t1,p1),t2),merge(p2,t3));
        }

        void cut(int x,int y){
            int p1=e[x][y],p2=e[y][x];
            e[x].erase(e[x].find(y));
            e[y].erase(e[y].find(x));
            pair<int,int> t1=split(p1,0);
            if(getRoot(p2)==t1.first){
                pair<int,int> t2=split(p2,0);
                merge(t2.first,t1.second);
            }else{
                pair<int,int> t2=split(p2,0);
                merge(t1.first,t2.second);
            }
            if(a[p1].tag){
                if(e[y].empty()) tagp[y]=-1;
                else{
                    tagp[y]=get_p(y);
                    a[tagp[y]].tag=1;
                    for(int i=tagp[y];i;i=a[i].p) pushup(i);
                }
            }
            if(a[p2].tag){
                if(e[x].empty()) tagp[x]=-1;
                else{
                    tagp[x]=get_p(x);
                    a[tagp[x]].tag=1;
                    for(int i=tagp[x];i;i=a[i].p) pushup(i);
                }
            }
            delnode(p1); delnode(p2);
        }

        void setTag(int x,int v){
            if(v==1){
                int p=get_p(x);
                if(!p){
                    tagp[x]=-1;
                    return;
                }
                x=tagp[x]=p;
            }else{
                int t=tagp[x];
                tagp[x]=0;
                x=t;
                if(x==-1){
                    tagp[x]=0;
                    return;
                }
            }
            a[x].tag=v;
            while(x){
                pushup(x);
                x=a[x].p;
            }
        }

        int get_size(int x){
            if(!x) return 1;
            return a[getRoot(x)].size/2+1;
        }

        void dfs1(int x,int i);
        bool dfs2(int x,int i,int py);
    }*s[MAXL];

    int l;
    map<int,int> gl[MAXN];
    set<int> g[MAXL][MAXN];
    pair<int,int> res_e;

    void init(int n){
        for(l=0;(1<<l)<n;l++);
        for(int i=0;i<=l;i++) s[i]=new ETT(n);
    }

    bool link(int x,int y){
        gl[x][y]=0;
        gl[y][x]=0;
        bool flagx=g[0][x].empty(),flagy=g[0][y].empty();
        g[0][x].insert(y);
        g[0][y].insert(x);
        int px=s[0]->get_p(x),py=s[0]->get_p(y);
        if(!px || !py || s[0]->getRoot(px)!=s[0]->getRoot(py)){
            s[0]->link(x,y);
            if(flagx) s[0]->setTag(x,1);
            if(flagy) s[0]->setTag(y,1);
            return 1;
        }
        if(flagx) s[0]->setTag(x,1);
        if(flagy) s[0]->setTag(y,1);
        return 0;
    }

    void level_up(int x,int y,int curl){
        if(g[curl][x].empty()) s[curl]->setTag(x,0);
        if(g[curl][y].empty()) s[curl]->setTag(y,0);
        gl[x][y]++;
        gl[y][x]++;
        bool flagx=g[curl+1][x].empty(),flagy=g[curl+1][y].empty();
        g[curl+1][x].insert(y);
        g[curl+1][y].insert(x);
        int px=s[curl+1]->get_p(x),py=s[curl+1]->get_p(y);
        if(!px || !py || s[curl+1]->getRoot(px)!=s[curl+1]->getRoot(py)) s[curl+1]->link(x,y);
        if(flagx) s[curl+1]->setTag(x,1);
        if(flagy) s[curl+1]->setTag(y,1);
    }

    void gao_1(int x,int i){
        auto it=g[i][x].begin();
        while(it!=g[i][x].end()){
            if(s[i]->e[x].count(*it)){
                int y=*it;
                g[i][y].erase(g[i][y].find(x));
                it=g[i][x].erase(it);
                level_up(x,y,i);
            }else it++;
        }
    }

    bool gao_2(int x,int i,int py){
        auto it=g[i][x].begin();
        while(it!=g[i][x].end()){
            if(s[i]->getRoot(s[i]->get_p(*it))==py){
                for(int k=0;k<=i;k++)
                    s[k]->link(x,*it);
                res_e=make_pair(x,*it);
                return 1;
            }else{
                int y=*it;
                g[i][y].erase(g[i][y].find(x));
                it=g[i][x].erase(it);
                level_up(x,y,i);
            }
        }
        return 0;
    }

    pair<int,int> cut(int x,int y){
        int curl=gl[x][y];
        gl[x].erase(gl[x].find(y));
        gl[y].erase(gl[y].find(x));
        g[curl][x].erase(g[curl][x].find(y));
        g[curl][y].erase(g[curl][y].find(x));
        if(s[0]->e[x].count(y)){
            for(int i=0;i<=curl;i++) s[i]->cut(x,y);
            for(int i=curl;i>=0;i--){
                int px=s[i]->get_p(x),py=s[i]->get_p(y);
                int sx=s[i]->get_size(px),sy=s[i]->get_size(py);
                if(sy<sx) swap(x,y),swap(sx,sy),swap(px,py);
                py=s[i]->getRoot(py);
                px=s[i]->getRoot(px);
                if(px){
                    s[i]->dfs1(px,i);
                    if(s[i]->dfs2(px,i,py)) return res_e;
                }else{
                    gao_1(x,i);
                    if(gao_2(x,i,py)) return res_e;
                }
            }
            return make_pair(0,0);
        }
        return make_pair(-1,0);
    }

    bool check(){
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++){
                bool flag=0;
                for(int k=l;k>=0;k--)
                    if(s[k]->e[i].count(j)) flag=1;
                    else if(flag) return 0;
            }
        return 1;
    }

    void ETT::dfs1(int x,int i){
        if(!a[x].stag) return;
        if(a[x].tag) gao_1(a[x].pt,i);
        if(a[a[x].lc].stag) dfs1(a[x].lc,i);
        if(a[a[x].rc].stag) dfs1(a[x].rc,i);
    }

    bool ETT::dfs2(int x,int i,int py){
        if(!a[x].stag) return 0;
        if(a[x].tag) if(gao_2(a[x].pt,i,py)) return 1;
        return dfs2(a[x].lc,i,py) || dfs2(a[x].rc,i,py);
    }
}

void link(int x,int y){ Duliu::link(x,y); }
void cut(int x,int y){ Duliu::cut(x,y); }

bool query(int x,int y){
    using namespace Duliu;
    if(x==y) return 1;
    int px=s[0]->get_p(x),py=s[0]->get_p(y);
    return px && py && s[0]->getRoot(px)==s[0]->getRoot(py);
}

int main(){
#ifdef DEBUG
    freopen("122.in","r",stdin);
#endif
    n=read(); m=read(); online=1;
    Duliu::init(n);
    int lastans=0;
    while(m--){
        int t=read(),x=read(),y=read();
        if(online){
            x^=lastans;
            y^=lastans;
        }
        if(t==0) link(x,y);
        if(t==1) cut(x,y);
        if(t==2){
            if(query(x,y)){
                lastans=x;
                puts("Y");
            }else{
                lastans=y;
                puts("N");
            }
        }
    }
    return 0;
}

查看详细内容