2^64内基于米勒罗宾素性测试的素性证明算法

无心人

对n进行基b的米勒罗宾素性测试定义为MR(n, b)
算法步骤
1、首先进行基 2 的 MR(n, 2)，
          如果测试失败，返回 n 是合数，结束
          否则，继续 2
2、对 n 求 s = $|log_2(log_2(n))|$
3、找到最小正整数 $t_0$ 满足 $t_0^s \ge n$
4、对于大于等于 $t_0$ 的连续素数 t 执行 s + 3 次 MR(n, t)，
          如果某次测试失败，返回 n 是合数，结束
          如果全部通过，返回 n 是素数，结束

此时，只有 214021506061375741 失败
如果改成 $t_0$ 开始连续整数，则有2个失败
9450361108305253
1578072131184991861

nyy

不要搞素数判定了，浪费时间！
有现成的BPSW算法，还不够好吗？

nyy

用你的算法试试这个大合数

2887148238050771212671429597130393991977609459279722700926516024197432\
3037991527331163289831446392259419778031109293496555784189494417409338\
0561511397999942154241693397290542371100275104208013496673175515285922\
6962916775325475044445856101949404200039904432116776619949629539250452\
6987193290703735640322737012784538991261203092448414947289768854060249\
76768122077071687938121709811322297802059565867

https://bbs.emath.ac.cn/forum.ph ... 8&fromuid=14149
如何分解这个变态的卡米歇尔强伪素数？

无心人

汇编写的MR(n, b)

1. ;bool mrtest(uint64_t p, uint64_t b)
   
2. ;win64下，前4个参数如果是整数或者指针保存到rcx, rdx, r8, r9
   
3. ;如果是浮点数，保存到xmm0, xmm1, xmm2, xmm3
   
4. ;如果是混合参数，则按照对应位置选择寄存器
   
5. ;多于4个的参数从右到左压栈，同时栈顶为前4个参数保留32字节影子空间
   
6. ;函数内可以自由使用rax, rcx, rdx, r8, r9, r10, r11, xmm0-xmm5
   
7. ;其他寄存器需要使用时候，则要保证在退出函数时恢复值
   
8. ;返回值保存在rax内，如果是浮点数保存在xmm0
   
9. ;linux, FreeBSD, MacOS等64位系统则遵循另外一个规则
   
10. ;前6个参数保存到rdi, rsi, rdx, rcx, r8, r9,
    
11. ;浮点数则保存到xmm0-xmm5
    
12. ;多于6个从右到左依次压栈，没有影子空间，
    
13. ;函数内可以自由使用rax, rcx, rdx, rsi, rdi, r8, r9, r10, r11
    
14. ;返回值保存到rax或者xmm0
    
15. USE64
    
16. section .bss
    
17. section .data
    
18. section .text
    
19. 
    
20. GLOBAL mrtest
    
21. 
    
22. mrtest:
    
23. %ifidn __OUTPUT_FORMAT__, win64
    
24.   mov r8, rcx
    
25.   mov r9, rdx
    
26. %else
    
27.   mov r8, rdi
    
28.   mov r9, rsi
    
29. %endif
    
30. ;r8 = p, r9 = b
    
31.   mov r10, r8
    
32.   mov rax, 1
    
33.   movq xmm0, rax
    
34.   movq xmm3, r8
    
35. ;xmm3 = p-1
    
36.   psubq xmm3, xmm0
    
37.   pxor xmm2, xmm2
    
38. ;xmm2 = 1
    
39.   paddq xmm2, xmm0
    
40.   sub r10, 1
    
41. ;r10 = p - 1
    
42.   bsf rcx, r10
    
43. ;rcx = k
    
44.   shr r10, cl
    
45. ;r10 = m, p - 1 = m * 2^k
    
46.   mov r11, 1
    
47.   cmp r10, 1
    
48.   je jmp1
    
49. loop0:
    
50.   cmp r10, 0
    
51.   je jmp0
    
52.   test r10, 1
    
53.   jz  loop1
    
54. ;r11 = r11 * r9 % r8
    
55.   mov rax, r11
    
56.   mul r9
    
57.   div r8
    
58.   mov r11, rdx
    
59. loop1:
    
60. ;r9 = r9 * r9 % r8
    
61.   mov rax, r9
    
62.   mul r9
    
