——曾晓奇
关于素数的算法是信息学竞赛和程序设计竞赛中常考的数论知识，在这里我跟大家讲一下寻找一定范围内素数的几个算法。看了以后相信
对大家一定有帮助。
正如大家都知道的那样，一个数 n 如果是合数，那么它的所有的因子不超过sqrt(n)--n的开方,那么我们可以用这个性质用最直观的方法
来求出小于等于n的所有的素数。
num = 0;
for(i=2; i<=n; i++)
{ for(j=2; j<=sqrt(i); j++)
if( j%i==0 ) break;
if( j>sqrt(i) ) prime[num++] = i; //这个prime[]是int型，跟下面讲的不同。
}
这就是最一般的求解n以内素数的算法。复杂度是o(nsqrt(n))，如果n很小的话，这种算法（其实这是不是算法我都怀疑，没有水平。当然没
接触过程序竞赛之前我也只会这一种求n以内素数的方法。-_-~）不会耗时很多
但是当n很大的时候，比如n=10000000时，nsqrt(n)>30000000000,数量级相当大。在一般的机子它不是一秒钟跑不出结果，它是好几分钟都跑不
出结果，这可不是我瞎掰的，想锻炼耐心的同学不妨试一试~。。。。
在程序设计竞赛中就必须要设计出一种更好的算法要求能在几秒钟甚至一秒钟之内找出n以内的所有素数。于是就有了素数筛法。
（我表达得不清楚的话不要骂我，见到我的时候扁我一顿我不说一句话。。。）
素数筛法是这样的：
1开一个大的bool型数组prime[]，大小就是n+1就可以了先把所有的下标为奇数的标为true,下标为偶数的标为false
2然后：
for( i=3; i<=sqrt(n); i+=2 )
{ if(prime[i])
for( j=i+i; j<=n; j+=i ) prime[j]=false;
}
3最后输出bool数组中的值为true的单元的下标，就是所求的n以内的素数了。
原理很简单，就是当i是质(素)数的时候，i的所有的倍数必然是合数。如果i已经被判断不是质数了，那么再找到i后面的质数来把这个质
数的倍数筛掉。
一个简单的筛素数的过程：n=30。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
第 1 步过后2 4 28 30这15个单元被标成false,其余为true。
第 2 步开始：
i=3; 由于prime[3]=true, 把prime[6], [9], [12], [15], [18], [21], [24], [27], [30]标为false
i=4; 由于prime[4]=false,不在继续筛法步骤。
i=5; 由于prime[5]=true, 把prime[10],[15],[20],[25],[30]标为false
i=6>sqrt(30)算法结束。
第 3 步把prime[]值为true的下标输出来：
for(i=2; i<=30; i++)
if(prime[i]) printf("%d ",i);
结果是 2 3 5 7 11 13 17 19 23 29
这就是最简单的素数筛选法，对于前面提到的10000000内的素数，用这个筛选法可以大大的降低时间复杂度。把一个只见黑屏的算法
优化到立竿见影，一下就得到结果。关于这个算法的时间复杂度，我不会描述，没看到过类似的记载。只知道算法书上如是说：前几年比
较好的算法的复杂度为o(n),空间复杂度为o(n^(1/2)/logn)另外还有时间复杂度为o(n/logn),但空间复杂度为O(n/(lognloglogn))的算法。
我水平有限啦，自己分析不来。最有说服力的就是自己上机试一试。下面给出这两个算法的程序：
//最普通的方法：
#include<stdioh>
#include<mathh>
#define N 10000001
int prime[N];
int main()
{
int i, j, num = 0;
for(i=2; i<N; i++)
{ for(j=2; j<=sqrt(i); j++)
if( j%i==0 ) break;
if( j>sqrt(i) ) prime[num++] = i;
}
for(i=2; i<100; i++) //由于输出将占用太多io时间，所以只输出2-100内的素数。可以把100改为N
if( prime[i] )printf("%d ",i);
return 0;
}
//用了筛法的方法：
#include<stdioh>
#include<mathh>
#define N 10000001
bool prime[N];
int main()
{
int i, j;
for(i=2; i<N; i++)
if(i%2) prime[i]=true;
else prime[i]=false;
for(i=3; i<=sqrt(N); i++)
{ if(prime[i])
for(j=i+i; j<N; j+=i) prime[i]=false;
}
for(i=2; i<100; i++)//由于输出将占用太多io时间，所以只输出2-100内的素数。可以把100改为N
if( prime[i] )printf("%d ",i);
return 0;
}
装了vc的同学上机跑一下这两个程序试一试。这个差别，绝对是天上地下。前面那个程序绝对是n分钟黑屏的说。
另外，对于这样的筛法，还可以进一步优化，就是bool型数组里面只存奇数不存偶数。如定义prime[N],则0表示
3，1表示5，2表示7，3表示9。如果prime[0]为true,则表示3时素数。prime[3]为false意味着9是合数。
这样的优化不是简单的减少了一半的循环时间，比如按照原始的筛法，数组的下标就对应数。则在计算30以内素
数的时候3个步骤加起来走了15个单位时间。但是用这样的优化则是这样：
则由于只存3 5 7 9 11 13 15 17 19 21 23 25 27 29，只需要14个单元
第 1 步 把14个单元赋为true (每个单元代表的数是2i+3,如第0单元代表3，第1单元代表5)
第 2 步开始：
i=0; 由于prime[0]=true, 把 [3], [6], [9], [12]标为false
i=1; 由于prime[1]=true, 把 [6], [11]标为false
i=2 2i+3>sqrt(30)算法结束。
这样优化以后总共只走6个单位时间。
当n相当大以后这样的优化效果就更加明显，效率绝对不仅仅是翻倍。
出了这样的优化以外，另外在每一次用当前已得出的素数筛选后面的数的时候可以一步跳到已经被判定不是素数的
数后面，这样就减少了大量的重复计算。（比如我们看到的，i=0与i=1时都标了[6],这个就是重复的计算。）
我们可以发现一个规律，那就是3（即i=0）是从下标为[3]的开始筛的，5（即i=1）是从下标为[11]开始筛的（因为[6]
已经被3筛过了）。然后如果n很大的话，继续筛。7（i=2）本来应该从下标为[9]开始筛，但是由于[9]被筛过了，而
[16]也已经被5（i=1）筛过了。于是7（i=2）从[23]（就是223+3=49）开始筛。
于是外围循环为i时，内存循环的筛法是从 i+(2i+3)(i+1)即i(2i+6)+3开始筛的。
这个优化也对算法复杂度的降低起到了很大的作用。
相比于一般的筛法，加入这两个优化后的筛法要高效很多。高兴去的同学可以试着自己编写程序看一看效率。我这里
有程序，需要的可以向我要。不懂得也可以问我。
上面的素数筛法是所有程序设计竞赛队员都必须掌握的，而后面加了两个优化的筛法是效率很高的算法，是湖南大学
huicpc39同学设计的（可能是学来的，也可能是自创的。相当强悍）。在数量级更大的情况下就可以发现一般筛法和
优化后的筛法的明显区别。
另外，台湾的ACMTino同学也给我介绍了他的算法：a是素数，则下一个起点是aa,把后面的所有的aa+2ia筛掉。
这上面的所有的素数筛选的算法都可以再进一步化为二次筛选法，就是欲求n以内的素数，就先把sqrt(n)内的素数求
出来，用已经求得的素数来筛出后面的合数。
我把一般的筛选法的过程详细的叙述了一遍，应该都懂了吧？后面的优化过程及不同的方法，能看懂最好。不是很难的。
相关知识：
最大公约数只有1和它本身的数叫做质数（素数）——这个应该知道吧？-_-b
至今为止，没有任何人发现素数的分布规律，也没有人能用一个公式计算出所有的素数。关于素数的很多的有趣的性质或者科学家的努力
我不在这里多说，大家有兴趣的话可以到百度或google搜一下。我在下面列出了一个网址，上面只有个大概。更多的知识需要大家一点一点
地动手收集。
http://wwwscitomcomcn/discovery/universe/home01html
1高斯猜测，n以内的素数个数大约与n/ln(n)相当，或者说，当n很大时，两者数量级相同。这就是著名的素数定理。
2十七世纪费马猜测，2的2^n次方+1，n=0，1，2…时是素数，这样的数叫费马素数，可惜当n=5时，2^32+1就不是素数，
至今也没有找到第六个费马素数。
318世纪发现的最大素数是2^31-1，19世纪发现的最大素数是2^127-1，20世纪末人类已知的最大素数是2^859433-1，用十进制表示，这是一个258715位的数字。
4孪生素数猜想：差为2的素数有无穷多对。目前知道的最大的孪生素数是1159142985×2^2304－1和1159142985×2^2304＋1。
5歌德巴赫猜想：大于2的所有偶数均是两个素数的和，大于5的所有奇数均是三个素数之和。其中第二个猜想是第一个的自然推论，因此歌德巴赫猜想又被称为1+1问题。我国数学家陈景润证明了1+2，即所有大于2的偶数都是一个素数和只有两个素数因数的合数的和。国际上称为陈氏定理。
本程序由国外的Vulture大哥编写，并公布了源码，这个是他95年的一个作品。真的很佩服他！太强了。好好学习啊！！！！
;编译方法： 1 tasm 3dasm
; 2 tlink 3dobj
; 3 exe2bin 3dexe 3dcom
; Assembler Program By Vulture ;
; 3D-system example Use the following formulas to rotate a point: ;
; ;
; Rotate around x-axis ;
; YT = Y COS(xang) - Z SIN(xang) / 256 ;
; ZT = Y SIN(xang) Z COS(xang) / 256 ;
; Y = YT ;
; Z = ZT ;
; ;
; Rotate around y-axis ;
; XT = X COS(yang) - Z SIN(yang) / 256 ;
; ZT = X SIN(yang) Z COS(yang) / 256 ;
; X = XT ;
; Z = ZT ;
; ;
; Rotate around z-axis ;
; XT = X COS(zang) - Y SIN(zang) / 256 ;
; YT = X SIN(zang) Y COS(zang) / 256 ;
; X = XT ;
; Y = YT ;
; ;
; Divide by 256 coz we have multiplyd our sin values with 256 too ;
; This example isn’t too fast right now but it’ll work just fine ;
; ;
; Current Date: 6-9-95 Vulture ;
; ;
IDEAL ; Ideal mode
P386 ; Allow 80386 instructions
JUMPS ; Tasm handles out of range jumps (rulez!:))
SEGMENT CODE ; Code segment starts
ASSUME cs:code,ds:code ; Let cs and ds point to code segment
ORG 100h ; Make a COM file
START: ; Main program
mov ax,0013h ; Init vga
int 10h
mov ax,cs
mov ds,ax ; ds points to codesegment
mov ax,0a000h
mov es,ax ; es points to vga [Page]
lea si,[Palette] ; Set palette
mov dx,3c8h
xor al,al
out dx,al
mov dx,3c9h
mov cx,1893
repz outsb
; === Set some variables ===
mov [DeltaX],1 ; Initial speed of rotation
mov [DeltaY],1 ; Change this and watch what
mov [DeltaZ],1 ; happens It’s fun!
mov [Xoff],256
mov [Yoff],256 ; Used for calculating vga-pos
mov [Zoff],300 ; Distance from viewer
MainLoop:
call MainProgram ; Yep do it all ;-)
in al,60h ; Scan keyboard
cmp al,1 ; Test on ESCAPE
jne MainLoop ; Continue if not keypressed
; === Quit to DOS ===
mov ax,0003h ; Back to textmode
int 10h
lea dx,[Credits]
mov ah,9
int 21h
mov ax,4c00h ; Return control to DOS
int 21h ; Call DOS interrupt
; === Sub-routines ===
PROC WaitVrt ; Waits for vertical retrace to reduce \"snow\"
mov dx,3dah
Vrt:
in al,dx
test al,8
jnz Vrt ; Wait until Verticle Retrace starts
NoVrt:
in al,dx
test al,8
jz NoVrt ; Wait until Verticle Retrace ends
ret ; Return to main program
ENDP WaitVrt
PROC UpdateAngles
; Calculates new x,y,z angles
; to rotate around
mov ax,[XAngle] ; Load current angles
mov bx,[YAngle]
mov cx,[ZAngle]
add ax,[DeltaX] ; Add velocity
and ax,11111111b ; Range from 0255
mov [XAngle],ax ; Update X
add bx,[DeltaY] ; Add velocity
and bx,11111111b ; Range from 0255
mov [YAngle],bx ; Update Y
add cx,[DeltaZ] ; Add velocity
and cx,11111111b ; Range from 0255
mov [ZAngle],cx ; Update Z
ret
ENDP UpdateAngles
PROC GetSinCos
; Needed : bx=angle (0255)
; Returns: ax=Sin bx=Cos
push bx ; Save angle (use as pointer)
shl bx,1 ; Grab a word so bx=bx2 [Page]
mov ax,[SinCos bx] ; Get sine
pop bx ; Restore pointer into bx
push ax ; Save sine on stack
add bx,64 ; Add 64 to get cosine
and bx,11111111b ; Range from 0255
shl bx,1 ; 2 coz it’s a word
mov ax,[SinCos bx] ; Get cosine
mov bx,ax ; Save it bx=Cos
pop ax ; Restore ax=Sin
ret
ENDP GetSinCos
PROC SetRotation
; Set sine & cosine of x,y,z
mov bx,[XAngle] ; Grab angle
call GetSinCos ; Get the sine&cosine
mov [Xsin],ax ; Save sin
mov [Xcos],bx ; Save cos
mov bx,[Yangle]
call GetSinCos
mov [Ysin],ax
mov [Ycos],bx
mov bx,[Zangle]
call GetSinCos
mov [Zsin],ax
mov [Zcos],bx
ret
ENDP SetRotation
PROC RotatePoint ; Rotates the point around x,y,z
; Gets original x,y,z values
; This can be done elsewhere
movsx ax,[Cube si] ; si = X (movsx coz of byte)
mov [X],ax
movsx ax,[Cube si 1] ; si 1 = Y
mov [Y],ax
movsx ax,[Cube si 2] ; si 2 = Z
mov [Z],ax
; Rotate around x-axis
; YT = Y COS(xang) - Z SIN(xang) / 256
; ZT = Y SIN(xang) Z COS(xang) / 256
; Y = YT
; Z = ZT
mov ax,[Y]
mov bx,[XCos]
imul bx ; ax = Y Cos(xang)
mov bp,ax
mov ax,[Z]
mov bx,[XSin]
imul bx ; ax = Z Sin(xang)
sub bp,ax ; bp = Y Cos(xang) - Z Sin(xang)
sar bp,8 ; bp = Y Cos(xang) - Z Sin(xang) / 256
mov [Yt],bp
mov ax,[Y]
mov bx,[XSin]
imul bx ; ax = Y Sin(xang)
mov bp,ax
mov ax,[Z]
mov bx,[XCos]
imul bx ; ax = Z Cos(xang)
add bp,ax ; bp = Y SIN(xang) Z COS(xang) [Page]
sar bp,8 ; bp = Y SIN(xang) Z COS(xang) / 256
mov [Zt],bp
mov ax,[Yt] ; Switch values
mov [Y],ax
mov ax,[Zt]
mov [Z],ax
; Rotate around y-axis
; XT = X COS(yang) - Z SIN(yang) / 256
; ZT = X SIN(yang) Z COS(yang) / 256
; X = XT
; Z = ZT
mov ax,[X]
mov bx,[YCos]
imul bx ; ax = X Cos(yang)
mov bp,ax
mov ax,[Z]
mov bx,[YSin]
imul bx ; ax = Z Sin(yang)
sub bp,ax ; bp = X Cos(yang) - Z Sin(yang)
sar bp,8 ; bp = X Cos(yang) - Z Sin(yang) / 256
mov [Xt],bp
mov ax,[X]
mov bx,[YSin]
imul bx ; ax = X Sin(yang)
mov bp,ax
mov ax,[Z]
mov bx,[YCos]
imul bx ; ax = Z Cos(yang)
add bp,ax ; bp = X SIN(yang) Z COS(yang)
sar bp,8 ; bp = X SIN(yang) Z COS(yang) / 256
mov [Zt],bp
mov ax,[Xt] ; Switch values
mov [X],ax
mov ax,[Zt]
mov [Z],ax
; Rotate around z-axis
; XT = X COS(zang) - Y SIN(zang) / 256
; YT = X SIN(zang) Y COS(zang) / 256
; X = XT
; Y = YT
mov ax,[X]
mov bx,[ZCos]
imul bx ; ax = X Cos(zang)
mov bp,ax
mov ax,[Y]
mov bx,[ZSin]
imul bx ; ax = Y Sin(zang)
sub bp,ax ; bp = X Cos(zang) - Y Sin(zang)
sar bp,8 ; bp = X Cos(zang) - Y Sin(zang) / 256
mov [Xt],bp
mov ax,[X]
mov bx,[ZSin]
imul bx ; ax = X Sin(zang)
mov bp,ax
mov ax,[Y]
mov bx,[ZCos]
imul bx ; ax = Y Cos(zang)
add bp,ax ; bp = X SIN(zang) Y COS(zang) [Page]
sar bp,8 ; bp = X SIN(zang) Y COS(zang) / 256
mov [Yt],bp
mov ax,[Xt] ; Switch values
mov [X],ax
mov ax,[Yt]
mov [Y],ax
ret
ENDP RotatePoint
PROC ShowPoint
; Calculates screenposition and
; plots the point on the screen
mov ax,[Xoff] ; XoffX / Z Zoff = screen x
mov bx,[X]
imul bx
mov bx,[Z]
add bx,[Zoff] ; Distance
idiv bx
add ax,[Mx] ; Center on screen
mov bp,ax
mov ax,[Yoff] ; YoffY / Z Zoff = screen y
mov bx,[Y]
imul bx
mov bx,[Z]
add bx,[Zoff] ; Distance
idiv bx
add ax,[My] ; Center on screen
mov bx,320
imul bx
add ax,bp ; ax = (y320) x
mov di,ax
mov ax,[Z] ; Get color from Z
add ax,100d ; (This piece of code could be improved)
mov [byte ptr es:di],al ; Place a dot with color al
mov [Erase si],di ; Save position for erase
ret
ENDP ShowPoint
PROC MainProgram
call UpdateAngles ; Calculate new angles
call SetRotation ; Find sine & cosine of those angles
xor si,si ; First 3d-point
mov cx,MaxPoints
ShowLoop:
call RotatePoint ; Rotates the point using above formulas
call ShowPoint ; Shows the point
add si,3 ; Next 3d-point
loop ShowLoop
call WaitVrt ; Wait for retrace
xor si,si ; Starting with point 0
xor al,al ; Color = 0 = black
mov cx,MaxPoints
Deletion:
mov di,[Erase si] ; di = vgapos old point
mov [byte ptr es:di],al ; Delete it
add si,3 ; Next point
loop Deletion
ret
ENDP MainProgram
; === DATA ===
Credits DB 13,10,\"Code by Vulture / Outlaw Triad\",13,10,\"$\"
Label SinCos Word ; 256 values
dw 0,6,13,19,25,31,38,44,50,56 [Page]
dw 62,68,74,80,86,92,98,104,109,115
dw 121,126,132,137,142,147,152,157,162,167
dw 172,177,181,185,190,194,198,202,206,209
dw 213,216,220,223,226,229,231,234,237,239
dw 241,243,245,247,248,250,251,252,253,254
dw 255,255,256,256,256,256,256,255,255,254
dw 253,252,251,250,248,247,245,243,241,239
dw 237,234,231,229,226,223,220,216,213,209
dw 206,202,198,194,190,185,181,177,172,167
dw 162,157,152,147,142,137,132,126,121,115
dw 109,104,98,92,86,80,74,68,62,56
dw 50,44,38,31,25,19,13,6,0,-6
dw -13,-19,-25,-31,-38,-44,-50,-
56,-62,-68
dw -74,-80,-86,-92,-98,-104,-109,-115,-121,-126
dw -132,-137,-142,-147,-152,-157,-162,-167,-172,-177
dw -181,-185,-190,-194,-198,-202,-206,-209,-213,-216
dw -220,-223,-226,-229,-231,-234,-237,-239,-241,-243
dw -245,-247,-248,-250,-251,-252,-253,-254,-255,-255
dw -256,-256,-256,-256,-256,-255,-255,-254,-253,-252
dw -251,-250,-248,-247,-245,-243,-241,-239,-237,-234
dw -231,-229,-226,-223,-220,-216,-213,-209,-206,-202
dw -198,-194,-190,-185,-181,-177,-172,-167,-162,-157
dw -152,-147,-142,-137,-132,-126,-121,-115,-109,-104
dw -98,-92,-86,-80,-74,-68,-62,-56,-50,-44
dw -38,-31,-25,-19,-13,-6
Label Cube Byte ; The 3d points
c = -35 ; 5x5y5z (=125) points
rept 5
b = -35
rept 5
a = -35
rept 5
db a,b,c
a = a 20
endm
b = b 20
endm
c = c 20
endm
Label Palette Byte ; The palette to use
db 0,0,0 ; 633 gray-tint
d = 63
rept 63
db d,d,d
db d,d,d
db d,d,d
d = d - 1
endm
X DW ; X variable for formula
Y DW
Z DW
Xt DW ; Temporary variable for x
Yt DW
Zt DW
XAngle DW 0 ; Angle to rotate around x
YAngle DW 0
ZAngle DW 0
DeltaX DW ; Amound Xangle is increased each time
DeltaY DW
DeltaZ DW
Xoff DW
Yoff DW
Zoff DW ; Distance from viewer
XSin DW ; Sine and cosine of angle to rotate around
XCos DW
YSin DW [Page]
YCos DW
ZSin DW
ZCos DW
Mx DW 160 ; Middle of the screen
My DW 100
MaxPoints EQU 125 ; Number of 3d Points
Erase DW MaxPoints DUP () ; Array for deletion screenpoints
ENDS CODE ; End of codesegment
END START ; The definite end :)
; You may use this code in your own productions but
; give credit where credit is due Only lamers steal
; code so try to create your own 3d-engine and use
; this code as an example
; Thanx must go to Arno Brouwer and Ash for releasing
; example sources
;
; Ciao dudoz,
;
; Vulture / Outlaw Triad
欢迎分享，转载请注明来源：品搜搜测评网
微信扫一扫
支付宝扫一扫
