math/big: add Baillie-PSW test to (*Int).ProbablyPrime

After x.ProbablyPrime(n) passes the n Miller-Rabin rounds,
add a Baillie-PSW test before declaring x probably prime.

Although the provable error bounds are unchanged, the empirical
error bounds drop dramatically: there are no known inputs
for which Baillie-PSW gives the wrong answer. For example,
before this CL, big.NewInt(443*1327).ProbablyPrime(1) == true.
Now it is (correctly) false.

The new Baillie-PSW test is two pieces: an added Miller-Rabin
round with base 2, and a so-called extra strong Lucas test.
(See the references listed in prime.go for more details.)
The Lucas test takes about 3.5x as long as the Miller-Rabin round,
which is close to theoretical expectations.

name                              time/op
ProbablyPrime/Lucas             2.91ms ± 2%
ProbablyPrime/MillerRabinBase2   850µs ± 1%
ProbablyPrime/n=0               3.75ms ± 3%

The speed of prime testing for a prime input does get slower:

name                  old time/op  new time/op   delta
ProbablyPrime/n=1    849µs ± 1%   4521µs ± 1%  +432.31%   (p=0.000 n=10+9)
ProbablyPrime/n=5   4.31ms ± 3%   7.87ms ± 1%   +82.70%  (p=0.000 n=10+10)
ProbablyPrime/n=10  8.52ms ± 3%  12.28ms ± 1%   +44.11%  (p=0.000 n=10+10)
ProbablyPrime/n=20  16.9ms ± 2%   21.4ms ± 2%   +26.35%   (p=0.000 n=9+10)

However, because the Baillie-PSW test is only added when the old
ProbablyPrime(n) would return true, testing composites runs at
the same speed as before, except in the case where the result
would have been incorrect and is now correct.

In particular, the most important use of this code is for
generating random primes in crypto/rand. That use spends
essentially all its time testing composites, so it is not
slowed down by the new Baillie-PSW check:

name                  old time/op  new time/op   delta
Prime                104ms ±22%    111ms ±16%      ~     (p=0.165 n=10+10)

Thanks to Serhat Şevki Dinçer for CL 20170, which this CL builds on.

Fixes #13229.

Change-Id: Id26dde9b012c7637c85f2e96355d029b6382812a
Reviewed-on: https://go-review.googlesource.com/30770
Run-TryBot: Russ Cox <rsc@golang.org>
TryBot-Result: Gobot Gobot <gobot@golang.org>
Reviewed-by: Robert Griesemer <gri@golang.org>

src/crypto/rand/util_test.go

@@ -7,7 +7,9 @@ package rand_test
 import (
 	"crypto/rand"
 	"math/big"
+	mathrand "math/rand"
 	"testing"
+	"time"
 )
 
 // https://golang.org/issue/6849.
@@ -63,3 +65,10 @@ func TestIntNegativeMaxPanics(t *testing.T) {
 	b := new(big.Int).SetInt64(int64(-1))
 	testIntPanics(t, b)
 }
+
+func BenchmarkPrime(b *testing.B) {
+	r := mathrand.New(mathrand.NewSource(time.Now().UnixNano()))
+	for i := 0; i < b.N; i++ {
+		rand.Prime(r, 1024)
+	}
+}

src/math/big/prime_test.go

