From 65e3a80e843c245f067b0355de9747ad22fc3563 Mon Sep 17 00:00:00 2001 From: Sander Schobers Date: Sat, 21 Dec 2019 13:44:19 +0100 Subject: [PATCH] Added primes package. --- primes/cache.go | 158 ++++++++++++++++++++++++++++++++++++++++++ primes/cache_test.go | 65 +++++++++++++++++ primes/factor.go | 92 ++++++++++++++++++++++++ primes/factor_test.go | 36 ++++++++++ primes/ints.go | 81 ++++++++++++++++++++++ primes/ints_test.go | 36 ++++++++++ primes/iterator.go | 42 +++++++++++ primes/math.go | 23 ++++++ primes/math_test.go | 17 +++++ primes/primes.go | 16 +++++ primes/wheel.go | 58 ++++++++++++++++ primes/wheel_test.go | 32 +++++++++ 12 files changed, 656 insertions(+) create mode 100644 primes/cache.go create mode 100644 primes/cache_test.go create mode 100644 primes/factor.go create mode 100644 primes/factor_test.go create mode 100644 primes/ints.go create mode 100644 primes/ints_test.go create mode 100644 primes/iterator.go create mode 100644 primes/math.go create mode 100644 primes/math_test.go create mode 100644 primes/primes.go create mode 100644 primes/wheel.go create mode 100644 primes/wheel_test.go diff --git a/primes/cache.go b/primes/cache.go new file mode 100644 index 0000000..994418f --- /dev/null +++ b/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) +} diff --git a/primes/cache_test.go b/primes/cache_test.go new file mode 100644 index 0000000..e0c4340 --- /dev/null +++ b/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) + } +} diff --git a/primes/factor.go b/primes/factor.go new file mode 100644 index 0000000..4e035bf --- /dev/null +++ b/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 +} diff --git a/primes/factor_test.go b/primes/factor_test.go new file mode 100644 index 0000000..0c8e59b --- /dev/null +++ b/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"))) +} diff --git a/primes/ints.go b/primes/ints.go new file mode 100644 index 0000000..6b995d1 --- /dev/null +++ b/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 +} diff --git a/primes/ints_test.go b/primes/ints_test.go new file mode 100644 index 0000000..2bee17e --- /dev/null +++ b/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)) + } +} diff --git a/primes/iterator.go b/primes/iterator.go new file mode 100644 index 0000000..8f4a8f4 --- /dev/null +++ b/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 +} diff --git a/primes/math.go b/primes/math.go new file mode 100644 index 0000000..aed72e4 --- /dev/null +++ b/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 +} diff --git a/primes/math_test.go b/primes/math_test.go new file mode 100644 index 0000000..98cf802 --- /dev/null +++ b/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) + } +} diff --git a/primes/primes.go b/primes/primes.go new file mode 100644 index 0000000..9134fe8 --- /dev/null +++ b/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) +} diff --git a/primes/wheel.go b/primes/wheel.go new file mode 100644 index 0000000..8ffda4e --- /dev/null +++ b/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 +} diff --git a/primes/wheel_test.go b/primes/wheel_test.go new file mode 100644 index 0000000..845980d --- /dev/null +++ b/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) +}