-module(viva_math@matrix).
-compile([no_auto_import, nowarn_unused_vars, nowarn_unused_function, nowarn_nomatch, inline]).
-define(FILEPATH, "src/viva_math/matrix.gleam").
-export([mat2_zero/0, mat2_identity/0, mat2_add/2, mat2_scale/2, mat2_mul/2, mat2_transpose/1, mat2_determinant/1, mat2_trace/1, mat2_inverse/1, mat2_rotation/1, mat3_zero/0, mat3_identity/0, mat3_add/2, mat3_scale/2, mat3_mul/2, mat3_transpose/1, mat3_determinant/1, mat3_trace/1, mat3_frobenius/1, mat3_inverse/1, mat3_mul_vec3/2, mat3_rot_x/1, mat3_rot_y/1, mat3_rot_z/1, mat2_eigenvalues/1, mat3_symmetric_eigenvalues/1, mat3_diagonal/3, mat4_identity/0, mat4_zero/0, mat4_translation/3, mat4_scale/3, matn_from_rows/1, matn_zeros/2, matn_identity/1, matn_transpose/1, matn_mul/2, matn_add/2, matn_scale/2, matn_trace/1, matn_frobenius/1, matn_mul_vec/2]).
-export_type([mat2/0, mat3/0, mat4/0, mat_n/0]).

-if(?OTP_RELEASE >= 27).
-define(MODULEDOC(Str), -moduledoc(Str)).
-define(DOC(Str), -doc(Str)).
-else.
-define(MODULEDOC(Str), -compile([])).
-define(DOC(Str), -compile([])).
-endif.

?MODULEDOC(
    " Small dense matrices (2×2, 3×3, 4×4) and a generic `MatN`.\n"
    "\n"
    " **Scope**: geometry, dynamics, small-system linear algebra. For batched\n"
    " or large-scale work (GEMM, BLAS, eigendecomposition) use `viva_tensor`,\n"
    " which dispatches to MKL/CUDA.\n"
    "\n"
    " ## Conventions\n"
    "\n"
    " - Row-major storage in the generic `MatN` (rows × cols).\n"
    " - Fixed-size matrices use named fields `m11..mNN`.\n"
    " - All constructors clamp to safe values where appropriate.\n"
).

-type mat2() :: {mat2, float(), float(), float(), float()}.

-type mat3() :: {mat3,
        float(),
        float(),
        float(),
        float(),
        float(),
        float(),
        float(),
        float(),
        float()}.

-type mat4() :: {mat4, list(list(float()))}.

-type mat_n() :: {mat_n, integer(), integer(), list(list(float()))}.

-file("src/viva_math/matrix.gleam", 25).
-spec mat2_zero() -> mat2().
mat2_zero() ->
    {mat2, +0.0, +0.0, +0.0, +0.0}.

-file("src/viva_math/matrix.gleam", 29).
-spec mat2_identity() -> mat2().
mat2_identity() ->
    {mat2, 1.0, +0.0, +0.0, 1.0}.

-file("src/viva_math/matrix.gleam", 33).
-spec mat2_add(mat2(), mat2()) -> mat2().
mat2_add(A, B) ->
    {mat2,
        erlang:element(2, A) + erlang:element(2, B),
        erlang:element(3, A) + erlang:element(3, B),
        erlang:element(4, A) + erlang:element(4, B),
        erlang:element(5, A) + erlang:element(5, B)}.

-file("src/viva_math/matrix.gleam", 37).
-spec mat2_scale(mat2(), float()) -> mat2().
mat2_scale(M, S) ->
    {mat2,
        erlang:element(2, M) * S,
        erlang:element(3, M) * S,
        erlang:element(4, M) * S,
        erlang:element(5, M) * S}.

