primes.go (6763B)
1 // Copyright (c) 2014 The mathutil Authors. All rights reserved. 2 // Use of this source code is governed by a BSD-style 3 // license that can be found in the LICENSE file. 4 5 package mathutil // import "modernc.org/mathutil" 6 7 import ( 8 "math" 9 ) 10 11 // IsPrimeUint16 returns true if n is prime. Typical run time is few ns. 12 func IsPrimeUint16(n uint16) bool { 13 return n > 0 && primes16[n-1] == 1 14 } 15 16 // NextPrimeUint16 returns first prime > n and true if successful or an 17 // undefined value and false if there is no next prime in the uint16 limits. 18 // Typical run time is few ns. 19 func NextPrimeUint16(n uint16) (p uint16, ok bool) { 20 return n + uint16(primes16[n]), n < 65521 21 } 22 23 // IsPrime returns true if n is prime. Typical run time is about 100 ns. 24 func IsPrime(n uint32) bool { 25 switch { 26 case n&1 == 0: 27 return n == 2 28 case n%3 == 0: 29 return n == 3 30 case n%5 == 0: 31 return n == 5 32 case n%7 == 0: 33 return n == 7 34 case n%11 == 0: 35 return n == 11 36 case n%13 == 0: 37 return n == 13 38 case n%17 == 0: 39 return n == 17 40 case n%19 == 0: 41 return n == 19 42 case n%23 == 0: 43 return n == 23 44 case n%29 == 0: 45 return n == 29 46 case n%31 == 0: 47 return n == 31 48 case n%37 == 0: 49 return n == 37 50 case n%41 == 0: 51 return n == 41 52 case n%43 == 0: 53 return n == 43 54 case n%47 == 0: 55 return n == 47 56 case n%53 == 0: 57 return n == 53 // Benchmarked optimum 58 case n < 65536: 59 // use table data 60 return IsPrimeUint16(uint16(n)) 61 default: 62 mod := ModPowUint32(2, (n+1)/2, n) 63 if mod != 2 && mod != n-2 { 64 return false 65 } 66 blk := &lohi[n>>24] 67 lo, hi := blk.lo, blk.hi 68 for lo <= hi { 69 index := (lo + hi) >> 1 70 liar := liars[index] 71 switch { 72 case n > liar: 73 lo = index + 1 74 case n < liar: 75 hi = index - 1 76 default: 77 return false 78 } 79 } 80 return true 81 } 82 } 83 84 // IsPrimeUint64 returns true if n is prime. Typical run time is few tens of µs. 85 // 86 // SPRP bases: http://miller-rabin.appspot.com 87 func IsPrimeUint64(n uint64) bool { 88 switch { 89 case n%2 == 0: 90 return n == 2 91 case n%3 == 0: 92 return n == 3 93 case n%5 == 0: 94 return n == 5 95 case n%7 == 0: 96 return n == 7 97 case n%11 == 0: 98 return n == 11 99 case n%13 == 0: 100 return n == 13 101 case n%17 == 0: 102 return n == 17 103 case n%19 == 0: 104 return n == 19 105 case n%23 == 0: 106 return n == 23 107 case n%29 == 0: 108 return n == 29 109 case n%31 == 0: 110 return n == 31 111 case n%37 == 0: 112 return n == 37 113 case n%41 == 0: 114 return n == 41 115 case n%43 == 0: 116 return n == 43 117 case n%47 == 0: 118 return n == 47 119 case n%53 == 0: 120 return n == 53 121 case n%59 == 0: 122 return n == 59 123 case n%61 == 0: 124 return n == 61 125 case n%67 == 0: 126 return n == 67 127 case n%71 == 0: 128 return n == 71 129 case n%73 == 0: 130 return n == 73 131 case n%79 == 0: 132 return n == 79 133 case n%83 == 0: 134 return n == 83 135 case n%89 == 0: 136 return n == 89 // Benchmarked optimum 137 case n <= math.MaxUint16: 138 return IsPrimeUint16(uint16(n)) 139 case n <= math.MaxUint32: 140 return ProbablyPrimeUint32(uint32(n), 11000544) && 141 ProbablyPrimeUint32(uint32(n), 31481107) 142 case n < 105936894253: 143 return ProbablyPrimeUint64_32(n, 2) && 144 ProbablyPrimeUint64_32(n, 1005905886) && 145 ProbablyPrimeUint64_32(n, 1340600841) 146 case n < 31858317218647: 147 return ProbablyPrimeUint64_32(n, 2) && 148 ProbablyPrimeUint64_32(n, 642735) && 149 ProbablyPrimeUint64_32(n, 553174392) && 150 ProbablyPrimeUint64_32(n, 3046413974) 151 case n < 3071837692357849: 152 return ProbablyPrimeUint64_32(n, 2) && 153 ProbablyPrimeUint64_32(n, 75088) && 154 ProbablyPrimeUint64_32(n, 642735) && 155 ProbablyPrimeUint64_32(n, 203659041) && 156 ProbablyPrimeUint64_32(n, 3613982119) 157 default: 158 return ProbablyPrimeUint64_32(n, 2) && 159 ProbablyPrimeUint64_32(n, 325) && 160 ProbablyPrimeUint64_32(n, 9375) && 161 ProbablyPrimeUint64_32(n, 28178) && 162 ProbablyPrimeUint64_32(n, 450775) && 163 ProbablyPrimeUint64_32(n, 9780504) && 164 ProbablyPrimeUint64_32(n, 1795265022) 165 } 166 } 167 168 // NextPrime returns first prime > n and true if successful or an undefined value and false if there 169 // is no next prime in the uint32 limits. Typical run time is about 2 µs. 170 func NextPrime(n uint32) (p uint32, ok bool) { 171 switch { 172 case n < 65521: 173 p16, _ := NextPrimeUint16(uint16(n)) 174 return uint32(p16), true 175 case n >= math.MaxUint32-4: 176 return 177 } 178 179 n++ 180 var d0, d uint32 181 switch mod := n % 6; mod { 182 case 0: 183 d0, d = 1, 4 184 case 1: 185 d = 4 186 case 2, 3, 4: 187 d0, d = 5-mod, 2 188 case 5: 189 d = 2 190 } 191 192 p = n + d0 193 if p < n { // overflow 194 return 195 } 196 197 for { 198 if IsPrime(p) { 199 return p, true 200 } 201 202 p0 := p 203 p += d 204 if p < p0 { // overflow 205 break 206 } 207 208 d ^= 6 209 } 210 return 211 } 212 213 // NextPrimeUint64 returns first prime > n and true if successful or an undefined value and false if there 214 // is no next prime in the uint64 limits. Typical run time is in hundreds of µs. 215 func NextPrimeUint64(n uint64) (p uint64, ok bool) { 216 switch { 217 case n < 65521: 218 p16, _ := NextPrimeUint16(uint16(n)) 219 return uint64(p16), true 220 case n >= 18446744073709551557: // last uint64 prime 221 return 222 } 223 224 n++ 225 var d0, d uint64 226 switch mod := n % 6; mod { 227 case 0: 228 d0, d = 1, 4 229 case 1: 230 d = 4 231 case 2, 3, 4: 232 d0, d = 5-mod, 2 233 case 5: 234 d = 2 235 } 236 237 p = n + d0 238 if p < n { // overflow 239 return 240 } 241 242 for { 243 if ok = IsPrimeUint64(p); ok { 244 break 245 } 246 247 p0 := p 248 p += d 249 if p < p0 { // overflow 250 break 251 } 252 253 d ^= 6 254 } 255 return 256 } 257 258 // FactorTerm is one term of an integer factorization. 259 type FactorTerm struct { 260 Prime uint32 // The divisor 261 Power uint32 // Term == Prime^Power 262 } 263 264 // FactorTerms represent a factorization of an integer 265 type FactorTerms []FactorTerm 266 267 // FactorInt returns prime factorization of n > 1 or nil otherwise. 268 // Resulting factors are ordered by Prime. Typical run time is few µs. 269 func FactorInt(n uint32) (f FactorTerms) { 270 switch { 271 case n < 2: 272 return 273 case IsPrime(n): 274 return []FactorTerm{{n, 1}} 275 } 276 277 f, w := make([]FactorTerm, 9), 0 278 for p := 2; p < len(primes16); p += int(primes16[p]) { 279 if uint(p*p) > uint(n) { 280 break 281 } 282 283 power := uint32(0) 284 for n%uint32(p) == 0 { 285 n /= uint32(p) 286 power++ 287 } 288 if power != 0 { 289 f[w] = FactorTerm{uint32(p), power} 290 w++ 291 } 292 if n == 1 { 293 break 294 } 295 } 296 if n != 1 { 297 f[w] = FactorTerm{n, 1} 298 w++ 299 } 300 return f[:w] 301 } 302 303 // PrimorialProductsUint32 returns a slice of numbers in [lo, hi] which are a 304 // product of max 'max' primorials. The slice is not sorted. 305 // 306 // See also: http://en.wikipedia.org/wiki/Primorial 307 func PrimorialProductsUint32(lo, hi, max uint32) (r []uint32) { 308 lo64, hi64 := int64(lo), int64(hi) 309 if max > 31 { // N/A 310 max = 31 311 } 312 313 var f func(int64, int64, uint32) 314 f = func(n, p int64, emax uint32) { 315 e := uint32(1) 316 for n <= hi64 && e <= emax { 317 n *= p 318 if n >= lo64 && n <= hi64 { 319 r = append(r, uint32(n)) 320 } 321 if n < hi64 { 322 p, _ := NextPrime(uint32(p)) 323 f(n, int64(p), e) 324 } 325 e++ 326 } 327 } 328 329 f(1, 2, max) 330 return 331 }