数学研发论坛

 找回密码
 欢迎注册
查看: 2790|回复: 116

[原创] 10^28开始的10^14个整数的素性概率性测试算法实践报告

[复制链接]
发表于 2020-11-26 16:11:32 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
打算在[10^28, 10^28+10^14]这个范围内
测试,基2伪素数PSP(2),基2强伪素数SPSP(2)的分布
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2020-11-26 16:18:29 | 显示全部楼层
  1. println("==========")
  2. @time for i in range(1, step = 2, stop = 9999999)
  3.     test = b + i
  4.     if fermatProbablePrimeTest(test, 2)
  5.         if millerRabinProbablePrimeTest(test, 2)
  6.             if !isprime(test)
  7.                 println("$test is SPSP(2)")
  8.             end
  9.         else
  10.             println("$test is PSP(2)")
  11.         end
  12.     end
  13. end
  14. println("==========")
复制代码

用这种方式,测试的每10^7需要34.87s
10^14需要4035天,还需要优化~~~

点评

现在可以优化到~11s多CPU有望在半年解决战斗. ispime调包方法太成问题了  发表于 2020-11-30 03:21
先筛就失去了意义,楼主的目的就是查找这些伪素数  发表于 2020-11-27 07:48
可以先拿小素数筛一下  发表于 2020-11-26 16:30
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2020-11-26 16:31:45 | 显示全部楼层
下面打算用筛法压缩下平均时间~~·否则太漫长了~~
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-11-26 16:34:17 | 显示全部楼层
无心人 发表于 2020-11-26 16:18
用这种方式,测试的每10^7需要34.87s
10^14需要4035天,还需要优化~~~

是Julia吗?
能放一下完整代码吗?

没有完整代码优化起来并不容易……
特别是大家CPU不一样的情况下,我们很难判定究竟是优化还是CPU起了作用
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-11-26 17:54:41 | 显示全部楼层
10^28 就必须的吗?如果不是,建议减少到 2^64 范围内,这样就能在机器的字长内进行运算,而不需要采用 BigInt 了。
但这样要降到 10^19 以下,可能与你的需要有较大的差距。

======================================================
实际测试两者似相差不大,可能即使降到 10^19 以下,但实际计算中用到的数还是远超出 2^64 范围了。
另外,我实验表明,julia 中,就普遍性而言,isprime() 远比计算 powermod(2, n-1, n) 要快得多。

  1. 文件:a.jl
  2. using Primes
  3. t = 0
  4. function test()
  5.         a = 10^18+1
  6.         for i=2:2:10^7
  7.                 if isprime(a+i)
  8.                         global t += 1
  9.                 end
  10.         end
  11. end

  12. @time test()
  13. println(t)       

  14. #=====================================================
  15. $ julia a.jl
  16.   7.429383 seconds (240.78 k allocations: 3.674 MiB)
  17. 241295
  18. =====================================================#

  19. 文件:b.jl
  20. using Primes
  21. t = 0
  22. function test(a)
  23.         a = big(10)^28+1
  24.         for i=2:2:10^7
  25.                 if isprime(a+i)
  26.                         global t += 1
  27.                 end
  28.         end
  29. end

  30. @time test()
  31. println(t)       

  32. #=====================================================
  33. $ julia a.jl                                                         
  34.   8.237764 seconds (25.31 M allocations: 831.328 MiB, 4.50% gc time)
  35. 155444                                                               
  36. =====================================================#
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2020-11-27 22:49:23 | 显示全部楼层
uk702 发表于 2020-11-26 17:54
10^28 就必须的吗?如果不是,建议减少到 2^64 范围内,这样就能在机器的字长内进行运算,而不需要采用 Big ...

2^64以内的数据,可以从网上查到,不需要我自己找的~
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2020-11-27 22:55:46 | 显示全部楼层
我打算这么做筛法
首先不考虑所有偶数

对某个素数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)

对小于某个值的所有小素数都进性这个操作,然后再测试未标注的整数即可
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-11-28 07:43:39 | 显示全部楼层
本帖最后由 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 个核同时计算。
或者找超算中心,你们的课题有一些的军事价值,还是有这种可能的。
而且,你们的题目做并行计算其实简单没难度。

