洛谷 P3726 [AH2017/HNOI2017]抛硬币 exgcd+扩展Lucas定理

题目描述
小 A 和小 B 是一对好朋友，他们经常一起愉快的玩耍。最近小 B 沉迷于**师手游，天天刷本，根本无心搞学习。但是已经入坑了几个月，却一次都没有抽到 SSR，让他非常怀疑人生。勤勉的小 A 为了劝说小 B 早日脱坑，认真学习，决定以抛硬币的形式让小 B 明白他是一个彻彻底底的非洲人，从而对这个游戏绝望。两个人同时抛 b 次硬币，如果小 A 的正面朝上的次数大于小 B 正面朝上的次数，则小 A 获胜。

但事实上，小 A 也曾经沉迷过拉拉游戏，而且他一次 UR 也没有抽到过，所以他对于自己的运气也没有太大把握。所以他决定在小 B 没注意的时候作弊，悄悄地多抛几次硬币，当然，为了不让小 B 怀疑，他不会抛太多次。现在小 A 想问你，在多少种可能的情况下，他能够胜过小 B 呢？由于答案可能太大，所以你只需要输出答案在十进制表示下的最后 k 位即可。

输入输出格式
输入格式：
有多组数据，对于每组数据输入三个数a,b,k，分别代表小A抛硬币的次数，小B抛硬币的次数，以及最终答案保留多少位整数。

输出格式：
对于每组数据，输出一个数，表示最终答案的最后 k 位为多少，若不足 k 位以 0 补全。

输入输出样例
输入样例#1：
2 1 9
3 2 1
输出样例#1：
000000004
6

说明
对于第一组数据，当小A抛2次硬币，小B抛1次硬币时，共有4种方案使得小A正面朝上的次数比小B多。
(01,0), (10,0), (11,0), (11,1)
对于第二组数据，当小A抛3次硬币，小B抛2次硬币时，共有16种方案使得小A正面朝上的次数比小B多。
(001,00), (010,00), (100,00), (011,00), (101,00), (110,00), (111,00), (011,01), (101,01), (110,01),(111,01), (011,10), (101,10), (110,10), (111,10), (111,11).

数据范围
10%的数据满足a,b≤20；
30%的数据满足a,b≤100；
70%的数据满足a,b≤100000，其中有20%的数据满足a=b；
100%的数据满足$1\le a,b\le 10^{15},b\le a\le b+10000,1\le k\le 9$，数据组数小于等于10。

分析：
假设一种方案为$s$，另一种方案$s^{\prime }$为$s$的每一位取反。
假设$a=b$，那么$s$为胜利状态，那么$s^{\prime }$必然为失败状态，只要减去平局状态再除以2即可。即$\frac{2^{a+b}-\left(\genfrac{}{}{0ex}{}{a+b}{a}\right)}{2}$。
对于$a>b$，假设$s$为平局或失败状态，$s^{\prime }$必然为胜利状态。还有一些胜利状态$s$对应胜利状态$s^{\prime }$。
因为$s$和$s^{\prime }$一一对应，显然$s^{\prime }$所代表的胜利状态恰好为$2^{a+b-1}$种，剩下只要加上$s$中的胜利状态个数。

现在这一类中$1$的个数满足，
$numa>numb$
$a-numa>b-numb$
所以
$a-b>numa-numb$
枚举其中一个的$1$的个数，那么这一个部分的答案为
$ans={\sum }_{i=0}^b{\sum }_{j=1}^{a-b-1}\left(\genfrac{}{}{0ex}{}{b}{i}\right)\ast \left(\genfrac{}{}{0ex}{}{a}{i+j}\right)={\sum }_{i=0}^b{\sum }_{j=1}^{a-b-1}\left(\genfrac{}{}{0ex}{}{b}{b-i}\right)\ast \left(\genfrac{}{}{0ex}{}{a}{i+j}\right)$
后面相当于直接从$a+b$中选$j$个，枚举$i$只是为了判断一个子集选多少个。
也就是求
$\frac{{\sum }_{i=1}^{a-b-1}\left(\genfrac{}{}{0ex}{}{a+b}{b+j}\right)}{2}$
总的方案数是$ans=\frac{2^{a+b}+{\sum }_{i=1}^{a-b-1}\left(\genfrac{}{}{0ex}{}{a+b}{b+j}\right)}{2}$

我们可以直接把模数设为$10^9$，后面只需要保留相应位数即可。
因为模数不是质数，我们可以先把模数因式分解成$mod=p_1^{a_1}\ast p_2^{a_2}\ast ...\ast p_m^{a_m}$的形式，对于每一个$p_i^{a_i}$次方分别求解，然后再中国剩余定理合并答案。

考虑怎样求$\left(\genfrac{}{}{0ex}{}{n}{m}\right)mod{p}^k$，其中$p$为质数。
我们把$n!$拆成$a\ast p^b$的形式，然后再计算。
其中$b$为$1$到$n$所有数包含$p$这个素因子的个数的和，即$b=\frac{n}{p}+\frac{n}{p^2}+\frac{n}{p^3}+....$
显然$n!=1\ast 2\ast 3\ast ....\ast n$
我们把是$p$的倍数的放在一边，就是
$n!=({\prod }_{i=1}^{p-1}i)\ast ({\prod }_{i=p+1}^{2\ast p-1}i)\ast ...\ast p\ast (2\ast p)\ast (3\ast p)\ast ....$
$=({\prod }_{i=1}^{p-1}i)\ast ({\prod }_{i=p+1}^{2\ast p-1}i)\ast ...\ast p^d\ast (n/p)!$

