关于伪素数与素性测试的一些事情

无心人

继续刷楼

mathematica

无心人 发表于 2019-9-22 09:43
基于Jacobi和的分圆域素性检验算法

我觉得你可真够无聊的,其实
baillie psw算法已经够好的了,
https://bbs.emath.ac.cn/forum.ph ... 776&fromuid=865
这儿有现成的代码!

无心人

mathematica 发表于 2019-9-29 09:19
我觉得你可真够无聊的,其实
baillie psw算法已经够好的了,
https://bbs.emath.ac.cn/forum.php?mod=vie ...

你可能不知道BPSW的算法里，除了Miller-Rabin的另一个算法，是跟我说的Frobenius类似的算法，你觉得，对它的优化，值不值得去讨论？
你实现那个，是用F-L序列实现的，而对F-L改进的算法是我说的基于二次域的\(1+\sqrt{c}\) 的，实现起来简单多了

mathematica

无心人 发表于 2019-9-29 11:21
你可能不知道BPSW的算法里，除了Miller-Rabin的另一个算法，是跟我说的Frobenius类似的算法，你觉得， ...

你说的我不懂,但是我知道bpsw目前来说还没发现反例,真的很优秀!

mathematica

无心人 发表于 2019-9-29 11:21
你可能不知道BPSW的算法里，除了Miller-Rabin的另一个算法，是跟我说的Frobenius类似的算法，你觉得， ...

一个2为底的miller rabin+一个lucas动参数的测试,
两者联合起来都已经很难找到反例了，
如果多几个miller rabin测试+几个lucas动参数，更难找到反例。
两者的因数分解的模的形式不一样，所以壬就考虑联合起来，
但是没想到很难找到反例！

无心人

mathematica 发表于 2019-10-1 09:56
一个2为底的miller rabin+一个lucas动参数的测试,
两者联合起来都已经很难找到反例了，
如果多几个mill ...

那个二次形式的Frobenius测试，就是你说的lucas在二次域的变形，原理是一样的，你可以看我给的wiki链接，但是，有个重要区别是，整数递归形式的测试，是存在很多伪素数例子的，而二次域的形式，10^8内我没找到伪素数，你可以说整数形式已经足够了，没有反例，但是这只是因为两个测试的伪素数没有发现叠加部分而已，而二次域形式的直接来自二次域的Fermat小定理，有论文分析过出现伪素数的可能，条件要远远的比整数上的fermat小定理的伪素数形式严格，甚至只要排除n是完全平方，二次域形式的Frobenius是可以独立做测试的

无心人

素性测试的基于GMP的C源代码看大数分解程序YAFU的源代码
解压后
头文件在./include/mpz_prp_prime.h
源文件在./top/mpz_prp_prime.c
这是从头文件里扒下来的它实现的测试名字

1. /* ******************************************************************
   
2. * mpz_prp: (also called a Fermat pseudoprime)
   
3. * A "pseudoprime" to the base a is a composite number n such that,
   
4. * (a,n)=1 and a^(n-1) = 1 mod n
   
5. * ******************************************************************/
   
6. int mpz_prp(mpz_t n, mpz_t a);
   
7. 
   
8. /* *************************************************************************
   
9. * mpz_euler_prp: (also called a Solovay-Strassen pseudoprime)
   
10. * An "Euler pseudoprime" to the base a is an odd composite number n with,
    
11. * (a,n)=1 such that a^((n-1)/2)=(a/n) mod n [(a/n) is the Jacobi symbol]
    
12. * *************************************************************************/
    
13. int mpz_euler_prp(mpz_t n, mpz_t a);
    
14. 
    
15. /* *********************************************************************************************
    
16. * mpz_sprp: (also called a Miller-Rabin pseudoprime)
    
17. * A "strong pseudoprime" to the base a is an odd composite n = (2^r)*s+1 with s odd such that
    
18. * either a^s == 1 mod n, or a^((2^t)*s) == -1 mod n, for some integer t, with 0 <= t < r.
    
19. * *********************************************************************************************/
    
20. int mpz_sprp(mpz_t n, mpz_t a);
    
21. 
    
22. /* *************************************************************************
    
23. * mpz_fibonacci_prp:
    
24. * A "Fibonacci pseudoprime" with parameters (P,Q), P > 0, Q=+/-1, is a
    
25. * composite n for which V_n == P mod n
    
26. * [V is the Lucas V sequence with parameters P,Q]
    
27. * *************************************************************************/
    
