Added primes package.

primes/cache.go

@@ -0,0 +1,158 @@
+package primes
+
+import "opslag.de/schobers/geom/ints"
+
+var primes = &cache{}
+
+const (
+	primesPerPage int = 8192
+)
+
+// Cache represents the cached state of some pre-calculated prime numbers.
+type Cache interface {
+	IsPrime(n int64) bool
+	OffsetOfPrime(n int64) int
+	Prime(n int) int64
+}
+
+// NewCache creates a new prime number cache.
+func NewCache() Cache {
+	return &cache{}
+}
+
+type cache struct {
+	pages []page
+}
+
+type page struct {
+	primes []int64
+	last   int64
+}
+
+func (c *cache) OffsetOfPrime(n int64) int {
+	p := 0
+	for {
+		if len(c.pages) == p {
+			c.generatePage()
+		}
+		if c.pages[p].last >= n {
+			for i := 0; i < primesPerPage; i++ {
+				if c.pages[p].primes[i] >= n {
+					return i + p*primesPerPage
+				}
+			}
+		}
+		p++
+	}
+}
+
+func (c *cache) Prime(n int) int64 {
+	i := n - 1
+	p := i / primesPerPage
+	for p >= len(c.pages) {
+		c.generatePage()
+	}
+	offset := i % primesPerPage
+	return c.pages[p].primes[offset]
+}
+
+func (c *cache) IsPrime(n int64) bool {
+	if inRange, isPrime := c.isKnownPrime(n); inRange {
+		return isPrime
+	}
+	upper := ints.Sqrt64(n)
+	it := NewIteratorFromCache(c)
+	for it.Next() {
+		p := it.Get()
+		if p > upper {
+			return true
+		}
+		if n%p == 0 {
+			return false
+		}
+	}
+	panic("out of primes")
+}
+
+func (c *cache) generatePage() {
+	pageI := len(c.pages)
+	c.pages = append(c.pages, page{})
+	var p = &c.pages[pageI]
+	p.primes = make([]int64, primesPerPage)
+	var i int
+	var last int64
+	if 0 == pageI {
+		copy(p.primes, seed)
+		i = 20
+		last = p.primes[19]
+	} else {
+		last = c.pages[pageI-1].last
+	}
+
+	wheel := newWheel6(c)
+	wheel.resetTo(last)
+
+	for ; i < primesPerPage; i++ {
+		prime := wheel.next()
+		for !c.IsPrime(prime) {
+			prime = wheel.next()
+		}
+		p.primes[i] = prime
+	}
+	p.last = p.primes[primesPerPage-1]
+}
+
+func (c *cache) isKnownPrime(n int64) (bool, bool) {
+	if 0 == len(c.pages) {
+		c.generatePage()
+	}
+	for _, p := range c.pages {
+		if n > p.last {
+			continue
+		}
+		l := len(p.primes)
+		offset := 0
+		shift := l >> 1
+		for shift > 0 {
+			pp := p.primes[offset+shift-1]
+			if n == pp {
+				return true, true
+			}
+			if n > pp {
+				offset += shift
+			}
+			shift >>= 1
+		}
+		if n == p.last {
+			return true, true
+		}
+		return true, false
+	}
+	return false, false
+}
+
+// Generate n pre-calculates the first n prime numbers for the default cache.
+func Generate(n int) {
+	primes.Prime(n)
+}
+
+// N returns the first n prime numbers.
+func N(n int) []int64 {
+	Generate(n)
+
+	ps := make([]int64, n)
+	var i = 0
+	for _, p := range primes.pages {
+		var toCopy = primesPerPage
+		if n-i < primesPerPage {
+			toCopy = n - i
+		}
+		i += copy(ps[i:], p.primes[:toCopy])
+	}
+	return ps
+}
+
+// Nth returns the nth prime number.
+func Nth(n int) int64 {
+	return primes.Prime(n)
+}

primes/cache_test.go

