SPOJ Prime Generator [PRIME1]

https://www.spoj.pl/problems/PRIME1/

首先想到的肯定是全部素数打表出来，但是1e9的规模，显然是TLE的。

于是，之后整个下午学习了Miller-Rabin素数测试法。

 

方法1：对于[a, b]内的每个数，进行MR判定，素数输出。特别注意的地方，由于MR是基于概率的，所以通常要多试几次，我尝试过只能在T>=3次的情况下AC，T小于3次，都不行。（据无来源消息，重复T次，出错的概率是2^(-T)，所以……）

代码如下：

#include <stdio.h>
#include <iostream>
#include <stdlib.h>
#include <time.h>

using namespace std;
typedef long long LL;

LL PowerMod(LL a, LL b, LL n)
{
    if (b == 1) return a % n;
    LL x = PowerMod(a, b >> 1, n);
    x = x * x % n;
    if (b & 1) x = x * a % n;
    return x;
}

bool Witness(LL a, LL n)
{
    LL m = n - 1;
    int q = 0;
    while((m&1) == 0)
    {
        q ++;
        m >>= 1;
    }
    LL x = PowerMod(a, m, n);
    if (x == 1 || x == n-1) return false;
    while(q --)
    {
        x = x * x % n;
        if (x == n-1) return false;
    }
    return true;
}

bool MillerRabin(int n, int T = 3)
{
    if (n < 2) return false;
    if (n == 2) return true;
    if ((n&1) == 0) return false;
    while(T --)
    {
        int a = rand() % (n-2) + 1;
        if (Witness(a, n)) return false;
    }
    return true;
}

int main()
{
    srand(unsigned(time(NULL)));
    int t;
    scanf("%d", &t);
    for (int i = 0; i < t; i ++)
    {
        int a, b;
        scanf("%d%d", &a, &b);
        i == 0 ? : puts("");
        while(a <= b)
        {
            if (MillerRabin(a))
            {
                printf("%d\n", a);
            }
            a ++;
        }
    }
    return 0;
}

 

代码跑出来之后，看下时间是1.80s，很慢。看看其他，都是1s内的，于是Google其他方法。

找到这个：http://www.algorithmist.com/index.php/SPOJ_PRIME1#Explanation

里面有解释。我翻译一下：

这题使用的是标准筛法的一个修改版本。如果用普通的筛法，对[2, n]内的每个数进行排查，明显不是高效的，会使用很多的时间和空间。然而，我们可以发现，[0, n]内没有一个合数的因子大于floor(sqrt(n))。所以，我们只需要把[2, sqrt(n)]，即[2, 31622]以内的素数筛出来。然后，对于每个询问[a, b]，使用预先筛好的素数进行第二次筛选，最后得到[a, b]内的素数，输出。

我开始没有明白最后一句，以为得到[2, sqrt(n)]的素数之后，对于每个询问[a, b]内的数n，只需要暴力测试一下。于是我写了之后，交上去，跑了2s+，比MR还慢。最后再细读Explanation，才明白。修改之后，交上去，跑了0.08ms。

代码如下：

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>

#define M 31625
#define N 100005

using namespace std;

bool Prime[M];
int PrimeList[3402];
int numP;
bool that[N];

void makePrime()
{
    int p = 2;
    while(p * p < M)
    {
        for (int i = p*p; i < M; i += p) Prime[i] = true;
        while(Prime[++p]);
    }
    for (int i = 2; i < M; i ++)
    {
        if (Prime[i] == false)
        {
            PrimeList[numP ++] = i;
        }
    }
//    cout << numP << endl;
}

int main()
{
    makePrime();
    int t;
    scanf("%d", &t);
    for (int i = 0; i < t; i ++)
    {
        i == 0 ? : puts("");
        int a, b;
        scanf("%d%d", &a, &b);
        if (a == 1) a++;

        memset(that, 0, sizeof(that));
        for (int i = 0; i < numP; i ++)
        {
            for (int j = PrimeList[i] * int(a / PrimeList[i]); j <= b; j += PrimeList[i])
            {
                if (j >= a)
                {
                    that[j - a] = true; // 把[a, b]内的合数筛掉
                }
            }
        }
        for (int i = a; i <= b; i ++)
        {
            if (that[i-a] == 0)
            {
                printf("%d\n", i);
            }
        }
    }
    return 0;
}