runewalsh

fuck

Jun 4th, 2022 (edited)
370
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 5.76 KB | None | 0 0
  1. // Вместо miller_rabin_probably_prime(tries=так, чтобы точно хватило) использовать небольшой, фиксированный набор оснований,
  2. // гарантирующий отсутствие ложноположительных срабатываний на всём диапазоне входных данных.
  3. #define USE_PREDEFINED_MILLER_RABIN_BASES
  4.  
  5. typedef unsigned int uint;
  6. typedef unsigned long long ulonglong;
  7.  
  8. // Делит n на fuck, пока делится. Записывает частное назад в n и возвращает число делений. n должно изначально делиться на fuck!
  9. uint extract_fucks(uint* n, uint fuck)
  10. {
  11.     uint nfuck = 0;
  12.     do { nfuck++, *n /= fuck; } while (*n % fuck == 0);
  13.     return nfuck;
  14. }
  15.  
  16. const uint small_fucks[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53 };
  17. // После 2 и 3, остальные простые числа идут начиная с 5 с шагами 2-4.
  18. const uint first_large_fuck = 59, first_large_fuck_d42 = first_large_fuck % 6 == 5 ? 2 : 4;
  19.  
  20. // Возвращает число делителей n, с учётом 1 и n.
  21. // Работает через факторизацию перебором.
  22. uint ndivs_ver1(uint n)
  23. {
  24.     if (n <= 1) return 2 - n;
  25.     uint r = 1, n_left = n, sqr = (uint)sqrt(n_left);
  26.    
  27.     for (uint fuck : small_fucks)
  28.     {
  29.         if (fuck > sqr) return 2 * r;
  30.         if (n_left % fuck == 0)
  31.         {
  32.             r *= 1 + extract_fucks(&n_left, fuck);
  33.             if (n_left == 1) return r;
  34.             sqr = (uint)sqrt(n_left);
  35.         }
  36.     }
  37.  
  38.     for (uint fuck = first_large_fuck, dfuck = first_large_fuck_d42; fuck <= sqr; fuck += dfuck, dfuck = 6 - dfuck)
  39.         if (n_left % fuck == 0)
  40.         {
  41.             r *= 1 + extract_fucks(&n_left, fuck);
  42.             if (n_left == 1) return r;
  43.             sqr = (uint)sqrt(n_left);
  44.         }
  45.     return 2 * r;
  46. }
  47.  
  48. // Номер старшего выставленного бита. Не определён для 0.
  49. int floor_log2(uint x)
  50. {
  51. #if defined(_MSC_VER)
  52.     unsigned long index;
  53.     _BitScanReverse(&index, x);
  54.     return index;
  55. #elif defined(__GNUC__)
  56.     return CHAR_BIT * sizeof(x) - 1 - __builtin_clz(x);
  57. #else
  58.     int r = -1;
  59.     while (x > 0) r++, x >>= 1;
  60.     return r;
  61. #endif
  62. }
  63.  
  64. // Номер младшего выставленного бита. Не определён для 0.
  65. int ffs(uint x)
  66. {
  67. #if defined(_MSC_VER)
  68.     unsigned long index;
  69.     _BitScanForward(&index, x);
  70.     return index;
  71. #elif defined(__GNUC__)
  72.     return __builtin_ctz(x);
  73. #else
  74.     for (int i = 0; i < CHAR_BIT * sizeof(x); i++)
  75.         if (x >> i & 1) return i;
  76.     return -1;
  77. #endif
  78. }
  79.  
  80. // b ** e % m.
  81. uint pow_mod(uint b, uint e, uint m)
  82. {
  83.     if (e == 0) return 1;
  84.     uint r = b;
  85.     for (int i = floor_log2(e) - 1; i >= 0; i--)
  86.     {
  87.         r = (ulonglong)r * r % m;
  88.         if (e >> i & 1) r = (ulonglong)r * b % m;
  89.     }
  90.     return r;
  91. }
  92.  
  93. // Вспомогательная функция для miller_rabin_probably_prime, раунд теста Миллера-Рабина.
  94. bool miller_rabin_round(uint n, uint y, uint q, uint k)
  95. {
  96.     y = pow_mod(y, q, n);
  97.     if (y == 1 || y == n - 1) return true;
  98.     for (;;)
  99.     {
  100.         if (!--k) return false;
  101.         y = (ulonglong)y * y % n;
  102.         if (y == n - 1) return true;
  103.         if (y <= 1) return false;
  104.     }
  105. }
  106.  
  107. #ifndef USE_PREDEFINED_MILLER_RABIN_BASES
  108. // Стохастический ответ, является ли число n простым. n должно быть нечётным > 2.
  109. // Отрицательный ответ всегда верен, максимальная вероятность ложноположительного ответа — 1/4**tries.
  110. bool miller_rabin_probably_prime(uint n, uint tries)
  111. {
  112.     uint k = ffs(n - 1), q = (n - 1) >> k;
  113.     if (!miller_rabin_round(n, 2, q, k)) return false;
  114.     for (uint itry = 0; itry + 1 < tries; itry++)
  115.     {
  116.         uint y = itry * itry + itry + 41;
  117.         if (y >= n - 1) return true;
  118.         if (!miller_rabin_round(n, y, q, k)) return false;
  119.     }
  120.     return true;
  121. }
  122. #endif
  123.  
  124. // Вспомогательная функция для ndivs_ver2, проверяющая, является ли n простым.
  125. // Не называется просто is_prime, потому что может в другой раз иметь особые требования к n.
  126. bool ndivs_ver2_test_prime(uint n)
  127. {
  128. #ifdef USE_PREDEFINED_MILLER_RABIN_BASES
  129.     // Версия miller_rabin_probably_prime с заготовленными основаниями.
  130.     // https ://miller-rabin.appspot .com
  131.     // Эти покрывают диапазон до 15,4 млрд., то есть 32-битный uint.
  132.     static const uint bases[] = { 15, 176006322, 4221622697 };
  133.     uint k = ffs(n - 1), q = (n - 1) >> k;
  134.     for (uint y : bases)
  135.         if (!miller_rabin_round(n, y, q, k)) return false;
  136.     return true;
  137. #else
  138.     return miller_rabin_probably_prime(n, 24);
  139. #endif
  140. }
  141.  
  142. // Возвращает число делителей n, с учётом 1 и n.
  143. // Работает через факторизацию перебором с дополнительной проверкой на простоту.
  144. uint ndivs_ver2(uint n)
  145. {
  146.     if (n <= 1) return 2 - n;
  147.     uint r = 1, n_left = n, sqr = (uint)sqrt(n_left);
  148.  
  149.     for (uint fuck : small_fucks)
  150.     {
  151.         if (fuck > sqr) return 2 * r;
  152.         if (n_left % fuck == 0)
  153.         {
  154.             r *= 1 + extract_fucks(&n_left, fuck);
  155.             if (n_left == 1) return r;
  156.             sqr = (uint)sqrt(n_left);
  157.         }
  158.     }
  159.  
  160.     if (first_large_fuck > sqr || ndivs_ver2_test_prime(n_left)) return 2 * r;
  161.     for (uint fuck = first_large_fuck, dfuck = first_large_fuck_d42; fuck <= sqr; fuck += dfuck, dfuck = 6 - dfuck)
  162.         if (n_left % fuck == 0)
  163.         {
  164.             r *= 1 + extract_fucks(&n_left, fuck);
  165.             if (n_left == 1) return r;
  166.             sqr = (uint)sqrt(n_left);
  167.             if (fuck + dfuck > sqr || ndivs_ver2_test_prime(n_left)) return 2 * r;
  168.         }
  169.     return 2 * r;
  170. }
Add Comment
Please, Sign In to add comment