[CF1067E]Random Forest Rank

Link

给定一个森林,每条边有\frac 12的概率被删去,随机删掉一些边,求剩下的图的临界矩阵的秩的期望值。

Solution

满秩的矩阵的行列式不为0。结合行列式的定义|A|=\sum_{\pi_i}(-1)^{\delta(\pi_i)}\prod a_{\pi_{ij}}。对于一个行列式不为零的邻接矩阵,其必要条件是\exist \pi_i,\text{s.t.}\prod a_{\pi_{ij}}\not=0
\pi_i循环分解,循环内相邻的点都有边连接。然而,对于一个不存在环的森林,这样的循环大小只能为2——即一条边连接的两个点。因此,求出这个森林的最大匹配,设为m,那么就能选出一个大小为2m的点集,其邻接矩阵满足\exist \pi_i,\text{s.t.}\prod a_{\pi_{ij}}\not=0的条件。
另外,选出的这个点集,完美匹配是唯一的,因此,只会有一个\pi_i满足上式子。
由此,问题转化为求期望最大匹配。
然后就不会了 先坑着吧

Code

[loj 6053]简单的函数

Link

Solution

p xor 1并不是积性函数
但是除了p=2的时候,p^1=p-1
p和1都是积性函数,就可以筛出来然后加上2把p=2的情况修正好。

Code

const int P=1e9+7,inv2=500000004;

int Add(int x,int const &y) {
    return (x+=y)>=P?x-P:x;
}

int Minus(int x,int const &y) {
    return (x-=y)<0?x+P:x;
}

int Mul(long long x,int const &y) {
    return x*y%P;
}


const int N=2e5,XN=N+11;

int prime[XN],pcnt,psum[XN],fsum[XN];
void Prep() {
    static bool notPrime[XN];
    for(int i=2;i<=N;++i) {
        if(!notPrime[i])
            prime[++pcnt]=i;
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0)
                break;
        }
    }
    for(int i=1;i<=pcnt;++i)
        psum[i]=Add(psum[i-1],prime[i]);
}

int lim,psz;long long n;
struct Identifier {
    int id[2][XN],cnt;
    int &operator [](long long x) {
        int &res=x<=lim?id[0][x]:id[1][n/x];
        if(res==0)
            res=++cnt;
        return res;
    }
}id;


int F(int p,int t) {
    return p^t;
}

int g[XN*2];

int H(long long n,int m) {
    if(n<=1 || m>psz)
        return 0;
    int res=Minus(g[id[n]],fsum[m-1]);
    for(int i=m;i<=psz && (long long)prime[i]*prime[i]<=n;++i) {
        long long pt=prime[i],pt1=pt*prime[i];
        for(int t=1;pt1<=n;++t,pt=pt1,pt1*=prime[i])
            res=Add(res,Add(Mul(F(prime[i],t),H(n/pt,i+1)),F(prime[i],t+1)));
    }
    return res;
}

int main() {
    Prep();
    static long long kp[XN*2];
    int kpc=0;
    fin>>n;
    lim=sqrt(n)+0.5,psz=std::upper_bound(prime+1,prime+1+pcnt,lim)-prime;
    static int g1[XN*2],g2[XN*2];
    for(long long l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        kp[++kpc]=n/l;
        g1[id[n/l]]=Minus((n/l%P)*(n/l%P+1)%P*inv2%P,1);
        g2[id[n/l]]=Minus(n/l%P,1);
    }
    for(int j=1;j<=psz;++j)
        for(int i=1;i<=kpc && (long long)prime[j]*prime[j]<=kp[i];++i) {
            g1[id[kp[i]]]=Minus(g1[id[kp[i]]],Mul(prime[j],Minus(g1[id[kp[i]/prime[j]]],psum[j-1])));
            g2[id[kp[i]]]=Minus(g2[id[kp[i]]],Minus(g2[id[kp[i]/prime[j]]],j-1));
        }
    for(int i=1;i<=psz;++i)
        fsum[i]=Add(fsum[i-1],prime[i]^1);
    for(int i=1;i<=kpc;++i)
        g[i]=Add(Minus(g1[i],g2[i]),kp[i]>=2?2:0);
    fout<<Add(H(n,1),1)<<'\n';
    return 0;
}

[ICPC 2018 Xuzhou Online]Easy Math

Link

Solution

\sum_{i=1}^n\mu(ix)=\mu(x)\sum_{i=1}^n\mu(i)[\gcd(i,x)=1]

