Проблема заключается в следующем.

Ввод: целое числоn

Выход: Наименьшее простое число больше, чем n.

Задача состоит в том, чтобы дать самый быстрый код для этого. Я протестирую код на значениях, начиная примерно с10^8 размера 10^200и удваивая его, пока это не займет более одной минуты на моем компьютере 10 секунд.

Код победы найдет следующее простое число для наибольшего размера ввода.

Для сравнения, простое сито, написанное на python, может найти следующее простое число больше 10^8 примерно за 20секунды.

Требование, чтобы я мог проверить это на моем компьютере Ubuntu 4 ГБ RAM, строгое. Весь код должен быть свободным (в обоих смыслах), и если он использует библиотеки, он также должен быть свободным и легко устанавливаемым. Любые сообщения о ложных простых числах немедленно дисквалифицируют представление.

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

Стол пока

• Python. Удивительная 357цифра простого числа 343239883006530485749095039954069660863471765007165270469723172959277159169882802606127982033072727748864815569574042901856099399985832190628701414555752857600000000000000000000000000000000000000002872284792758930912601189043411951050852357613658978971208596097634095500808832510259693761982135208603287199546795000697807728609476163156438356035166156820611была конечным числом менее 10 секунд с использованием кода, предоставленного primo. Кто-нибудь победит эту первую запись?

Python ~ 451 цифр

Это часть библиотеки, которую я написал для задачи о полупростой факторизации , с удалением ненужных функций. Он использует тест примитивности Baillie-PSW , который технически является вероятностным тестом, но на сегодняшний день нет известных псевдоприемников - и даже есть денежное вознаграждение, если вы можете его найти (или за предоставление доказательства того, что его не существует) ,

Изменить : я не понял, что Python имеет встроенное модульное возведение в степень. Замена собственной на встроенные приводит к увеличению производительности примерно на 33%.

my_math.py

# legendre symbol (a|m)
# note: returns m-1 if a is a non-residue, instead of -1
def legendre(a, m):
  return pow(a, (m-1) >> 1, m)

# strong probable prime
def is_sprp(n, b=2):
  d = n-1
  s = 0
  while d&1 == 0:
    s += 1
    d >>= 1

  x = pow(b, d, n)
  if x == 1 or x == n-1:
    return True

  for r in range(1, s):
    x = (x * x)%n
    if x == 1:
      return False
    elif x == n-1:
      return True

  return False

# lucas probable prime
# assumes D = 1 (mod 4), (D|n) = -1
def is_lucas_prp(n, D):
  P = 1
  Q = (1-D) >> 2

  # n+1 = 2**r*s where s is odd
  s = n+1
  r = 0
  while s&1 == 0:
    r += 1
    s >>= 1

  # calculate the bit reversal of (odd) s
  # e.g. 19 (10011) <=> 25 (11001)
  t = 0
  while s > 0:
    if s&1:
      t += 1
      s -= 1
    else:
      t <<= 1
      s >>= 1

  # use the same bit reversal process to calculate the sth Lucas number
  # keep track of q = Q**n as we go
  U = 0
  V = 2
  q = 1
  # mod_inv(2, n)
  inv_2 = (n+1) >> 1
  while t > 0:
    if t&1 == 1:
      # U, V of n+1
      U, V = ((U + V) * inv_2)%n, ((D*U + V) * inv_2)%n
      q = (q * Q)%n
      t -= 1
    else:
      # U, V of n*2
      U, V = (U * V)%n, (V * V - 2 * q)%n
      q = (q * q)%n
      t >>= 1

  # double s until we have the 2**r*sth Lucas number
  while r > 0:
      U, V = (U * V)%n, (V * V - 2 * q)%n
      q = (q * q)%n
      r -= 1

  # primality check
  # if n is prime, n divides the n+1st Lucas number, given the assumptions
  return U == 0

# primes less than 212
small_primes = set([
    2,  3,  5,  7, 11, 13, 17, 19, 23, 29,
   31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
   73, 79, 83, 89, 97,101,103,107,109,113,
  127,131,137,139,149,151,157,163,167,173,
  179,181,191,193,197,199,211])

