View source for Module:Q-odd-limit intervals

bit32 = require("bit32")
utils = require("Module:Utils")
ET = require("Module:ET")
limits = require("Module:Limits")
yesno = require("Module:yesno")

local p = {}

local PRIME_LIST = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61}

local function table_filter(t, thres)
	for i = 1, #t do
		if t[i] > thres then
			return {unpack(t, 1, i - 1)}
		end
	end
	return t
end

local function inner_product(val, monzo)
	local result = 0
	for i = 1, #val do
		result = result + val[i] * monzo[i]
	end
	return result
end

local function monzo2ratio(monzo, subgroup)
	local num = 1
	local den = 1
	for i = 1, #subgroup do
		if monzo[i] > 0 then
			num = num * subgroup[i]^monzo[i]
		elseif monzo[i] < 0 then
			den = den * subgroup[i]^(-monzo[i])
		end
	end
	return {num = num, den = den}
end

local function ratio2monzo(ratio, subgroup)
	local monzo = {}
	for i = 1, #subgroup do
		monzo[i] = 0
		while ratio.num % subgroup[i] == 0 do
			monzo[i] = monzo[i] + 1
			ratio.num = ratio.num / subgroup[i]
		end while ratio.den % subgroup[i] == 0 do
			monzo[i] = monzo[i] - 1
			ratio.den = ratio.den / subgroup[i]
		end
	end
	return monzo
end

local function monzo2cent(monzo, subgroup)
	local jip = {}
	for i = 1, #subgroup do
		jip[i] = 1200 * utils._log(subgroup[i], 2)
	end
	return inner_product(jip, monzo)
end

local function ratio_8ve_reduction(ratio)
	local oct = math.floor(utils._log(ratio.num / ratio.den, 2))
	if oct > 0 then
		ratio.den = ratio.den * 2^oct
	elseif oct < 0 then
		ratio.num = ratio.num * 2^(-oct)
	end
	return ratio
end

local function odd_limit_monzo_list_gen(limit)
	local monzo_list = {}
	local subgroup = table_filter(PRIME_LIST, limit)
	for num = 1, limit, 2 do
		for den = 1, num, 2 do
			if utils._gcd(num, den) == 1 then
				ratio = ratio_8ve_reduction({num = num, den = den})
				table.insert(monzo_list, ratio2monzo(ratio, subgroup))
			end
		end
	end
	return monzo_list
end

local function find_error(val, subgroup, monzo_list)
	local step_size = 1200 / val[1]
	local true_size
	local approx_size
	local error_list = {}
	for i = 1, #monzo_list do
		ratio = monzo2ratio(monzo_list[i], subgroup)
		comp = {2 * ratio.den, ratio.num}
		true_size = monzo2cent(monzo_list[i], subgroup)
		approx_size = step_size * inner_product(val, monzo_list[i])
		nearest_size = math.floor(true_size / step_size + 0.5)
		error_abs_direct = math.abs(nearest_size * step_size - true_size)
		error_abs_val = math.abs(approx_size - true_size)
		error_rel_direct = 100 * error_abs_direct / step_size
		error_rel_val = 100 * error_abs_val / step_size
		error_list[i] = {
			ratio = ratio, 
			comp = comp, 
			error_abs_direct = error_abs_direct,
			error_abs_val = error_abs_val, 
			error_rel_direct = error_rel_direct,
			error_rel_val = error_rel_val
		}
	end
	
	table.sort(error_list, function(a, b)
		return a.error_abs_val < b.error_abs_val
	end)
	
	return error_list
end

