http://www.cnblogs.com/BLADEVIL/p/3490321.html
http://www.cnblogs.com/zyfzyf/p/3997986.html
翻了翻题解，这两个合起来比较明白……
题意：求1~n!中与m!互质的数的数量(mod R)。
∵由欧几里得算法得gcd(a,b)=gcd(b,a%b)
∴gcd(a+b,b)=gcd(b,(a+b)%b)=gcd(b,a) 即 gcd(a,b)=gcd(a+b,b)
推广：gcd(a,b)=gcd(a+k*b,b)
∴若a与b互质，则a+k*b也与b互质
令n=5,m=4,则n!=120,m!=24
符合题意的数有以下：
5 29 53 77 101
7 31 55 79 103
11 35 59 83 107
13 37 61 85 109
17 41 65 89 113
19 43 67 91 115
23 47 71 95 119
∴若x(x<n!)与m!互质，则x+m!*k也和m!互质，根据我们刚才推出的性质，仔细观察，发现 x 在 mod m!的剩余系中恰好有n!/m!个数符合题意。
∴观察上表的第一列，我们只需考虑小于m!的数中与n!互质的数。
∴这样的数恰好有φ(m!)个
∴答案为φ(m!)*(n!/m!)
但是这样仍然超出了我们能承受的时间复杂度。
∴ans=φ(m!)*(n!/m!)
=m!*((p1-1)/p1)+((p2-1)/p2)+......+((pk-1)/pk)*n!/m! (p1......pk为m!的质因子，根据阶乘的定义，易证，m!的质因子恰好为<=m的所有质数)
=((p1-1)/p1)+((p2-1)/p2)+......+((pk-1)/pk)*n!
∴我们递推预处理出((p1-1)/p1)+((p2-1)/p2)+......+((pk-1)/pk)和n!即可O(1)地回答每个询问。
在预处理连续乘积的时候，我们要用到快速幂求乘法逆元。
综上，要预处理的是：
①maxm以内的素数
②maxm以内的素数的逆元
③1!~maxn!
④maxm以内的素数的((p1-1)/p1)+((p2-1)/p2)+......+((pk-1)/pk)
1 #include<cstdio> 2 #include<algorithm> 3 #include<iostream> 4 using namespace std; 5 typedef long long ll; 6 ll MOD; 7 int T,n[10001],m[10001],maxn,maxm,Ni[10000001],Fac[10000001],Muls[10000001]; 8 bool unPrime[10000001]; 9 void Shai_Prime() 10 { 11 unPrime[1]=true; 12 for(ll i=2;i<=maxm;i++) 13 for(ll j=i*i;j<=maxm;j+=i) unPrime[j]=true; 14 } 15 ll pow_mod(ll a,ll p,ll MOD) 16 { 17 if(!p) return 1; 18 ll ans=pow_mod(a,p>>1,MOD); 19 ans=ans*ans%MOD; 20 if((p&1)==1) ans=ans*a%MOD; 21 return ans; 22 } 23 void Init_Ni() 24 { 25 for(int i=1;i<=maxm;i++) 26 if(!unPrime[i]) 27 Ni[i]=(int)pow_mod(i,MOD-2,MOD); 28 } 29 void Init_Fac() 30 { 31 Fac[1]=1; 32 for(int i=2;i<=maxn;i++) Fac[i]=(int)((ll)Fac[i-1]*(ll)i%MOD); 33 } 34 void Init_Muls() 35 { 36 ll res=1; 37 for(int i=1;i<=maxm;i++) 38 { 39 if(!unPrime[i]) res=((res*(ll)(i-1))%MOD)*(ll)Ni[i]%MOD; 40 Muls[i]=res; 41 } 42 } 43 int main() 44 { 45 scanf("%d",&T); cin>>MOD; 46 for(int i=1;i<=T;i++) 47 { 48 scanf("%d%d",&n[i],&m[i]); 49 maxn=max(n[i],maxn); 50 maxm=max(m[i],maxm); 51 } 52 Shai_Prime(); Init_Ni(); Init_Fac(); Init_Muls(); 53 for(int i=1;i<=T;i++) 54 printf("%d ",(int)((ll)Muls[m[i]]*(ll)Fac[n[i]]%MOD)); 55 return 0; 56 }