Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- /*
- Program to quickly compute orders of integers modulo primes.
- See: https://mathoverflow.net/a/436684/496479
- You may use this code and its modifications for any purposes you like.
- */
- #include <iostream>
- #include <vector>
- #include <bitset>
- #include <iomanip>
- #include <algorithm>
- #include <math.h>
- using namespace std;
- typedef long long ll;
- typedef pair<ll, ll> pii;
- typedef vector<int> vi;
- typedef long double ld;
- const ll N = 4e4;
- bitset<N> pr; //Prime/composite info for n < N.
- vector<int> di[N+1]; //Vector procedure for maintaining information on divisors
- //Sieve
- void sieve() {
- for(int i = 2; i < N; ++i) {
- pr[i] = 1;
- }
- for(ll p = 2; p < N; ++p) {
- if(!pr[p]) continue;
- di[p].push_back(p);
- for(ll j = p*p; j < N; j+=p) {
- pr[j] = 0;
- }
- }
- }
- //Calculate a^e (mod p) in O(log e) time
- ll exp(ll a, ll e, ll p) {
- if(e == 0) return 1;
- if(e%2 == 0) {
- ll b = exp(a, e/2, p);
- return (b*b)%p;
- }
- ll b = exp(a, e-1, p);
- return (a*b)%p;
- }
- //Calculate q^{v_q(ord_p(a))}. It is given that v_q(p-1) = e.
- ll ordPart(ll a, ll p, ll q, ll e) {
- ll QU = (p-1)/exp(q, e, p);
- ll qu = exp(a, QU, p);
- ll A = e;
- ll ans = 1;
- while(qu != 1 && A > 0) {
- ans *= q;
- ans %= p;
- qu = exp(qu, q, p);
- --A;
- }
- return ans;
- }
- //Calculate ord_p(a)
- ll ord(ll a, ll p) {
- //Take a mod p
- a %= p;
- if(a < 0) {
- a += p;
- }
- ll ans = 1;
- ll n = p-1;
- for(auto j : di[(p-1)%N]) {
- int e = 0;
- while(n%j == 0) {
- ++e;
- n/=j;
- }
- //For each prime j dividing p-1, multiply answer by j-part ord_p(a).
- ans *= ordPart(a, p, j, e);
- }
- if(n > 1) {
- //p-1 has a large prime factor n.
- ans *= ordPart(a, p, n, 1);
- }
- return ans;
- }
- int main() {
- //Perform sieve computation
- sieve();
- //As an example, we compute the number of primes p for which a is a primitive root mod p (ord(a, p) = p-1).
- //Read integer a
- int a;
- cin >> a;
- ll am = 0; //Amount of primes such that a is a primitive root modulo p
- ll pi = 0; //Amount of primes p (encountered so far)
- int B = 1e7;
- //To make the output a bit cleaner
- cout << setprecision(8);
- cout << fixed;
- 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)
- //Print partial results in once in a while
- if(p%B == 0) {
- cout << p << " " << am << " " << pi << "\n";
- cout << (double) am/pi << "\n\n";
- }
- //Push the (small) prime factors of p-2 forward
- for(auto j : di[(p-2)%N]) {
- di[(p+j-2)%N].push_back(j);
- }
- di[(p-2)%N].clear();
- // If you run into memory problems, remember to shrink the vectors di[].
- // Frees up memory if the current length of di[] is much shorter than the maximum length it has been.
- // di(p-2)%N].shrink_to_fit();
- //If p is not a prime, just continue
- if(p < N) {
- if(!pr[p]) continue;
- } else {
- if(di[p%N].size() > 0) continue;
- }
- //p is prime
- ++pi;
- if(a >= p) {
- if(a%p == 0) continue;
- }
- //Check if a is a primitive root mod p
- //(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.)
- if(ord(a, p) == p-1) ++am;
- }
- }
Add Comment
Please, Sign In to add comment