View source for Module:Temperament data

local p = {}

local rat = require("Module:Rational")
local u = require("Module:Utils")

local function gcd(a, b)
    if type(a) == "number" and type(b) == "number" and a == math.floor(a) and b == math.floor(b) then
        if b == 0 then
            return a
        else
            return gcd(b, a % b) -- tail recursion
        end
    else
        error("Invalid argument to gcd (" .. tostring(a) .. ", " .. tostring(b) .. ")", 2)
    end
end

-- Linear algebra and RTT functions

local function matadd(a, b)
	local result = {}
	for i = 1, #a  do
		result[i] = {}
		for j = 1, #(b[1]) do
			result[i][j] = a[i][j] + b[i][j]
		end
	end
	return result
end

local function matsub(a, b)
	local result = {}
	for i = 1, #a  do
		result[i] = {}
		for j = 1, #(b[1]) do
			result[i][j] = a[i][j] - b[i][j]
		end
	end
	return result
end

local function matmul(a, b)
	local result = {}
	for i = 1, #a  do
		result[i] = {}
		for j = 1, #(b[1]) do
			result[i][j] = 0
			for k = 1, #(a[1]) do
				result[i][j] = result[i][j] + (a[i][k] * b[k][j])
			end
		end
	end
	return result
end


local function scalarmatmul(a, b)
	local result = {}
	for i = 1, #a  do
		result[i] = {}
		for j = 1, #(a[1]) do
			result[i][j] = a[i][j] * b
		end
	end
	return result
end

local function matinv(a)
	local xn = scalarmatmul(a, 1e-7)

	for i = 1, 75 do
		xn = matsub(scalarmatmul(xn, 2), matmul(xn, matmul(a, xn)))
	end
	
	return xn
end

local function transpose(a)
	local result = {}
	for i = 1, #a[1] do
		result[i] = {}
		for j = 1, #a do
			result[i][j] = a[j][i]
		end
	end
	return result
end

