HDU 4335 What is N?

#include <cstdio>
#include <cstring>
#include <vector>
#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <bitset>
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
typedef vector<int> VI;
typedef vector<VI> VVI;
const int maxn = 100000;
const int N = 11111;
const double eps = 1e-8;
bitset<maxn + 100> vis;
int pri[maxn], size;
void init() {
    int i, j;
    for (i = 2; i * i <= maxn; ++i)
        if (!vis[i])
            for (j = i * i; j <= maxn; vis[j] = 1, j += i)
                ;

    pri[size++] = 2;
    for (i = 3; i <= maxn; i += 2)
        if (!vis[i])
            pri[size++] = i;
}
ULL POW(ULL a, ULL b, ULL mod) {
    ULL res = 1;
    while (b) {
        if (b & 1)
            res = res * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return res;
}
int euler(int n) {
    int i, res = n;
    for (i = 0; i < size && pri[i] * pri[i] <= n; ++i) {
        if (n % pri[i] == 0) {
            res = res - res / pri[i];
            while (n % pri[i] == 0)
                n /= pri[i];
        }
    }
    if (n > 1)
        res = res - res / n;
    return res;
}
ULL remain[maxn];
ULL b, m;
int p;
int main() {
    ios::sync_with_stdio(false);
    init();
    int T, cas;
    ULL i;
    scanf("%d", &T);
    for (cas = 1; cas <= T; ++cas) {
        scanf("%I64u%d%I64u", &b, &p, &m);
        printf("Case #%d: ",cas);
        if (p == 1) {
            if (m == 18446744073709551615ULL)
                printf("18446744073709551616\n");
            else
                printf("%I64u\n", m + 1);
            continue;
        }
        ULL fac = 1;
        int sign = 0;
        ULL ans = 0, phi = euler(p);
        for (i = 0; i <= m && fac; ++i) {
            ans += POW(i, fac + sign * phi, p) == b;
            fac *= i + 1;
            sign = sign | (fac >= phi);
            fac %= phi;
        }
        int j;
        ULL cnt = 0;
        if (i <= m) {
            for (j = 0; j < p; ++j) {
                remain[j] = POW(i + j, phi, p);
                cnt += remain[j] == b;
            }
            ans += (m - i + 1) / p * cnt;
            int r = (m - i + 1) % p;
            for (j = 0; j < r; ++j) {
                ans += remain[j] == b;
            }
        }
        printf("%I64u\n",ans);
    }
    return 0;
}