[BZOJ2219] 数论之神 [原根&BSGS&原根&指标]

[$\mathfrak{L}\mathfrak{i}\mathfrak{n}\mathfrak{k}$]

---

整数模 n 乘法群

在 ${\mathbb{Z}}_n$ 群中，根据定义 $|{\mathbb{Z}}_n|=\varphi (n)$ 。

---

阶

$a$ 关于 $p$ 的阶 ${\delta }_p(a)$ 是满足 $a^x\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 的最小的 $x$ 。
也就是 $|<a>|$
同时根据群论拉格朗日定理 $|G|=[G:H]|H|$ 我们有 $|H|||G|$
（群论拉朗易得只需证等价，等价只需证自反对称传递）
那么 $|<a>||\varphi (p)$ 并且 $a^{|G|}=1$

---

原根

模 $p$ 的原根 $g$ 满足 ${\delta }_p(g)=\varphi (p)$ 。
$p$ 的所有原根乘起来 $\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$
$p$ 的所有原根加起来 $\equiv \mu (p-1)(\mathrm{m}\mathrm{o}\mathrm{d}p)$

模 $p$ 的原根个数为 $\varphi (\varphi (p))$ 。（证明：$|{\mathbb{Z}}_p|=\varphi (p)$ 并且生成元个数为 $\varphi (|{\mathbb{Z}}_p|)$ ）
（好多地方都写的上面这个证明，但是我实在是没有找到关于 ${\mathbb{Z}}_p$ 可逆元的说法。。）
（我感觉貌似意思是说原根能生成 ${\mathbb{Z}}_p$ 然后利用所有原根能互相得到的性质？）
（接着把不能作为生成元的与 $p$ 互质的数 $g^{x^{\prime }}$ 去掉）
（然后能作为生成元的 $x$ 也就是应该满足 $x$ 的倍数模 $\varphi (p)$ 能取到所有的 $x$ ）
（这种 $x$ 有 $\varphi (\varphi (p))$ 个）
或者也可以利用 ${\delta }_p(g^x)=\frac{{\delta }_p(g)}{(x,{\delta }_p(g))}$ 即 ${\delta }_p(g^x)=\frac{\varphi (g)}{(x,\varphi (g))}$ 且 $g_{other}=g^x$
又 ${\delta }_p(g^x)=\varphi (p)$ 的时候 $(x,\varphi (g))=1$ 即有 $\varphi (\varphi (g))$ 个 $g^x$ 为原根。得证。
（如果 $g$ 是 $p$ 的原根那么 $g^x$ 是 $p$ 的原根的充要条件为 $(x,\varphi (p))=1$

$g$ 是模 $p$ 的原根的充要条件为 $g^x(\mathrm{m}\mathrm{o}\mathrm{d}p)$，$x=0,1,\cdots ,{\delta }_p(g)$ 两两不相等。
也就是 $g^x\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 最小正整数解为 $x=\varphi (p)$

有原根的数只有 $1,2,4,p^a,2p^a$ ，其中 $p$ 是奇素数。
模 $p$ 有原根，意味着 ${\mathbb{Z}}_p$ 是一个循环群。

求原根：唯一分解 $\varphi (p)=\prod {p_i}^{a_i}$ ， 暴枚 $x\in [2,p)$
对每个 $x$ 判断 $x^{\frac{\varphi (p)}{p_i}}$ 里是否有一个数 $\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$ ，有就不是原根
求到一个原根就跑 $O(g(p)\omega (\varphi (p))\log (\varphi (p)))$
其中 $g(p)$ 是 $p$ 的最小正原根，$\omega ()$ 是素因数个数函数（不重复计算）

$g(p)=O(p^{1/4+ϵ})$ ， $\omega (p-1)\le p(\sqrt{n})+1=O(\frac{\sqrt{n}}{\ln \sqrt{n}})$
其中 $p(n)$ 指不大于 $n$ 的素数个数
我们暂且无视那个不好搞的 $\varphi ()$
然后可以取一个玄学上界 $O(\frac{\sqrt{n}}{\ln \sqrt{n}}{p}^{1/4+ϵ}\log p)$

另外：上面提到的素数快速幂可以提出去
那么复杂度变成 O ( max ⁡ { g ( p ) ω ( p ) , ω ( p ) log ⁡ ( p ) } ) O(\max\{g(p)\omega(p),\omega(p)\log(p)\}) O(max{g(p)ω(p),ω(p)log(p)})
然后大概大概地估计一下差不多就是 $O(\frac{\sqrt{n}}{\ln \sqrt{n}}\log p)$ （？

