Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- // Вместо miller_rabin_probably_prime(tries=так, чтобы точно хватило) использовать небольшой, фиксированный набор оснований,
- // гарантирующий отсутствие ложноположительных срабатываний на всём диапазоне входных данных.
- #define USE_PREDEFINED_MILLER_RABIN_BASES
- typedef unsigned int uint;
- typedef unsigned long long ulonglong;
- // Делит n на fuck, пока делится. Записывает частное назад в n и возвращает число делений. n должно изначально делиться на fuck!
- uint extract_fucks(uint* n, uint fuck)
- {
- uint nfuck = 0;
- do { nfuck++, *n /= fuck; } while (*n % fuck == 0);
- return nfuck;
- }
- const uint small_fucks[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53 };
- // После 2 и 3, остальные простые числа идут начиная с 5 с шагами 2-4.
- const uint first_large_fuck = 59, first_large_fuck_d42 = first_large_fuck % 6 == 5 ? 2 : 4;
- // Возвращает число делителей n, с учётом 1 и n.
- // Работает через факторизацию перебором.
- uint ndivs_ver1(uint n)
- {
- if (n <= 1) return 2 - n;
- uint r = 1, n_left = n, sqr = (uint)sqrt(n_left);
- for (uint fuck : small_fucks)
- {
- if (fuck > sqr) return 2 * r;
- if (n_left % fuck == 0)
- {
- r *= 1 + extract_fucks(&n_left, fuck);
- if (n_left == 1) return r;
- sqr = (uint)sqrt(n_left);
- }
- }
- for (uint fuck = first_large_fuck, dfuck = first_large_fuck_d42; fuck <= sqr; fuck += dfuck, dfuck = 6 - dfuck)
- if (n_left % fuck == 0)
- {
- r *= 1 + extract_fucks(&n_left, fuck);
- if (n_left == 1) return r;
- sqr = (uint)sqrt(n_left);
- }
- return 2 * r;
- }
- // Номер старшего выставленного бита. Не определён для 0.
- int floor_log2(uint x)
- {
- #if defined(_MSC_VER)
- unsigned long index;
- _BitScanReverse(&index, x);
- return index;
- #elif defined(__GNUC__)
- return CHAR_BIT * sizeof(x) - 1 - __builtin_clz(x);
- #else
- int r = -1;
- while (x > 0) r++, x >>= 1;
- return r;
- #endif
- }
- // Номер младшего выставленного бита. Не определён для 0.
- int ffs(uint x)
- {
- #if defined(_MSC_VER)
- unsigned long index;
- _BitScanForward(&index, x);
- return index;
- #elif defined(__GNUC__)
- return __builtin_ctz(x);
- #else
- for (int i = 0; i < CHAR_BIT * sizeof(x); i++)
- if (x >> i & 1) return i;
- return -1;
- #endif
- }
- // b ** e % m.
- uint pow_mod(uint b, uint e, uint m)
- {
- if (e == 0) return 1;
- uint r = b;
- for (int i = floor_log2(e) - 1; i >= 0; i--)
- {
- r = (ulonglong)r * r % m;
- if (e >> i & 1) r = (ulonglong)r * b % m;
- }
- return r;
- }
- // Вспомогательная функция для miller_rabin_probably_prime, раунд теста Миллера-Рабина.
- bool miller_rabin_round(uint n, uint y, uint q, uint k)
- {
- y = pow_mod(y, q, n);
- if (y == 1 || y == n - 1) return true;
- for (;;)
- {
- if (!--k) return false;
- y = (ulonglong)y * y % n;
- if (y == n - 1) return true;
- if (y <= 1) return false;
- }
- }
- #ifndef USE_PREDEFINED_MILLER_RABIN_BASES
- // Стохастический ответ, является ли число n простым. n должно быть нечётным > 2.
- // Отрицательный ответ всегда верен, максимальная вероятность ложноположительного ответа — 1/4**tries.
- bool miller_rabin_probably_prime(uint n, uint tries)
- {
- uint k = ffs(n - 1), q = (n - 1) >> k;
- if (!miller_rabin_round(n, 2, q, k)) return false;
- for (uint itry = 0; itry + 1 < tries; itry++)
- {
- uint y = itry * itry + itry + 41;
- if (y >= n - 1) return true;
- if (!miller_rabin_round(n, y, q, k)) return false;
- }
- return true;
- }
- #endif
- // Вспомогательная функция для ndivs_ver2, проверяющая, является ли n простым.
- // Не называется просто is_prime, потому что может в другой раз иметь особые требования к n.
- bool ndivs_ver2_test_prime(uint n)
- {
- #ifdef USE_PREDEFINED_MILLER_RABIN_BASES
- // Версия miller_rabin_probably_prime с заготовленными основаниями.
- // https ://miller-rabin.appspot .com
- // Эти покрывают диапазон до 15,4 млрд., то есть 32-битный uint.
- static const uint bases[] = { 15, 176006322, 4221622697 };
- uint k = ffs(n - 1), q = (n - 1) >> k;
- for (uint y : bases)
- if (!miller_rabin_round(n, y, q, k)) return false;
- return true;
- #else
- return miller_rabin_probably_prime(n, 24);
- #endif
- }
- // Возвращает число делителей n, с учётом 1 и n.
- // Работает через факторизацию перебором с дополнительной проверкой на простоту.
- uint ndivs_ver2(uint n)
- {
- if (n <= 1) return 2 - n;
- uint r = 1, n_left = n, sqr = (uint)sqrt(n_left);
- for (uint fuck : small_fucks)
- {
- if (fuck > sqr) return 2 * r;
- if (n_left % fuck == 0)
- {
- r *= 1 + extract_fucks(&n_left, fuck);
- if (n_left == 1) return r;
- sqr = (uint)sqrt(n_left);
- }
- }
- if (first_large_fuck > sqr || ndivs_ver2_test_prime(n_left)) return 2 * r;
- for (uint fuck = first_large_fuck, dfuck = first_large_fuck_d42; fuck <= sqr; fuck += dfuck, dfuck = 6 - dfuck)
- if (n_left % fuck == 0)
- {
- r *= 1 + extract_fucks(&n_left, fuck);
- if (n_left == 1) return r;
- sqr = (uint)sqrt(n_left);
- if (fuck + dfuck > sqr || ndivs_ver2_test_prime(n_left)) return 2 * r;
- }
- return 2 * r;
- }
Add Comment
Please, Sign In to add comment