题解 P3726 【[AH2017/HNOI2017]抛硬币】

· · 题解

过气数论魔法……

看了下网上的扩展lucas教程好像写的不是非常清楚……

反正我一开始是绕了不少弯之后才明白怎么写的……

所以就用这道题当作扩展lucas的板子吧……

本题题解

首先我们把小A抛硬币的胜负情况和小B抛硬币的胜负情况看成两个01序列，然后我们把这两个01序列接起来就形成了一个长度为a+b的01序列

现在我们将原问题转化为了对01序列计数的问题

以下设小a的序列中有W_{a}个1，小b的序列中有W_{b}个1

然后我们先考虑一个比较简单的情况a==b

那么我们假设现在有一个胜的01序列，那么一定满足条件W_{a}>W_{b}，由a==b立即推得a-W_{a}>b-W_{b}，因此我们证明了任意一个胜的01序列01反转之后一定是一个负的序列

而显然平的序列01反转之后依然是平的01序列,那么我们其实只需要计算出总方案，之后减去平的方案数，最后除以2即可了……

那么我们计算平的01序列可以直接暴力枚举这个01序列中有多少1来计数，式子长这样

\sum_{i=1}^{a}(C_{a}^{i})^{2}=C_{2a}^{a}

这个式子的证明下面会讲

下面我们来讨论略微复杂一点的情形，a>b

我们还是考虑01反转这个套路，但这时我们发现当W_{a}>W_{b}时我们只能推出-W_{a}<-W_{b}但是由于a>b因此我们无法比较a-W_{a}与b-W_{b}的大小，此时的情况略微的有点辣手……

但是我们反过来考虑呢？考虑一个负或者平的序列，当W_{a}<=W_{b}的时候我们可以推出-W_{a}>=-W_{b}又因为a>b因此我们可以立即推出a-W_{a}>b-W_{b}换句话说，一个负或者平的01序列反转之后一定是一个胜的01序列

因此我们发现只有部分胜的01序列反转之后还是胜的01序列，换句话来讲，只有这部分01序列没有“配对”，因此如果我们在总方案里加上这部分01序列，那么所有的胜的01序列均可以配对，此时只需要将方案数除2就是胜的01序列方案

那么我们需要找到一个条件使得胜的01序列反转之后仍然是一个胜的01序列

我们可以列出一个不等式组

W_{a}>W_{b}且a-W_{a}>b-W_{b}

我们可以解出这样一个比较优美的形式

a-b>W_{a}-W_{b}>0

那么我们呢就可以通过枚举W_{b}以及W_{a}-W_{b}来计数了……式子大概长这样

\sum_{i=0}^{b}\sum_{j=1}^{a-b-1}C_{b}^{i}C_{a}^{i+j}

我们交换一下Σ，即swap(i,j)得到

\sum_{i=1}^{a-b-1}\sum_{j=0}^{b}C_{b}^{j}C_{a}^{i+j}

=\sum_{i=1}^{a-b-1}\sum_{j=0}^{b}C_{b}^{b-j}C_{a}^{i+j}

观察到b-j+i+j=b+i因此我们可以重写内层Σ

\sum_{i=1}^{a-b-1}\sum_{j+k=b+i}C_{b}^{j}C_{a}^{k}

这个形式被称为范德蒙德卷积，即以下等式恒成立

\sum_{i+j=k}C_{a}^{i}C_{b}^{j}=C_{a+b}^{k}

所以我们的式子被化成了

\sum_{i=1}^{a-b-1}C_{a+b}^{b+i}

(备注：根据范德蒙德卷积我们也可以证明

\sum_{i=1}^{a}(C_{a}^{i})^{2}=\sum_{i=1}^{a}C_{a}^{i}C_{a}^{a-i}=C_{2a}^{a}

这个式子是成立的)

然后我们加上一个方案数然后除以2的话我们可以得到最后的答案式子

若A==B

Ans=\frac{2^{A+B}-C_{2a}^{a}}{2}

