Miller-Rabin Strong Pseudoprimes

Question

Given a non-negative integer N, output the smallest odd positive integer that is a strong pseudoprime to all of the first N prime bases.

This is OEIS sequence A014233.

Test Cases (one-indexed)

1       2047
2       1373653
3       25326001
4       3215031751
5       2152302898747
6       3474749660383
7       341550071728321
8       341550071728321
9       3825123056546413051
10      3825123056546413051
11      3825123056546413051
12      318665857834031151167461
13      3317044064679887385961981

Test cases for N > 13 are not available because those values have not been found yet. If you manage to find the next term(s) in the sequence, be sure to submit it/them to OEIS!

Rules

• You may choose to take N as a zero-indexed or a one-indexed value.
• It's acceptable for your solution to only work for the values representable within your language's integer range (so up to N = 12 for unsigned 64-bit integers), but your solution must theoretically work for any input with the assumption that your language supports arbitrary-length integers.

---

Background

Any positive even integer x can be written in the form x = d*2^s where d is odd. d and s can be computed by repeatedly dividing n by 2 until the quotient is no longer divisible by 2. d is that final quotient, and s is the number of times 2 divides n.

If a positive integer n is prime, then Fermat's little theorem states:

In any finite field Z/pZ (where p is some prime), the only square roots of 1 are 1 and -1 (or, equivalently, 1 and p-1).

We can use these three facts to prove that one of the following two statements must be true for a prime n (where d*2^s = n-1 and r is some integer in [0, s)):

The Miller-Rabin primality test operates by testing the contrapositive of the above claim: if there is a base a such that both of the above conditions are false, then n is not prime. That base a is called a witness.

Now, testing every base in [1, n) would be prohibitively expensive in computation time for large n. There is a probabilistic variant of the Miller-Rabin test that only tests some randomly-chosen bases in the finite field. However, it was discovered that testing only prime a bases is sufficient, and thus the test can be performed in an efficient and deterministic manner. In fact, not all of the prime bases need to be tested - only a certain number are required, and that number is dependent on the size of the value being tested for primality.

If an insufficient number of prime bases are tested, the test can produce false positives - odd composite integers where the test fails to prove their compositeness. Specifically, if a base a fails to prove the compositeness of an odd composite number, that number is called a strong pseudoprime to base a. This challenge is about finding the odd composite numbers who are strong psuedoprimes to all bases less than or equal to the Nth prime number (which is equivalent to saying that they are strong pseudoprimes to all prime bases less than or equal to the Nth prime number).

Answer 1

C, 349 295 277 267 255 bytes

N,i;__int128 n=2,b,o,l[999];P(m){i<N&&P(m<2?l[i++]=n:n%m?m-1:n++);}main(r){for(P(scanf("%d",&N));r|!b;)for(++n,b=i=N;i--&&b;){for(b=n-1,r=0;~b&1;b/=2)++r;for(o=1;b--;o=o*l[i]%n);for(b=o==1;r--;o=o*o%n)b|=o>n-2;for(o=r=1;++o<n;r&=n%o>0);}printf("%llu",n);}

Takes 1-based input on stdin, e.g.:

echo "1" | ./millerRabin

---

It certainly isn't going to discover any new values in the sequence any time soon, but it gets the job done. UPDATE: now even slower!

• Slightly faster and shorter, with inspiration from Neil A's answer (a^(d*2^r) == (a^d)^(2^r))
• Significantly slower again after realising that all solutions to this challenge will be odd, so there is no need to explicitly enforce that we only check odd numbers.
• Now using GCC's __int128, which is shorter than unsigned long long while also working with larger numbers! Also on little-endian machines, the printf with %llu still works fine.

Less-minified

N,i;
__int128 n=2,b,o,l[999];
P(m){i<N&&P(m<2?l[i++]=n:n%m?m-1:n++);}
main(r){
    for(P(scanf("%d",&N));r|!b;)
        for(++n,b=i=N;i--&&b;){
            for(b=n-1,r=0;~b&1;b/=2)++r;
            for(o=1;b--;o=o*l[i]%n);
            for(b=o==1;r--;o=o*o%n)b|=o>n-2;
            for(o=r=1;++o<n;r&=n%o>0);
        }
    printf("%llu",n);
}

(Outdated) Breakdown

unsigned long long                  // Use the longest type available
n=2,N,b,i,d,p,o,                    // Globals (mostly share names with question)
l[999];                             // Primes to check (limited to 999, but if you
                                    // want it to be "unlimited", change to -1u)
m(){for(o=1;p--;o=o*l[i]%n);}       // Inefficiently calculates (l[i]^p) mod n

// I cannot possibly take credit for this amazing prime finder;
// See https://codegolf.stackexchange.com/a/5818/8927
P(m){i<N&&P(m<2?l[i++]=n:n%m?m-1:n++);}

