数论之神 god of number theory

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <ctime>
#include <algorithm>
#include <map>

#define uns unsigned
#define int64 long long
#ifdef WIN32
#define fmt64 "%I64d"
#else
#define fmt64 "%lld"
#endif
#define oo 0x13131313
#define REP(i, n) for (i = 0; i < (n); ++i)
#define maxn 40000
#define MAX (1 << 30)

using namespace std;

int pm[maxn / 5], ptot; bool np[maxn];
int fac[20], ind[20], tot;

void prime()
{
    int i, j;
    for (i = 2; i < maxn; ++i) {
        if (!np[i]) pm[ptot++] = i;
        for (j = 0; j < ptot && i * pm[j] < maxn; ++j) {
            np[i * pm[j]] = 1;
            if (!(i % pm[j])) break;
        }
    }
}

int divide(int n, int fact[], int index[])
{
    int i, tot = 0;
    REP(i, ptot) {
        if (pm[i] * pm[i] > n) break;
        if (!(n % pm[i])) {
            fact[tot] = pm[i], index[tot++] = 1;
            for (n /= pm[i]; !(n % pm[i]); n /= pm[i])
                ++index[tot - 1];
        }
    }
    if (n > 1) fact[tot] = n, index[tot++] = 1;
    return tot;
}

int fpm(int64 a, int b, int c = MAX)
{
    int res = 1;
    for (; b; b >>= 1, a = a * a % c)
        if (b & 1) res = res * a % c;
    return res;
}

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

int gcd(int a, int b) { for (int c; b; c = a, a = b, b = c % b); return a; }

int work(int A, int B, int p, int k)
{
    int i, P = fpm(p, k), res = 0;
    if (B %= P, !B) {
        ++res;
        for (i = (int)ceil((double)k / A); i < k; ++i) {
            int t = fpm(p, k - i) - 1;
            res += t - t / p;
        }
    } else {
        int d = 0;
        for (; !(B % p); B /= p) ++d, --k;
        if (d % A) return 0;

        int phi = P / p * (p - 1);
        int fact[20] = {0}, ind[20] = {0}, tot = divide(phi, fact, ind);
        REP(i, tot) fact[i] = phi / fact[i];
        int g;
        for (g = 2; ; ++g) {
            REP(i, tot)
                if (fpm(g, fact[i], P) == 1) break;
            if (i == tot) break;
        }
        res = gcd(phi, A);
        map<int, int> hash; int SZ = (int)sqrt(phi);
        int64 s, t = 1;
        REP(i, SZ) hash[t] = i, t = t * g % P;
        for (i = 0, s = 1; ; i += SZ, s = s * t % P) {
            int64 x, y;
            euclid(s, P, x, y);
            if (x < 0) x = (-(-x % P) + P);
            x %= P, x = x * B % P;
            if (hash.count(x)) {
                if ((i + hash[x]) % res) res = 0;
                break;
            }
        }
    }
    return res;
}

int main()
{
    freopen("3731.in", "r", stdin);
    freopen("3731.out", "w", stdout);

    int T, A, B, K, i;
    prime();
    for (scanf("%d", &T); T--; ) {
        scanf("%d%d%d", &A, &B, &K);
        tot = divide(K << 1 | 1, fac, ind);
        int ans = 1;
        REP(i, tot) {
            ans *= work(A, B, fac[i], ind[i]);
            if (!ans) break;
        }
        printf("%d\n", ans);
    }
}