点评

个人爱好做的项目  发表于 2020-11-28 10:52
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2020-11-28 10:53:49 | 显示全部楼层
  1. using Primes

  2. struct Quadratic{T <: Integer, C} <: Integer
  3.    re::T
  4.    ir::T
  5. end

  6. function qfmul(left:: Quadratic{T, C}, right:: Quadratic{T, C}) where {T <: Integer, C}
  7.     re = left.re * right.re + left.ir * right.ir * C
  8.     ir = left.re * right.ir + left.ir * right.re
  9.     return Quadratic{T, C}(re, ir)
  10. end

  11. function qfmod(left:: Quadratic{T, C}, r:: T) where {T <: Integer, C}
  12.     re = left.re % r
  13.     ir = left.ir % r
  14.     return Quadratic{T, C}(re, ir)
  15. end

  16. function qfmulmod(left:: Quadratic{T, C}, right:: Quadratic{T, C}, r::T) where {T <: Integer, C}
  17.     re = (left.re * right.re + left.ir * right.ir * C) % r
  18.     ir = (left.re * right.ir + left.ir * right.re) % r
  19.     return Quadratic{T, C}(re, ir)
  20. end

  21. function qfpowermod(x::Quadratic{T, C}, p::T , m::T) where {T <: Integer, C}
  22.     @assert p >= 0
  23.     p == 0 && return Quadratic{T, C}(1, 0)
  24.     #b = oftype(m,mod(x,m))  # this also checks for divide by zero
  25.     b = x
  26.     t = prevpow(2, p)
  27.     r = Quadratic{T, C}(1, 0)
  28.     while true
  29.         if p >= t
  30.             r = qfmulmod(r, b, m)
  31.             p -= t
  32.         end
  33.         t >>>= 1
  34.         t <= 0 && break
  35.         r = qfmulmod(r, r ,m)
  36.     end
  37.     return r
  38. end

  39. function jacobi(a, b)
  40.     tab2 = [0, 1, 0, -1, 0, -1, 0, 1]
  41.     if b == 0
  42.         if abs(a)==1
  43.             return 1
  44.         else
  45.             return 0
  46.         end
  47.     end
  48.     if iseven(a) && iseven(b)
  49.         return 0
  50.     end
  51.     v = trailing_zeros(b)
  52.     b = b >> v
  53.     if iseven(v)
  54.         k = 1
  55.     else
  56.         k = tab2[a & 7 + 1]
  57.     end
  58.     if b < 0
  59.         b = -b
  60.     end
  61.     if a < 0
  62.         k = -k
  63.     end
  64.     while a != 0
  65.         v = trailing_zeros(a)
  66.         a = a >> v
  67.         if isodd(v)
  68.             k *= tab2[b & 7 + 1]
  69.         end
  70.         if a & b & 3 == 3
  71.             k = -k
  72.         end
  73.         r = abs(a)
  74.         a = b % r
  75.         b = r
  76.     end
  77.     if b > 1
  78.         return 0
  79.     else
  80.         return k
  81.     end
  82. end

  83. function frobenius(n)
  84.     small = [
  85.                3,   5,  7, 11, 13, 17, 19, 23,
  86.               29,  31, 37,  41, 43, 47, 53, 59,
  87.               61,  67, 71, 73, 79, 83, 89, 97,
  88.               101,103,107,109,113,127,131,137
  89.             ]
  90.     find = false
  91.     global p = 0
  92.     for c in small
  93.         if jacobi(c, n) == -1
  94.             find = true
  95.             global p = c
  96.             break
  97.         end
  98.     end
  99.     if !find
  100.         c = Primes.nextprime(last(small)+1)
  101.         while jacobi(c, n) != -1
  102.             c = Primes.nextprime(c+1)
  103.         end
  104.         global p = c
  105.     end
  106.     #println("$n: $p")
  107.     b = Quadratic{typeof(n), p}(1, 1)
  108.     r = qfpowermod(b, n, n)
  109.     return (r.re == 1) && (r.ir == n - 1)
  110. end

  111. function fermatProbablePrimeTest(n, b)
  112.     if gcd(n, b) > 1
  113.         return false
  114.     end
  115.     a = powermod(b, (n-1)>>1, n)
  116.     return (a == 1) || (a == (n-1))
  117. end

  118. function millerRabinProbablePrimeTest( n, b )
  119.     k = trailing_zeros(n-1)
  120.     m = (n-1) >> k
  121.     a = powermod(b, m, n)
  122.     if (a == 1) || (a == (n-1))
  123.         return true
  124.     else
  125.         for i in 1:k-1
  126.             a = powermod(a, 2, n)
  127.             if a == n - 1
  128.                 return true
  129.             end
  130.         end
  131.     end
  132.     return false
  133. end

  134. global t = 0
  135. global b = big(10)^28

  136. println("==========")
  137. @time for i in range(1, step = 2, stop = 9999999)
  138.     test = b + i
  139. #    if isqrt(test)^2 == test
  140. #        continue
  141. #    end

  142.     if fermatProbablePrimeTest(test, 2)
  143.         if millerRabinProbablePrimeTest(test, 2)
  144.             #if frobenius(test)
  145.                 if !isprime(test)
  146.                     #println("$test is frobenius pseudoprime!")
  147.                 #else
  148.                     #println("$test is prime")
  149.                 #end
  150.             #else
  151.                 println("$test is SPSP(2)")
  152.             end
  153.         else
  154.             println("$test is PSP(2)")
  155.         end
  156.     end
  157. end
  158. println("==========")
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-11-28 12:21:41 | 显示全部楼层
本帖最后由 .·.·. 于 2020-11-28 12:45 编辑

