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

无心人 · 发表于 2019-9-22 08:49:12

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？欢迎注册

×

如果奇合数$ q $满足费马素性测试 $ b^{(q-1)/2} \equiv \pm1  \mod q$ 称为以b为基的Fermat伪素数（Fermat pseudoprime）,以下简称fpsp(b)

更进一步，如果奇合数$ q $，$ q - 1 = 2^k m, m$为奇数，且
满足下面两个条件之一，（米勒罗宾测试）
1、$ b^{m} \equiv \pm1  \mod q$
2、$ b^{2^im} \equiv \pm1  \mod q, i < k$
称为以b为基的强伪素数（strong pseudoprime）,以下简称spsp(b)

无心人 · 发表于 2019-9-22 09:02:35

有用的几个性质
1、若n是fpsp(b)，则gcd(n, b)=1
2、spsp(b)必然是fpsp(b)
3、如果n同时是是fpsp(b1), fpsp(b2)那么n必然是fpsp(b1b2) ,反之则不成立，比如最小的fpsp(6)是35，而35既不是fpsp(2)，也不是fpsp(3)
4、存在一类整数C，是所有满足gcd(C, x)=1的fpsp(x)，即卡米切尔数（Carmichael number ）
5、n以内卡米切尔数大于$ n^{2/7 }$
6、spsp(b)要比fpsp(b)少的多，很多卡米切尔数是fpsp(b)，但不是spsp(b)
7、简单的计算可以发现一个事实，n以内的卡米切尔数，spsp(b), fpsp(b)，数量依次递增

有人证明n以内卡米切尔数个数$ C(n) > n^{0.332} $，计算表明10^15的卡米切尔数有105212个，所以，$ C(n) > \sqrt[3]{n} , n >= 10^{15}$ ?
我猜想，n以内spsp(b)数量大于 $ \sqrt[3]{n} $，而小于 $ \sqrt{n} $

无心人 · 发表于 2019-9-22 09:29:31

GMP里的素性测试用了高度复合数210，原理未知，可能是一次性排除待测试数字是2，3，5，7的倍数的可能吧

理论计算表明如果n是奇数，通过一次米勒罗宾测试是合数的概率是1/4，费马素性测试合数的概率是1/2
2^64内的实际计算表明，真实概率远小于这个，总体上，通过费马素性测试的合数多于米勒罗宾测试，
所以实践中采用多次米勒罗宾测试来减少出错可能

实践计算表明，spsp(2)数量要少于spsp(b)， b > 2，一般一次测试采用b=2
同样，因为合数b并不能减少米勒罗宾测试出错概率，所以我倾向于采用素数的基b的米勒罗宾测试，即用b=2,3,5,7....

用ψn表示通过前n个素数的强伪素数测试的最小合数，则：
ψ1 =  2047 = 23 * 89
ψ2 =  1373653 = 829 * 1657
ψ3 =  25326001 = 2251 * 11251
ψ4 =  32150 31751 = 151 * 751 * 28351
ψ5 =  215 23028 98747 = 6763 * 10627 * 29947
ψ6 =  347 47496 60383 = 1303 * 16927 * 157543
ψ7 = ψ8 =  34155 00717 28321 = 10670053 * 32010157
ψ9 = ψ10 = ψ11 =  3825 12305 65464 13051 = 149491 * 747451 * 34233211

实践中，也可以采取随机整数基测试

如果n是大整数，二进制有k位，一次米勒罗宾测试出错的实际概率 $ < k^2*4^{2−\sqrt{k}} $ （Damg˚ard et al. 1993）
实际上k = 500, 概率小于 $ 4^{-28} $，当然我猜想，应该小于$ 4^{-125} $，即2-4次就足够坚信n是素数了（张振祥建议次数用6次）

最近对10^28以上10^12内整数的实际测试，未找到任何基2伪素数和强伪素数，后来按照2^64内伪素数出现概率分析，
10^28范围，强伪素数应该每10^18左右的数字才会有1个，远小于上面一段的概率公式的结果

无心人 · 发表于 2019-9-22 09:42:52

Julia实现Fermat素性测试源代码

