src/lwes_mochinum.erl

   1 % This is the MIT license.
   2 %
   3 % Copyright (c) 2007 Mochi Media, Inc.
   4 %
   5 % Permission is hereby granted, free of charge, to any person obtaining a
   6 % copy of this software and associated documentation files (the "Software"),
   7 % to deal in the Software without restriction, including without limitation
   8 % the rights to use, copy, modify, merge, publish, distribute, sublicense,
   9 % and/or sell copies of the Software, and to permit persons to whom the
  10 % Software is furnished to do so, subject to the following conditions:
  11 %
  12 % The above copyright notice and this permission notice shall be included
  13 % in all copies or substantial portions of the Software.
  14 %
  15 % THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  16 % OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  17 % FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  18 % THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
  19 % OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
  20 % ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE
  21 % OR OTHER DEALINGS IN THE SOFTWARE.
  22 %
  23
  24 %% @copyright 2007 Mochi Media, Inc.
  25 %% @author Bob Ippolito <bob@mochimedia.com>
  26
  27 %% @doc Useful numeric algorithms for floats that cover some deficiencies
  28 %% in the math module. More interesting is digits/1, which implements
  29 %% the algorithm from:
  30 %% http://www.cs.indiana.edu/~burger/fp/index.html
  31 %% See also "Printing Floating-Point Numbers Quickly and Accurately"
  32 %% in Proceedings of the SIGPLAN '96 Conference on Programming Language
  33 %% Design and Implementation.
  34
  35 -module(lwes_mochinum).
  36 -author("Bob Ippolito <bob@mochimedia.com>").
  37 -export([digits/1, frexp/1, int_pow/2, int_ceil/1]).
  38
  39 %% IEEE 754 Float exponent bias
  40 -define(FLOAT_BIAS, 1022).
  41 -define(MIN_EXP, -1074).
  42 -define(BIG_POW, 4503599627370496).
  43
  44 %% External API
  45
  46 %% @spec digits(number()) -> string()
  47 %% @doc  Returns a string that accurately represents the given integer or float
  48 %%       using a conservative amount of digits. Great for generating
  49 %%       human-readable output, or compact ASCII serializations for floats.
  50 digits(N) when is_integer(N) ->
  51     integer_to_list(N);
  52 digits(0.0) ->
  53     "0.0";
  54 digits(Float) ->
  55     {Frac1, Exp1} = frexp_int(Float),
  56     [Place0 | Digits0] = digits1(Float, Exp1, Frac1),
  57     {Place, Digits} = transform_digits(Place0, Digits0),
  58     R = insert_decimal(Place, Digits),
  59     case Float < 0 of
  60         true ->
  61             [$- | R];
  62         _ ->
  63             R
  64     end.
  65
  66 %% @spec frexp(F::float()) -> {Frac::float(), Exp::float()}
  67 %% @doc  Return the fractional and exponent part of an IEEE 754 double,
  68 %%       equivalent to the libc function of the same name.
  69 %%       F = Frac * pow(2, Exp).
  70 frexp(F) ->
  71     frexp1(unpack(F)).
  72
  73 %% @spec int_pow(X::integer(), N::integer()) -> Y::integer()
  74 %% @doc  Moderately efficient way to exponentiate integers.
  75 %%       int_pow(10, 2) = 100.
  76 int_pow(_X, 0) ->
  77     1;
  78 int_pow(X, N) when N > 0 ->
  79     int_pow(X, N, 1).
  80
  81 %% @spec int_ceil(F::float()) -> integer()
  82 %% @doc  Return the ceiling of F as an integer. The ceiling is defined as
  83 %%       F when F == trunc(F);
  84 %%       trunc(F) when F &lt; 0;
  85 %%       trunc(F) + 1 when F &gt; 0.
  86 int_ceil(X) ->
  87     T = trunc(X),
  88     case (X - T) of
  89         Pos when Pos > 0 -> T + 1;
  90         _ -> T
  91     end.
  92
  93
  94 %% Internal API
  95
  96 int_pow(X, N, R) when N < 2 ->
  97     R * X;
  98 int_pow(X, N, R) ->
  99     int_pow(X * X, N bsr 1, case N band 1 of 1 -> R * X; 0 -> R end).
 100
 101 insert_decimal(0, S) ->
 102     "0." ++ S;
 103 insert_decimal(Place, S) when Place > 0 ->
 104     L = length(S),
 105     case Place - L of
 106          0 ->
 107             S ++ ".0";
 108         N when N < 0 ->
 109             {S0, S1} = lists:split(L + N, S),
 110             S0 ++ "." ++ S1;
 111         N when N < 6 ->
 112             %% More places than digits
 113             S ++ lists:duplicate(N, $0) ++ ".0";
 114         _ ->
 115             insert_decimal_exp(Place, S)
 116     end;
 117 insert_decimal(Place, S) when Place > -6 ->
 118     "0." ++ lists:duplicate(abs(Place), $0) ++ S;
 119 insert_decimal(Place, S) ->
 120     insert_decimal_exp(Place, S).
 121
 122 insert_decimal_exp(Place, S) ->
 123     [C | S0] = S,
 124     S1 = case S0 of
 125              [] ->
 126                  "0";
 127              _ ->
 128                  S0
 129          end,
 130     Exp = case Place < 0 of
 131               true ->
 132                   "e-";
 133               false ->
 134                   "e+"
 135           end,
 136     [C] ++ "." ++ S1 ++ Exp ++ integer_to_list(abs(Place - 1)).
 137
 138
 139 digits1(Float, Exp, Frac) ->
 140     Round = ((Frac band 1) =:= 0),
 141     case Exp >= 0 of
 142         true ->
 143             BExp = 1 bsl Exp,
 144             case (Frac =/= ?BIG_POW) of
 145                 true ->
 146                     scale((Frac * BExp * 2), 2, BExp, BExp,
 147                           Round, Round, Float);
 148                 false ->
 149                     scale((Frac * BExp * 4), 4, (BExp * 2), BExp,
 150                           Round, Round, Float)
 151             end;
 152         false ->
 153             case (Exp =:= ?MIN_EXP) orelse (Frac =/= ?BIG_POW) of
 154                 true ->
 155                     scale((Frac * 2), 1 bsl (1 - Exp), 1, 1,
 156                           Round, Round, Float);
 157                 false ->
 158                     scale((Frac * 4), 1 bsl (2 - Exp), 2, 1,
 159                           Round, Round, Float)
 160             end
 161     end.
 162
 163 scale(R, S, MPlus, MMinus, LowOk, HighOk, Float) ->
 164     Est = int_ceil(math:log10(abs(Float)) - 1.0e-10),
 165     %% Note that the scheme implementation uses a 326 element look-up table
 166     %% for int_pow(10, N) where we do not.
 167     case Est >= 0 of
 168         true ->
 169             fixup(R, S * int_pow(10, Est), MPlus, MMinus, Est,
 170                   LowOk, HighOk);
 171         false ->
 172             Scale = int_pow(10, -Est),
 173             fixup(R * Scale, S, MPlus * Scale, MMinus * Scale, Est,
 174                   LowOk, HighOk)
 175     end.
 176
 177 fixup(R, S, MPlus, MMinus, K, LowOk, HighOk) ->
 178     TooLow = case HighOk of
 179                  true ->
 180                      (R + MPlus) >= S;
 181                  false ->
 182                      (R + MPlus) > S
 183              end,
 184     case TooLow of
 185         true ->
 186             [(K + 1) | generate(R, S, MPlus, MMinus, LowOk, HighOk)];
 187         false ->
 188             [K | generate(R * 10, S, MPlus * 10, MMinus * 10, LowOk, HighOk)]
 189     end.
 190
 191 generate(R0, S, MPlus, MMinus, LowOk, HighOk) ->
 192     D = R0 div S,
 193     R = R0 rem S,
 194     TC1 = case LowOk of
 195               true ->
 196                   R =< MMinus;
 197               false ->
 198                   R < MMinus
 199           end,
 200     TC2 = case HighOk of
 201               true ->
 202                   (R + MPlus) >= S;
 203               false ->
 204                   (R + MPlus) > S
 205           end,
 206     case TC1 of
 207         false ->
 208             case TC2 of
 209                 false ->
 210                     [D | generate(R * 10, S, MPlus * 10, MMinus * 10,
 211                                   LowOk, HighOk)];
 212                 true ->
 213                     [D + 1]
 214             end;
 215         true ->
 216             case TC2 of
 217                 false ->
 218                     [D];
 219                 true ->
 220                     case R * 2 < S of
 221                         true ->
 222                             [D];
 223                         false ->
 224                             [D + 1]
 225                     end
 226             end
 227     end.
 228
 229 unpack(Float) ->
 230     <<Sign:1, Exp:11, Frac:52>> = <<Float:64/float>>,
 231     {Sign, Exp, Frac}.
 232
 233 frexp1({_Sign, 0, 0}) ->
 234     {0.0, 0};
 235 frexp1({Sign, 0, Frac}) ->
 236     Exp = log2floor(Frac),
 237     <<Frac1:64/float>> = <<Sign:1, ?FLOAT_BIAS:11, (Frac-1):52>>,
 238     {Frac1, -(?FLOAT_BIAS) - 52 + Exp};
 239 frexp1({Sign, Exp, Frac}) ->
 240     <<Frac1:64/float>> = <<Sign:1, ?FLOAT_BIAS:11, Frac:52>>,
 241     {Frac1, Exp - ?FLOAT_BIAS}.
 242
 243 log2floor(Int) ->
 244     log2floor(Int, 0).
 245
 246 log2floor(0, N) ->
 247     N;
 248 log2floor(Int, N) ->
 249     log2floor(Int bsr 1, 1 + N).
 250
 251
 252 transform_digits(Place, [0 | Rest]) ->
 253     transform_digits(Place, Rest);
 254 transform_digits(Place, Digits) ->
 255     {Place, [$0 + D || D <- Digits]}.
 256
 257
 258 frexp_int(F) ->
 259     case unpack(F) of
 260         {_Sign, 0, Frac} ->
 261             {Frac, ?MIN_EXP};
 262         {_Sign, Exp, Frac} ->
 263             {Frac + (1 bsl 52), Exp - 53 - ?FLOAT_BIAS}
 264     end.
 265
 266 %%
 267 %% Tests
 268 %%
 269 -ifdef(TEST).
 270 -include_lib("eunit/include/eunit.hrl").
 271
 272 int_ceil_test() ->
 273     ?assertEqual(1, int_ceil(0.0001)),
 274     ?assertEqual(0, int_ceil(0.0)),
 275     ?assertEqual(1, int_ceil(0.99)),
 276     ?assertEqual(1, int_ceil(1.0)),
 277     ?assertEqual(-1, int_ceil(-1.5)),
 278     ?assertEqual(-2, int_ceil(-2.0)),
 279     ok.
 280
 281 int_pow_test() ->
 282     ?assertEqual(1, int_pow(1, 1)),
 283     ?assertEqual(1, int_pow(1, 0)),
 284     ?assertEqual(1, int_pow(10, 0)),
 285     ?assertEqual(10, int_pow(10, 1)),
 286     ?assertEqual(100, int_pow(10, 2)),
 287     ?assertEqual(1000, int_pow(10, 3)),
 288     ok.
 289
 290 digits_test() ->
 291     ?assertEqual("0",
 292                  digits(0)),
 293     ?assertEqual("0.0",
 294                  digits(0.0)),
 295     ?assertEqual("1.0",
 296                  digits(1.0)),
 297     ?assertEqual("-1.0",
 298                  digits(-1.0)),
 299     ?assertEqual("0.1",
 300                  digits(0.1)),
 301     ?assertEqual("0.01",
 302                  digits(0.01)),
 303     ?assertEqual("0.001",
 304                  digits(0.001)),
 305     ?assertEqual("1.0e+6",
 306                  digits(1000000.0)),
 307     ?assertEqual("0.5",
 308                  digits(0.5)),
 309     ?assertEqual("4503599627370496.0",
 310                  digits(4503599627370496.0)),
 311     %% small denormalized number
 312     %% 4.94065645841246544177e-324 =:= 5.0e-324
 313     <<SmallDenorm/float>> = <<0,0,0,0,0,0,0,1>>,
 314     ?assertEqual("5.0e-324",
 315                  digits(SmallDenorm)),
 316     ?assertEqual(SmallDenorm,
 317                  list_to_float(digits(SmallDenorm))),
 318     %% large denormalized number
 319     %% 2.22507385850720088902e-308
 320     <<BigDenorm/float>> = <<0,15,255,255,255,255,255,255>>,
 321     ?assertEqual("2.225073858507201e-308",
 322                  digits(BigDenorm)),
 323     ?assertEqual(BigDenorm,
 324                  list_to_float(digits(BigDenorm))),
 325     %% small normalized number
 326     %% 2.22507385850720138309e-308
 327     <<SmallNorm/float>> = <<0,16,0,0,0,0,0,0>>,
 328     ?assertEqual("2.2250738585072014e-308",
 329                  digits(SmallNorm)),
 330     ?assertEqual(SmallNorm,
 331                  list_to_float(digits(SmallNorm))),
 332     %% large normalized number
 333     %% 1.79769313486231570815e+308
 334     <<LargeNorm/float>> = <<127,239,255,255,255,255,255,255>>,
 335     ?assertEqual("1.7976931348623157e+308",
 336                  digits(LargeNorm)),
 337     ?assertEqual(LargeNorm,
 338                  list_to_float(digits(LargeNorm))),
 339     %% issue #10 - mochinum:frexp(math:pow(2, -1074)).
 340     ?assertEqual("5.0e-324",
 341                  digits(math:pow(2, -1074))),
 342     ok.
 343
 344 frexp_test() ->
 345     %% zero
 346     ?assertEqual({0.0, 0}, frexp(0.0)),
 347     %% one
 348     ?assertEqual({0.5, 1}, frexp(1.0)),
 349     %% negative one
 350     ?assertEqual({-0.5, 1}, frexp(-1.0)),
 351     %% small denormalized number
 352     %% 4.94065645841246544177e-324
 353     <<SmallDenorm/float>> = <<0,0,0,0,0,0,0,1>>,
 354     ?assertEqual({0.5, -1073}, frexp(SmallDenorm)),
 355     %% large denormalized number
 356     %% 2.22507385850720088902e-308
 357     <<BigDenorm/float>> = <<0,15,255,255,255,255,255,255>>,
 358     ?assertEqual(
 359        {0.99999999999999978, -1022},
 360        frexp(BigDenorm)),
 361     %% small normalized number
 362     %% 2.22507385850720138309e-308
 363     <<SmallNorm/float>> = <<0,16,0,0,0,0,0,0>>,
 364     ?assertEqual({0.5, -1021}, frexp(SmallNorm)),
 365     %% large normalized number
 366     %% 1.79769313486231570815e+308
 367     <<LargeNorm/float>> = <<127,239,255,255,255,255,255,255>>,
 368     ?assertEqual(
 369         {0.99999999999999989, 1024},
 370         frexp(LargeNorm)),
 371     %% issue #10 - mochinum:frexp(math:pow(2, -1074)).
 372     ?assertEqual(
 373        {0.5, -1073},
 374        frexp(math:pow(2, -1074))),
 375     ok.
 376
 377 -endif.