题目地址：洛谷：【P2480】[SDOI2010]古代猪文 – 洛谷、BZOJ：Problem 1951. — [Sdoi2010]古代猪文

题目描述

猪王国的文明源远流长，博大精深。
iPig在大肥猪学校图书馆中查阅资料，得知远古时期猪文文字总个数为N。当然，一种语言如果字数很多，字典也相应会很大。当时的猪王国国王考虑到如果修一本字典，规模有可能远远超过康熙字典，花费的猪力、物力将难以估量。故考虑再三没有进行这一项劳猪伤财之举。当然，猪王国的文字后来随着历史变迁逐渐进行了简化，去掉了一些不常用的字。
iPig打算研究古时某个朝代的猪文文字。根据相关文献记载，那个朝代流传的猪文文字恰好为远古时期的k分之一，其中k是N的一个正约数（可以是1和N）。不过具体是哪k分之一，以及k是多少，由于历史过于久远，已经无从考证了。
iPig觉得只要符合文献，每一种能整除N的k都是有可能的。他打算考虑到所有可能的k。显然当k等于某个定值时，该朝的猪文文字个数为N / k。然而从N个文字中保留下N / k个的情况也是相当多的。iPig预计，如果所有可能的k的所有情况数加起来为P的话，那么他研究古代文字的代价将会是G的P次方。
现在他想知道猪王国研究古代文字的代价是多少。由于iPig觉得这个数字可能是天文数字，所以你只需要告诉他答案除以999911659的余数就可以了。

题意简述

计算$G^{\sum_{d|N} \mathrm{C}^d_N} \bmod 999911659$。

输入输出格式

输入格式：
输入文件ancient.in有且仅有一行：两个数N、G，用一个空格分开。

输出格式：
输出文件ancient.out有且仅有一行：一个数，表示答案除以999911659的余数。

输入输出样例

输入样例#1：

4 2

输出样例#1：

2048

说明

数据规模
10%的数据中，1 <= N <= 50；
20%的数据中，1 <= N <= 1000；
40%的数据中，1 <= N <= 100000；
100%的数据中，1 <= G <= 1000000000，1 <= N <= 1000000000。

题解

首先，组合数一堆加起来幂次可能很大，这时我们可以应用欧拉降幂公式降低幂次。我们把答案变成了这样子
$$ G^{\sum_{d|N} \mathrm{C}^d_N \bmod \varphi(999911659) + \varphi(999911659)} $$
由于$999911659$是个质数，我们知道$\varphi(999911659)=999911658$，而这个数并不是质数，并不好办。我们考虑使用扩展Lucas定理，将$999911658$分解质因数为$2 \times 3 \times 4679 \times 35617$，分别对每个质因数求一次组合数，然后应用中国剩余定理解这个同余方程组，得到最终结果。
复杂度$O(\sqrt{n} \log n)$。

代码

// Code by KSkun, 2018/7
#include <cstdio>
#include <cctype>
#include <cstring>

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2)
        ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1; register char c = fgc();
    for(; !isdigit(c); c = fgc()) if(c == '-') neg = -1;
    for(; isdigit(c); c = fgc()) res = (res << 1) + (res << 3) + c - '0';
    return res * neg;
}

inline LL fpow(LL n, LL k, LL p) {
    LL res = 1;
    for(; k; k >>= 1) {
        if(k & 1) res = res * n % p;
        n = n * n % p;
    }
    return res;
}

const int MAXN = 100005, MO[] = {0, 2, 3, 4679, 35617, 999911659};

LL N, G, a[5], mul[MAXN], inv[MAXN];

inline void init(LL p) {
    mul[0] = mul[1] = inv[0] = inv[1] = 1;
    for(int i = 2; i < p; i++) {
        mul[i] = mul[i - 1] * i % p;
    }
    for(int i = 2; i < p; i++) {
        inv[i] = (-(p / i) * inv[p % i] % p + p) % p;
    }
    for(int i = 2; i < p; i++) {
        inv[i] = inv[i - 1] * inv[i] % p;
    }
}

LL C(LL n, LL m, LL p) {
    if(m > n) return 0;
    if(m == 0) return 1;
    if(n < p && m < p) return mul[n] * inv[m] * inv[n - m] % p;
    return C(n % p, m % p, p) * C(n / p, m / p, p) % p;
}

inline LL calSum(LL p) {
    init(p);
    LL res = 0;
    for(int i = 1; i * i <= N; i++) {
        if(N % i) continue;
        res = (res + C(N, i, p)) % p;
        if(i * i == N) continue;
        res = (res + C(N, N / i, p)) % p;
    }
    return res;
}

inline LL crt() {
    LL p = 1, res = 0;
    for(int i = 1; i <= 4; i++) {
        p *= MO[i];
    }
    for(int i = 1; i <= 4; i++) {
        LL m = p / MO[i];
        res = (res + a[i] * m * fpow(m, MO[i] - 2, MO[i])) % p;
    }
    return res;
}

int main() {
    N = readint(); G = readint();
    if(G % MO[5] == 0) {
        puts("0"); return 0;
    }
    for(int i = 1; i <= 4; i++) {
        a[i] = calSum(MO[i]);
    }
    printf("%lld", fpow(G, crt() + MO[5] - 1, MO[5]));
    return 0;
}