Miller test (deterministic primality test)

miller_test.cpp

//////////////////////////////////////////////////////////
// Miller test, v1.10 (c) 2019 by Jin X (jin_x@list.ru) //
//////////////////////////////////////////////////////////

// ИСТОРИЯ ИЗМЕНЕНИЙ
// =================
// v1.10 (11.03.2019):
// [!] Добавлено получение остатка через умножение, в результате скорость возросла
//     (примерно в 2 раза или даже более), однако код заметно усложнился.
// [+] Добавлено несколько макро-идентификаторов для управления исходным кодом.
//
// v1.00 (07.03.2019):
// [!] Самая первая версия.

#include <iostream>
#include <ctime>

using std::cout;

////////////////////////////////////////////////////////////////////////////////////////////////////

// PRIME_ALLOW_INTRIN: Разрешить использование интринсиков (только для архитектуры x86/x86-64).
// Интринсики дают прирост производительности, но не все компиляторы их поддерживают.
// При использовании интринсиков происходит подключение заголовочного файла intrin.h.
#define PRIME_ALLOW_INTRIN

// PRIME_ALLOW_INT128: Разрешить использование типа __int128, если он поддерживается.
// Тип __int128 даёт небольшой прирост производительности.
#define PRIME_ALLOW_INT128

// PRIME_ARCH_X86: Архитектура, допустимые значения: 1 (не 0) - x86/x86-64, 0 - другие.
// Если идентификатор не определён, производится попытка определить значение автоматически
// (но не для всех компиляторов и платформ, поэтому рекомендуется определить его).
//#define PRIME_ARCH_X86 1

// PRIME_CODE_SIZE: Разрядность кода, значение должно содержать кол-во бит: 32 или 64.
// Если идентификатор не определён, производится попытка определить значение автоматически
// (но не для всех компиляторов и платформ, поэтому рекомендуется определить его).
//#define PRIME_CODE_SIZE 64

////////////////////////////////////////////////////////////////////////////////////////////////////

// Автоматическое определение архитектур x86/x86-64 при отсутствии идентификатора PRIME_ARCH_X86.
#if !defined(PRIME_ARCH_X86)
  #if defined(__i386__) || defined(_M_IX32) || defined(__x86_64__) || defined(_M_X64)
    #define PRIME_ARCH_X86 1
  #else
    #define PRIME_ARCH_X86 0
  #endif
#endif

// Автоматическое определение разрядности кода при отсутствии идентификатора PRIME_CODE_SIZE.
#if !defined(PRIME_CODE_SIZE)
  #if defined(__x86_64__) || defined(_M_X64) || defined(__ia64__) || defined(_M_IA64) || \
      defined(__aarch64__) || defined(_M_ARM64) || defined(__ppc64__)
    #define PRIME_CODE_SIZE 64
  #else
    #define PRIME_CODE_SIZE 32
  #endif
#endif


// PRIME_IF_X64: Выполнение then на архитектуре x86-64.
// Либо остальной части (если она задана) в противном случае.
#if PRIME_ARCH_X86 != 0 && PRIME_CODE_SIZE == 64
  #define PRIME_IF_X64(then, ...) then
#else
  #define PRIME_IF_X64(then, ...) __VA_ARGS__
#endif

// PRIME_IF_INTRIN: Выполнение then при разрешённых интринсиках архитектур x86/x86-64.
// Либо остальной части (если она задана) в противном случае.
#if defined(PRIME_ALLOW_INTRIN) && PRIME_ARCH_X86 != 0
  #include <intrin.h>
  #define PRIME_IF_INTRIN(then, ...) then
#else
  #define PRIME_IF_INTRIN(then, ...) __VA_ARGS__
#endif

// PRIME_IF_INTRIN_X64: Выполнение then при разрешённых интринсиках архитектуры x86-64.
// Либо остальной части (если она задана) в противном случае.
#if defined(PRIME_ALLOW_INTRIN) && PRIME_ARCH_X86 != 0 && PRIME_CODE_SIZE == 64
  #define PRIME_IF_INTRIN_X64(then, ...) then
