Module:Rational: Difference between revisions

← Older edit Newer edit →

@@ Line 450: / Line 450: @@
 		end
 		if diagram_divisors >= divisors then
+			return false
+		end
+		diagram = u.next_young_diagram(diagram)
+	end
+	return true
+end
+-- check if an integer is superabundant
+function p.is_superabundant(a)
+	if type(a) == 'number' then
+		a = p.new(a)
+	end
+	if a.nan or a.inf or a.zero then
+		return false
+	end
+	-- negative numbers are not superabundant
+	if a.sign == -1 then
+		return false
+	end
+	-- non-integers are not superabundant
+	for factor, power in pairs(a) do
+		if type(factor) == 'number' then
+			if power < 0 then
+				return false
+			end
+		end
+	end
+	local last_power = 1/0
+	local max_prime = p.max_prime(a)
+	local divisor_sum = p.new(1)
+	for i = 2, max_prime do
+		if u.is_prime(i) then
+			-- factors must be the first N primes
+			if a[i] == nil then
+				return false
+			end
+			-- powers must form a non-increasing sequence
+			if a[i] > last_power then
+				return false
+			end
+			last_power = a[i]
+			divisor_sum = p.mul(
+				divisor_sum,
+				p.div(
+					p.sub(p.pow(i, a[i] + 1), 1),
+					i - 1
+				)
+			)
+		end
+	end
+	-- last_power may be >1 only for 1, 4, 36
+	if last_power > 1 then
+		return p.eq(a, 1) or p.eq(a, 4) or p.eq(a, 36)
+	end
+	-- now we actually check whether it is superabundant
+	local n, _m = p.as_pair(a)
+	-- precision is very important here
+	local log2_n = 0
+	local t = 1
+	while t * 2 <= n do
+		log2_n = log2_n + 1
+		t = t * 2
+	end
+	local SA_ratio = p.div(divisor_sum, a)
+	local diagram_size = log2_n
+	local diagram = {log2_n}
+	local primes = {2}
+	function eval_diagram(d)
+		while #d > #primes do
+			local i = primes[#primes] + 1
+			while not u.is_prime(i) do
+				i = i + 1
+			end
+			table.insert(primes, i)
+		end
+		local m = 1
+		for i = 1, #d do
+			for j = 1, d[i] do
+				m = m * primes[i]
+			end
+		end
+		return m
+	end
+	-- iterate factorisations of some composite integers <n
+	while diagram do
+		while eval_diagram(diagram) >= n do
+			-- reduce diagram size, preserve diagram width
+			if diagram_size <= #diagram then
+				diagram = nil
+				break
+			end
+			diagram_size = diagram_size - 1
+			diagram[1] = diagram_size - #diagram + 1
+			for i = 2, #diagram do
+				diagram[i] = 1
+			end
+		end
+		if diagram == nil then
+			break
+		end
+		local diagram_divisor_sum = 1
+		for i = 1, #diagram do
+			diagram_divisor_sum = p.mul(
+				diagram_divisor_sum,
+				p.div(
+					p.sub(p.pow(primes[i], diagram[i] + 1), 1),
+					primes[i] - 1
+				)
+			)
+		end
+		local diagram_SA_ratio = p.div(diagram_divisor_sum, eval_diagram(diagram))
+		if p.geq(diagram_SA_ratio, SA_ratio) then
 			return false
 		end