function fermatProbablePrimeTest(n, b)
if gcd(n, b) > 1
return false
end
a = powermod(b, (n-1)>>1, n)
return (a == 1) || (a == (n-1))
end
println("341 test: ", fermatProbablePrimeTest(341, 2))
println("561 test: ", fermatProbablePrimeTest(561, 2))
println("2047 test: ", fermatProbablePrimeTest(2047, 2))
println("I17 test: ", fermatProbablePrimeTest(11111111111111111, 2))
println("I19 test: ", fermatProbablePrimeTest(1111111111111111111, 2))
println("2^31-1 test: ", fermatProbablePrimeTest(2^31-1, 2))
println("2^67-1 test: ", fermatProbablePrimeTest(2^67-1, 2))

复制代码

341 test: true
561 test: true
2047 test: true
I17 test: false
I19 test: true
2^31-1 test: true
2^67-1 test: false

无心人 · 发表于 2019-9-22 09:42:54

Julia实现Miller Rabin强素性测试源代码

function millerRabinProbablePrimeTest( n, b )
#n - 1 = 2^k * m, m % 2 != 0
k = trailing_zeros(n-1)
m = (n-1) >> k
a = powermod(b, m, n)
if (a == 1) || (a == (n-1))
return true
else
for i in 1:k-1
a = powermod(a, 2, n)
if a == n - 1
return true
end
end
end
return false
end
println("====================")
println("Miller Rabin Probable Prime Test")
println("341 test: ", millerRabinProbablePrimeTest(341, 2))
println("561 test: ", millerRabinProbablePrimeTest(561, 2))
println("2047 test: ", millerRabinProbablePrimeTest(2047, 2))
println("I17 test: ", millerRabinProbablePrimeTest(11111111111111111, 2))
println("I19 test: ", millerRabinProbablePrimeTest(1111111111111111111, 2))
println("2^31-1 test: ", millerRabinProbablePrimeTest(2^31-1, 2))
println("2^67-1 test: ", millerRabinProbablePrimeTest(2^67-1, 2))
println("====================")

复制代码

====================
Miller Rabin Probable Prime Test
341 test: false
561 test: false
2047 test: true
I17 test: false
I19 test: true
2^31-1 test: true
2^67-1 test: false
====================

julia> millerRabinProbablePrimeTest(big(2)^607-1, 2)
true

julia> millerRabinProbablePrimeTest(big(2)^9941-1, 2)
true

julia> millerRabinProbablePrimeTest(big(2)^9947-1, 2)
false

@time millerRabinProbablePrimeTest(big(2)^21701-1, 2)
1.604359 seconds (1.97 k allocations: 809.485 KiB)
true

@time millerRabinProbablePrimeTest(big(2)^132049-1, 2)
128.894809 seconds (1.60 M allocations: 24.091 GiB, 0.56% gc time)
true

无心人 · 发表于 2019-9-22 09:42:58

二次域Frobenius素性测试，

有用广义斐波那契卢卡斯序列描述的，但是我更倾向于用二次域描述的算法，更简洁
/*
备注下，BPSW用的是递归序列形式的测试，参考
1、https://en.wikipedia.org/wiki/Frobenius_pseudoprime
2、https://en.wikipedia.org/wiki/Baillie%E2%80%93PSW_primality_test
是有很多伪素数的，
而这里的二次域形式的根据网页1提供的线索，伪素数出现可能极低，至少wiki并没有列出来
如果有，在n小于2^60情况下，至少有个低于20的素因子或者c>128（c很大的可能很低很低）
*/

n是待测试正整数，不是完全平方数
令$C=[-1,2,3,5,7....]$是-1和素数的集合，$c \in C$
定义二次域，$ \ZZ(\sqrt{c}) $ 上面的整数 $ z = a + b \sqrt{c}  $, 其中 $a, b \in \ZZ $，且Jacobi符号$ (\frac{c}{n}) = -1 $
定义$\overline{z}= a - b \sqrt{c} $

第一种形式是
令$c$是最小满足$ c \in C $, jacobi符号$ (\frac{c}{n}) = -1 $的数
$z=1 + \sqrt{c}$
$ z^n \equiv \overline{z}  \mod n $

第二种形式是
令$c$是满足$ c \in C $, jacobi符号$ (\frac{c}{n}) = -1 $的数，
$ z = a + b \sqrt{c},  gcd(a, n)=1, gcd(b,n)=1$
$ z^n \equiv \overline{z}  \mod n $