local function antitranspose(a)
	local result = {}
	for i = 1, #a[1] do
		result[i] = {}
		for j = 1, #a do
			result[i][j] = a[#a - j + 1][#a[1] - i + 1]
		end
	end
	return result
end


local function pseudoinv(a)
	return matmul(transpose(a), matinv(matmul(a, transpose(a))))
end

local function nullspace(mapping)
	local identity = {}
	for i = 1, #mapping[1] do
		identity[i] = {}
		for j = 1, #mapping[1] do
			if i == j then
				identity[i][j] = 1
			else
				identity[i][j] = 0
			end
		end
	end

	-- local w = {{0},{1},{0}}
	-- for i = 1, #mapping[1] do
	-- 	w[i] = {10}
	-- end

	return matsub(identity, matmul(pseudoinv(mapping), mapping))
end

local function unreduced_mapping_from_basis(comma_basis)
	return antitranspose(nullspace(antitranspose(comma_basis)))
end

local function get_reduced_mapping(comma_basis, preimage)
	a =  pseudoinv(matmul(unreduced_mapping_from_basis(comma_basis), preimage))
	for i = 1, #a do
		for j = 1, #a[1] do
			a[i][j] = math.floor(a[i][j] + 0.5) -- round each entry (they are not exact integers)
		end
	end
	return a
end

local function get_te_generator(subgroup, comma_basis, preimage)
	local v = get_reduced_mapping(comma_basis, preimage)
	local w = {}
	for i = 1, #subgroup do
		w[i] = {}
		for j = 1, #subgroup do
			if i == j then
				w[i][j] = math.log(2) / math.log(subgroup[i])
			else
				w[i][j] = 0
			end
		end
	end
	
	local jw = {{}}
	for i = 1, #subgroup do
		jw[1][i] = 1
	end
	local vw = matmul(v, w)
	local g = matmul(jw, pseudoinv(vw))
	return g
end


local function get_pote_generator(subgroup, comma_basis, preimage)
	local period = 1
	for i = 1, #subgroup do
		period = period * (subgroup[i]^preimage[i][1])
	end
	local te = get_te_generator(subgroup, comma_basis, preimage)
	local stretch_factor = te[1][1] * math.log(2) / math.log(period)
	return scalarmatmul(te, 1 / stretch_factor)
end



-- Parsing/display functions

local function int_to_subgroup_monzo(subgroup, x) 
	local result = {}
	local x2 = x
	for i = 1, #subgroup do
		result[i] = 0
		while true do
			x2 = x2 / subgroup[i]
			if x2 ~= math.floor(x2) then
				break
			end
			result[i] = result[i] + 1
		end
		x2 = x
	end
	
	return result
end


local function rat_to_subgroup_monzo(subgroup, x)
	local n, d = rat.as_pair(x)
	return matsub({int_to_subgroup_monzo(subgroup, n)}, {int_to_subgroup_monzo(subgroup, d)})[1]	
end

local function rat_list_to_matrix(subgroup, list)
	local result = {}
	for j = 1, #subgroup do
		result[j] = {}
	end
	
	for i = 1, #list do
		local smonzo = rat_to_subgroup_monzo(subgroup, list[i])
		for j = 1, #subgroup do
			result[j][i] = smonzo[j]
		end
	end
	
	return result
end

local function mysplit (inputstr, sep)
        if sep == nil then
    		sep = "%s"
        end
        local t = {}
        for str in string.gmatch(inputstr, "([^" .. sep .. "]+)") do
                table.insert(t, str)
        end
        return t
end

local function trim(x)
	local str = x
	str = str:gsub("%s+", "")
    str = string.gsub(str, "%s+", "")
    return str
end

function p.temperament_data(frame)
	local subgroup = mysplit(frame.args["subgroup"], ".")
	for i = 1, #subgroup do
		subgroup[i] = tonumber(subgroup[i])
	end
	local comma_matrix = mysplit(frame.args["comma_list"], ",")
	for i = 1, #comma_matrix do
		comma_matrix[i] = rat.parse(comma_matrix[i])
	end
	comma_matrix = rat_list_to_matrix(subgroup, comma_matrix)
	local unparsed_gens = mysplit(frame.args["generators"], ",")
	local generators =  mysplit(frame.args["generators"], ",")
	for i = 1, #generators do
		generators[i] = rat.parse(generators[i])
	end
	generators = rat_list_to_matrix(subgroup, generators)
	local mapping = get_reduced_mapping(comma_matrix, generators)
	local cte_generator = frame.args["cte_generator"]
	local pote_generator = get_pote_generator(subgroup, comma_matrix, generators)
	local result =  "[[Subgroup]]: " .. frame.args["subgroup"]
	result = result .. "\n\n[[Comma list]]: " .. frame.args["comma_list"]
	result = result .. "\n\n[[Mapping]]: [⟨"
	for i = 1, #mapping do
		for j = 1, #mapping[1] do
			result = result .. mapping[i][j] .. " "
		end
		result = result:sub(0,-2) .. "], ⟨"
	end	
	result = result:sub(0,-6) .. "]"
	
	result = result .. "\n\n: mapping generators: "
	for i = 1, #unparsed_gens do
		result = result .. "~" .. trim(unparsed_gens[i]) .. ", "
	end
	result = result:sub(0, -3)
	if cte_generator ~= "" then
		cte_generator = mysplit(cte_generator, ",")
		for i = 1, #cte_generator do
			cte_generator[i] = tonumber(cte_generator[i])
		end
		
		result = result .. "\n\n[[Optimal tuning]]s:\n* [[CTE]]:"
		
		for i = 1, #(cte_generator) do
			result = result .. "~" .. trim(unparsed_gens[i]) .. " = "
			if subgroup[1] == 2 and i == 1 then
				result = result .. "1\\1"
			elseif subgroup[1] == 3 and i == 1 then
				result = result .. "1\\1edt"
			else
				result = result .. u._round(cte_generator[i], 7)
			end
			result = result .. ", "
		end
		
		result = result:sub(0,-3)
	
		if subgroup[1] ~= 2 then
			result = result .. "\n* [[Lp tuning|POL2]]:"
		else
			result = result .. "\n* [[POTE]]:"
		end
	else
		if subgroup[1] ~= 2 then
			result = result .. "\n\n[[Optimal tuning]] ([[Lp tuning|POL2]]): "
		else
			result = result .. "\n\n[[Optimal tuning]] ([[POTE]]): "
		end
	end
	
	for i = 1, #(pote_generator[1]) do
		result = result .. "~" .. trim(unparsed_gens[i]) .. " = "
		if subgroup[1] == 2 and i == 1 then
			result = result .. "1\\1"
		elseif subgroup[1] == 3 and i == 1 then
			result = result .. "1\\1edt"
			else
			result = result .. u._round(pote_generator[1][i] * 1200, 7)
		end
		result = result .. ", "
	end
	result = result:sub(0,-3)
	
	return result
end

return p