若A>B

Ans=\frac{2^{A+B}+\sum_{i=1}^{A-B-1}C_{A+B}^{B+i}}{2}

然后问题来了，我们怎么计算C_{n}^{m}mod10^9?

如果mod的是一个不大的质数，那么我们可以通过lucas定理轻松的求解，但是如果p不是质数呢？此时我们需要一种被称之为扩展lucas的上古魔法来帮助我们求解这类问题

拓展Lucas

首先由于我们的模数是任意的，我们将模数进行唯一分解，分解成P^{k}p为素数的连乘积形式，此时我们只需要分别计算出C_{n}^{m}modp^{k}最后使用中国剩余定理解同余方程就行了

问题来了怎么求出C_{n}^{m}modp^{k}呢？

此时网上的教程会告诉你我们只需要求出N!modp^{k}最后除在一起就好了……

但是问题来了，分母有极大的概率是p^{k}的倍数，此时不存在乘法逆元

因此我们需要略微费点功夫，我们将N!表示为a×P^{b}的形式

然后计算除法的时候各个阶乘的a部分乘法逆元都存在，然后在b部分运行指数上的加减法就好了

最后判以下指数是否大于k如果大于k直接输出0即可，否则就合并两个部分的值就行了

最后的问题，如何求出N!modP^{k}的不含有因子p的部分呢？

我们举一个非常老套的例子，计算19!mod3^2

把阶乘一字拆开得到

1×2×3×4×5×6×7×8×9×10×11×12×13×14×15×16×17×18×19

然后提走所有含有因子3的部分一个3式子变成了

6!×3^{6}×1×2×4×5×7×8×10×11×13×14×16×17×19

然后我们呢将后边的部分膜上一个3^2得到了

6!×3^{6}×(1×2×4×5×7×8)^{2}×19

让我们来观察各个部分的特征，左侧的阶乘可以直接递归的去处理，3^{6}可以单独计算，而中间的部分其实是一个伪阶乘的形式，即去掉了3的倍数的阶乘，发现最长不超过p^{k}可以打表预处理，而最后剩下的19其实是一个未完的周期，已经在我们的表当中了，因此如果记得伪阶乘为B_{i}那么我们可以得到这样的递归式来求阶乘的有乘法逆元的部分

lucas(n,p^{k})=lucas(\lfloor \frac{n}{p}\rfloor ,p^{k})B_{p^{k}-1}^{\lfloor \frac{n}{p^{k}} \rfloor}B_{n\%p^{k}}

其中预处理是O(p^k)的，每次询问阶乘是O(log^{2}n)的(计入了快速幂的log)

最后使用中国剩余定理解同余方程即可得到C(n,m)modp了

最后一个问题，2在膜10^9的意义下乘法逆元是不存在的，我们该怎么除2呢？

我们发现如果a+b是奇数，那么我们需要计算的式子其实是杨辉三角中间对称的几列，那么我们只计算左边的几列或者右边的几列就可以除2了，但是在a+b为偶数的时候我们会发现剩下了一个C_{a+b}^{(a+b)/2}情况会有点辣手

但是我们注意到根据帕斯卡恒等式(杨辉三角)可以得到

C_{2a}^{a}=C_{2a-1}^{a}+C_{2a-1}^{a-1}=2C_{2a-1}^{a}

于是我们计算C_{a+b-1}^{(a+b)/2}就可以了

上代码~

