同余方程 【代码】

此文包含：

• 线性同余方程
• 线性同余方程组
• 高次同余方程

详见代码注释。

/*
同余方程系列三合一 

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;
	}
}