Why modular exponentiations of random, coprime numbers equals to 1/3?

41 Views Asked by At

Why is the average of modular exponentiations of random, coprime numbers equal to 1/3 on average? A similar question for N random (not coprime) numbers, but 'Why ... equal to 1/2 on average'? C++ code that demonstrates this:

#include <ctime>
#include <iostream>
#include <numeric>
#include <random>

using namespace std;

uint64_t modexp(uint64_t x, uint64_t y, uint64_t N) {
    if (y == 0)
        return 1;
    uint64_t z = modexp(x, y / 2, N);
    if (y % 2 == 0)
        return (z * z) % N;
    else
        return (x * z * z) % N;
}

#define COPRIME
// #define RANDOM

int main() {
    const uint64_t N = 1000000;
    default_random_engine dre;
    uniform_int_distribution ud(uint64_t(2), uint64_t(1) << 31);
    dre.seed(time(0));
    uint64_t a, n, c = 0;
    for (uint64_t i = 0; i < N; ++i) {
#ifdef COPRIME
        while (gcd(a = ud(dre), n = ud(dre)) != 1) {}
#endif
#ifdef RANDOM
        a = ud(dre);
        n = ud(dre);
#endif
        c += !modexp(a, n - 1, uint64_t(1) << 31);
    }
    cout << c / double(N);
}

Example outputs when #define COPRIME and N = 1000000: 0.33252 0.333335 0.333448; when #define RANDOM: 0.500039 0.500289 0.499618. This correlation appears to be more accurate when N is larger.

EDIT: Not 'modular exponentiations' but probability of them to equal to 0 is 1/2 or 1/3. ('c += !modexp(a, n - 1, uint64_t(1) << 31)' is the same as 'c += (modexp(a, n - 1, uint64_t(1) << 31)) == 0'))