-- tkz_elements_tkz.lua
-- date 2025/06/10
-- version 4.10c
-- 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.

-- tkz-settings and math

tkz = {}
tkz.nb_dec = 10         
tkz.epsilon = 0.1 ^ tkz.nb_dec 
tkz.dc = 2   
tkz.phi = (1 + math.sqrt(5)) / 2     -- golden number φ
tkz.invphi = (math.sqrt(5) - 1) / 2 function tkz.midpoint(a, b)
tkz.sqrtphi = math.sqrt(tkz.phi)  
  return (a + b) / 2
end

function tkz.reset_defaults()
  tkz.nb_dec  = 10
  tkz.epsilon = 0.1 ^ tkz.nb_dec
  tkz.dc      = 2
end

function tkz.set_nb_dec(n)
  tkz.nb_dec = n
  tkz.epsilon = 0.1 ^ n
end

function tkz.midpoints(...)
  return midpoints_(...)
end

function tkz.length(a, b)
  return point.abs(a - b)
end

function tkz.parabola(a, b, c) --bug local
	local xa, xb, xc, ya, yb, yc
	xa = a.re
	ya = a.im
	xb = b.re
	yb = b.im
	xc = c.re
	yc = c.im
 return parabola_(xa, ya, xb, yb, xc, yc)
end

function tkz.get_angle(a, b, c)
	return angle_normalize_(get_angle_(a, b, c))
end

function tkz.angle_normalize(a)
	return angle_normalize_(a)
end

function tkz.round(num, idp)
  idp = idp or 0
  local mult = 10 ^ idp
	return math.floor(num * mult + 0.5) / mult
end


function tkz.solve_linear_system(M, N)
	if not M or not N then
		return nil, "Error: Input matrices are nil."
	end

	local m, n = M.rows, M.cols

	if N.rows ~= m or N.cols ~= 1 then
		return nil, "Error: Incompatible dimensions between M and N."
	end

	local Aug = {}

	-- Construction of the augmented matrix
	for i = 1, m do
		Aug[i] = {}
		for j = 1, n do
			Aug[i][j] = M.set[i][j]
		end
		Aug[i][n + 1] = N.set[i][1] -- Add column vector
	end

	local MA = matrix:new(Aug)

	local rankM = M:rank()
	local rankMA = MA:rank()

	-- System compatibility check
	if rankM < rankMA then
		return nil, "Incompatible system (no solution)"
	elseif rankM < n then
		return nil, "Indeterminate system (infinite solutions)"
	end

	-- Résolution par Gauss-Jordan
	local sol = MA:gauss_jordan_rect()

	--
	local S = {}
	for i = 1, m do
		S[i] = { sol.set[i][n + 1] } -- n + 1 car n=3, on veut la colonne 4
	end

	return matrix:new(S)
end

function tkz.bisector(a, b, c)
  local i = in_center_(a, b, c)
  return line:new(a, intersection_ll_(a, i, b, c))
end

function tkz.altitude(a, b, c)
  local o, p
  -- Get the orthocenter (which is the point of concurrency of the altitudes)
  o = ortho_center_(a, b, c)

  -- Get the perpendicular projection of point 'a' onto the line defined by 'b' and 'c'
  p = projection_(b, c, a)

  -- Return the altitude, which is the line from point 'a' to the point 'p'
  return line:new(a, p)
end

function tkz.bisector_ext(a, b, c) -- n=1 swap n=2 swap 2
  local i, p
  -- Get the incenter of the triangle
  i = in_center_(a, b, c)

  -- Rotate the incenter by 90 degrees around point 'a' to compute the external bisector
  p = rotation_(a, math.pi / 2, i)

  -- Return the external bisector as the line passing through points 'a' and 'p'
  return line:new(a, p)
end

function tkz.solve(...)
    local params = {...}
    local np = #params 
    if np == 2 then
        local a = params[1]
        local b = params[2]
        return -b / a
    elseif np == 3 then
            return tkz.solve_quadratic_(table.unpack(params))
    elseif np == 4 then
       return tkz.solve_cubic_(table.unpack(params))
    else
        tex.error("solve: Invalid number of parameters (" .. np .. "). Expected 2, 3, or 4.")
    end
end

function tkz.solve_cx_quadratic(a, b, c)
	local d = b * b - 4 * a * c
	local dcx = point.sqrt(d)
	local root1 = (-b + dcx) / (2 * a)
	local root2 = (-b - dcx) / (2 * a)
	return root1, root2
end

