记录编号	536543	评测结果	AAAAAAAAAA
题目名称	[HSOI 2018] 简单题233	最终得分	100
用户昵称	..	是否通过	通过
代码语言	C++	运行时间	1.186 s
提交时间	2019-07-08 09:58:21	内存使用	13.66 MiB

/*
1.因为p不是质数，所以可以用唯一分解定理，解得出，连续素数的指数乘积
2.那么分解完后每一项，pi^ai之间两两互质，所以只要我们能得出每组c(m,n)%(pi^ai)的解再用中国剩余定理合并即可
3.首先我们除一个数，等于乘上这个数的倒数，就相当于乘上逆元，
那么显然我们要求出阶乘和阶乘的逆元。由于阶乘与模数当前不互质，所以极有可能不存在直接的逆元，需要我们对式子进行进一步的拆分。 
n!=p^(n/p)(向下取整)*（n/p）！向下取整*（p整除i的数的乘积） 
22! mod 3^2
原式(22!)=(1*2*3*4*5*6*7*8*9)*(10*11*12*13*14*15*16*17*18）*(19*20*21*22)
现在化简后的式子(22!)=(1*2*4*5*7*8)*(10*11*13*14*16*17)*(19*20*22)*36*[(1*2*3*4*5*6*7)==（n/p）向下取整的！] 
发现按照如下分组之后前(n/pi^ki) 组在mod pi^ki 后结果相同，因此只需要做一组然后快速幂即可.
所以
式子中的p^(n/p向下取整)可以直接快速幂求出，
（n/p）！向下取整考虑递归求解，唯一比较麻烦的是后面剩下的不能被pp整除的数
可以发现，剩下这些项的乘积有循环节，长度小于p^a。由于显然x≡x+p^a(modp^a)，
所以可以求出共有几个循环节，对第一个循环节暴力求一遍，然后用快速幂搞定。
对于最后剩余的不能构成完整循环节的几项，依然是暴力求一遍搞定。0! = 1
首先我们把给定的模数p
p进行质因数分解。
然后对于每个pi^ai
pi^ai，都求一遍组合数c(m,n)mod(pi^ai)然后同时用crt合并。
*/ 


#include <iostream>
#include<cstdio>

using namespace std;
typedef long long ll;
const ll p=23333333;
ll n, m;
//const ll maxn=23333333;

ll exgcd(ll a, ll b, ll &x, ll &y){//扩欧求逆元 
	if(!b){
		x = 1, y = 0;
		return a;
	}
	ll res = exgcd(b, a%b, x, y);
	ll t = x;
	x = y;
	y = t-a/b*y;
	return res;
}

ll qpow(ll a, ll n, ll mod){//快速幂 
	ll res = 1;
	while(n){
		if(n&1) res = (res*a) % mod;
		a = (a*a) % mod;
		n >>= 1;
	}
	return res;
}

ll inv(ll a, ll p){//逆元的处理 
	ll x, y;
	exgcd(a, p, x, y);
	if(x+p > p) return x;
	return x+p;
}

ll crt(ll n, ll mod){
	return n*(p/mod)%p*inv(p/mod, mod)%p;//中国剩余定理 
}

ll fac(ll n, ll p, ll k) {      //n! % k (k=p^a)
    if(!n) return 1;		//0! = 1
    ll ans = 1;
    for(int i = 2; i <= k; i++) {		//处理循环节内的数字(22%9中的(1*2*4*5*7*8))
        if(i%p) ans = ans*i % k;
    }
    ans = qpow(ans, n/k, k);		//所有循环节的乘积(1*2*4*5*7*8)和(10*11*13*14*16*17)mod出来是一样的所以用快速幂直接算出来 
    for(int i = 2; i <= n%k; i++) {		//剩下单独的项（19*20*22） 
        if(i%p) ans = ans*i % k;
    }
    return ans*fac(n/p, p, k)%k;		//递归求解
    /*
	这样就可以在可承受的时间复杂度内求出阶乘。
	但是，为了达到除法的效果，我们需要考虑p的次幂一共出现了多少次。
	根据可以知道只除去一个p时，n!内包括了(n/p)xiazheng个p！
	剩下的数字如果要可能存在p的次幂也可以归为一个阶乘，即
 	n/p 下整的阶乘 
	所以f(n)=f(n/p)+[n/p]xiazheng 
	*/ 
}

/*
求C(n,m)%pk
求不了逆元，因为要求分母和模数互质，
所以C（nm）%pk=
(n!/(p^a1))/((m!/(p^a2))*((n-m)!/(p^a3))*p^(a1-a2-a3)modP^k
*/
ll C(ll n, ll m, ll p, ll k){	//k = p^x
   	if(n < m) return 0;
   	//a是n!/(p^a1))，b是(m!/(p^a2)，c=(n-m)!/(p^a3)
	ll a = fac(n,p,k), b = fac(m,p,k), c = fac(n-m,p,k);
	ll cnt = 0;
	/*
	上面说到，要在求组合数时单独处理形如p^x的项，
	那么具体的操作就是：我们用一个cnt变量记录最后的式子中有多少个因子p，
	那么cnt就等于n的阶乘中含有因子p的数量减去m和n-m的阶乘中含有的p的数量。
	*/
	for(ll i = p; i <= n; i *= p) cnt += n/i;
	for(ll i = p; i <= m; i *= p) cnt -= m/i;
	for(ll i = p; i <= n-m; i *= p) cnt -= (n-m)/i;
	return a*inv(b, k)%k * inv(c, k)%k * qpow(p, cnt, k)%k;
}

ll exlucas(){
	ll t = p, ans = 0;
	for(ll i = 2; i*i <= t; i++){//质因数分解 
		if(t%i) continue;
		ll tmp = 1;
		while(t%i == 0){
			tmp *= i;
			t /= i;
		}
		ans = (ans+crt(C(n+m-2 ,m-1, i, tmp)/*先算出c（m+n，n）modpi^ai*/, tmp))%p;
	}
	//说明p没有被完全合并 
	if(t > 1) ans = (ans+crt(C(n+m-2, m-1, t, t), t))%p;
	return ans%p;
}

int main()
{
	    	freopen("hst.in","r",stdin);
    	freopen("hst.out","w",stdout);
	cin >> n >> m ;
	cout << exlucas() << endl;
	
	return 0;
}