-file("src/viva_math/matrix.gleam", 41).
-spec mat2_mul(mat2(), mat2()) -> mat2().
mat2_mul(A, B) ->
    {mat2,
        (erlang:element(2, A) * erlang:element(2, B)) + (erlang:element(3, A) * erlang:element(
            4,
            B
        )),
        (erlang:element(2, A) * erlang:element(3, B)) + (erlang:element(3, A) * erlang:element(
            5,
            B
        )),
        (erlang:element(4, A) * erlang:element(2, B)) + (erlang:element(5, A) * erlang:element(
            4,
            B
        )),
        (erlang:element(4, A) * erlang:element(3, B)) + (erlang:element(5, A) * erlang:element(
            5,
            B
        ))}.

-file("src/viva_math/matrix.gleam", 50).
-spec mat2_transpose(mat2()) -> mat2().
mat2_transpose(M) ->
    {mat2,
        erlang:element(2, M),
        erlang:element(4, M),
        erlang:element(3, M),
        erlang:element(5, M)}.

-file("src/viva_math/matrix.gleam", 55).
?DOC(" det(M) = m₁₁m₂₂ - m₁₂m₂₁.\n").
-spec mat2_determinant(mat2()) -> float().
mat2_determinant(M) ->
    (erlang:element(2, M) * erlang:element(5, M)) - (erlang:element(3, M) * erlang:element(
        4,
        M
    )).

-file("src/viva_math/matrix.gleam", 60).
?DOC(" Trace tr(M) = Σ mᵢᵢ.\n").
-spec mat2_trace(mat2()) -> float().
mat2_trace(M) ->
    erlang:element(2, M) + erlang:element(5, M).

-file("src/viva_math/matrix.gleam", 65).
?DOC(" Inverse of a 2×2 matrix. Errors if singular.\n").
-spec mat2_inverse(mat2()) -> {ok, mat2()} | {error, nil}.
mat2_inverse(M) ->
    D = mat2_determinant(M),
    case D =:= +0.0 of
        true ->
            {error, nil};

        false ->
            Inv_d = case D of
                +0.0 -> +0.0;
                -0.0 -> -0.0;
                Gleam@denominator -> 1.0 / Gleam@denominator
            end,
            {ok,
                {mat2,
                    erlang:element(5, M) * Inv_d,
                    +0.0 - (erlang:element(3, M) * Inv_d),
                    +0.0 - (erlang:element(4, M) * Inv_d),
                    erlang:element(2, M) * Inv_d}}
    end.

-file("src/viva_math/matrix.gleam", 82).
?DOC(" 2-D rotation matrix by `theta` radians.\n").
-spec mat2_rotation(float()) -> mat2().
mat2_rotation(Theta) ->
    C = math:cos(Theta),
    S = math:sin(Theta),
    {mat2, C, +0.0 - S, S, C}.

-file("src/viva_math/matrix.gleam", 106).
-spec mat3_zero() -> mat3().
mat3_zero() ->
    {mat3, +0.0, +0.0, +0.0, +0.0, +0.0, +0.0, +0.0, +0.0, +0.0}.

-file("src/viva_math/matrix.gleam", 110).
-spec mat3_identity() -> mat3().
mat3_identity() ->
    {mat3, 1.0, +0.0, +0.0, +0.0, 1.0, +0.0, +0.0, +0.0, 1.0}.

-file("src/viva_math/matrix.gleam", 114).
-spec mat3_add(mat3(), mat3()) -> mat3().
mat3_add(A, B) ->
    {mat3,
        erlang:element(2, A) + erlang:element(2, B),
        erlang:element(3, A) + erlang:element(3, B),
        erlang:element(4, A) + erlang:element(4, B),
        erlang:element(5, A) + erlang:element(5, B),
        erlang:element(6, A) + erlang:element(6, B),
        erlang:element(7, A) + erlang:element(7, B),
        erlang:element(8, A) + erlang:element(8, B),
        erlang:element(9, A) + erlang:element(9, B),
        erlang:element(10, A) + erlang:element(10, B)}.

-file("src/viva_math/matrix.gleam", 128).
-spec mat3_scale(mat3(), float()) -> mat3().
mat3_scale(M, S) ->
    {mat3,
        erlang:element(2, M) * S,
        erlang:element(3, M) * S,
        erlang:element(4, M) * S,
        erlang:element(5, M) * S,
        erlang:element(6, M) * S,
        erlang:element(7, M) * S,
        erlang:element(8, M) * S,
        erlang:element(9, M) * S,
        erlang:element(10, M) * S}.

