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)$

那么取$x-1$次，取不完的方案数也为：$C(x)$

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

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

原题堆编号是从$1$开始对 。这里变成 $0$

那么对于$i\in[0,k-1]$

取$t$次，在第$i$堆石子结束游戏的情况数：

首先，$[0,i-1]$这部分，被取了$x$次，都没取完，这部分的方案贡献为：

$$C(x+1)^i$$

其次，$[i+1,k-1]$这部分，被取了$x-1$次，也都没取完，贡献为：

$$C(x)^{k-i-1}$$

而，第$i$堆取了$x$次，取完了，这部分贡献为：

$$C(x)$$

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

$$C(x+1)^iC(x)^{k-i}$$

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

记，$w=\sum_{i=1}^me^i$

则：

$$ans[i]=\sum_{t=1}^{w}C(t+1)^iC(t)^{k-i}$$

下面计算$C(x)$

对于$n$的质因数分解

$$<P_1,P_2,P_3,...P_m>$$

对应指数为：

<e1,e2,e3,...em>

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

对于$C(x)$ , $x\in[1,w]$

$P_i$将在$x$次取石子后被取完。

设 $r_t$ 是第 $t$ 次取走$P_i$的数量

我们有： $$r_1+r_2+r_3+...+r_x=e_i$$

对于这个方程的非负整数解的数量为：

$$\binom{e_i+x-1}{x-1}$$

注：这个公式可以通过对$x-1$个 分隔符 ‘$|$’，$e_i$个$*$的排列数量得到，分隔符把区间分成了$x$块，每块区间中$*$的数量对应一个$r$，排列的数量等同于分隔符选位置的方法数。

记：

$$H(x)=\prod_{i=1}^{m}\binom{e_i+x-1}{x-1}$$

首先，特别说明，$H(x)\ne C(x)$

因为，原题要求，每次至少取走一个素因子，而$H(x)$中，可能会有某一次，一个也没有取的情况。

设$x$次取数中，至少有$y$次，

一个素因子也没有取的方法数为$D(x,y)$

则：

$$D(x,y)=\binom{x}{y}H(x-y)=\binom{x}{x-y}H(x-y)$$

应用下式，容斥原理：

对于集合$X_n\{1,2,3,...,n\}$,

对于集合S的子集,$A_1,A_2,A_3,..A_n$

它们对立相反事件为：$A_1^-,A_2^-,A_3^-,...A_n^-$

那么S中不属于$A_1,A_2,A_3,...A_n$的元素数量为：

$$\big|A_1^-\bigcap A_2^-\bigcap A_3^-..\bigcap A_n^-\big|=\sum_{L\subseteq\ X_n}(-1)^{n-|L|}\big|\bigcap_{i\notin L} A_i\big |$$

应用到上式：

对于所有位置都至少取一个素因子的方案数量，等于至少有y个位置不能取，不同y之间的容斥。

$$C(x)=\sum_{y=0}^{x}(-1)^{x-y}D(x,y)$$

$$C(x)=\sum_{y=0}^{x}(-1)^{x-y}\binom{x}{x-y}H(x-y)$$

$$C(x)=\sum_{y=0}^{x}(-1)^{x-y}\binom{x}{y}H(y)$$

$$C(x)=\sum_{y=0}^{x}(-1)^{x-y}\frac{x!}{y!(x-y)!}H(y)$$

所以：

$$\frac{C(x)}{x!}=\sum_{y=0}^x\frac{(-1)^{x-y}}{(x-y)!}\frac{H(y)}{y!}$$

令多项式$A(x)$有：

$$A(x)=\frac{(-1)^0}{0!}x^0+\frac{(-1)^1}{1!}x^1+...+\frac{(-1)^w}{w!}x^w$$

令多项式$B(x)$有：

$$B(x)=\frac{H(0)}{0!}x^0+\frac{H(1)}{1!}x^1+...+\frac{H(w)}{w!}x^w$$

那么有：

$$A(x)B(x)=\frac{C(0)}{0!}x^0+\frac{C(1)}{1!}x^1+...+\frac{C(w)}{w!}x^{2w}$$

利用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
 */