以上纯属口胡 如果出锅 本人概不负责（
（如果你推出来的跟我推出来的不一样 我觉得你可以相信自己）

---

指标

$g^r\equiv a(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 中的 $r$ 是以 $g$ 为底 $a$ 的一个指标。其中 $g$ 是 $p$ 的一个原根。
实际上的记号一般采用 $ind$ 或者 $I$
显然 $I_x(a)\in [0,\varphi (p)](\mathrm{m}\mathrm{o}\mathrm{d}p)$ 且 $(a,p)=1$
我们有 $I_x(ab)\equiv I_x(a)+I_x(b)(\mathrm{m}\mathrm{o}\mathrm{d}\varphi (p))$
我们还有 $I_x(a^y)\equiv y{I}_x(a)(\mathrm{m}\mathrm{o}\mathrm{d}\varphi (p))$
可以理解指标为“离散对数”？需要特别注意的就是取指标后模数变了

---

BSGS

求 $g^r\equiv a(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 中的 $r$ ，其中 $p$ 不需要是素数，但需要有 $(g,p)=1$ 。
时间复杂度 $O(\sqrt{p})$ 。
$g$ 不需要是原根。

怎么求呢？设 $m=\lceil \sqrt{p}\rceil$，$r=im-j$ 其中 $i\in [1,m],j\in [0,m)$
那么 $$g^{im}\equiv a{g}^j(\mathrm{m}\mathrm{o}\mathrm{d}p),r=im-j$$ 把 $a{g}^j$ 存表，然后枚举 $i$ 查表。

存表的时候：如果需要让 $r$ 最小，更新表的时候就可以直接覆盖。

某种程度上这样的表算是可以起到约等于指标表的作用？

---

ExBSGS

虽然这道题用不到（

当 $p$ 不是素数时，设 $d=(x,p)$ 。
假如 $d$ 不整除 $y$ ，那么当 $y=1$ 的时候 $x=0$ ，否则无解。
假如整除就按照唯一分解不停搞掉，直到变成正常 bsgs 的形式

---

接下来考虑这道题。
好像就是求 exbsgs 解的个数。那考虑按照套路把模数唯一分解吧？
但是 exbsgs 是求解，这里是求方案数，？？
思考能不能分解成多个“求 bsgs 解的个数”的问题

$x^r\equiv a(\mathrm{m}\mathrm{o}\mathrm{d}2k+1)$ 里，范围内解 $x$ 的个数
容易发现相当于求同余方程组 { x r ≡ a ( m o d p i a i ) } \{x^r\equiv a\pmod{{p_i}^{a_i}}\} {xr≡a(modpi​ai​)} 解的个数
也就等于每个同余方程 $x^r\equiv a(\mathrm{m}\mathrm{o}\mathrm{d}{p_i}^{a_i})$ 把所有解组合一下？
所以答案等于每个同余方程 $x^r\equiv a(\mathrm{m}\mathrm{o}\mathrm{d}{p_i}^{a_i})$ 解的个数乘起来
（当然，这里谈到的同余方程的解取值一律为 $[0,mod)$

---

然后考虑求 $x^r\equiv a(\mathrm{m}\mathrm{o}\mathrm{d}{p_i}^{a_i})$ 解的个数
这个问题好像不必也不好继续化了？
想法：大力分类讨论
那首先挑一个简单地情况特判当场赐死他awsl启发一下？

如果 $a\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{p_i}^{a_i})$ 显然有 $x\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{p_i}^{a_i})$
显然只需把 $x$ 唯一分解，提出 $p_i$ 那一项，判断指数乘 $r$ 是否不小于 $a_i$ 即可
也就是说 $x$ 唯一分解后的 $p_i$ 一项的指数 $u$ 应该 $\ge \lceil \frac{a_i}{r}\rceil$
现在所有的 $x$ 都可以被表示为 $k{p_i}^u$ 也就是说它们的 $gcd$ 是 ${p_i}^u$
所以只需要求 $[0,{p_i}^{a_i})$ 内 ${p_i}^u$ 倍数的个数也就等于 ${p_i}^{a_i-u}$