考虑计算

\begin{aligned}
&\sum_{i=1}^nf(i)\\
&f(i)=\mu(i)[\gcd(i,n)=1]
\end{aligned}

首先,

f(i)=\begin{cases}
1,&i=1\\
-1,&i=p{\ \rm and\ }p\not|x\\
0,&i=p^k{\ \rm or\ }p|x
\end{cases}

然后,分\gcd(i,x)的不同情况讨论,容易知道f(i)是积性函数。
因此可以使用min25筛法了。

Tips

还是不必拘泥于g的形式……

Code

const int N=5e4,XN=N+11;

int prime[XN],pcnt;
void Prep() {
    static bool notPrime[XN];
    for(int i=2;i<=N;++i) {
        if(!notPrime[i])
            prime[++pcnt]=i;
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0)
                break;
        }
    }
}

int n,lim,psz;long long x;
struct Identifier {
    int id[2][XN],cnt;
    int &operator [](int x) {
        int &res=x<=lim?id[0][x]:id[1][n/x];
        if(res==0)
            res=++cnt;
        return res;
    }
}id;

long long dv[100];int dt;
int Mu(long long x) {
    int res=1;long long cx=x;
    for(long long d=2;d*d<=cx && x!=1;++d)
        if(x%d==0) {
            x/=d;
            dv[++dt]=d;
            if(x%d==0)
                return 0;
            else
                res*=-1;
        }
    if(x!=1) {
        res*=-1;
        dv[++dt]=x;
    }
    return res;
}

int F(int p,int t) {
    return t==1 && x%p?-1:0;
}

int g[XN*2],fsum[XN];

long long H(int n,int m) {
    if(n<=1 || m>psz)
        return 0;
    long long res=g[id[n]]-fsum[m-1];
    for(int i=m;i<=psz && prime[i]*prime[i]<=n;++i) {
        long long pt=prime[i],pt1=pt*prime[i];
        for(int t=1;pt1<=n;++t,pt=pt1,pt1*=prime[i])
            res+=F(prime[i],t)*H(n/pt,i+1)+F(prime[i],t+1);
    }
    return res;
}

int main() {
    Prep();
    static int kp[XN*2],kpc;
    fin>>n>>x;
    int mun=Mu(x);
    if(mun==0)
        fout<<0<<'\n';
    else {
        lim=sqrt(n)+0.5,psz=std::upper_bound(prime+1,prime+1+pcnt,lim)-prime;
        for(int l=1,r;l<=n;l=r+1) {
            r=n/(n/l);
            g[id[kp[++kpc]=n/l]]=n/l-1;
        }
        for(int j=1;j<=psz;++j)
            for(int i=1;i<=kpc && prime[j]*prime[j]<=kp[i];++i)
                g[id[kp[i]]]-=g[id[kp[i]/prime[j]]]-(j-1);//(1-(j-Count(j)));
        for(int i=1;i<=psz;++i)
            fsum[i]=fsum[i-1]+F(prime[i],1);
        for(int i=1;i<=kpc;++i)
            g[i]=std::upper_bound(dv+1,dv+1+dt,kp[i])-dv-1-g[i];
        fout<<mun*(H(n,1)+1)<<'\n';
    }
    return 0;
}

SPOJ DIVCNTK

Link

Solution

使用min_25筛
由于求h的过程都只会用到[\frac ni]集合中的点,那么求g的过程中把[\frac ni]构成的集合求出来,对于这个集合的所有点分层递推(递推过程中需要用到的g也属于[\frac ni]集合),就可以求出g(i,|P|)
然后递归求h,并不需要记忆化。边界见代码。

Code

typedef unsigned long long qword;

const int N=2e5,XN=N+11;

int prime[XN],pcnt;
void Prep() {
    static bool notPrime[XN];
    for(int i=2;i<=N;++i) {
        if(!notPrime[i])
            prime[++pcnt]=i;
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0)
                break;
        }
    }
}

long long n,k;
int lim,psz;

struct Identifier {
    int id[2][XN*2],cnt;
    int &operator [](long long x) {
        int &res=x<=lim?id[0][x]:id[1][n/x];
        if(res==0)
            res=++cnt;
        return res;
    }
}id;

qword g[XN*2];

