|
| 1 | +// The MIT License (MIT) |
| 2 | +// |
| 3 | +// Copyright (c) 2015 Filippo Tampieri |
| 4 | +// |
| 5 | +// Permission is hereby granted, free of charge, to any person obtaining a copy |
| 6 | +// of this software and associated documentation files (the "Software"), to deal |
| 7 | +// in the Software without restriction, including without limitation the rights |
| 8 | +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell |
| 9 | +// copies of the Software, and to permit persons to whom the Software is |
| 10 | +// furnished to do so, subject to the following conditions: |
| 11 | +// |
| 12 | +// The above copyright notice and this permission notice shall be included in |
| 13 | +// all copies or substantial portions of the Software. |
| 14 | +// |
| 15 | +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |
| 16 | +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
| 17 | +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE |
| 18 | +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
| 19 | +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, |
| 20 | +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN |
| 21 | +// THE SOFTWARE. |
| 22 | + |
| 23 | +// Package primes provides simple functionality for working with prime numbers. |
| 24 | +// |
| 25 | +// Call Sieve(n) to generate all prime numbers less than or equal to n, |
| 26 | +// IsPrime(n) to test for primality, Coprime(a,b) to test for coprimality, |
| 27 | +// and Pi(n) to count (or estimate) the number of primes less than or equal to n. |
| 28 | +// |
| 29 | +// The algorithms used to implement the functions above are fairly simple; |
| 30 | +// they work well with relatively small primes, but they are definitely not |
| 31 | +// intended for work in cryptography or any application requiring really |
| 32 | +// large primes. |
| 33 | +// |
| 34 | +package primes |
| 35 | + |
| 36 | +import ( |
| 37 | + "fmt" |
| 38 | + "math" |
| 39 | + "sort" |
| 40 | +) |
| 41 | + |
| 42 | +// primes is a cache of the first few prime numbers |
| 43 | +var primes []int |
| 44 | + |
| 45 | +func init() { |
| 46 | + // Cache the first 1,229 prime numbers (i.e. all primes <= 10,000) |
| 47 | + // FIXME: should I change n to get exactly 1,024 primes? |
| 48 | + primes = Sieve(10000) |
| 49 | + // FIXME: remove print |
| 50 | + fmt.Printf("cache len=%d, cap=%d\n", len(primes), cap(primes)) |
| 51 | +} |
| 52 | + |
| 53 | +// Pi returns the number of primes less than or equal to n. |
| 54 | +// If ok is true, the result is correct; otherwise it is an estimate |
| 55 | +// computed as n/(log(n)-1) with an error proven smaller than 11%. |
| 56 | +// This estimate is based on the prime number theorem which states that the |
| 57 | +// number of primes not exceeding n is asymptotic to n/log(n). |
| 58 | +// Better approximations are known, but they are more complex. |
| 59 | +// See https://primes.utm.edu/howmany.html#1, |
| 60 | +// https://en.wikipedia.org/wiki/Prime_number_theorem, and |
| 61 | +// https://en.wikipedia.org/wiki/Prime-counting_function for details. |
| 62 | +func Pi(n int) (pi int, ok bool) { |
| 63 | + // If n is smaller than or equal to the largest cached prime, |
| 64 | + // we have an exact count |
| 65 | + if i := sort.SearchInts(primes, n); i < len(primes) { |
| 66 | + if n == primes[i] { |
| 67 | + // n is the prime at index i |
| 68 | + pi = i + 1 |
| 69 | + } else { |
| 70 | + // n is not prime and primes[j] < n for all j in [0,i) |
| 71 | + pi = i |
| 72 | + } |
| 73 | + ok = true |
| 74 | + } else { |
| 75 | + // n is larger than the largest prime in the cache; |
| 76 | + // use the estimate |
| 77 | + pi = int(float64(n) / (math.Log(float64(n)) - 1)) |
| 78 | + } |
| 79 | + return |
| 80 | +} |
| 81 | + |
| 82 | +// IsPrime is a primality test: it returns true if n is prime. |
| 83 | +// It uses trial division to check whether n can be divided exactly by any |
| 84 | +// number that is less than n. |
| 85 | +// The algorithm can be very slow on large n despite a number of optimizations: |
| 86 | +// - If n is relatively small, it is compared against a cached table of primes. |
| 87 | +// - Only numbers up to sqrt(n) are used to check for primality. |
| 88 | +// - n is first checked for divisibility by the primes in the cache and |
| 89 | +// only if the test is inconclusive, n is checked against more numbers. |
| 90 | +// - Only odd numbers greater than the last prime in the cache are checked |
| 91 | +// after that. |
| 92 | +// See https://en.wikipedia.org/wiki/Primality_test and |
| 93 | +// https://en.wikipedia.org/wiki/Trial_division for details. |
| 94 | +func IsPrime(n int) bool { |
| 95 | + pMax := primes[len(primes)-1] |
| 96 | + if n <= pMax { |
| 97 | + // If n is prime, it must be in the cache |
| 98 | + i := sort.SearchInts(primes, n) |
| 99 | + return n == primes[i] |
| 100 | + } |
| 101 | + max := int(math.Ceil(math.Sqrt(float64(n)))) |
| 102 | + // Check if n is divisible by any of the cached primes |
| 103 | + for _, p := range primes { |
| 104 | + if p > max { |
| 105 | + return true |
| 106 | + } |
| 107 | + if n%p == 0 { |
| 108 | + return false |
| 109 | + } |
| 110 | + } |
| 111 | + // When you run out of cached primes, check if n is divisible by |
| 112 | + // any odd number larger than the largest prime in the cache. |
| 113 | + for d := pMax + 2; d <= max; d += 2 { |
| 114 | + if n%d == 0 { |
| 115 | + return false |
| 116 | + } |
| 117 | + } |
| 118 | + return true |
| 119 | +} |
| 120 | + |
| 121 | +// Coprime is a coprimality test: it returns true if the only positive integer |
| 122 | +// that divides evenly both a and b is 1. |
| 123 | +// This function implements the division-based version of the Euclidean algorithm. |
| 124 | +// See https://en.wikipedia.org/wiki/Coprime_integers and |
| 125 | +// https://en.wikipedia.org/wiki/Euclidean_algorithm for details. |
| 126 | +func Coprime(a, b int) bool { |
| 127 | + // Set a to gcd(a,b) |
| 128 | + var t int |
| 129 | + for b != 0 { |
| 130 | + t = b |
| 131 | + b = a % b |
| 132 | + a = t |
| 133 | + } |
| 134 | + // By definition, a and b are coprime if gcd(a,b) == 1 |
| 135 | + return a == 1 |
| 136 | +} |
| 137 | + |
| 138 | +// Sieve returns a list of the prime numbers less than or equal to n. |
| 139 | +// If n is less than 2, it returns an empty list. |
| 140 | +// The function uses the sieve of Eratosthenes algorithm |
| 141 | +// https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes |
| 142 | +// with the following optimizations: |
| 143 | +// - The initial list of candidate primes includes odd numbers only. |
| 144 | +// - Given a prime p, only multiples of p greater than or equal to p*p need |
| 145 | +// to be marked off since smaller multiples of p have already been marked off |
| 146 | +// by then. |
| 147 | +// - The above also implies that the algorithm can terminate as soon as it finds |
| 148 | +// a prime p such that p*p is greater than n. |
| 149 | +// Sieve takes O(n) memory and runs in O(n log log n) time. |
| 150 | +func Sieve(n int) []int { |
| 151 | + if n < 2 { |
| 152 | + return []int{} |
| 153 | + } |
| 154 | + // a[i] == false ==> p=2*i+3 is a candidate prime |
| 155 | + // p in [3,n] ==> i in [0,(n-3)/2] |
| 156 | + length := 1 + (n-3)/2 |
| 157 | + a := make([]bool, length, length) |
| 158 | + // Start with number 3 and consider only odd numbers |
| 159 | + for i, p := 0, 3; p*p <= n; p += 2 { |
| 160 | + if !a[i] { |
| 161 | + // 2*i+1 is a prime number; mark off its multiples |
| 162 | + for j := (p*p - 3) / 2; j < length; j += p { |
| 163 | + a[j] = true |
| 164 | + } |
| 165 | + } |
| 166 | + i++ |
| 167 | + } |
| 168 | + // ps will store the computed primes; its initial capacity is based |
| 169 | + // an estimate of the prime-counting function pi(n) |
| 170 | + pi, _ := Pi(n) |
| 171 | + ps := make([]int, 1, pi) |
| 172 | + ps[0] = 2 |
| 173 | + for i := 0; i < length; i++ { |
| 174 | + if !a[i] { |
| 175 | + ps = append(ps, 2*i+3) |
| 176 | + } |
| 177 | + } |
| 178 | + return ps |
| 179 | +} |
0 commit comments