Способ найти ближайшее простое число к длинному целому без знака (32 бита в ширину) в C?

ОБНОВЛЕНИЕ 2: исправлены (грубым образом) некоторые ошибки, приводившие к неправильным ответам для малых n. Спасибо Бретту Хейлу за то, что заметил! Также добавлены некоторые утверждения для документирования некоторых предположений.

ОБНОВЛЕНИЕ: я закодировал это, и это кажется достаточно быстрым для ваших требований (решил 1000 случайных экземпляров из [2 ^ 29, 2 ^ 32-1] за ‹100 мс, на машине с частотой 2,2 ГГц - не строгий тест, но тем не менее убедительный).

Он написан на C++, так как именно в нем был мой ситовой код (из которого я адаптировался), но преобразование в C должно быть простым. Использование памяти также (относительно) невелико, что можно увидеть при осмотре.

Вы можете видеть, что из-за того, как вызывается функция, возвращаемое число является ближайшим простым числом, которое умещается в 32 бита, но на самом деле это одно и то же, поскольку простые числа около 2^32 равны 4294967291 и 4294967311.

Я пытался убедиться, что не будет ошибок из-за целочисленного переполнения (поскольку мы имеем дело с числами вплоть до UINT_MAX); надеюсь, я не ошибся там. Код можно было бы упростить, если бы вы хотели использовать 64-битные типы (или знали, что ваши числа будут меньше, чем 2^32-256), поскольку вам не придется беспокоиться о переносе в условиях цикла. Также эта идея масштабируется для больших чисел, если вы готовы вычислять/сохранять маленькие простые числа до необходимого предела.

Я также должен отметить, что сито малых простых чисел работает довольно быстро для этих чисел (4-5 мс по грубому измерению), поэтому, если вам особенно не хватает памяти, выполнимо запускать его каждый раз вместо сохранения малых простых чисел (вы в этом случае, вероятно, хотелось бы сделать массивы mark[] более эффективными с точки зрения пространства)

#include <iostream>
#include <cmath>
#include <climits>
#include <cassert>

using namespace std;

typedef unsigned int UI;

const UI MAX_SM_PRIME = 1 << 16;
const UI MAX_N_SM_PRIMES = 7000;
const UI WINDOW = 256;

void getSMPrimes(UI primes[]) {
  UI pos = 0;
  primes[pos++] = 2;

  bool mark[MAX_SM_PRIME / 2] = {false};
  UI V_SM_LIM = UI(sqrt(MAX_SM_PRIME / 2));
  for (UI i = 0, p = 3; i < MAX_SM_PRIME / 2; ++i, p += 2)
    if (!mark[i]) {
      primes[pos++] = p;
      if (i < V_SM_LIM)
        for (UI j = p*i + p + i; j < MAX_SM_PRIME/2; j += p)
          mark[j] = true;
      }
  }

UI primeNear(UI n, UI min, UI max, const UI primes[]) {
  bool mark[2*WINDOW + 1] = {false};

  if (min == 0) mark[0] = true;
  if (min <= 1) mark[1-min] = true;

  assert(min <= n);
  assert(n <= max);
  assert(max-min <= 2*WINDOW);

  UI maxP = UI(sqrt(max));
  for (int i = 0; primes[i] <= maxP; ++i) {
    UI p = primes[i], k = min / p;
    if (k < p) k = p;
    UI mult = p*k;
    if (min <= mult)
      mark[mult-min] = true;
    while (mult <= max-p) {
      mult += p;
      mark[mult-min] = true;
      }
    }

  for (UI s = 0; (s <= n-min) || (s <= max-n); ++s)
    if ((s <= n-min) && !mark[n-s-min])
      return n-s;
    else if ((s <= max-n) && !mark[n+s-min])
      return n+s;

  return 0;
  }

