【容斥原理】AOJ-557 Redraiment猜想2

Redraiment猜想2

Time Limit: 1000 ms   Case Time Limit: 1000 ms   Memory Limit: 64 MB

Description

redraiment在家极度无聊，于是找了张纸开始统计素数的个数。 
设函数f(n)返回从1->n之间素数的个数。 
redraiment发现: 

f(1) = 0 
f(10) = 4 
f(100) = 25 
... 

满足g(m) = 17 * m2 / 3 - 22 * m / 3 + 5 / 3 
其中m为n的位数。 
他很激动，是不是自己发现了素数分布的规律了！ 
请你设计一个程序，求出1->n范围内素数的个数，来验证redraiment是不是正确的，也许还可以得诺贝尔奖呢。^_^

此题题面与前面的《Redraiment猜想》类似，不过数据范围更大，你能想出更快的算法吗？

此题诞生的故事：
N年前Tatsuya开始从id=1002（新OJ的非A+B的第1题）开始向后做题，没想到就在1007这个经典的求素数相关问题上卡住了，想啊想，正好上某门课不久，某个理论正好用上，成功解决本题。后来成了OJ管理员看到了数据，悲剧，与题目不符啊，拿现在的话来说XX啊。于是去别的网上找原题，果然题目上数据错了，虽然觉得有些对不起之前按照题面做的同学，不过还是把题面改“正确”了。不过这样一来，浪费了一个思考的好机会，所以时隔N年，此题又再度出山了。

Input

输入包括多组数据。 
每组数据仅有一个整数n (1≤n≤10 ⁸)。 
输入以0结束

Output

对于每组数据输入，输出一行，为1至n(含)之间的素数的个数

Sample Input

Original

Transformed

1
10
65
100
0

Sample Output

Original

Transformed

0
4
18
25

————————————————————困惑的分割线————————————————————

思路：做完了AOJ的Redraiment猜想1之后，来挑战这道题。。。。。。可是要学会容斥原理根本不那么简单。

首先我们来线性筛选素数打出10000以内的素数表，正好用来计算10^8以内素数的个数。

以下是误导思路  要想看正解，跳到容斥原理部分

接下来是问题的核心：倘若使用线性筛法，要么会MLE要么会TLE一般情况下，都是MLE因为就算是bool型的vis[100000000]数组，都会超过内存限制。为此我还专门学习了更加节省空间开销的办法：用int来储存32个数字的标记。方法如下：

#define POS(a) (1 << (a & ((1<<5)-1)))
int vis[3125000]， prime[5761460];
for(int i = 2; i < MAXN; i++) {
        if(!(vis[i>>5] & POS(i)))   prime[cnt++] = i;
        for(int j = 0; j < cnt && prime[j]*i < MAXN; j++) {
            index = prime[j]*i;
            vis[index>>5] |= POS(index);
            if(i % prime[j] == 0)   break;
        }
        if(i == 50000000)   prime[0] = 2;
    }

要弄懂这个并不简单，记住就更难了。首先是：vis[i>>5]的含义，即vis[i / 32]，vis[ ]数组开了3125000正好是100000000^0.5，这10^8个数，我分成每32个数字储存在一个vis[ i ]当中。查看标记的时候，对vis[ i ]的第k位进行判断即可。这就是宏定义POS(a)的意义。((1<<5)-1)))即32-1 = 31。即11111。

用 a & 11111 得到a / 32的余数。注意！这是代替Mod的方法：

a & ((1<<k)-1) <--> a Mod 2^k

然后 1<<... 得到a的位置(position)，对该位置进行操作可以：1.判断是否有标记 2.做标记

完成了位压缩之后，提交发现TLE。

只能学习容斥原理了。只看2、3、5。首先画出Venn图：2的倍数，3的倍数，5的倍数。假设这三个圈都是有交集的。我们知道质数的倍数是合数。如何去除合数呢？现在用10减去2的倍数，3的倍数，5的倍数这些数的个数。我们发现2和3的交集减了两次，3和5的交集减了两次，2和5的交集减了两次。2、3、5的交集减了三次。那么，我们加上2、3的交集，2、5的交集，3、5的交集，这时，2、3、5的倍数又被加了回来，再把它减掉一次。我们得到了10-5-3-2+1+1+0-0=2。因为素数2、3、5也被我们减掉了，所以加回来2+3=5。但是1不是素数，减去。5-1=4。即正确答案。

要实现生成子集并且对子集进行操作，利用dfs。该dfs需要的参数是什么？容斥的方法是用哪几个素数的倍数来筛选，那么参数就是：几个（num_pri）；素数(index)；的倍数(f_prod)。必须有剪枝，否则时间复杂度难以接受。我们用sqrt(10^8)以内的素数实现容斥原理，那么对于大于所输入n的素数之积不要再dfs了。for()循环怎么写呢？

参考我即将学习的“子集生成”。

代码如下：

#include <stdio.h>
#include <math.h>
const int N = 10005;
bool vis[N];
int cnt, sqrn, n;
int p[1500];
void get_prime(){ //得到10000以内的素数表p[]
    int k = 0;
    vis[1] = 1;
    for(int i = 2; i < N; i++){
        if(!vis[i])   p[k++] = i;
        for(int j = 0; j < k && p[j]*i < N; j++){ //注意！不论此数是不是素数，都要用，不能写else
            vis[p[j]*i] = 1;
            if(i % p[j] == 0)   break;
        }
    }
}
//num_pri是当前使用的素数个数，f_prod是本层前k个素数的积，index是使用的素数的下标
void dfs(int num_pri, int f_prod, int index){//容斥原理dfs+剪枝
    for(int i = index; p[i] <= sqrn; i++){ //从当前使用的素数开始，用sqrt(n)之前的素数进行筛选
        if(1LL*f_prod*p[i] > n)    return ;//强力剪枝，再乘一个素数就会超过范围，注意int的溢出
        dfs(num_pri+1, f_prod*p[i], i+1);//！！注意，进入下一层dfs，而在此之前，要完成本层dfs，即先完成for（）循环
        if(num_pri & 1)//奇数个素数，此时意味着减去这些个合数
            cnt -= n / (f_prod * p[i]);
        else//偶数个素数，把多减去的合数的个数加回来
            cnt += n / (f_prod * p[i]);
        if(num_pri == 1)    cnt++;//单独的素数本身也被减去了，加回1
    }
}
int main(){
    get_prime();
    while(scanf("%d", &n), n){
        cnt = n;
        sqrn = sqrt((double)n);//sqrt函数参数必须是double
        dfs(1, 1, 0);
        printf("%d\n", cnt - 1);//1没有被筛掉
    }
	return 0;
}