[Sdoi2010]古代猪文 (卢卡斯定理,欧拉函数)
哇,这道题真的好好,让我这个菜鸡充分体会到卢卡斯和欧拉函数的强大!
先把题意抽象出来!就是计算这个东西。
p=999911659是素数,p-1=2*3*4679*35617
所以:这样只要求出
然后再快速乘法就行了。
那好,怎么做呢?
有模运算的性质得到 然后就是卢卡斯原理。
先把卢卡斯原理放这里:
void init(int mod){ //对mod取余后,一定小于mod,因此把mod的阶乘存起来就够用 f[0] = 1; for (int i = 1; i <= mod; i++){ f[i] = f[i - 1] * i % mod; } } void ex_gcd(LL a, LL b, LL& d, LL& x, LL& y) { if (!b) { d = a; x = 1; y = 0; } 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 == 1 ? (x + m) % m : -1; } LL Lucas(LL m, LL n, LL p){ LL res = 1; 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 = 35617; int N, G, fact[maxn + 5], mod = 999911658; int prime[5] = { 0, 2, 3, 4679, 35617 }, Mod[5]; void get_fact() { fact[0] = 1; for (int i = 1; i <= maxn; i++) fact[i] = (ll)fact[i - 1] * i%mod; } int ex_t; void exgcd(int a, int b, int &x, int &y) { if (!b) { x = 1; y = 0; 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 = 1, 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 ? 0 : fact[n - m], p) % p; } return ret; } ll pow(int x, int n) { int ans = 1; for (; n;n>>=1, x=(ll)x*x%mod) if (n & 1)ans = (ll)ans*x%mod; return ans; } int main() { scanf("%d%d", &N, &G); if (G % (mod + 1) == 0){ printf("0"); return 0; } get_fact(); //得到阶乘 for (int i = 1; i*i <= N; ++i) //枚举因子 { if (N%i == 0) { for (int j = 1; j <= 4; ++j)Mod[j] = (Mod[j] + calc(i, prime[j])) % prime[j]; if (i*i!=N) for (int j = 1; j <= 4; ++j)Mod[j] = (Mod[j] + calc(N / i, prime[j])) % prime[j]; } } int x, y, b = 0; for (int i = 1; i <= 4; 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 += 1; //printf("%d\n", b); printf("%lld", pow(G, b)); }