基于2^64以下的以2为底的伪素数的数据分析
在下面网站找到了所有的低于2^64的以2为底伪素数PSP(2)的表http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html
下载后进行了一些分析，详细结果，后面列出来
下载的是第3项，包含Carmichael数跟强伪素数分类，和因子分解的部分
首先是统计信息，共118968379个, 其中强伪素数31894015个; Carmichael数4279356个, 既是Carmichael数又是强伪素数的32596个
可以给gxqcn用 对这些数据进行了强伪素数检测，首先过滤出31894016组base 2的强伪素数，
再进行base 3, 5, 7的单独或者组合测试，结果如下，数字代表通过测试的强伪素数数量
Total: 31894016
base 2, 3 -- 1501720
base 2, 3, 5 -- 131157
base 2, 3, 5, 7 -- 16826
base 2, 3, 7 -- 171359
base 2, 5 -- 1330740
base 2, 5, 7 -- 126999
base 2, 7 -- 1363175
过滤出强伪素数的代码
#include <stdio.h>
#include <stdlib.h>
int main( void )
{
FILE * psp2, *S64;
char buf;
int i=0;
if ((psp2=fopen("psp2.txt", "rt"))!=NULL)
if ((S64 = fopen("S64.txt", "wt"))!=NULL)
{
while (!feof(psp2))
{
fgets(buf, 256, psp2);
if (buf == 'S')
{
fputs(buf, S64);
i ++;
if (i >= 1000)
{
i=0;
fflush(S64);
}
}
}
}
fclose(psp2);
fclose(S64);
}
基于linux 64位汇编的测试代码
m64.s
.data
.text
.global powMod, oddDivisorOf, strongPseudoPrimeTest
powMod:
#return rax = rdi ^ rsi mod rdx
#unsigned long long powMod(unsigned long long n, unsigned long long p, unsigned long long m)
movq %rdi, %rax
movq %rdx, %rcx
movq $1, %r8
movq $1, %r9
powModStart:
cmpq %r9, %rsi
jbepowModExit
testq %r9, %rsi
je powModSquare
movq %rax, %rdi
mulq %r8
divq %rcx
movq %rdi, %rax
movq %rdx, %r8
powModSquare:
mulq %rax
divq %rcx
movq %rdx, %rax
shrq $1, %rsi
jmp powModStart
powModExit:
mulq %r8
divq %rcx
movq %rdx, %rax
ret
oddDivisorOf:
#rdi = 2^n * rax, rax is odd number
bsf %rdi, %rcx
movq %rdi, %rax
shrq %cl, %rax
testq $1, %rdi
cmovne %rdi, %rax
ret
strongPseudoPrimeTest:
#retrun rdi base rsi Strong PseudoPrime Test
movl $0, %eax
testq $1, %rdi
je strongPseudoPrimeTestExit
movq %rdi, %xmm0
subq $1, %rdi
bsf %rdi, %rcx
movd %ecx, %xmm1
shrq %cl, %rdi
xchg %rsi, %rdi
movq %xmm0, %rdx
call powMod
cmpq $1, %rax
je strongPseudoPrimeTestExit
movd %xmm1, %ecx
movq %xmm0, %r8
movq %xmm0, %r9
subq $1, %r8
cmpq %r8, %rax
je strongPseudoPrimeTestCheck
subl $1, %ecx
cmpl $0, %ecx
je strongPseudoPrimeTestCheck
strongPseudoPrimeTestLoop:
mulq %rax
divq %r9
movq %rdx, %rax
cmpq %r8, %rax
je strongPseudoPrimeTestCheck
subl $1, %ecx
jne strongPseudoPrimeTestLoop
strongPseudoPrimeTestCheck:
movq $1, %rcx
movq $0, %rdx
cmpq %r8, %rax
cmove %rcx, %rax
cmovne %rdx, %rax
strongPseudoPrimeTestExit:
ret
S2.c
#include <stdio.h>
#include <stdlib.h>
extern int strongPseudoPrimeTest(unsigned long long n, unsigned long long base);
int main( void )
{
FILE *S64, *S23, *S25, *S27, *S235, *S237, *S257, *S2357;
char buf, C;
int C64=0, C23=0, C25=0, C27=0, C235=0, C237=0, C257=0, C2357=0;
int P3, P5, P7;
unsigned long long p;
if ((S64=fopen("S64.txt", "rt"))!=NULL)
{
if ((S23 = fopen("S23.txt", "wt"))!=NULL)
if ((S25 = fopen("S25.txt", "wt"))!=NULL)
if ((S27 = fopen("S27.txt", "wt"))!=NULL)
if ((S235 = fopen("S235.txt", "wt"))!=NULL)
if ((S237 = fopen("S237.txt", "wt"))!=NULL)
if ((S257 = fopen("S257.txt", "wt"))!=NULL)
if ((S2357 = fopen("S2357.txt", "wt"))!=NULL)
{
while (!feof(S64))
{
fgets(buf, 256, S64);
C64++;
buf=' ';
sscanf(buf+2, "%llu", &p);
P3=strongPseudoPrimeTest(p, 3);
P5=strongPseudoPrimeTest(p, 5);
P7=strongPseudoPrimeTest(p, 7);
if (P3)
{
C23++;
fputs(buf, S23);
if (P5)
{
C235++;
fputs(buf, S235);
if (P7)
{
C2357++;
fputs(buf, S2357);
}
}
if (P7)
{
C237++;
fputs(buf, S237);
}
}
if (P5)
{
C25++;
fputs(buf, S25);
if (P7)
{
C257++;
fputs(buf, S257);
}
}
if (P7)
{
C27++;
fputs(buf, S27);
}
}
}
}
fclose(S64);
fclose(S23);
fclose(S25);
fclose(S27);
fclose(S235);
fclose(S237);
fclose(S257);
fclose(S2357);
printf("Total: %llu\n", C64);
printf(" base 2, 3 -- %llu\n", C23);
printf(" base 2, 3, 5 -- %llu\n", C235);
printf(" base 2, 3, 5, 7 -- %llu\n", C2357);
printf(" base 2, 3, 7 -- %llu\n", C237);
printf(" base 2, 5 -- %llu\n", C25);
printf(" base 2, 5, 7 -- %llu\n", C257);
printf(" base 2, 7 -- %llu\n", C27);
return 0;
}
linux下编译命令
gcc -O3 -o S2 S2.c m64.s
粗略的算了下，每次测试大概7100时钟左右，不知道是否存在优化的可能
无心人 发表于 2016-1-12 16:53
对这些数据进行了强伪素数检测，首先过滤出31894016组base 2的强伪素数，
再进行base 3, 5, 7的单独或者组 ...
你真棒,那么大的数据都能搞得动!
页:
[1]