#else
  #define PRIME_IF_INTRIN_X64(then, ...) __VA_ARGS__
#endif

// PRIME_IF_INT128: Выполнение then при поддержке типа __int128.
// Либо остальной части (если она задана) в противном случае.
#if defined(__SIZEOF_INT128__) && defined(PRIME_ALLOW_INT128)
  #define PRIME_IF_INT128(then, ...) then
#else
  #define PRIME_IF_INT128(then, ...) __VA_ARGS__
#endif

// PRIME_IF_INT128_OR_CODE64: Выполнение then при поддержке типа __int128 и/или на любых
// 64-битных архитектурах. Либо остальной части (если она задана) в противном случае.
#if defined(__SIZEOF_INT128__) && PRIME_CODE_SIZE == 64
  #define PRIME_IF_INT128_OR_CODE64(then, ...) then
#else
  #define PRIME_IF_INT128_OR_CODE64(then, ...) __VA_ARGS__
#endif


// PRIME_FAST_MOD: Вычисления остатка от деления n на mod.
// Метод выбирается в зависимости от возможностей компилятора/платформы. При наличии возможности
// быстрого деления (получения остатка от деления) используются переменные mod_r, mod_shift
// (обратная величина делителя и величина сдвига) и функция fast_mod64.
#define PRIME_FAST_MOD(n, mod, mod_r, mod_shift) PRIME_IF_INT128_OR_CODE64( fast_mod64(n, mod, mod_r, mod_shift) , n % mod )

////////////////////////////////////////////////////////////////////////////////////////////////////

// Вернуть номер старшего бита числа n (начиная с 0). Для n == 0 возвращает 0.
// Используется функцией is_prime_miller.
inline int bsr(uint32_t n)
{
  int size = 16, shift = 8;
  uint32_t mask = 0xFFFF;
  int result = 0;

  while (n > 1) {
    uint32_t low = n >> size;
    if (low) {
      n = low;
      result += size;
    } else {
      n &= mask;
    }
    size >>= 1;
    mask >>= shift;
    shift >>= 1;
  }

  return result;
}

// Вернуть старшие 64 бита результата умножения двух 64-битных чисел.
// Используется функцией fast_mod64.
inline uint64_t mul64_high(uint64_t a, uint64_t b)
{
  uint64_t x0, x1, x2, x3;
  const uint64_t al = uint32_t(a), bl = uint32_t(b), ah = a >> 32, bh = b >> 32;
  
  x3 = al * bl;  // low
  x0 = ah * bh;  // high
  x2 = ah * bl;  // mid2
  x1 = al * bh;  // mid1
  
  x2 += x3 >> 32;  // no carry, max = (2**32-1)**2 + 2**32-1
  x0 += x2 >> 32;
  x1 += uint32_t(x2);  // still no carry
  
  return x0 + (x1 >> 32);
  //mul64_low = (x1 << 32) | uint32_t(x3); - если это кому-то интересно :)
}

// Вернуть остаток от деления n на mod без деления (через умножение, используя mod_r и mod_shift).
// Используется функциями modpow и is_prime_miller.
inline uint64_t fast_mod64(uint64_t n, uint32_t mod, uint64_t mod_r, uint32_t mod_shift)
{
  PRIME_IF_INTRIN_X64(
    uint64_t quot;
    _umul128(n, mod_r, &quot);
  ,
    PRIME_IF_INT128(
      uint64_t quot = (unsigned __int128)(n) * mod_r >> 64;
    ,
      uint64_t quot = mul64_high(n, mod_r);
    )
  )
  quot >>= mod_shift;
  return n - quot * mod;  // n % mod;
}

