（n+1）位十进制整数中的最小素数

.·.·. · 发表于 2020-11-9 18:01:56

mathe 发表于 2020-11-9 16:59
确定性素性测试有AKS算法

我记得pari/gp专门讨论过aks
似乎是，在可接受的范围内，aks并不是最快的算法

pari/gp自带一个primecert
虽然不懂原理但据说这个够快

无心人 · 发表于 2020-11-9 18:59:21

mathe 发表于 2020-11-9 16:59
确定性素性测试有AKS算法

这玩意不如概率性算法好使

uk702 · 发表于 2020-11-10 22:29:58

无心人发表于 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:32:59

本帖最后由 uk702 于 2020-11-10 22:36 编辑

uk702 发表于 2020-11-10 22:29
我觉得大概可以从这几个方向来减少计算量：
第一步，减少总体素性测试的量，比如说偶数就不用测试，逢 5 ...

尝试根据模 30 做一个表，过滤掉 2、3、5 的倍数，测试并没带来收益。
这样，步骤一基本已经到头了。

using Primes
function test(n)
a=[1, 7, 11, 13, 17, 19, 23, 29]
d=big(10)
for i=1:n
j = 0
# 根据初始余数进行校准
r = d % 30
b=a.-r
u = 1
while b[u] < 0
u+=1
end
notfound = true
while notfound
s = u
while s <= 8
if isprime(d+j+b[s])
println(i, " ", j+b[s])
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:04:02

本帖最后由 uk702 于 2020-11-11 07:56 编辑

综合上述步骤，小素因子（2、3、5、7、11、13）直接查表，大素因子用 isprime() 来判断，从 520 秒下降到 500 秒，似乎又快了 20 秒（多次测试差异比较大，另外的一台家用台式机，332 秒就完成计算）。
这个思路估计是做到头了，下一步再考虑并行计算。

using Primes
a=[gcd(i, 3003)==1 for i=1:3013]
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[p] == 1 && isprime(d+j+1)
println(i, " ", j+1)
@goto nexti
end
if a[q] == 1 && isprime(d+j+3)
println(i, " ", j+3)
@goto nexti
end
if a[l] == 1 && isprime(d+j+7)
println(i, " ", j+7)
@goto nexti
end
if a[m] == 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)

zeroieme · 发表于 2020-11-11 07:11:34

无心人发表于 2014-11-21 20:49
结果是

微小的bug
因为\(1^0+1=2\)，所以首行输出应当是0：1

uk702 · 发表于 2020-11-11 15:21:02

多核并行计算版本，16核的机器起16条线程，82秒就算出来了。

# set JULIA_NUM_THREADS=4
using Base.Threads
using Primes
const a=[gcd(i, 3003)==1 for i=1:3013];
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[p] == 1 && isprime(d+j+1)
return i, j+1
end
if a[q] == 1 && isprime(d+j+3)
return i, j+3
end
if a[l] == 1 && isprime(d+j+7)
return i, j+7
end
if a[m] == 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)

复制代码

uk702 · 发表于 2020-11-14 10:14:03

经阅读 Primes.jl 的代码并经过实际测试，isprime() 函数内置了小素数因子测试，可以确认，#14、#15 给出的优化其实是没有效果的（最多是几条指令执行时间的区别）。

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

[原创] （n+1）位十进制整数中的最小素数