【题解】Luogu-P5345 【XR-1】快乐肥宅

P5345 【XR-1】快乐肥宅

Preface

好题！！！

当然也很毒瘤。

突然感到屠龙勇士在这题面前

就是逊内！！！

Description

给定高次同余方程组

\[\begin{cases} k_1^x\equiv r_1\pmod {g_1}\\ k_2^x\equiv r_2\pmod {g_2}\\ \cdots\\ k_n^x\equiv r_n\pmod {g_n} \end{cases} \]

请求出该方程组的最小 非负 整数解，若 无解 或 最小解 \(\ge 10^9\)，输出 Impossible。

  • 对于 \(100\%\) 的数据，\(1 \le n \le 10^3\)，\(1 \le k_i, r_i \le g_i \le 10^7\)。

Solution

前置芝士：

  • exBSGS
  • exCRT

一、求出单个同余方程的解

对于 \(a^x\bmod p\)，结果应该分为一段“尾巴”和一段循环部分，例如 \(124^x\bmod 495616\)：

感谢兔队的图

这时对于一个正整数 \(b\)，\(a^x\equiv b\pmod p\) 的解有 \(3\) 种情况：

  • 唯一解：即 \(x\) 在尾巴上。例如 \(b=15376\) 时，\(x=2\)。
  • 无穷解：即 \(x\) 在循环上。例如 \(b=241664\) 时，\(x\equiv 3+5\pmod 5\) 且 \(x\ge 3+5\)，其中 \(3\) 表示 \(b\) 是循环中的第 \(3\) 个，\(+5\) 表示尾巴的长度，模数为 \(5\) 表示循环的长度，不等式是为了保证 \(x\) 在循环上。
  • 无解：例如 \(b=114514\) 时。

在 exBSGS 中，除到 \(\gcd(a,p)=1\) 时，先用 BSGS 求出 \(a^x\equiv b\pmod p\) 的最小非负整数解，再求出 \(a\) 模 \(p\) 的阶（\(\gcd(a,p)=1\) 时阶一定存在），那么阶就是循环长度。

下面的代码返回一个 \(\operatorname{pair}\)，第一个数表示最小解（如果无解则为 \(-1\)），第二个数表示阶（如果不存在则为 \(-1\)）。

与 exBSGS 模板一样，需要特判 \(b=1\) 或 \(p=1\) 的情况。

pair<int, int> exBSGS(int a, int b, int p)
{
	if (p == 1)
	{
		return make_pair(0, 1);
	}
	a %= p, b %= p;
	int g = exgcd(a, p), f = 1, k = 0;
	while (g > 1)
	{
		if (b % g)
		{
			if (b == 1)
			{
				return make_pair(0, -1);
			}
			return make_pair(-1, -1);
		}
		k++;
		b /= g, p /= g;
		f = (ll)f * (a / g) % p;
		g = exgcd(a, p);
		if (f == b)
		{
			if (g > 1)
			{
				return make_pair(k, -1);
			}
			int y = BSGS(a, 1, p);
			return make_pair(k, y);
		}
	}
	int x = BSGS(a, (ll)b * inv(f, p) % p, p);
	if (x == -1)
	{
		return make_pair(-1, -1);
	}
	int y = BSGS(a, 1, p);
	return make_pair(x % y + k, y);
}

在主函数中：

pair<int, int> res = exBSGS(k, r, g);
b[i] = res.first, a[i] = res.second;
if (b[i] == -1) // x 无解就肯定无解
{
    puts("Impossible");
    return 0;
}
else if (a[i] == -1) // 阶不存在，用一个 X 记录唯一的解
{
    X = b[i];
    continue;
}
mn = max(mn, b[i]); // x ≥ mn

二、合并解

如果出现了唯一解（即上面代码中的 \(X\ne -1\)）的情况，直接将这个 \(X\) 带入剩下的方程中检验即可。

否则直接用 exCRT 合并即可，注意要加上 \(\operatorname{lcm}(a_1,a_2,\dots,a_n)\) 的倍数使得答案不小于上面代码中的 \(mn\)。

我们设前 \((i-1)\) 个方程组合并后的解为 \(x\equiv B\pmod A\)，第 \(i\) 个方程为 \(x\equiv b\pmod a\)。

问题来了， \(10^3\) 个 \(10^7\)，合并后 \(A,B\) 都有可能爆 long long。

观察题目中的要求：

若无解或 最小解 \(\ge 10^9\)，输出 Impossible。

当 \(A>10^9\) 时，直接检测 \(B\) 是否满足 \(B\equiv b\pmod a\) 即可，如果不满足那么一定要加上若干个 \(A\)，那么一定超过了 \(10^9\)，不满足要求。

