R-萊曼原始性測試中的模量警告

我花了一點時間來攻擊萊曼性能測試的R實現。該功能的設計，我從http://davidkendal.net/articles/2011/12/lehmann-primality-test R-萊曼原始性測試中的模量警告

這裏借是我的代碼：

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

來源

2011-12-20 jdennison

當然，表示整數有一個問題。在R中，整數將正確地表示爲2^53 - 1，大約是9e15。即使對於小數字，術語y^((n-1)/2)也會超出這個範圍。你將不得不計算(y^((n-1)/2)) %% n不斷平方y並取模數。這對應於(n-1)/2的二進制表示。

即使是'真正'的數字理論程序也是這樣做的 - 請參閱維基百科關於「模數求冪」的條目。也就是說應該提到像R（或者Matlab和其他數字計算系統）這樣的程序可能不適合實現數論算法的適當環境，甚至可能不是具有小整數的遊戲場。

編輯：原包裝是不正確的你可以利用功能modpower（）封裝內「pracma」是這樣的：

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

來源

2011-12-20 20:09:08

感謝您的幫助。我不能說這對工作是至關重要的，但是這個練習完全打敗了我。但是現在我知道模數求冪，我可以死一個快樂的人。在研究之後，我使用pow（）函數重寫了python中的函數。我很高興R中有一個實現。 – jdennison 2011-12-21 15:12:19

該解決方案適用於某些數字。但是，當指數不是自然數時，modpower函數會爆炸。這裏是包的來源。指數必須是一個自然數：floor（k）== ceiling（k）。當n = 4時，（n-1）/ 2 = 1.5並且modpower功能失敗。 – jdennison 2011-12-21 17:25:13

如果你只是使用base R，我會選擇＃2b ......「R在這個方面很糟糕」。在R整數（你似乎沒有使用）被限制爲16位精度。超過這個限制，你會得到舍入誤差。你應該看看：package：gmp或者package：Brobdingnag。包：gmp有大整數和大有理類。

來源

2011-12-20 19:39:14

R-萊曼原始性測試中的模量警告

回答

相關問題