-file("src/viva_math/matrix.gleam", 142).
-spec mat3_mul(mat3(), mat3()) -> mat3().
mat3_mul(A, B) ->
    {mat3,
        ((erlang:element(2, A) * erlang:element(2, B)) + (erlang:element(3, A) * erlang:element(
            5,
            B
        )))
        + (erlang:element(4, A) * erlang:element(8, B)),
        ((erlang:element(2, A) * erlang:element(3, B)) + (erlang:element(3, A) * erlang:element(
            6,
            B
        )))
        + (erlang:element(4, A) * erlang:element(9, B)),
        ((erlang:element(2, A) * erlang:element(4, B)) + (erlang:element(3, A) * erlang:element(
            7,
            B
        )))
        + (erlang:element(4, A) * erlang:element(10, B)),
        ((erlang:element(5, A) * erlang:element(2, B)) + (erlang:element(6, A) * erlang:element(
            5,
            B
        )))
        + (erlang:element(7, A) * erlang:element(8, B)),
        ((erlang:element(5, A) * erlang:element(3, B)) + (erlang:element(6, A) * erlang:element(
            6,
            B
        )))
        + (erlang:element(7, A) * erlang:element(9, B)),
        ((erlang:element(5, A) * erlang:element(4, B)) + (erlang:element(6, A) * erlang:element(
            7,
            B
        )))
        + (erlang:element(7, A) * erlang:element(10, B)),
        ((erlang:element(8, A) * erlang:element(2, B)) + (erlang:element(9, A) * erlang:element(
            5,
            B
        )))
        + (erlang:element(10, A) * erlang:element(8, B)),
        ((erlang:element(8, A) * erlang:element(3, B)) + (erlang:element(9, A) * erlang:element(
            6,
            B
        )))
        + (erlang:element(10, A) * erlang:element(9, B)),
        ((erlang:element(8, A) * erlang:element(4, B)) + (erlang:element(9, A) * erlang:element(
            7,
            B
        )))
        + (erlang:element(10, A) * erlang:element(10, B))}.

-file("src/viva_math/matrix.gleam", 156).
-spec mat3_transpose(mat3()) -> mat3().
mat3_transpose(M) ->
    {mat3,
        erlang:element(2, M),
        erlang:element(5, M),
        erlang:element(8, M),
        erlang:element(3, M),
        erlang:element(6, M),
        erlang:element(9, M),
        erlang:element(4, M),
        erlang:element(7, M),
        erlang:element(10, M)}.

-file("src/viva_math/matrix.gleam", 161).
?DOC(" Determinant via cofactor expansion along the first row.\n").
-spec mat3_determinant(mat3()) -> float().
mat3_determinant(M) ->
    ((erlang:element(2, M) * ((erlang:element(6, M) * erlang:element(10, M)) - (erlang:element(
        7,
        M
    )
    * erlang:element(9, M))))
    - (erlang:element(3, M) * ((erlang:element(5, M) * erlang:element(10, M)) - (erlang:element(
        7,
        M
    )
    * erlang:element(8, M)))))
    + (erlang:element(4, M) * ((erlang:element(5, M) * erlang:element(9, M)) - (erlang:element(
        6,
        M
    )
    * erlang:element(8, M)))).

-file("src/viva_math/matrix.gleam", 170).
-spec mat3_trace(mat3()) -> float().
mat3_trace(M) ->
    (erlang:element(2, M) + erlang:element(6, M)) + erlang:element(10, M).

-file("src/viva_math/matrix.gleam", 239).
?DOC(" Frobenius norm √(Σ aᵢⱼ²) of a Mat3.\n").
-spec mat3_frobenius(mat3()) -> float().
mat3_frobenius(M) ->
    math:sqrt(
        ((((((((erlang:element(2, M) * erlang:element(2, M)) + (erlang:element(
            3,
            M
        )
        * erlang:element(3, M)))
        + (erlang:element(4, M) * erlang:element(4, M)))
        + (erlang:element(5, M) * erlang:element(5, M)))
        + (erlang:element(6, M) * erlang:element(6, M)))
        + (erlang:element(7, M) * erlang:element(7, M)))
        + (erlang:element(8, M) * erlang:element(8, M)))
        + (erlang:element(9, M) * erlang:element(9, M)))
        + (erlang:element(10, M) * erlang:element(10, M))
    ).

