Codeforces Round #181 (Div. 2) Problem C 组合数取模 逆元

部分转自这里

题意：用 a 和 b 组成 n 位的数字 (1 ≤ n ≤ 10^6₎  ，如果各位数字之和仍只由 a 和 b 组成，则记为“good number”，问一共有多少个good number？ 最终答案对1000000007 (10⁹ + 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;
}