[sdoi2010]古代猪文

写这道题首先srO qw。。。。题目就不贴了。。。大意：求 (G^(Σd|n C(n, d))) mod M。

首先发现模是一个质数，所以根据费马小定理可知指数只需模模减一即可。

然后又发现模减一太大，而且不是质数。。。然后发现质因子分解后都是比较小的数字。。果断中国剩余定理合并啊。。。

然后记得lyp讲过什么n的正约数个数均摊是O(logn)级的。。。果断暴枚约数啊。。。

然后是组合数取模。。。n太大，实在想不出什么好办法直接处理阶乘，好像可以用勒让德定理，但是觉得太麻烦了。。。于是用zfone讲的lucas定理做吧。。。

然后组合数就线性预处理阶乘的模和逆元，每次直接用好了。。。

然后时间就被虐爆了。。。不知道好多版去了。。。

代码写了一百多行。。。各种繁琐。。。脑残错误有一开始搞错了模。。。

#include <cstdio>#include <cmath>#include <cstdlib>#include <cstring>#include <ctime>#include <cctype>#include <string>#include <map>#include <set>#include <queue>#include <algorithm>#include <fstream>#include <iostream>using namespace std;#ifdef WIN32#define fmt64 "%I64d"#else#define fmt64 "%lld"#endif#define PI M_PI#define oo 0x13131313#define iter iterator#define PB push_back#define PO pop_back#define MP make_pair#define fst first#define snd second#define cstr(a) (a).c_str()#define FOR(i, j, k) for(i = (j); i <= (k); ++i)#define ROF(i, j, k) for(i = (j); i >= (k); --i)#define FER(e, d, u) for(e = d[u]; e; e = e->n)#define FRE(i, a) for(i = (a).begin(); i != (a).end(); ++i)typedef unsigned int uint;typedef long long int64;typedef unsigned long long uint64;typedef long double real;template<class T> inline bool minim(T &a, const T &b) {return b < a ? a = b, 1 : 0;}template<class T> inline bool maxim(T &a, const T &b) {return b > a ? a = b, 1 : 0;}template<class T> inline T sqr(const T &a) {return a * a;}#define maxn 100005#define maxm 35700int pt, pm[maxn / 5];bool np[maxn + 1];int ft, fac[20], cnt[20];int dt, d[maxn];void get_prime(){int i, j, k;FOR(i, 2, maxn) {if (!np[i]) pm[++pt] = i;FOR(j, 1, pt) {if ((k = i * pm[j]) > maxn) break;np[k] = 1;if (!(i % pm[j])) break;}}}void divide(int k){int i, j; j = ft = 0;FOR(i, 1, pt) {if (sqr(pm[i]) > k) break;if (!(k % pm[i])) {fac[++ft] = pm[i], cnt[ft] = 1;for (k /= pm[i]; !(k % pm[i]); k /= pm[i]) ++cnt[ft];}}if (k != 1) fac[++ft] = k, cnt[ft] = 1;}void dfs(int t, int rest, int now){if (!rest) {d[++dt] = now; return;}if (t > ft) return;dfs(t + 1, rest, now);for (int i = 1; i <= cnt[t]; ++i) {if (i > rest) break;dfs(t + 1, rest - i, now *= fac[t]);}}const int mm = 2 * 3 * 4679 * 35617;const int m[4] = {2, 3, 4679, 35617};const int M[4] = {499955829, 333303886, 213702, 28074};int inv[4];int f[4][maxm + 1], v[4][maxm + 1];int ex_gcd(int a, int b, int &x, int &y){if (!b) return x = 1, y = 0, a;int ret = ex_gcd(b, a % b, x, y), t;return t = x, x = y, y = t - (a / b) * y, ret;}int inverse(int a, int b){int x, y;ex_gcd(a, b, x, y);if ((x %= b) < 0) x += b;return x;}void prepare(){int i, j, k, l;FOR(l, 0, 3) {inv[l] = inverse(M[l], m[l]);f[l][0] = f[l][1] = v[l][0] = v[l][1] = 1;}FOR(l, 0, 3) FOR(i, 2, m[l] - 1) {f[l][i] = (f[l][i - 1] * i) % m[l];if (!np[i]) v[l][i] = inverse(i, m[l]);FOR(j, 1, pt) {if ((k = pm[j] * i) > m[l]) break;(v[l][k] = v[l][i] * v[l][pm[j]]) %= m[l];if (!(i % pm[j])) break;}}FOR(l, 0, 3) FOR(i, 2, m[l] - 1)(v[l][i] *= v[l][i - 1]) %= m[l];}int C(int a, int b, int l){return a < b ? 0 : (int64(f[l][a] * v[l][b]) * v[l][a - b]) % m[l];}int lucas(int a, int b, int l){int64 ret = 1;for (; b; a /= m[l], b /= m[l])(ret *= C(a % m[l], b % m[l], l)) %= m[l];return ret;}int main(){freopen("ancient.in", "r", stdin);freopen("ancient.out", "w", stdout);int i, j; int64 k; j = k = 0;int n; int64 G; cin >> n >> G;if (G == mm + 1) puts("0"), exit(0);get_prime(), divide(n);FOR(i, 1, ft) j += cnt[i];FOR(i, 0, j) dfs(1, i, 1);prepare();FOR(i, 1, dt) {int64 sum = 0;FOR(j, 0, 3) {int64 t = lucas(n, d[i], j);(sum += inv[j] * M[j] * t) %= mm;}(k += sum) %= mm;}int64 ans = 1;for (i = k; i; (G *= G) %= mm + 1, i >>= 1)if (i & 1) (ans *= G) %= mm + 1;cout << ans << endl;return 0;}