-file("src/viva_math/matrix.gleam", 180).
?DOC(
    " Inverse of a 3×3 matrix via adjugate / determinant.\n"
    "\n"
    " Errors when the matrix is singular or ill-conditioned: the threshold is\n"
    " `|det| < ε · ‖M‖_F³` where ε ≈ machine epsilon. Pure `det == 0.0` is too\n"
    " permissive — matrices with `det = 1e-20` would still pass and produce\n"
    " catastrophic blow-up after division.\n"
).
-spec mat3_inverse(mat3()) -> {ok, mat3()} | {error, nil}.
mat3_inverse(M) ->
    Det = mat3_determinant(M),
    Frob = mat3_frobenius(M),
    Tolerance = ((2.22e-16 * Frob) * Frob) * Frob,
    Det_abs = case Det < +0.0 of
        true ->
            +0.0 - Det;

        false ->
            Det
    end,
    case Det_abs =< Tolerance of
        true ->
            {error, nil};

        false ->
            Inv_det = case Det of
                +0.0 -> +0.0;
                -0.0 -> -0.0;
                Gleam@denominator -> 1.0 / Gleam@denominator
            end,
            {ok,
                {mat3,
                    ((erlang:element(6, M) * erlang:element(10, M)) - (erlang:element(
                        7,
                        M
                    )
                    * erlang:element(9, M)))
                    * Inv_det,
                    ((erlang:element(4, M) * erlang:element(9, M)) - (erlang:element(
                        3,
                        M
                    )
                    * erlang:element(10, M)))
                    * Inv_det,
                    ((erlang:element(3, M) * erlang:element(7, M)) - (erlang:element(
                        4,
                        M
                    )
                    * erlang:element(6, M)))
                    * Inv_det,
                    ((erlang:element(7, M) * erlang:element(8, M)) - (erlang:element(
                        5,
                        M
                    )
                    * erlang:element(10, M)))
                    * Inv_det,
                    ((erlang:element(2, M) * erlang:element(10, M)) - (erlang:element(
                        4,
                        M
                    )
                    * erlang:element(8, M)))
                    * Inv_det,
                    ((erlang:element(4, M) * erlang:element(5, M)) - (erlang:element(
                        2,
                        M
                    )
                    * erlang:element(7, M)))
                    * Inv_det,
                    ((erlang:element(5, M) * erlang:element(9, M)) - (erlang:element(
                        6,
                        M
                    )
                    * erlang:element(8, M)))
                    * Inv_det,
                    ((erlang:element(3, M) * erlang:element(8, M)) - (erlang:element(
                        2,
                        M
                    )
                    * erlang:element(9, M)))
                    * Inv_det,
                    ((erlang:element(2, M) * erlang:element(6, M)) - (erlang:element(
                        3,
                        M
                    )
                    * erlang:element(5, M)))
                    * Inv_det}}
    end.

-file("src/viva_math/matrix.gleam", 209).
?DOC(" Matrix × Vec3 product.\n").
-spec mat3_mul_vec3(mat3(), viva_math@vector:vec3()) -> viva_math@vector:vec3().
mat3_mul_vec3(M, V) ->
    {vec3,
        ((erlang:element(2, M) * erlang:element(2, V)) + (erlang:element(3, M) * erlang:element(
            3,
            V
        )))
        + (erlang:element(4, M) * erlang:element(4, V)),
        ((erlang:element(5, M) * erlang:element(2, V)) + (erlang:element(6, M) * erlang:element(
            3,
            V
        )))
        + (erlang:element(7, M) * erlang:element(4, V)),
        ((erlang:element(8, M) * erlang:element(2, V)) + (erlang:element(9, M) * erlang:element(
            3,
            V
        )))
        + (erlang:element(10, M) * erlang:element(4, V))}.

