R-Lehmann Primality Test 中的模量警告答案

【问题标题】：Modulus warning in R- Lehmann Primality TestR-Lehmann Primality Test 中的模量警告
【发布时间】：2011-12-20 19:15:05
【问题描述】：

我花了一点时间破解了 lehmann 素数测试的 R 实现。我借用http://davidkendal.net/articles/2011/12/lehmann-primality-test的功能设计

这是我的代码：

primeTest <- function(n, iter){
  a <- sample(1:(n-1), 1)
    lehmannTest <- function(y, tries){
    x <- ((y^((n-1)/2)) %% n)
    if (tries == 0) {
      return(TRUE)
            }else{          
      if ((x == 1) | (x == (-1 %% n))){
        lehmannTest(sample(1:(n-1), 1), (tries-1))
      }else{
    return(FALSE)
      }
    }
  }
  lehmannTest(a, iter)
}

primeTest(4, 50) # false
primeTest(3, 50) # true
primeTest(10, 50)# false
primeTest(97, 50) # gives false # SHOULD BE TRUE !!!! WTF

prime_test<-c(2,3,5,7,11,13,17 ,19,23,29,31,37)

for (i in 1:length(prime_test)) {
  print(primeTest(prime_test[i], 50))
}

对于小的素数，它可以工作，但是当我大约 30 左右时，我会收到一条看起来很糟糕的消息，并且该功能停止正常工作：

2: In lehmannTest(a, iter) : probable complete loss of accuracy in modulus

经过一些调查，我认为这与浮点转换有关。非常大的数字会四舍五入，因此 mod 函数会给出错误的响应。

现在是问题。

这是一个浮点问题吗？还是在我的实现中？
是否有纯粹的 R 解决方案，或者 R 在这方面做得不好？

谢谢

解决方案：

在获得了很好的反馈和一个小时的关于模幂算法的阅读之后，我有了一个解决方案。首先是制作我自己的模幂函数。基本思想是模乘允许您计算中间结果。您可以在每次迭代后计算 mod，因此永远不会得到一个淹没 16 位 R int 的巨大讨厌的数字。

modexp<-function(a, b, n){
    r = 1
    for (i in 1:b){
        r = (r*a) %% n
    }
    return(r)
}


primeTest <- function(n, iter){
   a <- sample(1:(n-1), 1)
    lehmannTest <- function(y, tries){
      x <- modexp(y, (n-1)/2, n)   
    if (tries == 0) {
      return(TRUE)
            }else{          
      if ((x == 1) | (x == (-1 %% n))){
        lehmannTest(sample(1:(n-1), 1), (tries-1))
        }else{
        return(FALSE)
         }
    }
  }
   if( n < 2 ){
     return(FALSE)
     }else if (n ==2) {
       return(TRUE)
       } else{
         lehmannTest(a, iter)
         }
}

primeTest(4, 50) # false
primeTest(3, 50) # true
primeTest(10, 50)# false
primeTest(97, 50) # NOW IS TRUE !!!!


prime_test<-c(5,7,11,13,17 ,19,23,29,31,37,1009)

for (i in 1:length(prime_test)) {
  print(primeTest(prime_test[i], 50))
}
#ALL TRUE

【问题讨论】：

标签： r primes

【解决方案1】：

当然，表示整数是有问题的。在 R 中，整数将正确表示为 2^53 - 1，大约为 9e15。并且y^((n-1)/2) 一词即使对于小数字也很容易超过这个值。您必须通过不断地对y 求平方并取模来计算(y^((n-1)/2)) %% n。这对应于(n-1)/2 的二进制表示。

即使是“实数”数论程序也是如此——参见维基百科关于“模幂运算”的条目。也就是说，应该提到像 R（或 Matlab 和其他数值计算系统）这样的程序可能不是实现数论算法的合适环境，甚至可能不是小整数的游戏场。

编辑：原始包不正确您可以像这样使用包 'pracma' 中的函数 modpower()：

primeTest <- function(n, iter){
  a <- sample(1:(n-1), 1)
    lehmannTest <- function(y, tries){
    x <- modpower(y, (n-1)/2, n)  # ((y^((n-1)/2)) %% n)
    if (tries == 0) {
      return(TRUE)
            }else{          
      if ((x == 1) | (x == (-1 %% n))){
        lehmannTest(sample(1:(n-1), 1), (tries-1))
      }else{
    return(FALSE)
      }
    }
  }
  lehmannTest(a, iter)
}

以下测试成功，因为 1009 是该集合中唯一的素数：

prime_test <- seq(1001, 1011, by = 2)
for (i in 1:length(prime_test)) {
    print(primeTest(prime_test[i], 50))
}
# FALSE FALSE FALSE FALSE TRUE  FALSE

【讨论】：

感谢您的帮助。我不能说这对工作很重要，但这项运动完全打败了我。但是现在我知道了模幂运算，我可以快乐地死去。在研究它之后，我使用 pow() 函数在 python 中重写了该函数。我很高兴在 R 中有一个实现。
此解决方案适用于某些数字。但是，当指数不是自然数时，modpower 函数会爆炸。这是包的来源。指数必须是自然数：floor(k) == ceiling(k)。当 n=4 时，(n-1)/2 = 1.5 并且 modpower 函数失败。

【解决方案2】：

如果你只是使用基础 R，我会选择 #2b...“R 不擅长这个”。在 R 中，整数（您似乎没有使用）被限制为 16 位精度。超过该限制，您将获得舍入错误。您可能应该查看：package:gmp 或 package:Brobdingnag。 Package:gmp 有大整数类和大有理类。

【讨论】：