Xenharmonic Wiki

Module:Rational: Difference between revisions

← Older edit
Plumtree (talk | contribs)
autopatrolled
1,976 edits
m as_FJS() changed
ArrowHead294 (talk | contribs)
14,512 edits
mNo edit summary
 
(34 intermediate revisions by 6 users not shown)
Line 1: Line 1:
local seq = require('Module:Sequence')
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
Line 9: Line 12:
return { nan = true }
return { nan = true }
elseif n == 0 then
elseif n == 0 then
return { zero = true, sign = u.signum(m) }
return { zero = true, sign = utils.signum(m) }
elseif m == 0 then
elseif m == 0 then
return { inf = true, sign = u.signum(n) }
return { inf = true, sign = utils.signum(n) }
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 32: 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
Line 52: Line 58:
local factor = 1
local factor = 1
local a = { sign = 1 }
local a = { sign = 1 }
for i in s:gmatch('[%-%d]+') do
for i in s:gmatch("%S+") do
local power = tonumber(i)
local power = tonumber(i)
if power == nil then return nil end
if power == nil then
return nil
end
 
-- find the next prime
-- find the next prime
factor = factor + 1
factor = factor + 1
while not u.is_prime(factor) do
while not utils.is_prime(factor) do
factor = factor + 1
factor = factor + 1
end
end
 
if power ~= 0 then
if power ~= 0 then
a[factor] = power
a[factor] = power
Line 79: Line 87:
table.insert(data, n)
table.insert(data, n)
local frac = p.from_continued_fraction(data)
local frac = p.from_continued_fraction(data)
if type(stop) == 'function' and stop(frac) then
if type(stop) == "function" and stop(frac) then
break
break
elseif type(stop) == 'number' and i >= stop then
elseif type(stop) == "number" and i >= stop then
break
break
end
end
table.insert(convergents, frac)
table.insert(convergents, frac)
x = x - n
x = x - n
if x == 0 then
break
end
x = 1 / x
x = 1 / x
i = i + 1
i = i + 1
Line 93: Line 104:


-- determine whether a rational number is a convergent or a semiconvergent to `x`
-- determine whether a rational number is a convergent or a semiconvergent to `x`
-- TODO: document how this works
function p.converges(a, x)
function p.converges(a, x)
local n_a, m_a = p.as_pair(a)
local _, m_a = p.as_pair(a)
local convergents = p.convergents(
local convergents = p.convergents(x, function(b)
x,
local _, m_b = p.as_pair(b)
function(b)
return m_b >= m_a * 10000
local n_b, m_b = p.as_pair(b)
end)
return m_b >= m_a * 10000
for _, b in ipairs(convergents) do
end
)
for i, b in ipairs(convergents) do
if p.eq(a, b) then
if p.eq(a, b) then
return 'convergent'
return "convergent"
end
end
end
end
 