否则怎么做？仍然有可能 $a={p_i}^t\cdot q$ 我们就讨论这个。
假如是，那么 $x^r\equiv q{p_i}^t(\mathrm{m}\mathrm{o}\mathrm{d}{p_i}^{a_i})$ 即 $x^r=(q+k{p_i}^{a_i-t}){p_i}^t$
有解需要 $r|t$ 不妨设 $t=yr$ ， $x=p^{y}b$ 那么 ${p_i}^{yr}{b}^r\equiv {p_i}^{yr}q(\mathrm{m}\mathrm{o}\mathrm{d}{p_i}^{a_i})$
也就是 $b^r\equiv q(\mathrm{m}\mathrm{o}\mathrm{d}{p_i}^{a_i-yr})$ …？貌似变成了另一种情况？
（注意这里 $q=0$ 还要分类讨论）
求这条同余方程的解 $b$ 的个数。可以按照下面的来，
然后还原回 $x$ 的方案数的时候要乘上 ${p_i}^{a_i-t}$ （很显然吧？）
（ $q=0$ 的时候就不用了）

最后的情况是 $(a,{p_i}^{a_i})=1$
就是求 $r$ 次剩余的解数，我们把它化成一次同余即可，用指标
得到 $r{I}_g(x)\equiv I_g(a)(\mathrm{m}\mathrm{o}\mathrm{d}\varphi ({p_i}^{a_i}))$
由于指标和原数一一对应那么方案数一样
同时有解条件 $(a,\varphi ({p_i}^{a_i}))|I_g(a)$ 此时恰有 $(a,\varphi ({p_i}^{a_i}))$ 个解。
就是这样啦。

中间有用到了原根跟指标。因为 $P$ 是奇数所以模数 $p_i^{a_i}$ 中的 $p_i$ 会是奇素数，就没有问题啦。
然后我们开始写代码……
……
……
……
然后最后复杂度勉勉强强也许大概能够卡过去
加油（？）

辣辣讲这道题可以不用原根
但是我比较鰯