第一种形式，是第二种形式，a, b都等于1的特例
如果n是素数，则两种形式都成立，
如果n是合数，能使第一种形式恒等式成立的是Frobenius伪素数，简称FPP
如果n是合数，能使第二种形式恒等式成立的是以a,b,c为参数的Frobenius伪素数，简称FPP(a,b,c)

第一种形式，10^8内没发现伪素数，经过了双向验证，
Frobenius验证通过的，测试是不是非素数，没有发现，
常规素性测试通过的，测试是不是Frobenius验证为非素数，没有发现

2019-10-7，暴力遍历了所有$2^{64}$以内的基2强伪素数，除了$ 1093^2, 3511^2$这两个会导致测试死循环外，其他都通过测试，确认非素数
2019-10-9，测试代码里更新了，如果满足Jacobi Symbol (p, n) = -1的最小奇素数p大于137，就检测n是不是完全平方数，如果是，就直接返回非素数，
测试表明，大多数非完全平方数，满足Jacobi Symbol (p, n) = -1的最小奇素数p小于60，大于60的比例很低

Julia源代码用Hecke库重写了--2021-12-2

using Hecke
function qfnegmod(F, x, n)
re = mod(ZZ(-coeff(x, 0)), n)
ir = mod(ZZ(-coeff(x, 1)), n)
return F([re, ir])
end
function qfconjmod(F, x, n)
re = coeff(x, 0)
ir = mod(ZZ(-coeff(x, 1)), n)
return F([re, ir])
end
function qfmulmod(F, x, y, n)
z = x * y
re = mod(ZZ(coeff(z, 0)), n)
ir = mod(ZZ(coeff(z, 1)), n)
return return F([re, ir])
end
function qfpowermod(F, x, p , n)
@assert p >= 0
p == 0 && return F([1, 0])
b = x
t = ZZ(prevpow(BigInt(2), BigInt(p)))
r = F([1, 0])
while true
if p >= t
r = qfmulmod(F, r, b, n)
p -= t
end
t >>= 1
t <= 0 && break
r = qfmulmod(F, r, r ,n)
end
return r
end
function frobenius(n)
small = [
-1, 2, 3, 5, 7, 11, 13, 17, 19, 23,
29, 31, 37, 41, 43, 47, 53, 59,
61, 67, 71, 73, 79, 83, 89, 97,
101,103,107,109,113,127
]
find = false
global p = 0
for c in small
if jacobi_symbol(ZZ(c), ZZ(n)) == -1
find = true
global p = c
break
end
end
if !find
#n is prefect square number?
sq = isqrt(n)
if sq * sq == n
return false
end
c = next_prime(last(small)+1)
while jacobi_symbol(ZZ(c), ZZ(n)) != -1
c = next_prime(c+1)
end
global p = c
end
#println("$n: $p")
F, d = quadratic_field(p)
b = F([1, 1])
r = qfpowermod(F, b, n, n)
nb = qfconjmod(F, b, n)
return r == nb
end
global b = big(10)^67
@time for i in range(3, step = 2, stop = 99999)
test = b + i
if isqrt(test)^2 == test
continue
end
if isprime(test)
if !frobenius(test)
println("frobenius error at $test")
end
end
if frobenius(test)
if !isprime(test)
println("frobenius error at $test")
end
end
end
println("end")

复制代码

详细的理论分析见附件论文

无心人 · 发表于 2019-9-22 09:43:03

Extended Quadratic Frobenius Primality Test
扩展的二次域Frobenius素性检测算法，简称EQFT

定义环$\RR(n, c) = \ZZ(x)/(n, x^2-c)$
环上的整数$z=ax+b, a, b \in \ZZ_n, \overline{z}=-ax+b, N(z)=\overline{z} \bullet z=b^2-ca^2$

