数学研发论坛

 找回密码
 欢迎注册
12
返回列表 发新帖
楼主: 无心人

[原创] (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算法

这玩意不如概率性算法好使
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-11-10 22:32:59 | 显示全部楼层
本帖最后由 uk702 于 2020-11-10 22:36 编辑
uk702 发表于 2020-11-10 22:29
我觉得大概可以从这几个方向来减少计算量:
第一步,减少总体素性测试的量,比如说偶数就不用测试,逢 5 ...


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


  1. using Primes
  2. function test(n)
  3.     a=[1, 7, 11, 13, 17, 19, 23, 29]

  4.     d=big(10)
  5.     for i=1:n
  6.         j = 0
  7.         
  8.         # 根据初始余数进行校准
  9.         r = d % 30
  10.         b=a.-r
  11.         
  12.         u = 1
  13.         while b[u] < 0
  14.             u+=1
  15.         end
  16.         
  17.         notfound = true
  18.         while notfound
  19.             s = u
  20.             while s <= 8
  21.                 if isprime(d+j+b[s])
  22.                     println(i, " ", j+b[s])
  23.                     notfound = false
  24.                     break
  25.                 end   
  26.                 s+=1
  27.             end
  28.             j+=30
  29.             u = 1
  30.         end
  31.         
  32.         d*=10
  33.     end   
  34. end
复制代码


@time test(1000)
#552.174426 seconds (2.27 M allocations: 207.799 MiB, 0.01% gc time)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-11-11 07:04:02 | 显示全部楼层
本帖最后由 uk702 于 2020-11-11 07:56 编辑

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

  1. using Primes

  2. a=[gcd(i, 3003)==1 for i=1:3013]

  3. function test(n)
  4.     d=big(100)
  5.     for i=2:n
  6.         j = 0
  7.         
  8.         r = Int64(d % 3003)
  9.         #println("r = ", r)
  10.         
  11.         while true
  12.             if r > 3003
  13.                 r -= 3003
  14.             end
  15.             
  16.             p, q, l, m = r+1, r+3, r+7, r+9
  17.             #println("p = ", p, ", q = ", q, ", l = ", l, ", m = ", m)
  18.             #println("p in a = ", p in a)
  19.             
  20.             if a[p] == 1 && isprime(d+j+1)
  21.                 println(i, " ", j+1)
  22.                 @goto nexti
  23.             end   

  24.             if a[q] == 1 && isprime(d+j+3)
  25.                 println(i, " ", j+3)
  26.                 @goto nexti
  27.             end   

  28.             if a[l] == 1 && isprime(d+j+7)
  29.                 println(i, " ", j+7)
  30.                 @goto nexti
  31.             end   

  32.             if a[m] == 1 && isprime(d+j+9)
  33.                 println(i, " ", j+9)
  34.                 @goto nexti
  35.             end

  36.             j+=10
  37.             r+=10
  38.         end

  39. @label nexti
  40.         d*=10
  41.     end
  42. end

  43. @time test(1000)
复制代码


#503.741822 seconds (2.69 M allocations: 165.689 MiB, 0.01% gc time)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-11-11 07:11:34 | 显示全部楼层

微小的bug
因为\(1^0+1=2\),所以首行输出应当是0:1
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-11-11 15:21:02 | 显示全部楼层
多核并行计算版本,16核的机器起16条线程,82秒就算出来了。

  1. # set JULIA_NUM_THREADS=4

  2. using Base.Threads
  3. using Primes

  4. const a=[gcd(i, 3003)==1 for i=1:3013];
  5. result = Channel{Tuple}(1000);

  6. function one(i, d)
  7.         j = 0
  8.         r = d % 3003
  9.        
  10.         while true
  11.                 if r > 3003
  12.                         r -= 3003
  13.                 end
  14.                
  15.                 p, q, l, m = r+1, r+3, r+7, r+9
  16.                
  17.                 if a[p] == 1 && isprime(d+j+1)
  18.            return i, j+1
  19.                 end   

  20.                 if a[q] == 1 && isprime(d+j+3)
  21.             return i, j+3
  22.                 end   

  23.                 if a[l] == 1 && isprime(d+j+7)
  24.             return i, j+7
  25.                 end   
  26.         
  27.                 if a[m] == 1 && isprime(d+j+9)
  28.             return i, j+9
  29.                 end

  30.                 j+=10
  31.                 r+=10
  32.         end
  33. end;

  34. t1=time();
  35. n = 1000;
  36. @threads for i=2:n
  37.     put!(result, one(i, big(10)^i))
  38. end
  39.    
  40. while n > 1
  41.    i, p = take!(result)
  42.    println(i, " ", p)
  43.    global n = n - 1
  44. end;

  45. println("time = ", time() - t1)
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-11-14 10:14:03 | 显示全部楼层
经阅读 Primes.jl 的代码并经过实际测试,isprime() 函数内置了小素数因子测试,可以确认,#14、#15 给出的优化其实是没有效果的(最多是几条指令执行时间的区别)。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2021-1-26 11:20 , Processed in 0.064358 second(s), 14 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表