// Возведение в степень по модулю: (base ** exp) % mod.
// В версии с параметрами mod_r и mod_shift выполняется деление через умножение
// с последующим вычислением остатка от деления.
// Используется функцией is_prime_miller.
PRIME_IF_INT128_OR_CODE64(
inline uint32_t modpow(uint32_t base, uint32_t exp, uint32_t mod, uint64_t mod_r, uint32_t mod_shift)
,
inline uint32_t modpow(uint32_t base, uint32_t exp, uint32_t mod)
)
{
  // При использовании функции только в тесте Миллера эту проверку можно опустить,
  // т.к. в тесте не используются такие значения
  //if (mod <= 1 /*|| exp < 0*/) { return 0; }

  uint64_t result = 1;
  uint64_t base_long = base;  // % mod делать необязательно
  bool moded = true;  // для уменьшения кол-ва вычислений остатков от делений

  // Оптимизированный классический алгоритм (с использованием флага moded и if в середине цикла)
  while (true) {
    if (exp & 1) {
      result *= base_long;
      moded = (result >= 0x100000000);
      if (moded) {
        result = PRIME_FAST_MOD(result, mod, mod_r, mod_shift);  // result %= mod
      }
    }
    exp >>= 1;
    if (exp <= 0) { break; }
    base_long *= base_long;
    if (base_long >= 0x100000000) {
      base_long = PRIME_FAST_MOD(base_long, mod, mod_r, mod_shift);  // base_long %= mod
    }
  }

  return uint32_t(moded ? result : PRIME_FAST_MOD(result, mod, mod_r, mod_shift));  // : result %= mod
}

// Тест Миллера (детерминированный тест простоты числа n).
// Проверен на всём диапазоне uint32_t.
// https://en.wikipedia.org/wiki/Miller–Rabin_primality_test#Deterministic_variants
// https://ru.wikipedia.org/wiki/Тест_Миллера_(теория_чисел)
bool is_prime_miller(uint32_t n)
{
  // Проверяем, чётность, что число больше 1, а также частный случай простого числа (2)
  if ((n & 1) == 0 || n <= 1) {
    return (n == 2);  // число 2 простое, остальные чётные числа составные
  }

  // Находим кол-во простых оснований, по которым необходимо сделать проверку
  static const uint32_t primes[5] = { 2, 3, 5, 7, 11 };
  static const uint32_t max_n[4] = { 2047, 1373653, 25326001, 3215031751 };
  int last_index = 0;
  while (last_index < 4 && n >= max_n[last_index]) { ++last_index; }

  // Находим d и r такие, что n = d * 2**r + 1
  uint32_t n_minus_1 = n - 1, d = n_minus_1;
  int r;
  for (r = 0; (d & 1) == 0; d >>= 1, ++r) { }

  // Для ускорения операций получения остатка от делений может использоваться деление через
  // умножение, для чего необходимо вычислить обратную величину делителя и величину сдвига.
  // Этот метод эффективен, т.к. в коде многократно вычисляется остаток от деления на одно
  // и то же число, и нам достаточно выполнить деление только в этом блоке кода. Поскольку
  // n всегда нечётное (и не может быть степенью двойки), переполнения при делении не будет.
  // Такой приём позволяет значительно ускорить работу функции (~ в 2 раза или даже более)!
PRIME_IF_INT128_OR_CODE64(
  uint32_t n_shift;
  PRIME_IF_INTRIN(
    unsigned long ln_shift;
    _BitScanReverse(&ln_shift, n);
    n_shift = ln_shift;
  ,
    n_shift = bsr(n);
  )
  PRIME_IF_INT128(
    uint64_t n_r = (((unsigned __int128)(1) << 64 << n_shift) + (n >> 1)) / n;
    uint64_t quot;
  ,  // 64-битный код
    std::lldiv_t div = std::lldiv(1llu << 32 << n_shift, n);
    uint64_t n_r = div.quot << 32;
    div = std::lldiv((div.rem << 32) + (n >> 1), n);
    n_r += div.quot;
  )
)

  // Если число n сильно псевдопростое по всем простым основаниям
  // до primes[last_index] включительно, то оно точно простое
  for (int index = 0; index <= last_index; ++index) {
    bool prime = false;
    uint64_t x = PRIME_IF_INT128_OR_CODE64( modpow(primes[index], d, n, n_r, n_shift) ,
                                            modpow(primes[index], d, n) );
    // Число n сильно псевдопростое по основанию p, если p ** (2 * d) mod n == 1...
    if (x == 1 || x == n_minus_1) { continue; }
    // ...или если p ** (2**j * d) mod n == n - 1 при j = 0..r-1 (но для 0 мы уже проверили выше)
    for (int i = 1; i < r; ++i) {
      x = PRIME_FAST_MOD(x * x, n, n_r, n_shift);  // x = x * x % n;
      if (x == n_minus_1) {
        prime = true;
        break;
      }
    }
    if (!prime) { return false; }  // число составное
  }

  return true;  // число простое
}

