336 lines
6.6 KiB
Go
336 lines
6.6 KiB
Go
// 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 (
|
|
"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.
|
|
//
|
|
//TODO rename to IsPrimeUint32
|
|
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.
|
|
//
|
|
//TODO rename to NextPrimeUint32
|
|
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
|
|
}
|