viva_math

-module(viva_math@quaternion).
-compile([no_auto_import, nowarn_unused_vars, nowarn_unused_function, nowarn_nomatch, inline]).
-define(FILEPATH, "src/viva_math/quaternion.gleam").
-export([identity/0, from_axis_angle/2, raw/4, mul/2, conjugate/1, magnitude/1, magnitude_squared/1, normalize/1, inverse/1, dot/2, rotate/2, nlerp/3, slerp/3, to_axis_angle/1, is_close/3]).
-export_type([quaternion/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(
    " Unit quaternions for 3-D rotation.\n"
    "\n"
    " A quaternion `w + xi + yj + zk` extends complex numbers to four\n"
    " dimensions. Unit quaternions provide a numerically stable,\n"
    " gimbal-lock-free representation of 3-D rotation that interpolates\n"
    " smoothly via SLERP.\n"
    "\n"
    " ## Conventions\n"
    "\n"
    " - `q.w` is the scalar (real) part.\n"
    " - `(q.x, q.y, q.z)` is the vector (imaginary) part.\n"
    " - Rotation by angle `θ` around unit axis `(ax, ay, az)`:\n"
    "   `q = cos(θ/2) + sin(θ/2)·(ax·i + ay·j + az·k)`.\n"
    " - Right-handed coordinates.\n"
    "\n"
    " ## When to use\n"
    "\n"
    " - Smooth interpolation between two orientations (`slerp`).\n"
    " - Composing rotations without matrix multiplication overhead.\n"
    " - Avoiding gimbal lock present in Euler-angle representations.\n"
).

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

-file("src/viva_math/quaternion.gleam", 35).
?DOC(" The identity quaternion (no rotation).\n").
-spec identity() -> quaternion().
identity() ->
    {quaternion, 1.0, +0.0, +0.0, +0.0}.

-file("src/viva_math/quaternion.gleam", 43).
?DOC(
    " Build a unit quaternion from an axis-angle representation.\n"
    "\n"
    " `axis` need not be normalised — this function normalises it. `theta` is\n"
    " the rotation angle in radians.\n"
).
-spec from_axis_angle(viva_math@vector:vec3(), float()) -> quaternion().
from_axis_angle(Axis, Theta) ->
    N = viva_math@vector:normalize(Axis),
    Half = Theta / 2.0,
    S = math:sin(Half),
    {quaternion,
        math:cos(Half),
        erlang:element(2, N) * S,
        erlang:element(3, N) * S,
        erlang:element(4, N) * S}.

-file("src/viva_math/quaternion.gleam", 51).
?DOC(" Build a quaternion from raw components without normalisation.\n").
-spec raw(float(), float(), float(), float()) -> quaternion().
raw(W, X, Y, Z) ->
    {quaternion, W, X, Y, Z}.

-file("src/viva_math/quaternion.gleam", 60).
?DOC(" Hamilton product `a · b`. Non-commutative.\n").
-spec mul(quaternion(), quaternion()) -> quaternion().
mul(A, B) ->
    {quaternion,
        (((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(2, A) * erlang:element(3, B)) + (erlang:element(3, A)
        * erlang:element(2, B)))
        + (erlang:element(4, A) * erlang:element(5, B)))
        - (erlang:element(5, A) * erlang:element(4, B)),
        (((erlang:element(2, A) * erlang:element(4, B)) - (erlang:element(3, A)
        * erlang:element(5, B)))
        + (erlang:element(4, A) * erlang:element(2, B)))
        + (erlang:element(5, A) * erlang:element(3, B)),
        (((erlang:element(2, A) * erlang:element(5, B)) + (erlang:element(3, A)
        * erlang:element(4, B)))
        - (erlang:element(4, A) * erlang:element(3, B)))
        + (erlang:element(5, A) * erlang:element(2, B))}.

-file("src/viva_math/quaternion.gleam", 70).
?DOC(" Quaternion conjugate `(w, -x, -y, -z)`. For unit quaternions this is the inverse.\n").
-spec conjugate(quaternion()) -> quaternion().
conjugate(Q) ->
    {quaternion,
        erlang:element(2, Q),
        +0.0 - erlang:element(3, Q),
        +0.0 - erlang:element(4, Q),
        +0.0 - erlang:element(5, Q)}.

-file("src/viva_math/quaternion.gleam", 75).
?DOC(" Magnitude / Euclidean norm of a quaternion.\n").
-spec magnitude(quaternion()) -> float().
magnitude(Q) ->
    math:sqrt(
        (((erlang:element(2, Q) * erlang:element(2, Q)) + (erlang:element(3, Q)
        * erlang:element(3, Q)))
        + (erlang:element(4, Q) * erlang:element(4, Q)))
        + (erlang:element(5, Q) * erlang:element(5, Q))
    ).

-file("src/viva_math/quaternion.gleam", 80).
?DOC(" Squared magnitude (cheaper than `magnitude` when only comparison is needed).\n").
-spec magnitude_squared(quaternion()) -> float().
magnitude_squared(Q) ->
    (((erlang:element(2, Q) * erlang:element(2, Q)) + (erlang:element(3, Q) * erlang:element(
        3,
        Q
    )))
    + (erlang:element(4, Q) * erlang:element(4, Q)))
    + (erlang:element(5, Q) * erlang:element(5, Q)).

-file("src/viva_math/quaternion.gleam", 85).
?DOC(" Normalise to unit length. Returns identity if the input is zero.\n").
-spec normalize(quaternion()) -> quaternion().
normalize(Q) ->
    M = magnitude(Q),
    case M =:= +0.0 of
        true ->
            identity();

        false ->
            {quaternion, case M of
                    +0.0 -> +0.0;
                    -0.0 -> -0.0;
                    Gleam@denominator -> erlang:element(2, Q) / Gleam@denominator
                end, case M of
                    +0.0 -> +0.0;
                    -0.0 -> -0.0;
                    Gleam@denominator@1 -> erlang:element(3, Q) / Gleam@denominator@1
                end, case M of
                    +0.0 -> +0.0;
                    -0.0 -> -0.0;
                    Gleam@denominator@2 -> erlang:element(4, Q) / Gleam@denominator@2
                end, case M of
                    +0.0 -> +0.0;
                    -0.0 -> -0.0;
                    Gleam@denominator@3 -> erlang:element(5, Q) / Gleam@denominator@3
                end}
    end.

-file("src/viva_math/quaternion.gleam", 94).
?DOC(" Inverse `conj(q) / |q|²`. For unit quaternions equals the conjugate.\n").
-spec inverse(quaternion()) -> quaternion().
inverse(Q) ->
    M2 = magnitude_squared(Q),
    case M2 =:= +0.0 of
        true ->
            identity();

        false ->
            C = conjugate(Q),
            {quaternion, case M2 of
                    +0.0 -> +0.0;
                    -0.0 -> -0.0;
                    Gleam@denominator -> erlang:element(2, C) / Gleam@denominator
                end, case M2 of
                    +0.0 -> +0.0;
                    -0.0 -> -0.0;
                    Gleam@denominator@1 -> erlang:element(3, C) / Gleam@denominator@1
                end, case M2 of
                    +0.0 -> +0.0;
                    -0.0 -> -0.0;
                    Gleam@denominator@2 -> erlang:element(4, C) / Gleam@denominator@2
                end, case M2 of
                    +0.0 -> +0.0;
                    -0.0 -> -0.0;
                    Gleam@denominator@3 -> erlang:element(5, C) / Gleam@denominator@3
                end}
    end.

-file("src/viva_math/quaternion.gleam", 106).
?DOC(" Quaternion-quaternion dot product (treating each as a 4-vector).\n").
-spec dot(quaternion(), quaternion()) -> float().
dot(A, B) ->
    (((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/quaternion.gleam", 118).
?DOC(
    " Rotate a 3-D vector `v` by unit quaternion `q`.\n"
    "\n"
    " Computes `q · v · q⁻¹` using the optimised Rodrigues-style formula that\n"
    " skips two of the quaternion products.\n"
).
-spec rotate(quaternion(), viva_math@vector:vec3()) -> viva_math@vector:vec3().
rotate(Q, V) ->
    S = erlang:element(2, Q),
    U = {vec3, erlang:element(3, Q), erlang:element(4, Q), erlang:element(5, Q)},
    U_dot_v = viva_math@vector:dot(U, V),
    U_dot_u = viva_math@vector:dot(U, U),
    Cross_uv = viva_math@vector:cross(U, V),
    Term1 = viva_math@vector:scale(U, 2.0 * U_dot_v),
    Term2 = viva_math@vector:scale(V, (S * S) - U_dot_u),
    Term3 = viva_math@vector:scale(Cross_uv, 2.0 * S),
    viva_math@vector:add(viva_math@vector:add(Term1, Term2), Term3).

-file("src/viva_math/quaternion.gleam", 141).
?DOC(
    " Linear interpolation (LERP) between two quaternions, then normalise.\n"
    "\n"
    " Cheaper than SLERP but yields non-uniform angular speed. Acceptable for\n"
    " small angular distances (< ~30°). For large or critical interpolations\n"
    " use `slerp`.\n"
).
-spec nlerp(quaternion(), quaternion(), float()) -> quaternion().
nlerp(A, B, T) ->
    T_clamped = viva_math@scalar:clamp_unit(T),
    B_signed = case dot(A, B) < +0.0 of
        true ->
            {quaternion,
                +0.0 - erlang:element(2, B),
                +0.0 - erlang:element(3, B),
                +0.0 - erlang:element(4, B),
                +0.0 - erlang:element(5, B)};

        false ->
            B
    end,
    One_minus_t = 1.0 - T_clamped,
    normalize(
        {quaternion,
            (One_minus_t * erlang:element(2, A)) + (T_clamped * erlang:element(
                2,
                B_signed
            )),
            (One_minus_t * erlang:element(3, A)) + (T_clamped * erlang:element(
                3,
                B_signed
            )),
            (One_minus_t * erlang:element(4, A)) + (T_clamped * erlang:element(
                4,
                B_signed
            )),
            (One_minus_t * erlang:element(5, A)) + (T_clamped * erlang:element(
                5,
                B_signed
            ))}
    ).

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

-file("src/viva_math/quaternion.gleam", 162).
?DOC(
    " Spherical linear interpolation (SLERP) — constant angular velocity.\n"
    "\n"
    " Falls back to `nlerp` when the angle between the quaternions is very\n"
    " small (`sin(θ) < 1e-6`), where SLERP becomes numerically unstable.\n"
).
-spec slerp(quaternion(), quaternion(), float()) -> quaternion().
slerp(A, B, T) ->
    T_clamped = viva_math@scalar:clamp_unit(T),
    Cos_theta = dot(A, B),
    {B_signed, Cos_theta_signed} = case Cos_theta < +0.0 of
        true ->
            {{quaternion,
                    +0.0 - erlang:element(2, B),
                    +0.0 - erlang:element(3, B),
                    +0.0 - erlang:element(4, B),
                    +0.0 - erlang:element(5, B)},
                +0.0 - Cos_theta};

        false ->
            {B, Cos_theta}
    end,
    case Cos_theta_signed > 0.9995 of
        true ->
            nlerp(A, B_signed, T_clamped);

        false ->
            Theta = acos_safe(Cos_theta_signed),
            Sin_theta = math:sin(Theta),
            Ratio_a = case Sin_theta of
                +0.0 -> +0.0;
                -0.0 -> -0.0;
                Gleam@denominator -> math:sin((1.0 - T_clamped) * Theta) / Gleam@denominator
            end,
            Ratio_b = case Sin_theta of
                +0.0 -> +0.0;
                -0.0 -> -0.0;
                Gleam@denominator@1 -> math:sin(T_clamped * Theta) / Gleam@denominator@1
            end,
            {quaternion,
                (Ratio_a * erlang:element(2, A)) + (Ratio_b * erlang:element(
                    2,
                    B_signed
                )),
                (Ratio_a * erlang:element(3, A)) + (Ratio_b * erlang:element(
                    3,
                    B_signed
                )),
                (Ratio_a * erlang:element(4, A)) + (Ratio_b * erlang:element(
                    4,
                    B_signed
                )),
                (Ratio_a * erlang:element(5, A)) + (Ratio_b * erlang:element(
                    5,
                    B_signed
                ))}
    end.

-file("src/viva_math/quaternion.gleam", 196).
?DOC(
    " Quaternion → axis-angle. Returns `(axis, angle)` with `axis` unit-length.\n"
    " When the rotation is identity returns `(z-axis, 0)` by convention.\n"
).
-spec to_axis_angle(quaternion()) -> {viva_math@vector:vec3(), float()}.
to_axis_angle(Q) ->
    Qn = normalize(Q),
    W_clamped = case erlang:element(2, Qn) of
        W when W > 1.0 ->
            1.0;

        W@1 when W@1 < -1.0 ->
            -1.0;

        W@2 ->
            W@2
    end,
    Angle = 2.0 * acos_safe(W_clamped),
    S = math:sqrt(1.0 - (W_clamped * W_clamped)),
    case S < 1.0e-9 of
        true ->
            {{vec3, +0.0, +0.0, 1.0}, +0.0};

        false ->
            {{vec3, case S of
                        +0.0 -> +0.0;
                        -0.0 -> -0.0;
                        Gleam@denominator -> erlang:element(3, Qn) / Gleam@denominator
                    end, case S of
                        +0.0 -> +0.0;
                        -0.0 -> -0.0;
                        Gleam@denominator@1 -> erlang:element(4, Qn) / Gleam@denominator@1
                    end, case S of
                        +0.0 -> +0.0;
                        -0.0 -> -0.0;
                        Gleam@denominator@2 -> erlang:element(5, Qn) / Gleam@denominator@2
                    end}, Angle}
    end.

-file("src/viva_math/quaternion.gleam", 212).
?DOC(" Approximate equality up to a tolerance.\n").
-spec is_close(quaternion(), quaternion(), float()) -> boolean().
is_close(A, B, Tol) ->
    (((gleam@float:absolute_value(erlang:element(2, A) - erlang:element(2, B))
    =< Tol)
    andalso (gleam@float:absolute_value(
        erlang:element(3, A) - erlang:element(3, B)
    )
    =< Tol))
    andalso (gleam@float:absolute_value(
        erlang:element(4, A) - erlang:element(4, B)
    )
    =< Tol))
    andalso (gleam@float:absolute_value(
        erlang:element(5, A) - erlang:element(5, B)
    )
    =< Tol).