显然${\prod }_{i=1}^{p-1}i={\prod }_{i=p^k+1}^{p^k+p-1}i(mod{p}^k)$
也就是存在循环节，我们需要把一个循环节算$n/(p^k)$次方，再把多出来的乘上就可以了，$p^d$已经被统计进$b$中了。
然后阶乘要预处理，这里的阶乘是除掉$p$的倍数的，后面就递归下去即可。
然后再计算组合数时，$p$的指数直接减掉，而$a$的部分可以求逆元，显然$a$与$p^k$互质，考虑exgcd。
这样就算出组合数了。
最后要把结果除以$2$，这些组合数可以发现是两两对称的，我们只需要统计一半。如果中间多了一个，这个一定是$\left(\genfrac{}{}{0ex}{}{a+b}{(a+b)/2}\right)$，这个数一定是$2$的倍数。

代码：

// luogu-judger-enable-o2
#include <iostream>
#include <cstdio>
#include <cmath>
#define LL long long

const int maxn=2e6+7;
const LL mod[2]={512,1953125};
const LL M=1e9;

using namespace std;

LL n,m,k,ans;
LL jc[maxn][2];

void prework() //预处理阶乘，注意不含p的倍数
{
    jc[0][0]=1;
    for (LL i=1;i<mod[0];i++)
    {
        if (i%2) jc[i][0]=jc[i-1][0]*(LL)i%mod[0];
            else jc[i][0]=jc[i-1][0];
    }
    jc[0][1]=1;
    for (LL i=1;i<mod[1];i++)
    {
        if (i%5) jc[i][1]=jc[i-1][1]*(LL)i%mod[1];
            else jc[i][1]=jc[i-1][1];
    }
}

LL ksm(LL x,LL y,LL k)
{
    if (y==0) return 1;
    LL c=ksm(x,y/2,k);
    c=(c*c)%k;
    if (y&1) c=(c*x)%k;
    return c;
}

LL exgcd(LL a,LL b,LL &x,LL &y)
{
    if (b==0)
    {
        x=1;
        y=0;
        return a;
    }
    LL d=exgcd(b,a%b,x,y);
    LL z=x;
    x=y,y=z-a/b*y;
    return d;
}

LL inv(LL a,LL mod)
{
    LL x,y;
    exgcd(a,mod,x,y);
    return (x%mod+mod)%mod;
}

LL multi(LL n,LL pi,LL pk) //求出不含p的部分的积
{
    if (!n) return 1; //递归结束条件
    LL c=ksm(jc[pk-1][pi!=2],n/pk,pk)*jc[n%pk][pi!=2]%pk; //循环节有n/pk个，每个都是满的，直接乘方后再乘上多出来的n%pk个
    return c*multi(n/pi,pi,pk)%pk; //递归处理
}

LL binom(LL n,LL m,LL pi,LL pk,bool flag)
{
    LL p=0;
    for (LL i=pi;i<=n;i*=pi) p+=n/i;
    for (LL i=pi;i<=m;i*=pi) p-=m/i;
    for (LL i=pi;i<=n-m;i*=pi) p-=(n-m)/i; //算出p的幂次
    if ((flag) && (pi==2)) p--; //要对组合数除以2的情况
    if (p>=9) return 0;
    LL a=multi(n,pi,pk),b=multi(m,pi,pk),c=multi(n-m,pi,pk);
    LL cs=a*inv(b,pk)%pk*inv(c,pk)%pk*ksm(pi,p,pk)%pk;
    if ((flag) && (pi==5)) cs=cs*inv(2,pk)%pk; 
    return cs;
}

LL getans(LL n,LL m,bool flag) //中国剩余定理合并，因为打了exgcd，可以利用上
{
    LL a=binom(n,m,2,mod[0],flag); 
    LL b=binom(n,m,5,mod[1],flag);
    LL x,y;
    exgcd(mod[0],mod[1],x,y);
    x=(x*(b-a)%M*mod[0]%M+a+M)%M;	
    return x;
}

void prin(LL x)
{
    int b[20],num=0;
    for (int i=1;i<=k;i++) b[i]=0;
    while (x)
    {
        b[++num]=x%10;
        x/=10;
    }
    for (int i=k;i>0;i--) printf("%d",b[i]);
    printf("\n");
}

int main()
{
    prework();
    while (scanf("%lld%lld%lld",&n,&m,&k)!=EOF)
    {		
        if (n==m) ans=(ksm(2,n+m-1,M)+M-getans(n+m,n,1))%M;
        else
        {
            ans=ksm(2,n+m-1,M);			
            for (LL i=(n+m)/2+1;i<=n-1;i++) ans=(ans+getans(n+m,i,0))%M;
            if ((n+m)%2==0) ans=(ans+getans(n+m,(n+m)/2,1))%M;
        }
        prin(ans);
    }
}