viva_math

-module(viva_math@calculus).
-compile([no_auto_import, nowarn_unused_vars, nowarn_unused_function, nowarn_nomatch, inline]).
-define(FILEPATH, "src/viva_math/calculus.gleam").
-export([forward_diff/3, backward_diff/3, central_diff/3, five_point_diff/3, second_derivative/3, gradient/3, trapezoid/4, simpson/4, romberg/4, norm_l2/1]).

-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(
    " Numerical calculus: differentiation and integration.\n"
    "\n"
    " Closed-form symbolic calculus is out of scope. This module provides\n"
    " finite-difference approximations for derivatives and quadrature rules\n"
    " for integrals.\n"
    "\n"
    " ## Order of accuracy\n"
    "\n"
    " | Function                     | Order   | Notes                       |\n"
    " | ---------------------------- | ------- | --------------------------- |\n"
    " | `forward_diff` / `backward`  | O(h)    | One-sided                   |\n"
    " | `central_diff`               | O(h²)   | Recommended default         |\n"
    " | `five_point_diff`            | O(h⁴)   | Five evaluations, very accurate |\n"
    " | `second_derivative`          | O(h²)   | Central stencil             |\n"
    " | `trapezoid` / `simpson`      | O(h²/h⁴) | Composite rules             |\n"
    " | `romberg`                    | adaptive | Richardson extrapolation   |\n"
).

-file("src/viva_math/calculus.gleam", 26).
?DOC(" Forward difference (f(x+h) - f(x)) / h. O(h) accurate.\n").
-spec forward_diff(fun((float()) -> float()), float(), float()) -> float().
forward_diff(F, X, H) ->
    case H of
        +0.0 -> +0.0;
        -0.0 -> -0.0;
        Gleam@denominator -> (F(X + H) - F(X)) / Gleam@denominator
    end.

-file("src/viva_math/calculus.gleam", 31).
?DOC(" Backward difference (f(x) - f(x-h)) / h. O(h) accurate.\n").
-spec backward_diff(fun((float()) -> float()), float(), float()) -> float().
backward_diff(F, X, H) ->
    case H of
        +0.0 -> +0.0;
        -0.0 -> -0.0;
        Gleam@denominator -> (F(X) - F(X - H)) / Gleam@denominator
    end.

-file("src/viva_math/calculus.gleam", 36).
?DOC(" Central difference (f(x+h) - f(x-h)) / (2h). O(h²) accurate.\n").
-spec central_diff(fun((float()) -> float()), float(), float()) -> float().
central_diff(F, X, H) ->
    case (2.0 * H) of
        +0.0 -> +0.0;
        -0.0 -> -0.0;
        Gleam@denominator -> (F(X + H) - F(X - H)) / Gleam@denominator
    end.

-file("src/viva_math/calculus.gleam", 43).
?DOC(
    " Five-point stencil O(h⁴): the gold standard for smooth functions.\n"
    "\n"
    " f'(x) ≈ (-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)) / (12h)\n"
).
-spec five_point_diff(fun((float()) -> float()), float(), float()) -> float().
five_point_diff(F, X, H) ->
    F_p2 = F(X + (2.0 * H)),
    F_p1 = F(X + H),
    F_m1 = F(X - H),
    F_m2 = F(X - (2.0 * H)),
    case (12.0 * H) of
        +0.0 -> +0.0;
        -0.0 -> -0.0;
        Gleam@denominator -> ((((+0.0 - F_p2) + (8.0 * F_p1)) - (8.0 * F_m1)) + F_m2)
        / Gleam@denominator
    end.

-file("src/viva_math/calculus.gleam", 54).
?DOC(
    " Second derivative via central stencil. O(h²).\n"
    "\n"
    " f''(x) ≈ (f(x+h) - 2f(x) + f(x-h)) / h²\n"
).
-spec second_derivative(fun((float()) -> float()), float(), float()) -> float().
second_derivative(F, X, H) ->
    case (H * H) of
        +0.0 -> +0.0;
        -0.0 -> -0.0;
        Gleam@denominator -> ((F(X + H) - (2.0 * F(X))) + F(X - H)) / Gleam@denominator
    end.

