PageRenderTime 206ms CodeModel.GetById 61ms app.highlight 79ms RepoModel.GetById 21ms app.codeStats 0ms

/src/mochinum.erl

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