viva_math

-module(viva_math@bench_precision).
-compile([no_auto_import, nowarn_unused_vars, nowarn_unused_function, nowarn_nomatch, inline]).
-define(FILEPATH, "src/viva_math/bench_precision.gleam").
-export([main/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(
    " Benchmarks for the `precision` module.\n"
    "\n"
    " Runs with `gleam run -m precision_bench`. Compares naive `list.fold`\n"
    " summation against Neumaier compensated summation across input sizes\n"
    " and condition numbers.\n"
).

-file("src/viva_math/bench_precision.gleam", 82).
-spec lcg(float()) -> float().
lcg(X) ->
    V = math:fmod((X * 16807.0) + 0.123456789, 2.0) - 1.0,
    V.

-file("src/viva_math/bench_precision.gleam", 71).
-spec build_loop(integer(), float(), list(float())) -> list(float()).
build_loop(N, Seed, Acc) ->
    case N =< 0 of
        true ->
            Acc;

        false ->
            Next = lcg(Seed),
            build_loop(N - 1, Next, [Next | Acc])
    end.

-file("src/viva_math/bench_precision.gleam", 67).
-spec build_random(integer(), float()) -> list(float()).
build_random(N, Seed) ->
    build_loop(N, Seed, []).

-file("src/viva_math/bench_precision.gleam", 11).
-spec main() -> nil.
main() ->
    gleam_stdlib:println(<<"\n=== viva_math precision benchmarks ===\n"/utf8>>),
    Sizes = [100, 1000, 10000, 100000],
    gleam_stdlib:println(
        <<"[1] Sum accuracy on adversarial input [1, 1e100, 1, -1e100] — exact = 2.0"/utf8>>
    ),
    Pathological = [1.0, 1.0e100, 1.0, -1.0e100],
    Naive = gleam@list:fold(Pathological, +0.0, fun(Acc, X) -> Acc + X end),
    Kahan_val = viva_math@precision:kahan_sum(Pathological),
    Neumaier_val = viva_math@precision:neumaier_sum(Pathological),
    Fsum_val = viva_math@precision:fsum(Pathological),
    gleam_stdlib:println(
        <<"  naive sum     -> "/utf8, (erlang:float_to_binary(Naive))/binary>>
    ),
    gleam_stdlib:println(
        <<"  kahan_sum     -> "/utf8,
            (erlang:float_to_binary(Kahan_val))/binary>>
    ),
    gleam_stdlib:println(
        <<"  neumaier_sum  -> "/utf8,
            (erlang:float_to_binary(Neumaier_val))/binary>>
    ),
    gleam_stdlib:println(
        <<"  fsum          -> "/utf8,
            (erlang:float_to_binary(Fsum_val))/binary>>
    ),
    gleam_stdlib:println(
        <<"\n[2] Sum throughput by input size (raw timings)"/utf8>>
    ),
    gleam@list:each(
        Sizes,
        fun(N) ->
            Xs = build_random(N, 0.123),
            T0 = erlang:monotonic_time(),
            _ = gleam@list:fold(Xs, +0.0, fun(Acc@1, X@1) -> Acc@1 + X@1 end),
            T1 = erlang:monotonic_time(),
            _ = viva_math@precision:neumaier_sum(Xs),
            T2 = erlang:monotonic_time(),
            _ = viva_math@precision:pairwise_sum(Xs),
            T3 = erlang:monotonic_time(),
            gleam_stdlib:println(
                <<<<<<<<<<<<<<<<"  n="/utf8,
                                                (erlang:integer_to_binary(N))/binary>>/binary,
                                            "  naive="/utf8>>/binary,
                                        (erlang:integer_to_binary(T1 - T0))/binary>>/binary,
                                    "ns  neumaier="/utf8>>/binary,
                                (erlang:integer_to_binary(T2 - T1))/binary>>/binary,
                            "ns  pairwise="/utf8>>/binary,
                        (erlang:integer_to_binary(T3 - T2))/binary>>/binary,
                    "ns"/utf8>>
            )
        end
    ),
    gleam_stdlib:println(
        <<"\n[3] Pébay moments — variance + skew + kurtosis in one pass"/utf8>>
    ),
    Xs@1 = build_random(10000, 0.456),
    T0@1 = erlang:monotonic_time(),
    M = viva_math@precision:moments_from_list(Xs@1),
    T1@1 = erlang:monotonic_time(),
    _ = viva_math@precision:moments_variance(M),
    _ = viva_math@precision:moments_skewness(M),
    _ = viva_math@precision:moments_excess_kurtosis(M),
    gleam_stdlib:println(
        <<<<"  10k samples, single-pass Welford+Pébay → "/utf8,
                (erlang:integer_to_binary(T1@1 - T0@1))/binary>>/binary,
            "ns"/utf8>>
    ),
    gleam_stdlib:println(<<"\nDone.\n"/utf8>>).