-file("src/viva_math/calculus.gleam", 92).
-spec perturb_loop(list(float()), integer(), float(), integer(), list(float())) -> list(float()).
perturb_loop(Point, Index, Delta, I, Acc) ->
    case Point of
        [] ->
            lists:reverse(Acc);

        [X | Rest] ->
            Val = case I =:= Index of
                true ->
                    X + Delta;

                false ->
                    X
            end,
            perturb_loop(Rest, Index, Delta, I + 1, [Val | Acc])
    end.

-file("src/viva_math/calculus.gleam", 88).
-spec perturb(list(float()), integer(), float()) -> list(float()).
perturb(Point, Index, Delta) ->
    perturb_loop(Point, Index, Delta, 0, []).

-file("src/viva_math/calculus.gleam", 69).
-spec gradient_loop(
    fun((list(float())) -> float()),
    list(float()),
    float(),
    integer(),
    integer(),
    list(float())
) -> list(float()).
gradient_loop(F, Point, H, I, N, Acc) ->
    case I >= N of
        true ->
            lists:reverse(Acc);

        false ->
            Plus = perturb(Point, I, H),
            Minus = perturb(Point, I, +0.0 - H),
            Partial = case (2.0 * H) of
                +0.0 -> +0.0;
                -0.0 -> -0.0;
                Gleam@denominator -> (F(Plus) - F(Minus)) / Gleam@denominator
            end,
            gradient_loop(F, Point, H, I + 1, N, [Partial | Acc])
    end.

-file("src/viva_math/calculus.gleam", 61).
?DOC(
    " Gradient of a multivariate function `f: List(Float) -> Float`.\n"
    "\n"
    " Uses central differences along each axis with step `h`.\n"
).
-spec gradient(fun((list(float())) -> float()), list(float()), float()) -> list(float()).
gradient(F, Point, H) ->
    gradient_loop(F, Point, H, 0, erlang:length(Point), []).

-file("src/viva_math/calculus.gleam", 260).
-spec int_to_float(integer()) -> float().
int_to_float(N) ->
    erlang:float(N).

-file("src/viva_math/calculus.gleam", 129).
-spec trapezoid_inner_sum(
    fun((float()) -> float()),
    float(),
    float(),
    integer(),
    integer()
) -> float().
trapezoid_inner_sum(F, A, H, I, N) ->
    case I >= N of
        true ->
            +0.0;

        false ->
            F(A + (int_to_float(I) * H)) + trapezoid_inner_sum(
                F,
                A,
                H,
                I + 1,
                N
            )
    end.

-file("src/viva_math/calculus.gleam", 118).
?DOC(
    " Composite trapezoid rule over n equal subintervals.\n"
    "\n"
    " ∫f ≈ h · ( ½f(a) + Σ f(a + i·h) + ½f(b) ),  h = (b - a) / n\n"
).
-spec trapezoid(fun((float()) -> float()), float(), float(), integer()) -> float().
trapezoid(F, A, B, N) ->
    case N < 1 of
        true ->
            +0.0;

        false ->
            H = case int_to_float(N) of
                +0.0 -> +0.0;
                -0.0 -> -0.0;
                Gleam@denominator -> (B - A) / Gleam@denominator
            end,
            Inner = trapezoid_inner_sum(F, A, H, 1, N),
            H * (((0.5 * F(A)) + Inner) + (0.5 * F(B)))
    end.

-file("src/viva_math/calculus.gleam", 163).
-spec simpson_loop(
    fun((float()) -> float()),
    float(),
    float(),
    integer(),
    integer(),
    float()
) -> float().
simpson_loop(F, A, H, I, N, Acc) ->
    case I >= N of
        true ->
            Acc;

        false ->
            X = A + (int_to_float(I) * H),
            Weight = case (I rem 2) =:= 0 of
                true ->
                    2.0;

                false ->
                    4.0
            end,
            simpson_loop(F, A, H, I + 1, N, Acc + (Weight * F(X)))
    end.

