10^28开始的10^14个整数的素性概率性测试算法实践报告
打算在这个范围内测试，基2伪素数PSP(2)，基2强伪素数SPSP(2)的分布 println("==========")
@time for i in range(1, step = 2, stop = 9999999)
test = b + i
if fermatProbablePrimeTest(test, 2)
if millerRabinProbablePrimeTest(test, 2)
if !isprime(test)
println("$test is SPSP(2)")
end
else
println("$test is PSP(2)")
end
end
end
println("==========")
用这种方式，测试的每10^7需要34.87s
10^14需要4035天，还需要优化~~~ 下面打算用筛法压缩下平均时间~~·否则太漫长了~~ 无心人 发表于 2020-11-26 16:18
用这种方式，测试的每10^7需要34.87s
10^14需要4035天，还需要优化~~~
是Julia吗？
能放一下完整代码吗？
没有完整代码优化起来并不容易……
特别是大家CPU不一样的情况下，我们很难判定究竟是优化还是CPU起了作用 10^28 就必须的吗？如果不是，建议减少到 2^64 范围内，这样就能在机器的字长内进行运算，而不需要采用 BigInt 了。
但这样要降到 10^19 以下，可能与你的需要有较大的差距。
======================================================
实际测试两者似相差不大，可能即使降到 10^19 以下，但实际计算中用到的数还是远超出 2^64 范围了。
另外，我实验表明，julia 中，就普遍性而言，isprime() 远比计算 powermod(2, n-1, n) 要快得多。
文件：a.jl
using Primes
t = 0
function test()
a = 10^18+1
for i=2:2:10^7
if isprime(a+i)
global t += 1
end
end
end
@time test()
println(t)
#=====================================================
$ julia a.jl
7.429383 seconds (240.78 k allocations: 3.674 MiB)
241295
=====================================================#
文件：b.jl
using Primes
t = 0
function test(a)
a = big(10)^28+1
for i=2:2:10^7
if isprime(a+i)
global t += 1
end
end
end
@time test()
println(t)
#=====================================================
$ julia a.jl
8.237764 seconds (25.31 M allocations: 831.328 MiB, 4.50% gc time)
155444
=====================================================# uk702 发表于 2020-11-26 17:54
10^28 就必须的吗？如果不是，建议减少到 2^64 范围内，这样就能在机器的字长内进行运算，而不需要采用 Big ...
2^64以内的数据，可以从网上查到，不需要我自己找的~ 我打算这么做筛法
首先不考虑所有偶数
对某个素数p，基2的次数是k，显然k｜p-1
那么如果mp是PSP(2), 显然m=nk+1
即，只需要测试所有p(nk+1)即可
比如，对于31，2^5=1 (mod 31），次数是5
素因子包含31的合数，如果是PSP(2)必然有形式31(5k+1)
那么我可以对31(5k+3)的整数加不是PSP(2)的标志
进一步，因为是奇数所以形式31(10k+3, 7, 9)可以都标注非PSP(2)
对小于某个值的所有小素数都进性这个操作，然后再测试未标注的整数即可 本帖最后由 uk702 于 2020-11-28 09:00 编辑
如果不涉及保密问题，建议将代码放出来大家一起研究改进。如果有保密问题，我想大家也能理解，毕竟是正在进行的科研项目。
建议认真研究一下 Primes.jl ，其中的 isprime() 确实很快，效率真的是高。
以 10^28 为例，这时平均每 65 个自然数里头就有一个真素数， 1000 个数里面，大约有 15 个是真素数。
即使我们用小因子排除法，能被直接排除的很难超过4/5，
也就是说，有800个数我们能够直接排除，
但不幸的是，这 800 个数中，isprime() 本身就不费吹灰之力，也就是说，基本上没有多少提升的余地。
还有 200 个数，包括其中15个真素数，我们难以判断它们是素数还是合数，如果这一部分无法改进，
那天量计算恐难以避免，
如果真的要验证 10^14 这个数量，按照目前的计算能力，基本上可以肯定，
假设是PC机的计算能力的话，必将是数以百天计、千天计的计算量，
建议你们尝试找计算中心，看能否搞上 1000 个核同时计算。
或者找超算中心，你们的课题有一些的军事价值，还是有这种可能的。
而且，你们的题目做并行计算其实简单没难度。 using Primes
struct Quadratic{T <: Integer, C} <: Integer
re::T
ir::T
end
function qfmul(left:: Quadratic{T, C}, right:: Quadratic{T, C}) where {T <: Integer, C}
re = left.re * right.re + left.ir * right.ir * C
ir = left.re * right.ir + left.ir * right.re
return Quadratic{T, C}(re, ir)
end
function qfmod(left:: Quadratic{T, C}, r:: T) where {T <: Integer, C}
re = left.re % r
ir = left.ir % r
return Quadratic{T, C}(re, ir)
end
function qfmulmod(left:: Quadratic{T, C}, right:: Quadratic{T, C}, r::T) where {T <: Integer, C}
re = (left.re * right.re + left.ir * right.ir * C) % r
ir = (left.re * right.ir + left.ir * right.re) % r
return Quadratic{T, C}(re, ir)
end
function qfpowermod(x::Quadratic{T, C}, p::T , m::T) where {T <: Integer, C}
@assert p >= 0
p == 0 && return Quadratic{T, C}(1, 0)
#b = oftype(m,mod(x,m))# this also checks for divide by zero
b = x
t = prevpow(2, p)
r = Quadratic{T, C}(1, 0)
while true
if p >= t
r = qfmulmod(r, b, m)
p -= t
end
t >>>= 1
t <= 0 && break
r = qfmulmod(r, r ,m)
end
return r
end
function jacobi(a, b)
tab2 =
if b == 0
if abs(a)==1
return 1
else
return 0
end
end
if iseven(a) && iseven(b)
return 0
end
v = trailing_zeros(b)
b = b >> v
if iseven(v)
k = 1
else
k = tab2
end
if b < 0
b = -b
end
if a < 0
k = -k
end
while a != 0
v = trailing_zeros(a)
a = a >> v
if isodd(v)
k *= tab2
end
if a & b & 3 == 3
k = -k
end
r = abs(a)
a = b % r
b = r
end
if b > 1
return 0
else
return k
end
end
function frobenius(n)
small = [
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,131,137
]
find = false
global p = 0
for c in small
if jacobi(c, n) == -1
find = true
global p = c
break
end
end
if !find
c = Primes.nextprime(last(small)+1)
while jacobi(c, n) != -1
c = Primes.nextprime(c+1)
end
global p = c
end
#println("$n: $p")
b = Quadratic{typeof(n), p}(1, 1)
r = qfpowermod(b, n, n)
return (r.re == 1) && (r.ir == n - 1)
end
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
function millerRabinProbablePrimeTest( n, b )
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
global t = 0
global b = big(10)^28
println("==========")
@time for i in range(1, step = 2, stop = 9999999)
test = b + i
# if isqrt(test)^2 == test
# continue
# end
if fermatProbablePrimeTest(test, 2)
if millerRabinProbablePrimeTest(test, 2)
#if frobenius(test)
if !isprime(test)
#println("$test is frobenius pseudoprime!")
#else
#println("$test is prime")
#end
#else
println("$test is SPSP(2)")
end
else
println("$test is PSP(2)")
end
end
end
println("==========")
本帖最后由 .·.·. 于 2020-11-28 12:45 编辑
无心人 发表于 2020-11-28 10:53
fermat可以省略
if gcd(n, b) > 1
return false
end
这一步
我本来觉得fermat可以省略，毕竟millerRabin做了fermat该做的一切
但算起来省略fermat的millerRabin反而更慢了。
原因是我们做了太多不必要的判断
这里，gcd其实是不必要的判断，大多数情况下gcd都是不必要的
BTW让我优化的话……我只会改分布式玩并行
using Distributed
addprocs(6)
@everywhere begin
using Primes
const b = big(10)^28
global function fermatProbablePrimeTest2(n)
a = powermod(2, (n-1)>>1, n)
return (a == 1) || (a == (n-1))
end
global function millerRabinProbablePrimeTest( n, b )
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
end
@time begin
x=@distributed for i in range(1, step = 2, stop = 9999999)
test = b + i
# if isqrt(test)^2 == test
# continue
# end
if fermatProbablePrimeTest2(test)
if millerRabinProbablePrimeTest(test, 2)
#if frobenius(test)
if !isprime(test)
#println("$test is frobenius pseudoprime!")
#else
#println("$test is prime")
#end
#else
println("$test is SPSP(2)")
end
else
println("$test is PSP(2)")
end
end
end
wait(x)
end