hdu4335(What is N?)(数论)

http://acm.hdu.edu.cn/showproblem.php?pid=4335

这是道数论的题目，会的人可以瞬秒。

（1）先要熟悉欧拉定理：gcd(n, p) = 1时，  n^(phi(p)) ≡ 1 (mod p)；

（2）各种条件下的n:

       1)  n! < phi(p) ；

       2）phi(p) <= n! < s!, s=min{x|x! mod phi(p)=0}；

       3）n! % phi(p) = 0；

       4)其他；

（3）其实没多少看得懂。

他人具体代码：

#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define MAXN 100000
using namespace std;

typedef unsigned long long int Int64;

int p[MAXN+5], pm[10000], phi[MAXN+5], idx = -1, MOD, B;

int Fac[MAXN+5];

Int64 M, ans;

void GetPrime()
{
    for (int i = 2; i <= MAXN; ++i) {
        if (!p[i]) {
            pm[++idx] = i;                
        }
        for (int j = 0; pm[j]*i <= MAXN; ++j) {
            p[pm[j]*i] = 1;
            if (i % pm[j] == 0) {
                break;    
            }
        }
    }
/*    printf("%d\n", idx);
    for (int i = 0; i <= 100; ++i) {
        printf("%d ", pm[i]);
    }*/
}

void Eular()
{
    phi[1] = 1;
    for (int i = 2; i <= MAXN; ++i) {
    //    printf("phi[%d] = %d\n", i-1, phi[i-1]);
    //    getchar();
        if (!p[i]) {
            phi[i] = i - 1;    
            continue;
        }
        for (int j = 0; pm[j]*pm[j] <= i; ++j) {
            if (i % pm[j] == 0) {
                if (i / pm[j] % pm[j] == 0) {
                    phi[i] = pm[j] * phi[i/pm[j]];
                }
                else {
                    phi[i] = phi[pm[j]] * phi[i/pm[j]];
                }
                break;
            }
        }
        
    }
}

bool _pow(Int64 a, Int64 b)
{
    Int64 ret = 1;
    while (b) {
        if (b & 1) {
            ret *= a;
            ret %= MOD;
        }
        a *= a;
        a %= MOD;
        b >>= 1;
    }
    return ret == B;
}

void deal()
{
    int LIM = -1, ptr;
    for (ptr = 0; ptr <= M; ++ptr) {
        // 第一段，小于phi[MOD]的一部分 
        if (!ptr) { Fac[ptr] = 1; }
        else { Fac[ptr] = Fac[ptr-1] * ptr; }
        if (Fac[ptr] >= phi[MOD]) {
            break;
        }
        if (_pow(ptr, Fac[ptr])) {
            ++ans;
        }
    }
    if (ptr > M) { return; } // 如果已经超过了M的话 
    for (int i = ptr; i <= M; ++i) {
        if (!i) { Fac[i] = 1 % phi[MOD]; }
        else { Fac[i] = (Fac[i-1] * i) % phi[MOD]; }
        if (Fac[i] == 0) {  // 如果找到了这个点 
            LIM = i;
            break;
        }
        if (_pow(i, Fac[i]+phi[MOD])) {
            ++ans;
        }
    }
    if (LIM == -1) { return; } // 如果没有找到一个是phi[MOD]倍数的点 
    int t = 0;
    if ((M - LIM+1) >= MOD) { // LIM 这个点是没有计算的，所以长度为M-LIM+1，这里判定剩余是否大于MOD 
        for (int i = LIM; i < LIM+MOD; ++i) {
            if (_pow(i%MOD, phi[MOD])) { // 由于此时已经是phi[MOD]的倍数，所以直接把phi[MOD]这个指数代入
                ++t;
            }
        }
        ans += (M - LIM+1) / (1LLU*MOD) * t;
    }
    int left = (M-LIM+1) % MOD; // 残余的余数 
    for (int i = LIM; i < LIM+left; ++i) {
        if (_pow(i%MOD, phi[MOD])) {
            ++ans;
        }
    }
}

int main()
{
    GetPrime();
    Eular();
    int T, ca = 0;
    scanf("%d", &T);
    while (T--) {
        ans = 0;
        scanf("%d %d %I64u", &B, &MOD, &M);
        if(M == 18446744073709551615LLU && MOD == 1 && B == 0) {
            // 由于b=0,p=1,所以可以直接得到M+1，由于溢出原因所以需要特判 
            printf("Case #%d: 18446744073709551616\n", ++ca);
        }
        else 
        {    
            deal(); 
            printf("Case #%d: %I64u\n", ++ca, ans);
        }
    }
    return 0;
}