Blame | Last modification | View Log | RSS feed
%% @copyright 2007 Mochi Media, Inc.
%% @author Bob Ippolito <bob@mochimedia.com>
%%
%% Permission is hereby granted, free of charge, to any person obtaining a
%% copy of this software and associated documentation files (the "Software"),
%% to deal in the Software without restriction, including without limitation
%% the rights to use, copy, modify, merge, publish, distribute, sublicense,
%% and/or sell copies of the Software, and to permit persons to whom the
%% Software is furnished to do so, subject to the following conditions:
%%
%% The above copyright notice and this permission notice shall be included in
%% all copies or substantial portions of the Software.
%%
%% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
%% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
%% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
%% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
%% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
%% FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
%% DEALINGS IN THE SOFTWARE.
%% @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(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::float()) -> {Frac::float(), Exp::float()}
%% @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 < 0;
%% trunc(F) + 1 when F > 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)),
ok.
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)),
ok.
digits_test() ->
?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)),
%% small denormalized number
%% 4.94065645841246544177e-324 =:= 5.0e-324
<<SmallDenorm/float>> = <<0,0,0,0,0,0,0,1>>,
?assertEqual("5.0e-324",
digits(SmallDenorm)),
?assertEqual(SmallDenorm,
list_to_float(digits(SmallDenorm))),
%% large denormalized number
%% 2.22507385850720088902e-308
<<BigDenorm/float>> = <<0,15,255,255,255,255,255,255>>,
?assertEqual("2.225073858507201e-308",
digits(BigDenorm)),
?assertEqual(BigDenorm,
list_to_float(digits(BigDenorm))),
%% small normalized number
%% 2.22507385850720138309e-308
<<SmallNorm/float>> = <<0,16,0,0,0,0,0,0>>,
?assertEqual("2.2250738585072014e-308",
digits(SmallNorm)),
?assertEqual(SmallNorm,
list_to_float(digits(SmallNorm))),
%% large normalized number
%% 1.79769313486231570815e+308
<<LargeNorm/float>> = <<127,239,255,255,255,255,255,255>>,
?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))),
ok.
frexp_test() ->
%% 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)),
%% small denormalized number
%% 4.94065645841246544177e-324
<<SmallDenorm/float>> = <<0,0,0,0,0,0,0,1>>,
?assertEqual({0.5, -1073}, frexp(SmallDenorm)),
%% large denormalized number
%% 2.22507385850720088902e-308
<<BigDenorm/float>> = <<0,15,255,255,255,255,255,255>>,
?assertEqual(
{0.99999999999999978, -1022},
frexp(BigDenorm)),
%% small normalized number
%% 2.22507385850720138309e-308
<<SmallNorm/float>> = <<0,16,0,0,0,0,0,0>>,
?assertEqual({0.5, -1021}, frexp(SmallNorm)),
%% large normalized number
%% 1.79769313486231570815e+308
<<LargeNorm/float>> = <<127,239,255,255,255,255,255,255>>,
?assertEqual(
{0.99999999999999989, 1024},
frexp(LargeNorm)),
%% issue #10 - mochinum:frexp(math:pow(2, -1074)).
?assertEqual(
{0.5, -1073},
frexp(math:pow(2, -1074))),
ok.
-endif.