确定性素性测试有AKS算法
我记得pari/gp专门讨论过aks
似乎是，在可接受的范围内，aks并不是最快的算法
pari/gp自带一个primecert
虽然不懂原理但据说这个够快 mathe 发表于 2020-11-9 16:59
确定性素性测试有AKS算法
这玩意不如概率性算法好使
无心人 发表于 2014-11-21 21:13
原来想做个大的，计算到10000位，后来发现，想证明那么大的数是素数很难，
而到1000位，实际证明也是相当 ...
我觉得大概可以从这几个方向来减少计算量：
第一步，减少总体素性测试的量，比如说偶数就不用测试，逢 5 结尾不用测试。这样，只需对以1、3、7、9 结尾的进行测试。
第二步，减少某个具体大整数的素性测试的计算量，这个太难，我感觉没什么方法可以写出比库函数 isprime() 还要高效的方法，但我们可以尝试先试除 3、7、11、13 等小素数，如果都不整除，再执行 isprime()。
第一步想法，使用 julia 语言，对前 n=1:1000 进行测试，耗时 521 秒，代码如下：
using Primes
function test(n)
d=big(10)
for i=1:n
j = 0
while true
if isprime(d+j+1)
println(i, " ", j+1)
@goto nexti
end
if isprime(d+j+3)
println(i, " ", j+3)
@goto nexti
end
if isprime(d+j+7)
println(i, " ", j+7)
@goto nexti
end
if isprime(d+j+9)
println(i, " ", j+9)
@goto nexti
end
j+=10
end
@label nexti
d*=10
end
end
@time test(1000) 本帖最后由 uk702 于 2020-11-10 22:36 编辑
uk702 发表于 2020-11-10 22:29
我觉得大概可以从这几个方向来减少计算量：
第一步，减少总体素性测试的量，比如说偶数就不用测试，逢 5 ...
尝试根据模 30 做一个表，过滤掉 2、3、5 的倍数，测试并没带来收益。
这样，步骤一基本已经到头了。
using Primes
function test(n)
a=
d=big(10)
for i=1:n
j = 0
# 根据初始余数进行校准
r = d % 30
b=a.-r
u = 1
while b < 0
u+=1
end
notfound = true
while notfound
s = u
while s <= 8
if isprime(d+j+b)
println(i, " ", j+b)
notfound = false
break
end
s+=1
end
j+=30
u = 1
end
d*=10
end
end
@time test(1000)
#552.174426 seconds (2.27 M allocations: 207.799 MiB, 0.01% gc time) 本帖最后由 uk702 于 2020-11-11 07:56 编辑
综合上述步骤，小素因子（2、3、5、7、11、13）直接查表，大素因子用 isprime() 来判断，从 520 秒下降到 500 秒，似乎又快了 20 秒（多次测试差异比较大，另外的一台家用台式机，332 秒就完成计算）。
这个思路估计是做到头了，下一步再考虑并行计算。
using Primes
a=
function test(n)
d=big(100)
for i=2:n
j = 0
r = Int64(d % 3003)
#println("r = ", r)
while true
if r > 3003
r -= 3003
end
p, q, l, m = r+1, r+3, r+7, r+9
#println("p = ", p, ", q = ", q, ", l = ", l, ", m = ", m)
#println("p in a = ", p in a)
if a == 1 && isprime(d+j+1)
println(i, " ", j+1)
@goto nexti
end
if a == 1 && isprime(d+j+3)
println(i, " ", j+3)
@goto nexti
end
if a == 1 && isprime(d+j+7)
println(i, " ", j+7)
@goto nexti
end
if a == 1 && isprime(d+j+9)
println(i, " ", j+9)
@goto nexti
end
j+=10
r+=10
end
@label nexti
d*=10
end
end
@time test(1000)
#503.741822 seconds (2.69 M allocations: 165.689 MiB, 0.01% gc time) 无心人 发表于 2014-11-21 20:49
结果是
微小的bug
因为\(1^0+1=2\)，所以首行输出应当是0：1 多核并行计算版本，16核的机器起16条线程，82秒就算出来了。
# set JULIA_NUM_THREADS=4
using Base.Threads
using Primes
const a=;
result = Channel{Tuple}(1000);
function one(i, d)
j = 0
r = d % 3003
while true
if r > 3003
r -= 3003
end
p, q, l, m = r+1, r+3, r+7, r+9
if a == 1 && isprime(d+j+1)
return i, j+1
end
if a == 1 && isprime(d+j+3)
return i, j+3
end
if a == 1 && isprime(d+j+7)
return i, j+7
end
if a == 1 && isprime(d+j+9)
return i, j+9
end
j+=10
r+=10
end
end;
t1=time();
n = 1000;
@threads for i=2:n
put!(result, one(i, big(10)^i))
end
while n > 1
i, p = take!(result)
println(i, " ", p)
global n = n - 1
end;
println("time = ", time() - t1)
经阅读 Primes.jl 的代码并经过实际测试，isprime() 函数内置了小素数因子测试，可以确认，#14、#15 给出的优化其实是没有效果的（最多是几条指令执行时间的区别）。
页:
1
[2]