Module:Rational

Module documentation[view] [edit] [history] [purge]

This module depends on the following other modules: Module: Sequence Module: Utils

Todo: add documentation

---

local p = {}

local seq = require("Module:Sequence")
local utils = require("Module:Utils")

-- enter a numerator n and denominator m
-- returns a table of prime factors
-- similar to a monzo, but the indices are prime numbers. 
function p.new(n, m)
	m = m or 1
	if n == 0 and m == 0 then
		return { nan = true }
	elseif n == 0 then
		return { zero = true, sign = utils.signum(m) }
	elseif m == 0 then
		return { inf = true, sign = utils.signum(n) }
	end
	local sign = utils.signum(n) * utils.signum(m)
	-- ensure n and m are positive
	n = n * utils.signum(n)
	m = m * utils.signum(m)
	-- factorize n and m separately
	local n_factors = utils.prime_factorization_raw(n)
	local m_factors = utils.prime_factorization_raw(m)
	local factors = n_factors
	factors.sign = sign
	-- subtract the factors of m from the factors of n
	for factor, power in pairs(m_factors) do
		factors[factor] = factors[factor] or 0
		factors[factor] = factors[factor] - power
		if factors[factor] == 0 then
			factors[factor] = nil -- clear the zeros
		end
	end
	return factors
end

-- copy a rational number
function p.copy(a)
	local b = {}
	for factor, power in pairs(a) do
		b[factor] = power
	end
	return b
end

-- create a rational number from continued fraction array
function p.from_continued_fraction(data)
	local val = p.new(1, 0)
	for i = #data, 1, -1 do
		val = p.add(data[i], p.inv(val))
	end
	return val
end

-- create a rational number from a string of whitespace-separated integers
function p.from_ket(s)
	local factor = 1
	local a = { sign = 1 }
	for i in s:gmatch("%S+") do
		local power = tonumber(i)
		if power == nil then
			return nil
		end

		-- find the next prime
		factor = factor + 1
		while not utils.is_prime(factor) do
			factor = factor + 1
		end

		if power ~= 0 then
			a[factor] = power
		end
	end
	return a
end

-- list convergents to `x` with a given stop condition
-- `stop` is either a number or a function of rational numbers
function p.convergents(x, stop)
	local convergents = {}
	local data = {}
	local i = 0
	while true do
		local n = math.floor(x)
		table.insert(data, n)
		local frac = p.from_continued_fraction(data)
		if type(stop) == "function" and stop(frac) then
			break
		elseif type(stop) == "number" and i >= stop then
			break
		end
		table.insert(convergents, frac)
		x = x - n
		if x == 0 then
			break
		end
		x = 1 / x
		i = i + 1
	end
	return convergents
end

-- determine whether a rational number is a convergent or a semiconvergent to `x`
-- TODO: document how this works
function p.converges(a, x)
	local _, m_a = p.as_pair(a)
	local convergents = p.convergents(x, function(b)
		local _, m_b = p.as_pair(b)
		return m_b >= m_a * 10000
	end)
	for _, b in ipairs(convergents) do
		if p.eq(a, b) then
			return "convergent"
		end
	end

	for i = 2, #convergents - 1 do
		local n_delta, m_delta = p.as_pair(convergents[i])
		local n_c, m_c = p.as_pair(convergents[i - 1])
		while true do
			n_c = n_c + n_delta
			m_c = m_c + m_delta
			local c = p.new(n_c, m_c)
			if p.as_table(c)[2] >= p.as_table(convergents[i + 1])[2] then
				break
			end
			if p.eq(a, c) then
				return "semiconvergent"
			end
		end
	end
	return false
end

