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'))