Module:Rational
Jump to navigation
Jump to search
local seq = require("Module:Sequence")
local utils = 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 = utils.signum(m) }
elseif m == 0 then
return { inf = true, sign = utils.signum(n) }
end
local sign = utils.signum(n) * utils.signum(m)
n = n * utils.signum(n)
m = m * utils.signum(m)
local n_factors = utils.prime_factorization_raw(n)
local m_factors = utils.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)
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`
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
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 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
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 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 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 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
-- 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
-- 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
-- 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 rational number in ket notation
-- NaN, infinity, zero values use special representations
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 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 not only_numbers and a.sign < 0 then
s = s .. "-"
end
-- 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 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