Xenharmonic Wiki

Module:Rational: Difference between revisions

← Older edit
Plumtree (talk | contribs)
autopatrolled, emailconfirmed
1,976 edits
mNo edit summary
ArrowHead294 (talk | contribs)
14,512 edits
mNo edit summary
 
(89 intermediate revisions by 6 users not shown)
Line 1: Line 1:
local u = require('Module:Utils')
local p = {}
local p = {}


-- construct a rational number n/m
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)
function p.new(n, m)
m = m or 1
m = m or 1
if n == 0 then
if n == 0 and m == 0 then
if m == 0 then
return { nan = true }
return { nan = true }
elseif n == 0 then
else
return { zero = true, sign = utils.signum(m) }
return { zero = true, sign = u.signum(m) }
elseif m == 0 then
end
return { inf = true, sign = utils.signum(n) }
else
if m == 0 then
return { inf = true, sign = u.signum(n) }
end
end
end
local sign = u.signum(n) * u.signum(m)
local sign = utils.signum(n) * utils.signum(m)
n = n * u.signum(n)
-- ensure n and m are positive
m = m * u.signum(m)
n = n * utils.signum(n)
local n_factors = u.prime_factorization_raw(n)
m = m * utils.signum(m)
local m_factors = u.prime_factorization_raw(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
local factors = n_factors
factors.sign = sign
factors.sign = sign
-- subtract the factors of m from the factors of n
for factor, power in pairs(m_factors) do
for factor, power in pairs(m_factors) do
factors[factor] = factors[factor] or 0
factors[factor] = factors[factor] or 0
factors[factor] = factors[factor] - power
factors[factor] = factors[factor] - power
if factors[factor] == 0 then
if factors[factor] == 0 then
factors[factor] = nil
factors[factor] = nil -- clear the zeros
end
end
end
end
Line 33: Line 36:
end
end


-- copy a rational number
function p.copy(a)
function p.copy(a)
b = {}
local b = {}
for factor, power in pairs(a) do
for factor, power in pairs(a) do
b[factor] = power
b[factor] = power
end
end
return b
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
end


-- multiply two rational numbers; integers are also allowed
-- multiply two rational numbers; integers are also allowed
function p.mul(a, b)
function p.mul(a, b)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
if type(b) == 'number' then
if type(b) == "number" then
b = p.new(b)
b = p.new(b)
end
end
Line 68: Line 247:
local c = p.copy(a)
local c = p.copy(a)
for factor, power in pairs(b) do
for factor, power in pairs(b) do
if type(factor) == 'number' then
if type(factor) == "number" then
c[factor] = c[factor] or 0
c[factor] = c[factor] or 0
c[factor] = c[factor] + power
c[factor] = c[factor] + power
Line 82: Line 261:
-- compute 1/a for a rational number a; integers are also allowed
-- compute 1/a for a rational number a; integers are also allowed
function p.inv(a)
function p.inv(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 98: Line 277:
end
end
-- regular case: not NaN, not infinity, not zero
-- regular case: not NaN, not infinity, not zero
b = {}
local b = {}
for factor, power in pairs(a) do
for factor, power in pairs(a) do
if type(factor) == 'number' then
if type(factor) == "number" then
b[factor] = -power
b[factor] = -power
end
end
Line 111: Line 290:
function p.div(a, b)
function p.div(a, b)
return p.mul(a, p.inv(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
end


-- add two rational numbers; integers are also allowed
-- add two rational numbers; integers are also allowed
function p.add(a, b)
function p.add(a, b)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
if type(b) == 'number' then
if type(b) == "number" then
b = p.new(b)
b = p.new(b)
end
end
 
-- special case: NaN
-- special case: NaN
if a.nan or b.nan then
if a.nan or b.nan then
Line 140: Line 393:
return { inf = true, sign = b.sign }
return { inf = true, sign = b.sign }
end
end
-- regular case: both not NaN, not infinities
-- special case: one is zero
if a.zero then
n_a, m_a = p.as_pair(a)
return p.copy(b)
n_b, m_b = p.as_pair(b)
end
if b.zero then
n = n_a * m_b + n_b * m_a
return p.copy(a)
m = m_a * m_b
end
-- regular case: both not NaN, not infinities, not zeros
return p.new(n, m)
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
end


Line 154: Line 427:
function p.sub(a, b)
function p.sub(a, b)
return p.add(a, p.mul(b, -1))
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
end


Line 198: Line 481:
-- determine whether a rational number is equal to another; integers are also allowed
-- determine whether a rational number is equal to another; integers are also allowed
function p.eq(a, b)
function p.eq(a, b)
local c = p.sub(a, b)
if type(a) == "number" then
return c.zero
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
end


-- determine whether a rational number is integer
-- determine whether a rational number is integer
function p.is_int(a)
function p.is_int(a)
if type(a) == 'number' then
if type(a) == "number" then
return true
return true
end
end
Line 214: Line 527:
end
end
for factor, power in pairs(a) do
for factor, power in pairs(a) do
if type(factor) == 'number' then
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
if power < 0 then
return false
return false
end
end
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
end
return true
return true
end
end


-- return the (n, m) pair as a Lua tuple
-- 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)
function p.as_pair(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 244: Line 1,359:
local m = 1
local m = 1
for factor, power in pairs(a) do
for factor, power in pairs(a) do
if type(factor) == 'number' then
if type(factor) == "number" then
if power > 0 then
if power > 0 then
n = n * (factor ^ power)
n = n * (factor ^ power)
else
else
m = m * (factor ^ (-power))
m = m * (factor ^ -power)
end
end
end
end
Line 256: Line 1,371:
end
end


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


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


Line 273: Line 1,389:
end
end


-- return a rational number in ket notation
-- return a rational number in subgroup ket notation
-- NaN, infinity, zero values use special representations
function p.as_subgroup_ket(a, frame)
function p.as_ket(a, frame)
if type(a) == "number" then
if type(a) == 'number' then
a = p.new(a)
a = p.new(a)
end
end
-- special case: NaN
if a.nan or a.inf or a.zero or a.sign < 0 then
if a.nan then
return "n/a"
return 'NaN'
end
end
-- special case: infinity
local factors = {}
if a.inf then
for factor, _ in pairs(a) do
local sign = '+'
if type(factor) == "number" then
if a.sign < 0 then
table.insert(factors, factor)
sign = '-'
end
end
return sign .. '∞'
end
end
-- special case: zero
table.sort(factors)
if a.zero then
local subgroup = "1"
return '0'
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
end
-- regular case: not NaN, not infinity, not zero
local s = ''
-- special cases
if a.sign < 0 then
if a.nan or a.inf or a.zero or a.sign < 0 then
s = s .. '-'
return "n/a"
end
end
-- regular case: positive finite ratio
local s = ""
-- preparing the argument
-- preparing the argument
local largest_prime = -1
local max_prime = p.max_prime(a)
for factor, power in pairs(a) do
local template_arg = ""
if type(factor) == 'number' then
local template_arg_without_trailing_zeros = ""
if factor > largest_prime then
local zeros_n = 0
largest_prime = factor
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
end
end
end
local template_arg = ''
if #template_arg == 0 then
for i = 2, largest_prime do
template_arg = "0"
if u.is_prime(i) then
end
if i > 2 then template_arg = template_arg .. ' ' end
if only_numbers then
template_arg = template_arg .. (a[i] or 0)
s = s .. template_arg
end
else
s = s .. frame:expandTemplate({
title = "Monzo",
args = { template_arg },
})
end
end
s = s .. frame:expandTemplate{
title = 'Monzo',
args = {template_arg}
}
return s
return s
end
end
Line 327: Line 1,487:
-- returns nil on failure
-- returns nil on failure
function p.parse(unparsed)
function p.parse(unparsed)
if type(unparsed) ~= 'string' then
if type(unparsed) ~= "string" then
return nil
return nil
end
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
-- rational form
local sign, n, _m, m = unparsed:match('^%s*(%-?)%s*(%d+)%s*(/%s*(%d+))%s*$')
local sign, n, _, m = unparsed:match("^%s*(%-?)%s*(%d+)%s*(/%s*(%d+))%s*$")
if n == nil then
if n == nil then
-- integer form
-- integer form
sign, n = unparsed:match('^%s*(%-?)%s*(%d+)%s*$')
sign, n = unparsed:match("^%s*(%-?)%s*(%d+)%s*$")
if n == nil then
if n == nil then
-- parsing failure
-- parsing failure
Line 340: Line 1,508:
else
else
m = 1
m = 1
if #n > max_length then
return nil
end
n = tonumber(n)
n = tonumber(n)
if #sign > 0 then
if #sign > 0 then
Line 346: Line 1,517:
end
end
else
else
if #n > max_length then
return nil
end
n = tonumber(n)
n = tonumber(n)
if #m > max_length then
return nil
end
m = tonumber(m)
m = tonumber(m)
if #sign > 0 then
if #sign > 0 then
Line 357: Line 1,534:
-- a version of as_ket() that can be {{#invoke:}}d
-- a version of as_ket() that can be {{#invoke:}}d
function p.ket(frame)
function p.ket(frame)
local unparsed = frame.args[1] or '1'
local unparsed = frame.args[1] or "1"
local result = ""
local a = p.parse(unparsed)
local a = p.parse(unparsed)
if a == nil then
if a == nil then
return '<span style="color:red;">Invalid rational number: ' .. unparsed .. '.</span>'
result = '{{error|Invalid rational number: ' .. unparsed .. ".}}"
else
result = p.as_ket(a, frame)
end
end
return p.as_ket(a, frame)
return frame:preprocess(result)
end
end
p.monzo = p.ket
p.monzo = p.ket


return p
return p
Retrieved from "https://en.xen.wiki/w/Module:Rational"