////////////////////////////////////////////////////////////////////////////////////////////////////

// Класс для замера времени выполнения кода
class TimeMeasure
{
  clock_t _start_time;
public:
  // Запустить таймер
  void start()
  {
    _start_time = clock();
  }
  // Остановить таймер и вывести время
  void stop_and_show()
  {
    double time = clock() - _start_time;
    time /= CLOCKS_PER_SEC;
    cout << "[Elapsed time: " << time << " seconds]" << "\n";
  }
};

////////////////////////////////////////////////////////////////////////////////////////////////////

int main()
{
  // Вывод первых 1000 чисел
  for (int n = 0; n < 1000; ++n) {
    if (is_prime_miller(n)) {
      cout << n << " ";
    }
  }
  cout << "\n\n";

  // Замер производительности
  TimeMeasure tm;
  tm.start();
  int count = 0;
  for (int n = 0; n <= 10000000; ++n) {
    if (is_prime_miller(n)) { ++count; }
  }
  tm.stop_and_show();
  cout << "count = " << count << "\n";

  return 0;
}

miller_test.py

""" Miller test, v1.00 (c) 2019 by Jin X (jin_x@list.ru) """

# Массив простых чисел (постепенно пополняется при использовании функции is_prime_miller)
prime_list = [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71 ]

# Для проверки числа, меньшего, чем ключ словаря miller_max_n, достаточно делать
# проверку по простым основаниям, не превышающим соответствующее этому ключу значение
miller_max_n = { 2047 : 2,
          1373653 : 3,
          25326001 : 5,
          3215031751 : 7,
          2152302898747 : 11,
          3474749660383 : 13,
          341550071728321 : 17,
          3825123056546413051 : 23,
          318665857834031151167461 : 37,
          3317044064679887385961981 : 41,
          1543267864443420616877677640751301 : 61 }

# Тест Миллера (детерминированный тест простоты числа n).
# Если fast = True, для больших чисел (более 10**36) выполняется ускоренный тест,
# с уменьшенным числом простых оснований, по которым ведётся проверка.
# Достоверность такого теста пока не доказана, но по меньшей мере, его можно
# использовать как более достоверную альтернативу теста Миллера-Рабина :)
# https://en.wikipedia.org/wiki/Miller–Rabin_primality_test#Deterministic_variants
# https://ru.wikipedia.org/wiki/Тест_Миллера_(теория_чисел)
def is_prime_miller(n, fast = False):

  from math import log, sqrt

  def is_prime_list(n):
    max_d = int(sqrt(n))
    for d in prime_list:
      if d > max_d:
        return True
      if not n % d:
        return False
    return True  # на всякий случай :)

  # Проверка числа n на составное по основанию p
  def composite(p, n, d, r):
    n_minus_1 = n - 1
    x = pow(p, d, n)
    # Число n сильно псевдопростое по основанию p, если p** (2 * d) mod n == 1...
    if x == 1 or x == n_minus_1:
      return False
    # ...или если p ** (2**i * d) mod n == n - 1 при i = 0..r-1 (но для 0 мы уже проверили выше)
    for i in range(1, r):
      x = x * x % n
      if x == n_minus_1:
        return False
    return True

  # Проверяем, чётность, что число больше 1, а также частный случай простого числа (2)
  if not n & 1 or n <= 1:
    return (n == 2)  # число 2 простое, остальные чётные числа составные

  # Находим кол-во простых оснований, по которым необходимо сделать проверку
  if n > 10**36:
    log_n = log(n)
    max_prime = log_n * log(log_n) / log(2) if fast else 2 * log_n * log_n
    max_prime = int(max_prime) - 1 | 1
  else:
    max_prime = 71
    for key in miller_max_n:
      if n < key:
        max_prime = miller_max_n[key]
        break

  # Находим d и r такие, что n = d * 2**r + 1
  d, r = n - 1, 0
  while not d & 1:
    d //= 2
    r += 1

  # Если число n сильно псевдопростое по всем простым основаниям
  # до max_prime включительно, то оно точно простое
  for prime in prime_list:
    if prime > max_prime:
      return True  # число простое
    if composite(prime, n, d, r):
      return False  # число составное

  # Добавляем недостающие простые числа в prime_list и проверяем по ним
  last_prime = prime_list[-1]
  if last_prime < max_prime:
    delta = 4 - last_prime % 6 // 2  # 2 или 4
    while last_prime < max_prime:
      last_prime += delta
      if is_prime_list(last_prime):
        prime_list.append(last_prime)
        if composite(prime, n, d, r):  # проверка числа
          return False  # число составное
      delta ^= 6  # 2 -> 4, 4 -> 2

  return True  # число простое

