哇，这道题真的好好，让我这个菜鸡充分体会到卢卡斯和欧拉函数的强大！
先把题意抽象出来！就是计算这个东西。
p=999911659是素数，p-1=2*3*4679*35617
所以：这样只要求出
然后再快速乘法就行了。
那好，怎么做呢？
有模运算的性质得到 然后就是卢卡斯原理。
先把卢卡斯原理放这里：
void init(int mod){ //对mod取余后，一定小于mod，因此把mod的阶乘存起来就够用
f[] = ;
for (int i = ; i <= mod; i++){
f[i] = f[i - ] * i % mod;
}
}
void ex_gcd(LL a, LL b, LL& d, LL& x, LL& y)
{
if (!b) { d = a; x = ; y = ; }
else{ ex_gcd(b, a%b, d, y, x); y -= x*(a / b); }
}
LL inv(LL a, LL m)
{
LL d, x, y;
ex_gcd(a, m, d, x, y);
return d == ? (x + m) % m : -;
}
LL Lucas(LL m, LL n, LL p){
LL res = ;
while (n && m){
LL n1 = n % p;
LL m1 = m % p;
res = res * f[n1] * inv(f[n1 - m1], p) * inv(f[m1], p) % p;
n /= p;
m /= p;
}
return (res % p + p) % p;
}
则：那么我们其实把它每个存起来Mod[1-4]
然后，就是要找一个值来代替Mod[1-4]。利用中国剩余定理！（哇，太难打了公式了）
什么这样做？因为能同时被2，3，4679，35617那么一定会被99991165同余，那么这个数就是
注意：坑！快速幂一定要加long long，找了3小时的bug
#include<cstdio>
using namespace std;
#define ll long long
const int maxn = ;
int N, G, fact[maxn + ], mod = ;
int prime[] = { , , , , }, Mod[]; void get_fact() { fact[] = ; for (int i = ; i <= maxn; i++) fact[i] = (ll)fact[i - ] * i%mod; } int ex_t; void exgcd(int a, int b, int &x, int &y) { if (!b) { x = ; y = ; return; } exgcd(b, a%b, x, y); ex_t = x; x = y; y = ex_t - (a / b)*y; } int inv(int a, int p) { int x, y; exgcd(a, p, x, y); return (x%p + p) % p; } int calc(int i, int p) { int ret = , x, y, n, m; for (x = N, y = i; y; x /= p, y /= p) { n = x%p; m = y%p; //卢卡斯定理+预处理阶乘+逆元 ret = (ll)ret*fact[n] % p*inv(fact[m], p) % p*inv(n<m ? : fact[n - m], p) % p; } return ret; } ll pow(int x, int n)
{
int ans = ;
for (; n;n>>=, x=(ll)x*x%mod)
if (n & )ans = (ll)ans*x%mod;
return ans;
} int main()
{ scanf("%d%d", &N, &G);
if (G % (mod + ) == ){ printf(""); return ; }
get_fact(); //得到阶乘
for (int i = ; i*i <= N; ++i) //枚举因子
{
if (N%i == )
{
for (int j = ; j <= ; ++j)Mod[j] = (Mod[j] + calc(i, prime[j])) % prime[j];
if (i*i!=N)
for (int j = ; j <= ; ++j)Mod[j] = (Mod[j] + calc(N / i, prime[j])) % prime[j];
}
}
int x, y, b = ;
for (int i = ; i <= ; i++) //中国剩余定理 { exgcd(mod / prime[i], prime[i], x, y); b = (b + (ll)Mod[i] % mod*(mod / prime[i]) % mod*x%mod) % mod; }
b = (b%mod + mod) % mod; mod += ;
//printf("%d\n", b);
printf("%lld", pow(G, b));
}