main(r){
    for(
        P(scanf("%llu",&N));        // Read & calculate desired number of primes
        r|!b;){                     // While we haven't got an answer:
        n=n+1|1;                    // Advance to next odd number
        for(b=1,i=N;i--&&b;){       // For each prime...
            for(d=n-1,r=0;~d&1;d/=2)++r; // Calculate d and r: d*2^r = n-1
            // Check if there exists any r such that a^(d*2^r) == -1 mod n
            // OR a^d == 1 mod n
            m(p=d);
            for(b=o==1;r--;b|=o==n-1)m(p=d<<r);
            // If neither of these exist, we have proven n is not prime,
            // and the outer loop will keep going (!b)
        }
        // Check if n is actually prime;
        // if it is, the outer loop will keep going (r)
        for(i=r=1;++i<n;r&=n%i!=0);
    }
    printf("%llu",n);               // We found a non-prime; print it & stop.
}

---

As mentioned, this uses 1-based input. But for n=0 it produces 9, which follows the related sequence https://oeis.org/A006945. Not any more; now it hangs on 0.

Should work for all n (at least until the output reaches 2^64) but is incredibly slow. I have verified it on n=0, n=1 and (after a lot of waiting), n=2.

Answer 2

Python 2, 633 465 435 292 282 275 256 247 bytes

0-indexed

Question your implementation and try something new

Converting from a function to a program saves some bytes somehow...

If Python 2 gives me a shorter way to do the same thing, I'm going to use Python 2. Division is by default integer, so an easier way to divide by 2, and print doesn't need parentheses.

n=input()
f=i=3
z={2}
a=lambda:min([i%k for k in range(2,i)])
while n:
 if a():z|={i};n-=1
 i+=2
while f:
 i+=2;d,s,f=~-i,0,a()
 while~d&1:d/=2;s+=1
 for y in z:
  x=y**d%i
  if x-1:
   for _ in[]*s:
    x*=x
    if~x%i<1:break
   else:f=1
print i

Try it online!

Python is notoriously slow compared to other languages.

Defines a trial division test for absolute correctness, then repeatedly applies the Miller-Rabin test until a pseudoprime is found.

Try it online!

EDIT: Finally golfed the answer

EDIT: Used min for the trial division primality test and changed it to a lambda. Less efficient, but shorter. Also couldn't help myself and used a couple of bitwise operators (no length difference). In theory it should work (slightly) faster.

EDIT: Thanks @Dave. My editor trolled me. I thought I was using tabs but it was being converted to 4 spaces instead. Also went through pretty much every single Python tip and applied it.

EDIT: Switched to 0-indexing, allows me to save a couple of bytes with generating the primes. Also rethought a couple of comparisons

EDIT: Used a variable to store the result of the tests instead of for/else statements.

EDIT: Moved the lambda inside the function to eliminate the need for a parameter.

EDIT: Converted to a program to save bytes

EDIT: Python 2 saves me bytes! Also I don't have to convert the input to int

Answer 3

Perl + Math::Prime::Util, 81 + 27 = 108 bytes

1 until!is_provable_prime(++$\)&&is_strong_pseudoprime($\,2..nth_prime($_));$_=""

Run with -lpMMath::Prime::Util=:all (27-byte penalty, ouch).

Explanation

It's not just Mathematica that has a builtin for basically anything. Perl has CPAN, one of the first large library repositories, and that has a huge collection of ready-made solutions for tasks like this one. Unfortunately, they aren't imported (or even installed) by default, meaning that it's basically never a good option to use them in code-golf, but when one of them happens to fit the problem perfectly…

We run through consecutive integers until we find one that's not prime, and yet a strong pseudoprime to all integer bases from 2 to the nth prime. The command-line option imports the library that contains the builtin in question, and also sets implicit input (to line-at-a-time; Math::Prime::Util has its own builtin bignum library that doesn't like newlines in its integers). This uses the standard Perl trick of using $\ (output line separator) as a variable in order to reduce awkward parses and allow the output to be generated implicitly.

Note that we need to use is_provable_prime here to request a deterministic, rather than probabilistic, prime test. (Especially given that a probabilistic prime test is likely using Miller-Rabin in the first place, which we can't expect to give reliable results in this case!)

Perl + Math::Prime::Util, 71 + 17 = 88 bytes, in collaboration with @Dada

1until!is_provable_prime(++$\)&is_strong_pseudoprime$\,2..n‌​th_prime$_}{

Run with -lpMntheory=:all (17 byte penalty).

This uses a few Perl golfing tricks that I either didn't know (apparently Math::Prime::Util has an abbreviation!), knew about but didn't think of using (}{ to output $\ implicitly once, rather than "$_$\" implicitly every line), or knew about but somehow managed to misapply (removing parentheses from function calls). Thanks to @Dada for pointing these out to me. Apart from that, it's identical.