local u = require('Module:Utils')
local p = {}
-- construct a rational number n/m
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 = u.signum(m) }
elseif m == 0 then
return { inf = true, sign = u.signum(n) }
end
local sign = u.signum(n) * u.signum(m)
n = n * u.signum(n)
m = m * u.signum(m)
local n_factors = u.prime_factorization_raw(n)
local m_factors = u.prime_factorization_raw(m)
local factors = n_factors
factors.sign = sign
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
end
end
return factors
end
-- copy a rational number
function p.copy(a)
b = {}
for factor, power in pairs(a) do
b[factor] = power
end
return b
end
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
-- 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
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 i = 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`
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 = -1/0
local pos_power = 1/0
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
local a = p.div(a, gcd)
local b = p.div(b, gcd)
n_a, m_a = p.as_pair(a)
n_b, m_b = p.as_pair(b)
n = n_a * m_b + n_b * m_a
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
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
-- check if an integer is highly composite
function p.is_highly_composite(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 highly composite
if 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 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]
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, _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)
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_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
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
diagram = u.next_young_diagram(diagram)
end
return true
end
-- find max prime in ket notation
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, power in pairs(a) do
if type(factor) == 'number' then
if factor > max_factor then
max_factor = factor
end
end
end
return max_factor
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 the (n, m) pair as a Lua tuple
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 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 ket notation
-- NaN, infinity, zero values use special representations
function p.as_ket(a, frame)
if type(a) == 'number' then
a = p.new(a)
end
-- special case: NaN
if a.nan then
return 'NaN'
end
-- special case: infinity
if a.inf then
local sign = '+'
if a.sign < 0 then
sign = '-'
end
return sign .. '∞'
end
-- special case: zero
if a.zero then
return '0'
end
-- regular case: not NaN, not infinity, not zero
local s = ''
if a.sign < 0 then
s = s .. '-'
end
-- preparing the argument
local max_prime = p.max_prime(a)
local template_arg = ''
for i = 2, max_prime do
if u.is_prime(i) then
if i > 2 then template_arg = template_arg .. ' ' end
template_arg = template_arg .. (a[i] or 0)
end
end
s = s .. frame:expandTemplate{
title = 'Monzo',
args = {template_arg}
}
return s
end
-- parse a rational number
-- returns nil on failure
function p.parse(unparsed)
if type(unparsed) ~= 'string' then
return nil
end
-- rational form
local sign, n, _m, 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
n = tonumber(n)
if #sign > 0 then
n = -n
end
end
else
n = tonumber(n)
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 a = p.parse(unparsed)
if a == nil then
return '<span style="color:red;">Invalid rational number: ' .. unparsed .. '.</span>'
end
return p.as_ket(a, frame)
end
p.monzo = p.ket
return p
.svg){kind=icon}
