[SDOI2008]沙拉公主的困惑

题目描述

大富翁国因为通货膨胀，以及假钞泛滥，政府决定推出一项新的政策：现有钞票编号范围为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！素质的数的数量对R取模后的值

输入输出样例

输入样例#1：

1 11
4 2

输出样例#1：

1

数据范围：
对于100%的数据，1 < = N , M < = 10000000
题解：
素数筛法，同余，欧拉函数，逆元

首先，我们来引出一个定理

如果a与b互质，那么 $b*k+a$ 也与b互质。证明和证明gcd的证明类似。

反过来，我们也可以用 $gcd$ 证明，

因为 $gcd(a,b)=1$ ，所以 $gcd(a\%b,b)=1$

因为 $a\%b=a-k*b$ ，故 $gcd(a-k*b,b)=1$ ,及 $a-k*b$ 与 $b$ 互质。

根据这个特性，并且 $n>=m$ ,所以可以将n！分成若干段，每段为m!，每一段中与m!互质的个数都是相等的且等于1到m！中与m！互质的个数

我们可以得到式子 $ans={frac{n!}{m!}*phi(m!)}$

进一步拆开，我们可以得到（假设p为m!的质因数，很容易可以知道，p就是所有小于m的素数，r为质因数个数）

$ans={frac{n!}{m!}*m!*frac{prod limits_{i=1}^{r}(p_i-1)}{prodlimits_{i=1}^{r}p_i } o ans=n!*frac{prod limits_{i=1}^{r}(p_i-1)}{prodlimits_{i=1}^{r}p_i } }$

因为 $ans$ 要 $mod R$ ,所以我们也要算1到m的逆元，在累乘 $prodlimits_{i=1}^{r}p_i$ ，乘的是 $p_i$ 的逆元。

本题常数很大，有如下解决方法：
1.用int数组，计算时转为long long在取余，否则空间超限
2.离线处理,将最大的m，n记下，事先将素数筛出O(maxm)
3.线性时间求出1~maxn间的逆元，再O(maxn)求出阶乘
4.用类似前缀和的思想(称为前缀积)，将x之前所有素数的逆元，pi-1的积存起，记为AA[x],BB[x]
5.求逆元用一层循环，其他(pi-1)，逆元的前缀积，阶乘用一层循环搞定，查询时直接调用。
B[]为阶乘，A[]为逆元，AA[]为(素数-1)的前缀积，BB[]为逆元的前缀积

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cstring>
 4 using namespace std;
 5 long long s;
 6 int B[10000011],A[10000011],AA[10000011],BB[10000011];
 7 bool vis[10000011];
 8 int prime[1000011],Mod,n[10001],m[10001],maxm,maxn,tot;
 9 int main()
10 {int i,j,l,T,k;
11     cin>>T>>Mod;
12     for (l=1;l<=T;l++)
13     {scanf("%d%d",&n[l],&m[l]);
14     if (m[l]>maxm) maxm=m[l];
15     if (n[l]>maxn) maxn=n[l];
16     }
17     for (i=2;i<=maxm;i++)
18     {
19         if (vis[i]==0) 
20         {
21             tot++;
22             prime[tot]=i;
23         }
24         for (j=1;j<=tot;j++)
25         {
26             if (i*prime[j]>maxm) break;
27             vis[i*prime[j]]=1;
28             if (i%prime[j]==0) break;
29         }
30     }
31      A[1]=1;
32     for (i=2;i<=prime[tot];i++)
33      A[i]=((Mod-(Mod/i))*(long long)A[Mod%i])%Mod;
34     B[1]=1;
35     k=1;
36      AA[1]=1;BB[1]=1;
37     for (i=2;i<=maxn;i++)
38          {
39             B[i]=((long long)B[i-1]*(long long)i)%Mod;
40             if (i<=maxm)
41             {
42             AA[i]=1;
43             BB[i]=1;
44             if (prime[k]==i)
45             {
46                 AA[i]=(prime[k]-1)%Mod;
47                 BB[i]=A[i]%Mod;
48                  k++;        
49             }
50             AA[i]=((long long)AA[i]*(long long)AA[i-1])%Mod;
51             BB[i]=((long long)BB[i]*(long long)BB[i-1])%Mod;
52             }
53          } 
54     for (l=1;l<=T;l++)
55     {
56          s=B[n[l]];
57          s=(s*(long long)AA[m[l]])%Mod,s=((s*(long long)BB[m[l]]+Mod))%Mod;
58     printf("%d
",(s+Mod)%Mod);
59     }
60 }