%% @copyright 2007 Mochi Media, Inc. %% @author Bob Ippolito %% @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 "). -export([digits/1, frexp/1, int_pow/2, int_ceil/1, test/0]). %% 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) -> {Frac, Exp} = frexp(Float), Exp1 = Exp - 53, Frac1 = trunc(abs(Frac) * (1 bsl 53)), [Place | Digits] = digits1(Float, Exp1, Frac1), R = insert_decimal(Place, [$0 + D || D <- 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 Neg when Neg < 0 -> T; 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, Exp, Frac}. frexp1({_Sign, 0, 0}) -> {0.0, 0}; frexp1({Sign, 0, Frac}) -> Exp = log2floor(Frac), <> = <>, {Frac1, -(?FLOAT_BIAS) - 52 + Exp}; frexp1({Sign, Exp, Frac}) -> <> = <>, {Frac1, Exp - ?FLOAT_BIAS}. log2floor(Int) -> log2floor(Int, 0). log2floor(0, N) -> N; log2floor(Int, N) -> log2floor(Int bsr 1, 1 + N). test() -> ok = test_frexp(), ok = test_int_ceil(), ok = test_int_pow(), ok = test_digits(), ok. test_int_ceil() -> 1 = int_ceil(0.0001), 0 = int_ceil(0.0), 1 = int_ceil(0.99), 1 = int_ceil(1.0), -1 = int_ceil(-1.5), -2 = int_ceil(-2.0), ok. test_int_pow() -> 1 = int_pow(1, 1), 1 = int_pow(1, 0), 1 = int_pow(10, 0), 10 = int_pow(10, 1), 100 = int_pow(10, 2), 1000 = int_pow(10, 3), ok. test_digits() -> "0" = digits(0), "0.0" = digits(0.0), "1.0" = digits(1.0), "-1.0" = digits(-1.0), "0.1" = digits(0.1), "0.01" = digits(0.01), "0.001" = digits(0.001), ok. test_frexp() -> %% zero {0.0, 0} = frexp(0.0), %% one {0.5, 1} = frexp(1.0), %% negative one {-0.5, 1} = frexp(-1.0), %% small denormalized number %% 4.94065645841246544177e-324 <> = <<0,0,0,0,0,0,0,1>>, {0.5, -1073} = frexp(SmallDenorm), %% large denormalized number %% 2.22507385850720088902e-308 <> = <<0,15,255,255,255,255,255,255>>, {0.99999999999999978, -1022} = frexp(BigDenorm), %% small normalized number %% 2.22507385850720138309e-308 <> = <<0,16,0,0,0,0,0,0>>, {0.5, -1021} = frexp(SmallNorm), %% large normalized number %% 1.79769313486231570815e+308 <> = <<127,239,255,255,255,255,255,255>>, {0.99999999999999989, 1024} = frexp(LargeNorm), ok.