| %% @copyright 2007 Mochi Media, Inc. |
| %% @author Bob Ippolito <bob@mochimedia.com> |
| |
| %% @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, 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: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). |
| |
| |
| 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 |
| <<SmallDenorm/float>> = <<0,0,0,0,0,0,0,1>>, |
| {0.5, -1073} = frexp(SmallDenorm), |
| %% large denormalized number |
| %% 2.22507385850720088902e-308 |
| <<BigDenorm/float>> = <<0,15,255,255,255,255,255,255>>, |
| {0.99999999999999978, -1022} = frexp(BigDenorm), |
| %% small normalized number |
| %% 2.22507385850720138309e-308 |
| <<SmallNorm/float>> = <<0,16,0,0,0,0,0,0>>, |
| {0.5, -1021} = frexp(SmallNorm), |
| %% large normalized number |
| %% 1.79769313486231570815e+308 |
| <<LargeNorm/float>> = <<127,239,255,255,255,255,255,255>>, |
| {0.99999999999999989, 1024} = frexp(LargeNorm), |
| ok. |
