【数论】【欧拉函数】【筛法求素数】【乘法逆元】【快速幂取模】bzoj2186 [Sdoi2008]沙拉公主的困惑

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 }