-file("src/viva_math/matrix.gleam", 218).
?DOC(" Rotation around the X axis.\n").
-spec mat3_rot_x(float()) -> mat3().
mat3_rot_x(Theta) ->
    C = math:cos(Theta),
    S = math:sin(Theta),
    {mat3, 1.0, +0.0, +0.0, +0.0, C, +0.0 - S, +0.0, S, C}.

-file("src/viva_math/matrix.gleam", 225).
?DOC(" Rotation around the Y axis.\n").
-spec mat3_rot_y(float()) -> mat3().
mat3_rot_y(Theta) ->
    C = math:cos(Theta),
    S = math:sin(Theta),
    {mat3, C, +0.0, S, +0.0, 1.0, +0.0, +0.0 - S, +0.0, C}.

-file("src/viva_math/matrix.gleam", 232).
?DOC(" Rotation around the Z axis.\n").
-spec mat3_rot_z(float()) -> mat3().
mat3_rot_z(Theta) ->
    C = math:cos(Theta),
    S = math:sin(Theta),
    {mat3, C, +0.0 - S, +0.0, S, C, +0.0, +0.0, +0.0, 1.0}.

-file("src/viva_math/matrix.gleam", 273).
?DOC(
    " Real eigenvalues of a 2×2 matrix via the characteristic quadratic.\n"
    "\n"
    " λ² - tr(M)·λ + det(M) = 0  →  λ = (tr ± √(tr² - 4·det)) / 2.\n"
    "\n"
    " Returns `Error` for matrices with non-real eigenvalues (negative\n"
    " discriminant). For symmetric matrices the discriminant is always\n"
    " non-negative, so this always succeeds.\n"
).
-spec mat2_eigenvalues(mat2()) -> {ok, {float(), float()}} | {error, nil}.
mat2_eigenvalues(M) ->
    Tr = mat2_trace(M),
    Det = mat2_determinant(M),
    Discriminant = (Tr * Tr) - (4.0 * Det),
    case Discriminant < +0.0 of
        true ->
            {error, nil};

        false ->
            S = math:sqrt(Discriminant),
            {ok, {(Tr + S) / 2.0, (Tr - S) / 2.0}}
    end.

-file("src/viva_math/matrix.gleam", 347).
-spec sort_three(float(), float(), float()) -> {float(), float(), float()}.
sort_three(A, B, C) ->
    case {A =< B, B =< C, A =< C} of
        {true, true, _} ->
            {A, B, C};

        {true, false, true} ->
            {A, C, B};

        {true, false, false} ->
            {C, A, B};

        {false, true, true} ->
            {B, A, C};

        {false, true, false} ->
            {B, C, A};

        {_, _, _} ->
            {C, B, A}
    end.

-file("src/viva_math/matrix.gleam", 361).
-spec acos_safe(float()) -> float().
acos_safe(X) ->
    math:acos(X).

-file("src/viva_math/matrix.gleam", 294).
?DOC(
    " Real eigenvalues of a **symmetric** 3×3 matrix via Smith's trigonometric\n"
    " closed form. The matrix must be (numerically) symmetric.\n"
    "\n"
    " Algorithm: characteristic cubic det(M - λI) = 0 in the form\n"
    " λ³ - tr·λ² + …  Then shift by p = tr/3 and solve the depressed cubic\n"
    " using trigonometric substitution.\n"
    "\n"
    " Reference: Smith (1961) \"Eigenvalues of a Symmetric 3×3 Matrix\".\n"
).
-spec mat3_symmetric_eigenvalues(mat3()) -> {ok, {float(), float(), float()}} |
    {error, nil}.
