# Primitive roots density of number

## Definition:

A number m is called a primitive root of a prime p the condition that the smallest integer k for which p dividies mk-1 is p-1

Given a tuple (a,b) of positive integers, return the fraction:

(number of primes p equal to or less than a which has b as a primitive root) divided by (number of primes equal to or less than a)

## Test Cases

Assuming Artin's Conjecture to be true, giving (A very big number, 2) should return a value around .3739558136...

The challenge is for the fastest code, not necessarily the shortest one.

## Timing

My crappy code takes 10 seconds for a million on my PC (intel i3 6th Gen, 8GB RAM, others I forgot), so your code should naturally beat that speed.

This is my first challenge. Feeback is appreciated. :)

• Darn why there's no latex in CG SE :/ Sep 29, 2017 at 7:04
• Type it here, and post it as images. That's what I did in this one. Sep 29, 2017 at 7:12
• Can you add some test cases, please? Sep 29, 2017 at 7:24
• codegolf.stackexchange.com/tags/fastest-code/info Sep 29, 2017 at 8:03
• Why do you have a laptop made by Volkswagen??? Sep 29, 2017 at 10:25

# C (gcc)

Around 0.6 sec for 1,000,000.

int modpow(long b,int e,long n){
long r=1;
while(e){
if(e&1) r=r*b%n;
b=b*b%n;
e>>=1;
}
return r;
}
float f(int a,int b){
int ln=1;
for(int p=2;p<=a;p*=ln,ln++);
int* primes = malloc((a+1)*(sizeof(int)));
int** factors = malloc((a+1)*sizeof(int*));
primes[0] = 0;
primes[1] = 0;
for(int i=2;i<=a;i++){
primes[i] = 1;
factors[i] = malloc(ln*sizeof(int));
}
factors[1] = malloc(sizeof(int));
factors[1][1] = 0;
for(int p=2;p*p<=a;p++){
if(primes[p]){
for(int j=p*2;j<=a;j+=p){
primes[j] = 0;
factors[j][0]++;
factors[j][factors[j][0]] = p;
}
}
}
int count1 = 0;
int count2 = 0;
for(int p=2;p<=a;p++){
if(!primes[p]) continue;
count2++;
int n=p-1;
int m=p-1;
int prim=1;
for(int j=1;j<=factors[n][0];j++){
int q=factors[n][j];
if(modpow(b,n/q,p)==1){
prim=0;
break;
}
while(m%q==0){
m/=q;
}
}
if(m>1){
if(modpow(b,n/m,p)==1){
prim=0;
}
}
count1 += prim;
}
return (float)count1/count2;
}


Try it online!

• hey! +1 from me then! Sep 29, 2017 at 22:40
• @Jenny_mathy it's better now once I switch to C Sep 30, 2017 at 7:20

# Mathematica

(t=0;
f=PrimePi[#];
d=#2;
For[i=1,i<=f,i++,If[PrimitiveRoot[Prime[i],d]==d,t++]];
N[t/f])&[1000000,2]//AbsoluteTiming


Try it online Paste the code and press shift-enter

Takes 3.2 sec to find 1.000.000

• This should have PrimitiveRoot[Prime[i],d]==d to work in full generality. Otherwise, given input b=3, it will not count a prime that has 3 as a primitive root if it also has 2 as a primitive root, because it will find 2 first. Sep 29, 2017 at 14:43
• @MishaLavrov if you see the edit history you'll see that this was my first approach but then I thought it wasn't necessary. I'll change it back. thanks! Sep 29, 2017 at 15:05
• @Jenny_mathy hey we're on the same order of magnitude :P Sep 29, 2017 at 17:22
• I believe that NextPrime might be a bit faster than Prime. Sep 30, 2017 at 15:24

# C++11 + libop

1.23 sec on my slow laptop for 10,000,000.

#include <cstdint>
#include <vector>
#include <algorithm>
#include <iostream>
#include <utility>

#include "libop/op.h"

// Returns 2^-32 mod m, -m^(-1) mod 2^32
inline std::pair<uint64_t, uint64_t> mont_modinv32(uint64_t m) {
uint64_t a = 1ull << 31;
uint64_t u = 1;
uint64_t v = 0;

while (a > 0) {
a = a >> 1;
if ((u & 1) == 0) {
u = u >> 1; v = v >> 1;
} else {
u = ((u ^ m) >> 1) + (u & m);
v = (v >> 1) + (1ull << 31);
}
}

return std::make_pair(u, v);
}

// Returns (ab)R mod n given aR mod n, bR mod n, n and -n^(-1) mod R, with R = 2^32
inline uint64_t montmul32(uint64_t a, uint64_t b, uint64_t n, uint64_t nneginv) {
uint64_t T = a*b;
uint32_t m = T*nneginv; // m = T*-n^(-1) (mod 2^32)
uint64_t t = (T + m*n) >> 32;
return t >= n ? t - n : t;
}

double f(uint32_t a, uint32_t b) {
std::vector<uint32_t> primes;
op::primes_below(a + 1, std::back_inserter(primes));

int num_primitives = 0;

for (auto p : primes) {
if (p == 2) {
num_primitives += b == 1;
continue;
}

uint32_t s = p - 1;
uint32_t pneginv = mont_modinv32(p).second;
uint32_t montb = (uint64_t(b) << 32) % p;
uint64_t mont1 = (1ull << 32) % p;

uint32_t pows[32] = {montb};
for (int i = 1; s >> i; ++i) {
montb = montmul32(montb, montb, p, pneginv);
pows[i] = montb;
}

bool primitive_root = true;
for (auto& kv : op::factorization(s)) {
uint32_t e = s / kv.first;

uint64_t r = mont1;
for (int i = 0; e; ++i) {
if (e & 1) r = montmul32(r, pows[i], p, pneginv);
e >>= 1;
}

if (r == mont1) {
primitive_root = false;
break;
}
}

num_primitives += primitive_root;
}

return num_primitives / double(primes.size());
}