local function approx(steps, subgroup, monzo_list, t_title, tag, approx_t)
	local t_body_direct = {}
	local errlist_direct = {}
	local t_body_val = {}
	local errlist_val = {}
	local val = {}
	local bgclr, result

	for i = 1, #subgroup do
		val[i] = utils._round_dec(steps*utils._log(subgroup[i], 2))
	end
	error_list = find_error(val, subgroup, monzo_list)
	for i = 1, #error_list do
		ratiocomp = string.format("%d/%d, %d/%d",
			error_list[i].ratio.num, error_list[i].ratio.den,
			2 * error_list[i].ratio.den, error_list[i].ratio.num)
		eav = error_list[i].error_abs_val
		erv = error_list[i].error_rel_val
		ead = error_list[i].error_abs_direct
		erd = error_list[i].error_rel_direct
		error_abs_val = string.format("%.3f", eav)
		error_rel_val = string.format("%.1f", erv)
		error_abs_direct = string.format("%.3f", ead)
		error_rel_direct = string.format("%.1f", erd)
		
		if bit32.band(error_list[i].ratio.den, error_list[i].ratio.den - 1) == 0
		and utils.table_contains(subgroup, error_list[i].ratio.num) then -- check power of 2 for den and prime for num
			ratiocomp = "'''" .. ratiocomp .. "'''"
			error_abs_val = "'''" .. error_abs_val .. "'''"
			error_rel_val = "'''" .. error_rel_val .. "'''"
			error_abs_direct = "'''" .. error_abs_direct .. "'''"
			error_rel_direct = "'''" .. error_rel_direct .. "'''"
			bgclr = ""
		end
		if error_list[i].error_rel_val > 50 then
			ratiocomp = "''" .. ratiocomp .. "''"
			error_abs_val = "''" .. error_abs_val .. "''"
			error_rel_val = "''" .. error_rel_val .. "''"
			error_abs_direct = "''" .. error_abs_direct .. "''"
			error_rel_direct = "''" .. error_rel_direct .. "''"
			bgclr = " style=\"background-color: #dddddd;\""
		else
			bgclr = ""
		end
		
		errlist_direct[i] = {
			bgclr = bgclr,
			ratiocomp = ratiocomp,
			ead = ead,
			erd = erd,
			error_abs_direct = error_abs_direct,
			error_rel_direct = error_rel_direct,
		}
		
		errlist_val[i] = {
			bgclr = bgclr,
			ratiocomp = ratiocomp,
			eav = eav,
			erv = erv,
			error_abs_val = error_abs_val,
			error_rel_val = error_rel_val,
		}
	end
	
	table.sort(errlist_direct, function(a, b)
		return a.ead < b.ead
	end)
	
	for i = 1, #error_list do
		t_body_direct[i] = string.format("|-%s\n| %s || %s || %s",
			errlist_direct[i].bgclr,
			errlist_direct[i].ratiocomp,
			errlist_direct[i].error_abs_direct,
			errlist_direct[i].error_rel_direct)
		t_body_val[i] = string.format("|-%s\n| %s || %s || %s",
			errlist_val[i].bgclr,
			errlist_val[i].ratiocomp,
			errlist_val[i].error_abs_val,
			errlist_val[i].error_rel_val)
	end
	
	local t_head_0 = "{| class=\"wikitable center-all mw-collapsible mw-collapsed sortable\"\n"
		.. "|+ style=\"font-size: 105%; white-space: nowrap;\" | "
	local t_head_1 = "|-\n"
		.. "! class=\"unsortable\" | Interval and complement"
		.. " !! Error (abs, [[Cent|&#162;]])"
		.. " !! Error (rel, [[Relative cent|%]])\n"
	local t_foot = "\n|}"
	
	if tag ~= "none" then
		tag_direct = " (direct approximation, even if inconsistent)"
		tag_val = " (patent val mapping)"
	else
		tag_direct = ""
		tag_val = ""
	end
	
	if approx_t == "direct" then
		result = t_head_0 .. t_title .. tag_direct .. "\n" .. t_head_1 .. table.concat(t_body_direct, "\n") .. t_foot
	elseif approx_t == "val" then
		result = t_head_0 .. t_title .. tag_val .. "\n" .. t_head_1 .. table.concat(t_body_val, "\n") .. t_foot
	else
		result = ""
	end
	
	return result
end

-- local function prec_by_equal (steps)
-- 	return math.floor (math.log (steps*1.9)/math.log (10))
-- end

function p.q_odd_limit_intervals(frame)
	local steps = tonumber(frame.args["steps"])
	local limit = math.max(tonumber(frame.args["limit"]), 2)
	local constcy = limits.consistency_limit(ET.parse(steps .. "edo"), false, 43)
    local debugg = yesno(frame.args["debug"])
	local note = frame.args["note"]
	local title = frame.args["title"]
	local header = string.lower(frame.args["header"]:gsub("%s+", ""))
	local apx = string.lower(frame.args["apx"]:gsub("%s+", ""))
	local tag
	
	if not utils.value_provided(title) then
		title = string.format("%d-odd-limit intervals in %dedo", limit, steps)
	end
	
	local subgroup = table_filter(PRIME_LIST, limit)
	local monzo_list = odd_limit_monzo_list_gen(limit)
	local result
	if header ~= "none" then
		if constcy >= limit then
			apx = "val"
		end
		result = string.format("The following %s how [[%d-odd-limit]] intervals are represented in %dedo. ",
			string.gsub("table show", "()", {[((apx == "direct" or apx == "val") and 11 or 6)] = "s"}), limit, steps)
			.. "Prime harmonics are in '''bold'''"
			.. (constcy >= limit and string.format(".\n\nAs %dedo is consistent in the %d-odd-limit, "
				.. "the mappings by direct approximation and through the patent val are identical.", steps, limit)
				or "; inconsistent intervals are in ''italics''.")
			.. "\n\n"
		if constcy >= limit then
			tag = "none"
		end
		if note ~= "" then
			result = result .. note .. "\n\n"
		end
		if steps == 12 then
			result = result .. "Note that, since the [[cent]] was defined in terms of 12edo, the absolute and relative errors for 12edo are identical.\n\n"
		end
	else
		tag = string.lower(frame.args["tag"]:gsub("%s+", ""))
		result = ""
	end
	
	if utils.value_provided(apx) then
		result = result .. approx(steps, subgroup, monzo_list, title, tag, apx)
	else
		result = result .. approx(steps, subgroup, monzo_list, title, tag, "direct") .. "\n\n" .. approx(steps, subgroup, monzo_list, title, tag, "val")
	end
	
	return frame:preprocess(debugg == true and "<syntaxhighlight lang=\"wikitext\">" .. result .. "</syntaxhighlight>" or result)
end

return p