Форум программистов, компьютерный форум, киберфорум
Наши страницы
Jin X
Войти
Регистрация
Восстановить пароль
Рейтинг: 5.00. Голосов: 4.

Проверка числа на простоту – ускоряемся!

Запись от Jin X размещена 28.02.2019 в 01:29
Обновил(-а) Jin X 05.03.2019 в 21:24

Проверка числа на простоту – ускоряемся!

Самый простой способ проверки числа N на простоту – проверить его делимость на все числа от 2 до корня из N.

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
// Тест простоты числа n
bool is_prime(unsigned int n)
{
  // Проверяем, что число больше единицы (простое всегда больше)
  if (n <= 1) { return false; }
 
  // Проверяем делимость на числа от 2 до корня из n:
  for (unsigned int div = 2, max_div = sqrt(n); div <= max_div; ++div) {
    if (n % div == 0) { return false; }
  }
 
  // Число простое
  return true;
}
Но этот код слишком медленный. Первое, что приходит на ум для ускорения – исключить из списка делителей чётные числа, сократив таким образом количество итераций вдвое (по сути же, нам важно проверить делимость только на другие простые числа, а простые числа не бывают чётными, кроме 2).

Проверим делимость числа на 2, а затем на все нечётные, начиная с 3:
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
bool is_prime2(unsigned int n)
{
  // Проверяем, что число больше 1, а также частный случай простого числа (2)
  if (n <= 2) {
    return (n == 2);
  }
 
  // Проверяем делимость на 2
  if ((n & 1) == 0) { return false; }
 
  // Проверяем делимость на нечётные числа, начиная с 3:
  // 3, 5, 7, 9, 11, 13, 15, 17, 19, 21...
  for (unsigned int div = 3, max_div = sqrt(n);
       div <= max_div;
       div += 2) {
    if (n % div == 0) { return false; }
  }
 
  // Число простое
  return true;
}
Надеюсь, никого не смущает запись n & 1 вместо n % 2? Да, конечно, умный компилятор сам заменит % на &, но я хочу указать это явно. Если вас это напрягает, замените

Этот код тоже далёк от совершенства, поскольку мы можем исключить числа, делящиеся не только на 2, но и на 3, сократив кол-во итераций ещё в 1.5 раза (а в сравнении с первым вариантом – в 3 раза). Чтобы это сделать, нам нужно увеличивать делитель то на 2, то на 4. Для этого будем проделывать с "дельтой" операцию xor 6. Проверим: 2 xor 6 = 4, 4 xor 6 = 2. Всё до гениальности просто! Цикл начнём с 5.

Let's do it:
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
// Тест простоты числа n
bool is_prime3(unsigned int n)
{
  // Проверяем, что число больше 1, а также частные случаи простых чисел (2 и 3)
  if (n <= 3) {
    return (n >= 2);
  }
 
  // Проверяем делимость на 2 и 3
  // p.s. При оптимизации нормальными компиляторами с опцией -O2 деление (взятие остатка от деления)
  // на константу будет выполняться быстрее, чем на переменную с заранее неизвестным значением
  if ((n & 1) == 0 || n % 3 == 0) { return false; }
 
  // Проверяем делимость на числа, не делящиеся на 2 и 3, начиная с 5
  // (с приращением то на 2, то на 4): 5, 7, 11, 13, 17, 19, 23, 25, 29, 31...
  for (unsigned int div = 5, delta = 2, max_div = sqrt(n);
       div <= max_div;
       div += delta, delta ^= 6) {  // delta ^= 6 меняет значение delta с 2 на 4 и наоборот
    if (n % div == 0) { return false; }
  }
 
  // Число простое
  return true;
}
На этом можно было бы остановиться (и для большинства случаев такого варианта будет достаточно). Но мы всё же продолжим издеваться над математикой!

Исключим числа, делящиеся на 5, получив прирост ещё примерно на 25% (или в 3.75 раза в сравнении с первым вариантом). Но здесь нам придётся использовать массив значений delta, т.к. двумя числами (как в предыдущем случае – 2 и 4) мы уже не обойдёмся. Нам понадобится 8 чисел: { 4, 2, 4, 2, 4, 6, 2, 6 }, которые будут повторяться. И начинать цикл мы будет с 7.