mat3_symmetric_eigenvalues(M) ->
    P1 = ((erlang:element(3, M) * erlang:element(3, M)) + (erlang:element(4, M)
    * erlang:element(4, M)))
    + (erlang:element(7, M) * erlang:element(7, M)),
    case P1 =:= +0.0 of
        true ->
            {ok,
                sort_three(
                    erlang:element(2, M),
                    erlang:element(6, M),
                    erlang:element(10, M)
                )};

        false ->
            Q = mat3_trace(M) / 3.0,
            P2 = ((((erlang:element(2, M) - Q) * (erlang:element(2, M) - Q)) + ((erlang:element(
                6,
                M
            )
            - Q)
            * (erlang:element(6, M) - Q)))
            + ((erlang:element(10, M) - Q) * (erlang:element(10, M) - Q)))
            + (2.0 * P1),
            P = math:sqrt(P2 / 6.0),
            Inv_p = case P of
                +0.0 -> +0.0;
                -0.0 -> -0.0;
                Gleam@denominator -> 1.0 / Gleam@denominator
            end,
            B = {mat3,
                Inv_p * (erlang:element(2, M) - Q),
                Inv_p * erlang:element(3, M),
                Inv_p * erlang:element(4, M),
                Inv_p * erlang:element(5, M),
                Inv_p * (erlang:element(6, M) - Q),
                Inv_p * erlang:element(7, M),
                Inv_p * erlang:element(8, M),
                Inv_p * erlang:element(9, M),
                Inv_p * (erlang:element(10, M) - Q)},
            R = mat3_determinant(B) / 2.0,
            R_clamped = case R of
                R@1 when R@1 > 1.0 ->
                    1.0;

                R@2 when R@2 < -1.0 ->
                    -1.0;

                _ ->
                    R
            end,
            Phi = acos_safe(R_clamped) / 3.0,
            Two_pi_over_3 = 2.094395102393195,
            Eig1 = Q + ((2.0 * P) * math:cos(Phi)),
            Eig3 = Q + ((2.0 * P) * math:cos(Phi + Two_pi_over_3)),
            Eig2 = ((3.0 * Q) - Eig1) - Eig3,
            {ok, sort_three(Eig1, Eig2, Eig3)}
    end.

-file("src/viva_math/matrix.gleam", 366).
?DOC(" Diagonal matrix.\n").
-spec mat3_diagonal(float(), float(), float()) -> mat3().
mat3_diagonal(A, B, C) ->
    {mat3, A, +0.0, +0.0, +0.0, B, +0.0, +0.0, +0.0, C}.

-file("src/viva_math/matrix.gleam", 378).
-spec mat4_identity() -> mat4().
mat4_identity() ->
    {mat4,
        [[1.0, +0.0, +0.0, +0.0],
            [+0.0, 1.0, +0.0, +0.0],
            [+0.0, +0.0, 1.0, +0.0],
            [+0.0, +0.0, +0.0, 1.0]]}.

-file("src/viva_math/matrix.gleam", 387).
-spec mat4_zero() -> mat4().
mat4_zero() ->
    {mat4, gleam@list:repeat([+0.0, +0.0, +0.0, +0.0], 4)}.

-file("src/viva_math/matrix.gleam", 392).
?DOC(" Translation homogeneous transform.\n").
-spec mat4_translation(float(), float(), float()) -> mat4().
mat4_translation(Tx, Ty, Tz) ->
    {mat4,
        [[1.0, +0.0, +0.0, Tx],
            [+0.0, 1.0, +0.0, Ty],
            [+0.0, +0.0, 1.0, Tz],
            [+0.0, +0.0, +0.0, 1.0]]}.

-file("src/viva_math/matrix.gleam", 402).
?DOC(" Non-uniform scale.\n").
-spec mat4_scale(float(), float(), float()) -> mat4().
mat4_scale(Sx, Sy, Sz) ->
    {mat4,
        [[Sx, +0.0, +0.0, +0.0],
            [+0.0, Sy, +0.0, +0.0],
            [+0.0, +0.0, Sz, +0.0],
            [+0.0, +0.0, +0.0, 1.0]]}.