#include<cstdio>
#include<algorithm>
using namespace std;typedef unsigned long long ll;const ll mod[2]={512,1953125};const ll bs[2]={2,5};
const ll M=1e9;//省事起见直接mod1e9了 
inline ll po(ll a,ll p,ll mod){ll r=1;for(;p;p>>=1,a=a*a%mod)if(p&1)r=r*a%mod;return r;}
inline void exgcd(ll a,ll& x,ll b,ll& y)//乘法逆元和快速幂 
{if(b==0){x=1;y=0;return;}exgcd(b,x,a%b,y);ll tp=x;x=y;y=tp-a/b*y;}
inline ll inv(ll a,ll b){ll x;ll y;exgcd(a,x,b,y);return (x+b)%b;}
ll bru[2][1953200];ll a;ll b;int k;int nb[15];
ll lucas(ll n,ll t)//一行lucas 
{return n?lucas(n/bs[t],t)*bru[t][n%mod[t]]%mod[t]*po(bru[t][mod[t]-1],n/mod[t],mod[t])%mod[t]:1;}
inline ll c(ll n,ll m,int t)
{
    ll d=0;for(ll i=n;i;i/=bs[t])d+=i/bs[t];
    for(ll i=m;i;i/=bs[t])d-=i/bs[t];
    for(ll i=n-m;i;i/=bs[t])d-=i/bs[t];if(d>=9){return 0;}//剪枝，我们直接暴力计算出指数然后判断是否是0 
    ll p1=lucas(n,t);ll p2=lucas(m,t);ll p3=lucas(n-m,t);//计算非互质部分 
    return po(bs[t],d,mod[t])*p1%mod[t]*inv(p2,mod[t])%mod[t]*inv(p3,mod[t])%mod[t];//乘上p的幂次 
}
inline void prit(ll x)//打印函数 
{
    for(int i=1;i<=10;i++)nb[i]=0;
    for(int i=1;x;i++,x/=10)nb[i]=x%10;
    for(int i=k;i>=1;i--)printf("%d",nb[i]);printf("\n");
}
int main()//打表伪阶乘 
{
    bru[0][0]=1;for(int i=1;i<mod[0];i++){bru[0][i]=(i%2)?i*bru[0][i-1]%mod[0]:bru[0][i-1];}
    bru[1][0]=1;for(int i=1;i<mod[1];i++){bru[1][i]=(i%5)?i*bru[1][i-1]%mod[1]:bru[1][i-1];}
    while(scanf("%lld%lld%d",&a,&b,&k)!=-1)
    {
        ll ret2=0;ll ret5=0;
        if(a==b)//特判a==b 
        {
            ret2=c(2*a-1,a,0);ret5=c(2*a-1,a,1);
            ret2=ret2*mod[1]%M*inv(mod[1]%mod[0],mod[0])%M;
            (ret2+=ret5*mod[0]%M*inv(mod[0],mod[1]))%=M;ret2=M-ret2;
            (ret2+=po(2,a+b-1,M))%=M;prit(ret2);//中国剩余定理解同余方程 
        }
        else 
        {
            if((a+b)%2)//如果a+b是奇数 
            {
                for(ll i=b+1;i<=(a+b)/2;i++)(ret2+=c(a+b,i,0))%=mod[0];
                for(ll i=b+1;i<=(a+b)/2;i++)(ret5+=c(a+b,i,1))%=mod[1];//只算一半 
                ret2=ret2*mod[1]%M*inv(mod[1]%mod[0],mod[0])%M;
                (ret2+=ret5*mod[0]%M*inv(mod[0],mod[1]))%=M;(ret2+=po(2,a+b-1,M))%=M;
                prit(ret2);
            }
            else //偶数的时候通过计算C(a+b-1,(a+b)/2)来将C(a+b,(a+b)/2)除2 
            {
                for(ll i=b+1;i<(a+b)/2;i++)(ret2+=c(a+b,i,0))%=mod[0];(ret2+=c(a+b-1,(a+b)/2,0))%=mod[0];
                for(ll i=b+1;i<(a+b)/2;i++)(ret5+=c(a+b,i,1))%=mod[1];(ret5+=c(a+b-1,(a+b)/2,1))%=mod[1];               
                ret2=ret2*mod[1]%M*inv(mod[1]%mod[0],mod[0])%M;
                (ret2+=ret5*mod[0]%M*inv(mod[0],mod[1]))%=M;(ret2+=po(2,a+b-1,M))%=M;
                prit(ret2);
            }
        }
    }return 0;//拜拜程序~ 
}