@@ -0,0 +1,65 @@
+package primes
+
+import (
+	"testing"
+
+	"github.com/stretchr/testify/assert"
+)
+
+func TestIsKnownPrime(t *testing.T) {
+	primes.Prime(100)
+	var known = []struct {
+		n       int64
+		inRange bool
+		isPrime bool
+	}{
+		{2, true, true}, // prime
+		{3, true, true}, // prime
+		{4, true, false},
+		{5, true, true}, // prime
+		{52, true, false},
+		{53, true, true}, // prime
+		{54, true, false},
+		{58, true, false},
+		{59, true, true}, // prime
+		{60, true, false},
+		{1000000, false, false}, // not cached
+	}
+	for _, k := range known {
+		inRange, isPrime := primes.isKnownPrime(k.n)
+		assert.Equal(t, k.inRange, inRange)
+		assert.Equal(t, k.isPrime, isPrime)
+	}
+}
+
+func TestIsPrime(t *testing.T) {
+	assert.True(t, IsPrime(84017))
+}
+
+func TestPrimeN(t *testing.T) {
+	assert.Equal(t, int64(2), NewCache().Prime(1))
+	assert.Equal(t, int64(541), NewCache().Prime(100))
+	assert.Equal(t, int64(1987), NewCache().Prime(300))
+	assert.Equal(t, int64(7919), NewCache().Prime(1000))
+	assert.Equal(t, int64(84017), NewCache().Prime(8192))
+	assert.Equal(t, int64(104729), NewCache().Prime(10000))
+}
+
+func TestIterator(t *testing.T) {
+	nth := 1
+	reference := NewCache()
+	assert.Equal(t, int64(104729), reference.Prime(10000))
+	c := NewCache()
+	it := NewIteratorFromCache(c)
+	for it.Next() && nth < 10000 {
+		assert.Equal(t, reference.Prime(nth), it.Get())
+		nth++
+	}
+}
+
+func BenchmarkGenerateCache(b *testing.B) {
+	var c = &cache{}
+	for i := 0; i < b.N; i++ {
+		c.Prime(5000)
+	}
+}

primes/factor.go

@@ -0,0 +1,92 @@
+package primes
+
+import "math/big"
+
+// Factor represents a factor with a base (Prime) and an exponent (Exponent)
+type Factor struct {
+	Prime    int64
+	Exponent int
+}
+
+// Factorize calculates the factors of n and returns it as a slice of Factor's (uses a default cache).
+func Factorize(n int64) []Factor {
+	return FactorizeCache(primes, n)
+}
+
+// FactorizeCache calculates the factors of n and returns it as a slice of Factor's and specifies the cache to use.
+func FactorizeCache(c Cache, n int64) []Factor {
+	if 1 > n {
+		return []Factor{}
+	}
+	var factors []Factor
+	it := NewIteratorFromCache(c)
+	for it.Next() {
+		p := it.Get()
+		var exp int
+		for n%p == 0 {
+			n /= p
+			exp++
+		}
+		if 0 < exp {
+			factors = append(factors, Factor{p, exp})
+		}
+		if n < p*p {
+			break
+		}
+	}
+	if 1 < n {
+		factors = append(factors, Factor{n, 1})
+	} else if 1 == n && 0 == len(factors) {
+		return []Factor{{1, 1}}
+	}
+	return factors
+}
+
+// MustBigIntToInt64 returns the int64 representation of the big int n. If n is too large the method panics.
+func MustBigIntToInt64(n *big.Int) int64 {
+	if !n.IsInt64() {
+		panic("n can't be represented as a int64")
+	}
+	return n.Int64()
+}
+
+// FactorizeBigInt calculates the factors of n and returns it as a slice of Factor's. If the base of a factor is larger than can be represented as a int64 the method will panic.
+func FactorizeBigInt(n *big.Int) []Factor {
+	n = big.NewInt(0).Set(n)
+	var one = big.NewInt(1)
+	if one.Cmp(n) > 0 {
+		return []Factor{}
+	}
+	var zero = big.NewInt(0)
+	var divModZero = func(p *big.Int) bool {
+		var quo, rem = big.NewInt(0), big.NewInt(0)
+		quo.QuoRem(n, p, rem)
+		if rem.Cmp(zero) == 0 {
+			n = quo
+			return true
+		}
+		return false
+	}
+	var factors []Factor
+	it := NewIterator()
+	for it.Next() {
+		var p = big.NewInt(it.Get())
+		var exp int
+		for divModZero(p) {
+			exp++
+		}
+		if 0 < exp {
+			factors = append(factors, Factor{it.Get(), exp})
+		}
+		p.Mul(p, p)
+		if n.Cmp(p) < 0 {
+			break
+		}
+	}
+	if one.Cmp(n) < 0 {
+		factors = append(factors, Factor{MustBigIntToInt64(n), 1})
+	} else if one.Cmp(n) == 0 && 0 == len(factors) {
+		return []Factor{{1, 1}}
+	}
+	return factors
+}