# pre-calced sieve of eratosthenes for n = 2, 3, 5, 7
indices = [
    1, 11, 13, 17, 19, 23, 29, 31, 37, 41,
   43, 47, 53, 59, 61, 67, 71, 73, 79, 83,
   89, 97,101,103,107,109,113,121,127,131,
  137,139,143,149,151,157,163,167,169,173,
  179,181,187,191,193,197,199,209]

# distances between sieve values
offsets = [
  10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6,
   6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4,
   2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6,
   4, 2, 4, 6, 2, 6, 4, 2, 4, 2,10, 2]

max_int = 2147483647

# an 'almost certain' primality check
def is_prime(n):
  if n < 212:
    return n in small_primes

  for p in small_primes:
    if n%p == 0:
      return False

  # if n is a 32-bit integer, perform full trial division
  if n <= max_int:
    i = 211
    while i*i < n:
      for o in offsets:
        i += o
        if n%i == 0:
          return False
    return True

  # Baillie-PSW
  # this is technically a probabalistic test, but there are no known pseudoprimes
  if not is_sprp(n): return False
  a = 5
  s = 2
  while legendre(a, n) != n-1:
    s = -s
    a = s-a
  return is_lucas_prp(n, a)

# next prime strictly larger than n
def next_prime(n):
  if n < 2:
    return 2
  # first odd larger than n
  n = (n + 1) | 1
  if n < 212:
    while True:
      if n in small_primes:
        return n
      n += 2

  # find our position in the sieve rotation via binary search
  x = int(n%210)
  s = 0
  e = 47
  m = 24
  while m != e:
    if indices[m] < x:
      s = m
      m = (s + e + 1) >> 1
    else:
      e = m
      m = (s + e) >> 1

  i = int(n + (indices[m] - x))
  # adjust offsets
  offs = offsets[m:]+offsets[:m]
  while True:
    for o in offs:
      if is_prime(i):
        return i
      i += o

Пример тестового скрипта:

from time import clock
from my_math import *

n = i = 317**79
while True:
  i *= 317
  time1 = clock()
  n, o = next_prime(i), n
  span = clock()-time1
  if span > 10:
    break
  print(len(str(n)), span)
print(o)

Коэффициент 317 был выбран, потому что это примерно квадратный корень 10000 , прибавляя примерно 2,5 цифры на одну итерацию (и потому что удвоение было слишком медленным, чтобы пройти). Выходные данные показывают текущее количество цифр и время.

Пример результатов:

201 0.13121248650317288
203 0.059535499623555505
206 0.9157767258129175
208 0.2583420518529589
211 0.15367400046653978
213 0.32343915218274955
216 1.3962866788935466
218 0.5986165839513125
221 0.973842206202185
223 2.346910291671148
...
428 0.932809896229827
431 4.345940056627313
433 9.511724255457068
436 6.089835998709333
438 1.3793498894412721
441 4.290633027381972
443 3.5102506044762833
446 3.1629148397352083
448 3.364759208223404
451 7.34668009481652
1551197868099891386459896063244381932060770425565921999885096817830297496627504652115239001983985153119775350914638552307445919773021758654815641382344720913548160379485681746575245251059529720935264144339378936233043585239478807971817857394193701584822359805681429741446927344534491412763713568490429195862973508863067230162660278070962484418979417980291904500349345162151774412157280412235743457342694749679453616265540134456421369622519723266737913

Весь код теперь совместим с Python 3.

C ++ с GMP: 567 цифр

Использует реализацию Миллера-Рабина в GMP. Это может вернуть ложноположительный результат, но удача на самом деле ударит его с вероятностью 2 ^ -200.

#include <gmp.h>
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>

double time() {
  struct timeval t;
  gettimeofday(&t, NULL);
  return t.tv_usec  * 1e-6 + t.tv_sec;
}

int main(int argc, char *argv[]) {
  mpz_t n, m;
  mpz_init_set_ui(n, 10);
  mpz_pow_ui(n, n, 200);
  mpz_init(m);
  for (int i = 0; true; i++, mpz_mul_ui(n, n, 2)) {
    double start = time();
    for (mpz_add_ui(m, n, 1); !mpz_millerrabin(m, 100); mpz_add_ui(m, m, 2)) ;
    double t = time() - start;
    gmp_printf("%d %Zd %f\n", i, m, t);
    if (t > 10.0) break;
  }
}