-file("src/viva_math/matrix.gleam", 422).
?DOC(
    " Build a MatN from a row-major list of rows.\n"
    "\n"
    " Errors if rows is empty, has empty rows, or rows of inconsistent length.\n"
).
-spec matn_from_rows(list(list(float()))) -> {ok, mat_n()} | {error, nil}.
matn_from_rows(Rows) ->
    case Rows of
        [] ->
            {error, nil};

        [First | _] ->
            N_cols = erlang:length(First),
            case N_cols =:= 0 of
                true ->
                    {error, nil};

                false ->
                    Consistent = gleam@list:all(
                        Rows,
                        fun(R) -> erlang:length(R) =:= N_cols end
                    ),
                    case Consistent of
                        false ->
                            {error, nil};

                        true ->
                            {ok, {mat_n, erlang:length(Rows), N_cols, Rows}}
                    end
            end
    end.

-file("src/viva_math/matrix.gleam", 442).
?DOC(" Zero matrix of given shape.\n").
-spec matn_zeros(integer(), integer()) -> mat_n().
matn_zeros(Rows, Cols) ->
    Row = gleam@list:repeat(+0.0, Cols),
    {mat_n, Rows, Cols, gleam@list:repeat(Row, Rows)}.

-file("src/viva_math/matrix.gleam", 589).
-spec range_loop(integer(), integer(), list(integer())) -> list(integer()).
range_loop(From, To, Acc) ->
    case From > To of
        true ->
            lists:reverse(Acc);

        false ->
            range_loop(From + 1, To, [From | Acc])
    end.

-file("src/viva_math/matrix.gleam", 585).
-spec range_int(integer(), integer()) -> list(integer()).
range_int(From, To) ->
    range_loop(From, To, []).

-file("src/viva_math/matrix.gleam", 448).
?DOC(" Identity matrix of size n.\n").
-spec matn_identity(integer()) -> mat_n().
matn_identity(N) ->
    Data = begin
        _pipe = range_int(0, N - 1),
        gleam@list:map(_pipe, fun(I) -> _pipe@1 = range_int(0, N - 1),
                gleam@list:map(_pipe@1, fun(J) -> case I =:= J of
                            true ->
                                1.0;

                            false ->
                                +0.0
                        end end) end)
    end,
    {mat_n, N, N, Data}.

-file("src/viva_math/matrix.gleam", 571).
-spec list_at(list(float()), integer()) -> {ok, float()} | {error, nil}.
list_at(Xs, Idx) ->
    case {Xs, Idx} of
        {[], _} ->
            {error, nil};

        {[X | _], 0} ->
            {ok, X};

        {[_ | Rest], N} ->
            list_at(Rest, N - 1)
    end.

-file("src/viva_math/matrix.gleam", 469).
-spec transpose_lists(list(list(float()))) -> list(list(float())).
transpose_lists(Rows) ->
    case Rows of
        [] ->
            [];

        [First | _] ->
            N = erlang:length(First),
            _pipe = range_int(0, N - 1),
            gleam@list:map(
                _pipe,
                fun(Col) ->
                    gleam@list:map(Rows, fun(Row) -> case list_at(Row, Col) of
                                {ok, V} ->
                                    V;

                                {error, _} ->
                                    +0.0
                            end end)
                end
            )
    end.

-file("src/viva_math/matrix.gleam", 464).
?DOC(" Transpose of a MatN. O(rows · cols).\n").
-spec matn_transpose(mat_n()) -> mat_n().
matn_transpose(M) ->
    Rows = transpose_lists(erlang:element(4, M)),
    {mat_n, erlang:element(3, M), erlang:element(2, M), Rows}.

-file("src/viva_math/matrix.gleam", 502).
-spec dot_lists(list(float()), list(float())) -> float().
dot_lists(A, B) ->
    _pipe = gleam@list:zip(A, B),
    gleam@list:fold(
        _pipe,
        +0.0,
        fun(Acc, Pair) ->
            Acc + (erlang:element(1, Pair) * erlang:element(2, Pair))
        end
    ).

