-- tkz_elements_functions_matrices.lua
-- date 2025/03/04
-- version 3.34c
-- Copyright 2025  Alain Matthes
-- This work may be distributed and/or modified under the
-- conditions of the LaTeX Project Public License, either version 1.3
-- of this license or (at your option) any later version.
-- The latest version of this license is in
-- http://www.latex-project.org/lppl.txt
-- and version 1.3 or later is part of all distributions of LaTeX
-- version 2005/12/01 or later.
-- This work has the LPPL maintenance status “maintained”.
-- The Current Maintainer of this work is Alain Matthes.

-- ----------------------------------------------------------------------------
function print_matrix(m, mstyle)
	local mstyle = (mstyle or "bmatrix")
	local m = (m.type == "matrix" and m.set or m)
	tex.sprint("$")
	tex.sprint("\\begin{" .. mstyle .. "}")
	for i = 1, #m do
		for j = 1, #m[1] do
			local x = m[i][j]
			local st = display(x)
			tex.sprint(st)
			if j < #m[1] then
				tex.sprint(" & ")
			end
		end
		tex.sprint("\\\\")
	end
	tex.sprint("\\end{" .. mstyle .. "}")
	tex.sprint("$")
end

function print_array(matrix)
	local mdata = (matrix.type == "matrix" and matrix.set or matrix)
	tex.sprint("\\{%")
	for i = 1, #mdata do
		local row = mdata[i]
		local row_str = "{\\{"
		for j = 1, #row do
			row_str = row_str .. " " .. tostring(row[j])
			if j < #row then
				row_str = row_str .. "  ,"
			end
		end
		if (i ~= #mdata) and (j ~= #row) then
			if i > 1 then
				row_str = row_str .. " \\}},"
			else
				row_str = row_str .. " \\}},"
			end
			tex.sprint(row_str)
		else
			if i > 1 then
				row_str = row_str .. " \\}}"
			else
				row_str = row_str .. " \\}}"
			end
			tex.sprint(row_str)
		end
	end
	tex.sprint("\\}")
end

function mul_matrix(A, B)
	local adata = (A.type == "matrix" and A.set or A)
	local bdata = (B.type == "matrix" and B.set or B)
	local C = {}
	for i = 1, #adata do
		C[i] = {}
		for j = 1, #bdata[1] do
			local num = adata[i][1] * bdata[1][j]
			for k = 2, #adata[1] do
				num = num + adata[i][k] * bdata[k][j]
			end
			C[i][j] = num
		end
	end
	return matrix:new(C)
end

function add_matrix(A, B)
	local adata = (A.type == "matrix" and A.set or A)
	local bdata = (B.type == "matrix" and B.set or B)
	local S = {}
	for i = 1, #adata do
		local T = {}
		S[i] = T
		for j = 1, #adata[1] do
			T[j] = adata[i][j] + bdata[i][j]
		end
	end
	return matrix:new(S)
end

function k_mul_matrix(n, A)
	local adata = (A.type == "matrix" and A.set or A)
	local S = {}
	for i = 1, #adata, 1 do
		local T = {}
		S[i] = T
		for j = 1, #adata[1], 1 do
			T[j] = n * adata[i][j]
		end
	end
	return matrix:new(S)
end

function transposeMatrix(A)
	local mdata = (A.type == "matrix" and A.set or A)
	local transposedMatrix = {}
	for i = 1, #mdata[1] do
		transposedMatrix[i] = {}
		for j = 1, #mdata do
			transposedMatrix[i][j] = mdata[j][i]
		end
	end
	return matrix:new(transposedMatrix)
end

-- Function to calculate the determinant of a square matrix
function determinant(A)
	local matrix = (A.type == "matrix" and A.set or A)
	if #matrix == #matrix[1] then
		local n = #matrix
		if n == 1 then
			return matrix[1][1] -- Base case for 1x1 matrix
		elseif n == 2 then
			return matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1] -- Base case for 2x2 matrix
		else
			local det = 0
			for j = 1, n do
				local minor = {}
				for i = 2, n do
					minor[i - 1] = {}
					for k = 1, n do
						if k ~= j then
							minor[i - 1][#minor[i - 1] + 1] = matrix[i][k]
						end
					end
				end
				det = det + ((-1) ^ (j + 1)) * matrix[1][j] * determinant(minor) -- Recursive call for larger matrices
			end
			return det
		end
	else
		return nil
	end
end

function check_square_matrix(A)
	local matrix = (A.type == "matrix" and A.set or A)
	if #matrix == #matrix[1] then
		return true
	else
		return false
	end
end

function id_matrix(n)
	local identityMatrix = {}
	for i = 1, n do
		identityMatrix[i] = {}
		for j = 1, n do
			if i == j then
				identityMatrix[i][j] = 1
			else
				identityMatrix[i][j] = 0
			end
		end
	end
	return matrix:new(identityMatrix)
end

function inverse_2x2(A)
	local m = (A.type == "matrix" and A.set or A)
	local a, b, c, d = m[1][1], m[1][2], m[2][1], m[2][2]
	local D = A.det
	if D == 0 then
		return nil -- La matrice n'est pas inversible
	else
		local inv = {}
		inv[1] = {}
		inv[1][1] = d / D
		inv[1][2] = -b / D
		inv[2] = {}
		inv[2][1] = -c / D
		inv[2][2] = a / D
		return matrix:new(inv)
	end
end

function adjugate_(A)
	local m = (A.type == "matrix" and A.set or A)
	if #m == 2 then
		local a, b, c, d = m[2][2], -m[1][2], -m[2][1], m[1][1]
		return matrix:new({ { a, b }, { c, d } })
	elseif #m == 3 then
		local a, b, c = m[1][1], m[1][2], m[1][3]
		local d, e, f = m[2][1], m[2][2], m[2][3]
		local g, h, i = m[3][1], m[3][2], m[3][3]
		return transposeMatrix(matrix:new({
			{ e * i - f * h, -(d * i - f * g), d * h - e * g },
			{ -(b * i - c * h), a * i - c * g, -(a * h - b * g) },
			{ b * f - c * e, -(a * f - c * d), a * e - b * d },
		}))
	else
		return nil
	end
end

function inverse_3x3(A)
	local D = A.det
	if D == 0 then
		return nil -- La matrice n'est pas inversible
	else
		local adj = adjugate_(A)
		local m = (adj.type == "matrix" and adj.set or adj)
		local inv = {}
		for i = 1, 3 do
			inv[i] = {}
			for j = 1, 3 do
				inv[i][j] = m[i][j] / D
			end
		end
		return matrix:new(inv)
	end
end

-- inverse only for 2x2 or 3x3 matrix
function inv_matrix(A)
	if A.det == 0 then
		tex.print("Matrix not inversible: det = 0")
		return nil
	else
		local M = (A.type == "matrix" and A.set or A)
		local n = #M
		if n == 2 then
			local m = (A.type == "matrix" and A.set or A)
			local a, b, c, d = m[1][1], m[1][2], m[2][1], m[2][2]
			local D = A.det
			local inv = {}
			inv[1] = {}
			inv[1][1] = d / D
			inv[1][2] = -b / D
			inv[2] = {}
			inv[2][1] = -c / D
			inv[2][2] = a / D
			return matrix:new(inv)
		else
			local D = A.det
			local adj = adjugate_(A)
			local m = (adj.type == "matrix" and adj.set or adj)
			local inv = {}
			for i = 1, 3 do
				inv[i] = {}
				for j = 1, 3 do
					inv[i][j] = m[i][j] / D
				end
			end
			return matrix:new(inv)
		end
	end
end

function diagonalize_(A)
	local m = (A.type == "matrix" and A.set or A)
	local trace = m[1][1] + m[2][2]
	local a, b = m[1][1], m[1][2]
	local det = A.det
	local D = trace * trace - 4 * det
	if D > 0 then
		local va1 = (trace + math.sqrt(D)) / 2
		local va2 = (trace - math.sqrt(D)) / 2
		return matrix:new({ { va1, 0 }, { 0, va2 } }), matrix:new({ { 1, 1 }, { (va1 - a) / b, (va2 - a) / b } })
	else
		local va1 = point(trace / 2, math.sqrt(-D) / 2)
		local va2 = point(trace / 2, -math.sqrt(-D) / 2)
		return matrix:new({ { va1, 0 }, { 0, va2 } }) --,
		--   matrix : new ({{1,1},{ (va1 - a )/b, (va2 - a)/b}})
	end
end

function isDiagonal_(A)
	local matrix = (A.type == "matrix" and A.set or A)
	if check_square_matrix(A) == true then
		for i = 1, #matrix do
			for j = 1, #matrix[1] do
				if i ~= j and matrix[i][j] ~= 0 then
					return false
				end
			end
		end
		return true
	else
		return false
	end
end

function isOrthogonal_(A)
	local m = (A.type == "matrix" and A.set or A)
	if (check_square_matrix(A) == true) and (A.det ~= 0) then
		local mT = transposeMatrix(A)
		local mI = inv_matrix(A)
		if mT == mI then
			return true
		else
			return false
		end
	else
		return false
	end
end

function homogenization_(A)
	local m = (A.type == "matrix" and A.set or A)
	if A.cols ~= 1 then
		return nil
	else
		local a, b, c
		a = m[1][1]
		b = m[2][1]
		c = 1
		return matrix:new({ { a }, { b }, { c } })
	end
end

function get_element_(A, i, j)
	local m = (A.type == "matrix" and A.set or A)
	if m[i] and m[i][j] then
		return m[i][j]
	end
end

function get_htm_point(A)
	local m = (A.type == "matrix" and A.set or A)
	return point:new(m[1][1], m[2][1])
end

function htm_apply_(A, z)
	local V = homogenization_(z.mtx)
	local W = A * V
	return get_htm_point(W)
end

function htm_apply_L_(A, obj)
	local x, y
	x = htm_apply_(A, obj.pa)
	y = htm_apply_(A, obj.pb)
	return line:new(x, y)
end

function htm_apply_C_(A, obj)
	local x, y
	x = htm_apply_(A, obj.center)
	y = htm_apply_(A, obj.through)
	return circle:new(x, y)
end

function htm_apply_T_(A, obj)
	local x, y, z
	x = htm_apply_(A, obj.pa)
	y = htm_apply_(A, obj.pb)
	z = htm_apply_(A, obj.pc)
	return triangle:new(x, y, z)
end

function htm_apply_Q_(A, obj)
	local x, y, z, t
	x = htm_apply_(A, obj.pa)
	y = htm_apply_(A, obj.pb)
	z = htm_apply_(A, obj.pc)
	t = htm_apply_(A, obj.pd)
	if obj.type == "square" then
		return square:new(x, y, z, t)
	elseif obj.type == "rectangle" then
		return rectangle:new(x, y, z, t)
	elseif obj.type == "parallelogram" then
		return parallelogram:new(x, y, z, t)
	elseif obj.type == "quadrilateral" then
		return quadrilateral:new(x, y, z, t)
	end
end