基于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个

mathematica

可以给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

无心人

过滤出强伪素数的代码

1. #include <stdio.h>
   
2. #include <stdlib.h>
   
3. 
   
4. 
   
5. int main( void )
   
6. {
   
7.   FILE * psp2, *S64;
   
8.   char buf[256];
   
9.   int i=0;
   
10.   if ((psp2=fopen("psp2.txt", "rt"))!=NULL)
    
11.     if ((S64 = fopen("S64.txt", "wt"))!=NULL)
    
12.     {
    
13.       while (!feof(psp2))
    
14.       {
    
15.         fgets(buf, 256, psp2);
    
16.         if (buf[1] == 'S')
    
17.         {
    
18.           fputs(buf, S64);
    
19.           i ++;
    
20.           if (i >= 1000)
    
21.           {
    
22.             i=0;
    
23.             fflush(S64);
    
24.           }
    
25.         }
    
26.       }
    
27.     }
    
28.   fclose(psp2);
    
29.   fclose(S64);
    
30. }
    

复制代码

无心人

基于linux 64位汇编的测试代码
m64.s

1. .data
   
2. .text
   
3. .global powMod, oddDivisorOf, strongPseudoPrimeTest
   
4. 
   
5. powMod:
   
6. #return rax = rdi ^ rsi mod rdx
   
7. #unsigned long long powMod(unsigned long long n, unsigned long long p, unsigned long long m)
   
8.         movq %rdi, %rax
   
9.         movq %rdx, %rcx
   
10.         movq $1, %r8
    
11.         movq $1, %r9
    
12. powModStart:
    
13.         cmpq %r9, %rsi
    
14.         jbe  powModExit
    
15.         testq %r9, %rsi
    
16.         je powModSquare
    
17.         movq %rax, %rdi
    
18.         mulq %r8
    
19.         divq %rcx
    
20.         movq %rdi, %rax
    
21.         movq %rdx, %r8
    
22. powModSquare:
    
23.         mulq %rax
    
24.         divq %rcx
    
25.         movq %rdx, %rax
    
26.         shrq $1, %rsi
    
27.         jmp powModStart
    
28. powModExit:
    
29.         mulq %r8
    
30.         divq %rcx
    
31.         movq %rdx, %rax
    
32.         ret
    
33. 
    
34. oddDivisorOf:
    
35. #rdi = 2^n * rax, rax is odd number
    
36.         bsf %rdi, %rcx
    
37.         movq %rdi, %rax
    
38.         shrq %cl, %rax
    
39.         testq $1, %rdi
    
40.         cmovne %rdi, %rax
    
41.         ret
    
42. 
    
43. strongPseudoPrimeTest:
    
44. #retrun rdi base rsi Strong PseudoPrime Test
    
45.         movl $0, %eax
    
46.         testq $1, %rdi
    
47.         je strongPseudoPrimeTestExit
    
48.         movq %rdi, %xmm0
    
49.         subq $1, %rdi
    
50.         bsf %rdi, %rcx
    
51.         movd %ecx, %xmm1
    
52.         shrq %cl, %rdi
    
53.         xchg %rsi, %rdi
    
54.         movq %xmm0, %rdx
    
55.         call powMod
    
56.         cmpq $1, %rax
    
57.         je strongPseudoPrimeTestExit
    
58.         movd %xmm1, %ecx
    
59.         movq %xmm0, %r8
    
60.         movq %xmm0, %r9
    
61.         subq $1, %r8
    
62.         cmpq %r8, %rax
    
63.         je strongPseudoPrimeTestCheck
    
64.         subl $1, %ecx
    
65.         cmpl $0, %ecx
    
66.         je strongPseudoPrimeTestCheck
    
67. strongPseudoPrimeTestLoop:
    
68.         mulq %rax
    
69.         divq %r9
    
70.         movq %rdx, %rax
    
71.         cmpq %r8, %rax
    
72.         je strongPseudoPrimeTestCheck
    
73.         subl $1, %ecx
    
74.         jne strongPseudoPrimeTestLoop
    
75. strongPseudoPrimeTestCheck:
    
76.         movq $1, %rcx
    
77.         movq $0, %rdx
    
78.         cmpq %r8, %rax
    
79.         cmove %rcx, %rax
    
80.         cmovne %rdx, %rax
    