-- attempt to identify the ratio as a simple S-expression
-- returns a table of matched expressions
function p.find_S_expression(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero or a.sign < 0 then
		return {}
	end
	if p.eq(a, 1) then
		return {}
	end
	local max_prime = p.max_prime(a)
	if seq.square_superparticulars[max_prime] == nil then
		return {}
	end
	local expressions = {}
	local superparticular_indices = {}
	local superparticular_ratios = {}
	for _, k_array in pairs(seq.square_superparticulars) do
		for _, k in ipairs(k_array) do
			if k <= 1000 then
				table.insert(superparticular_indices, k)

				local Sk_num = p.pow(p.new(k), 2)
				local Sk_den = p.mul(k - 1, k + 1)
				local Sk = p.div(Sk_num, Sk_den)
				superparticular_ratios[k] = Sk
			end
		end
	end

	-- is it Sk?
	for _, k in ipairs(superparticular_indices) do
		if p.eq(a, superparticular_ratios[k]) then
			table.insert(expressions, "S" .. k)
		end
	end

	-- is it Sk*S(k+1) or Sk/S(k+1) or Sk^2*S(k+1) or Sk*S(k+1)^2?
	for _, k in ipairs(superparticular_indices) do
		local r1 = superparticular_ratios[k]
		local r2 = superparticular_ratios[k + 1]
		if r1 and r2 then
			if p.eq(a, p.mul(r1, r2)) then
				table.insert(expressions, "S" .. k .. " × S" .. (k + 1))
			end
			if p.eq(a, p.div(r1, r2)) then
				table.insert(expressions, "S" .. k .. " / S" .. (k + 1))
			end
			if p.eq(a, p.mul(p.pow(r1, 2), r2)) then
				table.insert(expressions, "S" .. k .. "<sup>2</sup> × S" .. (k + 1))
			end
			if p.eq(a, p.mul(r1, p.pow(r2, 2))) then
				table.insert(expressions, "S" .. k .. " * S" .. (k + 1) .. "<sup>2</sup>")
			end
		end
	end

	-- is it Sk/S(k+2)?
	for _, k in ipairs(superparticular_indices) do
		local r1 = superparticular_ratios[k]
		local r2 = superparticular_ratios[k + 2]
		if r1 and r2 then
			if p.eq(a, p.div(r1, r2)) then
				table.insert(expressions, "S" .. k .. " / S" .. (k + 2))
			end
		end
	end

	-- is it S(k-1)*Sk*S(k+1)?
	for _, k in ipairs(superparticular_indices) do
		local r1 = superparticular_ratios[k - 1]
		local r2 = superparticular_ratios[k]
		local r3 = superparticular_ratios[k + 1]
		if r1 and r2 and r3 then
			if p.eq(a, p.mul(r1, p.mul(r2, r3))) then
				table.insert(expressions, "S" .. (k - 1) .. " × S" .. k .. " × S" .. (k + 1))
			end
		end
	end

	return expressions
end

-- multiply two rational numbers; integers are also allowed
function p.mul(a, b)
	if type(a) == "number" then
		a = p.new(a)
	end
	if type(b) == "number" then
		b = p.new(b)
	end
	-- special case: NaN
	if a.nan or b.nan then
		return { nan = true }
	end
	-- special case: infinities
	if (a.inf and not b.zero) or (b.inf and not a.zero) then
		return { inf = true, sign = a.sign * b.sign }
	end
	-- special case: infinity * zero
	if (a.inf and b.zero) or (b.inf and a.zero) then
		return { nan = true }
	end
	-- special case: zeros
	if a.zero or b.zero then
		return { zero = true, sign = a.sign * b.sign }
	end
	-- regular case: both not NaN, not infinities, not zeros
	local c = p.copy(a)
	for factor, power in pairs(b) do
		if type(factor) == "number" then
			c[factor] = c[factor] or 0
			c[factor] = c[factor] + power
			if c[factor] == 0 then
				c[factor] = nil
			end
		end
	end
	c.sign = a.sign * b.sign
	return c
end

-- compute 1/a for a rational number a; integers are also allowed
function p.inv(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	-- special case: NaN
	if a.nan then
		return { nan = true }
	end
	-- special case: infinity
	if a.inf then
		return { zero = true, sign = a.sign }
	end
	-- special case: zero
	if a.zero then
		return { inf = true, sign = a.sign }
	end
	-- regular case: not NaN, not infinity, not zero
	local b = {}
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			b[factor] = -power
		end
	end
	b.sign = a.sign
	return b
end

-- divide a rational number a by b; integers are also allowed
function p.div(a, b)
	return p.mul(a, p.inv(b))
end

-- compute a^b; b must be an integer
function p.pow(a, b)
	if type(a) == "number" then
		a = p.new(a)
	end
	if type(b) ~= "number" then
		return nil
	end
	if a.nan then
		return { nan = true }
	end
	if a.inf then
		if b == 0 then
			return { nan = true }
		elseif b > 0 then
			return { inf = true, sign = math.pow(a.sign, b) }
		else
			return { zero = true, sign = math.pow(a.sign, b) }
		end
	end
	if a.zero then
		if b == 0 then
			return p.new(1)
		elseif b > 0 then
			return { zero = true, sign = math.pow(a.sign, b) }
		else
			return { inf = true, sign = math.pow(a.sign, b) }
		end
	end
	local c = p.new(1)
	for _ = 1, math.abs(b) do
		if b > 0 then
			c = p.mul(c, a)
		else
			c = p.div(c, a)
		end
	end
	return c
end

-- compute a canonical representation of `a` modulo powers of `b`
-- TODO: describe the exact behavior
--       it seems bugged when the equave is a fraction
function p.modulo_mul(a, b)
	if type(a) == "number" then
		a = p.new(a)
	end
	if type(b) == "number" then
		b = p.new(b)
	end
	if a.nan or b.nan or a.inf or b.inf or a.zero or b.zero then
		return p.copy(a)
	end
	local neg_power = -math.huge
	local pos_power = math.huge
	for factor, power in pairs(b) do
		if type(factor) == "number" then
			if (power > 0 and (a[factor] or 0) >= 0) or (power < 0 and (a[factor] or 0) <= 0) then
				pos_power = math.min(pos_power, math.floor((a[factor] or 0) / power))
			else
				neg_power = math.max(neg_power, -math.ceil(math.abs(a[factor] or 0) / math.abs(power)))
			end
		end
	end
	local power = 0
	if neg_power ~= neg_power + 1 and neg_power < 0 then
		power = neg_power
	end
	if pos_power ~= pos_power + 1 and pos_power > 0 then
		power = pos_power
	end
	return p.div(a, p.pow(b, power))
end

-- add two rational numbers; integers are also allowed
function p.add(a, b)
	if type(a) == "number" then
		a = p.new(a)
	end
	if type(b) == "number" then
		b = p.new(b)
	end

	-- special case: NaN
	if a.nan or b.nan then
		return { nan = true }
	end
	-- special case: infinities
	if a.inf and b.inf then
		if a.sign == b.sign then
			return { inf = true, sign = a.sign }
		else
			return { nan = true }
		end
	end
	if a.inf then
		return { inf = true, sign = a.sign }
	end
	if b.inf then
		return { inf = true, sign = b.sign }
	end
	-- special case: one is zero
	if a.zero then
		return p.copy(b)
	end
	if b.zero then
		return p.copy(a)
	end
	-- regular case: both not NaN, not infinities, not zeros
	local gcd = { sign = 1 }
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			if math.min(power, b[factor] or 0) > 0 then
				gcd[factor] = math.min(power, b[factor])
			end
			if math.max(power, b[factor] or 0) < 0 then
				gcd[factor] = math.max(power, b[factor])
			end
		end
	end
	a = p.div(a, gcd)
	b = p.div(b, gcd)

	local n_a, m_a = p.as_pair(a)
	local n_b, m_b = p.as_pair(b)

	local n = n_a * m_b + n_b * m_a
	local m = m_a * m_b

	return p.mul(p.new(n, m), gcd)
end

-- substract a rational number from another; integers are also allowed
function p.sub(a, b)
	return p.add(a, p.mul(b, -1))
end

-- absolute value of a rational number; integers are also allowed
function p.abs(a)
	if a.nan then
		return { nan = true }
	end
	local b = p.copy(a)
	b.sign = 1
	return b
end

-- determine whether a rational number is less than another; integers are also allowed
function p.lt(a, b)
	local c = p.sub(a, b)
	if c.zero then
		return false
	else
		return c.sign == -1
	end
end

-- determine whether a rational number is less or equal to the other; integers are also allowed
function p.leq(a, b)
	local c = p.sub(a, b)
	if c.zero then
		return true
	else
		return c.sign == -1
	end
end

-- determine whether a rational number is greater than another; integers are also allowed
function p.gt(a, b)
	local c = p.sub(a, b)
	if c.zero then
		return false
	else
		return c.sign == 1
	end
end

-- determine whether a rational number is greater or equal to the other; integers are also allowed
function p.geq(a, b)
	local c = p.sub(a, b)
	if c.zero then
		return true
	else
		return c.sign == 1
	end
end

-- determine whether a rational number is equal to another; integers are also allowed
function p.eq(a, b)
	if type(a) == "number" then
		a = p.new(a)
	end
	if type(b) == "number" then
		b = p.new(b)
	end
	if a.nan or b.nan then
		return false
	end
	if a.inf and b.inf then
		return a.sign == b.sign
	end
	if a.inf or b.inf then
		return false
	end
	if a.zero and b.zero then
		return true
	end
	if a.zero or b.zero then
		return false
	end
	for factor, power in pairs(a) do
		if b[factor] ~= power then
			return false
		end
	end
	for factor, power in pairs(b) do
		if a[factor] ~= power then
			return false
		end
	end
	return true
end

-- determine whether a rational number is integer
function p.is_int(a)
	if type(a) == "number" then
		return true
	end
	if a.nan then
		return false
	end
	if a.inf then
		return false
	end
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			if power < 0 then
				return false
			end
		end
	end
	return true
end

-- determine whether a rational number lies within [1; equave)
function p.is_reduced(a, equave, large)
	equave = equave or 2
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero or a.sign < 0 then
		return false
	end
	if large then
		-- an approximation
		local cents = p.cents(a)
		local cents_max = p.cents(equave)
		return cents >= 0 and cents < cents_max
	else
		return p.geq(a, 1) and p.lt(a, equave)
	end
end

-- determine whether a rational number represents a harmonic. 
-- reduced: check for reduced harmonic instead. 
function p.is_harmonic(a, reduced, large)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero or a.sign < 0 then
		return false
	end
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			if factor == 2 and reduced then
				-- pass (ignore factors of 2 for reduced harmonic check)
			elseif power < 0 then
				return false
			end
		end
	end
	if reduced then
		return p.is_reduced(a, 2, large)
	end
	return true
end

-- determine whether a rational number represents a subharmonic. 
-- reduced: check for reduced subharmonic instead. 
function p.is_subharmonic(a, reduced, large)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero or a.sign < 0 then
		return false
	end
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			if factor == 2 and reduced then
				-- pass (ignore factors of 2 for reduced subharmonic check)
			elseif power > 0 then
				return false
			end
		end
	end
	if reduced then
		return p.is_reduced(a, 2, large)
	end
	return true
end

-- determine whether a rational number is an integer power of another rational number
function p.is_power(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return false
	end
	if p.eq(a, 1) or p.eq(a, -1) then
		return false
	end

	local total_power = nil
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			if total_power then
				total_power = utils._gcd(total_power, math.abs(power))
			else
				total_power = math.abs(power)
			end
		end
	end
	return total_power > 1
end

-- determine whether a rational number is superparticular
function p.is_superparticular(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return false
	end
	local n, m = p.as_pair(a)
	return n - m == 1
end

-- determine whether a rational number is a square superparticular
function p.is_square_superparticular(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero or a.sign < 0 then
		return false
	end
	-- check the numerator
	local k = { sign = 1 }
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			if power > 0 and power % 2 ~= 0 then
				return false
			elseif power > 0 then
				k[factor] = math.floor(power / 2 + 0.5)
			end
		end
	end
	-- check the denominator
	local den = p.mul(p.add(k, 1), p.sub(k, 1))
	return p.eq(a, p.div(p.pow(k, 2), den))
end

-- check if an integer is prime
function p.is_prime(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	
	-- nan, inf, zero, and negative numbers are not prime
	if a.nan or a.inf or a.zero or a.sign < 0 then
		return false
	end
	
	local flag = false -- flag for having exactly one prime factor
	for factor, power in pairs(a) do
		if type(factor) == "number" and power then
			if flag or power ~= 1 then
				return false
			else
				flag = true
			end
		end
	end
	return flag
end

-- check if an integer is highly composite
function p.is_highly_composite(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	
	-- nan, inf, zero, and negative numbers are not highly composite
	if a.nan or a.inf or a.zero or a.sign == -1 then
		return false
	end
	
	-- non-integers are not highly composite
	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)
	for i = 2, max_prime do
		if utils.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]
		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 highly composite
	local n, _ = 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)
	local diagram_size = log2_n
	local diagram = { log2_n }
	local primes = { 2 }

	local function eval_diagram(d)
		while #d > #primes do
			local i = primes[#primes] + 1
			while not utils.is_prime(i) do
				i = i + 1
			end
			table.insert(primes, i)
		end
		local m = 1
		for i = 1, #d do
			for _ = 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 = utils.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 utils.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, _ = 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 }

	local function eval_diagram(d)
		while #d > #primes do
			local i = primes[#primes] + 1
			while not utils.is_prime(i) do
				i = i + 1
			end
			table.insert(primes, i)
		end
		local m = 1
		for i = 1, #d do
			for _ = 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
		diagram = utils.next_young_diagram(diagram)
	end
	return true
end

-- Check if ratio is within an int limit; that is, neither its numerator nor
-- denominator exceed that limit.
function p.is_within_int_limit(a, lim)
	return p.int_limit(a) <= lim
end

-- Find integer limit of a ratio
-- For a ratio p/q, this is simply max(p, q)
function p.int_limit(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return nil
	end
	local a_copy = p.copy(a)
	local num, den = p.as_pair(a_copy)
	return math.max(num, den)
end

-- Find odd limit of a ratio
-- For a ratio p/q, this is simply max(p, q) where powers of 2 are ignored
function p.odd_limit(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return nil
	end
	local a_copy = p.copy(a)
	if a_copy[2] ~= nil then
		a_copy[2] = 0
	end
	local num, den = p.as_pair(a_copy)
	return math.max(num, den)
end

-- find max prime involved in the factorisation
-- (a.k.a. prime limit or harmonic class) of a rational number
function p.max_prime(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return nil
	end
	local max_factor = 0
	for factor, _ in pairs(a) do
		if type(factor) == "number" then
			if factor > max_factor then
				max_factor = factor
			end
		end
	end
	return max_factor
end

-- convert a rational number to its size in octaves
-- equal to log2 of the rational number
function p.log(a, base)
	base = base or 2
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.inf and a.sign > 0 then
		return 1 / 0
	end
	if a.nan or a.inf then
		return nil
	end
	if a.zero then
		return -1 / 0
	end
	if a.sign < 0 then
		return nil
	end
	local logarithm = 0
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			logarithm = logarithm + power * utils._log(factor, base)
		end
	end
	return logarithm
end

-- convert a rational number to its size in cents
function p.cents(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.sign < 0 then
		return nil
	end
	if a.inf and a.sign > 0 then
		return 1 / 0
	end
	if a.zero then
		return -1 / 0
	end

	local c = 0
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			c = c + power * utils.log2(factor)
		end
	end
	return c * 1200
end

-- convert a rational number (interpreted as an interval) into Hz
function p.hz(a, base)
	base = base or 440
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.sign < 0 then
		return nil
	end
	if a.zero then
		return 0
	end
	local log_hz = math.log(base)
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			log_hz = log_hz + power * math.log(factor)
		end
	end
	return math.exp(log_hz)
end

-- FJS: x = a * 2^n : x >= 1, x < 2
local function red(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return nil
	end
	local b = p.copy(a)

	-- start with an approximation
	local log2 = p.log(b)
	b = p.div(b, p.pow(2, math.floor(log2)))

	while p.lt(b, 1) do
		b = p.mul(b, 2)
	end
	while p.geq(b, 2) do
		b = p.div(b, 2)
	end
	return b
end

-- FJS: x = a * 2^n : x >= 1/sqrt(2), x < sqrt(2)
local function reb(a)
	local b = red(a)
	if p.geq(p.mul(b, b), 2) then
		b = p.div(b, 2)
	end
	return b
end

-- FJS: master algorithm
local function FJS_master(prime)
	prime = red(prime)
	local tolerance = p.new(65, 63)
	local fifth = p.new(3, 2)
	local k = 0
	while true do
		local a = red(p.pow(fifth, k))
		if math.abs(p.cents(p.div(prime, a))) < p.cents(tolerance) then
			return k
		end
		if k == 0 then
			k = 1
		elseif k > 0 then
			k = -k
		else
			k = -k + 1
		end
	end
end

-- FJS: formal comma
local function formal_comma(prime)
	local fifth_shift = FJS_master(prime)
	return reb(p.div(prime, p.pow(3, fifth_shift)))
end

-- FJS representation of a rational number
-- might be a bit incorrect
-- TODO: confirm correctness
function p.as_FJS(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return nil
	end
	local b = p.copy(a)
	local otonal = {}
	local utonal = {}
	for factor, power in pairs(a) do
		if type(factor) == "number" and factor > 3 then
			local comma = formal_comma(factor)
			b = p.div(b, p.pow(comma, power))
			if power > 0 then
				for _ = 1, power do
					table.insert(otonal, factor)
				end
			else
				for _ = 1, -power do
					table.insert(utonal, factor)
				end
			end
		end
	end
	table.sort(otonal)
	table.sort(utonal)

	local fifths = b[3] or 0

	local o = math.floor((fifths * 2 + 3) / 7)
	local num = fifths * 11 + (b[2] or 0) * 7
	if num >= 0 then
		num = num + 1
	else
		num = num - 1
		o = -o
	end

	local num_mod = (num - utils.signum(num)) % 7
	local letter = ""
	if (num_mod == 0 or num_mod == 3 or num_mod == 4) and o == 0 then
		letter = "P"
	elseif o == 1 then
		letter = "M"
	elseif o == -1 then
		letter = "m"
	else
		if o >= 0 then
			o = o - 1
		else
			o = o + 1
		end
		if o > 0 then
			while o > 0 do
				letter = letter .. "A"
				o = o - 2
			end
		else
			while o < 0 do
				letter = letter .. "d"
				o = o + 2
			end
		end
		if #letter >= 5 then
			letter = #letter .. letter:sub(1, 1)
		end
	end

	local FJS = letter .. num
	if #otonal > 0 then
		FJS = FJS .. "^{" .. table.concat(otonal, ",") .. "}"
	end
	if #utonal > 0 then
		FJS = FJS .. "_{" .. table.concat(utonal, ",") .. "}"
	end
	return FJS
end

-- determine log2 product complexity
function p.tenney_height(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return nil
	end
	local h = 0
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			h = h + math.abs(power) * utils.log2(factor)
		end
	end
	return h
end

-- determine log2 max complexity
function p.weil_height(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return nil
	end
	local h1 = p.tenney_height(a)
	local h2 = 0
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			h2 = h2 + power * utils.log2(factor)
		end
	end
	h2 = math.abs(h2)
	return h1 + h2
end

-- determine sopfr complexity
function p.wilson_height(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return nil
	end
	local h = 0
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			h = h + math.abs(power) * factor
		end
	end
	return h
end

-- determine product complexity
function p.benedetti_height(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return nil
	end
	local n, m = p.as_pair(a)
	if (math.log(n) + math.log(m)) / math.log(10) <= 15 then
		return n * m
	else
		-- it is going to be an overflow
		return nil
	end
end

-- determine the number of rational divisors
function p.divisors(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return 0
	end
	local d = 1
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			d = d * (math.abs(power) + 1)
		end
	end
	return d
end

-- determine whether the rational number is +- p/q, where p, q are primes OR 1
function p.is_prime_ratio(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero then
		return false
	end
	local n_factors = 0
	local m_factors = 0
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			if power > 0 then
				n_factors = n_factors + 1
			else
				m_factors = m_factors + 1
			end
		end
	end
	return n_factors <= 1 and m_factors <= 1
end

-- return prime factorisation of a rational number
function p.factorisation(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero or p.eq(a, 1) or p.eq(a, -1) then
		return "n/a"
	end
	local s = ""
	if a.sign < 0 then
		s = s .. "-"
	end
	local factors = {}
	for factor, _ in pairs(a) do
		if type(factor) == "number" then
			table.insert(factors, factor)
		end
	end
	table.sort(factors)
	for i, factor in ipairs(factors) do
		if i > 1 then
			s = s .. " × "
		end
		s = s .. factor
		if a[factor] ~= 1 then
			s = s .. "<sup>" .. a[factor] .. "</sup>"
		end
	end
	return s
end

-- return the subgroup generated by primes from a rational number's prime factorisation
function p.subgroup(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	if p.eq(a, 1) then
		return "1"
	end
	if a.nan or a.inf or a.zero or p.eq(a, -1) then
		return "n/a"
	end
	local s = ""
	local factors = {}
	for factor, _ in pairs(a) do
		if type(factor) == "number" then
			table.insert(factors, factor)
		end
	end
	table.sort(factors)
	for i, factor in ipairs(factors) do
		if i > 1 then
			s = s .. "."
		end
		s = s .. factor
	end
	if a.sign < 0 then
		s = "-1." .. s
	end
	return s
end

-- unpack rational as two return values (n, m)
function p.as_pair(a)
	if type(a) == "number" then
		a = p.new(a)
	end
	-- special case: NaN
	if a.nan then
		return 0, 0
	end
	-- special case: infinity
	if a.inf then
		return a.sign, 0
	end
	-- special case: zero
	if a.zero then
		return 0, a.sign
	end
	-- regular case: not NaN, not infinity, not zero
	local n = 1
	local m = 1
	for factor, power in pairs(a) do
		if type(factor) == "number" then
			if power > 0 then
				n = n * (factor ^ power)
			else
				m = m * (factor ^ -power)
			end
		end
	end
	n = n * a.sign
	return n, m
end

-- return a string ratio representation
function p.as_ratio(a, separator)
	separator = separator or "/"
	local n, m = p.as_pair(a)
	return ("%d%s%d"):format(n, separator, m)
end

-- return the {n, m} pair as a Lua table
function p.as_table(a)
	return { p.as_pair(a) }
end

-- return n / m as a float approximation
function p.as_float(a)
	local n, m = p.as_pair(a)
	return n / m
end

-- return a rational number in subgroup ket notation
function p.as_subgroup_ket(a, frame)
	if type(a) == "number" then
		a = p.new(a)
	end
	if a.nan or a.inf or a.zero or a.sign < 0 then
		return "n/a"
	end
	local factors = {}
	for factor, _ in pairs(a) do
		if type(factor) == "number" then
			table.insert(factors, factor)
		end
	end
	table.sort(factors)
	local subgroup = "1"
	if not p.eq(a, 1) then
		subgroup = table.concat(factors, ".")
	end

	local powers = {}
	for _, factor in ipairs(factors) do
		table.insert(powers, a[factor])
	end
	local template_arg = "0"
	if not p.eq(a, 1) then
		template_arg = table.concat(powers, " ")
	end

	return subgroup .. " " .. frame:expandTemplate({
		title = "Monzo",
		args = { template_arg },
	})
end

-- return a string of a rational number in monzo notation
-- calling Template: Monzo
function p.as_ket(a, frame, skip_many_zeros, only_numbers)
	if skip_many_zeros == nil then
		skip_many_zeros = true
	end
	only_numbers = only_numbers or false
	if type(a) == "number" then
		a = p.new(a)
	end
	
	-- special cases
	if a.nan or a.inf or a.zero or a.sign < 0 then
		return "n/a"
	end
	
	-- regular case: positive finite ratio
	local s = ""

	-- preparing the argument
	local max_prime = p.max_prime(a)
	local template_arg = ""
	local template_arg_without_trailing_zeros = ""
	local zeros_n = 0
	for i = 2, max_prime do
		if utils.is_prime(i) then
			if i > 2 then
				template_arg = template_arg .. " "
			end
			template_arg = template_arg .. (a[i] or 0)

			if (a[i] or 0) ~= 0 then
				if skip_many_zeros and zeros_n >= 4 then
					template_arg = template_arg_without_trailing_zeros
					if #template_arg > 0 then
						template_arg = template_arg .. " "
					end
					template_arg = template_arg .. "0<sup>" .. zeros_n .. "</sup> " .. (a[i] or 0)
				end
				zeros_n = 0
				template_arg_without_trailing_zeros = template_arg
			else
				zeros_n = zeros_n + 1
			end
		end
	end
	if #template_arg == 0 then
		template_arg = "0"
	end
	if only_numbers then
		s = s .. template_arg
	else
		s = s .. frame:expandTemplate({
			title = "Monzo",
			args = { template_arg },
		})
	end
	return s
end

-- parse a rational number
-- returns nil on failure
function p.parse(unparsed)
	if type(unparsed) ~= "string" then
		return nil
	end
	-- removing whitespaces
	unparsed = unparsed:gsub("%s", "")
	-- removing <br> and <br/> tags
	unparsed = unparsed:gsub("<br/?>", "")

	-- length limit: very long strings are not converted into Lua numbers correctly
	local max_length = 15

	-- rational form
	local sign, n, _, m = unparsed:match("^%s*(%-?)%s*(%d+)%s*(/%s*(%d+))%s*$")
	if n == nil then
		-- integer form
		sign, n = unparsed:match("^%s*(%-?)%s*(%d+)%s*$")
		if n == nil then
			-- parsing failure
			return nil
		else
			m = 1
			if #n > max_length then
				return nil
			end
			n = tonumber(n)
			if #sign > 0 then
				n = -n
			end
		end
	else
		if #n > max_length then
			return nil
		end
		n = tonumber(n)
		if #m > max_length then
			return nil
		end
		m = tonumber(m)
		if #sign > 0 then
			n = -n
		end
	end
	return p.new(n, m)
end

-- a version of as_ket() that can be {{#invoke:}}d
function p.ket(frame)
	local unparsed = frame.args[1] or "1"
	local result = ""
	
	local a = p.parse(unparsed)
	if a == nil then
		result = '{{error|Invalid rational number: ' .. unparsed .. ".}}"
	else
		result = p.as_ket(a, frame)
	end
	
	return frame:preprocess(result)
end
p.monzo = p.ket

return p