primes/factor_test.go

@@ -0,0 +1,36 @@
+package primes
+
+import (
+	"math/big"
+	"testing"
+
+	"github.com/stretchr/testify/assert"
+)
+
+func TestFactorize(t *testing.T) {
+	assert.Equal(t, []Factor{}, Factorize(-1))
+	assert.Equal(t, []Factor{}, Factorize(0))
+	assert.Equal(t, []Factor{{1, 1}}, Factorize(1))
+	assert.Equal(t, []Factor{{2, 1}, {5, 1}}, Factorize(10))
+	assert.Equal(t, []Factor{{3, 1}, {7, 2}}, Factorize(147))
+	assert.Equal(t, []Factor{{3, 1}, {7, 1}, {11, 1}, {13, 1}, {37, 1}}, Factorize(111111))
+	assert.Equal(t, []Factor{{7, 2}, {73, 1}, {127, 1}, {337, 1}, {92737, 1}, {649657, 1}}, Factorize(9223372036854775807))
+	// assert.Equal(t, []Factor{{11, 1}, {41, 1}, {101, 1}, {271, 1}, {3541, 1}, {9091, 1}, {27961, 1}}, Factorize(11111111111111111111))
+}
+
+func stringToBigInt(s string) *big.Int {
+	var i = big.NewInt(0)
+	i, _ = i.SetString(s, 10)
+	return i
+}
+
+func TestFactorizeBigInt(t *testing.T) {
+	assert.Equal(t, []Factor{}, FactorizeBigInt(big.NewInt(-1)))
+	assert.Equal(t, []Factor{}, FactorizeBigInt(big.NewInt(0)))
+	assert.Equal(t, []Factor{{1, 1}}, FactorizeBigInt(big.NewInt(1)))
+	assert.Equal(t, []Factor{{2, 1}, {5, 1}}, FactorizeBigInt(big.NewInt(10)))
+	assert.Equal(t, []Factor{{3, 1}, {7, 2}}, FactorizeBigInt(big.NewInt(147)))
+	assert.Equal(t, []Factor{{3, 1}, {7, 1}, {11, 1}, {13, 1}, {37, 1}}, FactorizeBigInt(big.NewInt(111111)))
+	assert.Equal(t, []Factor{{7, 2}, {73, 1}, {127, 1}, {337, 1}, {92737, 1}, {649657, 1}}, FactorizeBigInt(big.NewInt(9223372036854775807)))
+	assert.Equal(t, []Factor{{11, 1}, {41, 1}, {101, 1}, {271, 1}, {3541, 1}, {9091, 1}, {27961, 1}}, FactorizeBigInt(stringToBigInt("11111111111111111111")))
+}

primes/ints.go

