%%%-----------------------------------------------------------------------------
%%% @copyright 2007 Mochi Media, Inc.
%%% @doc Useful numeric algorithms for floats that cover some deficiencies
%%% in the math module.
%%%
%%% More interesting is digits/1, which implements the algorithm from:
%%% http://www.cs.indiana.edu/~burger/fp/index.html
%%% See also "Printing Floating-Point Numbers Quickly and Accurately"
%%% in Proceedings of the SIGPLAN '96 Conference on Programming Language
%%% Design and Implementation.
%%% module renamed due to conflicts with rabbitmq_common
%%%
%%% TODO: REMOVE THIS FILE WHEN {@link kz_util} IS FIXED!
%%%
%%% @author Bob Ippolito <bob@mochimedia.com>
%%% @end
%%%-----------------------------------------------------------------------------


-module(kz_mochinum).
-author("Bob Ippolito <bob@mochimedia.com>").
-export([digits/1, frexp/1, int_pow/2, int_ceil/1]).

%% IEEE 754 Float exponent bias
-define(FLOAT_BIAS, 1022).
-define(MIN_EXP, -1074).
-define(BIG_POW, 4503599627370496).

%% External API

-spec digits(number()) -> string().
%% @doc Returns a string that accurately represents the given integer or float
%%       using a conservative amount of digits. Great for generating
%%       human-readable output, or compact ASCII serializations for floats.
digits(N) when is_integer(N) ->
    integer_to_list(N);
digits(0.0) ->
    "0.0";
digits(Float) ->
    {Frac1, Exp1} = frexp_int(Float),
    [Place0 | Digits0] = digits1(Float, Exp1, Frac1),
    {Place, Digits} = transform_digits(Place0, Digits0),
    R = insert_decimal(Place, Digits),
    case Float < 0 of
        true ->
            [$- | R];
        _ ->
            R
    end.

-spec frexp(F::number()) -> {Frac::float(), Exp::integer()}.
%% @doc Return the fractional and exponent part of an IEEE 754 double,
%%       equivalent to the libc function of the same name.
%%       F = Frac * pow(2, Exp).
frexp(F) ->
    frexp1(unpack(F)).

-spec int_pow(X::integer(), N::integer()) -> Y::integer().
%% @doc Moderately efficient way to exponentiate integers.
%%       int_pow(10, 2) = 100.
int_pow(_X, 0) ->
    1;
int_pow(X, N) when N > 0 ->
    int_pow(X, N, 1).

-spec int_ceil(F::float()) -> integer().
%% @doc Return the ceiling of F as an integer. The ceiling is defined as
%%       F when F == trunc(F);
%%       trunc(F) when F &lt; 0;
%%       trunc(F) + 1 when F &gt; 0.
int_ceil(X) ->
    T = trunc(X),
    case (X - T) of
        Pos when Pos > 0 -> T + 1;
        _ -> T
    end.


%% Internal API

int_pow(X, N, R) when N < 2 ->
    R * X;
int_pow(X, N, R) ->
    int_pow(X * X, N bsr 1, case N band 1 of 1 -> R * X; 0 -> R end).

insert_decimal(0, S) ->
    "0." ++ S;
insert_decimal(Place, S) when Place > 0 ->
    L = length(S),
    case Place - L of
        0 ->
            S ++ ".0";
        N when N < 0 ->
            {S0, S1} = lists:split(L + N, S),
            S0 ++ "." ++ S1;
        N when N < 6 ->
            %% More places than digits
            S ++ lists:duplicate(N, $0) ++ ".0";
        _ ->
            insert_decimal_exp(Place, S)
    end;
insert_decimal(Place, S) when Place > -6 ->
    "0." ++ lists:duplicate(abs(Place), $0) ++ S;
insert_decimal(Place, S) ->
    insert_decimal_exp(Place, S).

insert_decimal_exp(Place, S) ->
    [C | S0] = S,
    S1 = case S0 of
             [] ->
                 "0";
             _ ->
                 S0
         end,
    Exp = case Place < 0 of
              true ->
                  "e-";
              false ->
                  "e+"
          end,
    [C] ++ "." ++ S1 ++ Exp ++ integer_to_list(abs(Place - 1)).


