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 < 0; 62%% trunc(F) + 1 when F > 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. 355