-file("src/viva_math/calculus.gleam", 146).
?DOC(
    " Composite Simpson's rule. n must be even and ≥ 2.\n"
    "\n"
    " ∫f ≈ (h/3) · (f₀ + 4f₁ + 2f₂ + 4f₃ + ... + 4fₙ₋₁ + fₙ)\n"
).
-spec simpson(fun((float()) -> float()), float(), float(), integer()) -> {ok,
        float()} |
    {error, nil}.
simpson(F, A, B, N) ->
    case {N < 2, (N rem 2) /= 0} of
        {true, _} ->
            {error, nil};

        {_, true} ->
            {error, nil};

        {_, _} ->
            H = case int_to_float(N) of
                +0.0 -> +0.0;
                -0.0 -> -0.0;
                Gleam@denominator -> (B - A) / Gleam@denominator
            end,
            Total = simpson_loop(F, A, H, 1, N, +0.0),
            {ok, (H / 3.0) * ((F(A) + F(B)) + Total)}
    end.

-file("src/viva_math/calculus.gleam", 264).
-spec int_pow(integer(), integer()) -> integer().
int_pow(Base, Exp) ->
    case Exp of
        0 ->
            1;

        _ ->
            Base * int_pow(Base, Exp - 1)
    end.

-file("src/viva_math/calculus.gleam", 241).
-spec combine_pairs(list(float()), integer()) -> list(float()).
combine_pairs(Row, K) ->
    case Row of
        [A, B | Rest] ->
            Factor = begin
                _pipe = int_pow(4, K),
                int_to_float(_pipe)
            end,
            Improved = case (Factor - 1.0) of
                +0.0 -> +0.0;
                -0.0 -> -0.0;
                Gleam@denominator -> ((Factor * B) - A) / Gleam@denominator
            end,
            [Improved | combine_pairs([B | Rest], K)];

        _ ->
            []
    end.

-file("src/viva_math/calculus.gleam", 230).
-spec romberg_extrapolate(list(float()), integer()) -> float().
romberg_extrapolate(Row, K) ->
    case Row of
        [] ->
            +0.0;

        [Single] ->
            Single;

        _ ->
            Next_row = combine_pairs(Row, K),
            romberg_extrapolate(Next_row, K + 1)
    end.

-file("src/viva_math/calculus.gleam", 212).
-spec build_row(
    fun((float()) -> float()),
    float(),
    float(),
    integer(),
    integer(),
    list(float())
) -> list(float()).
build_row(F, A, B, I, Levels, Acc) ->
    case I >= Levels of
        true ->
            lists:reverse(Acc);

        false ->
            N = int_pow(2, I),
            Estimate = trapezoid(F, A, B, N),
            build_row(F, A, B, I + 1, Levels, [Estimate | Acc])
    end.

-file("src/viva_math/calculus.gleam", 203).
-spec romberg_row(fun((float()) -> float()), float(), float(), integer()) -> list(float()).
romberg_row(F, A, B, Levels) ->
    build_row(F, A, B, 0, Levels, []).

-file("src/viva_math/calculus.gleam", 188).
?DOC(
    " Romberg integration via Richardson extrapolation.\n"
    "\n"
    " Repeatedly halves `h` and combines trapezoid estimates to cancel error\n"
    " terms. `levels` controls the depth (typically 5-8 suffices).\n"
).
-spec romberg(fun((float()) -> float()), float(), float(), integer()) -> float().
romberg(F, A, B, Levels) ->
    case Levels < 1 of
        true ->
            trapezoid(F, A, B, 1);

        false ->
            First_row = romberg_row(F, A, B, Levels),
            romberg_extrapolate(First_row, 1)
    end.

-file("src/viva_math/calculus.gleam", 272).
-spec norm_l2(list(float())) -> float().
norm_l2(Values) ->
    math:sqrt(gleam@list:fold(Values, +0.0, fun(Acc, X) -> Acc + (X * X) end)).