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