/deps/mochiweb/src/mochinum.erl

http://github.com/zotonic/zotonic · Erlang · 354 lines · 267 code · 32 blank · 55 comment · 1 complexity · 24017d3f66ddccb3aa2c9b72dd4025e4 MD5 · raw file

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