Module:No-1s limits

From Xenharmonic Wiki
Jump to navigation Jump to search

Documentation transcluded from /doc
Icon-Todo.png Todo: Documentation

local rat = require('Module:Rational')
local ET = require('Module:ET')
local p = {}

-- compute all positive ratios n/m with n and m <= q modulo powers of equave
-- previous: already computed ratios for q-1
function p.limit_modulo_equave(q, equave, previous)
	equave = equave or 2
	local ratios = {}
	if rat.eq(rat.modulo_mul(rat.new(q, 1), equave), 1) then return ratios end
	if previous then
		for n = 2, q do
			if not rat.eq(rat.modulo_mul(rat.new(n, 1), equave), 1) then
				local a = rat.new(n, q)
				a = rat.modulo_mul(a, equave)
				local a_key = rat.as_ratio(a)
				
				local b = rat.new(q, n)
				b = rat.modulo_mul(b, equave)
				local b_key = rat.as_ratio(b)
				
				if previous[a_key] == nil then
					ratios[a_key] = a
				end
				if previous[b_key] == nil then
					ratios[b_key] = b
				end
			end
		end
	else
		for n = 2, q do
			if not rat.eq(rat.modulo_mul(rat.new(n, 1), equave), 1) then
				for m = 2, q do
					if not rat.eq(rat.modulo_mul(rat.new(n, 1), equave), 1) then
						local a = rat.new(n, m)
						a = rat.modulo_mul(a, equave)
						local key = rat.as_ratio(a)
						ratios[key] = a
					end
				end
			end
		end
	end
	return ratios
end

-- compute q-integer limit
-- if a function `norm` and a number `max_norm` are provided, the output will be additionally restricted
function p.integer_limit(q, norm, max_norm)
	local check_norm = type(norm) == 'function' and type(max_norm) == 'number'
	local ratios = {}
	for n = 1, q do
		for m = 1, q do
			local a = rat.new(n, m)
			if not check_norm or norm(a) <= max_norm then
				local key = rat.as_ratio(a)
				ratios[key] = a
			end
		end
	end
	return ratios
end

-- check additive consistency for a set of ratios (modulo powers of equave):
--   approx(a*b) = approx(a) + approx(b) forall a, b: a, b, ab in ratios
-- `distinct`: whether distinct ratios are required to be mapped to distinct approximations
-- `previous`: already computed ratios for the previous iteraton
function p.additively_consistent(et, ratios, distinct, previous)
	distinct = distinct or false
	previous = previous or {}
	if distinct then
		local approx_set = {}
		for a_key, a in pairs(previous) do
			local a_approx = ET.approximate(et, rat.as_float(a)) % et.size
			if approx_set[a_approx] then
				if not rat.eq(rat.modulo_mul(rat.div(a, approx_set[a_approx]), et.equave), 1) then
					mw.log(a_key .. ' -> ' .. a_approx .. ': conflict!')
					return false
				end
			end
			approx_set[a_approx] = a
			mw.log(a_key .. ' -> ' .. a_approx)
		end
		for a_key, a in pairs(ratios) do
			local a_approx = ET.approximate(et, rat.as_float(a)) % et.size
			if approx_set[a_approx] then
				if not rat.eq(rat.modulo_mul(rat.div(a, approx_set[a_approx]), et.equave), 1) then
					mw.log(a_key .. ' -> ' .. a_approx .. ': conflict!')
					return false
				end
			end
			approx_set[a_approx] = a
			mw.log(a_key .. ' -> ' .. a_approx)
		end
	end
	if type(distinct) == 'number' then
		return true
	end
	local previous_ordered = {}
	for a_key, a in pairs(previous) do
		table.insert(previous_ordered, a)
	end
	local ratios_ordered = {}
	for a_key, a in pairs(ratios) do
		table.insert(ratios_ordered, a)
	end
	for i, a in ipairs(ratios_ordered) do
		local a_approx = ET.approximate(et, rat.as_float(a))
		for j, b in ipairs(previous_ordered) do
			local b_approx = ET.approximate(et, rat.as_float(b))
			
			local c = rat.mul(a, b)
			local c_approx = ET.approximate(et, rat.as_float(c))
			
			c = rat.modulo_mul(c, et.equave)
			local c_key = rat.as_ratio(c)
			if previous[c_key] or ratios[c_key] then
				if c_approx ~= a_approx + b_approx then
					mw.log('a = ' .. rat.as_ratio(a) .. '; b = ' .. rat.as_ratio(b) .. '; ab = ' .. c_key)
					mw.log(a_approx .. ' + ' .. b_approx .. ' != ' .. c_approx)
					return false
				end
			end
		end
		for j, b in ipairs(ratios_ordered) do
			if i <= j then
				local b_approx = ET.approximate(et, rat.as_float(b))
				
				local c = rat.mul(a, b)
				local c_approx = ET.approximate(et, rat.as_float(c))
				
				c = rat.modulo_mul(c, et.equave)
				local c_key = rat.as_ratio(c)
				if previous[c_key] or ratios[c_key] then
					if c_approx ~= a_approx + b_approx then
						mw.log('a = ' .. rat.as_ratio(a) .. '; b = ' .. rat.as_ratio(b) .. '; ab = ' .. c_key)
						mw.log(a_approx .. ' + ' .. b_approx .. ' != ' .. c_approx)
						return false
					end
				end
			end
		end
	end
	return true