28. int mpz_fibonacci_prp(mpz_t n, long int p, long int q);
    
29. 
    
30. /* *******************************************************************************
    
31. * mpz_lucas_prp:
    
32. * A "Lucas pseudoprime" with parameters (P,Q) is a composite n with D=P^2-4Q,
    
33. * (n,2QD)=1 such that U_(n-(D/n)) == 0 mod n [(D/n) is the Jacobi symbol]
    
34. * *******************************************************************************/
    
35. int mpz_lucas_prp(mpz_t n, long int p, long int q);
    
36. 
    
37. /* *********************************************************************************************
    
38. * mpz_stronglucas_prp:
    
39. * A "strong Lucas pseudoprime" with parameters (P,Q) is a composite n = (2^r)*s+(D/n), where
    
40. * s is odd, D=P^2-4Q, and (n,2QD)=1 such that either U_s == 0 mod n or V_((2^t)*s) == 0 mod n
    
41. * for some t, 0 <= t < r. [(D/n) is the Jacobi symbol]
    
42. * *********************************************************************************************/
    
43. int mpz_stronglucas_prp(mpz_t n, long int p, long int q);
    
44. 
    
45. /* *******************************************************************************************
    
46. * mpz_extrastronglucas_prp:
    
47. * Let U_n = LucasU(p,1), V_n = LucasV(p,1), and D=p^2-4.
    
48. * An "extra strong Lucas pseudoprime" to the base p is a composite n = (2^r)*s+(D/n), where
    
49. * s is odd and (n,2D)=1, such that either U_s == 0 mod n or V_s == +/-2 mod n, or
    
50. * V_((2^t)*s) == 0 mod n for some t with 0 <= t < r-1 [(D/n) is the Jacobi symbol]
    
51. * *******************************************************************************************/
    
52. int mpz_extrastronglucas_prp(mpz_t n, long int p);
    
53. 
    
54. /* ***********************************************************************************************
    
55. * mpz_selfridge_prp:
    
56. * A "Lucas-Selfridge pseudoprime" n is a "Lucas pseudoprime" using Selfridge parameters of:
    
57. * Find the first element D in the sequence {5, -7, 9, -11, 13, ...} such that Jacobi(D,n) = -1
    
58. * Then use P=1 and Q=(1-D)/4 in the Lucas pseudoprime test.
    
59. * Make sure n is not a perfect square, otherwise the search for D will only stop when D=n.
    
60. * ***********************************************************************************************/
    
61. int mpz_selfridge_prp(mpz_t n);
    
62. 
    
63. /* *********************************************************************************************************
    
64. * mpz_strongselfridge_prp:
    
65. * A "strong Lucas-Selfridge pseudoprime" n is a "strong Lucas pseudoprime" using Selfridge parameters of:
    
66. * Find the first element D in the sequence {5, -7, 9, -11, 13, ...} such that Jacobi(D,n) = -1
    
67. * Then use P=1 and Q=(1-D)/4 in the strong Lucase pseudoprime test.
    
68. * Make sure n is not a perfect square, otherwise the search for D will only stop when D=n.
    
69. * **********************************************************************************************************/
    
70. int mpz_strongselfridge_prp(mpz_t n);
    
71. 
    
72. /* **********************************************************************************
    
73. * mpz_bpsw_prp:
    
74. * A "Baillie-Pomerance-Selfridge-Wagstaff pseudoprime" is a composite n such that
    
75. * n is a strong pseudoprime to the base 2 and
    
76. * n is a Lucas pseudoprime using the Selfridge parameters.
    
77. * **********************************************************************************/
    
78. int mpz_bpsw_prp(mpz_t n);
    
79. 
    
80. /* ****************************************************************************************
    
81. * mpz_strongbpsw_prp:
    
82. * A "strong Baillie-Pomerance-Selfridge-Wagstaff pseudoprime" is a composite n such that
    
83. * n is a strong pseudoprime to the base 2 and
    
84. * n is a strong Lucas pseudoprime using the Selfridge parameters.
    
85. * ****************************************************************************************/
    
86. int mpz_strongbpsw_prp(mpz_t n);
    
87. 
    