int main() {
  UI primes[MAX_N_SM_PRIMES];
  getSMPrimes(primes);

  UI n;
  while (cin >> n) {
    UI win_min = (n >= WINDOW) ? (n-WINDOW) : 0;
    UI win_max = (n <= UINT_MAX-WINDOW) ? (n+WINDOW) : UINT_MAX;

    if (!win_min)
      win_max = 2*WINDOW;
    else if (win_max == UINT_MAX)
      win_min = win_max-2*WINDOW;

    UI p = primeNear(n, win_min, win_max, primes);
    cout << "found nearby prime " << p << " from window " << win_min << ' ' << win_max << '\n';
    }
  }

Вы можете просеять интервалы в этом диапазоне, если знаете простые числа до 2^16 (существует только 6542 ‹= 2^16; вам следует подняться немного выше, если само простое число может быть больше 2^32 - 1). Не обязательно самый быстрый способ, но очень простой, и более сложные методы первичного тестирования действительно подходят для гораздо больших диапазонов.

По сути, сделайте обычное сито Эратосфена, чтобы получить «маленькие» простые числа (скажем, первые 7000). Очевидно, вам нужно сделать это только один раз в начале программы, но это должно быть очень быстро.

Затем, предположим, что ваше «целевое» число равно «a», рассмотрим интервал [a-n/2, a+n/2) для некоторого значения n. Вероятно, n = 128 — разумное место для начала; вам может потребоваться попробовать соседние интервалы, если числа в первом составные.

Для каждого «маленького» простого числа p вычеркните его кратные в диапазоне, используя деление, чтобы найти, с чего начать. Одна оптимизация заключается в том, что вам нужно только начать вычеркивать кратные числа, начиная с p*p (это означает, что вы можете перестать рассматривать простые числа, как только p*p превысит интервал).

Большинство простых чисел, за исключением первых нескольких, будут иметь внутри интервала либо единицу, либо ноль, кратные; чтобы воспользоваться этим, вы можете предварительно игнорировать несколько первых простых чисел. Проще всего игнорировать все четные числа, но нередко игнорируются числа, кратные 2, 3 и 5; это оставляет целые числа, конгруэнтные 1, 7, 11, 13, 17, 19, 23 и 29 по модулю 30 (их восемь, которые хорошо сопоставляются с битами байта при просеивании большого диапазона).

...Там как-то пошло не так; в любом случае, как только вы обработали все маленькие простые числа (до p*p > a+n/2), вы просто ищете в интервале числа, которые вы не вычеркнули; так как вы хотите ближе всего начать искать там и искать наружу в обоих направлениях.

наибольшая разница между простыми числами в диапазоне до (2^32 - 1) равна (335 ). Существует (6542) простых чисел меньше (2 ^ 16), которые можно свести в таблицу и использовать для просеивания последовательных нечетных значений после одноразовой настройки. Очевидно, что только простые числа ‹= floor(sqrt(candidate)) должны быть проверены на наличие конкретного значения-кандидата.

В качестве альтернативы: детерминированный вариант Миллера-Рабина test с основаниями SPRP: {2, 7, 61} достаточно, чтобы доказать простоту 32-битного значения. Из-за сложности теста (требует возведения в степень и т. д.) я сомневаюсь, что он будет таким же быстрым для таких маленьких кандидатов.

Редактировать: На самом деле, если умножение/уменьшение можно сохранить до 32-битного возведения в степень (может потребоваться 64-битная поддержка), тест M-R может быть лучше. Первичные зазоры, как правило, будут намного меньше, что делает затраты на настройку сита чрезмерными. Без больших таблиц поиска и т. д. вы также можете получить преимущество за счет лучшей локализации кэша.

Кроме того: произведение простых чисел {2, 3, 5, 7, 11, 13, 17, 19, 23} = (223092870). Явно протестируйте любого кандидата в [2, 23]. Вычислить наибольший общий делитель: g = gcd(u, 223092870UL). Если (g != 1), кандидат составной. Если (g == 1 && u < (29 * 29)), кандидат (u > 23) определенно является простым. В противном случае переходите к более дорогим тестам. Одиночный тест gcd с использованием 32-битной арифметики очень дешев, и, согласно теореме Мертенса (?), он обнаружит ~ 68,4% всех нечетных составных чисел.