fermat可以省略
    if gcd(n, b) > 1
        return false
    end
这一步

我本来觉得fermat可以省略,毕竟millerRabin做了fermat该做的一切
但算起来省略fermat的millerRabin反而更慢了。

原因是我们做了太多不必要的判断
这里,gcd其实是不必要的判断,大多数情况下gcd都是不必要的

BTW让我优化的话……我只会改分布式玩并行
  1. using Distributed
  2. addprocs(6)
  3. @everywhere begin
  4.     using Primes
  5.     const b = big(10)^28

  6.     global function fermatProbablePrimeTest2(n)
  7.         a = powermod(2, (n-1)>>1, n)
  8.         return (a == 1) || (a == (n-1))
  9.     end

  10.     global function millerRabinProbablePrimeTest( n, b )
  11.         k = trailing_zeros(n-1)
  12.         m = (n-1) >> k
  13.         a = powermod(b, m, n)
  14.         if (a == 1) || (a == (n-1))
  15.             return true
  16.         else
  17.             for i in 1:k-1
  18.                 a = powermod(a, 2, n)
  19.                 if a == n - 1
  20.                     return true
  21.                 end
  22.             end
  23.         end
  24.         return false
  25.     end
  26. end


  27. @time begin
  28.     x=@distributed for i in range(1, step = 2, stop = 9999999)
  29.         test = b + i
  30.     #    if isqrt(test)^2 == test
  31.     #        continue
  32.     #    end

  33.         if fermatProbablePrimeTest2(test)
  34.             if millerRabinProbablePrimeTest(test, 2)
  35.                 #if frobenius(test)
  36.                     if !isprime(test)
  37.                         #println("$test is frobenius pseudoprime!")
  38.                     #else
  39.                         #println("$test is prime")
  40.                     #end
  41.                 #else
  42.                     println("$test is SPSP(2)")
  43.                 end
  44.             else
  45.                 println("$test is PSP(2)")
  46.             end
  47.         end
  48.     end
  49.     wait(x)
  50. end
复制代码


点评

不,加速不了powermod,但是可以减少测试次数,我说的筛法就是一个尝试,因为符合条件的合数n如果有素因子p,对2的次数如果是k,那么余因子一定是mk+1形式,即n=p(mk+1)我们可以通过筛法,预先把不符合的n筛掉  发表于 2020-11-29 09:01
主要计算量在 powermod() 和 isprime() 上,如果不能加速 powermod(),我准备投降了。  发表于 2020-11-29 08:48
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2021-1-20 15:08 , Processed in 0.063342 second(s), 18 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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