题解 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}$就可以了 上代码~ ```C #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;//拜拜程序~ } ```