Тест на примитивность MillerRabin в C#

Question

Добро пожаловать. Я пытаюсь выполнить тест MillerRabin для проверки того, является ли большое заданное число простым. Вот мой код:

public static bool MillerRabinTest(BigInteger number) 
     { 

      BigInteger d; 
      var n = number - 1; 
      var s = FindK(n, out d); 

      BigInteger a = 2; 
      BigInteger y = Calc(a, d, number); //a^d mod number 
      if (y != BigInteger.One && y != n) 
      { 
       for (var r = 1; r <= s - 1; r++) 
       { 
        y = Calc(y, 2, number); 
        if (y == 1) 
         return false; 
       } 

       if (y != n) 
        return false; 
      } 
      return true; //it is probably prime 
     } 

Он отлично подходит для небольших Bigintegers. Но если мои программы должны анализировать числа, содержащие более 16 бит, программа зависает. Например, после успешной проверки, если число является простым, программа внезапно не реагирует. Я не понимаю, как это возможно. Если он проверил одно большое число, у него не должно возникнуть проблемы с проверкой другого. Даже отладчик не помогает, потому что step options исчезают. При необходимости я могу предоставить больше кода функций. Выше функция работает правильно для небольших чисел.

EDIT. Изменение моей функции модуляции для BigInteger.ModPow помогло. К сожалению, теперь для больших чисел, более 3000 бит, он никогда не возвращает простое число, что является весьма невозможным. Или действительно цифры prme трудно узнать?

Answer 1

Ну, это занимает около 5 секунд на моей рабочей станции (ядро 3.2GHz i5, IA64 .Net 4.5) для проверки того простым для чисел равна 2**3000:

public static class PrimeExtensions { 
    // Random generator (thread safe) 
    private static ThreadLocal<Random> s_Gen = new ThreadLocal<Random>(
    () => { 
     return new Random(); 
     } 
    ); 

    // Random generator (thread safe) 
    private static Random Gen { 
     get { 
     return s_Gen.Value; 
     } 
    } 

    public static Boolean IsProbablyPrime(this BigInteger value, int witnesses = 10) { 
     if (value <= 1) 
     return false; 

     if (witnesses <= 0) 
     witnesses = 10; 

     BigInteger d = value - 1; 
     int s = 0; 

     while (d % 2 == 0) { 
     d /= 2; 
     s += 1; 
     } 

     Byte[] bytes = new Byte[value.ToByteArray().LongLength]; 
     BigInteger a; 

     for (int i = 0; i < witnesses; i++) { 
     do { 
      Gen.NextBytes(bytes); 

      a = new BigInteger(bytes); 
     } 
     while (a < 2 || a >= value - 2); 

     BigInteger x = BigInteger.ModPow(a, d, value); 
     if (x == 1 || x == value - 1) 
      continue; 

     for (int r = 1; r < s; r++) { 
      x = BigInteger.ModPow(x, 2, value); 

      if (x == 1) 
      return false; 
      if (x == value - 1) 
      break; 
     } 

     if (x != value - 1) 
      return false; 
     } 

     return true; 
    } 
    } 

испытаний и тестов

BigInteger value = BigInteger.Pow(2, 3217) - 1; // Mersenne prime number (2.5e968) 

    Stopwatch sw = new Stopwatch(); 

    sw.Start(); 

    Boolean isPrime = value.IsProbablyPrime(10); 

    sw.Stop(); 

    Console.Write(isPrime ? "probably prime" : "not prime"); 
    Console.WriteLine(); 
    Console.Write(sw.ElapsedMilliseconds);