# Miller-Rabin Strong Pseudoprimes

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

• Sandbox post (now deleted) – Mego May 30 '17 at 22:42
• Is an algorithm that test all odd value from 1 to result for strong pseudo-primality allowed by the rules? – user202729 Jun 3 '17 at 8:41
• @user202729 I don't see why it wouldn't be. What would make you think it is? – Mego Jun 3 '17 at 8:41
• I would suggest making this a fastest-code question because most answers will simply be brute force. – Neil A. Jun 5 '17 at 6:35
• @NeilA. I disagree that this would be better as fastest code. While it's true that answers will almost certainly be brute force (since another algorithm hasn't yet been developed and I don't expect PPCG to do so), code golf is much simpler, has a much lower barrier for entry (since the submitters can score their own solutions), doesn't require me to run and score every solution (and deal with the exorbitant runtimes), and the problem is interesting enough as a golf challenge. – Mego Jun 5 '17 at 9:19

# C, 349295277267 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.

• I make a breakthrough on my solution, and then you just one-up me... Nice! – Neil A. Jun 6 '17 at 20:19
• @NeilA. Sorry! I was playing with shorter int types before you posted your update. I'm sure you'll find 2 bytes somewhere though; this is turning out to be surprisingly competitive considering it's 2 different non-golfing languages :D – Dave Jun 6 '17 at 20:29

# Python 2, 633465435292282275256 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: 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

• +1 for the way you handled a^(d*2^r) mod n! – Dave Jun 5 '17 at 20:50
• Are you aware that you can use single-space (or single-tab) indentation in Python to save… a whole lot of bytes, actually – Dave Jun 5 '17 at 21:45
• @Dave: This is using 1 tab per indentation level – Neil A. Jun 5 '17 at 21:52
• I think your IDE is messing with you and saving spaces while telling you it's using tabs; when I replace them for single spaces, I get a byte count of just 311 bytes! Try it online! – Dave Jun 5 '17 at 21:55
• @Dave: Ok, that's weird, thanks, I'll update the answer. – Neil A. Jun 5 '17 at 22:22

# 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 , 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. • Of course a golf-ish language comes and beats the rest. Well done! – Neil A. Jun 8 '17 at 7:14 • You can use ntheory instead of Math::Prime::Util. Also, }{ instead of ;$_="" should be fine. And you can omit the space after 1 and the parenthesis of a few function calls. Also, & works instead of &&. That should give 88 bytes: perl -Mntheory=:all -lpe '1until!is_provable_prime(++$\)&is_strong_pseudoprime$\,2..nth_prime\$_}{' – Dada Jun 8 '17 at 14:48
• I'd completely forgotten about }{. (Strangely enough, I remembered the parenthesis thing but it's been a while since I golfed in Perl and couldn't remember the rules for leaving it out.) I didn't know about the ntheory abbreviation at all, though. – user62131 Jun 8 '17 at 20:07