source: code/trunk/vendor/modernc.org/mathutil/primes.go@ 822

Last change on this file since 822 was 822, checked in by yakumo.izuru, 22 months ago

Prefer immortal.run over runit and rc.d, use vendored modules
for convenience.

Signed-off-by: Izuru Yakumo <yakumo.izuru@…>
File size: 6.6 KB
RevLine 
[822]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 // import "modernc.org/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.
24func 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
87func 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.
170func 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.
215func 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.
259type FactorTerm struct {
260 Prime uint32 // The divisor
261 Power uint32 // Term == Prime^Power
262}
263
264// FactorTerms represent a factorization of an integer
265type 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.
269func 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
307func 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}
Note: See TracBrowser for help on using the repository browser.