[NOIP2012 提高组] 同余方程

[NOIP2012 提高组] 同余方程

求关于 x 的同余方程 ax ≡ 1 (mod b) 的最小正整数解.
(a, b为正整数)

首先我们立刻就能想到三种方法（bushi

1.扩展欧几里得

我们知道，扩欧实际上是求解方程 ax + by = gcd(a, b) 的过程，那么这两个方程之间有什么关系呢？
因为 by 是 b 的倍数，所以不难证明:

ax + by ≡ 1 (mod b) ①

加了by跟没加一样
我们设 m = ax + by，因为 a 是 gcd(a, b) 的倍数，b 也是 gcd(a, b) 的倍数，
所以 ax，by 也都是 gcd(a, b)的倍数，那么 (ax + by) mod gcd(a, b) = 0.
即 m mod gcd(a, b) = 0.

因为在方程中 m = 1，所以 1 mod gcd(a, b) = 0，即 a 与 b互质.
因此我们只要求解方程 ax + by = gcd(a, b) 就OK啦（大大滴OK

直接上代码:

#include <iostream>
using namespace std;
typedef long long ll;

void exgcd(ll a, ll b, ll& x, ll& y) {
	if(b == 0) {x = 1; y = 0;}
	else {exgcd(b, a%b, y, x); y -= (a/b)*x;}
}
int main() {
	ll a, b, x, y;
	cin >> a >> b;
	exgcd(a, b, x, y);
	cout << (x + b) % b; //避免出现 x<0 的情况
	return 0;
}

对于这道题来讲，扩展欧几里得的时间复杂度应该是最低的，为 O(log n).
实际测试中，扩欧是这道题的正解.

2.费马小定理

首先给出费马小定理：

若 p 为素数，且 a 与 p 互质，则有
a^{p-1} ≡ 1 (mod p)

而题中的 ax ≡ 1 (mod b) 实际上可以与其合并：
∵a^{b-1} ≡ 1 (mod b), ax ≡ 1 (mod b)
∴a^{b-1} ≡ ax (mod b)
∴x ≡ a^{b-2} (mod b)
所以 x 的值只需要求出 a^{b-2} mod b 的值即可，时间复杂度也是 O(log n) 的（快速幂），但是这种方法只适用于 b 为素数的情况，因此只能拿到部分分数.

代码：

#include <iostream>
using namespace std;
typedef long long ll;
ll fpm(ll a, ll b, ll mod) {
	a %= mod;
	ll res = 1;
	while(b) {
		if(b&1) res = (res * a) % mod;
		a = (a*a) % mod;
		b >>= 1;
	}
	return res;
}
int main() {
	ll a, b;
	cin >> a >> b;
	cout << fpm(a, b-2, b);
	return 0;
}

实际测试中，这种写法只能拿到可怜的 20 分.

3.欧拉函数

欧拉定理：若正整数 a, p 互质，则有：
a^φ(n) ≡ 1 (mod n)

同费马小定理的解法，我们可以证明：
∵a^φ(b) ≡ 1 (mod b), ax ≡ 1 (mod b)
∴a^φ(b) ≡ ax (mod b)
∴x ≡ a^{φ(b)-1} (mod b)
所以可以求出 a^{φ(b)-1} mod b 的值来求出 x.
这里我们写一下欧拉函数：

ll phi(ll x) {
	ll res = x;
	for(int i = 2; i*i <= x; i++) {
		if(x % i == 0) res = res * (i-1) / i;
		while(x % i == 0) x /= i;
	}
	if(x > 0) res = res * (x-1) / x;
	return res;
}

整体来讲，欧拉函数的时间复杂度是 O(√n) 的.
所以最终时间复杂度应该是 O(√n + log n) （欧拉函数+快速幂）.

#include <iostream>
using namespace std;
typedef long long ll;
ll phi(ll x) {
	ll res = x;
	for(int i = 2; i*i <= x; i++) {
		if(x%i == 0) res = res * (i-1) / i;
		while(x%i == 0) x /= i;
	}
	if(x > 0) res = res * (x-1) / x;
	return res;
}
ll fpm(ll a, ll b, ll mod) {
	a %= mod;
    ll res = 1;
    while(b) {
        if(b&1) res = (res * a) % mod;
        a = (a * a) % mod;
        b >>= 1;
    }
    return (res + b) % mod;
}
int main() {
	ll a, b;
	cin >> a >> b;
	cout << fpm(a, phi(b)-1, b);
	return 0;
}

而很遗憾，这种算法的最终分数能拿到 70 分，所以不推荐.

所以说欧几里得 yyds 啊
点赞支持，感激不尽（磕头磕头磕头