算法1，Extended Quadratic Frobenius Test (EQFTac)
前置测试：
输入，奇数$n\geq5$
1: 如果n被不大于13素数整除，返回n是合数
2: 如果n是完全平方数，返回n是合数
3: 找到最小的c，jacobi符号$(\frac{c}{n})=-1$，输出c
第二部分测试：
输入$n, c, r_3 \in \{1\} \cup \{\xi \in R(n, c) | \phi_3(\xi)=0\}, r_4 \in \{1, -1\} \cup \{\xi \in R(n, c) | \phi_4(\xi)=0\}$
其中$\phi_3(x)=x^2+x+1, \phi_4(x)=x^2+1 $
$u, v$满足$n^2-1=2^u3^vq, (q, 6)=1$
4、随机选择$z \in R(n,c)^*$满足$(\frac{N(z)}{n})=1$
5、如果$z^n \neq \overline{z}$或者$ z^{(n^2-1)/2} \neq 1$，返回n是合数
6、如果$ z^{3^vq} \neq 1 $并且$ z^{2^i3^vq} \neq -1 $ 对所有$i = 0, ..., u-2$，返回n是合数
7、如果6找到的$i_0 \geq 1$, 满足 $ z^{2^{i_0}3^vq}=-1 $ (最多只能有一个值)
令$R_4(z)=z^{2^{i_0 - 1}3^vq} $否则，令$ R_4(z)=z^{3^vq}(=\pm 1) $
如果$ r_4 \neq \pm 1 $并且$ R_4(z)  \notin \{\pm 1, \pm r_4 \} $ ，返回n是合数
8、如果$ z^{2^uq} \neq 1 $并且$ \phi_3(z^{2^u3^iq}) \neq 0 $ 对所有$i = 0, ..., v-1$，返回n是合数
9、如果8找到的$i_0 \geq 1$, 满足 $ \phi_3(z^{2^u3^{i_0}q}) = 0 $ ,  (最多只能有一个值)
令$R_3(z)=z^{2^u3^{i_0}q} $,否则，令$ R_3(z)= 1 $
如果$ r_3 \neq 1 $并且$ R_3(z)  \notin \{1, r_3^{\pm 1} \} $ ，返回n是合数
10、如果$ r_3 = 1 $ 并且 $ R_3(z) \neq 1 $
   则令 $ s_3 = R_3(z) $ 否则 $ s_3 = r_3 $
   如果 $ r_4 = \pm 1 $ 并且 $ R_4(z) \neq \pm 1 $
   则令  $ s_4 = R_4(z) $ 否则 $ s_4 = r_4 $
   返回n可能是素数，和$ s_3, s_4 $

算法2，Extended Quadratic Frobenius Test (EQFTwc)
输入，奇数$ n \geq 5$
输出, n是合数或者可能素数， $ c \in \ZZ_n, r_{24} \in R(n, c)^*, (\frac{c}{n}) = -1, \phi_{24}(r_{24}) = 0 $
1、如果，n能被2，3整除，返回n是合数
2、如果，n是完全平方数，或者完全立方数，返回n是合数
3、找到小c，满足$ (\frac{c}{n}) = -1 $
4、计算, $ r \in R(n, c) $满足$ r^2 + r + 1 = 0 $, 可能返回n是合数
5、a:  如果 $ n \mod 3 \equiv 1 $ 随机选择 $ z \in R(n, c)^* $ 且 $ (\frac{N(z)}{n}) = -1, res_3(z) \neq 1 $
   b:  否则 $ n \mod 3 \equiv 2 $
            重复
                  做 Miller-Rabbin 测试，可能返回n是合数
                  随机选择 $ z \in R(n, c)^* $ 满足 $ (\frac{N(z)}{n}) = -1 $ 计算 $ res_3(z) $
            直到，或者  Miller-Rabbin返回n是合数，或者找到$ z $满足 $ res_3(z) \neq 1 $
6、如果 $ z^n \neq \overline{z} $ 返回n是合数
7、令 $ r_{24} = z^{(n^2-1)/24}$ 如果 $ r_{24}^8 \neq r^{\pm 1} $ 或者 $ r_{24}^{12} \neq - 1 $返回n是合数
8、返回n是可能素数，$c, r_{24} $

子序列迭代
输入 $ n, c, r_{24} $ 满足条件如上
9、随机选择$ z \in R(n, c)^* $
10、如果 $ z^n \neq \overline{z} $ 返回n是合数
11、如果 $ z^{(n^2-1)/24} \notin \{ r_{24}^i | i = 0, ...., 23 \}$返回n是合数
12、返回，n是可能素数

这个翻译的比较崩溃，说实话，有点看不懂加迷惑，好在下面的算法要比这个优秀，完全可以实现出来，这个就不实现了
原版论文已上传，有问题的看论文吧

无心人 · 发表于 2019-9-22 09:43:06

