原題地址這題真是太神犇了……可以讓人完全搞懂?dāng)?shù)論同余部分的全部內(nèi)容……
題解……由于虹貓大神已經(jīng)在空間里寫得很詳細(xì)了，所以就不腫么寫了囧……
主要說一下一些難想的和容易搞疵的地方：
（1）中國剩余定理的那個(gè)推論（多個(gè)同余方程的模數(shù)互質(zhì)，則整個(gè)方程組在小于所有模數(shù)之積的范圍內(nèi)的解數(shù)等于各個(gè)方程解數(shù)之積）其實(shí)是很強(qiáng)大的，不光對線性同余方程有用，對這種非線性的同余方程也有用，只需要所有方程都滿足：若模數(shù)為MOD，則a是解當(dāng)且僅當(dāng)(a+MOD)是解……本題顯然滿足，因此，只要在按質(zhì)因數(shù)拆分后求出各個(gè)方程的解數(shù)，再相乘即可（本沙茶就是這里木有想起來，結(jié)果看了虹貓的題解囧……）；
（2）對于余數(shù)不為0且和模數(shù)不互質(zhì)的情況要特別注意（這個(gè)好像很多標(biāo)程都疵了，比如虹貓給的標(biāo)程，不過數(shù)據(jù)弱，讓它們過了囧），首先必須是余數(shù)含p（p為該方程模數(shù)的質(zhì)因數(shù)）因子的個(gè)數(shù)j是a的倍數(shù)（也就是余數(shù)是p^a的倍數(shù)）才能有解，然后，當(dāng)有解時(shí)，轉(zhuǎn)化為解必須是p^(j/a)的倍數(shù)以及x/(p^(j/a))滿足一個(gè)模數(shù)指數(shù)為原來指數(shù)減j的方程，這里需要注意，這個(gè)新方程的解數(shù)乘以p^(j-j/a)才是原來方程的解數(shù)！！道理很簡單，因?yàn)槟?shù)除以了p^j，而x只除以了p^(j/a)……可以用一組數(shù)據(jù)檢驗(yàn)：3 330750 6643012，結(jié)果是135而不是15；
（3）原根只能暴力求（不過最小原根都很小，1000以內(nèi)的所有質(zhì)數(shù)最小原根最大只有19……），但在求的時(shí)候有一個(gè)小的優(yōu)化：首先p的原根也是p的任意整數(shù)次方的原根，然后求p的原根時(shí)，將(p-1)的非自身因數(shù)（預(yù)先求出）遞減排序，這樣可以比較快地排除不合法解；
（4）求逆元時(shí)一定要注意，如果得到的逆元是負(fù)數(shù)，要轉(zhuǎn)化為正數(shù)，另外要取模；
（5）BSGS的時(shí)候一定要注意去重，在保留重復(fù)元素的情況下即使使用另一種二分查找也會(huì)疵的；
（6）數(shù)組不要開小了；
代碼：
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <algorithm>
using namespace std;
#define re(i, n) for (int i=0; i<n; i++)
#define re1(i, n) for (int i=1; i<=n; i++)
#define re2(i, l, r) for (int i=l; i<r; i++)
#define re3(i, l, r) for (int i=l; i<=r; i++)
#define rre(i, n) for (int i=n-1; i>=0; i--)
#define rre1(i, n) for (int i=n; i>0; i--)
#define rre2(i, r, l) for (int i=r-1; i>=l; i--)
#define rre3(i, r, l) for (int i=r; i>=l; i--)
#define ll long long
const int MAXN = 110, MAXP = 50010, INF = ~0U >> 2;
int P_LEN, _P[MAXP + 1], P[MAXP + 1];
int A, B, M, n, DS[MAXN], DK[MAXN], R[MAXN], KR[MAXP], res;
struct sss {
int v, No;
bool operator< (sss s0) const {return v < s0.v || v == s0.v && No < s0.No;}
} Z[MAXP];
void prepare0()
{
P_LEN = 0; int v0;
for (int i=2; i<=MAXP; i++) {
if (!_P[i]) P[P_LEN++] = _P[i] = i; v0 = _P[i] <= MAXP / i ? _P[i] : MAXP / i;
for (int j=0; j<P_LEN && P[j]<=v0; j++) _P[i * P[j]] = P[j];
}
}
void prepare()
{
n = 0; int M0 = M;
re(i, P_LEN) if (!(M0 % P[i])) {
DS[n] = P[i]; DK[n] = 1; M0 /= P[i]; while (!(M0 % P[i])) {DK[n]++; M0 /= P[i];} n++;
if (M0 == 1) break;
}
if (M0 > 1) {DS[n] = M0; DK[n++] = 1;}
int x;
re(i, n) {
x = 1; re(j, DK[i]) x *= DS[i];
R[i] = B % x;
}
}
ll pow0(ll a, int b, ll MOD)
{
if (b) {ll _ = pow0(a, b >> 1, MOD); _ = _ * _ % MOD; if (b & 1) _ = _ * a % MOD; return _;} else return 1;
}
void exgcd(int a, int b, int &x, int &y)
{
if (b) {
int _x, _y; exgcd(b, a % b, _x, _y);
x = _y; y = _x - a / b * _y;
} else {x = 1; y = 0;}
}
int gcd(int a, int b)
{
int r = 0; while (b) {r = a % b; a = b; b = r;} return a;
}
void solve()
{
int x, y; res = 1;
re(i, n) if (!R[i]) {
if (DK[i] < A) x = 1; else x = (DK[i] - 1) / A + 1;
re2(j, x, DK[i]) res *= DS[i];
} else if (!(R[i] % DS[i])) {
x = 0; while (!(R[i] % DS[i])) {R[i] /= DS[i]; x++;}
if (x % A) {res = 0; return;} else {
DK[i] -= x; y = x / A;
re2(j, y, x) res *= DS[i];
}
}
int phi, m0, m1, KR_len, _r, v0, _left, _right, _mid, T; bool FF;
re(i, n) if (R[i]) {
x = DS[i] - 1; KR_len = 0;
for (int j=2; j*j<=x; j++) if (!(x % j)) {
KR[KR_len++] = j;
if (j * j < x) KR[KR_len++] = x / j;
}
KR[KR_len++] = 1;
re2(j, 2, DS[i]) {
FF = 1;
rre(k, KR_len) {
_r = (int) pow0(j, KR[k], DS[i]);
if (_r == 1) {FF = 0; break;}
}
if (FF) {x = j; break;}
}
phi = DS[i] - 1; re2(j, 1, DK[i]) phi *= DS[i]; v0 = phi / (DS[i] - 1) * DS[i];
m0 = (int) ceil(sqrt(phi) - 1e-10);
Z[0].v = 1; Z[0].No = 0; re2(j, 1, m0) {Z[j].v = (ll) Z[j - 1].v * x % v0; Z[j].No = j;}
_r = (ll) Z[m0 - 1].v * x % v0; sort(Z, Z + m0);
m1 = 1; re2(j, 1, m0) if (Z[j].v > Z[j - 1].v) Z[m1++] = Z[j];
exgcd(_r, v0, x, y); if (x < 0) x += v0; y = R[i];
re(j, m0) {
_left = 0; _right = m1 - 1; T = -1;
while (_left <= _right) {
_mid = _left + _right >> 1;
if (y == Z[_mid].v) {T = j * m0 + Z[_mid].No; break;}
else if (y < Z[_mid].v) _right = _mid - 1; else _left = _mid + 1;
}
if (T >= 0) break; else y = (ll) y * x % v0;
}
x = gcd(A, phi); if (T % x) {res = 0; break;} else res *= x;
}
}
int main()
{
int tests;
scanf("%d", &tests);
prepare0();
re(testno, tests) {
scanf("%d%d%d", &A, &B, &M); M += M + 1; B %= M;
if (!A) {
if (B == 1) res = M; else res = 0;
} else {
prepare();
solve();
}
printf("%d\n", res);
}
return 0;
}