---

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cctype>
#include<cmath>
#include<ctime>
#include<algorithm>
#include<iostream>
#include<map>
using namespace std;
const long long HM=192617;
int T;
long long A, B, K, P;
bool notPrime[35005];
long long PrimeList[35005];
long long pi[35005];
long long ai[35005];
int Prep[35005];
long long pw[35005];
struct Node
{
	long long Index, Value;
	int Next;
	Node(long long a = -1, long long b = -1, int c = -1)
	{
		Index = a; Value = b; Next = c;
	}
	bool IsEmpty()
	{
		return (Index==-1)&&(Value==-1)&&(Next==-1);
	}
	void Clear()
	{
		Index = -1;
		Value = -1;
		Next = -1;
	}
} MMP[192626];
int GetHashPos(long long OriVal)
{
	return (int) (OriVal%HM);
}
long long qpow(long long bas, long long ind)
{
	long long ret = 1;
	while (ind)
	{
		if (ind & 1) 
		{
			ret *= bas;
		}
		bas *= bas;
		ind >>= 1;
	}
	return ret;
}
long long qpowm(long long bas, long long ind, long long m)
{
	bas %= m;
	long long ret = 1;
	while (ind)
	{
		if (ind & 1) ret *= bas, ret %= m;
		bas *= bas, bas %= m;
		ind >>= 1;
	}
	return ret;
}
#define adj(x) ((x>=HM)?(x-HM):x)
void HashInsert(long long Val, long long Ind)
{
	int HashPos = GetHashPos(Val);
	while (MMP[HashPos].Next != -1) HashPos = MMP[HashPos].Next;
	if (!MMP[HashPos].IsEmpty())
	{
		while (!MMP[MMP[HashPos].Next].IsEmpty()) MMP[HashPos].Next=adj(1+MMP[HashPos].Next);
		HashPos = MMP[HashPos].Next;
	}
	Prep[++Prep[0]] = HashPos; 
	MMP[HashPos] = Node(Ind, Val);
}
long long HashQuery(long long Val)
{
	int HashPos = GetHashPos(Val); 
	while (HashPos != -1 && MMP[HashPos].Value != Val) HashPos = MMP[HashPos].Next;
	if (HashPos==-1 || MMP[HashPos].IsEmpty()) return -1;
	return MMP[HashPos].Index;
}
void PrePrime()
{
	notPrime[0] = notPrime[1] = 1;
	for (register int j, t, i = 2; i <= 34999; ++i)
	{
		if (!notPrime[i]) PrimeList[++PrimeList[0]] = 1ll*i;
		j = 1;
		while (true)
		{
			if ((j >= PrimeList[0]) || ((t = i * PrimeList[++j]) > 34999)) break;
			notPrime[t] = 1;
		}
	}
}
void Split(long long x)
{
	for (register int i = 1; i <= PrimeList[0] && PrimeList[i] * PrimeList[i] <= x; ++i)
	{
		if (x % PrimeList[i] == 0)
		{
			pi[++pi[0]] = PrimeList[i];
			ai[pi[0]] = 0;
			while (x % PrimeList[i] == 0)
			{
				++ai[pi[0]];
				x /= PrimeList[i];
			}
			pw[pi[0]] = qpow(pi[pi[0]],ai[pi[0]]);
		}
	}
	if (x > 1)
	{
		pi[++pi[0]] = x;
		ai[pi[0]] = 1;
		pw[pi[0]] = x;
	}
}
long long Solve0(const int& tmp)
{
	return qpow(1ll*pi[tmp], 1ll*ai[tmp]-1ll*ceil(1.0*ai[tmp]/A));
}
void exgcd(long long a, long long b, long long& d, long long& x, long long& y)
{
	if (!b)
	{
		d = a;
		x = 1;
		y = 0;
		return;
	}
	exgcd(b, a%b, d, y, x);
	y -= x * (a / b);
}
long long gcd(long long a, long long b)
{
	return !b?a:gcd(b,a%b);
}
long long BSGS(long long g, long long a, long long p)
{
	static long long sqr;
	
	while (Prep[0]) MMP[Prep[Prep[0]--]].Clear();
	sqr = 1ll * ceil(sqrt(1.0*p));
	for (register long long tem = a, i = 0; i < sqr; ++i, tem = (tem * g) % p)
	{
		HashInsert(tem, i);
	}
	for (register long long tem, i = 1, t = qpowm(g, sqr, p), s = t; i <= sqr; ++i, s *= t, s %= p)
	{
		tem = HashQuery(s);
		if(tem != -1) return i * sqr - tem;
	}
	return 0;
}
int TemL[35005];
void PMSplit(int x)
{
	TemL[0] = 0;
	for (register int i = 1; i <= PrimeList[0] && PrimeList[i] * PrimeList[i] <= x; ++i)
	{
		if (x % PrimeList[i] == 0)
		{
			TemL[++TemL[0]] = PrimeList[i];
			while (x % PrimeList[i] == 0) x /= PrimeList[i];
		}
	}
	if (x > 1) TemL[++TemL[0]] = x;
}
long long GetPrimRoot(long long p, long long phi)
{
	PMSplit(phi);
	bool fai_led = 0;
	for (register long long x = 2; x < p; ++x)
	{
		fai_led = 0;
		for (register int j = 1; j <= TemL[0]; ++j) 
		{
			if (qpowm(x, phi/TemL[j], p) == 1)
			{
				fai_led = 1;
				break;
			}
		}
		if (!fai_led) return x;
	}
	return 0;
}
long long Solve(const int& tmp)
{
	static long long t, q, g, pp, phi, I, Ret;
	
	pp = qpow(pi[tmp], ai[tmp]);
	q = B;
	if (q%pw[tmp]==0) return Solve0(tmp);
	t = 0;
	while (q%pi[tmp] == 0) 
	{
		q /= pi[tmp], ++t;
		pw[tmp] /= pi[tmp];
	}
	if (t % A) return 0;
	phi = pw[tmp] - pw[tmp]/pi[tmp];
	g = GetPrimRoot(pw[tmp], phi);
	I = BSGS(g, q, pw[tmp]);
	Ret = gcd(A, phi);
	if (I % Ret) return 0;
	return Ret*qpow(pi[tmp],t-t/A);
}
int main()
{
	long long Ans;
	
	PrePrime();
	scanf("%d", &T);
	while (T--)
	{
		scanf("%lld%lld%lld", &A, &B, &K);
		P = 2ll * K + 1;
		pi[0] = 0;
		Split(P);
		Ans = 1;
		for (int i = 1; i <= pi[0]; ++i)
		{
			Ans *= Solve(i);
			if (!Ans) break;
		}
		printf("%lld\n", Ans);
	}	
	return 0;
}

DATA_GENERATOR
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<sstream>
#include<iostream>
#include<cmath>
#include<ctime>
using namespace std;
int seed;
stringstream ss;
int Gen(int X)
{
	return (int)(1.0*rand()*X/RAND_MAX);
}
int main(int argc, char**argv)
{
	if (argc>1)
	{
		ss.clear();
		ss<<argv[1];
		ss>>seed;
	}
	srand(seed);
	int T=1;
	cout<<T<<endl;
	while(T--)
	{
		cout<<Gen(1e9)<<" "<<Gen(1e9)<<" "<<Gen(5e8)<<endl;
	}
	return 0;
}