Add radix type to output

This type measures which kind of classes each radix is a part of:
- Colossally Abundant (OEIS: A004490; factor score better than every
  other number if you account for size be dividing by a certain power of
  the number)
- Superabundant (OEIS: A004394; factor score better than every smaller
  number)
- Ordered-Exponent (OEIS: A025487; exponents in prime factorization go
  down as you get to bigger primes, and no prime is skipped)
- Practical (OEIS: A005153; factors can sum to any number below the
  original number without duplication)

Each of these groups is a subset of the next, so only the most specific
label is reported.

The purpose of this program is to give you useful info to help you
determine which radices are the best, and these categories give a rough,
quantitative measure of how useful a radix's factors are:
- Practical is approximately the minimum requirement for a worthwhile
  radix.  Non-practical radices above ~16 are probably terrible to use.
- Ordered-Exponent radices act like local maxima - you can't get any
  better (smaller) without changing the "shape" (exponents) of your prime
  factorization.
- Superabundant radices are the best radices below the next
  superabundant number (e.g. 12 is the best radix below 24).
- Colossally abundant radices are, in some sense, the best radices out of
  all numbers.

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
+}