ZJT's Blog

大変に気分がいい~

[自选题 #104] [UR #10] 列队

2017-12-21 0:412018-9-13 23:51
集训队作业自选题抽象代数DP容斥群论

传送门

这道题成功让我认识到了我的抽代基本等于没学。。。

题目问的是给定群$G$到$S_n$之间有多少种单同态,其中$S_n$是$n$阶对称群。

先考虑怎么把单同态的限制去掉。根据抽代中的一些理论,$\phi:G_1\to G_2$为单同态当且仅当$\ker \phi=\{e\}$,即$G_1$中只有单位元能映射到$G_2$的单位元。所以,我们只要枚举$G_1$除了$\{e\}$以外的所有正规子群$H$,递归算出满足$\ker \phi=H$的$\phi$有多少种(这等价于求$G_1/H$到$G_2$的单同态数目)。这时如果我们能求出$G_1\to G_2$的同态数目,就能直接容斥得到单同态数目了。

考虑怎么计算$G\to S_n$的同态数目。我们知道,$G\to S_n$的一个同态,对应了群$G$在$\{1,\dots,n\}$上的一个作用。显然,这个作用的轨道之间是独立的,我们可以把每个轨道单独拿出来考虑。

记$\{1,\cdots,n\}$中的某个轨道为$R$,$\phi:G\to S_{|R|}$为同态,则根据抽代中的结论,$|R|$整除$|G|$,且对于任意$r\in R$,$r$的稳定子群$G_r$的陪集和整个轨道$R$一一对应。可以证明(证明过程参见wyy题解),在$1$的稳定子群$G_1$确定的情况下(这里$1$也可以换成别的元素),同态的数目为$(|R|-1)!$,所以总的同态数目为$(|R|-1)!\times($$G$的大小为$\frac{|H|}{|R|}$的子群数$)$。

我们先预处理出$G$的每种大小的子群分别有多少个,然后就能DP啦。设$f_i$表示$G\to S_i$的同态数,我们考虑决策包含$1$的轨道,记$c_k=(k-1)!\times($$G$的大小为$\frac{|H|}k$的子群数$)$,则:$$f_i=\sum_{j=1}^ic_jf_{i-j}\binom{i-1}{j-1}$$

这样我们就完成了对$G\to S_n$的同态数目的统计。找子群可以用广搜,枚举加入的元素,然后跑出包含已有元素的最小子群。最后,容斥的时候有一个小技巧,根据群的第三同构定理,$G$的商群的正规子群一定对应了一个$G$的一个正规子群。所以我们只要在$G$的正规子群里面进行容斥就行了,这样就不用算商群了。总复杂度大概是$O(n_N(n_G+n|G|))$,其中$n_N$为正规子群数,$n_G$为子群数。由于$n_N$和$n_G$都很小,就能直接跑过去啦。

#include <cstdio>
#include <cstring>
#include <set>
#include <queue>
#include <vector>
#define MAXN 40
#define MAXM 1010
#define LL long long
using namespace std;

const LL P=998244353;

int n,m;
int g[MAXN][MAXN],inv[MAXN],unit;
set<int> s_sg;
vector<int> sg,nsg;
LL c[MAXM][MAXM],fac[MAXM];
LL h[MAXM];

int extend(int s,int x){
    static queue<int> Q;
    s|=1<<(x-1);
    while(1){
        bool flag=0;
        for(int i=1;i<=n;i++)
            if(s&(1<<(i-1)))
                for(int j=1;j<=n;j++)
                    if(s&(1<<(j-1)))
                        if((~s)&(1<<(g[i][j]-1))){
                            flag=1;
                            s|=1<<(g[i][j]-1);
                        }
        if(!flag) return s;
    }
}

bool check_normal(int s){
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            if((s&(1<<(j-1))) && ((~s)&(1<<(g[g[i][j]][inv[i]]-1))))
                return 0;
    return 1;
}

void bfs(){
    static queue<int> Q;
    for(int i=1;i<=n;i++){
        int x=extend(0,i);
        if(!s_sg.count(x)){
            s_sg.insert(x);
            Q.push(x);
        }
    }
    while(!Q.empty()){
        int x=Q.front();
        Q.pop();
        for(int i=1;i<=n;i++){
            int y=extend(x,i);
            if(!s_sg.count(y)){
                s_sg.insert(y);
                Q.push(y);
            }
        }
    }
    for(set<int>::iterator it=s_sg.begin();it!=s_sg.end();it++){
        if(check_normal(*it)) nsg.push_back(*it);
        sg.push_back(*it);
    }
}

void dp(int x){
    static int cnt[MAXN];
    memset(cnt,0,sizeof cnt);
    int s=nsg[x],sz=0;
    for(int l=1<<30;l;l>>=1) if(s&l) sz++;
    sz=n/sz;
    for(int i=0;i<sg.size();i++)
        if((sg[i]&s)==s){
            int temp=0;
            for(int l=1<<30;l;l>>=1) if(sg[i]&l) temp++;
            cnt[temp/(n/sz)]++;
        }
    static LL f[MAXM];
    f[0]=1;
    for(int i=1;i<=m;i++){
        f[i]=0;
        for(int j=1;j<=i && j<=sz;j++)
            if(sz%j==0)
                f[i]=(f[i]+c[i-1][j-1]*f[i-j]%P*cnt[sz/j]%P*fac[j-1])%P;
    }
    h[x]=f[m];
    for(int i=x+1;i<nsg.size();i++)
        if((nsg[i]&s)==s) h[x]=(h[x]-h[i])%P;
}

void solve(){
    s_sg.clear();
    sg.clear();
    nsg.clear();
    scanf("%d%d",&m,&n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            scanf("%d",&g[i][j]);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            for(int k=1;k<=n;k++)
                if(g[i][g[j][k]]!=g[g[i][j]][k]){
                    puts("0");
                    return;
                }
    for(int i=1;i<=n;i++)
        if(g[i][i]==i) unit=i;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            if(g[i][j]==unit)
                inv[i]=j;
    bfs();
    for(int i=nsg.size()-1;i>=0;i--) 
        dp(i);
    printf("%lld\n",(h[0]+P)%P);
}

void init(){
    for(int i=0;i<MAXM;i++){
        c[i][0]=1;
        for(int j=1;j<=i;j++)
            c[i][j]=(c[i-1][j-1]+c[i-1][j])%P;
    }
    fac[0]=1;
    for(int i=1;i<MAXM;i++) fac[i]=fac[i-1]*i%P;
}

int main(){
#ifdef DEBUG
    freopen("104.in","r",stdin);
#endif
    init();
    int T;
    scanf("%d",&T);
    while(T--) solve();
}

查看详细内容

集训队作业之AGC005-D

2017-12-11 17:492017-12-11 17:49
集训队作业AtCoder容斥DP

传送门

先考虑容斥,对于每个$i$处理出有$i$个位置不满足条件的方案数。发现要处理模$2k$同余的点之间的冲突关系,直接对于每个同余类dp出方案数(每个位置要么$|a_i-i|\ne k$,要么$a_i=i+k$,要么$a_i=i-k$,相邻两个位置不能前一个$a_i=i+k$,后一个$a_i=i-k$)。我们把这个DP的结果表示成多项式,那整体的多项式就是每个同余类的多项式卷积起来。由于$n$只有$2000$,直接暴力卷积就行了。(其实可以做到$O(n\log n)$?)

#include <cstdio>
#include <cstring>
#define MAXN 2010
#define LL long long

const LL P=924844033;

struct Poly{
    LL a[MAXN];
    int len;
}g[MAXN][2],res;

int n,k;
LL fac[MAXN];

void add(Poly &x,Poly &y){
    if(x.len>y.len) y.len=x.len;
    for(int i=0;i<=x.len;i++) y.a[i]=(y.a[i]+x.a[i])%P;
}

void add2(Poly &x,Poly &y){
    if(x.len+1>y.len) y.len=x.len+1;
    for(int i=0;i<=x.len;i++) y.a[i+1]=(y.a[i+1]+x.a[i])%P;
}

void mul(Poly &x,Poly &y){
    static LL temp[MAXN];
    for(int i=0;i<=x.len+y.len;i++) temp[i]=0;
    for(int i=0;i<=x.len;i++)
        for(int j=0;j<=y.len;j++)
            temp[i+j]=(temp[i+j]+x.a[i]*y.a[j])%P;
    y.len+=x.len;
    for(int i=0;i<=y.len;i++)
        y.a[i]=temp[i];
}

void gao(int x){
    int p1=x,p2=x+(n-x)/(2*k)*(2*k);
    bool flag1=(p1-k>=1),flag2=(p2+k<=n);
    g[p2][1].len=1;
    g[p2][1].a[1]=1;
    if(flag2){
        g[p2][0].len=1;
        g[p2][0].a[0]=g[p2][0].a[1]=1;
    }else{
        g[p2][0].len=0;
        g[p2][0].a[0]=1;
    }
    for(int i=p2-2*k;i>=p1;i-=2*k){
        int i2=i+2*k;
        add(g[i2][0],g[i][0]);
        add2(g[i2][0],g[i][0]);
        add2(g[i2][0],g[i][1]);
        add(g[i2][1],g[i][0]);
        add2(g[i2][1],g[i][1]);
    }
    if(flag1) add(g[p1][1],g[p1][0]);
}

int main(){
#ifdef DEBUG
    freopen("D.in","r",stdin);
#endif
    fac[0]=1;
    for(int i=1;i<MAXN;i++) fac[i]=fac[i-1]*i%P;
    scanf("%d%d",&n,&k);
    res.a[0]=1;
    for(int i=1;i<=2*k && i<=n;i++){
        gao(i);
        mul(g[i][0],res);
    }
    LL ans=0;
    for(int i=0;i<=n;i++){
        LL delta=(i&1)?(-1):1;
        ans=(ans+delta*res.a[i]*fac[n-i])%P;
    }
    ans=(ans+P)%P;
    printf("%lld\n",ans);
}

查看详细内容

集训队作业之AGC018-E

2017-11-23 18:552017-11-23 18:55
集训队作业AtCoder容斥组合数学FFT

传送门

看了题就先瞎推了一下,不知道怎么就推成了一个卷积的形式,然后一看模数$10^9+7$,范围是$10^6$,感觉AtCoder的机子那么牛逼应该可以强行任意模数FFT跑过去。然后就写了一个胡乱卷积,结果爆精度了。。。double改成long double之后又T了。。。然后尝试各种卡常+卡精度技巧,结果一点效果都没有。然后又仔细想了一下,发现这题随便$O(n)$。。。心情复杂.jpg

首先把起点、终点的限制容斥一下,问题转化成了在一个大框里选两个点作为起点和终点,在里面的一个小框里面选一个中间点,要求起点->中间点->终点路径数的和。对于某个中间点,设它到大框左端的距离是$d_x$,到顶端的距离是$d_y$,则在它左上方选一个起点走到它的路径总数就是:

$$ \begin{aligned} \sum_{x=0}^{d_x}\sum_{y=0}^{d_y}\binom{x+y}{x}&=\sum_{x=0}^{d_x}\binom{x+d_y+1}{x+1}\\ &=\sum_{x=0}^{d_x+1}\binom{x+d_y+1}{x+1}-1\\ &=\binom{d_x+d_y+2}{d_x+1}-1 \end{aligned} $$

后面那个$-1$在容斥的时候会被抵消掉,我们可以不用管他。

然而直接枚举中间点是$O(n^2)$的,为了减少这个复杂度,我想到了两个做法,一个是一开始那个爆精度的傻逼卷积做法(过不了),一个是比较简单的线性做法(能过)。

先讲一下那个过不了的弱智做法,就是考虑每条起点到终点的路径,如果穿过了中间那个小框,那么这条路径的贡献就是路径跟小框的交集的点数(这里面的每个点都能作为中间点产生$1$的贡献)。穿过小框的路径一定是从左边或者上方进入小框,从右边或者下方离开小框。我们根据入口、出口的位置,分左右、左下、上右、上下四种情况考虑,发现其实是一个卷积。然而这题模数是$10^9+7$,用按$\sqrt P$分段的那种任意模数FFT的话,double会爆精度,long double会T。。。这还是在AtCoder的机器上,在我的机子上double都要跑10s。。。

线性的做法也是类似,考虑每条路径穿过中间那个框的点数,但是把中间那个框容斥了一下,小框变成了大框左上角的一部分,这样就只用考虑出口了。由于随便选起点的路径总数能表示成一个组合数$\binom{d_x+d_y+2}{d_x+1}$,所以可以直接把起点看成固定点,终点也是一样。起点固定的情况下,起点->出口经过的点数也是确定的,所以直接枚举出口统计一下就行了,复杂度是$O(n)$(虽然因为容斥带了一个64的常数)。

代码($O(n)$,能过):

#include <cstdio>
#include <cstring>
#include <cassert>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#define MAXN 2100000
#define LL long long
#define y1 zjtsb_y1
using namespace std;

const LL P=1000000007;

LL fac[MAXN],invfac[MAXN];

LL getPow(LL x,LL y){
    LL res=1;
    while(y){
        if(y&1) res=res*x%P;
        x=x*x%P;
        y>>=1;
    }
    return res;
}

LL getC(int x,int y){
    if(x<y) return 0;
    return fac[x]*invfac[y]%P*invfac[x-y]%P;
}

void init(){
    fac[0]=1;
    for(int i=1;i<MAXN;i++) fac[i]=fac[i-1]*i%P;
    invfac[MAXN-1]=getPow(fac[MAXN-1],P-2);
    for(int i=MAXN-2;i>=0;i--) invfac[i]=invfac[i+1]*(i+1)%P;
}

LL calc(int x1,int x2,int y1,int y2){
    if(!x1 || !x2 || !y1 || !y2) return 0;
    LL res=0;
    for(int i=1;i<=x1;i++)
        res=(res+getC(y1+i-2,i-1)*getC(y2+x1+x2-i-1,y2-1)%P*(y1+i-1))%P;
    for(int i=1;i<=y1;i++)
        res=(res+getC(x1+i-2,i-1)*getC(x2+y1+y2-i-1,x2-1)%P*(x1+i-1))%P;
    return res;
}

int x1,x2,x3,x4,x5,x6;
int y1,y2,y3,y4,y5,y6;

LL gao(int sx,int sy,int ex,int ey){
    LL res=0;
    res+=calc(x4-sx+1,ex-x4,y4-sy+1,ey-y4);
    res-=calc(x3-sx,ex-x3+1,y4-sy+1,ey-y4);
    res-=calc(x4-sx+1,ex-x4,y3-sy,ey-y3+1);
    res+=calc(x3-sx,ex-x3+1,y3-sy,ey-y3+1);
    return (res%P+P)%P;
}

int main(){
#ifdef DEBUG
    freopen("E.in","r",stdin);
#endif
    scanf("%d%d%d%d%d%d",&x1,&x2,&x3,&x4,&x5,&x6);
    scanf("%d%d%d%d%d%d",&y1,&y2,&y3,&y4,&y5,&y6);
    init();
    LL ans=0;
    for(int i=0;i<16;i++){
        int c=0;
        int sx=x1-1,sy=y1-1,ex=x6+1,ey=y6+1;
        if(i&1){
            c++;
            sx=x2;
        }
        if(i&2){
            c++;
            sy=y2;
        }
        if(i&4){
            c++;
            ex=x5;
        }
        if(i&8){
            c++;
            ey=y5;
        }
        if(c&1) ans-=gao(sx,sy,ex,ey);
        else ans+=gao(sx,sy,ex,ey);
    }
    ans=(ans%P+P)%P;
    printf("%lld\n",ans);
    return 0;
}

代码($O(n\log n)$,爆精度):

#include <cstdio>
#include <cstring>
#include <cassert>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#define MAXN 2100000
#define LL long long
#define y1 zjtsb_y1
using namespace std;

const int MAXW=2097152;
const LL P=1000000007;
const long double PI=acos(-1.0);

namespace FFT{
    struct cplx{
        long double r,i;
        cplx(long double _r=0,long double _i=0):r(_r),i(_i){}
        friend cplx operator+(cplx x,cplx y){ return cplx(x.r+y.r,x.i+y.i); }
        friend cplx operator-(cplx x,cplx y){ return cplx(x.r-y.r,x.i-y.i); }
        friend cplx operator*(cplx x,cplx y){ return cplx(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r); }
    }wn[MAXW];

    void init(){
        for(int i=0;i<MAXW;i++) wn[i]=cplx(cos(2*PI/MAXW*i),sin(2*PI/MAXW*i));
    }

    void fft(cplx *a,int len,int flag){
        static int rev[MAXN],revlen;
        if(revlen!=len){
            revlen=len;
            for(int i=1;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?(len>>1):0);
        }
        for(int i=0;i<len;i++)
            if(i<rev[i])
                swap(a[i],a[rev[i]]);
        for(int l=2;l<=len;l<<=1){
            int l2=l>>1;
            for(int i=0;i<len;i+=l)
                for(int j=0;j<l2;j++){
                    cplx t1=a[i+j],t2=a[i+j+l2]*wn[MAXW/l*j];
                    a[i+j]=t1+t2;
                    a[i+j+l2]=t1-t2;
                }
        }
        if(flag==-1){
            for(int i=0;i<len;i++) a[i].r/=len;
            for(int i=1;i<len;i++)
                if(i<len-i) swap(a[i],a[len-i]);
        }
    }
}

LL fac[MAXN],invfac[MAXN];
int n,m;
int x1,x2,x3,x4,x5,x6;
int y1,y2,y3,y4,y5,y6;
LL f1[MAXN],f2[MAXN],f3[MAXN],f4[MAXN];

LL getPow(LL x,LL y){
    LL res=1;
    while(y){
        if(y&1) res=res*x%P;
        x=x*x%P;
        y>>=1;
    }
    return res;
}

LL getC(int x,int y){
    if(x<y) return 0;
    return fac[x]*invfac[y]%P*invfac[x-y]%P;
}

LL calcG(int x,int y){
    return getC(x+y+2,x+1)-1;
}

void init(){
    fac[0]=1;
    for(int i=1;i<MAXN;i++) fac[i]=fac[i-1]*i%P;
    invfac[MAXN-1]=getPow(fac[MAXN-1],P-2);
    for(int i=MAXN-2;i>=0;i--) invfac[i]=invfac[i+1]*(i+1)%P;
    FFT::init();
}

void gaoF(){
    for(int i=1;i<=m;i++){
        int x=x3-1,y=y3-1+i;
        f1[i]=((calcG(x-x1,y-y1)-calcG(x-x2-1,y-y1)-calcG(x-x1,y-y2-1)+calcG(x-x2-1,y-y2-1))%P+P)%P;
    }
    for(int i=1;i<=n;i++){
        int x=x3-1+i,y=y3-1;
        f2[i]=((calcG(x-x1,y-y1)-calcG(x-x2-1,y-y1)-calcG(x-x1,y-y2-1)+calcG(x-x2-1,y-y2-1))%P+P)%P;
    }
    for(int i=1;i<=m;i++){
        int x=x4+1,y=y3-1+i;
        f3[i]=((calcG(x6-x,y6-y)-calcG(x5-x-1,y6-y)-calcG(x6-x,y5-y-1)+calcG(x5-x-1,y5-y-1))%P+P)%P;
    }
    for(int i=1;i<=n;i++){
        int x=x3-1+i,y=y4+1;
        f4[i]=((calcG(x6-x,y6-y)-calcG(x5-x-1,y6-y)-calcG(x6-x,y5-y-1)+calcG(x5-x-1,y5-y-1))%P+P)%P;
    }
}

void mul(LL *a,LL *b,LL *c,int l1,int l2){
    using namespace FFT;
    const LL M=32000;
    static cplx t1[MAXN],t2[MAXN];
    static LL s1[MAXN],s2[MAXN],s3[MAXN];
    int sizew;
    for(sizew=1;sizew<=l1+l2;sizew<<=1);
    for(int i=0;i<sizew;i++) t1[i]=t2[i]=cplx();
    for(int i=0;i<=l1;i++) t1[i].r=(a[i]/M)+(a[i]%M);
    for(int i=0;i<=l2;i++) t2[i].r=(b[i]/M)+(b[i]%M);
    fft(t1,sizew,1); fft(t2,sizew,1);
    for(int i=0;i<sizew;i++) t1[i]=t1[i]*t2[i];
    fft(t1,sizew,-1);
    for(int i=0;i<sizew;i++) s1[i]=t1[i].r+0.5;

    for(int i=0;i<sizew;i++) t1[i]=t2[i]=cplx();
    for(int i=0;i<=l1;i++) t1[i].r=a[i]/M;
    for(int i=0;i<=l2;i++) t2[i].r=b[i]/M;
    fft(t1,sizew,1); fft(t2,sizew,1);
    for(int i=0;i<sizew;i++) t1[i]=t1[i]*t2[i];
    fft(t1,sizew,-1);
    for(int i=0;i<sizew;i++) s2[i]=t1[i].r+0.5;

    for(int i=0;i<sizew;i++) t1[i]=t2[i]=cplx();
    for(int i=0;i<=l1;i++) t1[i].r=a[i]%M;
    for(int i=0;i<=l2;i++) t2[i].r=b[i]%M;
    fft(t1,sizew,1); fft(t2,sizew,1);
    for(int i=0;i<sizew;i++) t1[i]=t1[i]*t2[i];
    fft(t1,sizew,-1);
    for(int i=0;i<sizew;i++) s3[i]=t1[i].r+0.5;

    for(int i=0;i<sizew;i++){
        s1[i]%=P;
        s2[i]%=P;
        s3[i]%=P;
    }
    for(int i=0;i<sizew;i++) s1[i]-=s2[i]+s3[i];
    for(int i=0;i<sizew;i++) c[i]=(s3[i]+M*s1[i]+M*M%P*s2[i])%P;
}

LL gao(){
    static LL t1[MAXN],t2[MAXN],t3[MAXN];
    LL res=0;
    //f1*f3
    for(int i=1;i<=m;i++) t1[i]=f1[i],t2[i]=f3[m-i+1];
    mul(t1,t2,t3,m,m);
    for(int i=0;i<=m-1;i++){
        LL t=t3[m+1-i];
        res=(res+t*(n+i)%P*getC(n-1+i,i))%P;
    }
    //f2*f4
    for(int i=1;i<=n;i++) t1[i]=f2[i],t2[i]=f4[n-i+1];
    mul(t1,t2,t3,n,n);
    for(int i=0;i<=n-1;i++){
        LL t=t3[n+1-i];
        res=(res+t*(m+i)%P*getC(m-1+i,i))%P;
    }
    //f1*f4
    for(int i=0;i<m;i++) t1[i]=f1[m-i]*invfac[i]%P;
    for(int i=0;i<n;i++) t2[i]=f4[i+1]*invfac[i]%P;
    mul(t1,t2,t3,m-1,n-1);
    for(int i=0;i<=n+m-2;i++){
        LL t=t3[i];
        res=(res+t*fac[i]%P*(i+1))%P;
    }
    //f2*f3
    for(int i=0;i<n;i++) t1[i]=f2[n-i]*invfac[i]%P;
    for(int i=0;i<m;i++) t2[i]=f3[i+1]*invfac[i]%P;
    mul(t1,t2,t3,n-1,m-1);
    for(int i=0;i<=n+m-2;i++){
        LL t=t3[i];
        res=(res+t*fac[i]%P*(i+1))%P;
    }
    return res;
}

int main(){
#ifdef DEBUG
    freopen("E.in","r",stdin);
#endif
    scanf("%d%d%d%d%d%d",&x1,&x2,&x3,&x4,&x5,&x6);
    scanf("%d%d%d%d%d%d",&y1,&y2,&y3,&y4,&y5,&y6);
    n=x4-x3+1;
    m=y4-y3+1;
    init();
    gaoF();
    LL ans=gao();
    printf("%lld\n",(ans%P+P)%P);
    return 0;
}

查看详细内容

集训队作业之AGC013-E

2017-11-2 17:72018-1-2 17:21
集训队作业AtCoderDP容斥线性递推

传送门

感觉这题非常不错啊,中间的思路挺妙的(虽然做的时候卡了一会儿)。

先考虑没有标记的情况,我们来搞一个dp,设$f_n$为总长为$n$的所有方案的美丽度之和,则我们有转移:$f_n=\sum k^2f_{n-k}$。考虑用矩乘算这个东西,我们记$$s_2=\sum_{k=1}^n f_k(n-k)^2,s_1=\sum_{k=1}^n f_k(n-k),s_0=\sum_{k=1}^n f_k$$,则每次$n:=n+1$时,所有$(n-k)^2$会变成$(n-k)^2+2(n-k)+1$,这时拿$s_1,s_0$去更新一下$s_2$就行了,同理也可以用$s_0$去更新$s_1$。新的$f_{n+1}$也可以从这几个东西里面得到,直接加到$s_0$里面去就行了。这个东西是个常系数线性递推,所以可以直接矩乘,这样我们就可以在$O(\log n)$的时间里面算出没有标记的一段的贡献了。

考虑有标记的情况,我们还是dp,设$f_i$表示完全覆盖左端点到第$i$个标记这一段的所有方案的美丽度之和(这里第$i$个之前的标记不能作为正方形的端点,但第$i$个标记必须作为最右边的正方形的右端点)。我们考虑容斥,用$g(len)$表示上面矩乘算出来的长度为$len$的答案,则$f_i=g(X_i)-\sum_{j<i}f_jg(X_i-X_j)$,这里$j$是枚举第一个不合法的标记的位置。直接这样容斥dp,维护一下矩阵的前缀乘积和前缀乘积的逆,就可以做到$O(M^2)$的复杂度。

为了优化这个做法,我们考虑每个$f_j$对后面所有$f_i$的贡献。当 $i:=i+1$ 时,贡献由$g(X_i-X_j)$变成了$g(X_{i+1}-X_j)$。由于$g(x)$可以表示成某个矩阵$T^x$中的某个位置的值,所以$g(X_{i+1}-X_j)$可以由$T^{X_i-X_j}$乘上$T^{X_{i+1}-X_i}$得到。由于每次 $i:=i+1$ ,所有$j$乘上的矩阵都是$T^{X_{i+1}-X_i}$,再根据矩阵乘法的分配律,我们可以直接带着矩阵来做,每次转移的时候都乘上$T^{X_{i+1}-X_i}$,就可以在$O(m\log n)$的时间里完成dp了。

#include <cstdio>
#include <cstring>
#define MAXN 100010
#define LL long long

const LL P=1000000007;

struct Mat{
    struct Data{
        LL a[3][3];
    }*data;

    Mat(){
        data=new Data;
        memset(data->a,0,sizeof data->a);
    }

    void release(){ delete data; }

    friend void mul(Mat &x,Mat y){
        static Mat res;
        for(int i=0;i<3;i++) 
            for(int j=0;j<3;j++){
                res.data->a[i][j]=0;
                for(int k=0;k<3;k++)
                    res.data->a[i][j]+=x.data->a[i][k]*y.data->a[k][j];
                res.data->a[i][j]%=P;
            }
        for(int i=0;i<3;i++)
            for(int j=0;j<3;j++)
                x.data->a[i][j]=res.data->a[i][j];
    }
}trans;

Mat getPow(Mat _x,int y){
    Mat res;
    for(int i=0;i<3;i++) res.data->a[i][i]=1;
    static Mat x;
    for(int i=0;i<3;i++) for(int j=0;j<3;j++) x.data->a[i][j]=_x.data->a[i][j];
    while(y){
        if(y&1) mul(res,x);
        mul(x,x);
        y>>=1;
    }
    return res;
}

void init(){
    trans.data->a[0][0]=1; trans.data->a[0][1]=1; trans.data->a[0][2]=1;
    trans.data->a[1][0]=0; trans.data->a[1][1]=1; trans.data->a[1][2]=2;
    trans.data->a[2][0]=1; trans.data->a[2][1]=1; trans.data->a[2][2]=2;
}

LL getG(int x){
    Mat temp=getPow(trans,x);
    return temp.data->a[0][2];
}

int n,m;
int p[MAXN];
LL f[MAXN];

void gao(){
    Mat temp;
    f[1]=getG(p[1]);
    temp.data->a[0][0]=f[1];
    for(int i=2;i<=m;i++){
        mul(temp,getPow(trans,p[i]-p[i-1]));
        f[i]=(getG(p[i])-temp.data->a[0][2]+P)%P;
        temp.data->a[0][0]=(temp.data->a[0][0]+f[i])%P;
    }
}

int main(){
#ifdef DEBUG
    freopen("E.in","r",stdin);
#endif
    scanf("%d%d",&n,&m);
    for(int i=1;i<=m;i++) scanf("%d",p+i);
    p[++m]=n;
    init();
    gao();
    printf("%lld\n",f[m]);
    return 0;
}

查看详细内容

集训队作业之ARC058-D

2017-10-25 14:252017-10-25 14:25
集训队作业AtCoder容斥组合数学

传送门

考虑容斥掉每个不符合条件的方案,枚举进入非法区域时的纵坐标,然后几个组合数随便乘一乘减一减就行了。

#include <cstdio>
#include <cstring>
#define MAXN 1000010
#define LL long long

const LL P=1000000007;
int n,m,a,b;
LL fac[MAXN],invfac[MAXN];

LL getPow(LL x,LL y){
    LL res=1;
    while(y){
        if(y&1) res=res*x%P;
        x=x*x%P;
        y>>=1;
    }
    return res;
}

void init(){
    fac[0]=1;
    for(int i=1;i<MAXN;i++) fac[i]=fac[i-1]*i%P;
    invfac[MAXN-1]=getPow(fac[MAXN-1],P-2);
    for(int i=MAXN-2;i>=0;i--) invfac[i]=invfac[i+1]*(i+1)%P;
}

LL getC(int x,int y){
    return fac[x]*invfac[y]%P*invfac[x-y]%P;
}

LL g(int x,int y){ return getC(x+y-2,y-1); }

int main(){
    init();
    scanf("%d%d%d%d",&n,&m,&a,&b);
    LL ans=g(n,m);
    for(int i=1;i<=b;i++)
        ans=(ans-g(n-a,i)*g(a,m-i+1)%P+P)%P;
    printf("%lld\n",ans);
}

查看详细内容