@@ -6,7 +6,9 @@ package big
 
 import (
 	"fmt"
+	"strings"
 	"testing"
+	"unicode"
 )
 
 var primes = []string{
@@ -48,28 +50,96 @@ var composites = []string{
 	"6084766654921918907427900243509372380954290099172559290432744450051395395951",
 	"84594350493221918389213352992032324280367711247940675652888030554255915464401",
 	"82793403787388584738507275144194252681",
+
+	// Arnault, "Rabin-Miller Primality Test: Composite Numbers Which Pass It",
+	// Mathematics of Computation, 64(209) (January 1995), pp. 335-361.
+	"1195068768795265792518361315725116351898245581", // strong pseudoprime to prime bases 2 through 29
+	// strong pseudoprime to all prime bases up to 200
+	`
+     80383745745363949125707961434194210813883768828755814583748891752229
+      74273765333652186502336163960045457915042023603208766569966760987284
+       0439654082329287387918508691668573282677617710293896977394701670823
+        0428687109997439976544144845341155872450633409279022275296229414984
+         2306881685404326457534018329786111298960644845216191652872597534901`,
+
+	// Extra-strong Lucas pseudoprimes. https://oeis.org/A217719
+	"989",
+	"3239",
+	"5777",
+	"10877",
+	"27971",
+	"29681",
+	"30739",
+	"31631",
+	"39059",
+	"72389",
+	"73919",
+	"75077",
+	"100127",
+	"113573",
+	"125249",
+	"137549",
+	"137801",
+	"153931",
+	"155819",
+	"161027",
+	"162133",
+	"189419",
+	"218321",
+	"231703",
+	"249331",
+	"370229",
+	"429479",
+	"430127",
+	"459191",
+	"473891",
+	"480689",
+	"600059",
+	"621781",
+	"632249",
+	"635627",
+
+	"3673744903",
+	"3281593591",
+	"2385076987",
+	"2738053141",
+	"2009621503",
+	"1502682721",
+	"255866131",
+	"117987841",
+	"587861",
+
+	"6368689",
+	"8725753",
+	"80579735209",
+	"105919633",
+}
+
+func cutSpace(r rune) rune {
+	if unicode.IsSpace(r) {
+		return -1
+	}
+	return r
 }
 
 func TestProbablyPrime(t *testing.T) {
 	nreps := 20
 	if testing.Short() {
-		nreps = 1
+		nreps = 3
 	}
 	for i, s := range primes {
 		p, _ := new(Int).SetString(s, 10)
-		if !p.ProbablyPrime(nreps) {
+		if !p.ProbablyPrime(nreps) || !p.ProbablyPrime(1) || !p.ProbablyPrime(0) {
 			t.Errorf("#%d prime found to be non-prime (%s)", i, s)
 		}
 	}
 
 	for i, s := range composites {
+		s = strings.Map(cutSpace, s)
 		c, _ := new(Int).SetString(s, 10)
-		if c.ProbablyPrime(nreps) {
+		if c.ProbablyPrime(nreps) || c.ProbablyPrime(1) || c.ProbablyPrime(0) {
 			t.Errorf("#%d composite found to be prime (%s)", i, s)
 		}
-		if testing.Short() {
-			break
-		}
 	}
 
 	// check that ProbablyPrime panics if n <= 0
@@ -77,7 +147,7 @@ func TestProbablyPrime(t *testing.T) {
 	for _, n := range []int{-1, 0, 1} {
 		func() {
 			defer func() {
-				if n <= 0 && recover() == nil {
+				if n < 0 && recover() == nil {
 					t.Fatalf("expected panic from ProbablyPrime(%d)", n)
 				}
 			}()
@@ -90,11 +160,55 @@ func TestProbablyPrime(t *testing.T) {
 
 func BenchmarkProbablyPrime(b *testing.B) {
 	p, _ := new(Int).SetString("203956878356401977405765866929034577280193993314348263094772646453283062722701277632936616063144088173312372882677123879538709400158306567338328279154499698366071906766440037074217117805690872792848149112022286332144876183376326512083574821647933992961249917319836219304274280243803104015000563790123", 10)
-	for _, rep := range []int{1, 5, 10, 20} {
-		b.Run(fmt.Sprintf("Rep=%d", rep), func(b *testing.B) {
+	for _, n := range []int{0, 1, 5, 10, 20} {
+		b.Run(fmt.Sprintf("n=%d", n), func(b *testing.B) {
 			for i := 0; i < b.N; i++ {
-				p.ProbablyPrime(rep)
+				p.ProbablyPrime(n)
 			}
 		})
 	}
+
+	b.Run("Lucas", func(b *testing.B) {
+		for i := 0; i < b.N; i++ {
+			p.abs.probablyPrimeLucas()
+		}
+	})
+	b.Run("MillerRabinBase2", func(b *testing.B) {
+		for i := 0; i < b.N; i++ {
+			p.abs.probablyPrimeMillerRabin(1, true)
+		}
+	})
+}
+
+func TestMillerRabinPseudoprimes(t *testing.T) {
+	testPseudoprimes(t, "probablyPrimeMillerRabin",
+		func(n nat) bool { return n.probablyPrimeMillerRabin(1, true) && !n.probablyPrimeLucas() },
+		// https://oeis.org/A001262
+		[]int{2047, 3277, 4033, 4681, 8321, 15841, 29341, 42799, 49141, 52633, 65281, 74665, 80581, 85489, 88357, 90751})
+}
+
+func TestLucasPseudoprimes(t *testing.T) {
+	testPseudoprimes(t, "probablyPrimeLucas",
+		func(n nat) bool { return n.probablyPrimeLucas() && !n.probablyPrimeMillerRabin(1, true) },
+		// https://oeis.org/A217719
+		[]int{989, 3239, 5777, 10877, 27971, 29681, 30739, 31631, 39059, 72389, 73919, 75077})
+}
+
+func testPseudoprimes(t *testing.T, name string, cond func(nat) bool, want []int) {
+	n := nat{1}
+	for i := 3; i < 100000; i += 2 {
+		n[0] = Word(i)
+		pseudo := cond(n)
+		if pseudo && (len(want) == 0 || i != want[0]) {
+			t.Errorf("%s(%v, base=2) = %v, want false", name, i)
+		} else if !pseudo && len(want) >= 1 && i == want[0] {
+			t.Errorf("%s(%v, base=2) = false, want true", name, i)
+		}
+		if len(want) > 0 && i == want[0] {
+			want = want[1:]
+		}
+	}
+	if len(want) > 0 {
+		t.Fatalf("forgot to test %v", want)
+	}
 }