同餘方程 【代碼】

此文包含：

• 線性同餘方程
• 線性同餘方程組
• 高次同餘方程

詳見代碼註釋。

/*
同餘方程系列三合一 

2021.01.02 - 2021.xx.xx

Author: zimujun 
*/

/*
Content

exgcd

線性同餘方程

線性同餘方程組 

高次同餘方程
*/

namespace Basic {
	template <typename Temp> inline void read(Temp & res) {
		Temp fh = 1; res = 0; char ch = getchar();
		for(; !isdigit(ch); ch = getchar()) if(ch == '-') fh = -1;
		for(; isdigit(ch); ch = getchar()) res = (res << 3) + (res << 1) + (ch ^ '0');
		res = res * fh;
	}
	template <typename Temp> inline void Checkmax(Temp & num, Temp comp) {if(comp > num) num = comp;}
	template <typename Temp> inline void Checkmin(Temp & num, Temp comp) {if(comp < num) num = comp;}
}
using namespace Basic;

namespace Math {
	template <typename Temp> inline Temp gcd_(Temp A, Temp B) {Temp C; while(B) C = B, B = A % B, A = C; return A;}
	template <typename Temp> inline Temp lcm_(Temp A, Temp B) {return A / gcd_<Temp>(A, B) * B;}
	/*
	該 GCD 函數是解決同餘方程問題的一個有用工具。 
	*/
	template <typename Temp> Temp gcd(Temp A, Temp B, Temp & X, Temp & Y) {
		if(B == 0) {X = 1; Y = 0; return A;}
		LL res = gcd(B, A % B, X, Y);
		LL Z = Y; Y = X - A / B * Y; X = Z;
		return res;
	}
	template <typename Temp> inline Temp qpow(Temp A, Temp P, Temp Mod) {Temp res = 1; while(P) res *= (P & 1 ? A : 1), res %= Mod, A = A * A % Mod, P >>= 1; return res;}
}

namespace equation {
	/*
	Description: 這是一個用於解決線性同餘方程問題的函數 
	Algorithm used: 歐幾里得算法 
	Step:
		1. 將同餘方程 ax = b (mod m) 轉化爲二元一次不定方程 ax + (-y)m = b； 
		2. 使用歐幾里得算法求出不定方程 ax + (-y)m = gcd(a, m) 的一個特解 x_0； 
		3. 令 x = x_0 * (b / gcd(a, m)) 那麼 x 就是該同餘方程的一個特解。通解爲 x + tm。 
	*/
	template <typename Temp> inline void solve(Temp A, Temp B, Temp M, Temp & X) {
		Temp D = Math::gcd_<Temp>(A, M);
		if(B % D) {X = -1; return;}
		Temp Y; D = Math::gcd<Temp>(A, M, X, Y);
		X = B / D * X;
		X = (X % M + M) % M;
	}
}

namespace equationset {
	/*
	Description: 這是一個用來解線性同餘方程組的函數。 
	Algorithms used: 中國剩餘定理及擴展（ex CRT） 
	Steps:
		1. 求出第一個線性同餘方程的解； 
		2. 計算前 k 個同餘方程的模數的 lcm，設爲 M； 
		3. 設前 k 個同餘方程的解爲 x, 那麼 (x + tM) 就是這些方程構成的方程組的通解。代入第 k + 1 個方程，整理獲得另外一個同餘方程並求解； 
		4. 重複 2，3 操做。 
	*/
	template <typename Temp> inline void solve(int number, Temp A[], Temp B[], Temp M[], Temp & X) {
		if(number >= 1) equation::solve(A[1], B[1], M[1], X); Temp Ms = M[1]; if(X == -1) return;
		for(register int i = 2; i <= number; ++i) {
			Temp tA = A[i] * Ms, tB = B[i] - A[i] * X, tM = M[i];
			tB = (tB % tM + tM) % tM; Temp resT;
			equation::solve(tA, tB, tM, resT); 
			if(resT == -1) return; X += resT * Ms;
			Ms = Math::lcm_<Temp>(Ms, M[i]);
			X = (X % Ms + Ms) % Ms; 
		}
	}
}

namespace HDCE {
	/*
	Description: 這是一個用來解決高次同餘方程的函數。 
	Algorithms used:BSGS(Baby_Step_Giant_Step)
	Step:
		設原方程爲 a^x = b (mod m)
		1. 記 t = ceil(sqrt(m))， 將 x 表示爲 it - j 的形式，那麼原方程化爲 (a^t)^i = b * a^j (mod m)，其中 j in [0, t - 1]； 
		2. 計算 j = 0 ~ (t - 1) 時分別的 b * a^j mod m 的值， 存入一個哈希表。
		3. 枚舉 i, 計算 i = 0 ~ t 時 (a^t)^i 的值，從哈希表中找到答案。 
	*/
	template <typename Temp> inline void BSGS(Temp A, Temp B, Temp M, Temp & X) {
		std::map<Temp, Temp> idx; idx.clear(); B = (B % M + M) % M; Temp t = sqrt(M) + 1; Temp tB = B % M; X = 0x7fffffff;
		for(register Temp i = 0; i < t; ++i) {
			/*std::cout << tB << " ";*/ idx[tB] = i + 1, tB = tB * A % M;
		}
		Temp tA = Math::qpow(A, t, M), tT = 1;
		/*
		for(register int i = 0; i < M; ++i) std::cout << idx[i] - 1 << " ";
		puts("");
		
		std::cout << tA << std::endl;
		*/
		if(tA == 0) {X = B == 0 ? 1 : -1; return;}
		for(register Temp i = 0; i <= t; ++i) {
			Temp j = idx[tT];
			//printf("%d ^ %d ^ %d = %d    ", A, t, i, tT);
			if(j) {
				j--; Temp answer = i * t - j;
				if(answer < 0) continue;
				Checkmin(X, answer);
			}
			tT = tT * tA % M;
		}
		if(X == 0x7fffffff) X = -1;
	}
}