Xenharmonic Wiki

Module:Rational: Difference between revisions

← Older editNewer edit →
Ganaram inukshuk (talk | contribs)
autopatrolled
4,810 edits
Added odd-limit function, as recommended by fredg999
Sintel (talk | contribs)
autopatrolled
1,547 edits
merge changes from dev
Line 1: Line 1:
local seq = require('Module:Sequence')
local seq = require("Module:Sequence")
local u = require('Module:Utils')
local utils = require("Module:Utils")
local p = {}
local p = {}


Line 9: Line 9:
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)
n = n * utils.signum(n)
m = m * u.signum(m)
m = m * utils.signum(m)
local n_factors = u.prime_factorization_raw(n)
local n_factors = utils.prime_factorization_raw(n)
local m_factors = u.prime_factorization_raw(m)
local m_factors = utils.prime_factorization_raw(m)
local factors = n_factors
local factors = n_factors
factors.sign = sign
factors.sign = sign
Line 32: Line 32:
-- 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 52:
local factor = 1
local factor = 1
local a = { sign = 1 }
local a = { sign = 1 }
for i in s:gmatch('%S+') 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 81:
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
Line 96: Line 98:


-- 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 122: Line 122:
end
end
if p.eq(a, c) then
if p.eq(a, c) then
return 'semiconvergent'
return "semiconvergent"
end
end
end
end
Line 132: Line 132:
-- 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 148: Line 148:
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.mul(k - 1, k + 1)
local Sk_den = p.mul(k - 1, k + 1)
Line 160: Line 160:
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 174: Line 174:
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 194: Line 194:
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 206: Line 206:
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 216: Line 216:
-- 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 241: Line 241:
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 255: Line 255:
-- 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 271: Line 271:
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 288: Line 288:
-- 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 316: Line 316:
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 328: Line 328:
-- compute a canonical representation of `a` modulo powers of `b`
-- compute a canonical representation of `a` modulo powers of `b`
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 337: Line 337:
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 364: Line 360:
-- 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 399: Line 395:
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 408: Line 404:
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 477: Line 473:
-- 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 513: Line 509:
-- 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 523: Line 519:
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 535: Line 531:
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 553: Line 549:
-- determine whether a rational number represents a harmonic
-- determine whether a rational number represents a harmonic
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 560: Line 556:
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 power < 0 then
-- do nothing
elseif power < 0 then
return false
return false
end
end
Line 576: Line 570:
-- determine whether a rational number represents a subharmonic
-- determine whether a rational number represents a subharmonic
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 583: Line 577:
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 power > 0 then
-- do nothing
elseif power > 0 then
return false
return false
end
end
Line 599: Line 591:
-- 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 611: Line 603:
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 = u._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 624: Line 616:
-- 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 636: Line 628:
-- 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
Line 643: Line 635:
end
end
-- check the numerator
-- check the numerator
local k = {sign = 1}
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
Line 660: Line 652:
-- 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
Line 672: Line 664:
-- 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 678: Line 670:
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 697: Line 689:
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 708: Line 700:
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 724: Line 716:
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 730: Line 722:
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 755: Line 747:
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 762: Line 754:
-- 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 774: Line 766:
-- 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 780: Line 772:
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 794: Line 786:
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 807: Line 793:
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 818: Line 804:
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 834: Line 820:
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 840: Line 826:
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 860: Line 846:
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 872: Line 853:
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 880: Line 861:
-- (a.k.a. prime limit or harmonic class) of a rational number
-- (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 868:
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 900: Line 881:
-- For a ratio p/q, this is simply max(p, q) where powers of 2 are ignored
-- For a ratio p/q, this is simply max(p, q) where powers of 2 are ignored
function p.odd_limit(a)
function p.odd_limit(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 918: Line 899:
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 928: Line 909:
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 935: Line 916:
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 * u._log(factor, base)
logarithm = logarithm + power * utils._log(factor, base)
end
end
end
end
Line 944: Line 925:
-- convert a rational number to its size in cents
-- convert a rational number to its size in cents
function p.cents(a)
function p.cents(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 951: Line 932:
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.zero then
if a.zero then
return -1/0
return -1 / 0
end
end


local c = 0
local c = 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
c = c + power * u._log(factor, base)
c = c + power * utils.log2(factor)
end
end
end
end
Line 969: Line 950:
function p.hz(a, base)
function p.hz(a, base)
base = base or 440
base = base or 440
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 980: Line 961:
local log_hz = math.log(base)
local log_hz = math.log(base)
for factor, power in pairs(a) do
for factor, power in pairs(a) do
if type(factor) == 'number' then
if type(factor) == "number" then
log_hz = log_hz + power * math.log(factor)
log_hz = log_hz + power * math.log(factor)
end
end
Line 989: Line 970:
-- 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 996: Line 977:
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 1,048: Line 1,029:
-- 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 1,059: Line 1,041:
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,075: Line 1,057:
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,086: Line 1,068:
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,103: Line 1,085:
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
end
if #letter >= 5 then
if #letter >= 5 then
letter = (#letter) .. letter:sub(1, 1)
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,129: Line 1,111:
-- 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,137: Line 1,119:
local h = 0
local h = 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
h = h + math.abs(power) * u._log(factor, 2)
h = h + math.abs(power) * utils.log2(factor)
end
end
end
end
Line 1,146: Line 1,128:
-- determine log2 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,152: Line 1,134:
return nil
return nil
end
end
local h1 = p.tenney_height (a)
local h1 = p.tenney_height(a)
local h2 = 0
local h2 = 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
h2 = h2 + power * u._log(factor, 2)
h2 = h2 + power * utils.log2(factor)
end
end
end
end
h2 = math.abs (h2)
h2 = math.abs(h2)
return h1 + h2
return h1 + h2
end
end
Line 1,165: Line 1,147:
-- determine sopfr complexity
-- determine sopfr complexity
function p.wilson_height(a)
function p.wilson_height(a)
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
Line 1,173: Line 1,155:
local h = 0
local h = 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
h = h + math.abs(power) * factor
h = h + math.abs(power) * factor
end
end
Line 1,182: Line 1,164:
-- 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,199: Line 1,181:
-- 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,207: Line 1,189:
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,216: Line 1,198:
-- 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,225: Line 1,207:
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,238: Line 1,220:
-- 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,256: Line 1,238:
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,267: Line 1,251:
-- 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,285: Line 1,269:
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,315: Line 1,301:
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,329: Line 1,315:
-- 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 ('%d%s%d'):format(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,347: Line 1,333:
-- 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


Line 1,387: Line 1,373:
end
end
only_numbers = only_numbers or false
only_numbers = only_numbers or false
if type(a) == 'number' then
if type(a) == "number" then
a = p.new(a)
a = p.new(a)
end
end
-- special case: NaN
-- special case: NaN
if a.nan then
if a.nan then
return 'NaN'
return "NaN"
end
end
-- special case: infinity
-- special case: infinity
if a.inf then
if a.inf then
local sign = '+'
local sign = "+"
if a.sign < 0 then
if a.sign < 0 then
sign = '-'
sign = "-"
end
end
return sign .. '∞'
return sign .. "∞"
end
end
-- special case: zero
-- special case: zero
if a.zero then
if a.zero then
return '0'
return "0"
end
end
-- regular case: not NaN, not infinity, not zero
-- regular case: not NaN, not infinity, not zero
 
local s = ''
local s = ""
if not only_numbers and a.sign < 0 then
if not only_numbers and a.sign < 0 then
s = s .. '-'
s = s .. "-"
end
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 skip_many_zeros and 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,440: Line 1,426:
end
end
if #template_arg == 0 then
if #template_arg == 0 then
template_arg = '0'
template_arg = "0"
end
end
if only_numbers then
if only_numbers then
s = s .. template_arg
s = s .. template_arg
else
else
s = s .. frame:expandTemplate{
s = s .. frame:expandTemplate({
title = 'Monzo',
title = "Monzo",
args = {template_arg}
args = { template_arg },
}
})
end
end
return s
return s
Line 1,456: Line 1,442:
-- 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,503: Line 1,489:
-- 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 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>'
return '<span style="color:red;">Invalid rational number: ' .. unparsed .. ".</span>"
end
end
return p.as_ket(a, frame)
return p.as_ket(a, frame)
Retrieved from "https://en.xen.wiki/w/Module:Rational"