全部合并完之后记得也要判断 \(B\) 是否在 \(10^9\) 内。

Code

//18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#include <map>
#include <cmath>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;

ll x, y;

ll exgcd(ll a, ll b)
{
	if (!b)
	{
		x = 1, y = 0;
		return a;
	}
	ll Gcd = exgcd(b, a % b);
	ll tmp = x;
	x = y;
	y = tmp - a / b * y;
	return Gcd;
}

ll A, B;

int merge(ll a, ll b)
{
	ll Gcd = exgcd(A, a), c = B - b;
	if (c % Gcd)
	{
		return -1;
	}
	c /= Gcd;
	x = (x * c % a + a) % a;
	ll Lcm = A / Gcd * a;
	B = (B - A * x % Lcm + Lcm) % Lcm;
	A = Lcm;
	return 1;
}

int div_ceil(int a, int b)
{
	return (a - 1) / b + 1;
}

const int MAXN = 1e3 + 5;
const int D = 1e9;

int a[MAXN], b[MAXN];

int exCRT(int n, int mn)
{
	A = a[1], B = b[1];
	for (int i = 2; i <= n; i++)
	{
		if (A > D && (B - b[i]) % a[i])
		{
			return -1;
		}
		else if (A <= D && merge(a[i], b[i]) == -1)
		{
			return -1;
		}
	}
	if (B < mn)
	{
		B += div_ceil(mn - B, A) * A;
	}
	if (B > D)
	{
		return -1;
	}
	return B;
}

int inv(int a, int p)
{
	exgcd(a, p);
	x = (x % p + p) % p;
	return x;
}

int phi(int p)
{
	int ans = p;
	for (int i = 2; i * i <= p; i++)
	{
		if (p % i == 0)
		{
			ans = ans / i * (i - 1);
			while (p % i == 0)
			{
				p /= i;
			}
		}
	}
	if (p != 1)
	{
		ans = ans / p * (p - 1);
	}
	return ans;
}

int BSGS(int a, int b, int p)
{
	map<int, int> hash;
	int t = ceil(sqrt(phi(p))), aj = 1;
	for (int j = 0; j < t; j++)
	{
		int val = (ll)b * aj % p;
		hash[val] = j;
		aj = (ll)aj * a % p;
	}
	a = aj;
	int ai = 1;
	for (int i = 1; i <= t; i++)
	{
		ai = (ll)ai * a % p;
		if (hash.find(ai) != hash.end())
		{
			int j = hash[ai];
			if (i * t - j >= 0)
			{
				return i * t - j;
			}
		}
	}
	return -1;
}

pair<int, int> exBSGS(int a, int b, int p)
{
	if (p == 1)
	{
		return make_pair(0, 1);
	}
	a %= p, b %= p;
	int g = exgcd(a, p), f = 1, k = 0;
	while (g > 1)
	{
		if (b % g)
		{
			if (b == 1)
			{
				return make_pair(0, -1);
			}
			return make_pair(-1, -1);
		}
		k++;
		b /= g, p /= g;
		f = (ll)f * (a / g) % p;
		g = exgcd(a, p);
		if (f == b)
		{
			if (g > 1)
			{
				return make_pair(k, -1);
			}
			int y = BSGS(a, 1, p);
			return make_pair(k, y);
		}
	}
	int x = BSGS(a, (ll)b * inv(f, p) % p, p);
	if (x == -1)
	{
		return make_pair(-1, -1);
	}
	int y = BSGS(a, 1, p);
	return make_pair(x % y + k, y);
}

int main()
{
	int n;
	scanf("%d", &n);
	int mn = 0, X = -1;
	for (int i = 1; i <= n; i++)
	{
		int k, g, r;
		scanf("%d%d%d", &k, &g, &r);
		pair<int, int> res = exBSGS(k, r, g);
		b[i] = res.first, a[i] = res.second;
		if (b[i] == -1)
		{
			puts("Impossible");
			return 0;
		}
		else if (a[i] == -1)
		{
			X = b[i];
			continue;
		}
		mn = max(mn, b[i]);
	}
	if (X != -1)
	{
		for (int i = 1; i <= n; i++)
		{
			if ((a[i] == -1 && X != b[i]) || (a[i] != -1 && (X < b[i] || (X - b[i]) % a[i])))
			{
				puts("Impossible");
				return 0;
			}
		}
		printf("%d\n", X);
		return 0;
	}
	int ans = exCRT(n, mn);
	if (ans == -1)
	{
		puts("Impossible");
	}
	else
	{
		printf("%d\n", ans);
	}
	return 0;
}

References

posted @ 2021-12-11 10:36  mango09  阅读(60)  评论(0)  编辑  收藏  举报
-->