HDU 1576 -- A/B (总结乘法逆元的几种求法)

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=1576

A/B

Time Limit: 1000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 7264    Accepted Submission(s): 5774

Problem Description

要求(A/B)%9973，但由于A很大，我们只给出n(n=A%9973)(我们给定的A必能被B整除，且gcd(B,9973) = 1)。

 

 

Input

数据的第一行是一个T，表示有T组数据。
每组数据有两个数n(0 <= n < 9973)和B(1 <= B <= 10^9)。

 

 

Output

对应每组数据输出(A/B)%9973。

 

 

Sample Input

2

1000 53

87 123456789

 

 

Sample Output

7922

6060

 

 

题意：给出A%9973和B的值，求(A/B)%9973。

解析：(A  /  B) % p = (A * inv(B) ) % p = (A % p * inv(B) % p) % p，其中p为模数9973, inv(B)为B关于模p的逆元。

 

同余式：

如果两个正整数a和b之差能被n整除，那么我们就说a和b对模n同余，记作：a≡b (mod n)。 
a≡b(mod n)等价于a与b分别用n去除，余数相同。

 

乘法逆元：

如果ax≡1 (mod p),且gcd(a,p)=1（a与p互质，此为一个数有逆元的充要条件，此时逆元唯一存在），则称a关于模p的乘法逆元为x。

逆元的含义：模n意义下，1个数a如果有逆元x，那么除以a相当于乘以x。

 

求解逆元的方法：

1. 扩展欧几里得算法：

已知整数a、b，扩展欧几里得算法可以在求得a、b的最大公约数的同时，能找到整数x、y（其中一个很可能是负数），使它们满足贝祖等式ax + by = gcd(a, b)。

当a关于模b的逆元存在，有gcd(a, b) == 1，即扩展欧几里得算法可求得x, y满足ax + by = 1, 两边同时余b，

ax % b + by % b = 1 % b 

ax % b = 1 % b

ax ≡ 1(mod b)

所以x是a的模b乘法逆元，同理y是b的模a乘法逆元。

算法时间复杂度：O(logn)

模板代码：

int ex_gcd(int a, int b, int &x, int &y) {  // 函数返回gcd(a, b)
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    int r = ex_gcd(b, a % b, y, x);
    y -= (a / b) * x;
    return r;
}

int main() {
    int a, b, x, y;
    cin >> a >> b;  // 求a关于模b的逆元
    cout << (ex_gcd(a, b, x, y) == 1 ? (x % b + b) % b : -1) << endl;  // -1表示逆元不存在

    return 0;
}

 

2. 费马小定理：

内容：假如a是一个整数，p是一个质数，那么a^p - a是p的倍数，可以表示为

a^p ≡ a (mod p)

如果a不是p的倍数(即gcd(a, p) == 1)，这个定理也可以写成

a^p-1 ≡ 1 (mod p)

变形得a * a^p-2 ≡ 1 (mod p)，所以a关于模p的逆元x = a^p-2 (mod p)，用快速幂模可快速求之。

算法时间复杂度：O(logn)

模板代码：

LL pow_mod(LL a, LL b, LL p) {    //a的b次方取模p 
    LL ret = 1;
    while (b) {
        if(b & 1) ret = (ret * a) % p;
        a = (a * a) % p;
        b >>= 1;
    }
    return ret;
}
LL Fermat(LL a, LL p) {   //费马小定理求a关于b的逆元 
        return pow_mod(a, p-2, p);
}

 

3. 欧拉定理：

欧拉函数：对正整数n，欧拉函数是小于n的正整数中与n互质的数的数目（φ(n) = n(1 – 1/p₁)(1 – 1/p₂)…(1 – 1/p_m), p_n为n的所有质因数,  φ(1) = 1）。

欧拉定理：若gcd(a, p) = 1，则a^φ(p) ≡ 1 (mod p)。

即 a*a^(φ(p)−1) ≡ 1(mod p)，那么a关于模p的逆元x = a^(φ(p)−1) (mod p)。

(当p为素数的时候φ(p)=p-1,则φ(p)-1=p-2可以看出欧拉定理是费马小定理的推广)

费马小定理要求模数p为素数，欧拉定理则没有此要求，但是似乎还有个问题？如何判断a是否有逆元呢?

再求一次gcd判断是否互质吗？这还不如直接用扩展欧几里得算法呢。

可以由逆元性质直接检验是否为逆元，看求出的值x与a相乘模p是否为1即可。

算法时间复杂度：O(logn)