-file("src/viva_math/matrix.gleam", 488).
?DOC(" Matrix product A·B. Errors on shape mismatch.\n").
-spec matn_mul(mat_n(), mat_n()) -> {ok, mat_n()} | {error, nil}.
matn_mul(A, B) ->
    case erlang:element(3, A) =:= erlang:element(2, B) of
        false ->
            {error, nil};

        true ->
            Bt = matn_transpose(B),
            New_rows = gleam@list:map(
                erlang:element(4, A),
                fun(Row_a) ->
                    gleam@list:map(
                        erlang:element(4, Bt),
                        fun(Row_b) -> dot_lists(Row_a, Row_b) end
                    )
                end
            ),
            {ok, {mat_n, erlang:element(2, A), erlang:element(3, B), New_rows}}
    end.

-file("src/viva_math/matrix.gleam", 508).
?DOC(" Add two matrices. Errors on shape mismatch.\n").
-spec matn_add(mat_n(), mat_n()) -> {ok, mat_n()} | {error, nil}.
matn_add(A, B) ->
    case (erlang:element(2, A) =:= erlang:element(2, B)) andalso (erlang:element(
        3,
        A
    )
    =:= erlang:element(3, B)) of
        false ->
            {error, nil};

        true ->
            New_rows = begin
                _pipe = gleam@list:zip(
                    erlang:element(4, A),
                    erlang:element(4, B)
                ),
                gleam@list:map(
                    _pipe,
                    fun(Pair) ->
                        _pipe@1 = gleam@list:zip(
                            erlang:element(1, Pair),
                            erlang:element(2, Pair)
                        ),
                        gleam@list:map(
                            _pipe@1,
                            fun(P) ->
                                erlang:element(1, P) + erlang:element(2, P)
                            end
                        )
                    end
                )
            end,
            {ok, {mat_n, erlang:element(2, A), erlang:element(3, A), New_rows}}
    end.

-file("src/viva_math/matrix.gleam", 524).
?DOC(" Scalar multiplication.\n").
-spec matn_scale(mat_n(), float()) -> mat_n().
matn_scale(M, S) ->
    New_rows = gleam@list:map(
        erlang:element(4, M),
        fun(Row) -> gleam@list:map(Row, fun(X) -> X * S end) end
    ),
    {mat_n, erlang:element(2, M), erlang:element(3, M), New_rows}.

-file("src/viva_math/matrix.gleam", 537).
-spec trace_loop(list(list(float())), integer(), float()) -> float().
trace_loop(Rows, Idx, Acc) ->
    case Rows of
        [] ->
            Acc;

        [Row | Rest] ->
            V = case list_at(Row, Idx) of
                {ok, X} ->
                    X;

                {error, _} ->
                    +0.0
            end,
            trace_loop(Rest, Idx + 1, Acc + V)
    end.

-file("src/viva_math/matrix.gleam", 530).
?DOC(" Trace (sum of diagonal). Defined for square matrices.\n").
-spec matn_trace(mat_n()) -> {ok, float()} | {error, nil}.
matn_trace(M) ->
    case erlang:element(2, M) =:= erlang:element(3, M) of
        false ->
            {error, nil};

        true ->
            {ok, trace_loop(erlang:element(4, M), 0, +0.0)}
    end.

-file("src/viva_math/matrix.gleam", 551).
?DOC(" Frobenius norm: √(Σ aᵢⱼ²).\n").
-spec matn_frobenius(mat_n()) -> float().
matn_frobenius(M) ->
    Sum_sq = gleam@list:fold(
        erlang:element(4, M),
        +0.0,
        fun(Acc, Row) ->
            Acc + gleam@list:fold(Row, +0.0, fun(S, X) -> S + (X * X) end)
        end
    ),
    math:sqrt(Sum_sq).

-file("src/viva_math/matrix.gleam", 560).
?DOC(" Matrix × column-vector product. Returns the resulting vector as a list.\n").
-spec matn_mul_vec(mat_n(), list(float())) -> {ok, list(float())} | {error, nil}.
matn_mul_vec(M, V) ->
    case erlang:element(3, M) =:= erlang:length(V) of
        false ->
            {error, nil};

        true ->
            {ok,
                gleam@list:map(
                    erlang:element(4, M),
                    fun(Row) -> dot_lists(Row, V) end
                )}
    end.