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

Your task
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. :)
C (gcc)

Around 0.6 sec for 1,000,000.



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


Try it online!