88. /* **********************************************************************************
    
89. * APR-CL (also known as APRT-CLE) is a prime proving algorithm developed by:
    
90. * L. Adleman, C. Pomerance, R. Rumely, H. Cohen, and A. K. Lenstra
    
91. * APRT-CLE = Adleman-Pomerance-Rumely Test Cohen-Lenstra Extended version
    
92. * You can find all the details of this implementation in the Cohen & Lenstra paper:
    
93. *    H. Cohen and A. K. Lenstra, "Implementation of a new primality test",
    
94. *    Math. Comp., 48 (1987) 103--121
    
95. *
    
96. * This C/GMP version is a conversion of Dario Alpern's Java based APRT-CLE code
    
97. * His code was based on Yuji Kida's UBASIC code
    
98. *
    
99. * Based on APRT-CLE Written by Dario Alejandro Alpern (Buenos Aires - Argentina)
    
100. * From: Last updated September 10th, 2011. See http://www.alpertron.com.ar/ECM.HTM
     
101. *
     
102. * With improvements based on Jason Moxham's APRCL v1.15 code, from 2003/01/01
     
103. * *********************************************************************************/
     
104. 
     
105. int mpz_aprcl(mpz_t N);
     
106. int mpz_aprtcle(mpz_t N, int verbose);
     

复制代码

mathematica

无心人 发表于 2019-10-7 09:46
素性测试的基于GMP的C源代码看大数分解程序YAFU的源代码
解压后
头文件在./include/mpz_prp_prime.h

在我看来，baillie psw算法足够了，不自信的话，就再加几个miller随机的测试，不纠结！

无心人

更新了SQFT测试的源代码，打算用Julia的Hecke包重写所有代码，Hecko包会包装GMP的一些函数

无心人

补充下 BPSW 测试

### （一）、Fibonacci 序列素性测试

定义 Fibonacci 序列，$F_0= 0， F_1 = F_2 = 1，F_{k+2} = F_{k+1} + F_k$，素数 $n > 5$， jacobi 符号 $J(\frac{n}{5})=j$，则 $F_{n-j} \equiv 0 (mod\ \ n)$。

利用上面的结果，假如 $n$ 是待测试整数，

1、若 $5 | n$ 则返回 $n$ 是合数。

2、否则若 $n \equiv \pm 1 (mod\ \ 5)$ 则 $j = 1$，若 $n \equiv \pm 2 (mod\ \ 5)$ 则 $j = -1$。

3、若 $F_{n-j} \equiv 0 (mod\ \ n)$，则返回 $n$ 可能是素数，否则返回 $n$ 是合数。

若合数 $n$ 通过上述测试，称为 Fibonacci 伪素数，简称 fpp，fpp 是无限多的。

julia示例代码

1. 
   
2. using Hecke
   
3. 
   
4. function fibonacciTest(n)
   
5.     r = n % 5
   
6.     r == 0 && return false
   
7.     if r == 1 || r == 4
   
8.         j = 1
   
9.     end
   
10.     if r == 2 || r == 3
    
11.         j = -1
    
12.     end
    
13.     return fibonacci(n-j) % n == 0
    
14. end
    
15. 
    
16. for i in range(ZZ(7), stop=ZZ(9999), step=2)
    
17.     if fibonacciTest(i)
    
18.         if !isprime(i)
    
19.             println("Fibonacci Pseudoprime: $i")
    
20.         end
    
21.     end
    
22.     if isprime(i)
    
23.         if !fibonacciTest(i)
    
24.             println("Fibonacci Primality test error at: $i")
    
25.         end
    
26.     end
    
27. end
    

复制代码

输出
Fibonacci Pseudoprime: 323
Fibonacci Pseudoprime: 377
Fibonacci Pseudoprime: 1891
Fibonacci Pseudoprime: 3827
Fibonacci Pseudoprime: 4181
Fibonacci Pseudoprime: 5777
Fibonacci Pseudoprime: 6601
Fibonacci Pseudoprime: 6721
Fibonacci Pseudoprime: 8149

对 10000 内大于 5 奇数测试表明，存在 9 个奇伪素数。如果不限定 $n$ 是奇数，则相应存在偶伪素数，最小的一个是 8539786。

另外，该代码是示例性质代码，并没考虑优化，序列的计算没有在模 $n$ 环上进行，当 $n$ 很大时候，可能会内存溢出、计算时间超长等。