R-Lehmann Primality Test 中的模数警告 [英] Modulus warning in R- Lehmann Primality Test

问题描述

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

I spent a little time hacking an R implementation of the lehmann primality test. The function design I borrowed from 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 左右，我就会收到一条看起来很糟糕的消息并且该函数停止正常工作:

For small primes it works but as soon as i get around ~30, i get a bad looking message and the function stops working correctly:

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

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

After some investigating i believe it has to do with floating point conversions. Very large numbers are rounded so that the mod function gives a bad response.

现在是问题.

1. 这是浮点问题吗?还是在我的实现中?
2. 是否有纯粹的 R 解决方案，或者 R 只是在这方面做得不好?

谢谢

解决方案:

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

After the great feedback and a hour reading about modular exponentiation algorithms i have a solution. first it is to make my own modular exponentiation function. The basic idea is that modular multiplication allows you calculate intermediate results. you can calculate the mod after each iteration, thus never getting a giant nasty number that swamps the 16-bit 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 中，整数将被正确表示为 2^53 - 1，大约为 9e15.并且术语 y^((n-1)/2) 即使对于小数字也很容易超过它.您必须通过不断平方 y 并取模数来计算 (y^((n-1)/2)) %% n .这对应于 (n-1)/2 的二进制表示.

Of course there is a problem with representing integers. In R integers will be represented correctly up to 2^53 - 1 which is about 9e15. And the term y^((n-1)/2) will exceed that even for small numbers easily. You will have to compute (y^((n-1)/2)) %% n by continually squaring y and taking the modulus. That corresponds to the binary representation of (n-1)/2.

即使是真实的"数论程序也是这样做的——请参阅维基百科关于模取幂"的条目.也就是说，应该提到的是，像 R(或 Matlab 和其他数值计算系统)这样的程序可能不是实现数论算法的合适环境，甚至可能不是小整数的竞技场.

Even the 'real' number theory programs do it like that -- see Wikipedia's entry on "modular exponentiation". That said it should be mentioned that programs like R (or Matlab and other systems for numerical computing) may not be a proper environment for implementing number theory algorithms, probably not even as playing fields with small integers.

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

The original package was incorrect You could utilize the function modpower() in package 'pracma' like this:

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 是该集合中唯一的质数:

The following test is successful as 1009 is the only prime in this set:

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

这篇关于R-Lehmann Primality Test 中的模数警告的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！