Барабанная дробь
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
// Тест простоты числа n
bool is_prime5(unsigned int n)
{
  // Проверяем, что число больше 1, а также частные случаи простых чисел (2, 3 и 5)
  if (n <= 5) {
    return (n >= 2 && n != 4);
  }
 
  // Проверяем делимость на 2, 3 и 5
  // p.s. При оптимизации нормальными компиляторами с опцией -O2 деление (взятие остатка от деления)
  // на константу будет выполняться быстрее, чем на переменную с заранее неизвестным значением
  if ((n & 1) == 0 || n % 3 == 0 || n % 5 == 0) { return false; }
 
  // Проверяем делимость на остальные числа, не делящиеся на 2, 3 и 5, начиная с 7
  // (с приращением, взятым из массива delta): 7, 11, 13, 17, 19, 23, 29, 31, 37...
  static const unsigned char delta[8] = { 4, 2, 4, 2, 4, 6, 2, 6 };  // массив приращений
  for (unsigned int div = 7, idx = 0, max_div = sqrt(n);  // idx - индекс в массиве delta
    div <= max_div;
    div += delta[idx], idx = idx+1 & 7) {  // увеличиваем индекс на 1 с остатком от деления на 8
    if (n % div == 0) { return false; }
  }
 
  // Число простое
  return true;
}
Да, здесь, как и раньше, нам придётся проверять делимость на 2, 3 и 5 вручную. Но это даже к лучшему: при оптимизации нормальными компиляторами с опцией -O2 деление (или взятие остатка от деления) на константу будет выполняться быстрее, чем на переменную с заранее неизвестным значением (как в цикле).

Кроме того, я специально использовал проверку вида ...(n % 3 == 0 && n != 3)..., а не вынес if (n == 2 || n == 3 || n == 5) { return true; } в начало, т.к. такой вариант должен работать быстрее (правда, тесты хоть и показывают прирост, но мизерный).

Хотите ещё?

Можно, это даст ещё приращение на 10-15%, но для исключения чисел, делящихся на 7 массив delta будет состоять уже из 48 чисел (цикл начинается с 11):
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
// Тест простоты числа n
bool is_prime7(unsigned int n)
{
  // Проверяем, что число больше 1, а также частные случаи простых чисел (2, 3, 5 и 7)
  if (n <= 10) {
    return (n == 2 || n == 3 || n == 5 || n == 7);
  }
 
  // Проверяем делимость на 2, 3, 5 и 7
  // p.s. При оптимизации нормальными компиляторами с опцией -O2 деление (взятие остатка от деления)
  // на константу будет выполняться быстрее, чем на переменную с заранее неизвестным значением
  if ((n & 1) == 0 || n % 3 == 0 || n % 5 == 0 || n % 7 == 0) { return false; }
 
  // Проверяем делимость на остальные числа, не делящиеся на 2, 3, 5 и 7, начиная с 11
  // (с приращением, взятым из массива delta): 11, 13, 17, 19, 23, 29, 31, 37, 41, 43...
  static const unsigned char delta[48] = {
    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, 10 };  // массив приращений
 
  for (unsigned int div = 11, idx = -1, max_div = sqrt(n);  // idx - индекс в массиве delta
    div <= max_div;
    div += delta[idx], idx = (idx == 47 ? 0 : idx + 1)) {  // увеличиваем индекс с остатком от деления на 48
    if (n % div == 0) { return false; }
  }
 
  // Число простое
  return true;
}

