/* This script is part of radix_info.
   Copyright (C) 2023  Adrien Hopkins

   This program is free software: you can redistribute it and/or modify
   it under the terms of version 3 of the GNU General Public License
   as published by the Free Software Foundation.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <https://www.gnu.org/licenses/>.
*/

package factors

// The set of all colossaly abundant numbers that are small enough
// to be stored in a uint64.
// The first 15 are also the first 15 superior highly composites.
// Number source: The Online Encyclopedia of Integer Sequences
// can be found at the URL https://oeis.org/A004490
var colossallyAbundantNums = [...]uint64{
	2, 6, 12, 60, 120, 360, 2520, 5040, 55440,
	720720, 1441440, 4324320, 21621600, 367567200, 6983776800,
	160626866400, 321253732800, 9316358251200, 288807105787200,
	2021649740510400, 6064949221531200, 224403121196654400,
	9200527969062830400}

// Number source: The Online Encyclopedia of Integer Sequences
// can be found at the URL https://oeis.org/A004394
var superabundantNums = [...]uint64{
	1, 2, 4, 6, 12, 24, 36, 48, 60, 120, 180, 240, 360, 720, 840, 1260, 1680,
	2520, 5040, 10080, 15120, 25200, 27720, 55440, 110880, 166320, 277200,
	332640, 554400, 665280, 720720, 1441440, 2162160, 3603600, 4324320,
	7207200, 8648640, 10810800, 21621600, 36756720, 61261200, 73513440,
	122522400, 147026880, 183783600, 367567200, 698377680, 735134400,
	1102701600, 1163962800, 1396755360, 2327925600, 2793510720, 3491888400,
	6983776800, 13967553600, 20951330400, 27935107200, 41902660800,
	48886437600, 80313433200, 160626866400, 321253732800, 481880599200,
	642507465600, 963761198400, 1124388064800, 1927522396800, 2248776129600,
	3373164194400, 4497552259200, 4658179125600, 6746328388800, 9316358251200,
	13974537376800, 18632716502400, 27949074753600, 32607253879200,
	55898149507200, 65214507758400, 97821761637600, 130429015516800,
	144403552893600, 195643523275200, 288807105787200, 433210658680800,
	577614211574400, 866421317361600, 1010824870255200, 1732842634723200,
	2021649740510400, 3032474610765600, 4043299481020800, 6064949221531200,
	10685862914126400, 12129898443062400, 21371725828252800, 24259796886124800,
	30324746107656000, 32057588742379200, 37400520199442400, 64115177484758400,
	74801040398884800, 112201560598327200, 149602080797769600,
	224403121196654400, 448806242393308800, 897612484786617600,
	1122015605983272000, 1346418727179926400, 1533421328177138400,
	2244031211966544000, 3066842656354276800, 4600263984531415200,
	6133685312708553600, 9200527969062830400, 18401055938125660800,
}

// A classification of numbers, based on how many factors they have.
// Each type is a subset of the next
// (except [Practical] and [NoneOfTheAbove]),
// so a number is only counted as its most exclusive type.
type CompositenessType 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 CompositenessType = 0xC0
	// A number whose factor score is higher than any smaller number.
	// All superabundant numbers have ordered exponents.
	Superabundant CompositenessType = 0xA0
	// 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 CompositenessType = 0x80
	// 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 CompositenessType = 0x60
	// None of the above types.  This is the zero value of CompositenessType.
	None CompositenessType = 0x00
)

// Type determines the [CompositenessType] of a number.
func Type(n uint) CompositenessType {
	if contains(colossallyAbundantNums[:], uint64(n)) {
		return ColossallyAbundant
	} else if contains(superabundantNums[:], uint64(n)) {
		return Superabundant
	} else if exponentsOrdered(n) {
		return OrderedExponent
	} else if practical(n) {
		return Practical
	} else {
		return None
	}
}

// invariant: p must be odd and >= 3
func nextPrime(p uint) uint {
	possiblePrime := p + 2
	for !isPrime(possiblePrime) {
		possiblePrime += 2
	}
	return possiblePrime
}

func exponentsOrdered(n uint) bool {
	if n <= 2 {
		return true
	} else if n%4 != 0 && n%6 != 0 {
		return false
	}

	pf := PrimeFactorize(n)
	maxPrime := maxUints(pf.Primes())

	for prime, prevPrime := uint(3), uint(2); prime <= maxPrime; {
		if pf.Exponent(prime) > pf.Exponent(prevPrime) {
			return false
		} else if pf.Exponent(prime) == 0 && prime < maxPrime {
			return false
		}

		prime, prevPrime = nextPrime(prime), prime
	}

	return true
}

func practical(n uint) bool {
	if n <= 2 {
		return true
	} else if n%4 != 0 && n%6 != 0 {
		return false
	}

	pf := PrimeFactorize(uint(n))
	primes := sortUints(pf.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
}