2017-03-26

bzoj

[bzoj 2219] 数论之神

给定正整数A, B, K, 求在上的解的个数, T组数据. (T ≤ 1000, 1 ≤ A, B ≤ 10^9, 1 ≤ K ≤ 5* 10^8) 我的做法好像不太主流, 后面会讲解两种方法.

背景知识

中国剩余定理

设, 两两互质, 则中的元素与中的元素一一对应.

原根

拉格朗日定理 设是有限群的子群, 则.

有限群中, 一个元素不断搞自己会生成一个子群, 称是的生成元, 满足的最小正整数是在中的阶, 记作. 可以证明, 且. 结合拉格朗日定理, 有.

下面考虑模乘法群. 首先, 由上面子群有关性质, 可以导出欧拉定理.

若, 则是整个群的生成元, 又称原根. 对所有奇素数和所有正整数, 有原根当且仅当 $Missing \left or extra \rightn \in \left\\{2, 4, p^e, 2p^e \right\\}$ .

不是原根 <=> 是的真约数 <=> 存在的质因子使得 <=> 存在的质因子使得. 由此, 得到一个不那么暴力的检验原根的算法. 从小到大枚举, 检验之, 即可求得一个原根. 据说由于原根一般不会太大所以这样就可以了......

离散对数

当, 可以使用大步小步算法解这个关于的方程.

设, , 方程转化为

使用meet-in-middle的策略, 枚举, 在哈希表中储存方程右边的所有取值, 再枚举, 到哈希表中查询. 平均时间复杂度.

同余式的一个性质

若, 则的公约数也是的约数. 再根据取模运算的分配律, 可以做除法: .

模线性方程

设. 是这个数某个排列的份复制. 于是上述方程有解当且仅当, 且如果有解, 则在上有个解.

两种做法公用的部分

首先将做质因数分解: , 则 (中国剩余定理). 考虑求这样的方程的解个数: , 是奇素数 (所以存在原根). 最后把答案乘起来即可.

先看一种特殊情形: , 存在原根, 互质.

求这个方程在上解的个数至少有两种方法.

取的一个原根, 设 (是唯一的), 则上述方程转化为, 对大步小步算法稍加改写, 以求得解的数目. 由于这里规定了的范围, 所以模的最后一个周期直接枚举.
如果和不互质, 则无解. 还是取原根, 设, 方程转化为. 用大步小步算法求出. 在上的解与原方程的解一一对应. 根据模线性方程相关知识, 该方程有解当且仅当, 如果有解, 解的数目是.

做法一

来自我自己.

若与互质, 用之前所述方法求之.

若与不互质, 则. 设, 分类讨论. - . 方程左边等于0. 若, 则有个解, 否则无解. - . 根据上面所说的同余式的一个性质, , 如果不满足, 则无解; 否则, 方程的左右两边和模同时约去: , 这个问题和原问题有相同的形式, 可递归求解. 但是取值范围变了. 递归下去, . 根据周期性, 只需将答案乘上.

问题顺利解决. 写代码的时候, 我把在最前面单独讨论. , 于是是范围内所有含因子的数, 共个.

做法二

参考 - Po姐的题解 - 【BZOJ 2219】【超详细题解】数论之神 Regina8023.

思路大体相同, 中国剩余定理.

做法一讨论了与互质, 与不互质两种情况. 事实上, 与的关系已经决定了与是否互质.

若, 直接返回.

若, 设, 则. 根据同余式的一个性质, 如果有解, 则, 并且由于右边不含, 左边也不能含, 故方程等价于其中与互质.

和做法一中的递归类似, 这个方程解的个数要乘上才是我们要的答案. 用公用部分所述方法求之.

感想

这道题的综合性比较强, 涉及了不少初等数论的知识. 细心分类讨论, 关注定义域的变化.

#include <bits/stdc++.h>
#define Rep(i, a, b) for (int i = (a), i##_ = (b); i < i##_; ++i)
#define For(i, a, b) for (int i = (a), i##_ = (b); i <= i##_; ++i)
#define Down(i, a, b) for (int i = (a), i##_ = (b); i >= i##_; --i)

using namespace std;
typedef long long ll;
typedef pair<int, int> ii;
typedef vector<ii> Vec;
typedef Vec::iterator iter;

struct Hash {
	const static int M = 38321;
	vector<ii> h[M];
	int cnt, used[M];

	void insert(int v)
	{
		int k = v % M;
		if (h[k].empty()) used[cnt++] = k;
		for (iter x = h[k].begin(); x != h[k].end(); ++x)
			if (x->first == v) {
				++x->second;
				return;
			}
		h[k].push_back(ii(v, 1));
	}

	int operator[](int v)
	{
		int k = v % M;
		for (iter x = h[k].begin(); x != h[k].end(); ++x)
			if (x->first == v)
				return x->second;
		return 0;
	}
	
	void clear()
	{
		while (cnt) h[used[--cnt]].clear();
	}
} H;

void decompose(int n, Vec& v)
{
	for (int d = 2; d*d <= n; ++d) {
		if (n % d == 0) {
			ii t(1, d);
			do {
				t.first *= d;
				n /= d;
			} while (n % d == 0);
			v.push_back(t);
		}
	}
	if (n > 1)
		v.push_back(ii(n, n));
}

ll fp(ll x, ll n)
{
	ll y = 1;
	for (; n; n >>= 1, x *= x)
		if (n & 1)
			y *= x;
	return y;
}

ll fpm(ll x, ll n, ll m)
{
	ll y = 1;
	for (; n; n >>= 1, (x *= x) %= m)
		if (n & 1)
			(y *= x) %= m;
	return y;
}

int bsgs(ll a, ll y, ll n, ll r)
{
	H.clear();
	int m = round(sqrt(r)), cnt = 0;
	ll t = y;
	For (i, 1, m) // y*a^i, i=1,2,...,m
		H.insert((t *= a) %= n);
	ll d = fpm(a, m, n);
	t = 1;
	For (i, 1, r/m) // (a^m)^i, i=1,2,...,floor(r/m)
		cnt += H[(t *= d) %= n];
	Rep (i, r/m*m, r) // a^i, i=m*floor(r/m),...,r-1
		cnt += t == y, (t *= a) %= n;
	return cnt;
}

int primitive_root(int a, int p, int r)
{
	Vec v;
	decompose(r, v);
	Rep (g, 2, a) if (g % p) {
		bool ok = true;
		for (iter x = v.begin(); x != v.end(); ++x)
			if (fpm(g, r/x->second, a) == 1) {
				ok = false;
				break;
			}
		if (ok) return g;
	}
	assert(0);
	return -1;
}

int solve(int a, int b, int m, int p, int k)
{
	if (!b) return fp(p, k-(k+a-1)/a);
	int phi = m/p*(p-1), g = primitive_root(m, p, phi), ans = bsgs(fpm(g, a%phi, m), b, m, phi);
	if (a < k) {
		int t = fp(p, a), m1 = m / t;
		if (b % t == 0)
			ans += t / p * solve(a, b/t%m1, m1, p, k-a);
	}
	return ans;
}

int main()
{
	int T, K, A, B;
	scanf("%d", &T);
	while (T--) {
		scanf("%d%d%d", &A, &B, &K);
		Vec v;
		decompose(2*K + 1, v);
		int ans = 1;
		for (iter x = v.begin(); x != v.end(); ++x) {
			int m = x->first, p = x->second, k = 1;
			for (int t = p; t < m; t *= p) ++k;
			ans *= solve(A, B % m, m, p, k);
		}
		printf("%d\n", ans);
	}
	return 0;
}