bzoj 2186: [Sdoi2008]沙拉公主的困惑

Time Limit: 10 Sec Memory Limit: 259 MB
Submit: 2463 Solved: 820
[Submit][Status][Discuss]

Description

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

Input

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

Output

共T行，对于每一对N，M，输出1至N！中与M！素质的数的数量对R取模后的值

Sample Input

1 11
4 2

Sample Output

1
数据范围：
对于100%的数据，1 < = N , M < = 10000000

题解：

　　这题要求在(1,N!)里与M!互质的数的个数，由于N!>M!，所以答案分成两部分，一部分是用欧拉函数phi(M!),表示1~M!里与M!互质的数的个数，还有一部分是(M!+1,N!)里的一部分数也有可能和M!互质。但是这部分数字比较难处理，现在由辗转相除法的原理可知：如果gcd(x,M!)==1，则一定有gcd(x+M!,M!)==1，所以(M!+1,N!)里的和M!互质的数一定是x+kM! (k∈Z)。所以答案就是：phi(M!)*(N!/M!)%R

　　现在的问题就是N!，M!都很大，如何求上面的式子。我们让N!%R单独算，把phi(M!)/(M!)放在一起，令f(M)=phi(M!)/M!，根据phi(x)=x*(1-1/p1)*(1-1/p2)*…*(1-1/pk) ，可以上下同时消掉M!，f(m)=(p1-1)/p1*(p2-1)/p2*…其中pi为不大于M的质数(因为由欧拉函数知pi是素因子，而M!=1*2*3*4***M，所以小于M的素数就是素因子)，所以对于f(i)，如果i是质数f(i)=f(i-1)*(i-1)/i，否则f(i)=f(i-1)。

　　还有一个问题，就是f(i)=f(i-1)*(i-1)/i 在模R意义下存在除法，这就有可能出现分数，这就需要用到乘法逆元，把除法转化成乘法。具体讲解：乘法逆元用扩展欧几里得求解乘法逆元

　　根据以上关系式可以预处理f(1)~f(10^7)。每次询问只需要输出f(M)*N!%R即可。

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 typedef long long LL;
 4 const int maxn=1e7+5;
 5 int N,M,P,T;
 6 LL fac[maxn],ans[maxn];
 7 bool f[maxn];
 8 inline void exgcd(int a,int b,int &x,int &y){
 9     if(b==0){
10         x=1; y=0; 
11         return ;
12     }
13     exgcd(b,a%b,x,y);
14     int t=x;
15     x=y; y=t-a/b*x;
16 }
17 inline int getinv(int a){
18     int x=0,y=0;
19     exgcd(a,P,x,y);
20     return (x%P+P)%P;
21 }
22 int main(){
23     scanf("%d%d",&T,&P);
24     fac[1]=1;
25     for(int i=2;i<=1e7;i++) fac[i]=fac[i-1]*i%P;//预处理 N!%R 
26     ans[1]=1;
27     for(int i=2;i<=1e7;i++){
28         if(f[i]==false){// 是素数 
29             ans[i]=ans[i-1]*(i-1)%P*getinv(i)%P;
30             for(int j=2;j<=1e7/i;j++) f[i*j]=true;//筛素数 
31         }
32         else ans[i]=ans[i-1];
33     }
34     while(T--){
35         scanf("%d%d",&N,&M);
36         printf("%lld
",ans[M]*fac[N]%P);
37     }
38     return 0;
39 }