一种无双筛现象的乘筛质数法

最新推荐文章于 2025-05-27 21:52:10 发布

harriszhuang

最新推荐文章于 2025-05-27 21:52:10 发布

阅读量289

点赞数

分类专栏：数论文章标签：质数筛法算法 Basic

本文链接：https://blog.youkuaiyun.com/harriszhuang/article/details/100515435

版权

数论专栏收录该内容

1 篇文章

订阅专栏

以往的质数筛都采用除筛法，通过逐个测试 $[3,N]$ 对 $[2, \sqrt N]$ 的可除性来确认一个数是不是质数。
众所周知，计算机除法是一系列减法的组合，运算速度较慢。
当计算机具有足够大的内存时，可以采用乘筛法，方法是先定义一个一维旗标数组，当做一条有限长的纸带，它的下标 $[1,N]$ 的每个元素对应一个自然数的状态。从 $I = 2$ 开始逐个筛去因数的倍数 $I \times J (J \geqslant 2)$ ，将被筛数所在的旗标置1（打孔），表示已经被验证为合数，直到因数 $I = [2, \sqrt N]$ 全部被筛完，这时纸带上未被打孔的就是质数。然后枚举旗标数组，将未被标记为合数的数字作为质数显示出来，例如：

最简单的办法是将 $[2, \sqrt N]$ 全部筛一遍，例如当纸带的长度是 $N=16$ 时，可将 $I=\{ 2,3,4 \}$ 全部筛一遍，这种方法叫做盲筛，即不管被筛数是否被筛过，只要满足因数的整数倍条件就进行筛除操作。这种筛法是完备有效的，但当我们采用手工操作时就会发现，某些数字至少被筛了两遍，出现了双筛，例如当 $N=16$ 时：

筛去1：

筛去2的倍数： $2 \times 2 = 4$ ， $2 \times 3 = 6$ ， $2 \times 4 = 8$ ， $2 \times 5 = 10$ ， $2 \times 6 = 12$ ， $2 \times 7 = 14$ ， $2 \times 8 = 16$

筛去3的倍数： $3 \times 2 = 6$ （已经被筛）， $3 \times 3 = 9$ ， $3 \times 4 = 12$ （已经被筛）， $3 \times 5 = 15$

筛去4的倍数： $4 \times 2 = 8$ （已经被筛）， $4 \times 3 = 12$ （已经被筛）， $4 \times 4 = 16$ （已经被筛）

（未完待续）

#INCLUDE "VBCOMPAT.BI"

'定义纸带长度
CONST PAPER_LEN AS INTEGER = 121

'定义合数指示旗标
DIM FLAG() AS INTEGER
DIM FLAG_TMP() AS INTEGER

'定义循环变量
DIM I AS INTEGER
DIM J AS INTEGER
DIM N AS INTEGER
DIM A AS DOUBLE
DIM B AS DOUBLE
DIM K AS INTEGER

REDIM FLAG(PAPER_LEN * 2 + 1)
REDIM FLAG_TMP(PAPER_LEN * 2 + 1)

'划掉第一个数
FLAG(1) = 1

N = PAPER_LEN

PRINT "PRIME v0.01"
PRINT "PAPER_LEN = "; N
PRINT "SQR(N) = "; SQR(N)

A = NOW
PRINT "BEGIN = "; FORMAT(A, "yy/mm/dd hh:mm:ss")

'划掉偶数
J = 2 
I = 2
K = 1
DO WHILE I < N
  I = I + J
  FLAG(I) = 1
  K = K + 1
LOOP

FOR I=3 TO SQR(N) STEP 2
  IF FLAG(I) = 0 THEN
    FOR J=I TO N/I STEP 2
      IF FLAG(J) = 0 THEN
        FLAG_TMP(I*J) = 1
        K = K + 1
      END IF
    NEXT J
    FOR J=I TO N STEP 2
      IF FLAG_TMP(J) = 1 THEN
        FLAG(J) = 1 : FLAG_TMP(J) = 0
      END IF
    NEXT J
  END IF
NEXT I

FOR I=1 TO N
  IF FLAG(I) = 0 THEN
    PRINT I,
  END IF
NEXT I

B = NOW
PRINT "OPERATIONS = "; K
PRINT "END = "; FORMAT(B, "yy/mm/dd hh:mm:ss")

END