for i = 2, #convergents - 1 do
for i = 2, #convergents - 1 do
local n_delta, m_delta = p.as_pair(convergents[i])
local n_delta, m_delta = p.as_pair(convergents[i])
Line 119: Line 128:
end
end
if p.eq(a, c) then
if p.eq(a, c) then
return 'semiconvergent'
return "semiconvergent"
end
end
end
end
Line 129: Line 138:
-- returns a table of matched expressions
-- returns a table of matched expressions
function p.find_S_expression(a)
function p.find_S_expression(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 145: Line 154:
local superparticular_indices = {}
local superparticular_indices = {}
local superparticular_ratios = {}
local superparticular_ratios = {}
for prime, k_array in pairs(seq.square_superparticulars) do
for _, k_array in pairs(seq.square_superparticulars) do
for i, k in ipairs(k_array) do
for _, k in ipairs(k_array) do
if k <= 1000 then
if k <= 1000 then
table.insert(superparticular_indices, k)
table.insert(superparticular_indices, k)
 
local Sk_num = p.pow(p.new(k), 2)
local Sk_num = p.pow(p.new(k), 2)
local Sk_den = p.sub(Sk_num, 1)
local Sk_den = p.mul(k - 1, k + 1)
local Sk = p.div(Sk_num, Sk_den)
local Sk = p.div(Sk_num, Sk_den)
superparticular_ratios[k] = Sk
superparticular_ratios[k] = Sk
Line 157: Line 166:
end
end
end
end
 
-- is it Sk?
-- is it Sk?
for _, k in ipairs(superparticular_indices) do
for _, k in ipairs(superparticular_indices) do
if p.eq(a, superparticular_ratios[k]) then
if p.eq(a, superparticular_ratios[k]) then
table.insert(expressions, 'S' .. k)
table.insert(expressions, "S" .. k)
end
end
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?
-- 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
for _, k in ipairs(superparticular_indices) do
Line 171: Line 180:
if r1 and r2 then
if r1 and r2 then
if p.eq(a, p.mul(r1, r2)) then
if p.eq(a, p.mul(r1, r2)) then
table.insert(expressions, 'S' .. k .. ' × S' .. (k + 1))
table.insert(expressions, "S" .. k .. " × S" .. (k + 1))
end
end
if p.eq(a, p.div(r1, r2)) then
if p.eq(a, p.div(r1, r2)) then
table.insert(expressions, 'S' .. k .. ' / S' .. (k + 1))
table.insert(expressions, "S" .. k .. " / S" .. (k + 1))
end
end
if p.eq(a, p.mul(p.pow(r1, 2), r2)) then
if p.eq(a, p.mul(p.pow(r1, 2), r2)) then
table.insert(expressions, 'S' .. k .. '<sup>2</sup> × S' .. (k + 1))
table.insert(expressions, "S" .. k .. "<sup>2</sup> × S" .. (k + 1))
end
end
if p.eq(a, p.mul(r1, p.pow(r2, 2))) then
if p.eq(a, p.mul(r1, p.pow(r2, 2))) then
table.insert(expressions, 'S' .. k .. ' * S' .. (k + 1) .. '<sup>2</sup>')
table.insert(expressions, "S" .. k .. " * S" .. (k + 1) .. "<sup>2</sup>")
end
end
end
end
end
end
 
-- is it Sk/S(k+2)?
-- is it Sk/S(k+2)?
for _, k in ipairs(superparticular_indices) do
for _, k in ipairs(superparticular_indices) do
Line 191: Line 200:
if r1 and r2 then
if r1 and r2 then
if p.eq(a, p.div(r1, r2)) then
if p.eq(a, p.div(r1, r2)) then
table.insert(expressions, 'S' .. k .. ' / S' .. (k + 2))
table.insert(expressions, "S" .. k .. " / S" .. (k + 2))
end
end
end
end
end
end
 
-- is it S(k-1)*Sk*S(k+1)?
-- is it S(k-1)*Sk*S(k+1)?
for _, k in ipairs(superparticular_indices) do
for _, k in ipairs(superparticular_indices) do
Line 203: Line 212:
if r1 and r2 and r3 then
if r1 and r2 and r3 then
if p.eq(a, p.mul(r1, p.mul(r2, 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))
table.insert(expressions, "S" .. (k - 1) .. " × S" .. k .. " × S" .. (k + 1))
end
end
end
end
end
end
 
return expressions
return expressions
end
end
Line 213: Line 222:
-- 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 238: 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 252: 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 268: 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 285: 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 313: 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 324: 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 334: 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 361: 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 396: 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 405: 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 474: 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)
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 510: Line 517:
-- 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 520: 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 532: Line 539:
function p.is_reduced(a, equave, large)
function p.is_reduced(a, equave, large)
equave = equave or 2
equave = equave or 2
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 548: Line 555:
end
end


-- determine whether a rational number represents a harmonic
-- determine whether a rational number represents a harmonic.
-- reduced: check for reduced harmonic instead.
function p.is_harmonic(a, reduced, large)
function p.is_harmonic(a, reduced, large)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 557: Line 565:
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 factor == 2 and reduced then
if factor == 2 and reduced then
-- do nothing
-- pass (ignore factors of 2 for reduced harmonic check)
elseif power < 0 then
elseif power < 0 then
return false
return false
Line 571: Line 579:
end
end


-- determine whether a rational number represents a subharmonic
-- determine whether a rational number represents a subharmonic.
-- reduced: check for reduced subharmonic instead.
function p.is_subharmonic(a, reduced, large)
function p.is_subharmonic(a, reduced, large)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 580: Line 589:
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 factor == 2 and reduced then
if factor == 2 and reduced then
-- do nothing
-- pass (ignore factors of 2 for reduced subharmonic check)
elseif power > 0 then
elseif power > 0 then
return false
return false
Line 596: Line 605:
-- determine whether a rational number is an integer power of another rational number
-- determine whether a rational number is an integer power of another rational number
function p.is_power(a)
function p.is_power(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 605: Line 614:
return false
return false
end
end
local function gcd(x, y)
 
if x < y then
x, y = y, x
end
while y > 0 do
x, y = y, x % y
end
return x
end
local total_power = nil
local total_power = nil
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 total_power then
if total_power then
total_power = gcd(total_power, math.abs(power))
total_power = utils._gcd(total_power, math.abs(power))
else
else
total_power = math.abs(power)
total_power = math.abs(power)
Line 630: Line 630:
-- determine whether a rational number is superparticular
-- determine whether a rational number is superparticular
function p.is_superparticular(a)
function p.is_superparticular(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 642: Line 642:
-- determine whether a rational number is a square superparticular
-- determine whether a rational number is a square superparticular
function p.is_square_superparticular(a)
function p.is_square_superparticular(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
if a.nan or a.inf or a.zero then
if a.nan or a.inf or a.zero or a.sign < 0 then
return false
return false
end
end
-- check the numerator
-- check the numerator
local k = { 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 power > 0 and power % 2 ~= 0 then
if power > 0 and power % 2 ~= 0 then
return false
return false
elseif power > 0 then
k[factor] = math.floor(power / 2 + 0.5)
end
end
end
end
end
end
return p.is_superparticular(a)
-- 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)
a = p.new(a)
end
end
if a.nan or a.inf or a.zero then
return false
-- nan, inf, zero, and negative numbers are not highly composite
end
if a.nan or a.inf or a.zero or a.sign == -1 then
-- negative numbers are not highly composite
if a.sign == -1 then
return false
return false
end
end
-- non-integers are not highly composite
-- non-integers are not highly composite
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 679: Line 707:
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)
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 698: Line 727:
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
-- now we actually check whether it is highly composite
local n, _m = p.as_pair(a)
local n, _ = p.as_pair(a)
 
-- precision is very important here
-- precision is very important here
local log2_n = 0
local log2_n = 0
Line 709: Line 738:
t = t * 2
t = t * 2
end
end
 
local divisors = p.divisors(a)
local divisors = p.divisors(a)
local diagram_size = log2_n
local diagram_size = log2_n
local diagram = {log2_n}
local diagram = { log2_n }
local primes = {2}
local primes = { 2 }
 
function eval_diagram(d)
local function eval_diagram(d)
while #d > #primes do
while #d > #primes do
local i = primes[#primes] + 1
local i = primes[#primes] + 1
while not u.is_prime(i) do
while not utils.is_prime(i) do
i = i + 1
i = i + 1
end
end
Line 725: Line 754:
local m = 1
local m = 1
for i = 1, #d do
for i = 1, #d do
for j = 1, d[i] do
for _ = 1, d[i] do
m = m * primes[i]
m = m * primes[i]
end
end
Line 731: Line 760:
return m
return m
end
end
 
-- iterate factorisations of some composite integers <n
-- iterate factorisations of some composite integers <n
while diagram do
while diagram do
Line 756: Line 785:
return false
return false
end
end
diagram = u.next_young_diagram(diagram)
diagram = utils.next_young_diagram(diagram)
end
end
return true
return true
Line 763: Line 792:
-- check if an integer is superabundant
-- check if an integer is superabundant
function p.is_superabundant(a)
function p.is_superabundant(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 775: Line 804:
-- non-integers are not superabundant
-- 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 781: 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)
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 795: Line 824:
end
end
last_power = a[i]
last_power = a[i]
divisor_sum = p.mul(
divisor_sum = p.mul(divisor_sum, p.div(p.sub(p.pow(i, a[i] + 1), 1), i - 1))
divisor_sum,
p.div(
p.sub(p.pow(i, a[i] + 1), 1),
i - 1
)
)
end
end
end
end
Line 808: 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 superabundant
-- now we actually check whether it is superabundant
local n, _m = p.as_pair(a)
local n, _ = p.as_pair(a)
 
-- precision is very important here
-- precision is very important here
local log2_n = 0
local log2_n = 0
Line 819: Line 842:
t = t * 2
t = t * 2
end
end
 
local SA_ratio = p.div(divisor_sum, a)
local SA_ratio = p.div(divisor_sum, a)
local diagram_size = log2_n
local diagram_size = log2_n
local diagram = {log2_n}
local diagram = { log2_n }
local primes = {2}
local primes = { 2 }
 
function eval_diagram(d)
local function eval_diagram(d)
while #d > #primes do
while #d > #primes do
local i = primes[#primes] + 1
local i = primes[#primes] + 1
while not u.is_prime(i) do
while not utils.is_prime(i) do
i = i + 1
i = i + 1
end
end
Line 835: Line 858:
local m = 1
local m = 1
for i = 1, #d do
for i = 1, #d do
for j = 1, d[i] do
for _ = 1, d[i] do
m = m * primes[i]
m = m * primes[i]
end
end
Line 841: Line 864:
return m
return m
end
end
 
-- iterate factorisations of some composite integers <n
-- iterate factorisations of some composite integers <n
while diagram do
while diagram do
Line 861: Line 884:
local diagram_divisor_sum = 1
local diagram_divisor_sum = 1
for i = 1, #diagram do
for i = 1, #diagram do
diagram_divisor_sum = p.mul(
diagram_divisor_sum =
diagram_divisor_sum,
p.mul(diagram_divisor_sum, p.div(p.sub(p.pow(primes[i], diagram[i] + 1), 1), primes[i] - 1))
p.div(
p.sub(p.pow(primes[i], diagram[i] + 1), 1),
primes[i] - 1
)
)
end
end
local diagram_SA_ratio = p.div(diagram_divisor_sum, eval_diagram(diagram))
local diagram_SA_ratio = p.div(diagram_divisor_sum, eval_diagram(diagram))
Line 873: Line 891:
return false
return false
end
end
diagram = u.next_young_diagram(diagram)
diagram = utils.next_young_diagram(diagram)
end
end
return true
return true
end
end


-- find max prime involved in the factorisation of a rational number
-- 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 887: 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 897: Line 953:
end
end


-- compute log2 of a rational number
-- convert a rational number to its size in octaves
-- equal to log2 of the rational number
function p.log(a, base)
function p.log(a, base)
base = base or 2
base = base or 2
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
if a.inf and a.sign > 0 then
if a.inf and a.sign > 0 then
return 1/0
return 1 / 0
end
end
if a.nan or a.inf then
if a.nan or a.inf then
Line 910: Line 967:
end
end
if a.zero then
if a.zero then
return -1/0
return -1 / 0
end
end
if a.sign < 0 then
if a.sign < 0 then
Line 917: Line 974:
local logarithm = 0
local logarithm = 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
logarithm = logarithm + power * math.log(factor) / math.log(base)
logarithm = logarithm + power * utils._log(factor, base)
end
end
end
end
return logarithm
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
end


-- FJS: x = a * 2^n : x >= 1, x < 2
-- FJS: x = a * 2^n : x >= 1, x < 2
local function red(a)
local function red(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 933: Line 1,035:
end
end
local b = p.copy(a)
local b = p.copy(a)
 
-- start with an approximation
-- start with an approximation
local log2 = p.log(b)
local log2 = p.log(b)
b = p.div(b, p.pow(2, math.floor(log2)))
b = p.div(b, p.pow(2, math.floor(log2)))
 
while p.lt(b, 1) do
while p.lt(b, 1) do
b = p.mul(b, 2)
b = p.mul(b, 2)
Line 985: Line 1,087:
-- FJS representation of a rational number
-- FJS representation of a rational number
-- might be a bit incorrect
-- might be a bit incorrect
-- TODO: confirm correctness
function p.as_FJS(a)
function p.as_FJS(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 996: Line 1,099:
local utonal = {}
local utonal = {}
for factor, power in pairs(a) do
for factor, power in pairs(a) do
if type(factor) == 'number' and factor > 3 then
if type(factor) == "number" and factor > 3 then
local comma = formal_comma(factor)
local comma = formal_comma(factor)
b = p.div(b, p.pow(comma, power))
b = p.div(b, p.pow(comma, power))
if power > 0 then
if power > 0 then
for i = 1, power do
for _ = 1, power do
table.insert(otonal, factor)
table.insert(otonal, factor)
end
end
else
else
for i = 1, -power do
for _ = 1, -power do
table.insert(utonal, factor)
table.insert(utonal, factor)
end
end
Line 1,012: Line 1,115:
table.sort(otonal)
table.sort(otonal)
table.sort(utonal)
table.sort(utonal)
 
local fifths = b[3] or 0
local fifths = b[3] or 0
 
local o = math.floor((fifths * 2 + 3) / 7)
local o = math.floor((fifths * 2 + 3) / 7)
local num = fifths * 11 + (b[2] or 0) * 7
local num = fifths * 11 + (b[2] or 0) * 7
Line 1,023: Line 1,126:
o = -o
o = -o
end
end
 
local num_mod = (num - u.signum(num)) % 7
local num_mod = (num - utils.signum(num)) % 7
local letter = ''
local letter = ""
if (num_mod == 0 or num_mod == 3 or num_mod == 4) and o == 0 then
if (num_mod == 0 or num_mod == 3 or num_mod == 4) and o == 0 then
letter = 'P'
letter = "P"
elseif o == 1 then
elseif o == 1 then
letter = 'M'
letter = "M"
elseif o == -1 then
elseif o == -1 then
letter = 'm'
letter = "m"
else
else
if o >= 0 then
if o >= 0 then
Line 1,040: Line 1,143:
if o > 0 then
if o > 0 then
while o > 0 do
while o > 0 do
letter = letter .. 'A'
letter = letter .. "A"
o = o - 2
o = o - 2
end
end
else
else
while o < 0 do
while o < 0 do
letter = letter .. 'd'
letter = letter .. "d"
o = o + 2
o = o + 2
end
end
end
if #letter >= 5 then
letter = #letter .. letter:sub(1, 1)
end
end
end
end
 
local FJS = letter .. num
local FJS = letter .. num
if #otonal > 0 then
if #otonal > 0 then
FJS = FJS .. '^{' .. table.concat(otonal, ',') .. '}'
FJS = FJS .. "^{" .. table.concat(otonal, ",") .. "}"
end
end
if #utonal > 0 then
if #utonal > 0 then
FJS = FJS .. '_{' .. table.concat(utonal, ',') .. '}'
FJS = FJS .. "_{" .. table.concat(utonal, ",") .. "}"
end
end
return FJS
return FJS
Line 1,063: Line 1,169:
-- determine log2 product complexity
-- determine log2 product complexity
function p.tenney_height(a)
function p.tenney_height(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 1,069: Line 1,175:
return nil
return nil
end
end
local n, m = p.as_pair(a)
local h = 0
return math.log(n * m) / math.log(2)
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
end


-- determine max complexity
-- determine log2 max complexity
function p.weil_height(a)
function p.weil_height(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 1,081: Line 1,192:
return nil
return nil
end
end
local n, m = p.as_pair(a)
local h1 = p.tenney_height(a)
return math.max(n, m)
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
end


-- determine product complexity
-- determine product complexity
function p.benedetti_height(a)
function p.benedetti_height(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 1,094: Line 1,229:
end
end
local n, m = p.as_pair(a)
local n, m = p.as_pair(a)
return n * m
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 1,107: 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 1,116: 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 1,125: 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 1,134: Line 1,274:
end
end
return n_factors <= 1 and m_factors <= 1
return n_factors <= 1 and m_factors <= 1
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
-- convert a rational number to cents
function p.cents(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 nil
end
local c = 0
for factor, power in pairs(a) do
if type(factor) == 'number' then
c = c + power * math.log(factor)
end
end
return c * 1200 / math.log(2)
end
end


-- return prime factorisation of a rational number
-- return prime factorisation of a rational number
function p.factorisation(a)
function p.factorisation(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
if a.nan or a.inf or a.zero or p.eq(a, 1) or p.eq(a, -1) then
if a.nan or a.inf or a.zero or p.eq(a, 1) or p.eq(a, -1) then
return 'n/a'
return "n/a"
end
end
local s = ''
local s = ""
if a.sign < 0 then
if a.sign < 0 then
s = s .. '-'
s = s .. "-"
end
end
local factors = {}
local factors = {}
for factor, power in pairs(a) do
for factor, _ in pairs(a) do
if type(factor) == 'number' then
if type(factor) == "number" then
table.insert(factors, factor)
table.insert(factors, factor)
end
end
Line 1,194: Line 1,296:
table.sort(factors)
table.sort(factors)
for i, factor in ipairs(factors) do
for i, factor in ipairs(factors) do
if i > 1 then s = s .. ' × ' end
if i > 1 then
s = s .. " × "
end
s = s .. factor
s = s .. factor
if a[factor] ~= 1 then
if a[factor] ~= 1 then
s = s .. '<sup>' .. a[factor] .. '</sup>'
s = s .. "<sup>" .. a[factor] .. "</sup>"
end
end
end
end
Line 1,205: Line 1,309:
-- return the subgroup generated by primes from a rational number's prime factorisation
-- return the subgroup generated by primes from a rational number's prime factorisation
function p.subgroup(a)
function p.subgroup(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
if p.eq(a, 1) then
if p.eq(a, 1) then
return '1'
return "1"
end
end
if a.nan or a.inf or a.zero or p.eq(a, -1) then
if a.nan or a.inf or a.zero or p.eq(a, -1) then
return 'n/a'
return "n/a"
end
end
local s = ''
local s = ""
local factors = {}
local factors = {}
for factor, power in pairs(a) do
for factor, _ in pairs(a) do
if type(factor) == 'number' then
if type(factor) == "number" then
table.insert(factors, factor)
table.insert(factors, factor)
end
end
Line 1,223: Line 1,327:
table.sort(factors)
table.sort(factors)
for i, factor in ipairs(factors) do
for i, factor in ipairs(factors) do
if i > 1 then s = s .. '.' end
if i > 1 then
s = s .. "."
end
s = s .. factor
s = s .. factor
end
end
if a.sign < 0 then
if a.sign < 0 then
s = '-1.' .. s
s = "-1." .. s
end
end
return s
return s
end
end


-- return the (n, m) pair as a Lua tuple
-- 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 1,253: 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 1,267: 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 1,285: Line 1,391:
-- return a rational number in subgroup ket notation
-- return a rational number in subgroup ket notation
function p.as_subgroup_ket(a, frame)
function p.as_subgroup_ket(a, frame)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
if a.nan or a.inf or a.zero or a.sign < 0 then
if a.nan or a.inf or a.zero or a.sign < 0 then
return 'n/a'
return "n/a"
end
end
local factors = {}
local factors = {}
for factor, power in pairs(a) do
for factor, _ in pairs(a) do
if type(factor) == 'number' then
if type(factor) == "number" then
table.insert(factors, factor)
table.insert(factors, factor)
end
end
end
end
table.sort(factors)
table.sort(factors)
local subgroup = '1'
local subgroup = "1"
if not p.eq(a, 1) then
if not p.eq(a, 1) then
subgroup = table.concat(factors, '.')
subgroup = table.concat(factors, ".")
end
end
 
local powers = {}
local powers = {}
for i, factor in ipairs(factors) do
for _, factor in ipairs(factors) do
table.insert(powers, a[factor])
table.insert(powers, a[factor])
end
end
local template_arg = '0'
local template_arg = "0"
if not p.eq(a, 1) then
if not p.eq(a, 1) then
template_arg = table.concat(powers, ' ')
template_arg = table.concat(powers, " ")
end
end
 
return subgroup .. ' ' .. frame:expandTemplate{
return subgroup .. " " .. frame:expandTemplate({
title = 'Monzo',
title = "Monzo",
args = {template_arg}
args = { template_arg },
}
})
end
end


-- return a rational number in ket notation
-- return a string of a rational number in monzo notation
-- NaN, infinity, zero values use special representations
-- calling Template: Monzo
function p.as_ket(a, frame)
function p.as_ket(a, frame, skip_many_zeros, only_numbers)
if type(a) == 'number' then
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)
a = p.new(a)
end
end
-- special case: NaN
if a.nan then
-- special cases
return 'NaN'
if a.nan or a.inf or a.zero or a.sign < 0 then
return "n/a"
end
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 = ''
-- regular case: positive finite ratio
if a.sign < 0 then
local s = ""
s = s .. '-'
 
end
-- 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 template_arg_without_trailing_zeros = ""
local zeros_n = 0
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
if i > 2 then
template_arg = template_arg .. ' '
template_arg = template_arg .. " "
end
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 (a[i] or 0) ~= 0 then
if zeros_n >= 4 then
if skip_many_zeros and zeros_n >= 4 then
template_arg = template_arg_without_trailing_zeros
template_arg = template_arg_without_trailing_zeros
if #template_arg > 0 then
if #template_arg > 0 then
template_arg = template_arg .. ' '
template_arg = template_arg .. " "
end
end
template_arg = template_arg .. '0<sup>' .. zeros_n .. '</sup> ' .. (a[i] or 0)
template_arg = template_arg .. "0<sup>" .. zeros_n .. "</sup> " .. (a[i] or 0)
end
end
zeros_n = 0
zeros_n = 0
Line 1,374: Line 1,471:
end
end
if #template_arg == 0 then
if #template_arg == 0 then
template_arg = '0'
template_arg = "0"
end
if only_numbers then
s = s .. template_arg
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 1,386: 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
-- removing whitespaces
unparsed = unparsed:gsub('%s', '')
unparsed = unparsed:gsub("%s", "")
-- removing <br> and <br/> tags
-- removing <br> and <br/> tags
unparsed = unparsed:gsub('<br/?>', '')
unparsed = unparsed:gsub("<br/?>", "")
 
-- length limit: very long strings are not converted into Lua numbers correctly
-- length limit: very long strings are not converted into Lua numbers correctly
local max_length = 15
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 1,433: 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"