Предупреждение о модуле в тесте R-Lehmann Primality - PullRequest
Новая
       1

Предупреждение о модуле в тесте R-Lehmann Primality

3 голосов
/ 20 декабря 2011

Я потратил немного времени на взлом 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

После некоторых исследований я считаю, что это связано с преобразованиями с плавающей запятой. Очень большие числа округляются, так что функция мода дает плохой ответ.

Теперь вопросы.

  1. Это проблема с плавающей запятой? или в моей реализации?
  2. Есть ли чисто R решение или R просто плох в этом?

Спасибо

Решение:

После отличной обратной связи и часового прочтения о алгоритмах модульного возведения в степень у меня есть решение. Во-первых, это сделать мою собственную модульную функцию возведения в степень. Основная идея заключается в том, что модульное умножение позволяет вычислять промежуточные результаты. Вы можете вычислить мод после каждой итерации, таким образом, никогда не получая гигантское противное число, которое забивает 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

Ответы [ 2 ]

5 голосов
/ 21 декабря 2011

Конечно, есть проблема с представлением целых чисел. В 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
2 голосов
/ 20 декабря 2011

Если вы просто используете базу R, я бы выбрал # 2b ... "R это плохо". В R целые числа (которые вы, похоже, не используете) ограничены 16-битной точностью. Выше этого предела вы получите ошибки округления. Вероятно, вы должны посмотреть на: package: gmp или package: Brobdingnag. Пакет: gmp имеет большие целочисленные и большие рациональные классы.

...