qword H(long long n,int m) {
    if(n<=1 || m>psz)
        return 0;
    qword res=(k+1)*(g[id[n]]-m+1);
    for(int i=m;i<=psz && (long long)prime[i]*prime[i]<=n;++i) {
        long long pt=prime[i];
        for(int t=1;pt*prime[i]<=n;++t,pt*=prime[i])
            res+=H(n/pt,i+1)*(k*t+1)+(k*(t+1)+1);//id[n/pt]
    }
    return res;
}

void Work() {
    lim=sqrt(n)+0.5;
    id.cnt=0;
    for(int i=0;i<=lim;++i)
        id.id[0][i]=id.id[1][i]=0;
    psz=std::lower_bound(prime+1,prime+1+pcnt,lim)-prime;
    static long long kp[XN*2];int kpc=0;
    for(long long l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        g[id[kp[++kpc]=n/l]]=n/l-1;
    }
    for(int j=1;j<=psz;++j)
        for(int i=1;i<=kpc && (long long)prime[j]*prime[j]<=kp[i];++i)
            g[id[kp[i]]]-=g[id[kp[i]/prime[j]]]-(j-1);
    std::cout<<H(n,1)+1<<'\n';
}

int main() {
    Prep();
    int T;fin>>T;
    while(T--) {
        fin>>n>>k;
        Work();
    }
    return 0;
}

胡扯min_25筛

记质数集合为P

Problem

以低于线性的复杂度求一类积性函数的和,函数需满足
f(p),f(p^k)(p\in P)很好求。

Solution

{\rm Min}(x)x的最小质因子。

g

g(n,i)=\sum f(x)[ {x\le n {\ \rm and\ } (x \in P {\ \rm or\ } {\rm Min}(x)>P_i) } ]
n < P_i^2时,g(n,i)=g(n,i-1),因为最小的{\rm Min}(x)= P_i的合数就是P_i^2
否则,考虑从g(n,i-1)中减去{\rm Min}(x)= P_i的合数的贡献,也就是\sum f(x)[{x\le n,{\rm Min}(x)=P_i}],提出f(P_i),式子变成了f(P_i)(g([\frac n{P_i}],i-1)-\sum_{j=1}^{i-1}f(P_j))
\therefore
g(n,i)= \begin{cases} g(n,i-1)&,n < P_i^2\\ g(n,i-1)-f(P_i)(g([\frac n{P_i}],i-1)-\sum_{j=1}^{i-1}f(P_j))&,\text{otherwise} \end{cases}

h

h(n,i)=\sum_{j=1}^nf(j)[{\rm Min}(j)\ge P_i]
对于h(n,i)n以内质数的贡献是g(n,|P|)-\sum_{j=1}^{i-1}f(P_j)
现在考虑合数的贡献。
枚举合数的最小质因子P_j,那么最小质因子为P_j的合数的贡献就是\sum_{t\ge 1,P_j^{t+1}\le n}(h([\frac n{P_j^t}],j+1)f(P_j^t)+f(P_j^{t+1}))
因此
h(n,i)=g(n,|P|)-\sum_{j=1}^{i-1}f(P_j)+\sum_{j\ge i}\sum_{t\ge 1,P_j^{t+1}\le n}(h([\frac n{P_j^t}],j+1)f(P_j^t)+f(P_j^{t+1}))
答案就是h(n,1)了。

Tips

具体实现的时候,g不一定求\sum f(x)[ {x\le n {\ \rm and\ } (x \in P {\ \rm or\ } {\rm Min}(x)>P_i) } ],求出的g只要能使得上式可以被很快求出即可。

[hdu6223][Shenyang Onsite 2017]Infinite Fraction Path

Link

每个点有一个出边,每个点上有一个字符。求字典序最大的长度为n的路径。

Solution

广义后缀数组!预处理一下一个点往后跳2^i步到的点,然后魔改一下后缀数组就好了。

Code

const int XN=150000+11,XLOGN=18;

