人生又能几多愁,恰是一江东水向东流
暴力求解
- 一次枚举n的每个约数d计算C(n, d) % mod, mod = 9999116598, 不为质数因此就不能用lucas和逆元了,只能用阶乘分解求出所有质因子,最后用快速幂求解
时间复杂度1600 * nlogn 肯定超时
这个东西真JB烦
题意:给你一个a, b, p, 求C(a, b) % p 的值 1e9 >= a, b, p 其中p不为质数
exlucas
- 一般这种题肯定能把p分解质因数 假定为 p = a1 ^b1 * a2^b2 * … * an^bn ,an <= 1e6
- 第一种所有指数全部为1 那么我们一次算出C(a, b) % ai 的值,最后用一遍中国剩余定理即可,为什么是对的?设x = C(a, b), 那么 x = k*p + r, 答案就是r 我们用 x依次mod上p的每个质因子等价于rmod每个质因子,因为ai两两互质,并且周期为p, 因此用中国定理求出的答案一定是对的。
- 本题中由于要求出很多个C(n, d) % ai再求和因此可以把每次的值加起来最后只计算一遍中国剩余定理即可。详情看代码。
- 时间复杂度 1600 * logplogn
- 对于处理数值比较大的情况一般建议define int long long main 函数前面写上signed
#include <iostream>
#include <cstring>
#include <vector>
#define int long long
using namespace std;
typedef pair<int, int> PII;
typedef long long LL;
const int mod = 999911659, N = 50010;
int primes[N], cnt;
bool st[N];
vector<PII> nums, seq;
int a[5], f[N];
void init(int n)
{
for(int i = 2; i <= n; i ++)
{
if(!st[i])primes[cnt ++] = i;
for(int j = 0; primes[j] <= n / i; j ++)
{
st[primes[j] * i] = true;
if(i % primes[j] == 0)break;
}
}
}
int qmi(int a, int b, int p)
{
int res = 1;
while(b)
{
if(b & 1)res = (LL)res * a % p;
a = (LL)a * a % p;
b >>= 1;
}
return res;
}
int Eulor(int n)
{
int res = n;
for(int i = 0; primes[i] <= n / primes[i]; i ++)
if(n % primes[i] == 0)
{
while(n % primes[i] == 0)
n /= primes[i];
res = res / primes[i] * (primes[i] - 1);
}
if(n > 1)res = res / n * (n - 1);
return res;
}
void work(int n, vector<PII> &nums)
{
for(int i = 0; primes[i] <= n / primes[i]; i ++)
{
int p = primes[i];
if(n % p == 0)
{
int s = 0;
while(n % p == 0)
{
n /= p;
s ++;
}
nums.push_back({primes[i], s});
}
}
if(n > 1)nums.push_back({n, 1});
}
void get(int p)
{
f[0] = 1;
for(int i = 1; i < p; i ++)f[i] = (LL)f[i - 1] * i % p;
}
int C(int a, int b, int p)
{
if(a < b)return 0;
return (LL)f[a] * qmi(f[b], p - 2, p) % p * qmi(f[a - b], p - 2, p) % p;
}
int lucas(int a, int b, int p)
{
if(a < p && b < p)return C(a, b, p);
return (LL)C(a % p, b % p, p) * lucas(a / p, b / p, p) % p;
}
vector<int> dir;
void dfs(int u, int sum)
{
if(u == seq.size())
{
dir.push_back(sum);
return ;
}
int p = 1;
for(int i = 0; i <= seq[u].second; i ++)
{
dfs(u + 1, sum * p);
p *= seq[u].first;
}
}
void exgcd(int a, int b, int &x, int &y)
{
if(!b)
{
x = 1;
y = 0;
return ;
}
exgcd(b, a % b, x, y);
int temp = x;
x = y;
y = temp - a / b * y;
}
int crt(int M)
{
int res = 0;
for(int i = 0; i < 4; i ++)
{
int mi = M / nums[i].first;
int x, y;
exgcd(mi, nums[i].first, x, y);
res = (res + mi * x * a[i]) % M;
}
return (res + M) % M;
}
signed main()
{
int n, q;
cin >> n >> q;
init(N - 1);
int fmod = Eulor(mod);
work(fmod, nums);
work(n, seq);
dfs(0, 1);
if(q%(mod)==0)
{
cout<<0;
return 0;
}
for(int i = 0; i < 4; i ++)
{
get(nums[i].first);
for(int j = 0; j < dir.size(); j ++)
{
a[i] = (a[i] + lucas(n, dir[j], nums[i].first)) % nums[i].first;
}
}
int t = crt(fmod);
cout << qmi(q, t, mod) << endl;
return 0;
}
对于一般情况下p^a 期中a> 1的可以自行查看exlucas