1 files changed, 136 insertions, 0 deletions
diff --git a/factors/type.go b/factors/type.go
new file mode 100644
index 0000000..e587ebb
--- /dev/null
+++ b/factors/type.go
@@ -0,0 +1,136 @@
+package factors
+
+import "slices"
+
+// The set of all colossaly abundant numbers that are small enough
+// to be stored in a uint32.
+// These numbers are also the first 14 superior highly composites.
+var colossallyAbundantNums = [...]uint32{
+	2, 6, 12, 60, 120, 360, 2520, 5040, 55440,
+	720720, 1441440, 4324320, 21621600, 367567200}
+
+// A cache containing all superabundant numbers up to a certain value
+var sanCache = []uint32{1}
+
+// A type of number (as determined by its factors)
+type NumberType byte
+
+const (
+	// A number whose factor score is higher than any other,
+	// if you adjust for size by dividing by some power of the number
+	// (different powers yield different best numbers).
+	// All colossally abundant numbers are also superabundant.
+	ColossallyAbundant NumberType = 0x84
+	// A number whose factor score is higher than any smaller number.
+	// All superabundant numbers have ordered exponents.
+	Superabundant NumberType = 0x83
+	// A number whose prime factorization exponents stay the same or decrease
+	// as you go from smaller to larger primes.
+	// All of these numbers are also practical.
+	OrderedExponent NumberType = 0x82
+	// A number whose factors can sum to any smaller number without duplication.
+	// All practical numbers besides 1 and 2 are divisible by 4 or 6.
+	Practical NumberType = 0x81
+	// None of the above types
+	NotPractical NumberType = 0x80
+)
+
+func Type(n uint32) NumberType {
+	if slices.Contains(colossallyAbundantNums[:], n) {
+		return ColossallyAbundant
+	} else if isSAN(n) {
+		return Superabundant
+	} else if exponentsOrdered(n) {
+		return OrderedExponent
+	} else if practical(n) {
+		return Practical
+	} else {
+		return NotPractical
+	}
+}
+
+func isSAN(n uint32) bool {
+	if n <= 2 {
+		return true
+	} else if n % 4 != 0 && n % 6 != 0 {
+		return false
+	}
+
+	cachedMax := sanCache[len(sanCache)-1]
+	maxScore := Score(uint(cachedMax))
+	nScore := Score(uint(n))
+	for i := cachedMax + 1; i <= n; i++ {
+		iScore := Score(uint(i))
+		if iScore > maxScore {
+			sanCache = append(sanCache, i)
+			maxScore = iScore
+			if maxScore > nScore {
+				return false
+			}
+		}
+	}
+
+	_, found := slices.BinarySearch(sanCache, n)
+	return found
+}
+
+func exponentsOrdered(n uint32) bool {
+	if n <= 2 {
+		return true
+	} else if n % 4 != 0 && n % 6 != 0 {
+		return false
+	}
+
+	pf := PrimeFactorize(uint(n))
+	maxPrime := slices.Max(pf.Primes())
+
+	lastPrime := uint(2)
+	for p := uint(3); p <= maxPrime; p += 2 {
+		if PrimeFactorize(p).Exponent(p) == 1 {
+			if pf.Exponent(p) > pf.Exponent(lastPrime) {
+				return false
+			} else if pf.Exponent(p) == 0 && p < maxPrime {
+				return false
+			}
+			lastPrime = p
+		}
+	}
+
+	return true
+}
+
+func practical(n uint32) bool {
+	if n <= 2 {
+		return true
+	} else if n % 4 != 0 && n % 6 != 0 {
+		return false
+	}
+
+	pf := PrimeFactorize(uint(n))
+	primes := pf.Primes()
+	slices.Sort(primes)
+
+	// algorithm from Wikipedia
+	for i := 0; i < pf.Size()-1; i++ {
+		factorSumUptoP := uint(1)
+		for j := 0; j <= i; j++ {
+			pj := primes[j]
+			ej := pf.Exponent(pj)
+			factorSumUptoP *= (uintpow(pj, ej+1) - 1) / (pj - 1)
+		}
+
+		if primes[i+1] > 1+factorSumUptoP {
+			return false
+		}
+	}
+
+	return true
+}
+
+func uintpow(a, b uint) uint {
+	result := uint(1)
+	for i := uint(0); i < b; i++ {
+		result *= a
+	}
+	return result
+}