################################################################################

if __name__ == '__main__':
  
  from time import time
  
  def test(n, fast, correct, out_n):
    result_str = { True: 'Prime', False: 'Composite' }
    result_char = { True: '+', False: '-' }
    correct_str = { True: 'correct!', False: 'WRONG !!!' }
    result = is_prime_miller(n, fast)
    if out_n:
      print(f'{n} : {result_str[result]} ({correct_str[result==correct]})')
    else:
      print(result_char[result], end='')
  
  # Гигантское простое число Бэлла (взято из https://primes.utm.edu/primes/page.php?id=68825)
  giant = 93074010508593618333390737086787716686243541818180888374481714163260295083047249735772767217242983302809460000977534180207136267393493843565165244803160631468142310993331802362163022244196698256971348116848613873420484831145693123069522851464977918752008824682368616575925219950816540476437905078816850961411910271333677033484348376563240090830825402255766897769433849999097072825800009667464269609454294177226626947696038219790689936750529759962907392539903707559051710294385488230399560792411004667389541080494330079820301708455191991301691046952805571310530593988681849292364334866037961772856489592631353387119588912930029104717108771889838640129969152323291892894440849105218746915343681917125357252705291568270795558661856036155452894768239527100657674295540074496397200222788808730895367566036148544193391411249930638389221671262805683765652287931774606379446475019704511886714677938481928235170230529782702798198071512791277092287199903418895014552315221588641882843663507160334370905336400728991988818738787812880619370936400288087393153584463921744028247514969200055227428251657073540384943632840780163289976513840862197690165020628387737123529546265113690620682595841836568269584610937749225174156188272089587049752042320435999387410834519188271208577279962713103121177117736246578720020849615208138290410547459210639313481871280957180797480718884005158357183543875383933081640708958426808475868980596063805203682907241542158286516965706176501691352009055980316953619941361963586164642762338959226194401591549258894070494114321789561253423874645767485673563010694616837670389191021116993326818985640677682311168596513135927575792933795897024983955212555699886813758658727223213225641249054291854713271825236198768288749577317015750567399596468873488940152346191448776708760260862506238255653851554400298770502418390469037927740196990922130058457344538461597268140533944714325634938884545914139335512028740689585456916586292846456683229763623263845961927185120666686527368190661471902546889836939907242929408922820078908112593178663177685220525268101383971283991711468187276352738607911284318470208480309880183721371226886168592890172034997639285024092759687525920453573640538777106302852351553054823775813450680545320959747676102527895283952321119565456131914284468837228528467270883859016854414206071054324686179724452435704506155435210031925383788518515000655319634148229746788564381020601053143272002310437607878237640840123305186361012402650803349859965081202294515347182519213721197040915413263249473539740781608786907273923065127196445526332443113542957189094428043671781635432417130645135281343627517544700098433529452127971455501702330453614487341357174977756767117068294356318437149385243962447271471217433312351733571281240293461665450829761335596591586210391312306966997773285594528746729346018039023115838145677882687881461094746947227827301981448949646394994521319966023728976458814160934241105908961406982946398028919136619104849916909167562570774608666807683843343671704615600840389419697202833796957720397144242132316427467001808219485482456195736463596111707485715440237594384459161928415836077852237526665117484048997247449275300584446550437546119923676017959462712581196976718470946270331842562972612728361652280030892982127111700793139354703946990580256780069816918913085639153753984131390468635302755088889211367474268779633561838363953846659135906229513613392267266814066735012127403702413192430883159093465043866796901005656737875251268652602552279882927553368913466086109551491491194789740567321879833676289767448527915219283173310873247606366284048111931661775107155492303602877951956085944593035383609134038714354896277016656832069877486297785906138808032199478041298446201040604532790336475388815136237012786896885113359098836781297206766936652939896086173174247047512753636594129042150203504101570396300673678733698741829932012118685174241471375329706399365116190852969742529040156369682356942527779947968734604866975128678544682420679340574499610003971519618440516937488305570847146550169657044672599070091974948167777874966859847710614503223797831808567277923192280169244649411883773069412702823257335336231095595728578499725749559860832154398909223367972897328124149325331610872832467828260416056336532880372285600592198569202140923550067859152359527561994066586118059826706360048699302434574822753790260890021587852244744860631603541791133619526793181058747657744949691213833240225173150813910220165362694840945664849801045511812349221956400744150897341864900858146835458095842131465890241777553970152159303024640985017795077066167761624753595558912106613186934139850039921207148473490125950756237204950896105932772825720856228894276720028932601843556239266022962890064643933856257123158863877848475127602406307792588595211703986239432550726691852534321783665511020801887786921172059700879026067285947359724825015339666089008466190352993533122497861078203104758579483775098260468006288071847722122648015073671866043577071451595398504208603083753549564030349663312899219714349854025049280827748345322558570768577632728362536039844811568342727553594911765165120515649387783590754617053981056103148444168056157453359284719489933160382315541998163668080430170189604432196012778454138678438477287617931797344197371492016925294293920435571230125441856377635430992383031641860714124008747101118765506416821101339484788073020563327420998407651377506476214794114648110563115377659800929141531398790434303752231082196206662416811943482851130965329365467377939976152662541912142094277951396022236086513406015986868896776219944649030822428734842983667486613881443827330912807645197812117423741459096670510731379104448554954932379389781897741126092080837369956421197209552100624952564640377427157704473125717924200923183638927529737419861594670798574290635462323151713380155501912570516986807486808850943703276543637772748688396151956171036621635411599055947803963844239542311343491664913052539486422164914206103599970496248266628264533316702997562527874219219368930835654767520116901325628008414215779935053300527454338095585904167331992744976932786018872605231337251856220370815003679385877026959961549970646582412358020805216001234519254278426271698879069271009841356762362776355079842217812045389559364458319251495759620406596953099358508401125247456305868316439298039940486165859806267316803843769250751290909537174237439362050585335391392533062430185171682628621992342320221279940728553592588293913252900669869333819652288398187260320155987924359450740883499506500844712757333449557083885253703080601131
  
  # Простые числа (взяты из https://ru.wikipedia.org/wiki/Список_простых_чисел, повторяющиеся исключены)
  numbers = (
      2, 5, 877, 27644437, 35742549198872617291353508656626642567, 359334085968622831041960188598043661065388726959079837,  # числа Бэлла
      7, 47, 223, 3967, 16127, 1046527, 16769023, 1073676287, 68718952447, 274876858367, 4398042316799, 1125899839733759, 18014398241046527, 1298074214633706835075030044377087,  # простые числа Кэрола
      31, 127, 8191, 131071, 524287, 2147483647, 2305843009213693951, 618970019642690137449562111, 162259276829213363391578010288127, 170141183460469231731687303715884105727,  # простые числа Марсенна
      41, 239, 9369319, 63018038201, 489133282872437279, 19175002942688032928599, 123426017006182806728593424683999798008235734137469123231828679,  # простые числа Ньюмена-Шэнкса-Уильямса
      13, 89, 233, 1597, 28657, 514229, 433494437, 2971215073, 99194853094755497, 1066340417491710595814572169, 19134702400093278081449423917,  # простые числа Фибоначчи
      11, 37, 101, 9091, 9901, 333667, 909091, 99990001, 999999000001, 9999999900000001, 909090909090909091, 1111111111111111111, 11111111111111111111111, 900900900900990990990991,  # уникальные
      23, 719, 5039, 39916801, 479001599, 87178291199, 10888869450418352160768000001, 265252859812191058636308479999999, 263130836933693530167218012159999999, 8683317618811886495518194401279999999  # факториальные
    )

  # Большие простые числа (взяты из http://rosettacode.org/wiki/Miller–Rabin_primality_test)
  large = (
    35201546659608842026088328007565866231962578784643756647773109869245232364730066609837018108561065242031153677,
    10513733234846849736873637829838635104309714688896631127438692162131857778044158273164093838937083421380041997,
    24684249032065892333066123534168930441269525239006410135714283699648991959894332868446109170827166448301044689,
    76921421106760125285550929240903354966370431827792714920086011488103952094969175731459908117375995349245839343,
    32998886283809577512914881459957314326750856532546837122557861905096711274255431870995137546575954361422085081,
    30925729459015376480470394633122420870296099995740154747268426836472045107181955644451858184784072167623952123,
    14083359469338511572632447718747493405040362318205860500297736061630222431052998057250747900577940212317413063,
    10422980533212493227764163121880874101060283221003967986026457372472702684601194916229693974417851408689550781,
    36261430139487433507414165833468680972181038593593271409697364115931523786727274410257181186996611100786935727,
    15579763548573297857414066649875054392128789371879472432457450095645164702139048181789700140949438093329334293
  )

  # Огромные числа (простое только одно: 2**1279-1)
  huge = {
    643808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153: True,
    743808006803554439230129854961492699151386107534013432918073439524138264842370630061369715394739134090922937332590384720397133335969549256322620979036686633213903952966175107096769180017646161851573147596390153: False,
    2**1279-3: False, 2**1279-2: False, 2**1279-1: True, 2**1279+1: False
  }
  
  # Тесты скорости небольших чисел
  
  print('\nFast test:', end=' ')
  start = time()
  for n in numbers:
    test(n, True, True, False)
  elapsed1 = time()-start
  print(f'\nElapsed time: {elapsed1:.3f} sec')

  print('\nSlow test:', end=' ')
  start = time()
  for n in numbers:
    test(n, False, True, False)
  elapsed2 = time()-start
  print(f'\nElapsed time: {elapsed2:.3f} sec')
  print(f'Speed ratio = {elapsed2/elapsed1:.3f}')

  # Тесты скорости больших чисел

  print('\nFast test of large numbers:')
  start = time()
  for n in large:
    test(n - 2, True, False, True)
    test(n, True, True, True)
  elapsed1 = time()-start
  print(f'Elapsed time: {elapsed1:.3f} sec')
  
  print('\nSlow test of large numbers:')
  start = time()
  for n in large:
    test(n - 2, False, False, True)
    test(n, False, True, True)
  elapsed2 = time()-start
  print(f'Elapsed time: {elapsed2:.3f} sec')
  print(f'Speed ratio = {elapsed2/elapsed1:.3f}')

  # Тесты скорости огромных чисел

  print('\nFast test of HUGE numbers:')
  start = time()
  for n in huge:
    test(n, True, huge[n], True)
  elapsed = time()-start
  print(f'Elapsed time: {elapsed:.3f} sec')
  
  # Тесты скорости гигантского числа Бэлла

  print('\nFast test of G*I*A*N*T numbers:')
  start = time()
  test(giant-2, True, False, True)
  elapsed = time()-start
  print(f'Elapsed time: {elapsed:.3f} sec')

  ''' Слишком медленно, чтобы запускать в работу :)
  start = time()
  test(giant, True, True, True)
  elapsed = time()-start
  print(f'Elapsed time: {elapsed:.3f} sec')
  '''

  print(f'\nPrimes in prime_list: {len(prime_list)}')

