我花了一点时间破解莱曼素性测试的 R 实现。我借鉴的功能设计http://davidkendal.net/articles/2011/12/lehmann-primality-test 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 只是在这方面表现不佳?
Thanks
解决方案:
经过良好的反馈和一个小时的有关模幂算法的阅读后,我有了一个解决方案。首先是制作我自己的模幂函数。基本思想是模乘允许您计算中间结果。您可以在每次迭代后计算 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