/deps/mochiweb/src/mochinum.erl

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