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.

3 files changed, 205 insertions, 5 deletions

diff --git a/factors/factors_test.go b/factors/factors_test.go
index db6e098..12d4384 100644
--- a/factors/factors_test.go
+++ b/factors/factors_test.go
@@ -87,7 +87,7 @@ var factorScoreCases = map[uint]float64{
 
 func TestFactorScore(t *testing.T) {
 	// factors.Score is accurate enough that we can test for exact floats!
-	tableTest(t, Score, factorScoreCases, stdEquals[float64], "Score")
+	tableTest(t, Score, factorScoreCases, stdEquals, "Score")
 }
 
 var basicRankCases = map[uint]string{
@@ -97,7 +97,7 @@ var basicRankCases = map[uint]string{
 }
 
 func TestBasicRank(t *testing.T) {
-	tableTest(t, BasicRank, basicRankCases, stdEquals[string], "BasicRank")
+	tableTest(t, BasicRank, basicRankCases, stdEquals, "BasicRank")
 }
 
 var gcdCases = map[struct{ a, b uint }]uint{
@@ -114,7 +114,7 @@ func gcdTest(c struct{ a, b uint }) uint {
 }
 
 func TestGCD(t *testing.T) {
-	tableTest(t, gcdTest, gcdCases, stdEquals[uint], "gcd")
+	tableTest(t, gcdTest, gcdCases, stdEquals, "gcd")
 }
 
 var mtcCases = map[uint]uint{
@@ -127,7 +127,60 @@ var mtcCases = map[uint]uint{
 }
 
 func TestMTC(t *testing.T) {
-	tableTest(t, MTC, mtcCases, stdEquals[uint], "MTC")
+	tableTest(t, MTC, mtcCases, stdEquals, "MTC")
+}
+
+var sanCases = map[uint32]bool{
+	1: true, 2: true, 4: true, 6: true, 12: true, 24: true, 36: true,
+	48: true, 60: true, 120: true, 180: true, 240: true, 360: true,
+	720: true, 840: true, 1260: true, 1680: true, 2520: true, 5040: true,
+	10080: true, 15120: true, 25200: true, 27720: true, 55440: true,
+	3: false, 8: false, 13: false, 18: false, 20: false, 30: false,
+	31: false, 72: false, 107: false, 135: false, 144: false, 300: false,
+	1728: false, 4000: false, 6912: false, 7560: false,
+}
+
+func TestSAN(t *testing.T) {
+	tableTest(t, isSAN, sanCases, stdEquals, "isSAN")
+}
+
+var expOrdCases = map[uint32]bool{
+	1: true, 2: true, 4: true, 5: false, 6: true, 12: true, 18: false,
+	20: false, 30: true, 44: false, 60: true, 70: false, 96: true,
+	101: false, 120: true, 144: true, 770: false, 6912: true, 10010: false,
+}
+
+func TestExponentsOrdered(t *testing.T) {
+	tableTest(t, exponentsOrdered, expOrdCases, stdEquals, "exponentsOrdered")
+}
+
+var practicalCases = map[uint32]bool{
+	1: true, 2: true, 4: true, 6: true, 8: true, 12: true, 16: true,
+	18: true, 20: true, 24: true, 28: true, 30: true, 32: true, 36: true,
+	48: true, 60: true, 120: true, 180: true, 240: true, 360: true,
+	720: true, 840: true, 1260: true, 1680: true, 2520: true, 5040: true,
+	10080: true, 15120: true, 25200: true, 27720: true, 55440: true,
+	3: false, 10: false, 27: false, 44: false, 45: false, 68: false,
+	70: false, 114: false, 121: false, 770: false, 1729: false, 10010: false,
+}
+
+func TestPractical(t *testing.T) {
+	tableTest(t, practical, practicalCases, stdEquals, "practical")
+}
+
+var typeCases = map[uint32]NumberType {
+	2: ColossallyAbundant, 3: NotPractical, 4: Superabundant,
+	6: ColossallyAbundant, 8: OrderedExponent, 10: NotPractical,
+	12: ColossallyAbundant, 18: Practical, 20: Practical, 24: Superabundant,
+	28: Practical, 30: OrderedExponent, 31: NotPractical, 34: NotPractical,
+	60: ColossallyAbundant, 70: NotPractical, 90: Practical,
+	120: ColossallyAbundant, 240: Superabundant, 720: Superabundant,
+	770: NotPractical, 1729: NotPractical, 6912: OrderedExponent,
+	10010: NotPractical, 10080: Superabundant,
+}
+
+func TestType(t *testing.T) {
+	tableTest(t, Type, typeCases, stdEquals, "Type")
 }
 
 // to be used as the equal paramater for tableTest
diff --git a/factors/prime_factorization.go b/factors/prime_factorization.go
index 69d0e08..f91a2cf 100644
--- a/factors/prime_factorization.go
+++ b/factors/prime_factorization.go
@@ -45,7 +45,7 @@ func (factors PrimeFactorization) Exponent(p uint) uint {
 	return factors.exponents[p]
 }
 
-// ExponentMap creates and returns a map mapping primes to their exponents.
+// ExponentMap returns a map mapping primes to their exponents.
 func (factors PrimeFactorization) ExponentMap() map[uint]uint {
 	exponentMap := make(map[uint]uint, len(factors.exponents))
 	for p, e := range factors.exponents {
@@ -54,6 +54,17 @@ func (factors PrimeFactorization) ExponentMap() map[uint]uint {
 	return exponentMap
 }
 
+// Primes returns a slice containing all of the primes with nonzero exponents
+func (factors PrimeFactorization) Primes() []uint {
+	primes := make([]uint, 0, len(factors.exponents))
+	for p := range factors.exponents {
+		if factors.exponents[p] > 0 {
+			primes = append(primes, p)
+		}
+	}
+	return primes
+}
+
 // Size returns how many primes in this factorization have nonzero exponents.
 func (factors PrimeFactorization) Size() int {
 	return len(factors.exponents)
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
+}