function tkz.solve_quadratic_(a, b, c)
	local root1, root2, delta, sqrtdelta
	if (type(a) == "number") and (type(b) == "number") and (type(c) == "number") then
		delta = b * b - 4 * a * c
		if math.abs(delta) < tkz.epsilon then
			delta = 0
		end
		if delta < 0 then
			root1, root2 = false, false --solve_cx_quadratic(a, b, c)
		elseif delta == 0 then
			root1 = -b / (2 * a)
			root2 = -b / (2 * a)
		else
			sqrtdelta = math.sqrt(delta)

			root1 = (-b + sqrtdelta) / (2 * a)
			root2 = (-b - sqrtdelta) / (2 * a)
		end
	else
		root1, root2 = tkz.solve_cx_quadratic(a, b, c)
	end
	return root1, root2 -- Two real roots
end

function tkz.solve_cubic_(a, b, c, d)
    local p = (3 * a * c - b * b) / (3 * a * a)
    local q = (2 * b * b * b - 9 * a * b * c + 27 * a * a * d) / (27 * a * a * a)

    local delta = q * q / 4 + p * p * p / 27
    local epsilon = 1e-10 -- Tolérance pour les comparaisons de nombres flottants

    local function cubert(x)
        return x >= 0 and x ^ (1 / 3) or -(-x) ^ (1 / 3)
    end

    if delta > epsilon then
        local u = cubert((-q + math.sqrt(delta)) / 2)
        local v = cubert((-q - math.sqrt(delta)) / 2)
        local x1 = u + v - b / (3 * a)

        return x1, nil, nil
    elseif math.abs(delta) <= epsilon then
        local u = cubert(-q / 2)
        local x1 = 2 * u - b / (3 * a)
        local x2 = -u - b / (3 * a)

        return x1, x2, x2
    else
        local r = math.sqrt(-p * p * p / 27)
        local t = -q / (2 * r)

        if t > 1 then t = 1 end
        if t < -1 then t = -1 end

        local theta = math.acos(t)
        local x1 = 2 * cubert(r) * math.cos(theta / 3) - b / (3 * a)
        local x2 = 2 * cubert(r) * math.cos((theta + 2 * math.pi) / 3) - b / (3 * a)
        local x3 = 2 * cubert(r) * math.cos((theta + 4 * math.pi) / 3) - b / (3 * a)

        return x1, x2, x3
    end
end

-- straight line from a through b
function tkz.line_coefficients(xa, ya, xb, yb)
	return (yb - ya) / (xb - xa), (xb * ya - xa * yb) / (xb - xa)
end

function tkz.angle_between_vectors(a, b, c, d)
  return tkz_angle_between_vectors_(a, b, c, d)
end

function tkz.barycenter(...)
	return barycenter_(...)
end

function tkz.is_linear(z1, z2, z3)
	return math.abs((z2 - z1) ^ (z3 - z1)) < tkz.epsilon
end

function tkz.is_ortho(z1, z2, z3)
	return math.abs((z2 - z1) .. (z3 - z1)) < tkz.epsilon
end

function tkz.is_zero(x)
	return is_zero_(x)
end

function tkz.set_zero(x)
  return set_zero_(x)
end

function tkz.residue(x)
  return residue_(x)
end

function tkz.is_integer(x)
  return type(x) == "number"
     and x == x               -- exclude NaN
     and x < math.huge        -- exclude +inf
     and -x < math.huge       -- exclude -inf
     and is_zero_(x - tkz_round_(x))
end

function tkz.near_integer(x)
	return is_zero_(x % 1)
end


function tkz.orient2d(x1, y1, x2, y2, x3, y3)
    return (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)
end

function tkz.cramer(...)
  local t = {...}
  local n = #t
  if n == 4 then
    local a1, a2, b1, b2 = table.unpack(t)
    return a1 * b2 - a2 * b1
  elseif n == 9 then
    local a1, a2, a3, b1, b2, b3, c1, c2, c3 = table.unpack(t)
    return a1 * b2 * c3 + a3 * b1 * c2 + a2 * b3 * c1
         - a3 * b2 * c1 - a1 * b3 * c2 - a2 * b1 * c3
  else
    error("cramer: invalid number of arguments (expected 4 or 9)")
  end
end

function tkz.dot_product(a, b, c)
	return (b - a) .. (c - a)
end

function tkz.round(num, idp)
	return tkz_round_(num, idp)
end

function tkz.display(z)
  return tkz_display_(z)
end

function tkz.poncelet_point(a, b, c, d)
  local e1 = euler_center_(a, b, c)
  local e2 = euler_center_(a, c, d)
  local m = midpoint_(a, c)
  local x, y = intersection_cc_(e1, m, e2, m)

  if not x or not y then
    error("Poncelet construction failed: no intersection.")
  end

  if x == m then
    return y
  else
    return x
  end
end

function orthopole(a, b, c, l)
  local ap, bp, cp = l:projection(a, b, c)
  local bpp = self.ca:projection(bp)
  local app = self.bc:projection(ap)
  local la = line(ap, app)
  local lb = line(bp, bpp)
  return intersection(la, lb)
end
return tkz