end

-- find additive consistency limit
-- returns nil when at least `max_n`
-- `distinct`: whether distinct ratios are required to be mapped to distinct approximations
-- - if an integer, it is the regular consistency limit already known
function p.consistency_limit(et, distinct, max_n)
	if et.size == 0 then
		-- the answer is known already
		return '∞'
	end
	max_n = max_n or 1/0
	distinct = distinct or false
	local n = 1
	local last_n = 1
	local previous = {}
	while true do
		if type(distinct) == 'number' and n > distinct then
			return last_n
		end
		local ratios = p.limit_modulo_equave(n, et.equave, previous)
		for key, ratio in pairs(ratios) do
			mw.log('step ' .. n .. ': ' .. key)
		end
		if next(ratios) ~= nil then
			local consistent = p.additively_consistent(et, ratios, distinct, previous)
			if not consistent then
				mw.log('Not consistent at step ' .. n .. ', returning ' .. last_n)
				return last_n
			end
			for key, ratio in pairs(ratios) do
				previous[key] = ratio
			end
			last_n = n
		end
		n = n + 1
		if n > max_n then
			return nil
		end
	end
end

-- find vals generate no-1s consistency intervals
-- returns string with wart notation
-- limit must meet no-1s odd-limit consistency
function p.vals(et, limit)
	local primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43}
	local wart = 'abcdefghijklmn'
	local patent_val = {}
	local wart_dir = {}
	local minprime
	local val = {}
	local rets = {}
	for i = 1, #primes do
		patent_val[i] = ET.approximate(et, primes[i])
		if ET.cents(et, patent_val[i]) < math.log(primes[i]) / math.log(2) * 1200 then
			wart_dir[i] = -1
		else
			wart_dir[i] = 1
		end
	end
	if et.equave == 2 then
		minprime = 3
	else
		minprime = 2
	end
	local previ = 1
	for i = 1, #primes do
		local ip = primes[i]
		if ip > limit then break end
		if ip == et.equave then
			val[i] = patent_val[i]
		elseif ip == minprime then
			val[i] = patent_val[i]
			previ = i
		else
			val[i] = val[previ] + ET.approximate(et, ip / primes[previ])
			previ = i
		end
	end

	rets[1] = ''
	if et.equave ~= 2 then
		local i = table.find(primes, et.equave)
		if i then
			rets[1] = rets[1] .. string.sub(wart, i, i)
		else
			rets[1] = rets[1] .. 'q'
		end
	end
	rets[1] = rets[1] .. et.size
	if limit < minprime * minprime then
		local nonpatent = false
		local binded = false
		local bias = 0
		for i = 1, #val do
			local ip = primes[i]
			if ip == et.equave then
			elseif ip == minprime then
				bias = wart_dir[i]
			else
				if val[i] ~= patent_val[i] then
					nonpatent = true
				elseif wart_dir[i] ~= bias then
					binded = true
				end
			end
		end
		if nonpatent or not binded then
			rets[2] = '' .. rets[1]
			for i = 1, #val do
				local ip = primes[i]
				local difference = (val[i] - patent_val[i]) * wart_dir[i]
				if difference > 0 then
					rets[1] = rets[1] .. string.rep(string.sub(wart, i, i), difference * 2)
				elseif difference < 0 then
					rets[1] = rets[1] .. string.rep(string.sub(wart, i, i), -difference * 2 - 1)
				end
			end
			for i = 1, #val do
				local ip = primes[i]
				if ip ~= et.equave then val[i] = val[i] - bias end
				local difference = (val[i] - patent_val[i]) * wart_dir[i]
				if difference > 0 then
					rets[2] = rets[2] .. string.rep(string.sub(wart, i, i), difference * 2)
				elseif difference < 0 then
					rets[2] = rets[2] .. string.rep(string.sub(wart, i, i), -difference * 2 - 1)
				end
			end
			if bias > 0 then rets = {rets[2], rets[1]} end
		end
	else
		for i = 1, #val do
			local ip = primes[i]
			if ip * minprime <= limit then
				if val[i] ~= patent_val[i] then
					return '(error)'
				end
			end
			local difference = (val[i] - patent_val[i]) * wart_dir[i]
			if difference > 0 then
				rets[1] = rets[1] .. string.rep(string.sub(wart, i, i), difference * 2)
			elseif difference < 0 then
				rets[1] = rets[1] .. string.rep(string.sub(wart, i, i), -difference * 2 - 1)
			end
		end
	end
	return table.concat (rets, ", ")
end

function p.for_small_edos (frame)
	local lines = tonumber (frame.args['lines'])
	local t_body = {}
	local val = {}
	for i = 1, lines do
		local et = ET.parse('' .. i .. 'edo')
		local limit = p.consistency_limit(et, false, 43)
		local vals = p.vals(et, limit)
		t_body[i] = string.format ("|-\n| %s\n| %s\n| %s", i, limit, vals)
	end

	return "{| class=\"wikitable center-all mw-collapsible mw-collapsed\"\n" ..
	"|-\n" ..
	"! EDO\n" ..
	"! No-1s consistency limit\n" ..
	"! Associated vals\n" ..
	table.concat (t_body, "\n") ..
	"\n|}"
end

return p