int par[XN][XLOGN];
char v[XN];
void Work() {
    int n;fin>>n>>(v+1);
    for(int i=1;i<=n;++i)
        par[i][0]=((long long)(i-1)*(i-1)+1)%n+1;
    for(int b=1;(1<<b)<=n;++b)
        for(int i=1;i<=n;++i)
            par[i][b]=par[par[i][b-1]][b-1];
    static int temp[2][XN],cnt[XN],*x=temp[0],*y=temp[1],sa[XN];
    int m=256;
    std::fill(cnt+1,cnt+1+m,0);
    for(int i=1;i<=n;++i)
        cnt[x[i]=v[i]]++;
    std::partial_sum(cnt+1,cnt+1+m,cnt+1);
    for(int i=1;i<=n;++i)
        sa[cnt[x[i]]--]=i;
    for(int b=0;(1<<b)<n;++b) {
        std::fill(cnt+1,cnt+1+m,0);
        for(int i=1;i<=n;++i)
            cnt[x[par[i][b]]]++;
        std::partial_sum(cnt+1,cnt+1+m,cnt+1);
        for(int i=1;i<=n;++i)
            y[cnt[x[par[i][b]]]--]=i;
        std::fill(cnt+1,cnt+1+m,0);
        for(int i=1;i<=n;++i)
            cnt[x[i]]++;
        std::partial_sum(cnt+1,cnt+1+m,cnt+1);
        for(int i=n;i>=1;--i)
            sa[cnt[x[y[i]]]--]=y[i];
        std::swap(x,y);
        int p=x[sa[1]]=1;
        for(int i=2;i<=n;++i)
            x[sa[i]]=y[sa[i-1]]==y[sa[i]] && y[par[sa[i-1]][b]]==y[par[sa[i]][b]]?p:++p;
        if((m=p)==n)
            break;
    }

    for(int i=1,p=sa[n];i<=n;++i,p=par[p][0])
        fout<<v[p];
    fout<<'\n';
}

int main() {
    int T;fin>>T;
    for(int ks=1;ks<=T;++ks) {
        printf("Case #%d: ",ks);
        Work();
    }
    return 0;
}

[cf1051F]The shortest statement

Link

n个点,m(m-n< 20)条边的连通图,多组询问两点的最短路。

Solution

求出图的任意一个生成树,把不在生成树上的边所连接的点都做一遍最短路。
如果两点最短路在生成树上,那么求出树上最短路就得到答案。
否则,考虑floyd算法,选取点u,用d[x][u]+d[u][y]更新答案。选取的点u只需选取那些不在生成树上的边连接的点,因为如果生成树上的某点到点x在树上的距离不是最短路,那么最短路一定经过某个不在生成树上的边连接的点,因此只用这些点松弛就行了。

[zoj4053][icpc 2018 Qingdao Online]Couleur

Link

Solution

假设将区间[l,r]中间的p ban掉,并且知道[l,r]的逆序堆数,容易计算出[l,p-1][p+1,r]之间,以及他们分别和p产生的逆序数。
用启发式的思想,每次将这两个区间中小的那个遍历,在另一个区间的主席树中查询,最终复杂度为O(n\log^2n)

Tips

赛场上没有过这道题,是因为偷懒直接把问题一般化,得到O(n\sqrt {n\log n})的做法,最终自闭。

Code

const int XN=1e5+11;

struct SegTree {
    struct Node {
        Node *son[2];
        long long cnt;

        Node():cnt(0) {
            son[0]=son[1]=null;
        }

        Node(void*):cnt(0) {
            son[0]=son[1]=this;
        }

        void *operator new(size_t flag) {
            static Node *pool=(Node*)malloc(XN*20*sizeof(Node)),*mem=pool++;
            return flag?mem++:mem=pool;
        }

        void Up() {
            cnt=son[0]->cnt+son[1]->cnt;
        }

    }*root[XN];
    int n,L,R;
    static Node *null;

    SegTree(int L,int R):n(0),L(L),R(R) {
        root[0]=null;
    }

    void Append(int v) {
        ++n;
        Modify(root[n]=root[n-1],L,R,v);
    }

    long long Query(int lbd,int rbd,int l,int r) {
        return lbd<=rbd && l<=r?Query(root[lbd-1],root[rbd],L,R,l,r):0;
    }

    static void Modify(Node *&pos,int L,int R,int goal) {
        pos=new Node(*pos);
        if(L==R)
            pos->cnt++;
        else {
            int M=(L+R)/2;
            if(goal<=M)
                Modify(pos->son[0],L,M,goal);
            else
                Modify(pos->son[1],M+1,R,goal);
            pos->Up();
        }
    }

    static long long Query(Node *pl,Node *pr,int L,int R,int l,int r) {
        if(L==l && R==r)
            return pr->cnt-pl->cnt;
        else {
            int M=(L+R)/2;
            if(r<=M)
                return Query(pl->son[0],pr->son[0],L,M,l,r);
            else if(M+1<=l)
                return Query(pl->son[1],pr->son[1],M+1,R,l,r);
            else
                return Query(pl->son[0],pr->son[0],L,M,l,M)
                      +Query(pl->son[1],pr->son[1],M+1,R,M+1,r);
        }
    }
}seg(-1,-1);