@@ -0,0 +1,81 @@
+package primes
+
+import "opslag.de/schobers/geom/ints"
+
+// GreatestCommonDivisor returns the greatest common divisor of a and b.
+func GreatestCommonDivisor(a, b int64) int64 { return GCD(a, b) }
+
+// GCD returns the greatest common divisor of a and b.
+func GCD(a, b int64) int64 {
+	factors := Factorize(a)
+	var i int
+	var gcd int64 = 1
+	for _, f := range Factorize(b) {
+		for ; i < len(factors) && factors[i].Prime < f.Prime; i++ {
+		}
+		if i < len(factors) && factors[i].Prime == f.Prime {
+			var exponent = ints.Min(f.Exponent, factors[i].Exponent)
+			gcd *= ints.Pow64(f.Prime, int64(exponent))
+			i++
+		}
+	}
+	return gcd
+}
+
+// LeastCommonMultiple returns the least common multiple of a and b.
+func LeastCommonMultiple(a, b int64) int64 { return LCM(a, b) }
+
+// LCM returns the least common multiple of a and b.
+func LCM(a, b int64) int64 {
+	factors := Factorize(a)
+	var i int
+	var lcm int64 = 1
+	for _, f := range Factorize(b) {
+		for i < len(factors) {
+			var p = factors[i].Prime
+			if p < f.Prime {
+				lcm *= p
+				i++
+			} else {
+				break
+			}
+		}
+		var exponent = f.Exponent
+		if i < len(factors) && factors[i].Prime == f.Prime {
+			exponent = ints.Max(exponent, factors[i].Exponent)
+			i++
+		}
+		lcm *= ints.Pow64(f.Prime, int64(exponent))
+	}
+	for ; i < len(factors); i++ {
+		f := factors[i]
+		lcm *= ints.Pow64(f.Prime, int64(f.Exponent))
+	}
+	return lcm
+}
+
+// RadInt returns the radical of integer n.
+func RadInt(n int) int {
+	if 1 == n {
+		return 1
+	}
+	factors := Factorize(int64(n))
+	var r int = 1
+	for _, f := range factors {
+		r *= int(f.Prime)
+	}
+	return r
+}
+
+// RadInt64 returns the radical of integer n.
+func RadInt64(n int64) int64 {
+	if 1 == n {
+		return 1
+	}
+	factors := Factorize(n)
+	var r int64 = 1
+	for _, f := range factors {
+		r *= f.Prime
+	}
+	return r
+}

primes/ints_test.go

@@ -0,0 +1,36 @@
+package primes
+
+import (
+	"testing"
+
+	"github.com/stretchr/testify/assert"
+)
+
+func TestLCM(t *testing.T) {
+	var tests = []struct {
+		Expected int64
+		A, B     int64
+	}{
+		{2 * 3 * 5 * 7, 2 * 3 * 5, 2 * 3 * 7},
+		{2 * 2 * 3 * 5 * 7, 2 * 3 * 5, 2 * 2 * 7},
+	}
+	for _, test := range tests {
+		assert.Equal(t, test.Expected, LCM(test.A, test.B))
+		assert.Equal(t, test.Expected, LCM(test.B, test.A))
+	}
+}
+
+func TestGCD(t *testing.T) {
+	var tests = []struct {
+		Expected int64
+		A, B     int64
+	}{
+		{2 * 3, 2 * 3 * 5, 2 * 3 * 7},
+		{2, 2 * 3 * 5, 2 * 2 * 7},
+		{1, 3 * 5, 2 * 2 * 7},
+	}
+	for _, test := range tests {
+		assert.Equal(t, test.Expected, GCD(test.A, test.B))
+		assert.Equal(t, test.Expected, GCD(test.B, test.A))
+	}
+}

primes/iterator.go

@@ -0,0 +1,42 @@
+package primes
+
+// Iterator represents an iterator for iterating over numbers (prime numbers).
+type Iterator interface {
+	Next() bool
+	Get() int64
+	GoTo(p int64) bool
+	Nth() int
+}
+
+// NewIterator creates a new iterator for prime numbers (uses a default cache).
+func NewIterator() Iterator {
+	return &iterator{primes, -1}
+}
+
+// NewIteratorFromCache creates a new iterator for prime numbers and specifies the cache to use.
+func NewIteratorFromCache(c Cache) Iterator {
+	return &iterator{c, -1}
+}
+
+type iterator struct {
+	cache  Cache
+	offset int
+}
+
+func (i *iterator) Next() bool {
+	i.offset++
+	return true
+}
+
+func (i *iterator) GoTo(p int64) bool {
+	i.offset = i.cache.OffsetOfPrime(p)
+	return true
+}
+
+func (i *iterator) Get() int64 {
+	return i.cache.Prime(i.Nth())
+}
+
+func (i *iterator) Nth() int {
+	return i.offset + 1
+}

primes/math.go

@@ -0,0 +1,23 @@
+package primes
+
+// DecimalPeriod tries to find the decimal period (repetition) of the decimals of 1/n. Panics if n is not a prime number.
+func DecimalPeriod(n int) int {
+	return DecimalPeriod64(int64(n))
+}
+
+// DecimalPeriod64 tries to find the decimal period (repetition) of the decimals of 1/n. Panics if n is not a prime number.
+func DecimalPeriod64(n int64) int {
+	if !primes.IsPrime(n) {
+		panic("not a prime")
+	}
+	if n <= 5 {
+		return 1
+	}
+	length := 1
+	a := int64(1)
+	for a > 0 {
+		a = (10*a + 1) % n
+		length++
+	}
+	return length
+}

