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