[数论][SDOI]沙拉公主的疑惑

题目传送门:https://www.luogu.org/problemnew/show/P2155

---

题目描述

大富翁国因为通货膨胀，以及假钞泛滥，政府决定推出一项新的政策：现有钞票编号范围为1到N的阶乘，但是，政府只发行编号与M!互质的钞票。房地产第一大户沙拉公主决定预测一下大富翁国现在所有真钞票的数量。现在，请你帮助沙拉公主解决这个问题，由于可能张数非常大，你只需计算出对R取模后的答案即可。R是一个质数。

输入格式：
第一行为两个整数T，R。R<=10^9+10，T<=10000，表示该组中测试数据数目，R为模 后面T行，每行一对整数N，M，见题目描述 m<=n

输出格式：
共T行，对于每一对N，M(1<=N,M<=10000000)，输出1至N！中与M！素质的数的数量对R取模后的值

---

题解

本题目的大意就是求1到N!中有多少个数与M!互质。
由于互质，我们很容易就想到欧拉函数，所以本题的解就可以被分为两部分：
(1,M!)这个区间中与M!互质的数+(M!,N!)这个区间中与M!互质的数
第一个部分很容易求，就是φ(M!)
第二个部分就比较麻烦了，首先，根据GCD(a,b)=GCD(b,a%b)辗转相除的原理可知：
假如GCD(x,M!)==1,那么GCD(x+K*M!,M!)==0.(K为整数)
所以，在(M!,N!)这个区间里与M!互质的数就是x+K*M!
而K可以有多大呢？因为N>M,所以N!与M!肯定是成倍数关系，所以K=N!/M!
所以最后两个区间的答案相加可以得到$φ(M!)*K\%R$.
将K代换得

$$φ(M!)*\frac{N!}{M!}\%R$$
现在关键就是如何求这个表达式的值。
我们先把表达式的运算顺序改变一下: $$(N!)\%R*\frac{φ(M!)}{M!}$$
根据欧拉函数的运算式 $$φ(x)=x*(1-\frac{1}{p_1})(1-\frac{1}{p_2})…(1-\frac{1}{p_k})$$ ，(φM!)/M!中M!可以被抵消，所以还剩后面的(1-1/p1)(1-1/p2)…(1-1/pk)。(其中p1,p2,p3…..pk为x的质因子)
将除号放到括号外得 $$f(m)=\frac{p_1-1}{p_1}*\frac{p_2-1}{p_2}.....*\frac{p_k-1}{p_k}$$
由于 $M!=1*2*3*...*M$
所以小于M的质数就是它所有的质因子。
所以对于f(i)来说，当i为质数， $f(i)=f(i-1)*\frac{i-1}i\%R$
否则， $f(i)=f(i-1)$.
因为这里出现了模R意义下的除法，我们需要用到乘法逆元。
最后，转换为f(i)=f(i-1)*(i-1)再乘i的逆元%R。
(如何求逆元就不再赘述了，如果不会的话，自行百度 awa)

可是对如果每个询问都处理一次，肯定超时，所以在处理询问前先进行预处理，算出每个N!模R的值，用线性筛求出所有素数，用拓展欧几里得定理求出所有素数在%R意义下的逆元，最后预处理出所有的f(i),这样，对于每一次询问就可以在O(1)的时间内完成。

---

代码

#include<cmath>
#include<ctime>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdio>
#define inf 1000000000ll
#define N 10000000
#define ll long long
using namespace std;
ll T,R,n,m,cnt;
ll fac[N+5],inverse[N+5],prime[N+5],ans[N+5];
bool mark[N+5];
void exgcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0)
    {
        x=1;
        y=0;
        return;
    }
    exgcd(b,a%b,x,y);
    ll t=x;
    x=y;
    y=t-a/b*y;
}
ll getinverse(ll t)
{
    ll x,y;
    exgcd(t,R,x,y);
    return (x%R+R)%R;
}
void pre()
{
    fac[1]=1;
    for(ll i=2;i<=N;i++) fac[i]=(ll)fac[i-1]*i%R;
    inverse[1]=1;
    for(ll i=2;i<=N;i++)
    {
        if(!mark[i]) 
        {
            prime[++cnt]=i;
            inverse[i]=getinverse(i);
        }
        for(ll j=1;prime[j]*i<=N&&j<=cnt;j++)
        {
            mark[prime[j]*i]=1;
            if(i%prime[j]==0) break;
        }
    }
    ans[1]=1;
    for(ll i=2;i<=N;i++)
    {
        ans[i]=ans[i-1];
        if(!mark[i]) ans[i]=(ll)ans[i]*(i-1)%R*inverse[i]%R;
    }
}
int main()
{
    scanf("%lld%lld",&T,&R);
    pre();
    while(T--)
    {
        scanf("%lld%lld",&n,&m);
        printf("%lld\n",(ll)(fac[n]*ans[m])%R);
    }
    return 0;
}