BZOJ3738 [Ontak2013]Kapitał

题目：
给出三个数字(N,M,K)。求(C_{N+M}^N)去掉所有末尾的(0)后对(10^K)取模的结果.
(1leq N,Mleq 10^15,1leq kleq 9)

分析：
把(10^K)分解为(2^K 5^K)
后面CRT合并就好了
两个都是质数，处理方法相同，下面直接用质数(P)代替
考虑把(X=C_{N+M}^N)表示为(aP^b)的形式，其中(gcd(a,P)=1)
我们将(a)对(P^K)取模
这是一个类似于ExLucas的递归，我们先预处理一下(P^K)以内的前缀积（去除(P)的倍数）
这样可以(O(P^K))预处理，(O(logn))查询
于是有了两个等式：
(frac{X}{2^{b_1}} equiv a_1(mod 2^K))
(frac{X}{5^{b_2}} equiv a_2(mod 5^K))
我们先要把后缀(0)去掉，假设有(m)个，相当于
将(a_1)乘上((5^{-1})^{m})（模(2^K)意义下的），(b_1)减去(m)
将(a_2)乘上((2^{-1})^{m})（模(5^K)意义下的），(b_2)减去(m)
把(a_{1}2^{b_1})和(a_{2}5^{b_2})CRT合并就是答案

#include<cstdio>
#include<cstring>
#include<cmath>
#include<iostream>
#include<queue>
#include<algorithm>
#include<vector>
#include<set>

#define maxn 500005
#define INF 0x3f3f3f3f
#define MOD 998244353

using namespace std;

inline long long getint()
{
	long long num=0,flag=1;char c;
	while((c=getchar())<'0'||c>'9')if(c=='-')flag=-1;
	while(c>='0'&&c<='9')num=num*10+c-48,c=getchar();
	return num*flag;
}

long long n,m,k;
long long T,P,Pk;
long long ans0,ans1,ans;
long long s[1953130];
inline void exgcd(long long p,long long q,long long &x,long long &y)
{
	if(!q){x=1,y=0;return;}
	exgcd(q,p%q,y,x),y-=p/q*x;
}
inline long long ksm(long long num,long long k,long long P)
{
	long long ret=1;
	for(;k;k>>=1,num=num*num%P)if(k&1)ret=ret*num%P;
	return ret;
}
long long inv(long long p,long long q){long long x,y;exgcd(p,q,x,y);return (x%q+q)%q;}
struct node{
	long long a,b;
	node(){}node(long long x,long long y){a=x,b=y;}
	node operator *(node x){return node(a*x.a%Pk,b+x.b);}
	node operator /(node x){return node(a*inv(x.a,Pk)%Pk,b-x.b);}
}st[2];

inline node cal(long long n){return n?node(s[n%Pk]*ksm(s[Pk],n/Pk,Pk)%Pk,n/P)*cal(n/P):node(1,0);}
inline void init(){s[0]=1;for(int i=1;i<Pk;i++)if(i%P)s[i]=s[i-1]*i%Pk;else s[i]=s[i-1];s[Pk]=s[Pk-1];}

int main()
{
	n=getint(),m=getint(),k=getint();
	T=1;while(k--)T*=10;
	P=2,Pk=512,init();
	st[0]=cal(n+m)/cal(n)/cal(m);
	P=5,Pk=1953125,init();
	st[1]=cal(n+m)/cal(n)/cal(m);
	int tmp=min(st[0].b,st[1].b);
	while(tmp--)Pk=512,st[0]=st[0]/node(5,1),Pk=1953125,st[1]=st[1]/node(2,1);
	P=2,Pk=512,ans0=st[0].a*ksm(P,st[0].b,Pk)%Pk;
	P=5,Pk=1953125,ans1=st[1].a*ksm(P,st[1].b,Pk)%Pk;
	ans=(1953125ll*inv(1953125,512)%T*ans0%T+512ll*inv(512,1953125)%T*ans1%T)%T;
	while(ans*10<T)putchar('0'),T/=10;
	return printf("%lld",ans),0;
}