模板代码：

#include <iostream>
using namespace std;
typedef long long LL;
LL euler(LL n) {  // 欧拉函数
    LL res = n, i;
    for (i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            res = res / i * (i - 1);
            while (n % i == 0) n /= i;
        }
    }
    if (n != 1) res = res / n * (n - 1);
    return res;
}
LL pow_mod(LL a, LL b, LL p) {  // 快速幂模
    LL ret = 1;
    while (b) {
        if(b & 1) ret = (ret * a) % p;
        a = (a * a) % p;
        b >>= 1;
    }
    return ret;
}

int main() {
    LL a, p, inv;
    cin >> a >> p;
    inv = pow_mod(a, euler(p) - 1, p); // inv为a关于模p的逆元
    if (a * inv % p == 1) cout << inv << endl;  // 由逆元性质检验逆元是否存在
    else cout << -1 << endl;  // 不存在输出-1

    return 0;
}

 

4. O(n)求1~n逆元表

在模质数p下，求1~n逆元(n < p)

inv(a) = (p - p / a) * inv(p % a) % p

证明：

设x = p % a,y = p / a
于是有 x + y * a = p
(x + y * a) % p = 0
移项得 x % p = (-y) * a % p
x * inv(a) % p = (-y) % p
inv(a) = (p - y) * inv(x) % p
于是 inv(a) = (p - p / a) * inv(p % a) % p

模板代码：

#include<cstdio>
const int N = 200000 + 5;
const int MOD = (int)1e9 + 7;
int inv[N];
void init(){
    inv[1] = 1;
    for(int i = 2; i < N; i ++){
        inv[i] = (MOD - MOD / i) * 1ll * inv[MOD % i] % MOD;
    }
}
int main(){
    init();
}

 

阶乘逆元求法

fac[0]=fac[1]=1;

for(int i=2;i<=MAXN;i++)fac[i]=fac[i-1]*i%mod;

inv[MAXN]=qpow(fac[MAXN],mod-2);

for(int i=MAXN-1;i>=0;i--)inv[i]=inv[i+1]*(i+1)%mod;

 

inv(x) 是积性函数

证明：设 x 是 a 关于 p 的逆元，y 是 b 关于 p 的逆元，即 xa mod p = yb mod p = 1，则

xayb mod p = 1

(ab)*(xy) mod p = 1

即 xy 是 ab 关于模 p 的逆元，即 inv(ab) = xy = inv(a)*inv(b)。

 

总结：

1. 逆元求解一般利用扩欧。

2. 当模数p为质数的时候直接使用费马小定理。

3. p非质数使用欧拉函数(一般不用)。

 

HDU 1576 -- A/B   AC代码：

#include <iostream>
#define MOD 9973
using namespace std;
typedef long long LL;
LL qmod(LL a, LL b) {   // 快速幂模
    LL res = 1;
    while (b) {
        if (b & 1) res = res * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return res;
}

int main() {
    int t; cin >> t;
    while (t--) {
        int n, b; cin >> n >> b;
        cout << 1ll * n * qmod(1ll * b, MOD - 2) % MOD << endl;  // 费马小定理求b的逆元
    }
    return 0;
}

进阶题：HDU 5976 -- Detachment (贪心+逆元表+前缀和、积)：

#include <stdio.h>
const int maxn = 1e5 + 10;
const int MOD = 1e9 + 7;
typedef long long LL;
LL mul[maxn], sum[maxn], inv[maxn];
void init() {
    mul[1] = 1;
    sum[1] = 0;
    inv[1] = 1;
    for (int i = 2; i < maxn; i++) {
        sum[i] = sum[i-1] + i;
        mul[i] = (i * mul[i-1]) % MOD;
        inv[i] = (MOD - MOD / i) * 1ll * inv[MOD % i] % MOD;
    }
}

int main() {
    int t, x;
    init();
    scanf("%d", &t);
    while (t--) {
        scanf("%d", &x);
        if (x == 1) { puts("1"); continue; }
        int l = 2, r = maxn, mid, idx;
        while (l <= r) {
            mid = (l + r) / 2;
            if (sum[mid] <= x) idx = mid, l = mid + 1;
            else r = mid - 1;
        }
        int rest = x - sum[idx];
        LL ans = 0;
        if (rest == idx) ans = (mul[idx] * inv[2] % MOD * (idx + 2)) % MOD;
        else ans = mul[idx+1] * inv[idx+1-rest] % MOD;
        printf("%I64d\n", ans);
    }
    return 0;
}