6
\$\begingroup\$

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. :)

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

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!

\$\endgroup\$
  • \$\begingroup\$ hey! +1 from me then! \$\endgroup\$ – J42161217 Sep 29 '17 at 22:40
  • \$\begingroup\$ @Jenny_mathy it's better now once I switch to C \$\endgroup\$ – Leaky Nun Sep 30 '17 at 7:20
2
\$\begingroup\$

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

\$\endgroup\$
  • \$\begingroup\$ 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. \$\endgroup\$ – Misha Lavrov Sep 29 '17 at 14:43
  • \$\begingroup\$ @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! \$\endgroup\$ – J42161217 Sep 29 '17 at 15:05
  • \$\begingroup\$ @Jenny_mathy hey we're on the same order of magnitude :P \$\endgroup\$ – Leaky Nun Sep 29 '17 at 17:22
  • \$\begingroup\$ I believe that NextPrime might be a bit faster than Prime. \$\endgroup\$ – LegionMammal978 Sep 30 '17 at 15:24
2
\$\begingroup\$

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());
}
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.