miller_test_v1.00.cpp

// Miller test, v1.00 (c) 2019 by Jin X (jin_x@list.ru)

#include <iostream>
#include <ctime>

using std::cout;

////////////////////////////////////////////////////////////////////////////////////////////////////

// Возведение в степень по модулю: (base ** exp) % mod.
// Используется функцией is_prime_miller.
unsigned int modpow(unsigned int base, unsigned int exp, unsigned int mod)
{
  // При использовании функции только в тесте Миллера эту проверку можно опустить,
  // т.к. в тесте не используются такие значения
  if (mod == 1 /*|| exp < 0*/) { return 0; }

  unsigned long long result = 1;
  unsigned long long base_long = base;  // % mod
  bool moded = true;  // для уменьшения кол-ва вычислений остатков от делений

  // Оптимизированный классический алгоритм (с использованием флага moded и if в середине цикла)
  while (true) {
    if (exp & 1) {
      result *= base_long;
      moded = (result >= 0x100000000LL);
      if (moded) { result %= mod; }
    }
    exp >>= 1;
    if (exp <= 0) { break; }
    base_long *= base_long;
    if (base_long >= 0x100000000LL) {
      base_long %= mod;
    }
  }

  return moded ? result : result % mod;
}

// Тест Миллера (детерминированный тест простоты числа n).
// Проверен на всём диапазоне unsigned int.
// https://en.wikipedia.org/wiki/Miller–Rabin_primality_test#Deterministic_variants
// https://ru.wikipedia.org/wiki/Тест_Миллера_(теория_чисел)
bool is_prime_miller(unsigned int n)
{
  // Проверяем, чётность, что число больше 1, а также частный случай простого числа (2)
  if ((n & 1) == 0 || n <= 1) {
    return (n == 2);  // число 2 простое, остальные чётные числа составные
  }

  // Находим кол-во простых оснований, по которым необходимо сделать проверку
  static const unsigned int primes[5] = { 2, 3, 5, 7, 11 };
  static const unsigned int max_n[4] = { 2047, 1'373'653, 25'326'001, 3'215'031'751 };
  int last_index = 0;
  while (last_index < 4 && n >= max_n[last_index]) { ++last_index; }

  // Находим d и r такие, что n = d * 2**r + 1
  unsigned int r, d = n - 1;
  for (r = 0; (d & 1) == 0; d >>= 1, ++r) { }

  // Если число n сильно псевдопростое по всем простым основаниям
  // до primes[last_index] включительно, то оно точно простое
  unsigned int n_minus_1 = n - 1;
  for (int i = 0; i <= last_index; ++i) {
    bool prime = false;
    unsigned long long x = modpow(primes[i], d, n);
    // Число сильно псевдопростое по основанию p, если p ** (2 * d) mod n == 1...
    if (x == 1 || x == n_minus_1) { continue; }
    // ...или если p ** (2**j * d) mod n == n - 1 при j = 0..r-1 (но для 0 мы уже проверили выше)
    for (int j = 1; j < r; ++j) {
      x = x * x % n;
      if (x == n_minus_1) {
        prime = true;
        break;
      }
    }
    if (!prime) { return false; }  // число составное
  }

  return true;  // число простое
}

////////////////////////////////////////////////////////////////////////////////////////////////////

// Класс для замера времени выполнения кода
class TimeMeasure
{
  clock_t _start_time;
public:
  // Запустить таймер
  void start()
  {
    _start_time = clock();
  }
  // Остановить таймер и вывести время
  void stop_and_show()
  {
    double time = clock() - _start_time;
    time /= CLOCKS_PER_SEC;
    cout << "[Elapsed time: " << time << " seconds]" << "\n";
  }
};

////////////////////////////////////////////////////////////////////////////////////////////////////

int main()
{
  // Вывод первых 1'000 чисел
  for (int n = 0; n < 1'000; ++n) {
    if (is_prime_miller(n)) {
      cout << n << " ";
    }
  }
  cout << "\n\n";

  // Замер производительности
  TimeMeasure tm;
  tm.start();
  int count = 0;
  for (int n = 0; n <= 10'000'000; ++n) {
    if (is_prime_miller(n)) { ++count; }
  }
  tm.stop_and_show();
  cout << "count = " << count << "\n";

  return 0;
}