luogu P5345 【XR-1】快乐肥宅

题目链接

观察形如 (k^xequiv rpmod g) 的式子，发现 (k^x) 在模 (g) 意义下一定有循环节，并且一定是一个 ( ho)。定义除环上点外其余点在“尾巴”上。

先考虑一个数 (i) 如果在环上会有什么特征（(d=gcd(i,g))）：

[i imes k^pequiv i pmod g\ frac i d imes k^pequiv frac i dpmod {frac g d}\ k^pequiv 1pmod {frac g d} ]

也就是 ((k,frac g d)=1)。据此可以得出“尾巴”的长度最多为 (log)。

暴力找到第一个在环上的点 (c)，假设 (c=k^t)，那么原式可以写为（(d=(c,g))）：

[k^xequiv rpmod g\ c imes k^{x-t}equiv rpmod g\ k^{x-t}equiv frac r d left(frac c d ight)^{-1}pmod {frac g d} ]

则可以用 BSGS 算出 (y=x-t) 以及 (k) 对 (frac g d) 的阶（假设为 (o)）。现在一个方程就被转化成了 (xequiv ypmod o)（(xgeq y)）。

接下来的问题就是合并这些方程。考虑直接用 excrt 合并，过程中如果模数 (ge 10^9) 就直接跳出并判断当前解是否对所以方程成立即可。

时间复杂度 (O(nsqrt g))。

代码：

#include<bits/stdc++.h>
#define int long long

using namespace std;

const int N = 1009;
int n, k[N], g[N], r[N], a[N], m[N], Max;
int o[N];
struct my_hash {
    static uint64_t splitmix64(uint64_t x) {
        x += 0x9e3779b97f4a7c15;
        x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
        x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
        return x ^ (x >> 31);
    }

    size_t operator()(uint64_t x) const {
        static const uint64_t FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
        return splitmix64(x + FIXED_RANDOM);
    }
};

unordered_map <int, int, my_hash> Map;

int gsc(int a, int b, int M)
{
	int res = 0;
	while (b)
	{
		if (b & 1)	
			res = (res + a) % M;
		b >>= 1, a = (a + a) % M;
	}
	return res;
}

int ksm(int a, int b, int M)
{
	int res = 1;
	while (b)
	{
		if (b & 1)
			res = res * a % M;
		b >>= 1, a = a * a % M;
	}
	return res;
}

void Exgcd(int a, int b, int &d, int &x, int &y)
{
	if (b == 0)
	{
		x = 1, y = 0, d = a;
		return;
	}
	Exgcd(b, a % b, d, x, y);
	int t = x;
	x = y, y = t - a / b * y;
}

int invs(int x, int M)
{
	int d, X, Y;
	Exgcd(x, M, d, X, Y);
	return (X % M + M) % M;
}

void init()
{
	scanf("%lld", &n);
	for (int i = 1; i <= n; i++)
		scanf("%lld %lld %lld", &k[i], &g[i], &r[i]), r[i] %= g[i];
}

int BSGS(int a, int b, int M)
{
	a %= M, b %= M;
	if (b == 1 || M == 1) return 0;
	Map.clear();
	int Sqr = ceil(sqrt(M)), j = 1;
	for (int i = 0; i < Sqr; i++)
		Map[b * j % M] = i, j = j * a % M;
	for (int now = j, i = 1; i <= Sqr; i++, now = now * j % M)
		if (Map.find(now) != Map.end())
			return i * Sqr - Map[now];
	return -1;
}

int _BSGS(int a, int b, int M)
{
	a %= M, b %= M;
	Map.clear();
	int Sqr = ceil(sqrt(M)), j = 1;
	for (int i = 0; i < Sqr; i++)
		Map[b * j % M] = i, j = j * a % M;
	for (int now = j, i = 1; i <= Sqr; i++, now = now * j % M)
		if (Map.find(now) != Map.end())
			return i * Sqr - Map[now];
	return -1;
}

int ExBSGS(int a, int b, int M)
{
	a %= M, b %= M;
	for (int i = 0; ; i++)
	{
		if (b == 1) return i;
		int d = __gcd(a, M);
		if (b % d != 0) return -1;
		if (d == 1)
		{
			int k = BSGS(a, b, M);
			if (k == -1) return -1;
			return k + i;
		}
		M /= d;
		if (M == 1) return i + 1;
		b = b / d * invs(a / d, M) % M;
	}
	return 0;
}

bool Check(int a[], int m[], int A)
{
	for (int i = 1; i <= n; i++)
		if (A % m[i] != a[i] % m[i])
			return 0;
	return 1;
}

int ExCRT(int a[], int m[])
{
	int lcm = m[1], A = a[1];
	for (int i = 2; i <= n; i++)
	{
		if (lcm >= 1e9)
		{
			if (A <= 1e9 && Check(a, m, A)) return A;
			return -1;
		}
		int X0, Y0, d, L = lcm, c = ((a[i] - A) % m[i] + m[i]) % m[i];
		Exgcd(lcm, m[i], d, X0, Y0);
		if (c % d != 0)
			return -1;
		lcm = lcm / __gcd(lcm, m[i]) * m[i];
		A = ((gsc(gsc(X0, c / d, m[i] / d), L, lcm) + A) % lcm + lcm) % lcm;
	}
	int res = (A + lcm) % lcm;
	if (res < Max)
		res = ((Max - res - 1) / lcm + 1) * lcm + res;
	return res;
}

void check_(int x)
{
	int flag = 1;
	for (int i = 1; i <= n; i++)
		if (ksm(k[i], x, g[i]) % g[i] != r[i] % g[i])
		{
			flag = 0;
			break;
		}
	if (flag) printf("%lld
", x), exit(0);
}

void check(int x)
{
	int flag = 1;
	for (int i = 1; i <= n; i++)
		if (ksm(k[i], x, g[i]) != r[i] % g[i])
		{
			flag = 0;
			break;
		}
	if (flag) printf("%lld
", x);
	else puts("Impossible");
	exit(0);
}

void work()
{
	for (int i = 1; i <= n; i++)
		if (ExBSGS(k[i], r[i], g[i]) == -1)
		{
			puts("Impossible");
			return;
		}
	for (int i = 1; i <= n; i++)
		if (__gcd(r[i], g[i]) != __gcd(r[i] * k[i], g[i]))
			check(ExBSGS(k[i], r[i], g[i]));
	for (int i = 1; i <= n; i++)
	{
		int X = 1, o = 0;
		while (__gcd(X, g[i]) != __gcd(X * k[i], g[i]))
			X = X * k[i] % g[i], o++;
		int d = __gcd(X, g[i]);
		if (r[i] % d != 0)
		{
			puts("Impossible");
			exit(0);
		}
		r[i] /= d, g[i] /= d;
		r[i] = r[i] * invs(X / d, g[i]) % g[i];
		m[i] = _BSGS(k[i], 1, g[i]);
		a[i] = BSGS(k[i], r[i], g[i]) + o;
		Max = max(Max, a[i]);
	}
	int ans = ExCRT(a, m);
	if (ans != -1)
		printf("%lld
", ans);
	else
		puts("Impossible");
}

signed main()
{
	init();
	work();
	return 0;
}