HDU 6036 Division Game

HDU 6036 Division Game

http://acm.hdu.edu.cn/showproblem.php?pid=6036

题目中需要的用到知识清单:

1: FFT mod P 意义下的实现

2:容斥原理与莫比乌斯反演

3:简单的组和数计算。

注:2,3可以参阅组合数学第2,5,6章节

<script type="math/tex; mode=display" id="MathJax-Element-3"></script>

————————————————————————————

注:这是题解思路

正题:

对于一堆石子。取之前。要是取之后的倍数。
相当于。将 n 质因数分解后,取走若干素因子。直到取完游戏结束。

令 , 取x次,恰好取完的方案数为: C(x)

那么取 x1 次,取不完的方案数也为: C(x)

这是因为,两种方案不同情况,一一对应。

<script type="math/tex; mode=display" id="MathJax-Element-71"></script>

原题堆编号是从 1 开始对 。这里变成 0
那么对于 i[0,k1]
t 次,在第i堆石子结束游戏的情况数:
首先, [0,i1] 这部分,被取了 x 次,都没取完,这部分的方案贡献为:

C(x+1)i

其次, [i+1,k1] 这部分,被取了 x1 次,也都没取完,贡献为:

C(x)ki1

而,第 i 堆取了x次,取完了,这部分贡献为:

C(x)

综上。取 x 次,在第i堆结束游戏的情况数量为:

C(x+1)iC(x)ki

记, ans[i] 是在第 i 堆结束的情况数量

记,w=mi=1ei

则:

ans[i]=t=1wC(t+1)iC(t)ki

下面计算 C(x)
对于 n 的质因数分解

<P1,P2,P3,...Pm>

对应指数为:
<e1,e2,e3,...em>
<script type="math/tex; mode=display" id="MathJax-Element-33"> </script>
对于 C(x) , x[1,w]
Pi 将在 x 次取石子后被取完。

rt 是第 t 次取走Pi的数量

我们有:
r1+r2+r3+...+rx=ei
对于这个方程的非负整数解的数量为:

(ei+x1x1)

注:这个公式可以通过对 x1 个 分隔符 ‘ | ’,ei 的排列数量得到,分隔符把区间分成了x块,每块区间中 的数量对应一个r,排列的数量等同于分隔符选位置的方法数。
记:

H(x)=i=1m(ei+x1x1)

首先,特别说明, H(x)C(x)
因为,原题要求,每次至少取走一个素因子,而 H(x) 中,可能会有某一次,一个也没有取的情况。
x 次取数中,至少有y次,
一个素因子也没有取的方法数为 D(x,y)
则:

D(x,y)=(xy)H(xy)=(xxy)H(xy)

应用下式,容斥原理:
对于集合 Xn{1,2,3,...,n} ,
对于集合S的子集, A1,A2,A3,..An
它们对立相反事件为: A1,A2,A3,...An
那么S中不属于 A1,A2,A3,...An 的元素数量为:

A1A2A3..An=L Xn(1)n|L|iLAi

应用到上式:
对于所有位置都至少取一个素因子的方案数量,等于至少有y个位置不能取,不同y之间的容斥。

C(x)=y=0x(1)xyD(x,y)

C(x)=y=0x(1)xy(xxy)H(xy)

C(x)=y=0x(1)xy(xy)H(y)

C(x)=y=0x(1)xyx!y!(xy)!H(y)

所以:

C(x)x!=y=0x(1)xy(xy)!H(y)y!

令多项式 A(x) 有:

A(x)=(1)00!x0+(1)11!x1+...+(1)ww!xw

令多项式 B(x) 有:

B(x)=H(0)0!x0+H(1)1!x1+...+H(w)w!xw

那么有:

A(x)B(x)=C(0)0!x0+C(1)1!x1+...+C(w)w!x2w

利用FFT计算出这个卷积。
我们也就得到了 C[]
注:由于代码写的有点搓,提交请使用G++
下面是代码:
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <math.h>
#define MAXN 100050
using namespace std;