digits1(Float, Exp, Frac) ->
    Round = ((Frac band 1) =:= 0),
    case Exp >= 0 of
        true ->
            BExp = 1 bsl Exp,
            case (Frac =/= ?BIG_POW) of
                true ->
                    scale((Frac * BExp * 2), 2, BExp, BExp,
                          Round, Round, Float);
                false ->
                    scale((Frac * BExp * 4), 4, (BExp * 2), BExp,
                          Round, Round, Float)
            end;
        false ->
            case (Exp =:= ?MIN_EXP)
                orelse (Frac =/= ?BIG_POW)
            of
                true ->
                    scale((Frac * 2), 1 bsl (1 - Exp), 1, 1,
                          Round, Round, Float);
                false ->
                    scale((Frac * 4), 1 bsl (2 - Exp), 2, 1,
                          Round, Round, Float)
            end
    end.

scale(R, S, MPlus, MMinus, LowOk, HighOk, Float) ->
    Est = int_ceil(math:log10(abs(Float)) - 1.0e-10),
    %% Note that the scheme implementation uses a 326 element look-up table
    %% for int_pow(10, N) where we do not.
    case Est >= 0 of
        true ->
            fixup(R, S * int_pow(10, Est), MPlus, MMinus, Est,
                  LowOk, HighOk);
        false ->
            Scale = int_pow(10, -Est),
            fixup(R * Scale, S, MPlus * Scale, MMinus * Scale, Est,
                  LowOk, HighOk)
    end.

fixup(R, S, MPlus, MMinus, K, LowOk, HighOk) ->
    TooLow = case HighOk of
                 true ->
                     (R + MPlus) >= S;
                 false ->
                     (R + MPlus) > S
             end,
    case TooLow of
        true ->
            [(K + 1) | generate(R, S, MPlus, MMinus, LowOk, HighOk)];
        false ->
            [K | generate(R * 10, S, MPlus * 10, MMinus * 10, LowOk, HighOk)]
    end.

generate(R0, S, MPlus, MMinus, LowOk, HighOk) ->
    D = R0 div S,
    R = R0 rem S,
    TC1 = case LowOk of
              true ->
                  R =< MMinus;
              false ->
                  R < MMinus
          end,
    TC2 = case HighOk of
              true ->
                  (R + MPlus) >= S;
              false ->
                  (R + MPlus) > S
          end,
    case TC1 of
        false ->
            case TC2 of
                false ->
                    [D | generate(R * 10, S, MPlus * 10, MMinus * 10,
                                  LowOk, HighOk)];
                true ->
                    [D + 1]
            end;
        true ->
            case TC2 of
                false ->
                    [D];
                true ->
                    case R * 2 < S of
                        true ->
                            [D];
                        false ->
                            [D + 1]
                    end
            end
    end.

unpack(Float) ->
    <<Sign:1, Exp:11, Frac:52>> = <<Float:64/float>>,
    {Sign, Exp, Frac}.

frexp1({_Sign, 0, 0}) ->
    {0.0, 0};
frexp1({Sign, 0, Frac}) ->
    Exp = log2floor(Frac),
    <<Frac1:64/float>> = <<Sign:1, ?FLOAT_BIAS:11, (Frac-1):52>>,
    {Frac1, -(?FLOAT_BIAS) - 52 + Exp};
frexp1({Sign, Exp, Frac}) ->
    <<Frac1:64/float>> = <<Sign:1, ?FLOAT_BIAS:11, Frac:52>>,
    {Frac1, Exp - ?FLOAT_BIAS}.

log2floor(Int) ->
    log2floor(Int, 0).

log2floor(0, N) ->
    N;
log2floor(Int, N) ->
    log2floor(Int bsr 1, 1 + N).


transform_digits(Place, [0 | Rest]) ->
    transform_digits(Place, Rest);
transform_digits(Place, Digits) ->
    {Place, [$0 + D || D <- Digits]}.


frexp_int(F) ->
    case unpack(F) of
        {_Sign, 0, Frac} ->
            {Frac, ?MIN_EXP};
        {_Sign, Exp, Frac} ->
            {Frac + (1 bsl 52), Exp - 53 - ?FLOAT_BIAS}
    end.

%%
%% Tests
%%
-ifdef(TEST).
-include_lib("eunit/include/eunit.hrl").

int_ceil_test_() ->
    [?_assertEqual(1, int_ceil(0.0001))
    ,?_assertEqual(0, int_ceil(0.0))
    ,?_assertEqual(1, int_ceil(0.99))
    ,?_assertEqual(1, int_ceil(1.0))
    ,?_assertEqual(-1, int_ceil(-1.5))
    ,?_assertEqual(-2, int_ceil(-2.0))
    ].

