1// Copyright (c) 2014 The mathutil Authors. All rights reserved.
2// Use of this source code is governed by a BSD-style
3// license that can be found in the LICENSE file.
4
5package mathutil
6
7import (
8	"bytes"
9	"fmt"
10	"math"
11	"math/big"
12	"math/rand"
13	"os"
14	"path"
15	"runtime"
16	"sort"
17	"strings"
18	"sync"
19	"testing"
20)
21
22func caller(s string, va ...interface{}) {
23	_, fn, fl, _ := runtime.Caller(2)
24	fmt.Fprintf(os.Stderr, "caller: %s:%d: ", path.Base(fn), fl)
25	fmt.Fprintf(os.Stderr, s, va...)
26	fmt.Fprintln(os.Stderr)
27	_, fn, fl, _ = runtime.Caller(1)
28	fmt.Fprintf(os.Stderr, "\tcallee: %s:%d: ", path.Base(fn), fl)
29	fmt.Fprintln(os.Stderr)
30}
31
32func dbg(s string, va ...interface{}) {
33	if s == "" {
34		s = strings.Repeat("%v ", len(va))
35	}
36	_, fn, fl, _ := runtime.Caller(1)
37	fmt.Fprintf(os.Stderr, "dbg %s:%d: ", path.Base(fn), fl)
38	fmt.Fprintf(os.Stderr, s, va...)
39	fmt.Fprintln(os.Stderr)
40}
41
42func TODO(...interface{}) string {
43	_, fn, fl, _ := runtime.Caller(1)
44	return fmt.Sprintf("TODO: %s:%d:\n", path.Base(fn), fl)
45}
46
47func use(...interface{}) {}
48
49// ============================================================================
50
51func r32() *FC32 {
52	r, err := NewFC32(math.MinInt32, math.MaxInt32, true)
53	if err != nil {
54		panic(err)
55	}
56
57	return r
58}
59
60var (
61	r64lo          = big.NewInt(math.MinInt64)
62	r64hi          = big.NewInt(math.MaxInt64)
63	_3             = big.NewInt(3)
64	MinIntM1       = MinInt
65	MaxIntP1       = MaxInt
66	MaxUintP1 uint = MaxUint
67)
68
69func init() {
70	MinIntM1--
71	MaxIntP1++
72	MaxUintP1++
73}
74
75func r64() *FCBig {
76	r, err := NewFCBig(r64lo, r64hi, true)
77	if err != nil {
78		panic(err)
79	}
80
81	return r
82}
83
84func benchmark1eN(b *testing.B, r *FC32) {
85	b.StartTimer()
86	for i := 0; i < b.N; i++ {
87		r.Next()
88	}
89}
90
91func BenchmarkFC1e3(b *testing.B) {
92	b.StopTimer()
93	r, _ := NewFC32(0, 1e3, false)
94	benchmark1eN(b, r)
95}
96
97func BenchmarkFC1e6(b *testing.B) {
98	b.StopTimer()
99	r, _ := NewFC32(0, 1e6, false)
100	benchmark1eN(b, r)
101}
102
103func BenchmarkFC1e9(b *testing.B) {
104	b.StopTimer()
105	r, _ := NewFC32(0, 1e9, false)
106	benchmark1eN(b, r)
107}
108
109func Test0(t *testing.T) {
110	const N = 10000
111	for n := 1; n < N; n++ {
112		lo, hi := 0, n-1
113		period := int64(hi) - int64(lo) + 1
114		r, err := NewFC32(lo, hi, false)
115		if err != nil {
116			t.Fatal(err)
117		}
118		if r.Cycle()-period > period {
119			t.Fatalf("Cycle exceeds 2 * period")
120		}
121	}
122	for n := 1; n < N; n++ {
123		lo, hi := 0, n-1
124		period := int64(hi) - int64(lo) + 1
125		r, err := NewFC32(lo, hi, true)
126		if err != nil {
127			t.Fatal(err)
128		}
129		if r.Cycle()-2*period > period {
130			t.Fatalf("Cycle exceeds 3 * period")
131		}
132	}
133}
134
135func Test1(t *testing.T) {
136	const (
137		N = 360
138		S = 3
139	)
140	for hq := 0; hq <= 1; hq++ {
141		for n := 1; n < N; n++ {
142			for seed := 0; seed < S; seed++ {
143				lo, hi := -n, 2*n
144				period := int64(hi - lo + 1)
145				r, err := NewFC32(lo, hi, hq == 1)
146				if err != nil {
147					t.Fatal(err)
148				}
149				r.Seed(int64(seed))
150				m := map[int]bool{}
151				v := make([]int, period, period)
152				p := make([]int64, period, period)
153				for i := lo; i <= hi; i++ {
154					x := r.Next()
155					p[i-lo] = r.Pos()
156					if x < lo || x > hi {
157						t.Fatal("t1.0")
158					}
159					if m[x] {
160						t.Fatal("t1.1")
161					}
162					m[x] = true
163					v[i-lo] = x
164				}
165				for i := lo; i <= hi; i++ {
166					x := r.Next()
167					if x < lo || x > hi {
168						t.Fatal("t1.2")
169					}
170					if !m[x] {
171						t.Fatal("t1.3")
172					}
173					if x != v[i-lo] {
174						t.Fatal("t1.4")
175					}
176					if r.Pos() != p[i-lo] {
177						t.Fatal("t1.5")
178					}
179					m[x] = false
180				}
181				for i := lo; i <= hi; i++ {
182					r.Seek(p[i-lo] + 1)
183					x := r.Prev()
184					if x < lo || x > hi {
185						t.Fatal("t1.6")
186					}
187					if x != v[i-lo] {
188						t.Fatal("t1.7")
189					}
190				}
191			}
192		}
193	}
194}
195
196func Test2(t *testing.T) {
197	const (
198		N = 370
199		S = 3
200	)
201	for hq := 0; hq <= 1; hq++ {
202		for n := 1; n < N; n++ {
203			for seed := 0; seed < S; seed++ {
204				lo, hi := -n, 2*n
205				period := int64(hi - lo + 1)
206				r, err := NewFC32(lo, hi, hq == 1)
207				if err != nil {
208					t.Fatal(err)
209				}
210				r.Seed(int64(seed))
211				m := map[int]bool{}
212				v := make([]int, period, period)
213				p := make([]int64, period, period)
214				for i := lo; i <= hi; i++ {
215					x := r.Prev()
216					p[i-lo] = r.Pos()
217					if x < lo || x > hi {
218						t.Fatal("t2.0")
219					}
220					if m[x] {
221						t.Fatal("t2.1")
222					}
223					m[x] = true
224					v[i-lo] = x
225				}
226				for i := lo; i <= hi; i++ {
227					x := r.Prev()
228					if x < lo || x > hi {
229						t.Fatal("t2.2")
230					}
231					if !m[x] {
232						t.Fatal("t2.3")
233					}
234					if x != v[i-lo] {
235						t.Fatal("t2.4")
236					}
237					if r.Pos() != p[i-lo] {
238						t.Fatal("t2.5")
239					}
240					m[x] = false
241				}
242				for i := lo; i <= hi; i++ {
243					s := p[i-lo] - 1
244					if s < 0 {
245						s = r.Cycle() - 1
246					}
247					r.Seek(s)
248					x := r.Next()
249					if x < lo || x > hi {
250						t.Fatal("t2.6")
251					}
252					if x != v[i-lo] {
253						t.Fatal("t2.7")
254					}
255				}
256			}
257		}
258	}
259}
260
261func benchmarkBig1eN(b *testing.B, r *FCBig) {
262	b.StartTimer()
263	for i := 0; i < b.N; i++ {
264		r.Next()
265	}
266}
267
268func BenchmarkFCBig1e3(b *testing.B) {
269	b.StopTimer()
270	hi := big.NewInt(0).SetInt64(1e3)
271	r, _ := NewFCBig(big0, hi, false)
272	benchmarkBig1eN(b, r)
273}
274
275func BenchmarkFCBig1e6(b *testing.B) {
276	b.StopTimer()
277	hi := big.NewInt(0).SetInt64(1e6)
278	r, _ := NewFCBig(big0, hi, false)
279	benchmarkBig1eN(b, r)
280}
281
282func BenchmarkFCBig1e9(b *testing.B) {
283	b.StopTimer()
284	hi := big.NewInt(0).SetInt64(1e9)
285	r, _ := NewFCBig(big0, hi, false)
286	benchmarkBig1eN(b, r)
287}
288
289func BenchmarkFCBig1e12(b *testing.B) {
290	b.StopTimer()
291	hi := big.NewInt(0).SetInt64(1e12)
292	r, _ := NewFCBig(big0, hi, false)
293	benchmarkBig1eN(b, r)
294}
295
296func BenchmarkFCBig1e15(b *testing.B) {
297	b.StopTimer()
298	hi := big.NewInt(0).SetInt64(1e15)
299	r, _ := NewFCBig(big0, hi, false)
300	benchmarkBig1eN(b, r)
301}
302
303func BenchmarkFCBig1e18(b *testing.B) {
304	b.StopTimer()
305	hi := big.NewInt(0).SetInt64(1e18)
306	r, _ := NewFCBig(big0, hi, false)
307	benchmarkBig1eN(b, r)
308}
309
310var (
311	big0 = big.NewInt(0)
312)
313
314func TestBig0(t *testing.T) {
315	const N = 7400
316	lo := big.NewInt(0)
317	hi := big.NewInt(0)
318	period := big.NewInt(0)
319	c := big.NewInt(0)
320	for n := int64(1); n < N; n++ {
321		hi.SetInt64(n - 1)
322		period.Set(hi)
323		period.Sub(period, lo)
324		period.Add(period, _1)
325		r, err := NewFCBig(lo, hi, false)
326		if err != nil {
327			t.Fatal(err)
328		}
329		if r.cycle.Cmp(period) < 0 {
330			t.Fatalf("Period exceeds cycle")
331		}
332		c.Set(r.Cycle())
333		c.Sub(c, period)
334		if c.Cmp(period) > 0 {
335			t.Fatalf("Cycle exceeds 2 * period")
336		}
337	}
338	for n := int64(1); n < N; n++ {
339		hi.SetInt64(n - 1)
340		period.Set(hi)
341		period.Sub(period, lo)
342		period.Add(period, _1)
343		r, err := NewFCBig(lo, hi, true)
344		if err != nil {
345			t.Fatal(err)
346		}
347		if r.cycle.Cmp(period) < 0 {
348			t.Fatalf("Period exceeds cycle")
349		}
350		c.Set(r.Cycle())
351		c.Sub(c, period)
352		c.Sub(c, period)
353		if c.Cmp(period) > 0 {
354			t.Fatalf("Cycle exceeds 3 * period")
355		}
356	}
357}
358
359func TestBig1(t *testing.T) {
360	const (
361		N = 120
362		S = 3
363	)
364	lo := big.NewInt(0)
365	hi := big.NewInt(0)
366	seek := big.NewInt(0)
367	for hq := 0; hq <= 1; hq++ {
368		for n := int64(1); n < N; n++ {
369			for seed := 0; seed < S; seed++ {
370				lo64 := -n
371				hi64 := 2 * n
372				lo.SetInt64(lo64)
373				hi.SetInt64(hi64)
374				period := hi64 - lo64 + 1
375				r, err := NewFCBig(lo, hi, hq == 1)
376				if err != nil {
377					t.Fatal(err)
378				}
379				r.Seed(int64(seed))
380				m := map[int64]bool{}
381				v := make([]int64, period, period)
382				p := make([]int64, period, period)
383				for i := lo64; i <= hi64; i++ {
384					x := r.Next().Int64()
385					p[i-lo64] = r.Pos().Int64()
386					if x < lo64 || x > hi64 {
387						t.Fatal("tb1.0")
388					}
389					if m[x] {
390						t.Fatal("tb1.1")
391					}
392					m[x] = true
393					v[i-lo64] = x
394				}
395				for i := lo64; i <= hi64; i++ {
396					x := r.Next().Int64()
397					if x < lo64 || x > hi64 {
398						t.Fatal("tb1.2")
399					}
400					if !m[x] {
401						t.Fatal("tb1.3")
402					}
403					if x != v[i-lo64] {
404						t.Fatal("tb1.4")
405					}
406					if r.Pos().Int64() != p[i-lo64] {
407						t.Fatal("tb1.5")
408					}
409					m[x] = false
410				}
411				for i := lo64; i <= hi64; i++ {
412					r.Seek(seek.SetInt64(p[i-lo64] + 1))
413					x := r.Prev().Int64()
414					if x < lo64 || x > hi64 {
415						t.Fatal("tb1.6")
416					}
417					if x != v[i-lo64] {
418						t.Fatal("tb1.7")
419					}
420				}
421			}
422		}
423	}
424}
425
426func TestBig2(t *testing.T) {
427	const (
428		N = 120
429		S = 3
430	)
431	lo := big.NewInt(0)
432	hi := big.NewInt(0)
433	seek := big.NewInt(0)
434	for hq := 0; hq <= 1; hq++ {
435		for n := int64(1); n < N; n++ {
436			for seed := 0; seed < S; seed++ {
437				lo64, hi64 := -n, 2*n
438				lo.SetInt64(lo64)
439				hi.SetInt64(hi64)
440				period := hi64 - lo64 + 1
441				r, err := NewFCBig(lo, hi, hq == 1)
442				if err != nil {
443					t.Fatal(err)
444				}
445				r.Seed(int64(seed))
446				m := map[int64]bool{}
447				v := make([]int64, period, period)
448				p := make([]int64, period, period)
449				for i := lo64; i <= hi64; i++ {
450					x := r.Prev().Int64()
451					p[i-lo64] = r.Pos().Int64()
452					if x < lo64 || x > hi64 {
453						t.Fatal("tb2.0")
454					}
455					if m[x] {
456						t.Fatal("tb2.1")
457					}
458					m[x] = true
459					v[i-lo64] = x
460				}
461				for i := lo64; i <= hi64; i++ {
462					x := r.Prev().Int64()
463					if x < lo64 || x > hi64 {
464						t.Fatal("tb2.2")
465					}
466					if !m[x] {
467						t.Fatal("tb2.3")
468					}
469					if x != v[i-lo64] {
470						t.Fatal("tb2.4")
471					}
472					if r.Pos().Int64() != p[i-lo64] {
473						t.Fatal("tb2.5")
474					}
475					m[x] = false
476				}
477				for i := lo64; i <= hi64; i++ {
478					s := p[i-lo64] - 1
479					if s < 0 {
480						s = r.Cycle().Int64() - 1
481					}
482					r.Seek(seek.SetInt64(s))
483					x := r.Next().Int64()
484					if x < lo64 || x > hi64 {
485						t.Fatal("tb2.6")
486					}
487					if x != v[i-lo64] {
488						t.Fatal("tb2.7")
489					}
490				}
491			}
492		}
493	}
494}
495
496func TestPermutations(t *testing.T) {
497	data := sort.IntSlice{3, 2, 1}
498	check := [][]int{
499		{1, 2, 3},
500		{1, 3, 2},
501		{2, 1, 3},
502		{2, 3, 1},
503		{3, 1, 2},
504		{3, 2, 1},
505	}
506	i := 0
507	for PermutationFirst(data); ; i++ {
508		if i >= len(check) {
509			t.Fatalf("too much permutations generated: %d > %d", i+1, len(check))
510		}
511
512		for j, v := range check[i] {
513			got := data[j]
514			if got != v {
515				t.Fatalf("permutation %d:\ndata: %v\ncheck: %v\nexpected data[%d] == %d, got %d", i, data, check[i], j, v, got)
516			}
517		}
518
519		if !PermutationNext(data) {
520			if i != len(check)-1 {
521				t.Fatal("permutations generated", i, "expected", len(check))
522			}
523			break
524		}
525	}
526}
527
528func TestIsPrime(t *testing.T) {
529	const p4M = 283146 // # of primes < 4e6
530	n := 0
531	for i := uint32(0); i <= 4e6; i++ {
532		if IsPrime(i) {
533			n++
534		}
535	}
536	t.Log(n)
537	if n != p4M {
538		t.Fatal(n)
539	}
540}
541
542func BenchmarkIsPrime(b *testing.B) {
543	b.StopTimer()
544	n := make([]uint32, b.N)
545	rng := rand.New(rand.NewSource(1))
546	for i := 0; i < b.N; i++ {
547		n[i] = rng.Uint32()
548	}
549	b.StartTimer()
550	for _, n := range n {
551		IsPrime(n)
552	}
553}
554
555func BenchmarkNextPrime(b *testing.B) {
556	b.StopTimer()
557	n := make([]uint32, b.N)
558	rng := rand.New(rand.NewSource(1))
559	for i := 0; i < b.N; i++ {
560		n[i] = rng.Uint32()
561	}
562	b.StartTimer()
563	for _, n := range n {
564		NextPrime(n)
565	}
566}
567
568func BenchmarkIsPrimeUint64(b *testing.B) {
569	const N = 1 << 16
570	b.StopTimer()
571	a := make([]uint64, N)
572	r := r64()
573	for i := range a {
574		a[i] = uint64(r.Next().Int64())
575	}
576	runtime.GC()
577	b.StartTimer()
578	for i := 0; i < b.N; i++ {
579		IsPrimeUint64(a[i&(N-1)])
580	}
581}
582
583func BenchmarkNextPrimeUint64(b *testing.B) {
584	b.StopTimer()
585	n := make([]uint64, b.N)
586	rng := rand.New(rand.NewSource(1))
587	for i := 0; i < b.N; i++ {
588		n[i] = uint64(rng.Int63())
589		if i&1 == 0 {
590			n[i] ^= 1 << 63
591		}
592	}
593	b.StartTimer()
594	for _, n := range n {
595		NextPrimeUint64(n)
596	}
597}
598
599func TestNextPrime(t *testing.T) {
600	const p4M = 283146 // # of primes < 4e6
601	n := 0
602	var p uint32
603	for {
604		p, _ = NextPrime(p)
605		if p >= 4e6 {
606			break
607		}
608		n++
609	}
610	t.Log(n)
611	if n != p4M {
612		t.Fatal(n)
613	}
614}
615
616func TestNextPrime2(t *testing.T) {
617	type data struct {
618		x  uint32
619		y  uint32
620		ok bool
621	}
622	tests := []data{
623		{0, 2, true},
624		{1, 2, true},
625		{2, 3, true},
626		{3, 5, true},
627		{math.MaxUint32, 0, false},
628		{math.MaxUint32 - 1, 0, false},
629		{math.MaxUint32 - 2, 0, false},
630		{math.MaxUint32 - 3, 0, false},
631		{math.MaxUint32 - 4, 0, false},
632		{math.MaxUint32 - 5, math.MaxUint32 - 4, true},
633	}
634
635	for _, test := range tests {
636		y, ok := NextPrime(test.x)
637		if ok != test.ok || ok && y != test.y {
638			t.Fatalf("x %d, got y %d ok %t, expected y %d ok %t", test.x, y, ok, test.y, test.ok)
639		}
640	}
641}
642
643func TestNextPrimeUint64(t *testing.T) {
644	const (
645		lo = 2000000000000000000
646		hi = 2000000000000100000
647		k  = 2346 // PrimePi(hi)-PrimePi(lo)
648	)
649	n := 0
650	p := uint64(lo) - 1
651	var ok bool
652	for {
653		p0 := p
654		p, ok = NextPrimeUint64(p)
655		if !ok {
656			t.Fatal(p0)
657		}
658
659		if p > hi {
660			break
661		}
662
663		n++
664	}
665	if n != k {
666		t.Fatal(n, k)
667	}
668}
669
670func TestISqrt(t *testing.T) {
671	for n := int64(0); n < 5e6; n++ {
672		x := int64(ISqrt(uint32(n)))
673		if x2 := x * x; x2 > n {
674			t.Fatalf("got ISqrt(%d) == %d, too big", n, x)
675		}
676		if x2 := x*x + 2*x + 1; x2 < n {
677			t.Fatalf("got ISqrt(%d) == %d, too low", n, x)
678		}
679	}
680	for n := int64(math.MaxUint32); n > math.MaxUint32-5e6; n-- {
681		x := int64(ISqrt(uint32(n)))
682		if x2 := x * x; x2 > n {
683			t.Fatalf("got ISqrt(%d) == %d, too big", n, x)
684		}
685		if x2 := x*x + 2*x + 1; x2 < n {
686			t.Fatalf("got ISqrt(%d) == %d, too low", n, x)
687		}
688	}
689	rng := rand.New(rand.NewSource(1))
690	for i := 0; i < 5e6; i++ {
691		n := int64(rng.Uint32())
692		x := int64(ISqrt(uint32(n)))
693		if x2 := x * x; x2 > n {
694			t.Fatalf("got ISqrt(%d) == %d, too big", n, x)
695		}
696		if x2 := x*x + 2*x + 1; x2 < n {
697			t.Fatalf("got ISqrt(%d) == %d, too low", n, x)
698		}
699	}
700}
701
702func TestSqrtUint64(t *testing.T) {
703	for n := uint64(0); n < 2e6; n++ {
704		x := SqrtUint64(n)
705		if x > math.MaxUint32 {
706			t.Fatalf("got Sqrt(%d) == %d, too big", n, x)
707		}
708		if x2 := x * x; x2 > n {
709			t.Fatalf("got Sqrt(%d) == %d, too big", n, x)
710		}
711		if x2 := x*x + 2*x + 1; x2 < n {
712			t.Fatalf("got Sqrt(%d) == %d, too low", n, x)
713		}
714	}
715	const H = uint64(18446744056529682436)
716	for n := H; n > H-2e6; n-- {
717		x := SqrtUint64(n)
718		if x > math.MaxUint32 {
719			t.Fatalf("got Sqrt(%d) == %d, too big", n, x)
720		}
721		if x2 := x * x; x2 > n {
722			t.Fatalf("got Sqrt(%d) == %d, too big", n, x)
723		}
724		if x2 := x*x + 2*x + 1; x2 < n {
725			t.Fatalf("got Sqrt(%d) == %d, too low", n, x)
726		}
727	}
728	rng := rand.New(rand.NewSource(1))
729	for i := 0; i < 2e6; i++ {
730		n := uint64(rng.Uint32())<<31 | uint64(rng.Uint32())
731		x := SqrtUint64(n)
732		if x2 := x * x; x2 > n {
733			t.Fatalf("got Sqrt(%d) == %d, too big", n, x)
734		}
735		if x2 := x*x + 2*x + 1; x2 < n {
736			t.Fatalf("got Sqrt(%d) == %d, too low", n, x)
737		}
738	}
739}
740
741func TestSqrtBig(t *testing.T) {
742	const N = 3e4
743	var n, lim, x2 big.Int
744	lim.SetInt64(N)
745	for n.Cmp(&lim) != 0 {
746		x := SqrtBig(&n)
747		x2.Mul(x, x)
748		if x.Cmp(&n) > 0 {
749			t.Fatalf("got sqrt(%s) == %s, too big", &n, x)
750		}
751		x2.Add(&x2, x)
752		x2.Add(&x2, x)
753		x2.Add(&x2, _1)
754		if x2.Cmp(&n) < 0 {
755			t.Fatalf("got sqrt(%s) == %s, too low", &n, x)
756		}
757		n.Add(&n, _1)
758	}
759	rng := rand.New(rand.NewSource(1))
760	var h big.Int
761	h.SetBit(&h, 1e3, 1)
762	for i := 0; i < N; i++ {
763		n.Rand(rng, &h)
764		x := SqrtBig(&n)
765		x2.Mul(x, x)
766		if x.Cmp(&n) > 0 {
767			t.Fatalf("got sqrt(%s) == %s, too big", &n, x)
768		}
769		x2.Add(&x2, x)
770		x2.Add(&x2, x)
771		x2.Add(&x2, _1)
772		if x2.Cmp(&n) < 0 {
773			t.Fatalf("got sqrt(%s) == %s, too low", &n, x)
774		}
775	}
776}
777
778func TestFactorInt(t *testing.T) {
779	chk := func(n uint64, f []FactorTerm) bool {
780		if n < 2 {
781			return len(f) == 0
782		}
783
784		for i := 1; i < len(f); i++ { // verify ordering
785			if t, u := f[i-1], f[i]; t.Prime >= u.Prime {
786				return false
787			}
788		}
789
790		x := uint64(1)
791		for _, v := range f {
792			if p := v.Prime; p < 0 || !IsPrime(uint32(v.Prime)) {
793				return false
794			}
795
796			for i := uint32(0); i < v.Power; i++ {
797				x *= uint64(v.Prime)
798				if x > math.MaxUint32 {
799					return false
800				}
801			}
802		}
803		return x == n
804	}
805
806	for n := uint64(0); n < 3e5; n++ {
807		f := FactorInt(uint32(n))
808		if !chk(n, f) {
809			t.Fatalf("bad FactorInt(%d): %v", n, f)
810		}
811	}
812	for n := uint64(math.MaxUint32); n > math.MaxUint32-12e4; n-- {
813		f := FactorInt(uint32(n))
814		if !chk(n, f) {
815			t.Fatalf("bad FactorInt(%d): %v", n, f)
816		}
817	}
818	rng := rand.New(rand.NewSource(1))
819	for i := 0; i < 13e4; i++ {
820		n := rng.Uint32()
821		f := FactorInt(n)
822		if !chk(uint64(n), f) {
823			t.Fatalf("bad FactorInt(%d): %v", n, f)
824		}
825	}
826}
827
828func TestFactorIntB(t *testing.T) {
829	const N = 3e5 // must be < math.MaxInt32
830	factors := make([][]FactorTerm, N+1)
831	// set up the divisors
832	for prime := uint32(2); prime <= N; prime, _ = NextPrime(prime) {
833		for n := int(prime); n <= N; n += int(prime) {
834			factors[n] = append(factors[n], FactorTerm{prime, 0})
835		}
836	}
837	// set up the powers
838	for n := 2; n <= N; n++ {
839		f := factors[n]
840		m := uint32(n)
841		for i, v := range f {
842			for m%v.Prime == 0 {
843				m /= v.Prime
844				v.Power++
845			}
846			f[i] = v
847		}
848		factors[n] = f
849	}
850	// check equal
851	for n, e := range factors {
852		g := FactorInt(uint32(n))
853		if len(e) != len(g) {
854			t.Fatal(n, "len", g, "!=", e)
855		}
856
857		for i, ev := range e {
858			gv := g[i]
859			if ev.Prime != gv.Prime {
860				t.Fatal(n, "prime", gv, ev)
861			}
862
863			if ev.Power != gv.Power {
864				t.Fatal(n, "power", gv, ev)
865			}
866		}
867	}
868}
869
870func BenchmarkISqrt(b *testing.B) {
871	b.StopTimer()
872	n := make([]uint32, b.N)
873	rng := rand.New(rand.NewSource(1))
874	for i := 0; i < b.N; i++ {
875		n[i] = rng.Uint32()
876	}
877	b.StartTimer()
878	for _, n := range n {
879		ISqrt(n)
880	}
881}
882
883func BenchmarkSqrtUint64(b *testing.B) {
884	b.StopTimer()
885	n := make([]uint64, b.N)
886	rng := rand.New(rand.NewSource(1))
887	for i := 0; i < b.N; i++ {
888		n[i] = uint64(rng.Uint32())<<32 | uint64(rng.Uint32())
889	}
890	b.StartTimer()
891	for _, n := range n {
892		SqrtUint64(n)
893	}
894}
895
896func benchmarkSqrtBig(b *testing.B, bits int) {
897	b.StopTimer()
898	n := make([]*big.Int, b.N)
899	rng := rand.New(rand.NewSource(1))
900	var nn, h big.Int
901	h.SetBit(&h, bits, 1)
902	for i := 0; i < b.N; i++ {
903		n[i] = nn.Rand(rng, &h)
904	}
905	runtime.GC()
906	b.StartTimer()
907	for _, n := range n {
908		SqrtBig(n)
909	}
910}
911
912func BenchmarkSqrtBig2e1e1(b *testing.B) {
913	benchmarkSqrtBig(b, 1e1)
914}
915
916func BenchmarkSqrtBig2e1e2(b *testing.B) {
917	benchmarkSqrtBig(b, 1e2)
918}
919
920func BenchmarkSqrtBig2e1e3(b *testing.B) {
921	benchmarkSqrtBig(b, 1e3)
922}
923
924func BenchmarkSqrtBig2e1e4(b *testing.B) {
925	benchmarkSqrtBig(b, 1e4)
926}
927
928func BenchmarkSqrtBig2e1e5(b *testing.B) {
929	benchmarkSqrtBig(b, 1e5)
930}
931
932func BenchmarkFactorInt(b *testing.B) {
933	b.StopTimer()
934	n := make([]uint32, b.N)
935	rng := rand.New(rand.NewSource(1))
936	for i := 0; i < b.N; i++ {
937		n[i] = rng.Uint32()
938	}
939	b.StartTimer()
940	for _, n := range n {
941		FactorInt(n)
942	}
943}
944
945func TestIsPrimeUint16(t *testing.T) {
946	for n := 0; n <= math.MaxUint16; n++ {
947		if IsPrimeUint16(uint16(n)) != IsPrime(uint32(n)) {
948			t.Fatal(n)
949		}
950	}
951}
952
953func BenchmarkIsPrimeUint16(b *testing.B) {
954	b.StopTimer()
955	n := make([]uint16, b.N)
956	rng := rand.New(rand.NewSource(1))
957	for i := 0; i < b.N; i++ {
958		n[i] = uint16(rng.Uint32())
959	}
960	b.StartTimer()
961	for _, n := range n {
962		IsPrimeUint16(n)
963	}
964}
965
966func TestNextPrimeUint16(t *testing.T) {
967	for n := 0; n <= math.MaxUint16; n++ {
968		p, ok := NextPrimeUint16(uint16(n))
969		p2, ok2 := NextPrime(uint32(n))
970		switch {
971		case ok:
972			if !ok2 || uint32(p) != p2 {
973				t.Fatal(n, p, ok)
974			}
975		case !ok && ok2:
976			if p2 < 65536 {
977				t.Fatal(n, p, ok)
978			}
979		}
980	}
981}
982
983func BenchmarkNextPrimeUint16(b *testing.B) {
984	b.StopTimer()
985	n := make([]uint16, b.N)
986	rng := rand.New(rand.NewSource(1))
987	for i := 0; i < b.N; i++ {
988		n[i] = uint16(rng.Uint32())
989	}
990	b.StartTimer()
991	for _, n := range n {
992		NextPrimeUint16(n)
993	}
994}
995
996/*
997
998From: http://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetKernighan
999
1000Counting bits set, Brian Kernighan's way
1001
1002unsigned int v; // count the number of bits set in v
1003unsigned int c; // c accumulates the total bits set in v
1004for (c = 0; v; c++)
1005{
1006  v &= v - 1; // clear the least significant bit set
1007}
1008
1009Brian Kernighan's method goes through as many iterations as there are set bits.
1010So if we have a 32-bit word with only the high bit set, then it will only go
1011once through the loop.
1012
1013Published in 1988, the C Programming Language 2nd Ed. (by Brian W. Kernighan
1014and Dennis M. Ritchie) mentions this in exercise 2-9. On April 19, 2006 Don
1015Knuth pointed out to me that this method "was first published by Peter Wegner
1016in CACM 3 (1960), 322. (Also discovered independently by Derrick Lehmer and
1017published in 1964 in a book edited by Beckenbach.)"
1018*/
1019func bcnt(v uint64) (c int) {
1020	for ; v != 0; c++ {
1021		v &= v - 1
1022	}
1023	return
1024}
1025
1026func TestPopCount(t *testing.T) {
1027	const N = 4e5
1028	maxUint64 := big.NewInt(0)
1029	maxUint64.SetBit(maxUint64, 64, 1)
1030	maxUint64.Sub(maxUint64, big.NewInt(1))
1031	rng := r64()
1032	for i := 0; i < N; i++ {
1033		n := uint64(rng.Next().Int64())
1034		if g, e := PopCountByte(byte(n)), bcnt(uint64(byte(n))); g != e {
1035			t.Fatal(n, g, e)
1036		}
1037
1038		if g, e := PopCountUint16(uint16(n)), bcnt(uint64(uint16(n))); g != e {
1039			t.Fatal(n, g, e)
1040		}
1041
1042		if g, e := PopCountUint32(uint32(n)), bcnt(uint64(uint32(n))); g != e {
1043			t.Fatal(n, g, e)
1044		}
1045
1046		if g, e := PopCount(int(n)), bcnt(uint64(uint(n))); g != e {
1047			t.Fatal(n, g, e)
1048		}
1049
1050		if g, e := PopCountUint(uint(n)), bcnt(uint64(uint(n))); g != e {
1051			t.Fatal(n, g, e)
1052		}
1053
1054		if g, e := PopCountUint64(n), bcnt(n); g != e {
1055			t.Fatal(n, g, e)
1056		}
1057
1058		if g, e := PopCountUintptr(uintptr(n)), bcnt(uint64(uintptr(n))); g != e {
1059			t.Fatal(n, g, e)
1060		}
1061	}
1062}
1063
1064var gcds = []struct{ a, b, gcd uint64 }{
1065	{8, 12, 4},
1066	{12, 18, 6},
1067	{42, 56, 14},
1068	{54, 24, 6},
1069	{252, 105, 21},
1070	{1989, 867, 51},
1071	{1071, 462, 21},
1072	{2 * 3 * 5 * 7 * 11, 5 * 7 * 11 * 13 * 17, 5 * 7 * 11},
1073	{2 * 3 * 5 * 7 * 7 * 11, 5 * 7 * 7 * 11 * 13 * 17, 5 * 7 * 7 * 11},
1074	{2 * 3 * 5 * 7 * 7 * 11, 5 * 7 * 7 * 13 * 17, 5 * 7 * 7},
1075	{2 * 3 * 5 * 7 * 11, 13 * 17 * 19, 1},
1076}
1077
1078func TestGCD(t *testing.T) {
1079	for i, v := range gcds {
1080		if v.a <= math.MaxUint16 && v.b <= math.MaxUint16 {
1081			if g, e := uint64(GCDUint16(uint16(v.a), uint16(v.b))), v.gcd; g != e {
1082				t.Errorf("%d: got gcd(%d, %d) %d, exp %d", i, v.a, v.b, g, e)
1083			}
1084			if g, e := uint64(GCDUint16(uint16(v.b), uint16(v.a))), v.gcd; g != e {
1085				t.Errorf("%d: got gcd(%d, %d) %d, exp %d", i, v.b, v.a, g, e)
1086			}
1087		}
1088		if v.a <= math.MaxUint32 && v.b <= math.MaxUint32 {
1089			if g, e := uint64(GCDUint32(uint32(v.a), uint32(v.b))), v.gcd; g != e {
1090				t.Errorf("%d: got gcd(%d, %d) %d, exp %d", i, v.a, v.b, g, e)
1091			}
1092			if g, e := uint64(GCDUint32(uint32(v.b), uint32(v.a))), v.gcd; g != e {
1093				t.Errorf("%d: got gcd(%d, %d) %d, exp %d", i, v.b, v.a, g, e)
1094			}
1095		}
1096		if g, e := GCDUint64(v.a, v.b), v.gcd; g != e {
1097			t.Errorf("%d: got gcd(%d, %d) %d, exp %d", i, v.a, v.b, g, e)
1098		}
1099		if g, e := GCDUint64(v.b, v.a), v.gcd; g != e {
1100			t.Errorf("%d: got gcd(%d, %d) %d, exp %d", i, v.b, v.a, g, e)
1101		}
1102	}
1103}
1104
1105func lg2(n uint64) (lg int) {
1106	if n == 0 {
1107		return -1
1108	}
1109
1110	for n >>= 1; n != 0; n >>= 1 {
1111		lg++
1112	}
1113	return
1114}
1115
1116func TestLog2(t *testing.T) {
1117	if g, e := Log2Byte(0), -1; g != e {
1118		t.Error(g, e)
1119	}
1120	if g, e := Log2Uint16(0), -1; g != e {
1121		t.Error(g, e)
1122	}
1123	if g, e := Log2Uint32(0), -1; g != e {
1124		t.Error(g, e)
1125	}
1126	if g, e := Log2Uint64(0), -1; g != e {
1127		t.Error(g, e)
1128	}
1129	const N = 1e6
1130	rng := r64()
1131	for i := 0; i < N; i++ {
1132		n := uint64(rng.Next().Int64())
1133		if g, e := Log2Uint64(n), lg2(n); g != e {
1134			t.Fatalf("%b %d %d", n, g, e)
1135		}
1136		if g, e := Log2Uint32(uint32(n)), lg2(n&0xffffffff); g != e {
1137			t.Fatalf("%b %d %d", n, g, e)
1138		}
1139		if g, e := Log2Uint16(uint16(n)), lg2(n&0xffff); g != e {
1140			t.Fatalf("%b %d %d", n, g, e)
1141		}
1142		if g, e := Log2Byte(byte(n)), lg2(n&0xff); g != e {
1143			t.Fatalf("%b %d %d", n, g, e)
1144		}
1145	}
1146}
1147
1148func TestBitLen(t *testing.T) {
1149	if g, e := BitLenByte(0), 0; g != e {
1150		t.Error(g, e)
1151	}
1152	if g, e := BitLenUint16(0), 0; g != e {
1153		t.Error(g, e)
1154	}
1155	if g, e := BitLenUint32(0), 0; g != e {
1156		t.Error(g, e)
1157	}
1158	if g, e := BitLenUint64(0), 0; g != e {
1159		t.Error(g, e)
1160	}
1161	if g, e := BitLenUintptr(0), 0; g != e {
1162		t.Error(g, e)
1163	}
1164	const N = 1e6
1165	rng := r64()
1166	for i := 0; i < N; i++ {
1167		n := uint64(rng.Next().Int64())
1168		if g, e := BitLenUintptr(uintptr(n)), lg2(uint64(uintptr(n)))+1; g != e {
1169			t.Fatalf("%b %d %d", n, g, e)
1170		}
1171		if g, e := BitLenUint64(n), lg2(n)+1; g != e {
1172			t.Fatalf("%b %d %d", n, g, e)
1173		}
1174		if g, e := BitLenUint32(uint32(n)), lg2(n&0xffffffff)+1; g != e {
1175			t.Fatalf("%b %d %d", n, g, e)
1176		}
1177		if g, e := BitLen(int(n)), lg2(uint64(uint(n)))+1; g != e {
1178			t.Fatalf("%b %d %d", n, g, e)
1179		}
1180		if g, e := BitLenUint(uint(n)), lg2(uint64(uint(n)))+1; g != e {
1181			t.Fatalf("%b %d %d", n, g, e)
1182		}
1183		if g, e := BitLenUint16(uint16(n)), lg2(n&0xffff)+1; g != e {
1184			t.Fatalf("%b %d %d", n, g, e)
1185		}
1186		if g, e := BitLenByte(byte(n)), lg2(n&0xff)+1; g != e {
1187			t.Fatalf("%b %d %d", n, g, e)
1188		}
1189	}
1190}
1191
1192func BenchmarkGCDByte(b *testing.B) {
1193	const N = 1 << 16
1194	type t byte
1195	type u struct{ a, b t }
1196	b.StopTimer()
1197	rng := r32()
1198	a := make([]u, N)
1199	for i := range a {
1200		a[i] = u{t(rng.Next()), t(rng.Next())}
1201	}
1202	b.StartTimer()
1203	for i := 0; i < b.N; i++ {
1204		v := a[i&(N-1)]
1205		GCDByte(byte(v.a), byte(v.b))
1206	}
1207}
1208
1209func BenchmarkGCDUint16(b *testing.B) {
1210	const N = 1 << 16
1211	type t uint16
1212	type u struct{ a, b t }
1213	b.StopTimer()
1214	rng := r32()
1215	a := make([]u, N)
1216	for i := range a {
1217		a[i] = u{t(rng.Next()), t(rng.Next())}
1218	}
1219	b.StartTimer()
1220	for i := 0; i < b.N; i++ {
1221		v := a[i&(N-1)]
1222		GCDUint16(uint16(v.a), uint16(v.b))
1223	}
1224}
1225
1226func BenchmarkGCDUint32(b *testing.B) {
1227	const N = 1 << 16
1228	type t uint32
1229	type u struct{ a, b t }
1230	b.StopTimer()
1231	rng := r32()
1232	a := make([]u, N)
1233	for i := range a {
1234		a[i] = u{t(rng.Next()), t(rng.Next())}
1235	}
1236	b.StartTimer()
1237	for i := 0; i < b.N; i++ {
1238		v := a[i&(N-1)]
1239		GCDUint32(uint32(v.a), uint32(v.b))
1240	}
1241}
1242
1243func BenchmarkGCDUint64(b *testing.B) {
1244	const N = 1 << 16
1245	type t uint64
1246	type u struct{ a, b t }
1247	b.StopTimer()
1248	rng := r64()
1249	a := make([]u, N)
1250	for i := range a {
1251		a[i] = u{t(rng.Next().Int64()), t(rng.Next().Int64())}
1252	}
1253	b.StartTimer()
1254	for i := 0; i < b.N; i++ {
1255		v := a[i&(N-1)]
1256		GCDUint64(uint64(v.a), uint64(v.b))
1257	}
1258}
1259
1260func BenchmarkLog2Byte(b *testing.B) {
1261	const N = 1 << 16
1262	type t byte
1263	b.StopTimer()
1264	rng := r32()
1265	a := make([]t, N)
1266	for i := range a {
1267		a[i] = t(rng.Next())
1268	}
1269	b.StartTimer()
1270	for i := 0; i < b.N; i++ {
1271		Log2Byte(byte(a[i&(N-1)]))
1272	}
1273}
1274
1275func BenchmarkLog2Uint16(b *testing.B) {
1276	const N = 1 << 16
1277	type t uint16
1278	b.StopTimer()
1279	rng := r32()
1280	a := make([]t, N)
1281	for i := range a {
1282		a[i] = t(rng.Next())
1283	}
1284	b.StartTimer()
1285	for i := 0; i < b.N; i++ {
1286		Log2Uint16(uint16(a[i&(N-1)]))
1287	}
1288}
1289
1290func BenchmarkLog2Uint32(b *testing.B) {
1291	const N = 1 << 16
1292	type t uint32
1293	b.StopTimer()
1294	rng := r32()
1295	a := make([]t, N)
1296	for i := range a {
1297		a[i] = t(rng.Next())
1298	}
1299	b.StartTimer()
1300	for i := 0; i < b.N; i++ {
1301		Log2Uint32(uint32(a[i&(N-1)]))
1302	}
1303}
1304
1305func BenchmarkLog2Uint64(b *testing.B) {
1306	const N = 1 << 16
1307	type t uint64
1308	b.StopTimer()
1309	rng := r64()
1310	a := make([]t, N)
1311	for i := range a {
1312		a[i] = t(rng.Next().Int64())
1313	}
1314	b.StartTimer()
1315	for i := 0; i < b.N; i++ {
1316		Log2Uint64(uint64(a[i&(N-1)]))
1317	}
1318}
1319func BenchmarkBitLenByte(b *testing.B) {
1320	const N = 1 << 16
1321	type t byte
1322	b.StopTimer()
1323	rng := r32()
1324	a := make([]t, N)
1325	for i := range a {
1326		a[i] = t(rng.Next())
1327	}
1328	b.StartTimer()
1329	for i := 0; i < b.N; i++ {
1330		BitLenByte(byte(a[i&(N-1)]))
1331	}
1332}
1333
1334func BenchmarkBitLenUint16(b *testing.B) {
1335	const N = 1 << 16
1336	type t uint16
1337	b.StopTimer()
1338	rng := r32()
1339	a := make([]t, N)
1340	for i := range a {
1341		a[i] = t(rng.Next())
1342	}
1343	b.StartTimer()
1344	for i := 0; i < b.N; i++ {
1345		BitLenUint16(uint16(a[i&(N-1)]))
1346	}
1347}
1348
1349func BenchmarkBitLenUint32(b *testing.B) {
1350	const N = 1 << 16
1351	type t uint32
1352	b.StopTimer()
1353	rng := r32()
1354	a := make([]t, N)
1355	for i := range a {
1356		a[i] = t(rng.Next())
1357	}
1358	b.StartTimer()
1359	for i := 0; i < b.N; i++ {
1360		BitLenUint32(uint32(a[i&(N-1)]))
1361	}
1362}
1363
1364func BenchmarkBitLen(b *testing.B) {
1365	const N = 1 << 16
1366	type t int
1367	b.StopTimer()
1368	rng := r64()
1369	a := make([]t, N)
1370	for i := range a {
1371		a[i] = t(rng.Next().Int64())
1372	}
1373	b.StartTimer()
1374	for i := 0; i < b.N; i++ {
1375		BitLen(int(a[i&(N-1)]))
1376	}
1377}
1378
1379func BenchmarkBitLenUint(b *testing.B) {
1380	const N = 1 << 16
1381	type t uint
1382	b.StopTimer()
1383	rng := r64()
1384	a := make([]t, N)
1385	for i := range a {
1386		a[i] = t(rng.Next().Int64())
1387	}
1388	b.StartTimer()
1389	for i := 0; i < b.N; i++ {
1390		BitLenUint(uint(a[i&(N-1)]))
1391	}
1392}
1393
1394func BenchmarkBitLenUintptr(b *testing.B) {
1395	const N = 1 << 16
1396	type t uintptr
1397	b.StopTimer()
1398	rng := r64()
1399	a := make([]t, N)
1400	for i := range a {
1401		a[i] = t(rng.Next().Int64())
1402	}
1403	b.StartTimer()
1404	for i := 0; i < b.N; i++ {
1405		BitLenUintptr(uintptr(a[i&(N-1)]))
1406	}
1407}
1408
1409func BenchmarkBitLenUint64(b *testing.B) {
1410	const N = 1 << 16
1411	type t uint64
1412	b.StopTimer()
1413	rng := r64()
1414	a := make([]t, N)
1415	for i := range a {
1416		a[i] = t(rng.Next().Int64())
1417	}
1418	b.StartTimer()
1419	for i := 0; i < b.N; i++ {
1420		BitLenUint64(uint64(a[i&(N-1)]))
1421	}
1422}
1423
1424func BenchmarkPopCountByte(b *testing.B) {
1425	const N = 1 << 16
1426	type t byte
1427	b.StopTimer()
1428	rng := r32()
1429	a := make([]t, N)
1430	for i := range a {
1431		a[i] = t(rng.Next())
1432	}
1433	b.StartTimer()
1434	for i := 0; i < b.N; i++ {
1435		PopCountByte(byte(a[i&(N-1)]))
1436	}
1437}
1438
1439func BenchmarkPopCountUint16(b *testing.B) {
1440	const N = 1 << 16
1441	type t uint16
1442	b.StopTimer()
1443	rng := r32()
1444	a := make([]t, N)
1445	for i := range a {
1446		a[i] = t(rng.Next())
1447	}
1448	b.StartTimer()
1449	for i := 0; i < b.N; i++ {
1450		PopCountUint16(uint16(a[i&(N-1)]))
1451	}
1452}
1453
1454func BenchmarkPopCountUint32(b *testing.B) {
1455	const N = 1 << 16
1456	type t uint32
1457	b.StopTimer()
1458	rng := r32()
1459	a := make([]t, N)
1460	for i := range a {
1461		a[i] = t(rng.Next())
1462	}
1463	b.StartTimer()
1464	for i := 0; i < b.N; i++ {
1465		PopCountUint32(uint32(a[i&(N-1)]))
1466	}
1467}
1468
1469func BenchmarkPopCount(b *testing.B) {
1470	const N = 1 << 16
1471	type t int
1472	b.StopTimer()
1473	rng := r64()
1474	a := make([]t, N)
1475	for i := range a {
1476		a[i] = t(rng.Next().Int64())
1477	}
1478	b.StartTimer()
1479	for i := 0; i < b.N; i++ {
1480		PopCount(int(a[i&(N-1)]))
1481	}
1482}
1483
1484func BenchmarkPopCountUint(b *testing.B) {
1485	const N = 1 << 16
1486	type t uint
1487	b.StopTimer()
1488	rng := r64()
1489	a := make([]t, N)
1490	for i := range a {
1491		a[i] = t(rng.Next().Int64())
1492	}
1493	b.StartTimer()
1494	for i := 0; i < b.N; i++ {
1495		PopCountUint(uint(a[i&(N-1)]))
1496	}
1497}
1498
1499func BenchmarkPopCountUintptr(b *testing.B) {
1500	const N = 1 << 16
1501	type t uintptr
1502	b.StopTimer()
1503	rng := r64()
1504	a := make([]t, N)
1505	for i := range a {
1506		a[i] = t(rng.Next().Int64())
1507	}
1508	b.StartTimer()
1509	for i := 0; i < b.N; i++ {
1510		PopCountUintptr(uintptr(a[i&(N-1)]))
1511	}
1512}
1513
1514func BenchmarkPopCountUint64(b *testing.B) {
1515	const N = 1 << 16
1516	type t uint64
1517	b.StopTimer()
1518	rng := r64()
1519	a := make([]t, N)
1520	for i := range a {
1521		a[i] = t(rng.Next().Int64())
1522	}
1523	b.StartTimer()
1524	for i := 0; i < b.N; i++ {
1525		PopCountUint64(uint64(a[i&(N-1)]))
1526	}
1527}
1528
1529func TestUintptrBits(t *testing.T) {
1530	switch g := UintptrBits(); g {
1531	case 32, 64:
1532		// ok
1533		t.Log(g)
1534	default:
1535		t.Fatalf("got %d, expected 32 or 64", g)
1536	}
1537}
1538
1539func BenchmarkUintptrBits(b *testing.B) {
1540	for i := 0; i < b.N; i++ {
1541		UintptrBits()
1542	}
1543}
1544
1545func TestModPowByte(t *testing.T) {
1546	data := []struct{ b, e, m, r byte }{
1547		{0, 1, 1, 0},
1548		{0, 2, 1, 0},
1549		{0, 3, 1, 0},
1550
1551		{1, 0, 1, 0},
1552		{1, 1, 1, 0},
1553		{1, 2, 1, 0},
1554		{1, 3, 1, 0},
1555
1556		{2, 0, 1, 0},
1557		{2, 1, 1, 0},
1558		{2, 2, 1, 0},
1559		{2, 3, 1, 0},
1560
1561		{2, 11, 23, 1}, // 23|M11
1562		{2, 11, 89, 1}, // 89|M11
1563		{2, 23, 47, 1}, // 47|M23
1564		{5, 3, 13, 8},
1565	}
1566
1567	for _, v := range data {
1568		if g, e := ModPowByte(v.b, v.e, v.m), v.r; g != e {
1569			t.Errorf("b %d e %d m %d: got %d, exp %d", v.b, v.e, v.m, g, e)
1570		}
1571	}
1572}
1573
1574func TestModPowUint16(t *testing.T) {
1575	data := []struct{ b, e, m, r uint16 }{
1576		{0, 1, 1, 0},
1577		{0, 2, 1, 0},
1578		{0, 3, 1, 0},
1579
1580		{1, 0, 1, 0},
1581		{1, 1, 1, 0},
1582		{1, 2, 1, 0},
1583		{1, 3, 1, 0},
1584
1585		{2, 0, 1, 0},
1586		{2, 1, 1, 0},
1587		{2, 2, 1, 0},
1588		{2, 3, 1, 0},
1589
1590		{2, 11, 23, 1},     // 23|M11
1591		{2, 11, 89, 1},     // 89|M11
1592		{2, 23, 47, 1},     // 47|M23
1593		{2, 929, 13007, 1}, // 13007|M929
1594		{4, 13, 497, 445},
1595		{5, 3, 13, 8},
1596	}
1597
1598	for _, v := range data {
1599		if g, e := ModPowUint16(v.b, v.e, v.m), v.r; g != e {
1600			t.Errorf("b %d e %d m %d: got %d, exp %d", v.b, v.e, v.m, g, e)
1601		}
1602	}
1603}
1604
1605func TestModPowUint32(t *testing.T) {
1606	data := []struct{ b, e, m, r uint32 }{
1607		{0, 1, 1, 0},
1608		{0, 2, 1, 0},
1609		{0, 3, 1, 0},
1610
1611		{1, 0, 1, 0},
1612		{1, 1, 1, 0},
1613		{1, 2, 1, 0},
1614		{1, 3, 1, 0},
1615
1616		{2, 0, 1, 0},
1617		{2, 1, 1, 0},
1618		{2, 2, 1, 0},
1619		{2, 3, 1, 0},
1620
1621		{2, 23, 47, 1},        // 47|M23
1622		{2, 67, 193707721, 1}, // 193707721|M67
1623		{2, 929, 13007, 1},    // 13007|M929
1624		{4, 13, 497, 445},
1625		{5, 3, 13, 8},
1626		{2, 500471, 264248689, 1},
1627		{2, 1000249, 112027889, 1},
1628		{2, 2000633, 252079759, 1},
1629		{2, 3000743, 222054983, 1},
1630		{2, 4000741, 1920355681, 1},
1631		{2, 5000551, 330036367, 1},
1632		{2, 6000479, 1020081431, 1},
1633		{2, 7000619, 840074281, 1},
1634		{2, 8000401, 624031279, 1},
1635		{2, 9000743, 378031207, 1},
1636		{2, 10000961, 380036519, 1},
1637		{2, 20000723, 40001447, 1},
1638	}
1639
1640	for _, v := range data {
1641		if g, e := ModPowUint32(v.b, v.e, v.m), v.r; g != e {
1642			t.Errorf("b %d e %d m %d: got %d, exp %d", v.b, v.e, v.m, g, e)
1643		}
1644	}
1645}
1646
1647func TestModPowUint64(t *testing.T) {
1648	data := []struct{ b, e, m, r uint64 }{
1649		{0, 1, 1, 0},
1650		{0, 2, 1, 0},
1651		{0, 3, 1, 0},
1652
1653		{1, 0, 1, 0},
1654		{1, 1, 1, 0},
1655		{1, 2, 1, 0},
1656		{1, 3, 1, 0},
1657
1658		{2, 0, 1, 0},
1659		{2, 1, 1, 0},
1660		{2, 2, 1, 0},
1661		{2, 3, 1, 0},
1662
1663		{2, 23, 47, 1},        // 47|M23
1664		{2, 67, 193707721, 1}, // 193707721|M67
1665		{2, 929, 13007, 1},    // 13007|M929
1666		{4, 13, 497, 445},
1667		{5, 3, 13, 8},
1668		{2, 500471, 264248689, 1}, // m|Me ...
1669		{2, 1000249, 112027889, 1},
1670		{2, 2000633, 252079759, 1},
1671		{2, 3000743, 222054983, 1},
1672		{2, 4000741, 1920355681, 1},
1673		{2, 5000551, 330036367, 1},
1674		{2, 6000479, 1020081431, 1},
1675		{2, 7000619, 840074281, 1},
1676		{2, 8000401, 624031279, 1},
1677		{2, 9000743, 378031207, 1},
1678		{2, 10000961, 380036519, 1},
1679		{2, 20000723, 40001447, 1},
1680		{2, 1000099, 1872347344039, 1},
1681
1682		{9223372036854775919, 9223372036854776030, 9223372036854776141, 7865333882915297658},
1683	}
1684
1685	for _, v := range data {
1686		if g, e := ModPowUint64(v.b, v.e, v.m), v.r; g != e {
1687			t.Errorf("b %d e %d m %d: got %d, exp %d", v.b, v.e, v.m, g, e)
1688		}
1689	}
1690}
1691
1692func TestModPowBigInt(t *testing.T) {
1693	data := []struct {
1694		b, e int64
1695		m    interface{}
1696		r    int64
1697	}{
1698		{0, 1, 1, 0},
1699		{0, 2, 1, 0},
1700		{0, 3, 1, 0},
1701
1702		{1, 0, 1, 0},
1703		{1, 1, 1, 0},
1704		{1, 2, 1, 0},
1705		{1, 3, 1, 0},
1706
1707		{2, 0, 1, 0},
1708		{2, 1, 1, 0},
1709		{2, 2, 1, 0},
1710		{2, 3, 1, 0},
1711
1712		{2, 23, 47, 1},        // 47|M23
1713		{2, 67, 193707721, 1}, // 193707721|M67
1714		{2, 929, 13007, 1},    // 13007|M929
1715		{4, 13, 497, 445},
1716		{5, 3, 13, 8},
1717		{2, 500471, 264248689, 1}, // m|Me ...
1718		{2, 1000249, 112027889, 1},
1719		{2, 2000633, 252079759, 1},
1720		{2, 3000743, 222054983, 1},
1721		{2, 4000741, 1920355681, 1},
1722		{2, 5000551, 330036367, 1},
1723		{2, 6000479, 1020081431, 1},
1724		{2, 7000619, 840074281, 1},
1725		{2, 8000401, 624031279, 1},
1726		{2, 9000743, 378031207, 1},
1727		{2, 10000961, 380036519, 1},
1728		{2, 20000723, 40001447, 1},
1729		{2, 100279, "11502865265922183403581252152383", 1},
1730
1731		{2, 7293457, "533975545077050000610542659519277030089249998649", 1},
1732	}
1733
1734	for _, v := range data {
1735		var m big.Int
1736		switch x := v.m.(type) {
1737		case int:
1738			m.SetInt64(int64(x))
1739		case string:
1740			m.SetString(x, 10)
1741		}
1742		b, e, r := big.NewInt(v.b), big.NewInt(v.e), big.NewInt(v.r)
1743		if g, e := ModPowBigInt(b, e, &m), r; g.Cmp(e) != 0 {
1744			t.Errorf("b %s e %s m %v: got %s, exp %s", b, e, m, g, e)
1745		}
1746	}
1747
1748	s := func(n string) *big.Int {
1749		i, ok := big.NewInt(0).SetString(n, 10)
1750		if !ok {
1751			t.Fatal(ok)
1752		}
1753
1754		return i
1755	}
1756
1757	if g, e := ModPowBigInt(
1758		s("36893488147419103343"), s("36893488147419103454"), s("36893488147419103565")), s("34853007610367449339"); g.Cmp(e) != 0 {
1759		t.Fatal(g, e)
1760	}
1761}
1762
1763func BenchmarkModPowByte(b *testing.B) {
1764	const N = 1 << 16
1765	b.StopTimer()
1766	type t struct{ b, e, m byte }
1767	a := make([]t, N)
1768	r := r32()
1769	for i := range a {
1770		a[i] = t{
1771			byte(r.Next() | 2),
1772			byte(r.Next() | 2),
1773			byte(r.Next() | 2),
1774		}
1775	}
1776	runtime.GC()
1777	b.StartTimer()
1778	for i := 0; i < b.N; i++ {
1779		v := a[i&(N-1)]
1780		ModPowByte(v.b, v.e, v.m)
1781	}
1782}
1783
1784func BenchmarkModPowUint16(b *testing.B) {
1785	const N = 1 << 16
1786	b.StopTimer()
1787	type t struct{ b, e, m uint16 }
1788	a := make([]t, N)
1789	r := r32()
1790	for i := range a {
1791		a[i] = t{
1792			uint16(r.Next() | 2),
1793			uint16(r.Next() | 2),
1794			uint16(r.Next() | 2),
1795		}
1796	}
1797	runtime.GC()
1798	b.StartTimer()
1799	for i := 0; i < b.N; i++ {
1800		v := a[i&(N-1)]
1801		ModPowUint16(v.b, v.e, v.m)
1802	}
1803}
1804
1805func BenchmarkModPowUint32(b *testing.B) {
1806	const N = 1 << 16
1807	b.StopTimer()
1808	type t struct{ b, e, m uint32 }
1809	a := make([]t, N)
1810	r := r32()
1811	for i := range a {
1812		a[i] = t{
1813			uint32(r.Next() | 2),
1814			uint32(r.Next() | 2),
1815			uint32(r.Next() | 2),
1816		}
1817	}
1818	runtime.GC()
1819	b.StartTimer()
1820	for i := 0; i < b.N; i++ {
1821		v := a[i&(N-1)]
1822		ModPowUint32(v.b, v.e, v.m)
1823	}
1824}
1825
1826func BenchmarkModPowUint64(b *testing.B) {
1827	const N = 1 << 16
1828	b.StopTimer()
1829	type t struct{ b, e, m uint64 }
1830	a := make([]t, N)
1831	r := r64()
1832	for i := range a {
1833		a[i] = t{
1834			uint64(r.Next().Int64() | 2),
1835			uint64(r.Next().Int64() | 2),
1836			uint64(r.Next().Int64() | 2),
1837		}
1838	}
1839	runtime.GC()
1840	b.StartTimer()
1841	for i := 0; i < b.N; i++ {
1842		v := a[i&(N-1)]
1843		ModPowUint64(v.b, v.e, v.m)
1844	}
1845}
1846
1847func BenchmarkModPowBigInt(b *testing.B) {
1848	const N = 1 << 16
1849	b.StopTimer()
1850	type t struct{ b, e, m *big.Int }
1851	a := make([]t, N)
1852	mx := big.NewInt(math.MaxInt64)
1853	mx.Mul(mx, mx)
1854	r, err := NewFCBig(big.NewInt(2), mx, true)
1855	if err != nil {
1856		b.Fatal(err)
1857	}
1858	for i := range a {
1859		a[i] = t{
1860			r.Next(),
1861			r.Next(),
1862			r.Next(),
1863		}
1864	}
1865	runtime.GC()
1866	b.StartTimer()
1867	for i := 0; i < b.N; i++ {
1868		v := a[i&(N-1)]
1869		ModPowBigInt(v.b, v.e, v.m)
1870	}
1871}
1872
1873func TestAdd128(t *testing.T) {
1874	const N = 1e5
1875	r := r64()
1876	var mm big.Int
1877	for i := 0; i < N; i++ {
1878		a, b := uint64(r.Next().Int64()), uint64(r.Next().Int64())
1879		aa, bb := big.NewInt(0).SetUint64(a), big.NewInt(0).SetUint64(b)
1880		mhi, mlo := AddUint128_64(a, b)
1881		m := big.NewInt(0).SetUint64(mhi)
1882		m.Lsh(m, 64)
1883		m.Add(m, big.NewInt(0).SetUint64(mlo))
1884		mm.Add(aa, bb)
1885		if m.Cmp(&mm) != 0 {
1886			t.Fatalf("%d\na %40d\nb %40d\ng %40s %032x\ne %40s %032x", i, a, b, m, m, &mm, &mm)
1887		}
1888	}
1889}
1890
1891func TestMul128(t *testing.T) {
1892	const N = 1e5
1893	r := r64()
1894	var mm big.Int
1895	f := func(a, b uint64) {
1896		aa, bb := big.NewInt(0).SetUint64(a), big.NewInt(0).SetUint64(b)
1897		mhi, mlo := MulUint128_64(a, b)
1898		m := big.NewInt(0).SetUint64(mhi)
1899		m.Lsh(m, 64)
1900		m.Add(m, big.NewInt(0).SetUint64(mlo))
1901		mm.Mul(aa, bb)
1902		if m.Cmp(&mm) != 0 {
1903			t.Fatalf("\na %40d\nb %40d\ng %40s %032x\ne %40s %032x", a, b, m, m, &mm, &mm)
1904		}
1905	}
1906	for i := 0; i < N; i++ {
1907		f(uint64(r.Next().Int64()), uint64(r.Next().Int64()))
1908	}
1909	for x := 0; x <= 1<<9; x++ {
1910		for y := 0; y <= 1<<9; y++ {
1911			f(math.MaxUint64-uint64(x), math.MaxUint64-uint64(y))
1912		}
1913	}
1914}
1915
1916func BenchmarkMul128(b *testing.B) {
1917	const N = 1 << 16
1918	b.StopTimer()
1919	type t struct{ a, b uint64 }
1920	a := make([]t, N)
1921	r := r64()
1922	for i := range a {
1923		a[i] = t{
1924			uint64(r.Next().Int64()),
1925			uint64(r.Next().Int64()),
1926		}
1927	}
1928	runtime.GC()
1929	b.StartTimer()
1930	for i := 0; i < b.N; i++ {
1931		v := a[i&(N-1)]
1932		MulUint128_64(v.a, v.b)
1933	}
1934}
1935
1936func BenchmarkMul128Big(b *testing.B) {
1937	const N = 1 << 16
1938	b.StopTimer()
1939	type t struct{ a, b *big.Int }
1940	a := make([]t, N)
1941	r := r64()
1942	for i := range a {
1943		a[i] = t{
1944			big.NewInt(r.Next().Int64()),
1945			big.NewInt(r.Next().Int64()),
1946		}
1947	}
1948	var x big.Int
1949	runtime.GC()
1950	b.StartTimer()
1951	for i := 0; i < b.N; i++ {
1952		v := a[i&(N-1)]
1953		x.Mul(v.a, v.b)
1954	}
1955}
1956
1957func TestIsPrimeUint64(t *testing.T) {
1958	f := func(lo, hi uint64, exp int) {
1959		got := 0
1960		for n := lo; n <= hi; {
1961			if IsPrimeUint64(n) {
1962				got++
1963			}
1964			n0 := n
1965			n++
1966			if n < n0 {
1967				break
1968			}
1969		}
1970		if got != exp {
1971			t.Fatal(lo, hi, got, exp)
1972		}
1973	}
1974
1975	// lo, hi, PrimePi(hi)-PrimePi(lo)
1976	f(0, 1e4, 1229)
1977	f(1e5, 1e5+1e4, 861)
1978	f(1e6, 1e6+1e4, 753)
1979	f(1e7, 1e7+1e4, 614)
1980	f(1e8, 1e8+1e4, 551)
1981	f(1e9, 1e9+1e4, 487)
1982	f(1e10, 1e10+1e4, 406)
1983	f(1e11, 1e11+1e4, 394)
1984	f(1e12, 1e12+1e4, 335)
1985	f(1e13, 1e13+1e4, 354)
1986	f(1e14, 1e14+1e4, 304)
1987	f(1e15, 1e15+1e4, 263)
1988	f(1e16, 1e16+1e4, 270)
1989	f(1e17, 1e17+1e4, 265)
1990	f(1e18, 1e18+1e4, 241)
1991	f(1e19, 1e19+1e4, 255)
1992	f(1<<64-1e4, 1<<64-1, 218)
1993}
1994
1995func TestProbablyPrimeUint32(t *testing.T) {
1996	f := func(n, firstFail uint32, primes []uint32) {
1997		for ; n <= firstFail; n += 2 {
1998			prp := true
1999			for _, a := range primes {
2000				if !ProbablyPrimeUint32(n, a) {
2001					prp = false
2002					break
2003				}
2004			}
2005			if prp != IsPrime(n) && n != firstFail {
2006				t.Fatal(n)
2007			}
2008		}
2009	}
2010	if !ProbablyPrimeUint32(5, 2) {
2011		t.Fatal(false)
2012	}
2013	if !ProbablyPrimeUint32(7, 2) {
2014		t.Fatal(false)
2015	}
2016	if ProbablyPrimeUint32(9, 2) {
2017		t.Fatal(true)
2018	}
2019	if !ProbablyPrimeUint32(11, 2) {
2020		t.Fatal(false)
2021	}
2022	// http://oeis.org/A014233
2023	f(5, 2047, []uint32{2})
2024	f(2047, 1373653, []uint32{2, 3})
2025	f(1373653, 25326001, []uint32{2, 3, 5})
2026}
2027
2028func BenchmarkProbablyPrimeUint32(b *testing.B) {
2029	const N = 1 << 16
2030	b.StopTimer()
2031	type t struct{ n, a uint32 }
2032	data := make([]t, N)
2033	r := r32()
2034	for i := range data {
2035		n := uint32(r.Next()) | 1
2036		if n <= 3 {
2037			n += 5
2038		}
2039		a := uint32(r.Next())
2040		if a <= 1 {
2041			a += 2
2042		}
2043		data[i] = t{n, a}
2044	}
2045	runtime.GC()
2046	b.StartTimer()
2047	for i := 0; i < b.N; i++ {
2048		v := data[i&(N-1)]
2049		ProbablyPrimeUint32(v.n, v.a)
2050	}
2051}
2052
2053func TestProbablyPrimeUint64_32(t *testing.T) {
2054	f := func(n, firstFail uint64, primes []uint32) {
2055		for ; n <= firstFail; n += 2 {
2056			prp := true
2057			for _, a := range primes {
2058				if !ProbablyPrimeUint64_32(n, a) {
2059					prp = false
2060					break
2061				}
2062			}
2063			if prp != IsPrimeUint64(n) && n != firstFail {
2064				t.Fatal(n)
2065			}
2066		}
2067	}
2068	if !ProbablyPrimeUint64_32(5, 2) {
2069		t.Fatal(false)
2070	}
2071	if !ProbablyPrimeUint64_32(7, 2) {
2072		t.Fatal(false)
2073	}
2074	if ProbablyPrimeUint64_32(9, 2) {
2075		t.Fatal(true)
2076	}
2077	if !ProbablyPrimeUint64_32(11, 2) {
2078		t.Fatal(false)
2079	}
2080	// http://oeis.org/A014233
2081	f(5, 2047, []uint32{2})
2082	f(2047, 1373653, []uint32{2, 3})
2083}
2084
2085func BenchmarkProbablyPrimeUint64_32(b *testing.B) {
2086	const N = 1 << 16
2087	b.StopTimer()
2088	type t struct {
2089		n uint64
2090		a uint32
2091	}
2092	data := make([]t, N)
2093	r := r32()
2094	r2 := r64()
2095	for i := range data {
2096		var n uint64
2097		for n <= 3 {
2098			n = uint64(r2.Next().Int64()) | 1
2099		}
2100		var a uint32
2101		for a <= 1 {
2102			a = uint32(r.Next())
2103		}
2104		data[i] = t{n, a}
2105	}
2106	runtime.GC()
2107	b.StartTimer()
2108	for i := 0; i < b.N; i++ {
2109		v := data[i&(N-1)]
2110		ProbablyPrimeUint64_32(v.n, v.a)
2111	}
2112}
2113
2114func TestProbablyPrimeBigInt_32(t *testing.T) {
2115	f := func(n0, firstFail0 uint64, primes []uint32) {
2116		n, firstFail := big.NewInt(0).SetUint64(n0), big.NewInt(0).SetUint64(firstFail0)
2117		for ; n.Cmp(firstFail) <= 0; n.Add(n, _2) {
2118			prp := true
2119			for _, a := range primes {
2120				if !ProbablyPrimeBigInt_32(n, a) {
2121					prp = false
2122					break
2123				}
2124			}
2125			if prp != IsPrimeUint64(n0) && n0 != firstFail0 {
2126				t.Fatal(n)
2127			}
2128			n0 += 2
2129		}
2130	}
2131	if !ProbablyPrimeBigInt_32(big.NewInt(5), 2) {
2132		t.Fatal(false)
2133	}
2134	if !ProbablyPrimeBigInt_32(big.NewInt(7), 2) {
2135		t.Fatal(false)
2136	}
2137	if ProbablyPrimeBigInt_32(big.NewInt(9), 2) {
2138		t.Fatal(true)
2139	}
2140	if !ProbablyPrimeBigInt_32(big.NewInt(11), 2) {
2141		t.Fatal(false)
2142	}
2143	// http://oeis.org/A014233
2144	f(5, 2047, []uint32{2})
2145	f(2047, 1373653, []uint32{2, 3})
2146}
2147
2148func BenchmarkProbablyPrimeBigInt_32(b *testing.B) {
2149	const N = 1 << 16
2150	b.StopTimer()
2151	type t struct {
2152		n *big.Int
2153		a uint32
2154	}
2155	data := make([]t, N)
2156	r := r32()
2157	r2 := r64()
2158	for i := range data {
2159		var n uint64
2160		for n <= 3 {
2161			n = uint64(r2.Next().Int64()) | 1
2162		}
2163		var a uint32
2164		for a <= 1 {
2165			a = uint32(r.Next())
2166		}
2167		data[i] = t{big.NewInt(0).SetUint64(n), a}
2168	}
2169	runtime.GC()
2170	b.StartTimer()
2171	for i := 0; i < b.N; i++ {
2172		v := data[i&(N-1)]
2173		ProbablyPrimeBigInt_32(v.n, v.a)
2174	}
2175}
2176
2177func TestProbablyPrimeBigInt(t *testing.T) {
2178	f := func(n0, firstFail0 uint64, primes []uint32) {
2179		n, firstFail := big.NewInt(0).SetUint64(n0), big.NewInt(0).SetUint64(firstFail0)
2180		for ; n.Cmp(firstFail) <= 0; n.Add(n, _2) {
2181			prp := true
2182			var a big.Int
2183			for _, a0 := range primes {
2184				a.SetInt64(int64(a0))
2185				if !ProbablyPrimeBigInt(n, &a) {
2186					prp = false
2187					break
2188				}
2189			}
2190			if prp != IsPrimeUint64(n0) && n0 != firstFail0 {
2191				t.Fatal(n)
2192			}
2193			n0 += 2
2194		}
2195	}
2196	if !ProbablyPrimeBigInt(big.NewInt(5), _2) {
2197		t.Fatal(false)
2198	}
2199	if !ProbablyPrimeBigInt(big.NewInt(7), _2) {
2200		t.Fatal(false)
2201	}
2202	if ProbablyPrimeBigInt(big.NewInt(9), _2) {
2203		t.Fatal(true)
2204	}
2205	if !ProbablyPrimeBigInt(big.NewInt(11), _2) {
2206		t.Fatal(false)
2207	}
2208	// http://oeis.org/A014233
2209	f(5, 2047, []uint32{2})
2210	f(2047, 1373653, []uint32{2, 3})
2211}
2212
2213var once2059 sync.Once
2214
2215func BenchmarkProbablyPrimeBigInt64(b *testing.B) {
2216	const N = 1 << 16
2217	b.StopTimer()
2218	once2059.Do(func() { b.Log("64 bit n, 64 bit a\n") })
2219	type t struct {
2220		n, a *big.Int
2221	}
2222	data := make([]t, N)
2223	r := r64()
2224	for i := range data {
2225		var n uint64
2226		for n <= 3 {
2227			n = uint64(r.Next().Int64()) | 1
2228		}
2229		var a uint64
2230		for a <= 1 {
2231			a = uint64(r.Next().Int64())
2232		}
2233		data[i] = t{big.NewInt(0).SetUint64(n), big.NewInt(0).SetUint64(uint64(a))}
2234	}
2235	runtime.GC()
2236	b.StartTimer()
2237	for i := 0; i < b.N; i++ {
2238		v := data[i&(N-1)]
2239		ProbablyPrimeBigInt(v.n, v.a)
2240	}
2241}
2242
2243var once2090 sync.Once
2244
2245func BenchmarkProbablyPrimeBigInt128(b *testing.B) {
2246	const N = 1 << 16
2247	b.StopTimer()
2248	once2090.Do(func() { b.Log("128 bit n, 128 bit a\n") })
2249	type t struct {
2250		n, a *big.Int
2251	}
2252	data := make([]t, N)
2253	r := r64()
2254	for i := range data {
2255		n := big.NewInt(0).SetUint64(uint64(r.Next().Int64()))
2256		n.Lsh(n, 64)
2257		n.Add(n, big.NewInt(0).SetUint64(uint64(r.Next().Int64())|1))
2258		a := big.NewInt(0).SetUint64(uint64(r.Next().Int64()))
2259		a.Lsh(a, 64)
2260		a.Add(a, big.NewInt(0).SetUint64(uint64(r.Next().Int64())))
2261		data[i] = t{n, a}
2262	}
2263	runtime.GC()
2264	b.StartTimer()
2265	for i := 0; i < b.N; i++ {
2266		v := data[i&(N-1)]
2267		ProbablyPrimeBigInt(v.n, v.a)
2268	}
2269}
2270
2271func TestQCmpUint32(t *testing.T) {
2272	const N = 6e4
2273	r := r32()
2274	var x, y big.Rat
2275	for i := 0; i < N; i++ {
2276		a, b, c, d := uint32(r.Next()), uint32(r.Next()), uint32(r.Next()), uint32(r.Next())
2277		x.SetFrac64(int64(a), int64(b))
2278		y.SetFrac64(int64(c), int64(d))
2279		if g, e := QCmpUint32(a, b, c, d), x.Cmp(&y); g != e {
2280			t.Fatal(a, b, c, d, g, e)
2281		}
2282	}
2283}
2284
2285func TestQScaleUint32(t *testing.T) {
2286	const N = 4e4
2287	r := r32()
2288	var x, y big.Rat
2289	var a uint64
2290	var b, c, d uint32
2291	for i := 0; i < N; i++ {
2292		for {
2293			b, c, d = uint32(r.Next()), uint32(r.Next()), uint32(r.Next())
2294			a = QScaleUint32(b, c, d)
2295			if a <= math.MaxInt64 {
2296				break
2297			}
2298		}
2299		x.SetFrac64(int64(a), int64(b))
2300		y.SetFrac64(int64(c), int64(d))
2301		if g := x.Cmp(&y); g < 0 {
2302			t.Fatal(a, b, c, d, g, "expexted 1 or 0")
2303		}
2304
2305		if a != 0 {
2306			x.SetFrac64(int64(a-1), int64(b))
2307			if g := x.Cmp(&y); g > 0 {
2308				t.Fatal(a, b, c, d, g, "expected -1 or 0")
2309			}
2310		}
2311	}
2312}
2313
2314var smalls = []uint32{2, 3, 5, 7, 11, 13, 17, 19, 23, 29}
2315
2316func isPrimorialProduct(t FactorTerms, maxp uint32) bool {
2317	if len(t) == 0 {
2318		return false
2319	}
2320
2321	pmax := uint32(32)
2322	for i, v := range t {
2323		if v.Prime != smalls[i] || v.Power > pmax || v.Power > maxp {
2324			return false
2325		}
2326		pmax = v.Power
2327	}
2328	return true
2329}
2330
2331func TestPrimorialProductsUint32(t *testing.T) {
2332	r := PrimorialProductsUint32(2*3*5*7*11*13*17*19+1, math.MaxUint32, 1)
2333	if len(r) != 1 {
2334		t.Fatal(len(r))
2335	}
2336
2337	if r[0] != 2*3*5*7*11*13*17*19*23 {
2338		t.Fatal(r[0])
2339	}
2340
2341	r = PrimorialProductsUint32(0, math.MaxUint32, math.MaxUint32)
2342	if g, e := len(r), 1679; g != e {
2343		t.Fatal(g, e)
2344	}
2345
2346	m := map[uint32]struct{}{}
2347	for _, v := range r {
2348		if _, ok := m[v]; ok {
2349			t.Fatal(v)
2350		}
2351
2352		m[v] = struct{}{}
2353	}
2354
2355	for lo := uint32(0); lo < 5e4; lo += 1e3 {
2356		hi := 1e2 * lo
2357		for max := uint32(0); max <= 32; max++ {
2358			m := map[uint32]struct{}{}
2359			for i, v := range PrimorialProductsUint32(lo, hi, max) {
2360				f := FactorInt(v)
2361				if v < lo || v > hi {
2362					t.Fatal(lo, hi, max, v)
2363				}
2364
2365				if _, ok := m[v]; ok {
2366					t.Fatal(i, lo, hi, max, v, f)
2367				}
2368
2369				m[v] = struct{}{}
2370				if !isPrimorialProduct(f, max) {
2371					t.Fatal(i, v)
2372				}
2373
2374				for _, v := range f {
2375					if v.Power > max {
2376						t.Fatal(i, v, f)
2377					}
2378				}
2379			}
2380		}
2381	}
2382}
2383
2384func BenchmarkPrimorialProductsUint32(b *testing.B) {
2385	for i := 0; i < b.N; i++ {
2386		PrimorialProductsUint32(0, math.MaxUint32, math.MaxUint32)
2387	}
2388}
2389
2390func powerizeUint32BigInt(b uint32, n *big.Int) (e uint32, p *big.Int) {
2391	p = big.NewInt(1)
2392	bb := big.NewInt(int64(b))
2393	for p.Cmp(n) < 0 {
2394		p.Mul(p, bb)
2395		e++
2396	}
2397	return
2398}
2399
2400func TestPowerizeUint32BigInt(t *testing.T) {
2401	var data = []struct{ b, n, e, p int }{
2402		{0, 10, 0, -1},
2403		{1, 10, 0, -1},
2404		{2, -1, 0, -1},
2405		{2, 0, 0, 1},
2406		{2, 1, 0, 1},
2407		{2, 2, 1, 2},
2408		{2, 3, 2, 4},
2409		{3, 0, 0, 1},
2410		{3, 1, 0, 1},
2411		{3, 2, 1, 3},
2412		{3, 3, 1, 3},
2413		{3, 4, 2, 9},
2414		{3, 8, 2, 9},
2415		{3, 9, 2, 9},
2416		{3, 10, 3, 27},
2417		{3, 80, 4, 81},
2418	}
2419
2420	var n big.Int
2421	for _, v := range data {
2422		b := v.b
2423		n.SetInt64(int64(v.n))
2424		e, p := PowerizeUint32BigInt(uint32(b), &n)
2425		if e != uint32(v.e) {
2426			t.Fatal(b, &n, e, p, v.e, v.p)
2427		}
2428
2429		if v.p < 0 {
2430			if p != nil {
2431				t.Fatal(b, &n, e, p, v.e, v.p)
2432			}
2433			continue
2434		}
2435
2436		if p.Int64() != int64(v.p) {
2437			t.Fatal(b, &n, e, p, v.e, v.p)
2438		}
2439	}
2440	const N = 1e5
2441	var nn big.Int
2442	for _, base := range []uint32{2, 3, 15, 17} {
2443		for n := 0; n <= N; n++ {
2444			nn.SetInt64(int64(n))
2445			ge, gp := PowerizeUint32BigInt(base, &nn)
2446			ee, ep := powerizeUint32BigInt(base, &nn)
2447			if ge != ee || gp.Cmp(ep) != 0 {
2448				t.Fatal(base, n, ge, gp, ee, ep)
2449			}
2450
2451			gp.Div(gp, big.NewInt(int64(base)))
2452			if gp.Sign() > 0 && gp.Cmp(&nn) >= 0 {
2453				t.Fatal(gp.Sign(), gp.Cmp(&nn))
2454			}
2455		}
2456	}
2457}
2458
2459func benchmarkPowerizeUint32BigInt(b *testing.B, base uint32, exp int) {
2460	b.StopTimer()
2461	var n big.Int
2462	n.SetBit(&n, exp, 1)
2463	b.StartTimer()
2464	for i := 0; i < b.N; i++ {
2465		PowerizeUint32BigInt(base, &n)
2466	}
2467}
2468
2469func BenchmarkPowerizeUint32BigInt_2_2e1e1(b *testing.B) {
2470	benchmarkPowerizeUint32BigInt(b, 2, 1e1)
2471}
2472
2473func BenchmarkPowerizeUint32BigInt_2_2e1e2(b *testing.B) {
2474	benchmarkPowerizeUint32BigInt(b, 2, 1e2)
2475}
2476
2477func BenchmarkPowerizeUint32BigInt_2_2e1e3(b *testing.B) {
2478	benchmarkPowerizeUint32BigInt(b, 2, 1e3)
2479}
2480
2481func BenchmarkPowerizeUint32BigInt_2_2e1e4(b *testing.B) {
2482	benchmarkPowerizeUint32BigInt(b, 2, 1e4)
2483}
2484
2485func BenchmarkPowerizeUint32BigInt_2_2e1e5(b *testing.B) {
2486	benchmarkPowerizeUint32BigInt(b, 2, 1e5)
2487}
2488
2489func BenchmarkPowerizeUint32BigInt_2_2e1e6(b *testing.B) {
2490	benchmarkPowerizeUint32BigInt(b, 2, 1e6)
2491}
2492
2493func BenchmarkPowerizeUint32BigInt_2_2e1e7(b *testing.B) {
2494	benchmarkPowerizeUint32BigInt(b, 2, 1e7)
2495}
2496
2497func BenchmarkPowerizeUint32BigInt_3_2e1e1(b *testing.B) {
2498	benchmarkPowerizeUint32BigInt(b, 3, 1e1)
2499}
2500
2501func BenchmarkPowerizeUint32BigInt_3_2e1e2(b *testing.B) {
2502	benchmarkPowerizeUint32BigInt(b, 3, 1e2)
2503}
2504
2505func BenchmarkPowerizeUint32BigInt_3_2e1e3(b *testing.B) {
2506	benchmarkPowerizeUint32BigInt(b, 3, 1e3)
2507}
2508
2509func BenchmarkPowerizeUint32BigInt_3_2e1e4(b *testing.B) {
2510	benchmarkPowerizeUint32BigInt(b, 3, 1e4)
2511}
2512
2513func BenchmarkPowerizeUint32BigInt_3_2e1e5(b *testing.B) {
2514	benchmarkPowerizeUint32BigInt(b, 3, 1e5)
2515}
2516
2517func BenchmarkPowerizeUint32BigInt_3_2e1e6(b *testing.B) {
2518	benchmarkPowerizeUint32BigInt(b, 3, 1e6)
2519}
2520
2521func BenchmarkPowerizeUint32BigInt_15_2e1e1(b *testing.B) {
2522	benchmarkPowerizeUint32BigInt(b, 15, 1e1)
2523}
2524
2525func BenchmarkPowerizeUint32BigInt_15_2e1e2(b *testing.B) {
2526	benchmarkPowerizeUint32BigInt(b, 15, 1e2)
2527}
2528
2529func BenchmarkPowerizeUint32BigInt_15_2e1e3(b *testing.B) {
2530	benchmarkPowerizeUint32BigInt(b, 15, 1e3)
2531}
2532
2533func BenchmarkPowerizeUint32BigInt_15_2e1e4(b *testing.B) {
2534	benchmarkPowerizeUint32BigInt(b, 15, 1e4)
2535}
2536
2537func BenchmarkPowerizeUint32BigInt_15_2e1e5(b *testing.B) {
2538	benchmarkPowerizeUint32BigInt(b, 15, 1e5)
2539}
2540
2541func BenchmarkPowerizeUint32BigInt_15_2e1e6(b *testing.B) {
2542	benchmarkPowerizeUint32BigInt(b, 15, 1e6)
2543}
2544
2545func BenchmarkPowerizeUint32BigInt_17_2e1e1(b *testing.B) {
2546	benchmarkPowerizeUint32BigInt(b, 17, 1e1)
2547}
2548
2549func BenchmarkPowerizeUint32BigInt_17_2e1e2(b *testing.B) {
2550	benchmarkPowerizeUint32BigInt(b, 17, 1e2)
2551}
2552
2553func BenchmarkPowerizeUint32BigInt_17_2e1e3(b *testing.B) {
2554	benchmarkPowerizeUint32BigInt(b, 17, 1e3)
2555}
2556
2557func BenchmarkPowerizeUint32BigInt_17_2e1e4(b *testing.B) {
2558	benchmarkPowerizeUint32BigInt(b, 17, 1e4)
2559}
2560
2561func BenchmarkPowerizeUint32BigInt_17_2e1e5(b *testing.B) {
2562	benchmarkPowerizeUint32BigInt(b, 17, 1e5)
2563}
2564
2565func BenchmarkPowerizeUint32BigInt_17_2e1e6(b *testing.B) {
2566	benchmarkPowerizeUint32BigInt(b, 17, 1e6)
2567}
2568
2569func TestPowerizeBigInt(t *testing.T) {
2570	var data = []struct{ b, n, e, p int }{
2571		{0, 10, 0, -1},
2572		{1, 10, 0, -1},
2573		{2, -1, 0, -1},
2574		{2, 0, 0, 1},
2575		{2, 1, 0, 1},
2576		{2, 2, 1, 2},
2577		{2, 3, 2, 4},
2578		{3, 0, 0, 1},
2579		{3, 1, 0, 1},
2580		{3, 2, 1, 3},
2581		{3, 3, 1, 3},
2582		{3, 4, 2, 9},
2583		{3, 8, 2, 9},
2584		{3, 9, 2, 9},
2585		{3, 10, 3, 27},
2586		{3, 80, 4, 81},
2587	}
2588
2589	var b, n big.Int
2590	for _, v := range data {
2591		b.SetInt64(int64(v.b))
2592		n.SetInt64(int64(v.n))
2593		e, p := PowerizeBigInt(&b, &n)
2594		if e != uint32(v.e) {
2595			t.Fatal(&b, &n, e, p, v.e, v.p)
2596		}
2597
2598		if v.p < 0 {
2599			if p != nil {
2600				t.Fatal(&b, &n, e, p, v.e, v.p)
2601			}
2602			continue
2603		}
2604
2605		if p.Int64() != int64(v.p) {
2606			t.Fatal(&b, &n, e, p, v.e, v.p)
2607		}
2608	}
2609	const N = 1e5
2610	var nn big.Int
2611	for _, base := range []uint32{2, 3, 15, 17} {
2612		b.SetInt64(int64(base))
2613		for n := 0; n <= N; n++ {
2614			nn.SetInt64(int64(n))
2615			ge, gp := PowerizeBigInt(&b, &nn)
2616			ee, ep := powerizeUint32BigInt(base, &nn)
2617			if ge != ee || gp.Cmp(ep) != 0 {
2618				t.Fatal(base, n, ge, gp, ee, ep)
2619			}
2620
2621			gp.Div(gp, &b)
2622			if gp.Sign() > 0 && gp.Cmp(&nn) >= 0 {
2623				t.Fatal(gp.Sign(), gp.Cmp(&nn))
2624			}
2625		}
2626	}
2627}
2628
2629func benchmarkPowerizeBigInt(b *testing.B, base uint32, exp int) {
2630	b.StopTimer()
2631	var bb, n big.Int
2632	n.SetBit(&n, exp, 1)
2633	bb.SetInt64(int64(base))
2634	b.StartTimer()
2635	for i := 0; i < b.N; i++ {
2636		PowerizeBigInt(&bb, &n)
2637	}
2638}
2639
2640func BenchmarkPowerizeBigInt_2_2e1e1(b *testing.B) {
2641	benchmarkPowerizeBigInt(b, 2, 1e1)
2642}
2643
2644func BenchmarkPowerizeBigInt_2_2e1e2(b *testing.B) {
2645	benchmarkPowerizeBigInt(b, 2, 1e2)
2646}
2647
2648func BenchmarkPowerizeBigInt_2_2e1e3(b *testing.B) {
2649	benchmarkPowerizeBigInt(b, 2, 1e3)
2650}
2651
2652func BenchmarkPowerizeBigInt_2_2e1e4(b *testing.B) {
2653	benchmarkPowerizeBigInt(b, 2, 1e4)
2654}
2655
2656func BenchmarkPowerizeBigInt_2_2e1e5(b *testing.B) {
2657	benchmarkPowerizeBigInt(b, 2, 1e5)
2658}
2659
2660func BenchmarkPowerizeBigInt_2_2e1e6(b *testing.B) {
2661	benchmarkPowerizeBigInt(b, 2, 1e6)
2662}
2663
2664func BenchmarkPowerizeBigInt_2_2e1e7(b *testing.B) {
2665	benchmarkPowerizeBigInt(b, 2, 1e7)
2666}
2667
2668func BenchmarkPowerizeBigInt_3_2e1e1(b *testing.B) {
2669	benchmarkPowerizeBigInt(b, 3, 1e1)
2670}
2671
2672func BenchmarkPowerizeBigInt_3_2e1e2(b *testing.B) {
2673	benchmarkPowerizeBigInt(b, 3, 1e2)
2674}
2675
2676func BenchmarkPowerizeBigInt_3_2e1e3(b *testing.B) {
2677	benchmarkPowerizeBigInt(b, 3, 1e3)
2678}
2679
2680func BenchmarkPowerizeBigInt_3_2e1e4(b *testing.B) {
2681	benchmarkPowerizeBigInt(b, 3, 1e4)
2682}
2683
2684func BenchmarkPowerizeBigInt_3_2e1e5(b *testing.B) {
2685	benchmarkPowerizeBigInt(b, 3, 1e5)
2686}
2687
2688func BenchmarkPowerizeBigInt_3_2e1e6(b *testing.B) {
2689	benchmarkPowerizeBigInt(b, 3, 1e6)
2690}
2691
2692func BenchmarkPowerizeBigInt_15_2e1e1(b *testing.B) {
2693	benchmarkPowerizeBigInt(b, 15, 1e1)
2694}
2695
2696func BenchmarkPowerizeBigInt_15_2e1e2(b *testing.B) {
2697	benchmarkPowerizeBigInt(b, 15, 1e2)
2698}
2699
2700func BenchmarkPowerizeBigInt_15_2e1e3(b *testing.B) {
2701	benchmarkPowerizeBigInt(b, 15, 1e3)
2702}
2703
2704func BenchmarkPowerizeBigInt_15_2e1e4(b *testing.B) {
2705	benchmarkPowerizeBigInt(b, 15, 1e4)
2706}
2707
2708func BenchmarkPowerizeBigInt_15_2e1e5(b *testing.B) {
2709	benchmarkPowerizeBigInt(b, 15, 1e5)
2710}
2711
2712func BenchmarkPowerizeBigInt_15_2e1e6(b *testing.B) {
2713	benchmarkPowerizeBigInt(b, 15, 1e6)
2714}
2715
2716func BenchmarkPowerizeBigInt_17_2e1e1(b *testing.B) {
2717	benchmarkPowerizeBigInt(b, 17, 1e1)
2718}
2719
2720func BenchmarkPowerizeBigInt_17_2e1e2(b *testing.B) {
2721	benchmarkPowerizeBigInt(b, 17, 1e2)
2722}
2723
2724func BenchmarkPowerizeBigInt_17_2e1e3(b *testing.B) {
2725	benchmarkPowerizeBigInt(b, 17, 1e3)
2726}
2727
2728func BenchmarkPowerizeBigInt_17_2e1e4(b *testing.B) {
2729	benchmarkPowerizeBigInt(b, 17, 1e4)
2730}
2731
2732func BenchmarkPowerizeBigInt_17_2e1e5(b *testing.B) {
2733	benchmarkPowerizeBigInt(b, 17, 1e5)
2734}
2735
2736func BenchmarkPowerizeBigInt_17_2e1e6(b *testing.B) {
2737	benchmarkPowerizeBigInt(b, 17, 1e6)
2738}
2739
2740func TestEnvelope(t *testing.T) {
2741	const prec = 1e-3
2742	type check struct {
2743		approx Approximation
2744		x, y   float64
2745	}
2746	data := []struct {
2747		points []float64
2748		checks []check
2749	}{
2750		{
2751			[]float64{0, 1},
2752			[]check{
2753				{Linear, 0, 0},
2754				{Linear, 0.25, 0.25},
2755				{Linear, 0.5, 0.5},
2756				{Linear, 0.75, 0.75},
2757				{Linear, 0.9999, 1},
2758			},
2759		},
2760		{
2761			[]float64{-1, 0},
2762			[]check{
2763				{Linear, 0, -1},
2764				{Linear, 0.25, -0.75},
2765				{Linear, 0.5, -0.5},
2766				{Linear, 0.75, -0.25},
2767				{Linear, 0.9999, 0},
2768			},
2769		},
2770		{
2771			[]float64{-1, 1},
2772			[]check{
2773				{Linear, 0, -1},
2774				{Linear, 0.25, -0.5},
2775				{Linear, 0.5, 0},
2776				{Linear, 0.75, 0.5},
2777				{Linear, 0.9999, 1},
2778			},
2779		},
2780		{
2781			[]float64{-1, 1, -2},
2782			[]check{
2783				{Linear, 0, -1},
2784				{Linear, 0.25, 0},
2785				{Linear, 0.5, 1},
2786				{Linear, 0.75, -0.5},
2787				{Linear, 0.9, -1.4},
2788				{Linear, 0.9999, -2},
2789			},
2790		},
2791		{
2792			[]float64{-1, 1},
2793			[]check{
2794				{Sinusoidal, 0, -1},
2795				{Sinusoidal, 0.25, -math.Sqrt2 / 2},
2796				{Sinusoidal, 0.5, 0},
2797				{Sinusoidal, 0.75, math.Sqrt2 / 2},
2798				{Sinusoidal, 0.9999, 1},
2799			},
2800		},
2801		{
2802			[]float64{-1, 1, -2},
2803			[]check{
2804				{Sinusoidal, 0, -1},
2805				{Sinusoidal, 1. / 8, -math.Sqrt2 / 2},
2806				{Sinusoidal, 2. / 8, 0},
2807				{Sinusoidal, 3. / 8, math.Sqrt2 / 2},
2808				{Sinusoidal, 4. / 8, 1},
2809				{Sinusoidal, 5. / 8, (3*math.Sqrt2 - 2) / 4},
2810				{Sinusoidal, 6. / 8, -0.5},
2811				{Sinusoidal, 7. / 8, (-3*math.Sqrt2 - 2) / 4},
2812				{Sinusoidal, 0.9999, -2},
2813			},
2814		},
2815	}
2816	for i, suite := range data {
2817		for j, test := range suite.checks {
2818			e, g := test.y, Envelope(test.x, suite.points, test.approx)
2819			d := math.Abs(e - g)
2820			if d > prec {
2821				t.Errorf(
2822					"i %d, j %d, x %v, e %v, g %v, d %v, prec %v",
2823					i, j, test.x, e, g, d, prec,
2824				)
2825			}
2826		}
2827	}
2828}
2829
2830func TestMaxInt(t *testing.T) {
2831	n := int64(MaxInt)
2832	if n != math.MaxInt32 && n != math.MaxInt64 {
2833		t.Error(n)
2834	}
2835
2836	t.Logf("64 bit ints: %t, MaxInt: %d", n == math.MaxInt64, n)
2837}
2838
2839func TestMinInt(t *testing.T) {
2840	n := int64(MinInt)
2841	if n != math.MinInt32 && n != math.MinInt64 {
2842		t.Error(n)
2843	}
2844
2845	t.Logf("64 bit ints: %t. MinInt: %d", n == math.MinInt64, n)
2846}
2847
2848func TestMaxUint(t *testing.T) {
2849	n := uint64(MaxUint)
2850	if n != math.MaxUint32 && n != math.MaxUint64 {
2851		t.Error(n)
2852	}
2853
2854	t.Logf("64 bit uints: %t. MaxUint: %d", n == math.MaxUint64, n)
2855}
2856
2857func TestMax(t *testing.T) {
2858	tests := []struct{ a, b, e int }{
2859		{MinInt, MinIntM1, MaxInt},
2860		{MinIntM1, MinInt, MaxInt},
2861		{MinIntM1, MinIntM1, MaxInt},
2862
2863		{MinInt, MinInt, MinInt},
2864		{MinInt + 1, MinInt, MinInt + 1},
2865		{MinInt, MinInt + 1, MinInt + 1},
2866
2867		{-1, -1, -1},
2868		{-1, 0, 0},
2869		{-1, 1, 1},
2870
2871		{0, -1, 0},
2872		{0, 0, 0},
2873		{0, 1, 1},
2874
2875		{1, -1, 1},
2876		{1, 0, 1},
2877		{1, 1, 1},
2878
2879		{MaxInt, MaxInt, MaxInt},
2880		{MaxInt - 1, MaxInt, MaxInt},
2881		{MaxInt, MaxInt - 1, MaxInt},
2882
2883		{MaxIntP1, MaxInt, MaxInt},
2884		{MaxInt, MaxIntP1, MaxInt},
2885		{MaxIntP1, MaxIntP1, MinInt},
2886	}
2887
2888	for _, test := range tests {
2889		if g, e := Max(test.a, test.b), test.e; g != e {
2890			t.Fatal(test.a, test.b, g, e)
2891		}
2892	}
2893}
2894
2895func TestMin(t *testing.T) {
2896	tests := []struct{ a, b, e int }{
2897		{MinIntM1, MinInt, MinInt},
2898		{MinInt, MinIntM1, MinInt},
2899		{MinIntM1, MinIntM1, MaxInt},
2900
2901		{MinInt, MinInt, MinInt},
2902		{MinInt + 1, MinInt, MinInt},
2903		{MinInt, MinInt + 1, MinInt},
2904
2905		{-1, -1, -1},
2906		{-1, 0, -1},
2907		{-1, 1, -1},
2908
2909		{0, -1, -1},
2910		{0, 0, 0},
2911		{0, 1, 0},
2912
2913		{1, -1, -1},
2914		{1, 0, 0},
2915		{1, 1, 1},
2916
2917		{MaxInt, MaxInt, MaxInt},
2918		{MaxInt - 1, MaxInt, MaxInt - 1},
2919		{MaxInt, MaxInt - 1, MaxInt - 1},
2920
2921		{MaxIntP1, MaxInt, MinInt},
2922		{MaxInt, MaxIntP1, MinInt},
2923		{MaxIntP1, MaxIntP1, MinInt},
2924	}
2925
2926	for _, test := range tests {
2927		if g, e := Min(test.a, test.b), test.e; g != e {
2928			t.Fatal(test.a, test.b, g, e)
2929		}
2930	}
2931}
2932
2933func TestUMax(t *testing.T) {
2934	tests := []struct{ a, b, e uint }{
2935		{0, 0, 0},
2936		{0, 1, 1},
2937		{1, 0, 1},
2938
2939		{10, 10, 10},
2940		{10, 11, 11},
2941		{11, 10, 11},
2942		{11, 11, 11},
2943
2944		{MaxUint, MaxUint, MaxUint},
2945		{MaxUint, MaxUint - 1, MaxUint},
2946		{MaxUint - 1, MaxUint, MaxUint},
2947		{MaxUint - 1, MaxUint - 1, MaxUint - 1},
2948
2949		{MaxUint, MaxUintP1, MaxUint},
2950		{MaxUintP1, MaxUint, MaxUint},
2951		{MaxUintP1, MaxUintP1, 0},
2952	}
2953
2954	for _, test := range tests {
2955		if g, e := UMax(test.a, test.b), test.e; g != e {
2956			t.Fatal(test.a, test.b, g, e)
2957		}
2958	}
2959}
2960
2961func TestUMin(t *testing.T) {
2962	tests := []struct{ a, b, e uint }{
2963		{0, 0, 0},
2964		{0, 1, 0},
2965		{1, 0, 0},
2966
2967		{10, 10, 10},
2968		{10, 11, 10},
2969		{11, 10, 10},
2970		{11, 11, 11},
2971
2972		{MaxUint, MaxUint, MaxUint},
2973		{MaxUint, MaxUint - 1, MaxUint - 1},
2974		{MaxUint - 1, MaxUint, MaxUint - 1},
2975		{MaxUint - 1, MaxUint - 1, MaxUint - 1},
2976
2977		{MaxUint, MaxUintP1, 0},
2978		{MaxUintP1, MaxUint, 0},
2979		{MaxUintP1, MaxUintP1, 0},
2980	}
2981
2982	for _, test := range tests {
2983		if g, e := UMin(test.a, test.b), test.e; g != e {
2984			t.Fatal(test.a, test.b, g, e)
2985		}
2986	}
2987}
2988
2989func TestMaxByte(t *testing.T) {
2990	tests := []struct{ a, b, e byte }{
2991		{0, 0, 0},
2992		{0, 1, 1},
2993		{1, 0, 1},
2994
2995		{10, 10, 10},
2996		{10, 11, 11},
2997		{11, 10, 11},
2998		{11, 11, 11},
2999
3000		{math.MaxUint8, math.MaxUint8, math.MaxUint8},
3001		{math.MaxUint8, math.MaxUint8 - 1, math.MaxUint8},
3002		{math.MaxUint8 - 1, math.MaxUint8, math.MaxUint8},
3003		{math.MaxUint8 - 1, math.MaxUint8 - 1, math.MaxUint8 - 1},
3004	}
3005
3006	for _, test := range tests {
3007		if g, e := MaxByte(test.a, test.b), test.e; g != e {
3008			t.Fatal(test.a, test.b, g, e)
3009		}
3010	}
3011}
3012
3013func TestMinByte(t *testing.T) {
3014	tests := []struct{ a, b, e byte }{
3015		{0, 0, 0},
3016		{0, 1, 0},
3017		{1, 0, 0},
3018
3019		{10, 10, 10},
3020		{10, 11, 10},
3021		{11, 10, 10},
3022		{11, 11, 11},
3023
3024		{math.MaxUint8, math.MaxUint8, math.MaxUint8},
3025		{math.MaxUint8, math.MaxUint8 - 1, math.MaxUint8 - 1},
3026		{math.MaxUint8 - 1, math.MaxUint8, math.MaxUint8 - 1},
3027		{math.MaxUint8 - 1, math.MaxUint8 - 1, math.MaxUint8 - 1},
3028	}
3029
3030	for _, test := range tests {
3031		if g, e := MinByte(test.a, test.b), test.e; g != e {
3032			t.Fatal(test.a, test.b, g, e)
3033		}
3034	}
3035}
3036
3037func TestMaxUint16(t *testing.T) {
3038	tests := []struct{ a, b, e uint16 }{
3039		{0, 0, 0},
3040		{0, 1, 1},
3041		{1, 0, 1},
3042
3043		{10, 10, 10},
3044		{10, 11, 11},
3045		{11, 10, 11},
3046		{11, 11, 11},
3047
3048		{math.MaxUint16, math.MaxUint16, math.MaxUint16},
3049		{math.MaxUint16, math.MaxUint16 - 1, math.MaxUint16},
3050		{math.MaxUint16 - 1, math.MaxUint16, math.MaxUint16},
3051		{math.MaxUint16 - 1, math.MaxUint16 - 1, math.MaxUint16 - 1},
3052	}
3053
3054	for _, test := range tests {
3055		if g, e := MaxUint16(test.a, test.b), test.e; g != e {
3056			t.Fatal(test.a, test.b, g, e)
3057		}
3058	}
3059}
3060
3061func TestMinUint16(t *testing.T) {
3062	tests := []struct{ a, b, e uint16 }{
3063		{0, 0, 0},
3064		{0, 1, 0},
3065		{1, 0, 0},
3066
3067		{10, 10, 10},
3068		{10, 11, 10},
3069		{11, 10, 10},
3070		{11, 11, 11},
3071
3072		{math.MaxUint16, math.MaxUint16, math.MaxUint16},
3073		{math.MaxUint16, math.MaxUint16 - 1, math.MaxUint16 - 1},
3074		{math.MaxUint16 - 1, math.MaxUint16, math.MaxUint16 - 1},
3075		{math.MaxUint16 - 1, math.MaxUint16 - 1, math.MaxUint16 - 1},
3076	}
3077
3078	for _, test := range tests {
3079		if g, e := MinUint16(test.a, test.b), test.e; g != e {
3080			t.Fatal(test.a, test.b, g, e)
3081		}
3082	}
3083}
3084
3085func TestMaxUint32(t *testing.T) {
3086	tests := []struct{ a, b, e uint32 }{
3087		{0, 0, 0},
3088		{0, 1, 1},
3089		{1, 0, 1},
3090
3091		{10, 10, 10},
3092		{10, 11, 11},
3093		{11, 10, 11},
3094		{11, 11, 11},
3095
3096		{math.MaxUint32, math.MaxUint32, math.MaxUint32},
3097		{math.MaxUint32, math.MaxUint32 - 1, math.MaxUint32},
3098		{math.MaxUint32 - 1, math.MaxUint32, math.MaxUint32},
3099		{math.MaxUint32 - 1, math.MaxUint32 - 1, math.MaxUint32 - 1},
3100	}
3101
3102	for _, test := range tests {
3103		if g, e := MaxUint32(test.a, test.b), test.e; g != e {
3104			t.Fatal(test.a, test.b, g, e)
3105		}
3106	}
3107}
3108
3109func TestMinUint32(t *testing.T) {
3110	tests := []struct{ a, b, e uint32 }{
3111		{0, 0, 0},
3112		{0, 1, 0},
3113		{1, 0, 0},
3114
3115		{10, 10, 10},
3116		{10, 11, 10},
3117		{11, 10, 10},
3118		{11, 11, 11},
3119
3120		{math.MaxUint32, math.MaxUint32, math.MaxUint32},
3121		{math.MaxUint32, math.MaxUint32 - 1, math.MaxUint32 - 1},
3122		{math.MaxUint32 - 1, math.MaxUint32, math.MaxUint32 - 1},
3123		{math.MaxUint32 - 1, math.MaxUint32 - 1, math.MaxUint32 - 1},
3124	}
3125
3126	for _, test := range tests {
3127		if g, e := MinUint32(test.a, test.b), test.e; g != e {
3128			t.Fatal(test.a, test.b, g, e)
3129		}
3130	}
3131}
3132
3133func TestMaxUint64(t *testing.T) {
3134	tests := []struct{ a, b, e uint64 }{
3135		{0, 0, 0},
3136		{0, 1, 1},
3137		{1, 0, 1},
3138
3139		{10, 10, 10},
3140		{10, 11, 11},
3141		{11, 10, 11},
3142		{11, 11, 11},
3143
3144		{math.MaxUint64, math.MaxUint64, math.MaxUint64},
3145		{math.MaxUint64, math.MaxUint64 - 1, math.MaxUint64},
3146		{math.MaxUint64 - 1, math.MaxUint64, math.MaxUint64},
3147		{math.MaxUint64 - 1, math.MaxUint64 - 1, math.MaxUint64 - 1},
3148	}
3149
3150	for _, test := range tests {
3151		if g, e := MaxUint64(test.a, test.b), test.e; g != e {
3152			t.Fatal(test.a, test.b, g, e)
3153		}
3154	}
3155}
3156
3157func TestMinUint64(t *testing.T) {
3158	tests := []struct{ a, b, e uint64 }{
3159		{0, 0, 0},
3160		{0, 1, 0},
3161		{1, 0, 0},
3162
3163		{10, 10, 10},
3164		{10, 11, 10},
3165		{11, 10, 10},
3166		{11, 11, 11},
3167
3168		{math.MaxUint64, math.MaxUint64, math.MaxUint64},
3169		{math.MaxUint64, math.MaxUint64 - 1, math.MaxUint64 - 1},
3170		{math.MaxUint64 - 1, math.MaxUint64, math.MaxUint64 - 1},
3171		{math.MaxUint64 - 1, math.MaxUint64 - 1, math.MaxUint64 - 1},
3172	}
3173
3174	for _, test := range tests {
3175		if g, e := MinUint64(test.a, test.b), test.e; g != e {
3176			t.Fatal(test.a, test.b, g, e)
3177		}
3178	}
3179}
3180
3181func TestMaxInt8(t *testing.T) {
3182	tests := []struct{ a, b, e int8 }{
3183		{math.MinInt8, math.MinInt8, math.MinInt8},
3184		{math.MinInt8 + 1, math.MinInt8, math.MinInt8 + 1},
3185		{math.MinInt8, math.MinInt8 + 1, math.MinInt8 + 1},
3186
3187		{-1, -1, -1},
3188		{-1, 0, 0},
3189		{-1, 1, 1},
3190
3191		{0, -1, 0},
3192		{0, 0, 0},
3193		{0, 1, 1},
3194
3195		{1, -1, 1},
3196		{1, 0, 1},
3197		{1, 1, 1},
3198
3199		{math.MaxInt8, math.MaxInt8, math.MaxInt8},
3200		{math.MaxInt8 - 1, math.MaxInt8, math.MaxInt8},
3201		{math.MaxInt8, math.MaxInt8 - 1, math.MaxInt8},
3202	}
3203
3204	for _, test := range tests {
3205		if g, e := MaxInt8(test.a, test.b), test.e; g != e {
3206			t.Fatal(test.a, test.b, g, e)
3207		}
3208	}
3209}
3210
3211func TestMinInt8(t *testing.T) {
3212	tests := []struct{ a, b, e int8 }{
3213		{math.MinInt8, math.MinInt8, math.MinInt8},
3214		{math.MinInt8 + 1, math.MinInt8, math.MinInt8},
3215		{math.MinInt8, math.MinInt8 + 1, math.MinInt8},
3216
3217		{-1, -1, -1},
3218		{-1, 0, -1},
3219		{-1, 1, -1},
3220
3221		{0, -1, -1},
3222		{0, 0, 0},
3223		{0, 1, 0},
3224
3225		{1, -1, -1},
3226		{1, 0, 0},
3227		{1, 1, 1},
3228
3229		{math.MaxInt8, math.MaxInt8, math.MaxInt8},
3230		{math.MaxInt8 - 1, math.MaxInt8, math.MaxInt8 - 1},
3231		{math.MaxInt8, math.MaxInt8 - 1, math.MaxInt8 - 1},
3232	}
3233
3234	for _, test := range tests {
3235		if g, e := MinInt8(test.a, test.b), test.e; g != e {
3236			t.Fatal(test.a, test.b, g, e)
3237		}
3238	}
3239}
3240
3241func TestMaxInt16(t *testing.T) {
3242	tests := []struct{ a, b, e int16 }{
3243		{math.MinInt16, math.MinInt16, math.MinInt16},
3244		{math.MinInt16 + 1, math.MinInt16, math.MinInt16 + 1},
3245		{math.MinInt16, math.MinInt16 + 1, math.MinInt16 + 1},
3246
3247		{-1, -1, -1},
3248		{-1, 0, 0},
3249		{-1, 1, 1},
3250
3251		{0, -1, 0},
3252		{0, 0, 0},
3253		{0, 1, 1},
3254
3255		{1, -1, 1},
3256		{1, 0, 1},
3257		{1, 1, 1},
3258
3259		{math.MaxInt16, math.MaxInt16, math.MaxInt16},
3260		{math.MaxInt16 - 1, math.MaxInt16, math.MaxInt16},
3261		{math.MaxInt16, math.MaxInt16 - 1, math.MaxInt16},
3262	}
3263
3264	for _, test := range tests {
3265		if g, e := MaxInt16(test.a, test.b), test.e; g != e {
3266			t.Fatal(test.a, test.b, g, e)
3267		}
3268	}
3269}
3270
3271func TestMinInt16(t *testing.T) {
3272	tests := []struct{ a, b, e int16 }{
3273		{math.MinInt16, math.MinInt16, math.MinInt16},
3274		{math.MinInt16 + 1, math.MinInt16, math.MinInt16},
3275		{math.MinInt16, math.MinInt16 + 1, math.MinInt16},
3276
3277		{-1, -1, -1},
3278		{-1, 0, -1},
3279		{-1, 1, -1},
3280
3281		{0, -1, -1},
3282		{0, 0, 0},
3283		{0, 1, 0},
3284
3285		{1, -1, -1},
3286		{1, 0, 0},
3287		{1, 1, 1},
3288
3289		{math.MaxInt16, math.MaxInt16, math.MaxInt16},
3290		{math.MaxInt16 - 1, math.MaxInt16, math.MaxInt16 - 1},
3291		{math.MaxInt16, math.MaxInt16 - 1, math.MaxInt16 - 1},
3292	}
3293
3294	for _, test := range tests {
3295		if g, e := MinInt16(test.a, test.b), test.e; g != e {
3296			t.Fatal(test.a, test.b, g, e)
3297		}
3298	}
3299}
3300
3301func TestMaxInt32(t *testing.T) {
3302	tests := []struct{ a, b, e int32 }{
3303		{math.MinInt32, math.MinInt32, math.MinInt32},
3304		{math.MinInt32 + 1, math.MinInt32, math.MinInt32 + 1},
3305		{math.MinInt32, math.MinInt32 + 1, math.MinInt32 + 1},
3306
3307		{-1, -1, -1},
3308		{-1, 0, 0},
3309		{-1, 1, 1},
3310
3311		{0, -1, 0},
3312		{0, 0, 0},
3313		{0, 1, 1},
3314
3315		{1, -1, 1},
3316		{1, 0, 1},
3317		{1, 1, 1},
3318
3319		{math.MaxInt32, math.MaxInt32, math.MaxInt32},
3320		{math.MaxInt32 - 1, math.MaxInt32, math.MaxInt32},
3321		{math.MaxInt32, math.MaxInt32 - 1, math.MaxInt32},
3322	}
3323
3324	for _, test := range tests {
3325		if g, e := MaxInt32(test.a, test.b), test.e; g != e {
3326			t.Fatal(test.a, test.b, g, e)
3327		}
3328	}
3329}
3330
3331func TestMinInt32(t *testing.T) {
3332	tests := []struct{ a, b, e int32 }{
3333		{math.MinInt32, math.MinInt32, math.MinInt32},
3334		{math.MinInt32 + 1, math.MinInt32, math.MinInt32},
3335		{math.MinInt32, math.MinInt32 + 1, math.MinInt32},
3336
3337		{-1, -1, -1},
3338		{-1, 0, -1},
3339		{-1, 1, -1},
3340
3341		{0, -1, -1},
3342		{0, 0, 0},
3343		{0, 1, 0},
3344
3345		{1, -1, -1},
3346		{1, 0, 0},
3347		{1, 1, 1},
3348
3349		{math.MaxInt32, math.MaxInt32, math.MaxInt32},
3350		{math.MaxInt32 - 1, math.MaxInt32, math.MaxInt32 - 1},
3351		{math.MaxInt32, math.MaxInt32 - 1, math.MaxInt32 - 1},
3352	}
3353
3354	for _, test := range tests {
3355		if g, e := MinInt32(test.a, test.b), test.e; g != e {
3356			t.Fatal(test.a, test.b, g, e)
3357		}
3358	}
3359}
3360
3361func TestMaxInt64(t *testing.T) {
3362	tests := []struct{ a, b, e int64 }{
3363		{math.MinInt64, math.MinInt64, math.MinInt64},
3364		{math.MinInt64 + 1, math.MinInt64, math.MinInt64 + 1},
3365		{math.MinInt64, math.MinInt64 + 1, math.MinInt64 + 1},
3366
3367		{-1, -1, -1},
3368		{-1, 0, 0},
3369		{-1, 1, 1},
3370
3371		{0, -1, 0},
3372		{0, 0, 0},
3373		{0, 1, 1},
3374
3375		{1, -1, 1},
3376		{1, 0, 1},
3377		{1, 1, 1},
3378
3379		{math.MaxInt64, math.MaxInt64, math.MaxInt64},
3380		{math.MaxInt64 - 1, math.MaxInt64, math.MaxInt64},
3381		{math.MaxInt64, math.MaxInt64 - 1, math.MaxInt64},
3382	}
3383
3384	for _, test := range tests {
3385		if g, e := MaxInt64(test.a, test.b), test.e; g != e {
3386			t.Fatal(test.a, test.b, g, e)
3387		}
3388	}
3389}
3390
3391func TestMinInt64(t *testing.T) {
3392	tests := []struct{ a, b, e int64 }{
3393		{math.MinInt64, math.MinInt64, math.MinInt64},
3394		{math.MinInt64 + 1, math.MinInt64, math.MinInt64},
3395		{math.MinInt64, math.MinInt64 + 1, math.MinInt64},
3396
3397		{-1, -1, -1},
3398		{-1, 0, -1},
3399		{-1, 1, -1},
3400
3401		{0, -1, -1},
3402		{0, 0, 0},
3403		{0, 1, 0},
3404
3405		{1, -1, -1},
3406		{1, 0, 0},
3407		{1, 1, 1},
3408
3409		{math.MaxInt64, math.MaxInt64, math.MaxInt64},
3410		{math.MaxInt64 - 1, math.MaxInt64, math.MaxInt64 - 1},
3411		{math.MaxInt64, math.MaxInt64 - 1, math.MaxInt64 - 1},
3412	}
3413
3414	for _, test := range tests {
3415		if g, e := MinInt64(test.a, test.b), test.e; g != e {
3416			t.Fatal(test.a, test.b, g, e)
3417		}
3418	}
3419}
3420
3421func TestPopCountBigInt(t *testing.T) {
3422	const N = 1e4
3423	rng := rand.New(rand.NewSource(42))
3424	lim := big.NewInt(0)
3425	lim.SetBit(lim, 1e3, 1)
3426	z := big.NewInt(0)
3427	m1 := big.NewInt(-1)
3428	for i := 0; i < N; i++ {
3429		z.Rand(rng, lim)
3430		g := PopCountBigInt(z)
3431		e := 0
3432		for bit := 0; bit < z.BitLen(); bit++ {
3433			if z.Bit(bit) != 0 {
3434				e++
3435			}
3436		}
3437		if g != e {
3438			t.Fatal(g, e)
3439		}
3440
3441		z.Mul(z, m1)
3442		if g := PopCountBigInt(z); g != e {
3443			t.Fatal(g, e)
3444		}
3445	}
3446}
3447
3448func benchmarkPopCountBigInt(b *testing.B, bits int) {
3449	lim := big.NewInt(0)
3450	lim.SetBit(lim, bits, 1)
3451	z := big.NewInt(0)
3452	z.Rand(rand.New(rand.NewSource(42)), lim)
3453	b.ResetTimer()
3454	for i := 0; i < b.N; i++ {
3455		PopCountBigInt(z)
3456	}
3457}
3458
3459func BenchmarkPopCountBigInt1e1(b *testing.B) {
3460	benchmarkPopCountBigInt(b, 1e1)
3461}
3462
3463func BenchmarkPopCountBigInt1e2(b *testing.B) {
3464	benchmarkPopCountBigInt(b, 1e2)
3465}
3466
3467func BenchmarkPopCountBigInt1e3(b *testing.B) {
3468	benchmarkPopCountBigInt(b, 1e3)
3469}
3470
3471func BenchmarkPopCountBigIbnt1e4(b *testing.B) {
3472	benchmarkPopCountBigInt(b, 1e4)
3473}
3474
3475func BenchmarkPopCountBigInt1e5(b *testing.B) {
3476	benchmarkPopCountBigInt(b, 1e5)
3477}
3478
3479func BenchmarkPopCountBigInt1e6(b *testing.B) {
3480	benchmarkPopCountBigInt(b, 1e6)
3481}
3482
3483func TestToBase(t *testing.T) {
3484	x := ToBase(big.NewInt(0), 42)
3485	e := []int{0}
3486	if g, e := len(x), len(e); g != e {
3487		t.Fatal(g, e)
3488	}
3489
3490	for i, g := range x {
3491		if e := e[i]; g != e {
3492			t.Fatal(i, g, e)
3493		}
3494	}
3495
3496	x = ToBase(big.NewInt(2047), 22)
3497	e = []int{1, 5, 4}
3498	if g, e := len(x), len(e); g != e {
3499		t.Fatal(g, e)
3500	}
3501
3502	for i, g := range x {
3503		if e := e[i]; g != e {
3504			t.Fatal(i, g, e)
3505		}
3506	}
3507
3508	x = ToBase(big.NewInt(-2047), 22)
3509	e = []int{-1, -5, -4}
3510	if g, e := len(x), len(e); g != e {
3511		t.Fatal(g, e)
3512	}
3513
3514	for i, g := range x {
3515		if e := e[i]; g != e {
3516			t.Fatal(i, g, e)
3517		}
3518	}
3519}
3520
3521func TestBug(t *testing.T) {
3522	if BitLenUint(MaxUint) != 64 {
3523		t.Logf("Bug reproducible only on 64 bit architecture")
3524		return
3525	}
3526
3527	_, err := NewFC32(MinInt, MaxInt, true)
3528	if err == nil {
3529		t.Fatal("Expected non nil err")
3530	}
3531}
3532
3533func poly(a ...int) string {
3534	var b bytes.Buffer
3535	for i, v := range a {
3536		p := len(a) - i - 1
3537		if v == 0 && p != 0 {
3538			continue
3539		}
3540
3541		if v == 0 && p == 0 && b.Len() != 0 {
3542			continue
3543		}
3544
3545		if av := abs(v); av == 1 && p != 0 {
3546			if b.Len() != 0 {
3547				if v == 1 {
3548					b.WriteByte('+')
3549				} else {
3550					b.WriteByte('-')
3551				}
3552			} else if v == -1 {
3553				b.WriteByte('-')
3554			}
3555		} else {
3556			switch {
3557			case b.Len() == 0:
3558				fmt.Fprintf(&b, "%d", v)
3559			default:
3560				fmt.Fprintf(&b, "%+d", v)
3561			}
3562		}
3563
3564		if p == 0 {
3565			continue
3566		}
3567
3568		if p == 1 {
3569			fmt.Fprintf(&b, "x")
3570			continue
3571		}
3572
3573		fmt.Fprintf(&b, "x^%d", p)
3574	}
3575	return b.String()
3576}
3577
3578func polyK(k int) string {
3579	switch {
3580	case k == -1:
3581		return "-"
3582	case k == 1:
3583		return ""
3584	default:
3585		return fmt.Sprint(k)
3586	}
3587}
3588
3589func TestQuadPolyDiscriminant(t *testing.T) {
3590	for i, test := range []struct {
3591		a, b, c, ds, d int
3592	}{
3593		{-1, -5, 6, 49, 7},
3594		{-1, 5, 6, 49, 7},
3595		{1, -5, -6, 49, 7},
3596		{1, 5, -6, 49, 7},
3597		{1, 5, 6, 1, 1},
3598		{2, 3, 5, -31, -1},
3599		{2, 7, 3, 25, 5},
3600		{3, 8, 5, 4, 2},
3601		{3, 9, 5, 21, -1},
3602		{4, 5, 1, 9, 3},
3603		{5, 3, 2, -31, -1},
3604	} {
3605		ds, d, err := QuadPolyDiscriminant(test.a, test.b, test.c)
3606		if err != nil {
3607			t.Fatal(i, err)
3608		}
3609
3610		if g, e := ds, test.ds; g != e {
3611			t.Fatal(i, g, e)
3612		}
3613
3614		if g, e := d, test.d; g != e {
3615			t.Fatal(i, g, e)
3616		}
3617	}
3618}
3619
3620func testQuadPolyFactors(t *testing.T, p1, q1, p2, q2, k, cases int) {
3621	a := k * p1 * p2
3622	b := k * (p1*q2 + q1*p2)
3623	c := k * (q1 * q2)
3624	con, f, err := QuadPolyFactors(a, b, c)
3625	if err != nil {
3626		t.Fatalf(
3627			"%d: %s(%s)(%s) = %s -> %v",
3628			cases, polyK(k), poly(p1, q1), poly(p2, q2), poly(a, b, c),
3629			err,
3630		)
3631	}
3632
3633	switch {
3634	case a == 0:
3635		if g, e := len(f), 1; g != e {
3636			t.Fatalf(
3637				"%d: %s(%s)(%s) = %s -> got %v factors, expected %v",
3638				cases, polyK(k), poly(p1, q1), poly(p2, q2), poly(a, b, c),
3639				g, e,
3640			)
3641		}
3642
3643		a2 := 0
3644		b2 := con * f[0].P
3645		c2 := con * f[0].Q
3646		if a != a2 || b != b2 || c != c2 {
3647			t.Fatalf(
3648				"%d: %s(%s)(%s) = %s -> %s(%s) = %s",
3649				cases, polyK(k), poly(p1, q1), poly(p2, q2), poly(a, b, c),
3650				polyK(con), poly(f[0].P, f[0].Q), poly(a2, b2, c2),
3651			)
3652		}
3653
3654		t.Logf(
3655			"%d: %s(%s)(%s) = %s -> %s(%s) = %s",
3656			cases, polyK(k), poly(p1, q1), poly(p2, q2), poly(a, b, c),
3657			polyK(con), poly(f[0].P, f[0].Q), poly(a2, b2, c2),
3658		)
3659	default:
3660		if g, e := len(f), 2; g != e {
3661			t.Fatalf(
3662				"%d: %s(%s)(%s) = %s -> got %v factors, expected %v",
3663				cases, polyK(k), poly(p1, q1), poly(p2, q2), poly(a, b, c),
3664				g, e,
3665			)
3666		}
3667
3668		a2 := con * f[0].P * f[1].P
3669		b2 := con * (f[0].P*f[1].Q + f[0].Q*f[1].P)
3670		c2 := con * f[0].Q * f[1].Q
3671		if a != a2 || b != b2 || c != c2 {
3672			dbg("", con, f)
3673			t.Fatalf(
3674				"%d: %s(%s)(%s) = %s -> %s(%s)(%s) = %s",
3675				cases, polyK(k), poly(p1, q1), poly(p2, q2), poly(a, b, c),
3676				polyK(con), poly(f[0].P, f[0].Q), poly(f[1].P, f[1].Q), poly(a2, b2, c2),
3677			)
3678		}
3679
3680		t.Logf(
3681			"%d: %s(%s)(%s) = %s -> %s(%s)(%s) = %s",
3682			cases, polyK(k), poly(p1, q1), poly(p2, q2), poly(a, b, c),
3683			polyK(con), poly(f[0].P, f[0].Q), poly(f[1].P, f[1].Q), poly(a2, b2, c2),
3684		)
3685	}
3686}
3687
3688func TestQuadPolyFactors(t *testing.T) {
3689	cases := 0
3690
3691	const N = 1e4
3692	mask := 1<<14 - 1
3693	if IntBits < 64 {
3694		mask = 1<<7 - 1
3695	}
3696	rng := rand.New(rand.NewSource(42))
3697	for i := 0; i < N; i++ {
3698		p1 := int(rng.Int63()) & mask
3699		q1 := int(rng.Int63()) & mask
3700		p2 := int(rng.Int63()) & mask
3701		q2 := int(rng.Int63()) & mask
3702		k := int(rng.Int63()) & mask
3703		testQuadPolyFactors(t, p1, q1, p2, q2, k, cases)
3704		cases++
3705	}
3706
3707	cons := []int{-1, 1}
3708	const lim = 7
3709	for p1 := -lim; p1 <= lim; p1++ {
3710		for q1 := -lim; q1 <= lim; q1++ {
3711			for p2 := -lim; p2 <= lim; p2++ {
3712				for q2 := -lim; q2 <= lim; q2++ {
3713					for _, k := range cons {
3714						testQuadPolyFactors(t, p1, q1, p2, q2, k, cases)
3715						cases++
3716					}
3717				}
3718			}
3719		}
3720	}
3721}
3722