SegTree::Node *SegTree::null=new SegTree::Node((void*)0);

int a[XN];

void Work() {
    SegTree::Node::operator new(0);
    int n;fin>>n;
    new(&seg) SegTree(1,n);
    long long last=0;
    std::set<int> unav;
    std::multiset<long long> S;
    for(int i=1;i<=n;++i) {
        fin>>a[i];
        last+=seg.Query(1,i-1,a[i]+1,n);
        seg.Append(a[i]);
    }
    unav.insert(0);
    unav.insert(n+1);
    static std::multiset<long long>::iterator its[XN];
    its[1]=S.insert(last);
    for(int lop=1;lop<=n;++lop) {
        fout<<last<<(lop==n?'\n':' ');
        long long p;fin>>p;p^=last;
        std::set<int>::iterator pt=unav.upper_bound(p);
        int r=(*pt)-1,l=(*--pt)+1;
        unav.insert(p);
        long long tot=*its[l];
        S.erase(its[l]);
        tot-=seg.Query(l,p-1,a[p]+1,n)+seg.Query(p+1,r,1,a[p]-1);
        if(p-l>r-p) {
            long long ic=0,lc=0;
            for(int i=p+1;i<=r;++i) {
                lc+=seg.Query(l,p-1,a[i]+1,n);
                ic+=seg.Query(p+1,i-1,a[i]+1,n);
            }
            its[l]=S.insert(tot-lc-ic);
            its[p+1]=S.insert(ic);
        } else {
            long long ic=0,rc=0;
            for(int i=l;i<=p-1;++i) {
                rc+=seg.Query(p+1,r,1,a[i]-1);
                ic+=seg.Query(l,i-1,a[i]+1,n);
            }
            its[l]=S.insert(ic);
            its[p+1]=S.insert(tot-rc-ic);
        }
        last=*S.rbegin();
    }
}

int main() {
    int T;fin>>T;
    while(T--)
        Work();
    return 0;
}

[zoj4056][ICPC 2018 Qingdao Online]Press the Button

Link

Solution

问题最后是求一个长度为t的数轴上有多少对相邻的关键点距离\ge v,关键点是所有ac的倍数的点。
假设a < c
如果a < v则不存在
否则
考虑c的倍数与它后面相邻的a的倍数的距离
根据根据鸽巢原理,至多c步,这个距离就会开始循环
然后把t分成三部分:
不在循环+循环部分+一段没循环完的部分
细节巨多写到脑壳疼
~正解是模拟到{\rm lcm}(a,b),直接判断,然后循环~

Code

const int XN=1e6+11;
void Work() {
    long long v,t,a,b,c,d;
    fin>>a>>b>>c>>d>>v>>t;
    long long Ans=(long long)((t/a)+1)*b+(long long)((t/c)+1)*d-1;
    if(a>v && c>v) {
        if(a>c)
            std::swap(a,c);
        static int us[XN],ut;
        static long long crec[XN],prec[XN];
        ++ut;long long cnt=0,pc=0,py=0;
        for(;;pc+=c) {
            long long x=(pc/a)*a,y=x+a; 
            cnt+=((pc-x)>v)+((y-pc)>v)+(x-py)/a;
            py=y;
            if(us[y-pc]==ut)
                break;
            us[y-pc]=ut;
            crec[y-pc]=cnt;
            prec[y-pc]=py;
        }
        if(t>prec[py-pc]) {
            long long circnt=cnt-crec[py-pc],cirlen=py-prec[py-pc];
            Ans-=circnt*((t-prec[py-pc])/cirlen)+crec[py-pc];
            py+=((t-prec[py-pc])/cirlen-1)*cirlen;
            pc=(py/c+1)*c;
        } else {
            py=pc=0;
        }
        if(pc<=t)
            for(;;pc+=c) {
                long long x=(pc/a)*a,y=x+a;
                Ans-=((pc-x)>v)+(x-py)/a;
                if(y<=t) {
                    Ans-=((y-pc)>v);
                }
                py=y;
                if(pc+c>t)
                    break;
            }
        if(py<t)
            Ans-=(t-py)/a;
    }
    fout<<Ans<<'\n';
}

int main() {
    int T;fin>>T;
    while(T--)
        Work();
    return 0;
}