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