-- tkz_elements_functions_maths.lua
-- date 2025/05/25
-- version 4.00c
-- 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 tkz_round_str(num, idp)
	return topoint(string.format("%." .. (idp or 0) .. "f", num))
end

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


function format_coord(x, decimals)
  decimals = decimals or 5
  local fmt = "%." .. tostring(decimals) .. "f"
  return string.format(fmt, checknumber(x))
end

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

function Cramer33(a1, a2, a3, b1, b2, b3, c1, c2, c3)
	return a1 * b2 * c3 + a3 * b1 * c2 + a2 * b3 * c1 - a3 * b2 * c1 - a1 * b3 * c2 - a2 * b1 * c3
end

function Cramer22(a1, a2, b1, b2)
	return a1 * b2 - a2 * b1
end

function det_(x1, y1, x2, y2, x3, y3)
    return (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)
end

function aligned(m, a, b)
	local z
	z = (b - a) / (m - b)
	if math.abs(z.im) < tkz_epsilon then
		return true
	else
		return false
	end
end

function islinear(z1, z2, z3)
	local dp
	dp = (z2 - z1) ^ (z3 - z1)
	if math.abs(dp) < tkz_epsilon then
		return true
	else
		return false
	end
end
is_linear = islinear

function isortho(z1, z2, z3) --modif
	local dp = (z2 - z1) .. (z3 - z1)
	if math.abs(dp) < tkz_epsilon then
		return true
	else
		return false
	end
end

is_ortho = isortho

function 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
	local D = (xa - xb) * (xa - xc) * (xb - xc)
	local A = (xc * (yb - ya) + xb * (ya - yc) + xa * (yc - yb)) / D
	local B = (xc * xc * (ya - yb) + xb * xb * (yc - ya) + xa * xa * (yb - yc)) / D
	local C = (xb * xc * (xb - xc) * ya + xc * xa * (xc - xa) * yb + xa * xb * (xa - xb) * yc) / D
	return A, B, C
end

function parabola_(xa, ya, xb, yb, xc, yc) -- added
	local D = (xa - xb) * (xa - xc) * (xb - xc)
	local A = (xc * (yb - ya) + xb * (ya - yc) + xa * (yc - yb)) / D
	local B = (xc * xc * (ya - yb) + xb * xb * (yc - ya) + xa * xa * (yb - yc)) / D
	local C = (xb * xc * (xb - xc) * ya + xc * xa * (xc - xa) * yb + xa * xb * (xa - xb) * yc) / D
	return A, B, C
end

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

function get_angle_(a, b, c)
	if (b == a) or (c == a) then
		tex.error("Points confondus dans get_angle_")
	end
	return point.arg((c - a) / (b - a))
end

function angle_normalize(a)
	return angle_normalize_(a)
end

function angle_normalize_(a)
	while a < 0 do
		a = a + 2 * math.pi
	end
	while a >= 2 * math.pi do
		a = a - 2 * math.pi
	end
	return a
end

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

-- real func
function is_integer(x)
  if type(x) ~= "number" then
    return false -- Not a number, so not an integer
  end
  if not math.isfinite(x) then
    return false -- NaN or infinity, not an integer
  end
  local rounded = tkz_round(x)
  return math.abs(x - rounded) < tkz_epsilon -- Use a small tolerance
end

function near_integer(x)
	local _, r = math.modf(x)
	if is_zero(r) then
		return true
	end
	return false
end

function residue(x)
  local _, frac = math.modf(x)
  return frac
end

function is_zero(x)
	return math.abs(x) < tkz_epsilon
end

function set_zero(x)
  return is_zero(x) and 0 or x
end




function checknumber(x)
  if type(x) == "number" then
    return string.format("%.12f", x)
  elseif type(x) == "string" and tonumber(x) then
    return string.format("%.12f", tonumber(x))
  else
    return x -- table ou invalide
  end
end

local safety_margin = 2
function format_number(number, dcpl)
  if type(number) == "table" then
    return number
  else
    local num_dcpl
    if dcpl ~= nil then
      num_dcpl = dcpl
    else
      num_dcpl = tkz_nb_dec + safety_margin -- Default value 
    end

    local format_string = string.format("%%.%df", num_dcpl)
    return string.format(format_string, number)
  end
end

function get_sign(number)
	local sgn
	if math.abs(number) < tkz_epsilon then
		sgn = ""
	elseif number > 0 then
		sgn = "+"
	else
		sgn = "-"
	end
	return sgn
end

function 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 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 solve_quadratic(table.unpack(params))
    elseif np == 4 then
       return solve_cubic(table.unpack(params))
    else
        tex.error("solve: Invalid number of parameters (" .. np .. "). Expected 2, 3, or 4.")
    end
end

function 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 = solve_cx_quadratic(a, b, c)
	end
	return root1, root2 -- Two real roots
end

function 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

function display_real(r)
	local format_string, format_string
	if r == nil then
		return ""
	else
		if near_integer(r) then
			r = math.round(r)
			format_string = string.format("%%.%df", 0)
		else
			format_string = string.format("%%.%df", tkz_dc)
		end
		local st = string.format(format_string, r)
		return st
	end
end

function display_imag(r)
	local sgn
	sgn = get_sign(r)
	r = math.abs(r)
	if math.abs(r - 1) < tkz_epsilon then
		r = nil
	elseif near_integer(r) then
		r = math.abs(math.round(r))
	end
	st = display_real(r)
	return sgn, st
end

function display(z)
	if type(z) == "number" then
		return display_real(z)
	else
		local real, imag
		real = z.re
		imag = z.im
		if is_zero(imag) then
			return display_real(real)
		else
			str = display_real(real)
			sgni, sti = display_imag(imag)
			if str == "0" then
				str = ""
				sgni = ""
			end
			local st = str .. sgni .. sti .. "i"
			return st
		end
	end
end

function get_v(L, p) --bug
	return report_(L.pa, L.pb, 1, p)
end

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

-- intersection line / parabola  good OCCS !! p paramètre de la parabole
--  y= x^{2} /2p  with occs S,i,j j same direction KF
-- m and n parameters of the line

function solve_para_line(p, m, n)
	return solve_quadratic(1, -2 * p * m, -2 * p * n)
end

function solve_hyper_line(a, b, m, p)
	local A = a ^ 2 / b ^ 2 - m ^ 2
	local B = -2 * m * p
	local C = a ^ 2 - p ^ 2
	return solve_quadratic(A, B, C)
end


function angle_between_vectors(a, b, c, d)
	-- Calcul des vecteurs
	local zab = b - a
	local zcd = d - c

	-- Angle entre les vecteurs en utilisant argument du rapport
	local theta = math.atan2(zab.im * zcd.re - zab.re * zcd.im, zab.re * zcd.re + zab.im * zcd.im)

	return theta -- L'angle en radians
end

table.getn = function(t) return #t end

function get_nb_elements(t)
  return table.getn(t)
end

function get_elements(t)
  return table.unpack(t)
end

function reverse(t)
  local n = #t
  for i = 1, math.floor(n / 2) do
    t[i], t[n - i + 1] = t[n - i + 1], t[i]
  end
  return t
end