Module:Rational: Difference between revisions

← Older edit Newer edit →

@@ Line 394: / Line 394: @@
 		return p.eq(a, 1) or p.eq(a, 4) or p.eq(a, 36)
 	end
 	-- now we actually check whether it is highly composite
 	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 divisors = p.divisors(a)
-	-- FIXME: inefficient; we should iterate factorisations of integers <n directly
+	local diagram_size = log2_n
-	for i = n - 1, 1, -1 do
+	local diagram = {log2_n}
-		if p.divisors(i) >= divisors then
+	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_divisors = 1
+		for i = 1, #diagram do
+			diagram_divisors = diagram_divisors * (diagram[i] + 1)
+		end
+		if diagram_divisors >= divisors then
 			return false
 		end
+		diagram = u.next_young_diagram(diagram)
 	end
 	return true