Для исключения чисел, делящихся ещё и на 11, массив delta будет состоять уже из 480 чисел (цикл начинается с 13).
Кликните здесь для просмотра всего текста
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
// Тест простоты числа n
bool is_prime11(unsigned int n)
{
  // Проверяем, что число больше 1, а также частные случаи простых чисел (2, 3, 5, 7 и 11)
  if (n <= 12) {
    return (n == 2 || n == 3 || n == 5 || n == 7 || n == 11);
  }
 
  // Проверяем делимость на 2, 3, 5, 7 и 11
  // p.s. При оптимизации нормальными компиляторами с опцией -O2 деление (взятие остатка от деления)
  // на константу будет выполняться быстрее, чем на переменную с заранее неизвестным значением
  if ((n & 1) == 0 || n % 3 == 0 || n % 5 == 0 || n % 7 == 0 || n % 11 == 0) { return false; }
 
  // Проверяем делимость на остальные числа, не делящиеся на 2, 3, 5, 7 и 11, начиная с 13
  // (с приращением, взятым из массива delta): 13, 17, 19, 23, 29, 31, 37, 41, 43, 47...
  static const unsigned char delta[480] = {
    4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 4,
    14, 4, 6, 2, 10, 2, 6, 6, 4, 2, 4, 6, 2, 10, 2, 4, 2, 12, 10, 2, 4, 2, 4, 6,
    2, 6, 4, 6, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 6, 8, 6, 10, 2, 4, 6,
    2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 6, 10, 2, 10, 2, 4, 2, 4, 6, 8, 4, 2, 4,
    12, 2, 6, 4, 2, 6, 4, 6, 12, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 10, 2,
    4, 6, 2, 6, 4, 2, 4, 2, 10, 2, 10, 2, 4, 6, 6, 2, 6, 6, 4, 6, 6, 2, 6, 4,
    2, 6, 4, 6, 8, 4, 2, 6, 4, 8, 6, 4, 6, 2, 4, 6, 8, 6, 4, 2, 10, 2, 6, 4,
    2, 4, 2, 10, 2, 10, 2, 4, 2, 4, 8, 6, 4, 2, 4, 6, 6, 2, 6, 4, 8, 4, 6, 8,
    4, 2, 4, 2, 4, 8, 6, 4, 6, 6, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2,
    10, 2, 10, 2, 6, 4, 6, 2, 6, 4, 2, 4, 6, 6, 8, 4, 2, 6, 10, 8, 4, 2, 4, 2,
    4, 8, 10, 6, 2, 4, 8, 6, 6, 4, 2, 4, 6, 2, 6, 4, 6, 2, 10, 2, 10, 2, 4, 2,
    4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 6, 6, 4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4,
    8, 4, 6, 2, 6, 6, 4, 2, 4, 6, 8, 4, 2, 4, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2,
    10, 2, 4, 6, 8, 6, 4, 2, 6, 4, 6, 8, 4, 6, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2,
    6, 6, 4, 6, 6, 2, 6, 6, 4, 2, 10, 2, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 10, 6,
    2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4, 2, 12, 6, 4, 6, 2, 4, 6, 2, 12, 4, 2, 4,
    8, 6, 4, 2, 4, 2, 10, 2, 10, 6, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2,
    10, 6, 8, 6, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 6, 4, 6, 2, 6, 4, 2,
    4, 2, 10, 12, 2, 4, 2, 10, 2, 6, 4, 2, 4, 6, 6, 2, 10, 2, 6, 4, 14, 4, 2, 4,
    2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 12, 2, 12 };  // массив приращений
 
  for (unsigned int div = 13, idx = 0, max_div = sqrt(n);  // idx - индекс в массиве delta
    div <= max_div;
    div += delta[idx], idx = (idx == 479 ? 0 : idx + 1)) {  // увеличиваем индекс с остатком от деления на 480
    if (n % div == 0) { return false; }
  }
 
  // Число простое
  return true;
}




Мой спич был бы не полным, если бы я не привёл результаты теста скорости.
Итак, запускаем проверку чисел от 0 до 10'000'000 на простоту на моём домашнем компьютере (Intel Core i5-2500K / Windows 10 x64).
Я использовал компиляторы GCC (g++) и Clang с ключом оптимизации -O2 для платформ x86 и x64.
Результаты между различными компиляторами и платформами оказались настолько близкими (с расхождением порядка 1%), что я решил не разделять их, а объединить в одну таблицу.
Здесь указано значение в секундах, доля и ускорение относительно первого варианта.

Вариант
Время
%
Ускорение
is_prime
7.50
100%
1x
is_prime2
3.80
51%
2.0x
is_prime3
2.52
34%
3.0x
is_prime5
2.02
27%
3.7x
is_prime7
1.75
23%
4.3x
is_prime11
1.60
21%
4.7x