81. strongPseudoPrimeTestExit:
    
82.         ret

复制代码

S2.c

1. #include <stdio.h>
   
2. #include <stdlib.h>
   
3. 
   
4. extern int strongPseudoPrimeTest(unsigned long long n, unsigned long long base);
   
5. 
   
6. int main( void )
   
7. {
   
8.   FILE *S64, *S23, *S25, *S27, *S235, *S237, *S257, *S2357;
   
9.   char buf[256], C;
   
10.   int C64=0, C23=0, C25=0, C27=0, C235=0, C237=0, C257=0, C2357=0;
    
11.   int P3, P5, P7;
    
12.   unsigned long long p;
    
13.   if ((S64=fopen("S64.txt", "rt"))!=NULL)
    
14.   {
    
15.     if ((S23 = fopen("S23.txt", "wt"))!=NULL)
    
16.     if ((S25 = fopen("S25.txt", "wt"))!=NULL)
    
17.     if ((S27 = fopen("S27.txt", "wt"))!=NULL)
    
18.     if ((S235 = fopen("S235.txt", "wt"))!=NULL)
    
19.     if ((S237 = fopen("S237.txt", "wt"))!=NULL)
    
20.     if ((S257 = fopen("S257.txt", "wt"))!=NULL)
    
21.     if ((S2357 = fopen("S2357.txt", "wt"))!=NULL)
    
22.     {
    
23.       while (!feof(S64))
    
24.       {
    
25.         fgets(buf, 256, S64);
    
26.         C64++;
    
27.         buf[1]=' ';
    
28.         sscanf(buf+2, "%llu", &p);
    
29.         P3=strongPseudoPrimeTest(p, 3);
    
30.         P5=strongPseudoPrimeTest(p, 5);
    
31.         P7=strongPseudoPrimeTest(p, 7);
    
32.         if (P3)
    
33.         {
    
34.                 C23++;
    
35.                 fputs(buf, S23);
    
36. 
    
37.                 if (P5)
    
38.                 {
    
39.                         C235++;
    
40.                         fputs(buf, S235);
    
41. 
    
42.                         if (P7)
    
43.                         {
    
44.                                 C2357++;
    
45.                                 fputs(buf, S2357);
    
46.                         }
    
47.                 }
    
48. 
    
49.                 if (P7)
    
50.                 {
    
51.                         C237++;
    
52.                         fputs(buf, S237);
    
53.                 }
    
54.         }
    
55. 
    
56.         if (P5)
    
57.         {
    
58.                 C25++;
    
59.                 fputs(buf, S25);
    
60. 
    
61.                 if (P7)
    
62.                 {
    
63.                         C257++;
    
64.                         fputs(buf, S257);
    
65.                 }
    
66.         }
    
67. 
    
68.         if (P7)
    
69.         {
    
70.                 C27++;
    
71.                 fputs(buf, S27);
    
72.         }
    
73. 
    
74.       }
    
75.     }
    
76.   }
    
77.   fclose(S64);
    
78.   fclose(S23);
    
79.   fclose(S25);
    
80.   fclose(S27);
    
81.   fclose(S235);
    
82.   fclose(S237);
    
83.   fclose(S257);
    
84.   fclose(S2357);
    
85.   printf("Total: %llu\n", C64);
    
86.   printf("   base 2, 3 -- %llu\n", C23);
    
87.   printf("       base 2, 3, 5 -- %llu\n", C235);
    
88.   printf("           base 2, 3, 5, 7 -- %llu\n", C2357);
    
89.   printf("       base 2, 3, 7 -- %llu\n", C237);
    
90.   printf("   base 2, 5 -- %llu\n", C25);
    
91.   printf("       base 2, 5, 7 -- %llu\n", C257);
    
92.   printf("   base 2, 7 -- %llu\n", C27);
    
93.   return 0;
    
94. }
    

复制代码

linux下编译命令
gcc -O3 -o S2 S2.c m64.s

粗略的算了下，每次测试大概7100时钟左右，不知道是否存在优化的可能

mathematica

无心人 发表于 2016-1-12 16:53
对这些数据进行了强伪素数检测，首先过滤出31894016组base 2的强伪素数，
再进行base 3, 5, 7的单独或者组 ...

你真棒,那么大的数据都能搞得动!