Codeforces Round #181 (Div. 2) Problem C 组合数取模 逆元
部分转自这里
题意:用 a 和 b 组成 n 位的数字 (1 ≤ n ≤ 106) ,如果各位数字之和仍只由 a 和 b 组成,则记为“good number”,问一共有多少个good number? 最终答案对1000000007 (109 + 7)取模。
分析:假设有k个a,那么b的个数为n-k , 各位数字之和为 s = k*a+(n-k)*b, 如果s为good number,那么一共有 C(n,k)种方法。
C(n,k) = n! / k! / (n-k)! ,加,减,乘的运算可以取余,但是除法不能直接取余,a/b%mod = a*reverse(b)%mod; 所以这里要对k! 和 (n-k)! 取逆元。
定理:a存在模p的乘法逆元的充要条件是gcd(a,p) = 1
证明:
充分性: 如果gcd(a,p) = 1,根据欧拉定理,a^φ(p) ≡ 1 mod p,因此 显然a^(φ(p)-1) mod p是a的模p乘法逆元。
必要性: 假设存在a模p的乘法逆元为b ab ≡ 1 mod p 则ab = kp +1 ,所以1 = ab - kp 因为gcd(a,p) = d 所以d | 1 所以d只能为1
求逆元:
方法1:
(p是mod) 首先这里要保证 gcd(b,p)=1
根据费马小定理 若gcd(b,p)=1 则 b^phi(p)== 1(mod p), a / b * 1==a /b * b^phi(p)== a * b^(phi(p)-1) (mod p) 这里要用到欧拉函数二分求幂
快速求a ^ b : O(logn)
把b按二进制展开为:b = p(n)*2^n + p(n-1)*2^(n-1) +…+ p(1)*2 + p(0) 其中p(i) 为 0 或 1
这样 a^b = a^( p(n)*2^n + p(n-1)*2^(n-1) +...+ p(1)*2 + p(0) )= a^(p(n)*2^n) * a^(p(n-1)*2^(n-1)) *...* a^(p(1)*2) * a^p(0)
p(i)=0的情况, a^(p(i) * 2^(i-1) ) = a^0 = 1,不用处理
p(i)=1: a^(2^i) = a^( 2^(i-1)*2 ) = ( a^( 2^(i-1) ) ) ^ 2
欧拉
const LL mod = (LL)1e9+7; LL A[1000005] = {1, 1}; bool judge(int n){ while( n ){ int t = n % 10; n /= 10; if(t != a && t != b)return false; } return true; } LL ans; int a,b,n; int euler1(int n){//求φ(n) int res=n; for(int i=2; i<= (int)sqrt( n*1.0 ); i++) if(n%i==0){//i是n的质因数 res=res/i*(i-1); while(n%i==0)n/=i; } if(n>1) res=res/n*(n-1);//n是质因数 return res; } //计算 a ^ b mod n int modexp(int a, int b, int n) { LL ans = 1, t = a; while(b) { if(b & 1) ans = ans * t % n; t = t * t % n; b >>= 1; } return ans; } int main(){ #ifndef ONLINE_JUDGE freopen("in.txt","r",stdin); #endif FOR(i, 2, 1000004) A[i] = A[i-1] * i % mod; int pow = euler1(mod) - 1; cin>>a>>b>>n; FOE(k,0,n) if(judge(k*a+(n-k)*b)){ //cout<<"AAAAAAAAAA"<<endl; LL t = A[n]; LL t1 = modexp(A[k], pow, mod); LL t2 = modexp(A[n-k], pow, mod); ans = (ans + ((t*t1)%mod * t2) )%mod; } cout<<ans<<endl; return 0; }
方法2:
若gcd(b,p)=1,则可根据扩展欧几里得算法得出 bx==1 (mod p) , a/b*(b*x)==ax (mod p)
扩展欧几里德算法,就是已知a、b,求一组解(x,y)使得a*x+b*y=1。这里求得的x即为a%b的逆元,y为b%a的逆元(方程两边都模上b)。
扩展gcd
const LL mod = (LL)1e9+7; LL A[1000005] = {1, 1}; int a,b,n; LL ans; bool judge(int n){ while( n ){ int t = n % 10; n /= 10; if(t != a && t != b)return false; } return true; } LL extgcd(LL a, LL b, LL & x, LL & y){ if (b == 0) { x=1; y=0; return a; } LL d = extgcd(b, a % b, x, y); LL t = x; x = y; y = t - a / b * y; return d; } int main(){ #ifndef ONLINE_JUDGE READ #endif FOR(i, 2, 1000004) A[i] = A[i-1] * i % mod; cin>>a>>b>>n; LL t,t1,t2,t3; FOE(k,0,n) if(judge(k*a+(n-k)*b)){ t = A[n]; extgcd(A[k], mod, t1, t3); extgcd(A[n-k], mod, t2, t3); t1 = (t1 % mod + mod) % mod; t2 = (t2 % mod + mod) % mod; ans = (ans + ((t*t1)%mod * t2) )%mod; } cout<<ans<<endl; return 0; }