Используйте ссылки в таблице выше, посмотреть и протестировать код на платформе Try It Online.

Исходники с тестовым кодом (более полным, чем на tio.run) прикреплены к статье в виде архива.

На этом всё.
Если вы знаете более быстрые варианты точного определения простоты числа, буду рад, если поделитесь в комментариях



update 04.03.2019: добавлены функции is_prime7, is_prime11 (в спойлерах), тип изменён на unsigned int. Архив перезалит.
Вложения
Тип файла: zip is_prime.zip (6.6 Кб, 29 просмотров)
Просмотров 457 Комментарии 6
Всего комментариев 6
Комментарии
  1. Старый комментарий
    Уважаемый Jin X,
    я решал эту задачу. И естественно нашел более быстрый и простой алгоритм. Суть его том, что в качестве делителей берутся только простые числа. Для этого объявляется массив нужного размера и заполняется простыми числами. Конкретно я брал массив длиной порядка 45000. Очевидно, что массив либо заполнялся с файла (простых чисел), либо был (инициировался) при загрузке программы.
    Вариант 2
    Вычислялись все простые числа до 2 000 000 000 и записывались в файл прямого доступа. Оставалось только проверить: находится ли данное число в файле?
    Запись от нтч размещена 28.02.2019 в 08:02 нтч на форуме
  2. Старый комментарий
    Аватар для Jin X
    нтч, брать значение из массива (или лучше битового набора – для 1 млрд чисел нужно ≈ 120 Мб памяти) – это очевидно самый быстрый способ. Для вычисления можно даже не заморачиваться с файлом, а использовать решето Эратосфена или решето Аткина, к примеру (правда, я их ещё не тестил их на скорость).

    Но такой расход памяти не всегда оправдан. К примеру, если вам нужно лишь несколько значений из заранее неизвестного (очень широкого) диапазона. Либо если вы не знаете заранее, понадобятся ли вам эти значения – смысл держать выделенные мегабайты памяти?

    Давайте в рамках этой статьи ограничимся разовым тестом числа на простоту без предварительных расчётов. Есть же всякие тесты Миллера и пр. Вопрос только в том, насколько они быстрые и не чересчур ли сложные?
    А что касается решёт, то когда у меня дойдут до них руки, я протестирую их скорость и напишу об этом отдельную статью
    p.s. По-хорошему, надо ещё разобраться со способами факторизации числа, нахождения списка делителей (кроме самых очевидных).
    Запись от Jin X размещена 28.02.2019 в 09:32 Jin X вне форума
    Обновил(-а) Jin X 28.02.2019 в 10:08
  3. Старый комментарий
    Запись от Croessmah размещена 02.03.2019 в 14:33 Croessmah вне форума
  4. Старый комментарий
    Запись от bedvit размещена 02.03.2019 в 17:22 bedvit вне форума
  5. Старый комментарий
    Аватар для Jin X
    Цитата:
    Спасибо, сохранил несколько ссылок для просмотра. А алгоритм из шапки работает с той же скоростью (потому что, по сути, на том же принципе основан), что и мой последний (is_prime5)
    Правда, на tio.run почему-то в 2 раза медленнее...

    Цитата:
    Значит, я сам того не зная, применил у себя этот wheel optimization
    Распараллеливание – это интересная идея, но полагаю, что числа должны быть достаточно большими + их должно быть несколько, иначе процесс запуска/остановки распараллеливания будет долгим.

    Кстати, я уже сделал тест Миллера (детерминированный, а не вероятностный Миллера-Рабина). Работает довольно быстро (лично у меня на компе прирост начинается с отметки примерно в полмиллиона, на tio.run эта отметка ниже).

    За решето ещё не брался
    Запись от Jin X размещена 04.03.2019 в 18:07 Jin X вне форума
  6. Старый комментарий
    Аватар для Jin X
    Добавил функции is_prime7, is_prime11 (см. в спойлерах), тип изменил на unsigned int.
    Архив перезалил.
    Запись от Jin X размещена 04.03.2019 в 20:04 Jin X вне форума
    Обновил(-а) Jin X 05.03.2019 в 21:24
 
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.