diff options
Diffstat (limited to 'vendor/modernc.org/mathutil/primes.go')
-rw-r--r-- | vendor/modernc.org/mathutil/primes.go | 331 |
1 files changed, 331 insertions, 0 deletions
diff --git a/vendor/modernc.org/mathutil/primes.go b/vendor/modernc.org/mathutil/primes.go new file mode 100644 index 00000000..ec1e5f5a --- /dev/null +++ b/vendor/modernc.org/mathutil/primes.go @@ -0,0 +1,331 @@ +// Copyright (c) 2014 The mathutil Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package mathutil // import "modernc.org/mathutil" + +import ( + "math" +) + +// IsPrimeUint16 returns true if n is prime. Typical run time is few ns. +func IsPrimeUint16(n uint16) bool { + return n > 0 && primes16[n-1] == 1 +} + +// NextPrimeUint16 returns first prime > n and true if successful or an +// undefined value and false if there is no next prime in the uint16 limits. +// Typical run time is few ns. +func NextPrimeUint16(n uint16) (p uint16, ok bool) { + return n + uint16(primes16[n]), n < 65521 +} + +// IsPrime returns true if n is prime. Typical run time is about 100 ns. +func IsPrime(n uint32) bool { + switch { + case n&1 == 0: + return n == 2 + case n%3 == 0: + return n == 3 + case n%5 == 0: + return n == 5 + case n%7 == 0: + return n == 7 + case n%11 == 0: + return n == 11 + case n%13 == 0: + return n == 13 + case n%17 == 0: + return n == 17 + case n%19 == 0: + return n == 19 + case n%23 == 0: + return n == 23 + case n%29 == 0: + return n == 29 + case n%31 == 0: + return n == 31 + case n%37 == 0: + return n == 37 + case n%41 == 0: + return n == 41 + case n%43 == 0: + return n == 43 + case n%47 == 0: + return n == 47 + case n%53 == 0: + return n == 53 // Benchmarked optimum + case n < 65536: + // use table data + return IsPrimeUint16(uint16(n)) + default: + mod := ModPowUint32(2, (n+1)/2, n) + if mod != 2 && mod != n-2 { + return false + } + blk := &lohi[n>>24] + lo, hi := blk.lo, blk.hi + for lo <= hi { + index := (lo + hi) >> 1 + liar := liars[index] + switch { + case n > liar: + lo = index + 1 + case n < liar: + hi = index - 1 + default: + return false + } + } + return true + } +} + +// IsPrimeUint64 returns true if n is prime. Typical run time is few tens of µs. +// +// SPRP bases: http://miller-rabin.appspot.com +func IsPrimeUint64(n uint64) bool { + switch { + case n%2 == 0: + return n == 2 + case n%3 == 0: + return n == 3 + case n%5 == 0: + return n == 5 + case n%7 == 0: + return n == 7 + case n%11 == 0: + return n == 11 + case n%13 == 0: + return n == 13 + case n%17 == 0: + return n == 17 + case n%19 == 0: + return n == 19 + case n%23 == 0: + return n == 23 + case n%29 == 0: + return n == 29 + case n%31 == 0: + return n == 31 + case n%37 == 0: + return n == 37 + case n%41 == 0: + return n == 41 + case n%43 == 0: + return n == 43 + case n%47 == 0: + return n == 47 + case n%53 == 0: + return n == 53 + case n%59 == 0: + return n == 59 + case n%61 == 0: + return n == 61 + case n%67 == 0: + return n == 67 + case n%71 == 0: + return n == 71 + case n%73 == 0: + return n == 73 + case n%79 == 0: + return n == 79 + case n%83 == 0: + return n == 83 + case n%89 == 0: + return n == 89 // Benchmarked optimum + case n <= math.MaxUint16: + return IsPrimeUint16(uint16(n)) + case n <= math.MaxUint32: + return ProbablyPrimeUint32(uint32(n), 11000544) && + ProbablyPrimeUint32(uint32(n), 31481107) + case n < 105936894253: + return ProbablyPrimeUint64_32(n, 2) && + ProbablyPrimeUint64_32(n, 1005905886) && + ProbablyPrimeUint64_32(n, 1340600841) + case n < 31858317218647: + return ProbablyPrimeUint64_32(n, 2) && + ProbablyPrimeUint64_32(n, 642735) && + ProbablyPrimeUint64_32(n, 553174392) && + ProbablyPrimeUint64_32(n, 3046413974) + case n < 3071837692357849: + return ProbablyPrimeUint64_32(n, 2) && + ProbablyPrimeUint64_32(n, 75088) && + ProbablyPrimeUint64_32(n, 642735) && + ProbablyPrimeUint64_32(n, 203659041) && + ProbablyPrimeUint64_32(n, 3613982119) + default: + return ProbablyPrimeUint64_32(n, 2) && + ProbablyPrimeUint64_32(n, 325) && + ProbablyPrimeUint64_32(n, 9375) && + ProbablyPrimeUint64_32(n, 28178) && + ProbablyPrimeUint64_32(n, 450775) && + ProbablyPrimeUint64_32(n, 9780504) && + ProbablyPrimeUint64_32(n, 1795265022) + } +} + +// NextPrime returns first prime > n and true if successful or an undefined value and false if there +// is no next prime in the uint32 limits. Typical run time is about 2 µs. +func NextPrime(n uint32) (p uint32, ok bool) { + switch { + case n < 65521: + p16, _ := NextPrimeUint16(uint16(n)) + return uint32(p16), true + case n >= math.MaxUint32-4: + return + } + + n++ + var d0, d uint32 + switch mod := n % 6; mod { + case 0: + d0, d = 1, 4 + case 1: + d = 4 + case 2, 3, 4: + d0, d = 5-mod, 2 + case 5: + d = 2 + } + + p = n + d0 + if p < n { // overflow + return + } + + for { + if IsPrime(p) { + return p, true + } + + p0 := p + p += d + if p < p0 { // overflow + break + } + + d ^= 6 + } + return +} + +// NextPrimeUint64 returns first prime > n and true if successful or an undefined value and false if there +// is no next prime in the uint64 limits. Typical run time is in hundreds of µs. +func NextPrimeUint64(n uint64) (p uint64, ok bool) { + switch { + case n < 65521: + p16, _ := NextPrimeUint16(uint16(n)) + return uint64(p16), true + case n >= 18446744073709551557: // last uint64 prime + return + } + + n++ + var d0, d uint64 + switch mod := n % 6; mod { + case 0: + d0, d = 1, 4 + case 1: + d = 4 + case 2, 3, 4: + d0, d = 5-mod, 2 + case 5: + d = 2 + } + + p = n + d0 + if p < n { // overflow + return + } + + for { + if ok = IsPrimeUint64(p); ok { + break + } + + p0 := p + p += d + if p < p0 { // overflow + break + } + + d ^= 6 + } + return +} + +// FactorTerm is one term of an integer factorization. +type FactorTerm struct { + Prime uint32 // The divisor + Power uint32 // Term == Prime^Power +} + +// FactorTerms represent a factorization of an integer +type FactorTerms []FactorTerm + +// FactorInt returns prime factorization of n > 1 or nil otherwise. +// Resulting factors are ordered by Prime. Typical run time is few µs. +func FactorInt(n uint32) (f FactorTerms) { + switch { + case n < 2: + return + case IsPrime(n): + return []FactorTerm{{n, 1}} + } + + f, w := make([]FactorTerm, 9), 0 + for p := 2; p < len(primes16); p += int(primes16[p]) { + if uint(p*p) > uint(n) { + break + } + + power := uint32(0) + for n%uint32(p) == 0 { + n /= uint32(p) + power++ + } + if power != 0 { + f[w] = FactorTerm{uint32(p), power} + w++ + } + if n == 1 { + break + } + } + if n != 1 { + f[w] = FactorTerm{n, 1} + w++ + } + return f[:w] +} + +// PrimorialProductsUint32 returns a slice of numbers in [lo, hi] which are a +// product of max 'max' primorials. The slice is not sorted. +// +// See also: http://en.wikipedia.org/wiki/Primorial +func PrimorialProductsUint32(lo, hi, max uint32) (r []uint32) { + lo64, hi64 := int64(lo), int64(hi) + if max > 31 { // N/A + max = 31 + } + + var f func(int64, int64, uint32) + f = func(n, p int64, emax uint32) { + e := uint32(1) + for n <= hi64 && e <= emax { + n *= p + if n >= lo64 && n <= hi64 { + r = append(r, uint32(n)) + } + if n < hi64 { + p, _ := NextPrime(uint32(p)) + f(n, int64(p), e) + } + e++ + } + } + + f(1, 2, max) + return +} |