Xenharmonic Wiki

Module:Rational: Difference between revisions

← Older edit
Plumtree (talk | contribs)
autopatrolled
1,976 edits
m is_highly_composite(): basic implementation
ArrowHead294 (talk | contribs)
14,512 edits
mNo edit summary
 
(79 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 35: Line 38:
-- copy a rational number
-- 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 69: 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 83: 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 99: 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 116: Line 294:
-- compute a^b; b must be an integer
-- compute a^b; b must be an integer
function p.pow(a, b)
function p.pow(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
return nil
return nil
end
end
Line 144: Line 322:
end
end
local c = p.new(1)
local c = p.new(1)
for i = 1, math.abs(b) do
for _ = 1, math.abs(b) do
if b > 0 then
if b > 0 then
c = p.mul(c, a)
c = p.mul(c, a)
Line 155: Line 333:


-- compute a canonical representation of `a` modulo powers of `b`
-- 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)
function p.modulo_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 165: Line 345:
return p.copy(a)
return p.copy(a)
end
end
local neg_power = -1/0
local neg_power = -math.huge
local pos_power = 1/0
local pos_power = math.huge
for factor, power in pairs(b) do
for factor, power in pairs(b) do
if type(factor) == 'number' then
if type(factor) == "number" then
if (power > 0 and (a[factor] or 0) >= 0) or (power < 0 and (a[factor] or 0) <= 0) 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,
pos_power = math.min(pos_power, math.floor((a[factor] or 0) / power))
math.floor((a[factor] or 0) / power)
)
else
else
neg_power = math.max(neg_power,
neg_power = math.max(neg_power, -math.ceil(math.abs(a[factor] or 0) / math.abs(power)))
-math.ceil(math.abs(a[factor] or 0) / math.abs(power))
)
end
end
end
end
Line 192: Line 368:
-- 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 227: Line 403:
local gcd = { sign = 1 }
local gcd = { sign = 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 math.min(power, b[factor] or 0) > 0 then
if math.min(power, b[factor] or 0) > 0 then
gcd[factor] = math.min(power, b[factor])
gcd[factor] = math.min(power, b[factor])
Line 236: Line 412:
end
end
end
end
local a = p.div(a, gcd)
a = p.div(a, gcd)
local b = p.div(b, gcd)
b = p.div(b, gcd)
 
n_a, m_a = p.as_pair(a)
local n_a, m_a = p.as_pair(a)
n_b, m_b = p.as_pair(b)
local n_b, m_b = p.as_pair(b)
 
n = n_a * m_b + n_b * m_a
local n = n_a * m_b + n_b * m_a
m = m_a * m_b
local m = m_a * m_b
 
return p.mul(p.new(n, m), gcd)
return p.mul(p.new(n, m), gcd)
end
end
Line 305: 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 321: 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
if power < 0 then
return false
return false
Line 328: Line 534:
end
end
return true
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
end


-- check if an integer is highly composite
-- check if an integer is highly composite
function p.is_highly_composite(a)
function p.is_highly_composite(a)
if type(a) == 'number' then
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)
a = p.new(a)
end
end
Line 338: Line 798:
return false
return false
end
end
-- negative numbers are not highly composite
-- negative numbers are not superabundant
if a.sign == -1 then
if a.sign == -1 then
return false
return false
end
end
-- non-integers are not highly composite
-- non-integers are not superabundant
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
return false
return false
Line 350: Line 810:
end
end
end
end
local last_power = 1/0
local last_power = 1 / 0
local max_prime = p.max_prime(a)
local max_prime = p.max_prime(a)
local divisor_sum = p.new(1)
for i = 2, max_prime do
for i = 2, max_prime do
if u.is_prime(i) then
if utils.is_prime(i) then
-- factors must be the first N primes
-- factors must be the first N primes
if a[i] == nil then
if a[i] == nil then
Line 363: Line 824:
end
end
last_power = a[i]
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
end
end
Line 369: Line 831:
return p.eq(a, 1) or p.eq(a, 4) or p.eq(a, 36)
return p.eq(a, 1) or p.eq(a, 4) or p.eq(a, 36)
end
end
-- now we actually check whether it is highly composite
 
local n, _m = p.as_pair(a)
-- now we actually check whether it is superabundant
local divisors = p.divisors(a)
local n, _ = p.as_pair(a)
-- FIXME: inefficient; we should iterate factorisations of integers <n directly
 
for i = n - 1, 1, -1 do
-- precision is very important here
if p.divisors(i) >= divisors then
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
return false
end
end
diagram = utils.next_young_diagram(diagram)
end
end
return true
return true
end
end


-- find max prime in ket notation
-- 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)
function p.max_prime(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 390: Line 943:
end
end
local max_factor = 0
local max_factor = 0
for factor, power in pairs(a) do
for factor, _ in pairs(a) do
if type(factor) == 'number' then
if type(factor) == "number" then
if factor > max_factor then
if factor > max_factor then
max_factor = factor
max_factor = factor
Line 398: Line 951:
end
end
return max_factor
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
end


-- determine the number of rational divisors
-- determine the number of rational divisors
function p.divisors(a)
function p.divisors(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 410: Line 1,247:
local d = 1
local d = 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
d = d * (math.abs(power) + 1)
d = d * (math.abs(power) + 1)
end
end
Line 419: Line 1,256:
-- determine whether the rational number is +- p/q, where p, q are primes OR 1
-- determine whether the rational number is +- p/q, where p, q are primes OR 1
function p.is_prime_ratio(a)
function p.is_prime_ratio(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 428: Line 1,265:
local m_factors = 0
local m_factors = 0
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_factors = n_factors + 1
n_factors = n_factors + 1
Line 439: Line 1,276:
end
end


-- return the (n, m) pair as a Lua tuple
-- 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 460: 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 474: Line 1,373:
-- return a string ratio representation
-- 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 490: 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 max_prime = p.max_prime(a)
local max_prime = p.max_prime(a)
local template_arg = ''
local template_arg = ""
local template_arg_without_trailing_zeros = ""
local zeros_n = 0
for i = 2, max_prime do
for i = 2, max_prime do
if u.is_prime(i) then
if utils.is_prime(i) then
if i > 2 then template_arg = template_arg .. ' ' end
if i > 2 then
template_arg = template_arg .. " "
end
template_arg = template_arg .. (a[i] or 0)
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
s = s .. frame:expandTemplate{
if #template_arg == 0 then
title = 'Monzo',
template_arg = "0"
args = {template_arg}
end
}
if only_numbers then
s = s .. template_arg
else
s = s .. frame:expandTemplate({
title = "Monzo",
args = { template_arg },
})
end
return s
return s
end
end
Line 537: 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 550: 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 556: 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 567: 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"