EQFT的简化实现在附件的论文里，
称为$Simplified \: Quadratic \: Frobenius \: test$，简称$SQFT$

下面算法中的$ \phi_3(x)=x^2+x+1, \phi_4(x)=x^2+1, \phi_8(x)=x^4+1$

一、前置算法 $MR2 (Miller-Rabin \: test \: with \: basis \: two \: or \: a \: small \: nonresidue)$
输入：奇数$n$
输出： $n$为合数，或者
         $ c, (\frac{c}{n}) = -1，\epsilon \in R(n, c),  \epsilon^4 = -1,  \phi_8(\epsilon) = 0 $
1、$ n \mod 4 = 3 $
计算 $ \alpha = 2^{(n-3)/4}  (\mod n) $
如果$ 2\alpha^2 \neq \pm 1  (\mod n) $ 输出$n$是合数
否则，输出$c = -1, \epsilon = \alpha + \alpha x $
2、$ n \mod 8 = 5 $
计算 $ \alpha = 2^{(n-1)/4}  (\mod n) $
如果$ \alpha^2 \neq - 1  (\mod n) $ 输出$n$是合数
否则，输出$c = 2, \epsilon = \frac{1 + \alpha}{2} x $

（这里的分数，在实际上的模 n 环 $R_n$ 中的计算可以转换成整数，即如果$\alpha$是偶数，
$\frac{1+\alpha}{2} \equiv \frac{1+\alpha+n}{2} (\mod n)$）
3、$ n \mod 8 = 1 $
如果$n$是完全平方数，输出$n$是合数
否则找到小的$c$满足$ (\frac{c}{n}) = - 1 $ //c从3开始
计算 $ \alpha = c^{(n-1)/8}  (\mod n) $
如果$ \alpha^4 \neq - 1  (\mod n) $ 输出$n$是合数
否则，输出$c = -1, \epsilon = \alpha $

二、SQFT循环 $SQFT_{round}$
输入：奇数$n$，$ c, (\frac{c}{n}) = -1，\epsilon \in R(n, c),  \epsilon^4 = -1,  \phi_8(\epsilon) = 0 $
输出：$n$是合数或者可能是素数
1、随机选择$ z \in R(n, c) $满足$ (\frac{N(z)}{n}) = -1 $
无心人备注：这个满足条件可以去掉，考虑1加上小素数的集合，$ z \in R(n, c) $ 表达成$a + bx $，那么$ a, b \in \{ q | q = 1 \quad or \quad q \in P, q < T\} $ 这里$T$可以取个固定上限比如100，则满足26*25=650次测试需求了
2、如果$ z^n \neq \overline{z} $ 输出$n$是合数
3、如果$ z^{(n^2-1)/8} \notin \{ \pm \epsilon, \pm \epsilon^3 \} $ 输出$n$是合数
无心人备注：经过用29实际测试，发现这里应该是 $ \{ \pm 1, \pm \epsilon, \pm \epsilon^2, \pm \epsilon^3 \} $
4、输出$n$可能是素数

三、SQFT测试 $SQFT: (Simplified \: Quadratic \: Frobenius \: test) $
输入： $n, n > 20 $, 测试次数$t$
输出：$n$是合数或者可能是素数
1、如果$n$被小于$20$素数整除，输出$n$是合数
2、调用算法$MR2$
3、如果$n$没被判定为合数，则$MR2$输出$c, \epsilon $
4、重复$t$次：
用$n, c, \epsilon $调用算法$ SQFT_{round} $，如果算法判定$n$是合数，则终止
5、输出$n$可能是素数

附加三次根测试的$SQFT3$算法
四、SQFT3循环 $ SQFT3_{round} $
输入：整数$ n, gcd(n, 6) = 1 $, 小整数$ c, (\frac{c}{n} ) = -1 $,
      值$ \epsilon, \epsilon^4 = -1 $和$ \epsilon_3 $满足 $ \epsilon_3  =1 $ 或者 $ \phi_3(\epsilon_3) = 0$,