typedef long long LL;
const LL P=985661441;
const LL G=3;

LL Pow(LL a,int b)
{
    LL tmp=1;
    while(b)
    {
        if(b&1)
            tmp=tmp*a%P;
        a=a*a%P;
        b>>=1;
    }
    return tmp;
}


void change(LL y[],int n)
{
    int b=n>>1,s=n-1;
    for(int i=1,j=n>>1;i<s;i++)
    {
        if(i<j)swap(y[i],y[j]);
        int k=b;
        while(j>=k)
        {
            j-=k;
            k>>=1;
        }
        j+=k;
    }
}

void NTT(LL y[],int len,int on)
{
    change(y,len);
    for(int h=2;h<=len;h<<=1)
    {
        LL wh=Pow(G,(P-1)/h);
        if(on<0)wh=Pow(wh,P-2);
        for(int i=0;i<len;i+=h)
        {
            LL w=1;
            int r=h>>1;
            for(int k=i,s=r+i;k<s;k++)
            {
                LL u=y[k];
                LL t=w*y[k+r]%P;
                y[k]=u+t;   if(y[k]>=P)y[k]-=P;
                y[k+r]=u-t; if(y[k+r]<0) y[k+r]+=P;
                w=w*wh%P;
            }
        }
    }
    if(on<0)
    {
        LL I=Pow((LL)len,P-2);
        for(int i=0;i<len;i++) y[i]=y[i]*I%P;
    }
}

LL f[MAXN*5];
LL E[20];
LL p[20];
LL B[MAXN*5];
LL A[MAXN*5];
LL F[MAXN*5];
LL IV[MAXN*5];
LL V[MAXN*5];
LL iv[MAXN];
LL ans[MAXN];
LL x;
int main ()
{
    int ca=1;
    F[0]=1;
    for(int i=1,size=5*MAXN;i<size;i++) F[i]=F[i-1]*i%P;
    for(int i=1,size=5*MAXN;i<size;i++) IV[i]=Pow(F[i],P-2);
    for(int i=1,size=5*MAXN;i<size;i++) V[i]=Pow((LL)i,P-2);
    int m,kk;
    while(scanf("%d %d",&m,&kk)==2)
    {
        x=0;
        for(int i=0;i<m;i++) scanf("%lld %lld",p+i,E+i);
        for(int i=0;i<m;i++) x+=E[i];
        B[1]=1;
        B[0]=0;
        x++;
        int len=1;
        while(len<2*x)len<<=1;

        for(int k=1;k<x;k++)
        {
            B[k+1]=B[k];
            for(int i=0;i<m;i++)
                B[k+1]=B[k+1]*(E[i]+k)%P*V[k]%P;
        }

        for(int k=1;k<x;k++)B[k]=B[k]*IV[k]%P;

        for(int k=(int)x;k<len;k++) B[k]=0;

        for(int k=0,on=1;k<x;k++,on=-on) A[k]=(IV[k]*(LL)on+P)%P;

        for(int k=(int)x;k<len;k++) A[k]=0;

        A[0]=1;
        NTT(B,len,1);
        NTT(A,len,1);

        for(int k=0;k<len;k++)  f[k]=(A[k]*B[k])%P;

        NTT(f,len,-1);

        for(int k=0;k<len;k++)  f[k]=f[k]*F[k]%P;
        f[x]=0;
        for(int i=1;i<x;i++) iv[i]=Pow(f[i],P-2);
        for(int k=0;k<(int)x;k++)   ans[k]=Pow(f[k],kk)%P;
        printf("Case #%d:",ca++);
        for(int i=0;i<kk;i++)
        {
            LL a=0;
            for(int k=1;k<(int)x;k++)
            {
                a+=ans[k];
                if(a>=P)a-=P;
                ans[k]=ans[k]*f[k+1]%P*iv[k]%P;
            }
            printf(" %lld",a);
        }
        printf("\n");
    }
    return 0;
}
/*
1 1
2 2
2 1
3 1
5 1
1 2
2 3
2 2
2 4
5 4
 */
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值