Находит простое число 10^200 * 2^1216 + 361(567 цифр) перед тем, как работать на моем медленном ноутбуке.

Perl с модулем GMP, 1300 цифр

Используя мой модуль Math :: Prime :: Util и его GMP-сервер . Точная точка пересечения будет зависеть от вашего компьютера и наличия у вас последней библиотеки GMP. Весь код бесплатный (модули на github и CPAN, а GMP находится в свободном доступе). Я запускал их на Ubuntu AWS, а также на настольной Ubuntu (и Fedora, и AIX, и NetBSD и т. Д.).

Код ядра находится в C и C + GMP. next_prime из MPU видит, что число слишком велико, и перенаправляет его на серверную часть GMP (или чистый код Perl, если серверная часть не установлена). Это преобразует строку в формат mpz и преобразует его обратно во входной тип объекта или Math :: BigInt. Сам по себе next_prime делает:

• мод 30 колесо
• отслеживает оставшийся мод 23 #, так что он может делать собственные модули для простых чисел до 23
• вероятный главный тест на вещах, которые проходят эти.

Возможное простое испытание:

• проверять крошечные делители, используя mpz_gcd_ui (в 64-битных двух из них проверяется до 101)
• проверять на наличие малых делителей, используя большие примитивы, вычисленные по одному. Это либо простые числа до 10 КБ, либо 40 КБ в зависимости от размера ввода.
• для значений, превышающих 2 ^ 1600, выполняется дальнейшее пробное деление с использованием сетки деревьев. Это может быть сделано более эффективно.
• наконец, проводится ES BPSW (тест Миллера-Рабина с основанием 2, за которым следует очень сильный тест Лукаса ).

Все, что до ES BPSW - это просто оптимизация, которую мы, конечно, хотим использовать для next_prime. next_prime также реализован в Perl с использованием модуля Math :: BigInt (в ядре с необязательным внутренним интерфейсом Pari и GMP). Это делает AES BPSW (как Pari), но не так оптимизирован.

Я подумал о достоинствах версии на основе частичного сита, используя диапазон, например, 2 достоинств. Я просто не уверен, будет ли это действительно лучше, так как большую часть времени мы будем делать ненужное просеивание, поскольку разрыв будет небольшим, а иногда для большого разрыва нам придется повторять это несколько раз.

Библиотека реализует ECPP (включая сертификаты), чтобы мы могли выполнить проверку результата, но 1200 цифр действительно слишком велико для крошечного набора по умолчанию включенных многочленов (есть метод для загрузки больших наборов - доказательства могут занять немного меньше времени) 15 минут, что немного быстрее, чем APR-CL Пари, но немного медленнее, чем mpz_aprcl от WraithX). Одним из недостатков ECPP по сравнению с APR-CL является то, что он имеет большую временную дисперсию, поэтому есть вероятность, что он превысит 10 секунд для некоторого числа до того, как будет достигнуто среднее время. С доказательством я думаю, что мы ограничены чем-то в диапазоне 400 цифр, если мы не разрешаем многопоточное программное обеспечение.

#!/usr/bin/env perl
use warnings;
use strict;
use Math::Prime::Util ":all";
use Math::Prime::Util::GMP;  # Barf if the backend isn't installed
use Time::HiRes qw(gettimeofday tv_interval);
use Math::GMP;

my $n = Math::GMP->new(10) ** 200;
while (1) {
  my $start = [gettimeofday];
  my $np = next_prime($n);
  my $sec = tv_interval($start);
  my $len = length($n);
  die "next_prime $len = +",$np-$n," in $sec seconds\n" if $sec > 10;
  warn "  next_prime $len = +",$np-$n," in $sec seconds\n";
  $n *= 10;
}

Я решил попробовать с той же последовательностью, которую использовал primo. Он достиг 1191 цифры, так как здесь мы пробили 18138. Я также протестировал код primo, используя последний my_math.py. Он достигает 630 цифр с последовательностью 10 ^ е и 641 с его последовательностью. Очень впечатляет компактный полностью Python-код без множества предварительных тестов.