63.   div r8
    
64.   mov r9, rdx
    
65.   shr r10, 1
    
66.   jmp loop0
    
67. jmp0:
    
68. ; xmm0 = (r11 == 1 || r11 == p-1)? -11 : 0
    
69.   movq xmm0, r11
    
70.   movq xmm1, r11
    
71.   pcmpeqq xmm0, xmm2;xmm0=r11 == 1?(-1):0
    
72.   pcmpeqq xmm1, xmm3;xmm1=r11 == (p-1)?(-1):0
    
73.   por xmm0, xmm1
    
74. ;if (rax = -xmm0 == 1) then exit
    
75.   movq rax, xmm0
    
76.   neg rax
    
77.   cmp rax, 1
    
78.   je exit0
    
79.   mov r9, r11
    
80. jmp1:
    
81.   sub rcx, 1
    
82.   jz exit0
    
83. loop3:
    
84.   mov rax, r9
    
85.   mul r9
    
86.   div r8
    
87.   mov r9, rdx
    
88. ;xmm1 = r9 == p-1? (-1):0
    
89.   movq xmm1, r9
    
90.   pcmpeqq xmm1, xmm3
    
91.   movq rax, xmm1
    
92.   neg rax
    
93.   cmp rax, 1
    
94.   je exit0
    
95.   sub rcx, 1
    
96.   jnz loop3
    
97. exit0:   
    
98.   ret
    

复制代码

测试代码

1. extern bool mrtest(uint64_t p, uint64_t b);
   
2. 
   
3. uint64_t ts[8][2] = {{257,2},  {255,2}, {2047, 2},  {121, 3},
   
4.      {10095020520187, 2}, {581130733, 3}, {55031295157, 234587},
   
5.      {55031295179, 234587}};
   
6. 
   
7. void testmrtest( void )
   
8. {
   
9.   for (uint64_t i = 0; i < 8; i ++)
   
10.     if (mrtest(ts[i][0], ts[i][1]))
    
11.       printf("mrtest(%lu, %lu) is true\n", ts[i][0], ts[i][1]);
    
12.     else
    
13.       printf("mrtest(%lu, %lu) is false\n", ts[i][0], ts[i][1]);
    
14. }
    
15. 
    
16. int main( void )
    
17. {  
    
18.     testmrtest();
    
19. }
    

复制代码

假设汇编名为 mr.asm, C代码 mrtest.c
linux下
nasm -f elf64 mr.asm -o mr.o
gcc mrtest.c mr.o -o mrtest
windows下
nasm -f win64 mr.asm
gcc mrtest.c mr.obj -o mrtest
即可生成测试程序

汇编代码也可以用yasm编译通过

nyy

无心人 发表于 2023-4-25 10:01
汇编写的MR(n, b)

其实我主要是为了搞明白64汇编怎么和C语言混合编程

你还不如直接从机器语言开始编程，混合编程我不会，我连汇编都不会，我只会简单弱智的mathematica

nyy

无心人 发表于 2023-4-25 10:01
汇编写的MR(n, b)

假设
n=2887148238050771212671429597130393991977609459279722700926516024197432\
3037991527331163289831446392259419778031109293496555784189494417409338\
0561511397999942154241693397290542371100275104208013496673175515285922\
6962916775325475044445856101949404200039904432116776619949629539250452\
6987193290703735640322737012784538991261203092448414947289768854060249\
76768122077071687938121709811322297802059565867

那么Log[2, Log[2, n]]，对2取两次对数，结果是10.3631，
但是这个n，能通过1到307的强伪素数的测试！

nyy

无心人 发表于 2023-4-25 10:01
汇编写的MR(n, b)

那你二次域给我看看呀

nyy

nyy 发表于 2023-4-25 13:31
假设
n=2887148238050771212671429597130393991977609459279722700926516024197432\
30379915273311632 ...

https://bbs.emath.ac.cn/forum.ph ... 3&fromuid=14149

用lucas U lucas V就行了！

无心人

各一次基2，3的MR
然后卢卡斯测试或者二次域测试，我有时间比较下哪个快

nyy

无心人 发表于 2023-4-25 16:08
各一次基2，3的MR
然后卢卡斯测试或者二次域测试，我有时间比较下哪个快

整天研究素数判定，难道你还想搞出个大新闻？？？？？？？