批改娘 10082. Fast Find Prime Number (OpenMP)

1. 題目描述
2. 輸入格式
3. 輸出格式
4. 輸入範例
5. 輸出範例
6. 提示
7. Solution
1. 7.1. 線性篩法
2. 7.2. 作弊打表

題目描述

找出 $[ L, \; R ]$ 區間內有多少個質數。

輸入格式

有多組測資，每組測資一行兩個整數 $L$, $R$，表示左右區間。

$1 \le L \le R \le 2^{31}-1$

輸出格式

對於每一組詢問輸出一行一個整數。

輸入範例

1 2	1000000 2000000 1 1000000000

輸出範例

1 2	70435 50847534

提示

Sieve of Eratosthenes

Solution

這一題沒有出得很好，後來想到可以作弊打表，目前仍沒有同學挑戰成功，略感傷心。

線性篩法

線性篩法 Accepted (1724 ms, 2816 KB)

由於只要統計 $2^{31}$ 以內的質數個數，只需要把 50000 以內的所有質數找到，利用這不到 5000 的質數進行區塊篩法即可。

由於 $[1, 2^{31}]$ 也是非常龐大的區間，可以利用 bitmask 的技術，將記憶體壓縮 32 倍，但這樣子平行度還是不高，因此若詢問區間 $[l, r]$ 有多少個質數，將區間切割成數個足夠大的 chunk，這裡決定 chunk 大約在 100000 到 1000000 都有不錯的效果。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <omp.h>
 
#define MAXL (50000>>5)+1
#define BLOCK (1<<17)
#define MAXN (BLOCK>>5)+1
#define GET1(x) (mark1[(x)>>5]>>((x)&31)&1)
#define SET1(x) (mark1[(x)>>5] |= 1<<((x)&31))
#define GET2(x) (mark2[(x)>>5]>>((x)&31)&1)
#define SET2(x) (mark2[(x)>>5] |= 1<<((x)&31))
int mark1[MAXL];
int P[5500], Pt = 0;
void sieve() {
    register int i, j, k;
    SET1(1);
    int n = 50000;
    for (i = 2; i <= n; i++) {
        if (!GET1(i)) {
            for (k = n/i, j = i*k; k >= i; k--, j -= i)
                SET1(j);
            P[Pt++] = i;
        }
    }
}
 
int local_sieve(int a, int b) {
    int sqr = (int) sqrt(b), gap = b - a;
    int mark2[MAXN] = {};
 
    for (int i = 0; i < Pt && P[i] <= sqr; i++) {
        unsigned p = P[i], base = a / p * p;
        while (base < a)    base += p;
        if (base == p)    base += p;
        for (unsigned j = base; j <= b; j += p)
            SET2(j - a);
    }
    if (a == 1)    SET2(0);
    int ret = gap + 1;
    for (int i = gap>>5; i >= 0; i--)
        ret -= __builtin_popcount(mark2[i]);
    return ret;
}
long long min(long long x, long long y) {
    return x < y ? x : y;
}
int main() {
    sieve();
    int l, r;
    while (scanf("%d %d", &l, &r) == 2) {
        int ret = 0;
        #pragma omp parallel for schedule(dynamic) reduction(+: ret)
        for (long long i = l; i <= r; i += BLOCK)
            ret += local_sieve(i, min(i+BLOCK-1, r));
        printf("%d\n", ret);
    }
    return 0;
}

作弊打表

由於問題設計只針對個數，還可以偷偷本地建表完成每一個區間的值，針對非完全包含的區間單獨計算即可。打表部分可以藉由上述的程序做到。

作弊打表 Accepted (622 ms, 384 KB)

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <omp.h>
 
#define MAXL (50000>>5)+1
#define BLOCK (1<<20)
#define MAXN (BLOCK>>5)+1
#define GET1(x) (mark1[(x)>>5]>>((x)&31)&1)
#define SET1(x) (mark1[(x)>>5] |= 1<<((x)&31))
#define GET2(x) (mark2[(x)>>5]>>((x)&31)&1)
#define SET2(x) (mark2[(x)>>5] |= 1<<((x)&31))
int preprocess[] = {82025,73586,70938,69398,68248,67307, <...>};
int mark1[MAXL];
int P[5500], Pt = 0;
void sieve() {
    register int i, j, k;
    SET1(1);
    int n = 50000;
    for (i = 2; i <= n; i++) {
        if (!GET1(i)) {
            for (k = n/i, j = i*k; k >= i; k--, j -= i)
                SET1(j);
            P[Pt++] = i;
        }
    }
}
 
int local_sieve(int a, int b) {
    int sqr = (int) sqrt(b), gap = b - a;
    int mark2[MAXN] = {};
 
    for (int i = 0; i < Pt && P[i] <= sqr; i++) {
        unsigned p = P[i], base = a / p * p;
        while (base < a)    base += p;
        if (base == p)    base += p;
        for (unsigned j = base; j <= b; j += p)
            SET2(j - a);
    }
    if (a == 1)    SET2(0);
    int ret = gap + 1;
    for (int i = gap>>5; i >= 0; i--)
        ret -= __builtin_popcount(mark2[i]);
    return ret;
}
 
long long min(long long x, long long y) {
    return x < y ? x : y;
}
long long max(long long x, long long y) {
    return x > y ? x : y;
}
int main() {
    sieve();
 
    int prepro_n = sizeof(preprocess) / sizeof(int);
    int l, r;
    while (scanf("%d %d", &l, &r) == 2) {
        int ret = 0; 
        for (int i = 0; i < prepro_n; i++) {
            long long tl = (long long) i * BLOCK + 1;
            long long tr = (long long) (i+1) * BLOCK;
            if (l <= tl && tr <= r) {
                ret += preprocess[i];
            } else {
                if (tl <= l && l <= tr) {
                    ret += local_sieve(l, min(r, tr));
                } else if (tl <= r && r <= tr) {
                    ret += local_sieve(max(tl, l), r);
                }
            }
        }
 
        printf("%d\n", ret);
    }
    return 0;
}

Morris' Blog