输出：$ n $是合数，或者是可能的素数，同时输出下面的值
      $ {\epsilon'}_3 $, 满足 $ {\epsilon'}_3  =1 $ 或者 $ \phi_3({\epsilon'}_3) = 0$,
      如果， $ \epsilon_3 \neq 1$，那么  ${\epsilon'}_3  = \epsilon_3^{\pm 1} $
1、随机选择$ z \in R(n, c) $满足$ (\frac{N(z)}{n}) = -1 $
无心人备注：这里应该跟上面一样，忽略掉对z的范的Jacobi符号要求
2、如果$ z^n \neq \overline{z} $ 输出$n$是合数
3、如果$ z^{(n^2-1)/8} \notin \{ \pm \epsilon, \pm \epsilon^3 \} $ 输出$n$是合数
无心人备注：这里跟上面的SQFTRound一致即可
4、置$ u = v_3(n^2 - 1), n^2 - 1 = 3^u r $  //这里的$v_3$看不出什么意思来~
5、令$ i = \min\{ j: 0 \leq j \leq u, z^{3^j r} = 1 \} $
6、如果$ i = 0$，输出$n$可能是素数，$ {\epsilon'}_3 = \epsilon_3 $
7、置$ {\epsilon'}_3 = z^{3^{i-1}r}  $
8、如果$ \epsilon_3 = 1 $ 且 $ \phi_3({\epsilon'}_3) \neq 0 $，输出n是合数
9、如果$ \epsilon_3 \neq 1 $ 且$ {\epsilon'}_3 \neq \epsilon_3^{\pm 1} $，输出n是合数
无心人备注：这里猜测和SQFTRound一样，是符合三次根，SQFTRound是八个八次根
10、输出$n$可能是素数和 ${\epsilon'}_3$

五、SQFT3测试 $ SQFT3 $
输入： $n, n > 200 $, 测试次数$t$
输出：$n$是合数或者可能是素数
1、如果 $n$被小于$200$素数整除，输出$n$ 是合数
2、调用$MR2$，如果判定 $n$是合数，结束
3、从$MR2$获得输出的$ c, \epsilon $，并且置 $ \epsilon_3 = 1 $
4、循环$t$次
用$ n, c, \epsilon, \epsilon_3 $调用$ SQFT3_{round} $
如果$ SQFT3_{round} $判定$n$是合数，结束
置$ \epsilon_3 = {\epsilon'}_3 $，这里  $ {\epsilon'}_3 $是算法$ SQFT3_{round} $的输出
5、输出$n$可能是素数

乘法实现：
$w, z \in R(n, c)， w = A_wx + B_w, z = A_zx + B_z $
$z · w = (m_1 + m_2)x + (cA_z + B_z)(A_w + B_w) − cm_1 − m_2 $
$ m_1 = A_zB_w , m_2 = B_zA_w $
这里其实是带入$x=\sqrt{c}$来计算的

Julia源代码，不包含SQFT3实现

using Hecke
function qfnegmod(F, x, n)
re = mod(ZZ(-coeff(x, 0)), n)
ir = mod(ZZ(-coeff(x, 1)), n)
return F([re, ir])
end
function qfconjmod(F, x, n)
re = coeff(x, 0)
ir = mod(ZZ(-coeff(x, 1)), n)
return F([re, ir])
end
function qfmulmod(F, x, y, n)
z = x * y
re = mod(ZZ(coeff(z, 0)), n)
ir = mod(ZZ(coeff(z, 1)), n)
return return F([re, ir])
end
function qfpowermod(F, x, p , n)
@assert p >= 0
p == 0 && return F([1, 0])
b = x
t = ZZ(prevpow(BigInt(2), BigInt(p)))
r = F([1, 0])
while true
if p >= t
r = qfmulmod(F, r, b, n)
p -= t
end
t >>= 1
t <= 0 && break
r = qfmulmod(F, r, r ,n)
end
return r
end
function makeEpsilonList(F, e, n)
list = []
one = F([1, 0])
append!(list, [one])
append!(list, [qfnegmod(F, one, n)])
append!(list, [e])
append!(list, [qfnegmod(F, e, n)])
e2 = qfmulmod(F, e, e, n)
append!(list, [e2])
append!(list, [qfnegmod(F, e2, n)])
e3 = qfmulmod(F, e2, e, n)
append!(list,[e3])
append!(list, [qfnegmod(F, e3, n)])
return list
end
function MR2(n)
if n % 4 == 3
alpha = powermod(ZZ(2), (n-3)>>2, n)
tmp = (2*alpha*alpha) % n
(tmp != 1) && (tmp != (n-1)) && return (false, 0, [])
F, t = quadratic_field(-1)
return (true, F, makeEpsilonList(F, F([alpha, alpha]), n))
end
if n % 8 == 5
alpha = powermod(ZZ(2), (n-1)>>2, n)
tmp = (alpha*alpha) % n
(tmp != (n-1)) && return (false, 0, [])
iseven(alpha) && (alpha+=n)
alpha = (alpha+1)>>1
F, t = quadratic_field(2)
return (true, F, makeEpsilonList(F, F([0, alpha]), n))
end
if n % 8 == 1
issquare(n) && return (false, 0, [])
c = ZZ(3)
while jacobi_symbol(c, n) != -1
c = next_prime(c)
end
alpha = powermod(c, (n-1)>>3, n)
tmp = powermod(alpha, 4, n)
(tmp != (n-1)) && return (false, 0, [])
F, t = quadratic_field(c)
return (true, F, makeEpsilonList(F, F([alpha, 0]), n))
end
end
function SQFTRound(F, x, n, list)
r = qfpowermod(F, x, n, n)
tmp = qfconjmod(F, x, n)
#println("1:$n: x:$x, r:$r, conj(x):$tmp")
r != tmp && return false
tmp = (n^2-1) >> 3
(td, tr) = divrem(tmp, n)
tmp1 = qfpowermod(F, r, td, n)
tmp2 = qfpowermod(F, x, tr, n)
r = qfmulmod(F, tmp1, tmp2, n)
#println("2: $n: $x, $r")
return in(r, list)
end
function SQFT(n, t)
Prm = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]
for p in Prm
(n % p == 0) && return (n == p)
end
Q = append!([1], Prm)
(r, F, list) = MR2(n)
if r
#println("MR2 True at $n: $F")
#println("E: $list")
i = 0
for x in Q
for y in Q
if x != y
if SQFTRound(F, F([x, y]), n, list)
i = i + 1
(i > t) && return true
else
return false
end
end
end
end
else
return false
end
end
########################################################
global b = ZZ(10)^67
@time for i in range(3, step = 2, stop = 19999)
test = b + i
isprime(test) && !SQFT(test, 1) && println("frobenius error at $test")
SQFT(test, 1) && !isprime(test) && println("frobenius error at $test")
end
println("end")

复制代码

对于可能是奇数的平方数的情况
如果n不能被3，5整除
n模240必然是[1, 49, 121, 169]之一

无心人 · 发表于 2019-9-22 09:43:17

基于Jacobi和的分圆域素性检验算法
假设需要证明$n$的素性
1、确定高度复合数$t$，$t$有小因子
比如，$t=5040=2^4\*3^2\*5\*7$
2、确定
$e(t)=2\prod_{q-1|t}q^{v_q(t)+1}$
其中$q$是素数，$v_q(t)$是$q$在$t$中出现的次数，保证$e(t)>\sqrt{n}$
比如$e(5040)=2^6\*3^3\*5^2\*7^2\*11\*13\*17\*19\*29\*31\*37\*41\*43\*61\*71\*73$
$\*113\*127\*181\*211\*241\*281\*337\*421\*631\*1009\*2521$

Primality Testing and Jacobi Sums.zip (920.6 KB, 下载次数: 6)

Primality Testing and Jacobi Sums.z02 (1000 KB, 下载次数: 1)

Primality Testing and Jacobi Sums.z01 (1000 KB, 下载次数: 1)

mpz_aprcl_v1.2.zip (108.4 KB, 下载次数: 0)

无心人 · 发表于 2019-9-27 15:08:33

素性测试从费马测试，到米勒罗宾测试，条件越来越严格，但是测试还是在$ \ZZ_n $中进行
而Frobenius为代表的一些测试，则把测试扩展到了代数数域的二次域甚至三次域等
Jacobi和的严格素性证明，则用到了分圆域，其实本质也是代数数域，只不过，比二次域或者三次域更复杂了

椭圆曲线则显得比较另类，用到了椭圆曲线的知识

其实这个跟整数分解一样，一开始是连分数，然后是二次域，最后是一般代数域，同样也用到了椭圆曲线的方法

账号		自动登录	找回密码
密码			欢迎注册

[讨论] 关于伪素数与素性测试的一些事情

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

点评

点评