primes/math_test.go

@@ -0,0 +1,17 @@
+package primes
+
+import (
+	"testing"
+
+	"github.com/stretchr/testify/assert"
+)
+
+func TestDecimalPeriod(t *testing.T) {
+	tests := []struct {
+		Value  int
+		Period int
+	}{{2, 1}, {3, 1}, {5, 1}, {7, 6}, {11, 2}, {13, 6}, {17, 16}, {41, 5}, {271, 5}}
+	for _, test := range tests {
+		assert.Equal(t, test.Period, DecimalPeriod(test.Value), "for p = %d", test.Value)
+	}
+}

primes/primes.go

@@ -0,0 +1,16 @@
+package primes
+
+var seed = []int64{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71}
+
+// FirstN returns the first n-amount of prime numbers.
+func FirstN(n int) []int64 {
+	if n <= 20 {
+		return seed
+	}
+	return nil
+}
+
+// IsPrime checks if n is prime number.
+func IsPrime(n int64) bool {
+	return primes.IsPrime(n)
+}

primes/wheel.go

@@ -0,0 +1,58 @@
+package primes
+
+type wheel struct {
+	cache    Cache
+	elements []int64
+	primes   []int64
+	size     int64
+	offset   int
+	n        int64
+}
+
+func newWheel6(cache Cache) *wheel {
+	return &wheel{cache, []int64{2, 4}, []int64{2, 3}, 6, 0, 5}
+}
+
+func (s *wheel) resetTo(n int64) {
+	nn := n + 1
+	s.n = (nn/s.size)*s.size - 1
+	s.offset = 0
+	for s.n+s.elements[s.offset] <= n {
+		s.next()
+	}
+}
+
+func (s *wheel) next() int64 {
+	s.n += s.elements[s.offset]
+	s.offset = (s.offset + 1) % len(s.elements)
+	return s.n
+}
+
+func (s *wheel) expand() *wheel {
+	prime := s.primes[len(s.primes)-1]
+	it := NewIteratorFromCache(s.cache)
+	for it.Next() && it.Get() <= prime {
+	}
+	prime = it.Get()
+	if s.n <= prime {
+		return nil
+	}
+	primes := append([]int64(nil), s.primes...)
+	primes = append(primes, prime)
+	var elements []int64
+	var v int64 = -1
+	var inc int64
+	for i := int64(0); i < prime; i++ {
+		for _, el := range s.elements {
+			v += el
+			inc += el
+			if 0 != v%prime {
+				elements = append(elements, inc)
+				inc = 0
+			}
+		}
+	}
+	exp := &wheel{s.cache, elements, primes, s.size * prime, 0, s.n}
+	exp.resetTo(s.n)
+	return exp
+}

primes/wheel_test.go

@@ -0,0 +1,32 @@
+package primes
+
+import (
+	"testing"
+
+	"github.com/stretchr/testify/assert"
+)
+
+func TestWheelCantBeExpandedYet(t *testing.T) {
+	s := newWheel6(primes)
+	assert.Nil(t, s.expand())
+}
+
+func TestWheelIsExpanded(t *testing.T) {
+	s6 := newWheel6(primes)
+	s6.resetTo(419)
+	assert.Equal(t, 0, s6.offset)
+
+	s30 := s6.expand()
+	assert.EqualValues(t, []int64{2, 6, 4, 2, 4, 2, 4, 6}, s30.elements)
+	assert.EqualValues(t, []int64{2, 3, 5}, s30.primes)
+	assert.Equal(t, int64(30), s30.size)
+	assert.Equal(t, 0, s30.offset)
+	assert.Equal(t, s6.n, s30.n)
+
+	s210 := s30.expand()
+	assert.EqualValues(t, []int64{2, 3, 5, 7}, s210.primes)
+	assert.Equal(t, 48, len(s210.elements))
+	assert.Equal(t, int64(210), s210.size)
+	assert.Equal(t, 0, s210.offset)
+	assert.Equal(t, s30.n, s210.n)
+}