[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 35700

int 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;
}