int_pow_test_() ->
    [?_assertEqual(1, int_pow(1, 1))
    ,?_assertEqual(1, int_pow(1, 0))
    ,?_assertEqual(1, int_pow(10, 0))
    ,?_assertEqual(10, int_pow(10, 1))
    ,?_assertEqual(100, int_pow(10, 2))
    ,?_assertEqual(1000, int_pow(10, 3))
    ].

digits_test_() ->
    %% small de-normalized number
    %% 4.94065645841246544177e-324 =:= 5.0e-324
    <<SmallDenorm/float>> = <<0,0,0,0,0,0,0,1>>,

    %% large de-normalized number
    %% 2.22507385850720088902e-308
    <<BigDenorm/float>> = <<0,15,255,255,255,255,255,255>>,

    %% small normalized number
    %% 2.22507385850720138309e-308
    <<SmallNorm/float>> = <<0,16,0,0,0,0,0,0>>,

    %% large normalized number
    %% 1.79769313486231570815e+308
    <<LargeNorm/float>> = <<127,239,255,255,255,255,255,255>>,

    [?_assertEqual("0",
                   digits(0))
    ,?_assertEqual("0.0",
                   digits(0.0))
    ,?_assertEqual("1.0",
                   digits(1.0))
    ,?_assertEqual("-1.0",
                   digits(-1.0))
    ,?_assertEqual("0.1",
                   digits(0.1))
    ,?_assertEqual("0.01",
                   digits(0.01))
    ,?_assertEqual("0.001",
                   digits(0.001))
    ,?_assertEqual("1.0e+6",
                   digits(1000000.0))
    ,?_assertEqual("0.5",
                   digits(0.5))
    ,?_assertEqual("4503599627370496.0",
                   digits(4503599627370496.0))
    ,?_assertEqual("5.0e-324",
                   digits(SmallDenorm))
    ,?_assertEqual(SmallDenorm,
                   list_to_float(digits(SmallDenorm)))
    ,?_assertEqual("2.225073858507201e-308",
                   digits(BigDenorm))
    ,?_assertEqual(BigDenorm,
                   list_to_float(digits(BigDenorm)))
    ,?_assertEqual("2.2250738585072014e-308",
                   digits(SmallNorm))
    ,?_assertEqual(SmallNorm,
                   list_to_float(digits(SmallNorm)))
    ,?_assertEqual("1.7976931348623157e+308",
                   digits(LargeNorm))
    ,?_assertEqual(LargeNorm,
                   list_to_float(digits(LargeNorm)))
     %% issue #10 - mochinum:frexp(math:pow(2, -1074)).
    ,?_assertEqual("5.0e-324",
                   digits(math:pow(2, -1074)))
    ].

frexp_test_() ->
    %% small de-normalized number
    %% 4.94065645841246544177e-324
    <<SmallDenorm/float>> = <<0,0,0,0,0,0,0,1>>,

    %% large de-normalized number
    %% 2.22507385850720088902e-308
    <<BigDenorm/float>> = <<0,15,255,255,255,255,255,255>>,

    %% small normalized number
    %% 2.22507385850720138309e-308
    <<SmallNorm/float>> = <<0,16,0,0,0,0,0,0>>,

    %% large normalized number
    %% 1.79769313486231570815e+308
    <<LargeNorm/float>> = <<127,239,255,255,255,255,255,255>>,

    %% zero
    [?_assertEqual({0.0, 0}, frexp(0.0))
     %% one
    ,?_assertEqual({0.5, 1}, frexp(1.0))
     %% negative one
    ,?_assertEqual({-0.5, 1}, frexp(-1.0))
    ,?_assertEqual({0.5, -1073}, frexp(SmallDenorm))
    ,?_assertEqual({0.99999999999999978, -1022}, frexp(BigDenorm))
    ,?_assertEqual({0.5, -1021}, frexp(SmallNorm))
    ,?_assertEqual({0.99999999999999989, 1024}, frexp(LargeNorm))
     %% issue #10 - mochinum:frexp(math:pow(2, -1074)).
    ,?_assertEqual({0.5, -1073}, frexp(math:pow(2, -1074)))
    ].

-endif.