Guest User

Untitled

a guest
Mar 6th, 2023
117
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 3.29 KB | Source Code | 0 0
  1. /*
  2. Program to quickly compute orders of integers modulo primes.
  3.  
  4. See: https://mathoverflow.net/a/436684/496479
  5.  
  6. You may use this code and its modifications for any purposes you like.
  7. */
  8.  
  9. #include <iostream>
  10. #include <vector>
  11. #include <bitset>
  12. #include <iomanip>
  13. #include <algorithm>
  14. #include <math.h>
  15.  
  16. using namespace std;
  17.  
  18. typedef long long ll;
  19. typedef pair<ll, ll> pii;
  20. typedef vector<int> vi;
  21. typedef long double ld;
  22.  
  23. const ll N = 4e4;
  24.  
  25. bitset<N> pr; //Prime/composite info for n < N.
  26.  
  27. vector<int> di[N+1]; //Vector procedure for maintaining information on divisors
  28.  
  29. //Sieve
  30. void sieve() {
  31.     for(int i = 2; i < N; ++i) {
  32.         pr[i] = 1;
  33.     }
  34.     for(ll p = 2; p < N; ++p) {
  35.         if(!pr[p]) continue;
  36.         di[p].push_back(p);
  37.         for(ll j = p*p; j < N; j+=p) {
  38.             pr[j] = 0;
  39.         }
  40.     }
  41.  
  42. }
  43.  
  44. //Calculate a^e (mod p) in O(log e) time
  45. ll exp(ll a, ll e, ll p) {
  46.     if(e == 0) return 1;
  47.     if(e%2 == 0) {
  48.         ll b = exp(a, e/2, p);
  49.         return (b*b)%p;
  50.     }
  51.     ll b = exp(a, e-1, p);
  52.     return (a*b)%p;
  53. }
  54.  
  55. //Calculate q^{v_q(ord_p(a))}. It is given that v_q(p-1) = e.
  56. ll ordPart(ll a, ll p, ll q, ll e) {
  57.     ll QU = (p-1)/exp(q, e, p);
  58.     ll qu = exp(a, QU, p);
  59.     ll A = e;
  60.     ll ans = 1;
  61.     while(qu != 1 && A > 0) {
  62.         ans *= q;
  63.         ans %= p;
  64.         qu = exp(qu, q, p);
  65.         --A;
  66.     }
  67.     return ans;
  68. }
  69.  
  70. //Calculate ord_p(a)
  71. ll ord(ll a, ll p) {
  72.     //Take a mod p
  73.     a %= p;
  74.     if(a < 0) {
  75.         a += p;
  76.     }
  77.  
  78.     ll ans = 1;
  79.  
  80.     ll n = p-1;
  81.     for(auto j : di[(p-1)%N]) {
  82.         int e = 0;
  83.         while(n%j == 0) {
  84.             ++e;
  85.             n/=j;
  86.         }
  87.         //For each prime j dividing p-1, multiply answer by j-part ord_p(a).
  88.         ans *= ordPart(a, p, j, e);
  89.     }
  90.  
  91.     if(n > 1) {
  92.         //p-1 has a large prime factor n.
  93.         ans *= ordPart(a, p, n, 1);
  94.     }
  95.     return ans;
  96. }
  97.  
  98.  
  99. int main() {
  100.     //Perform sieve computation
  101.     sieve();
  102.  
  103.     //As an example, we compute the number of primes p for which a is a primitive root mod p (ord(a, p) = p-1).
  104.     //Read integer a
  105.  
  106.     int a;
  107.     cin >> a;
  108.  
  109.     ll am = 0; //Amount of primes such that a is a primitive root modulo p
  110.     ll pi = 0; //Amount of primes p (encountered so far)
  111.  
  112.     int B = 1e7;
  113.  
  114.     //To make the output a bit cleaner
  115.     cout << setprecision(8);
  116.     cout << fixed;
  117.  
  118.  
  119.     for(ll p = 2; p < 1e9; ++p) { //(Note: current implementation starts to break down if p is too large (p*p > 2^64, i.e. p > 4e9), because of multiplication overflow in exp-function)
  120.    
  121.         //Print partial results in once in a while
  122.         if(p%B == 0) {
  123.             cout << p << " " << am << " " << pi << "\n";
  124.             cout << (double) am/pi << "\n\n";
  125.         }
  126.  
  127.         //Push the (small) prime factors of p-2 forward
  128.         for(auto j : di[(p-2)%N]) {
  129.             di[(p+j-2)%N].push_back(j);
  130.         }
  131.         di[(p-2)%N].clear();
  132.  
  133. //      If you run into memory problems, remember to shrink the vectors di[].
  134. //      Frees up memory if the current length of di[] is much shorter than the maximum length it has been.
  135. //      di(p-2)%N].shrink_to_fit();
  136.  
  137.  
  138.         //If p is not a prime, just continue
  139.         if(p < N) {
  140.             if(!pr[p]) continue;
  141.         } else {
  142.             if(di[p%N].size() > 0) continue;
  143.         }
  144.  
  145.         //p is prime
  146.         ++pi;
  147.         if(a >= p) {
  148.             if(a%p == 0) continue;
  149.         }
  150.  
  151.         //Check if a is a primitive root mod p
  152.         //(Note: the specific problem "is a primitive root modulo p" could be done slightly faster than computing ord(a. p) and checking whether it is p-1 by tweaking the code above.)
  153.         if(ord(a, p) == p-1) ++am;
  154.     }
  155.  
  156. }
Add Comment
Please, Sign In to add comment