From 0482a605f5f6782db77c7b8d334b9cb66f89d8c7 Mon Sep 17 00:00:00 2001 From: Marshall Roch Date: Mon, 3 Apr 2017 15:06:40 -0700 Subject: [PATCH] initial commit --- .gitignore | 8 + .merlin | 4 + .ocp-indent | 1 + CHANGES.md | 4 + LICENSE | 61 ++++++ PATENTS | 33 ++++ README.md | 29 +++ _tags | 4 + doc/api.odocl | 1 + doc/dev.odocl | 1 + opam | 22 +++ pkg/META | 7 + pkg/pkg.ml | 11 ++ src/cached_powers.c | 135 +++++++++++++ src/cached_powers.h | 22 +++ src/diy_fp.c | 72 +++++++ src/diy_fp.h | 33 ++++ src/dtoa.ml | 12 ++ src/dtoa.mli | 27 +++ src/dtoa.mllib | 1 + src/dtoa_stubs.c | 287 ++++++++++++++++++++++++++++ src/fast_dtoa.c | 425 +++++++++++++++++++++++++++++++++++++++++ src/fast_dtoa.h | 61 ++++++ src/ieee.c | 143 ++++++++++++++ src/ieee.h | 43 +++++ src/utils.h | 37 ++++ test/test.ml | 21 ++ tests/ecma_test.ml | 63 ++++++ tests/g_fmt_test.ml | 63 ++++++ tests/shortest_test.ml | 67 +++++++ tests/test.ml | 18 ++ 31 files changed, 1716 insertions(+) create mode 100644 .gitignore create mode 100644 .merlin create mode 100644 .ocp-indent create mode 100644 CHANGES.md create mode 100644 LICENSE create mode 100644 PATENTS create mode 100644 README.md create mode 100644 _tags create mode 100644 doc/api.odocl create mode 100644 doc/dev.odocl create mode 100644 opam create mode 100644 pkg/META create mode 100755 pkg/pkg.ml create mode 100644 src/cached_powers.c create mode 100644 src/cached_powers.h create mode 100644 src/diy_fp.c create mode 100644 src/diy_fp.h create mode 100644 src/dtoa.ml create mode 100644 src/dtoa.mli create mode 100644 src/dtoa.mllib create mode 100644 src/dtoa_stubs.c create mode 100644 src/fast_dtoa.c create mode 100644 src/fast_dtoa.h create mode 100644 src/ieee.c create mode 100644 src/ieee.h create mode 100644 src/utils.h create mode 100644 test/test.ml create mode 100644 tests/ecma_test.ml create mode 100644 tests/g_fmt_test.ml create mode 100644 tests/shortest_test.ml create mode 100644 tests/test.ml diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b48ed60 --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +_build +tmp +*~ +\.\#* +\#*# +*.install +*.native +*.byte \ No newline at end of file diff --git a/.merlin b/.merlin new file mode 100644 index 0000000..0e0495c --- /dev/null +++ b/.merlin @@ -0,0 +1,4 @@ +PKG bytes +S src +S test +B _build/** diff --git a/.ocp-indent b/.ocp-indent new file mode 100644 index 0000000..6c5b8f4 --- /dev/null +++ b/.ocp-indent @@ -0,0 +1 @@ +strict_with=always,match_clause=4,strict_else=never diff --git a/CHANGES.md b/CHANGES.md new file mode 100644 index 0000000..ca9a474 --- /dev/null +++ b/CHANGES.md @@ -0,0 +1,4 @@ +v0.1.0 2017-04-03 +-------------------------- + +First release. diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..d181d8d --- /dev/null +++ b/LICENSE @@ -0,0 +1,61 @@ +BSD License + +For ocaml-dtoa software + +Copyright (c) 2017-present, Facebook, Inc. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + * Neither the name Facebook nor the names of its contributors may be used to + endorse or promote products derived from this software without specific + prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +---- + +Based on double-conversion (https://github.com/google/double-conversion): + +Copyright 2006-2011, the V8 project authors. All rights reserved. +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + * Neither the name of Google Inc. nor the names of its + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/PATENTS b/PATENTS new file mode 100644 index 0000000..68b481e --- /dev/null +++ b/PATENTS @@ -0,0 +1,33 @@ +Additional Grant of Patent Rights Version 2 + +"Software" means the ocaml-dtoa software contributed by Facebook, Inc. + +Facebook, Inc. ("Facebook") hereby grants to each recipient of the Software +("you") a perpetual, worldwide, royalty-free, non-exclusive, irrevocable +(subject to the termination provision below) license under any Necessary +Claims, to make, have made, use, sell, offer to sell, import, and otherwise +transfer the Software. For avoidance of doubt, no license is granted under +Facebook’s rights in any patent claims that are infringed by (i) modifications +to the Software made by you or any third party or (ii) the Software in +combination with any software or other technology. + +The license granted hereunder will terminate, automatically and without notice, +if you (or any of your subsidiaries, corporate affiliates or agents) initiate +directly or indirectly, or take a direct financial interest in, any Patent +Assertion: (i) against Facebook or any of its subsidiaries or corporate +affiliates, (ii) against any party if such Patent Assertion arises in whole or +in part from any software, technology, product or service of Facebook or any of +its subsidiaries or corporate affiliates, or (iii) against any party relating +to the Software. Notwithstanding the foregoing, if Facebook or any of its +subsidiaries or corporate affiliates files a lawsuit alleging patent +infringement against you in the first instance, and you respond by filing a +patent infringement counterclaim in that lawsuit against that party that is +unrelated to the Software, the license granted hereunder will not terminate +under section (i) of this paragraph due to such counterclaim. + +A "Necessary Claim" is a claim of a patent owned by Facebook that is +necessarily infringed by the Software standing alone. + +A "Patent Assertion" is any lawsuit or other action alleging direct, indirect, +or contributory infringement or inducement to infringe any patent, including a +cross-claim or counterclaim. diff --git a/README.md b/README.md new file mode 100644 index 0000000..4f982a3 --- /dev/null +++ b/README.md @@ -0,0 +1,29 @@ +# ocaml-dtoa + +This library provides a function that converts OCaml floats into strings, using the efficient Grisu3 algorithm. + +The Grisu3 algorithm is described in ["Printing Floating-Point Numbers Quickly And Accurately with Integers"](http://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf) by Florian Loitsch. + +The implementation is adapted from [double-conversion](https://github.com/google/double-conversion). + +## Current Status + +Currently, this library exposes three functions: + +- `ecma_string_of_float : float -> string`: formats the float according to the ECMAScript specification's implementation of [Number.prototype.toString](https://tc39.github.io/ecma262/#sec-tostring-applied-to-the-number-type). Notably, the output of this function is valid JSON. + +- `shortest_string_of_float : float -> string`: formats the float as compactly as possible, for example returning `123e6` instead of `123000000` or `1.23e+08`. + +- `g_fmt : float -> string`: formats the float in the same way as David M. Gay's [`g_fmt`](http://www.netlib.org/fp/g_fmt.c). + +The underlying [`fast_dtoa()`](src/fast_dtoa.h) function computes the significand and exponent, which are formatted by the above functions in [`dtoa_stubs.c`](src/dtoa_stubs.c). It is a port of the [`double-conversion`](https://github.com/google/double-conversion) library from C++ to C. + +Many features of `double-conversion` are still missing. Patches are welcome! + +## License + +`ocaml-dtoa` is [BSD-licensed](LICENSE). We also provide an additional [patent grant](PATENTS). + +## Author + +`ocaml-dtoa` was created by Facebook for the [Flow](https://flow.org) project. diff --git a/_tags b/_tags new file mode 100644 index 0000000..a0d3945 --- /dev/null +++ b/_tags @@ -0,0 +1,4 @@ +true : bin_annot, safe_string, package(bytes) + + : include + : include \ No newline at end of file diff --git a/doc/api.odocl b/doc/api.odocl new file mode 100644 index 0000000..ff024d2 --- /dev/null +++ b/doc/api.odocl @@ -0,0 +1 @@ +Dtoa \ No newline at end of file diff --git a/doc/dev.odocl b/doc/dev.odocl new file mode 100644 index 0000000..ff024d2 --- /dev/null +++ b/doc/dev.odocl @@ -0,0 +1 @@ +Dtoa \ No newline at end of file diff --git a/opam b/opam new file mode 100644 index 0000000..e19e5ff --- /dev/null +++ b/opam @@ -0,0 +1,22 @@ +opam-version: "1.2" +maintainer: "Marshall Roch " +authors: ["Marshall Roch "] +homepage: "https://github.com/flowtype/ocaml-dtoa" +doc: "https://github.com/flowtype/ocaml-dtoa" +license: "BSD" +dev-repo: "https://github.com/flowtype/ocaml-dtoa.git" +bug-reports: "https://github.com/flowtype/ocaml-dtoa/issues" +tags: [] +available: [ ocaml-version >= "4.01.0"] +depends: +[ + "ocamlfind" {build} + "ocamlbuild" {build} + "topkg" {build & >= "0.9.0"} +] +depopts: [] +build: +[[ + "ocaml" "pkg/pkg.ml" "build" + "--pinned" "%{pinned}%" +]] diff --git a/pkg/META b/pkg/META new file mode 100644 index 0000000..f5a8290 --- /dev/null +++ b/pkg/META @@ -0,0 +1,7 @@ +description = "converts OCaml floats into strings, using the efficient Grisu3 algorithm." +version = "%%VERSION_NUM%%" +requires = "" +archive(byte) = "dtoa.cma" +archive(native) = "dtoa.cmxa" +plugin(byte) = "dtoa.cma" +plugin(native) = "dtoa.cmxs" diff --git a/pkg/pkg.ml b/pkg/pkg.ml new file mode 100755 index 0000000..b7f3173 --- /dev/null +++ b/pkg/pkg.ml @@ -0,0 +1,11 @@ +#!/usr/bin/env ocaml +#use "topfind" +#require "topkg" +open Topkg + +let () = + Pkg.describe "dtoa" + ~licenses:[ Pkg.std_file "LICENSE"; Pkg.std_file "PATENTS" ] + @@ fun c -> + Ok [ Pkg.mllib "src/dtoa.mllib"; + Pkg.test "test/test"; ] diff --git a/src/cached_powers.c b/src/cached_powers.c new file mode 100644 index 0000000..d7efda1 --- /dev/null +++ b/src/cached_powers.c @@ -0,0 +1,135 @@ +/** + * Copyright (c) 2017-present, Facebook, Inc. + * Copyright (c) 2010, the V8 project authors. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + */ + +#include +#include "cached_powers.h" +#include "utils.h" + +static const int kCachedPowersOffset = 348; // -1 * the first decimal_exponent. +static const double kD_1_LOG2_10 = 0.30102999566398114; // 1 / lg(10) +// Difference between the decimal exponents in the table above. +const int kDecimalExponentDistance = 8; +const int kMinDecimalExponent = -348; +const int kMaxDecimalExponent = 340; + +typedef struct { + uint64_t significand; + int16_t binary_exponent; + int16_t decimal_exponent; +} CachedPower; + +static const CachedPower kCachedPowers[] = { + {UINT64_2PART_C(0xfa8fd5a0, 081c0288), -1220, -348}, + {UINT64_2PART_C(0xbaaee17f, a23ebf76), -1193, -340}, + {UINT64_2PART_C(0x8b16fb20, 3055ac76), -1166, -332}, + {UINT64_2PART_C(0xcf42894a, 5dce35ea), -1140, -324}, + {UINT64_2PART_C(0x9a6bb0aa, 55653b2d), -1113, -316}, + {UINT64_2PART_C(0xe61acf03, 3d1a45df), -1087, -308}, + {UINT64_2PART_C(0xab70fe17, c79ac6ca), -1060, -300}, + {UINT64_2PART_C(0xff77b1fc, bebcdc4f), -1034, -292}, + {UINT64_2PART_C(0xbe5691ef, 416bd60c), -1007, -284}, + {UINT64_2PART_C(0x8dd01fad, 907ffc3c), -980, -276}, + {UINT64_2PART_C(0xd3515c28, 31559a83), -954, -268}, + {UINT64_2PART_C(0x9d71ac8f, ada6c9b5), -927, -260}, + {UINT64_2PART_C(0xea9c2277, 23ee8bcb), -901, -252}, + {UINT64_2PART_C(0xaecc4991, 4078536d), -874, -244}, + {UINT64_2PART_C(0x823c1279, 5db6ce57), -847, -236}, + {UINT64_2PART_C(0xc2109436, 4dfb5637), -821, -228}, + {UINT64_2PART_C(0x9096ea6f, 3848984f), -794, -220}, + {UINT64_2PART_C(0xd77485cb, 25823ac7), -768, -212}, + {UINT64_2PART_C(0xa086cfcd, 97bf97f4), -741, -204}, + {UINT64_2PART_C(0xef340a98, 172aace5), -715, -196}, + {UINT64_2PART_C(0xb23867fb, 2a35b28e), -688, -188}, + {UINT64_2PART_C(0x84c8d4df, d2c63f3b), -661, -180}, + {UINT64_2PART_C(0xc5dd4427, 1ad3cdba), -635, -172}, + {UINT64_2PART_C(0x936b9fce, bb25c996), -608, -164}, + {UINT64_2PART_C(0xdbac6c24, 7d62a584), -582, -156}, + {UINT64_2PART_C(0xa3ab6658, 0d5fdaf6), -555, -148}, + {UINT64_2PART_C(0xf3e2f893, dec3f126), -529, -140}, + {UINT64_2PART_C(0xb5b5ada8, aaff80b8), -502, -132}, + {UINT64_2PART_C(0x87625f05, 6c7c4a8b), -475, -124}, + {UINT64_2PART_C(0xc9bcff60, 34c13053), -449, -116}, + {UINT64_2PART_C(0x964e858c, 91ba2655), -422, -108}, + {UINT64_2PART_C(0xdff97724, 70297ebd), -396, -100}, + {UINT64_2PART_C(0xa6dfbd9f, b8e5b88f), -369, -92}, + {UINT64_2PART_C(0xf8a95fcf, 88747d94), -343, -84}, + {UINT64_2PART_C(0xb9447093, 8fa89bcf), -316, -76}, + {UINT64_2PART_C(0x8a08f0f8, bf0f156b), -289, -68}, + {UINT64_2PART_C(0xcdb02555, 653131b6), -263, -60}, + {UINT64_2PART_C(0x993fe2c6, d07b7fac), -236, -52}, + {UINT64_2PART_C(0xe45c10c4, 2a2b3b06), -210, -44}, + {UINT64_2PART_C(0xaa242499, 697392d3), -183, -36}, + {UINT64_2PART_C(0xfd87b5f2, 8300ca0e), -157, -28}, + {UINT64_2PART_C(0xbce50864, 92111aeb), -130, -20}, + {UINT64_2PART_C(0x8cbccc09, 6f5088cc), -103, -12}, + {UINT64_2PART_C(0xd1b71758, e219652c), -77, -4}, + {UINT64_2PART_C(0x9c400000, 00000000), -50, 4}, + {UINT64_2PART_C(0xe8d4a510, 00000000), -24, 12}, + {UINT64_2PART_C(0xad78ebc5, ac620000), 3, 20}, + {UINT64_2PART_C(0x813f3978, f8940984), 30, 28}, + {UINT64_2PART_C(0xc097ce7b, c90715b3), 56, 36}, + {UINT64_2PART_C(0x8f7e32ce, 7bea5c70), 83, 44}, + {UINT64_2PART_C(0xd5d238a4, abe98068), 109, 52}, + {UINT64_2PART_C(0x9f4f2726, 179a2245), 136, 60}, + {UINT64_2PART_C(0xed63a231, d4c4fb27), 162, 68}, + {UINT64_2PART_C(0xb0de6538, 8cc8ada8), 189, 76}, + {UINT64_2PART_C(0x83c7088e, 1aab65db), 216, 84}, + {UINT64_2PART_C(0xc45d1df9, 42711d9a), 242, 92}, + {UINT64_2PART_C(0x924d692c, a61be758), 269, 100}, + {UINT64_2PART_C(0xda01ee64, 1a708dea), 295, 108}, + {UINT64_2PART_C(0xa26da399, 9aef774a), 322, 116}, + {UINT64_2PART_C(0xf209787b, b47d6b85), 348, 124}, + {UINT64_2PART_C(0xb454e4a1, 79dd1877), 375, 132}, + {UINT64_2PART_C(0x865b8692, 5b9bc5c2), 402, 140}, + {UINT64_2PART_C(0xc83553c5, c8965d3d), 428, 148}, + {UINT64_2PART_C(0x952ab45c, fa97a0b3), 455, 156}, + {UINT64_2PART_C(0xde469fbd, 99a05fe3), 481, 164}, + {UINT64_2PART_C(0xa59bc234, db398c25), 508, 172}, + {UINT64_2PART_C(0xf6c69a72, a3989f5c), 534, 180}, + {UINT64_2PART_C(0xb7dcbf53, 54e9bece), 561, 188}, + {UINT64_2PART_C(0x88fcf317, f22241e2), 588, 196}, + {UINT64_2PART_C(0xcc20ce9b, d35c78a5), 614, 204}, + {UINT64_2PART_C(0x98165af3, 7b2153df), 641, 212}, + {UINT64_2PART_C(0xe2a0b5dc, 971f303a), 667, 220}, + {UINT64_2PART_C(0xa8d9d153, 5ce3b396), 694, 228}, + {UINT64_2PART_C(0xfb9b7cd9, a4a7443c), 720, 236}, + {UINT64_2PART_C(0xbb764c4c, a7a44410), 747, 244}, + {UINT64_2PART_C(0x8bab8eef, b6409c1a), 774, 252}, + {UINT64_2PART_C(0xd01fef10, a657842c), 800, 260}, + {UINT64_2PART_C(0x9b10a4e5, e9913129), 827, 268}, + {UINT64_2PART_C(0xe7109bfb, a19c0c9d), 853, 276}, + {UINT64_2PART_C(0xac2820d9, 623bf429), 880, 284}, + {UINT64_2PART_C(0x80444b5e, 7aa7cf85), 907, 292}, + {UINT64_2PART_C(0xbf21e440, 03acdd2d), 933, 300}, + {UINT64_2PART_C(0x8e679c2f, 5e44ff8f), 960, 308}, + {UINT64_2PART_C(0xd433179d, 9c8cb841), 986, 316}, + {UINT64_2PART_C(0x9e19db92, b4e31ba9), 1013, 324}, + {UINT64_2PART_C(0xeb96bf6e, badf77d9), 1039, 332}, + {UINT64_2PART_C(0xaf87023b, 9bf0ee6b), 1066, 340}, +}; + +void cached_power_for_binary_exponent_range( + int min_exponent, + int max_exponent, + diy_fp* power, + int* decimal_exponent) { + int kQ = kDiyFpSignificandSize; + double k = ceil((min_exponent + kQ - 1) * kD_1_LOG2_10); + int foo = kCachedPowersOffset; + int index = (foo + (int)(k) - 1) / kDecimalExponentDistance + 1; + ASSERT(0 <= index && index < (int)(ARRAY_SIZE(kCachedPowers))); + CachedPower cached_power = kCachedPowers[index]; + ASSERT(min_exponent <= cached_power.binary_exponent); + (void) max_exponent; // Mark variable as used. + ASSERT(cached_power.binary_exponent <= max_exponent); + diy_fp r = { cached_power.significand, cached_power.binary_exponent }; + *decimal_exponent = cached_power.decimal_exponent; + *power = r; +} diff --git a/src/cached_powers.h b/src/cached_powers.h new file mode 100644 index 0000000..844f54b --- /dev/null +++ b/src/cached_powers.h @@ -0,0 +1,22 @@ +/** + * Copyright (c) 2017-present, Facebook, Inc. + * Copyright (c) 2010, the V8 project authors. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + */ + +#ifndef OCAML_DTOA_CACHED_POWERS_H_ +#define OCAML_DTOA_CACHED_POWERS_H_ + +#include "diy_fp.h" + +void cached_power_for_binary_exponent_range( + int min_exponent, + int max_exponent, + diy_fp* power, + int* decimal_exponent); + +#endif // OCAML_DTOA_CACHED_POWERS_H_ diff --git a/src/diy_fp.c b/src/diy_fp.c new file mode 100644 index 0000000..bf97473 --- /dev/null +++ b/src/diy_fp.c @@ -0,0 +1,72 @@ +/** + * Copyright (c) 2017-present, Facebook, Inc. + * Copyright (c) 2010, the V8 project authors. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + */ + +#include + +#include "diy_fp.h" +#include "utils.h" + +const int kDiyFpSignificandSize = 64; +static const uint64_t kUint64MSB = UINT64_2PART_C(0x80000000, 00000000); + +diy_fp diy_fp_minus(diy_fp a, diy_fp b) { + diy_fp r; + assert(a.e == b.e); + assert(a.f >= b.f); + r.f = a.f - b.f; + r.e = a.e; + return r; +} + +diy_fp diy_fp_multiply(diy_fp x, diy_fp y) { + // Simply "emulates" a 128 bit multiplication. + // However: the resulting number only contains 64 bits. The least + // significant 64 bits are only used for rounding the most significant 64 + // bits. + diy_fp r; + const uint64_t kM32 = 0xFFFFFFFFU; + uint64_t a = x.f >> 32; + uint64_t b = x.f & kM32; + uint64_t c = y.f >> 32; + uint64_t d = y.f & kM32; + uint64_t ac = a * c; + uint64_t bc = b * c; + uint64_t ad = a * d; + uint64_t bd = b * d; + uint64_t tmp = (bd >> 32) + (ad & kM32) + (bc & kM32); + // By adding 1U << 31 to tmp we round the final result. + // Halfway cases will be round up. + tmp += 1U << 31; + r.f = ac + (ad >> 32) + (bc >> 32) + (tmp >> 32); + r.e = x.e + y.e + 64; + return r; +} + +diy_fp diy_fp_normalize(diy_fp x) { + diy_fp r; + assert(x.f != 0); + uint64_t significand = x.f; + int exponent = x.e; + + // This method is mainly called for normalizing boundaries. In general + // boundaries need to be shifted by 10 bits. We thus optimize for this case. + const uint64_t k10MSBits = UINT64_2PART_C(0xFFC00000, 00000000); + while ((significand & k10MSBits) == 0) { + significand <<= 10; + exponent -= 10; + } + while ((significand & kUint64MSB) == 0) { + significand <<= 1; + exponent--; + } + r.f = significand; + r.e = exponent; + return r; +} diff --git a/src/diy_fp.h b/src/diy_fp.h new file mode 100644 index 0000000..d52f0f8 --- /dev/null +++ b/src/diy_fp.h @@ -0,0 +1,33 @@ +/** + * Copyright (c) 2017-present, Facebook, Inc. + * Copyright (c) 2010, the V8 project authors. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + */ + +#ifndef OCAML_DTOA_DIY_FP_H_ +#define OCAML_DTOA_DIY_FP_H_ + +#include + +typedef struct diy_fp { + uint64_t f; + int e; +} diy_fp; + +extern const int kDiyFpSignificandSize; + +// Returns a - b. +// The exponents of both numbers must be the same and this must be bigger +// than other. The result will not be normalized. +diy_fp diy_fp_minus(diy_fp a, diy_fp b); + +// returns a * b; +diy_fp diy_fp_multiply(diy_fp a, diy_fp b); + +diy_fp diy_fp_normalize(diy_fp x); + +#endif // OCAML_DTOA_DIY_FP_H_ diff --git a/src/dtoa.ml b/src/dtoa.ml new file mode 100644 index 0000000..1483a56 --- /dev/null +++ b/src/dtoa.ml @@ -0,0 +1,12 @@ +(** + * Copyright (c) 2017-present, Facebook, Inc. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + *) + +external shortest_string_of_float: float -> string = "flow_shortest_string_of_float" +external ecma_string_of_float: float -> string = "flow_ecma_string_of_float" +external g_fmt: float -> string = "flow_g_fmt" diff --git a/src/dtoa.mli b/src/dtoa.mli new file mode 100644 index 0000000..73a5578 --- /dev/null +++ b/src/dtoa.mli @@ -0,0 +1,27 @@ +(*--------------------------------------------------------------------------- + Copyright (c) 2017 Marshall Roch. All rights reserved. + Distributed under the ISC license, see terms at the end of the file. + %%NAME%% %%VERSION%% + ---------------------------------------------------------------------------*) + +(** converts OCaml floats into strings, using the efficient Grisu3 algorithm. + + {e %%VERSION%% — {{:%%PKG_HOMEPAGE%% }homepage}} *) + +(** {1 Dtoa} *) + +(*--------------------------------------------------------------------------- + Copyright (c) 2017 Marshall Roch + + Permission to use, copy, modify, and/or distribute this software for any + purpose with or without fee is hereby granted, provided that the above + copyright notice and this permission notice appear in all copies. + + THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + ---------------------------------------------------------------------------*) diff --git a/src/dtoa.mllib b/src/dtoa.mllib new file mode 100644 index 0000000..dea0a3f --- /dev/null +++ b/src/dtoa.mllib @@ -0,0 +1 @@ +Dtoa diff --git a/src/dtoa_stubs.c b/src/dtoa_stubs.c new file mode 100644 index 0000000..eac8abd --- /dev/null +++ b/src/dtoa_stubs.c @@ -0,0 +1,287 @@ +/** + * Copyright (c) 2017-present, Facebook, Inc. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + */ + +#include +#include +#include +#include +#include +#include +#include +#include "fast_dtoa.h" + +typedef enum { + NO_FLAGS = 0, + LEADING_ZERO = 1, // 0.123 instead of .123 + PLUS_IN_EXPONENT = 2, // 1e+3 instead of 1e3 + PAD_EXPONENT = 4, // 1e+03 instead of 1e+3 +} flags; + +static int handle_special(double v, char *dst) { + int fp_class = fpclassify(v); + if (fp_class == FP_NAN) { + return snprintf(dst, 4, "NaN"); + } + if (fp_class == FP_ZERO) { + return snprintf(dst, 2, "0"); + } + if (fp_class == FP_INFINITE) { + if (v < 0) { + return snprintf(dst, 10, "-Infinity"); + } else { + return snprintf(dst, 9, "Infinity"); + } + } + return -1; // not special +} + +static int i_to_str(int val, char *str) { + int len, i; + char *s; + char *begin = str; + if (val < 0) { *str++ = '-'; val = -val; } + s = str; + + for(;;) + { + int ni = val / 10; + int digit = val - ni*10; + *s++ = (char)('0' + digit); + if (ni == 0) + break; + val = ni; + } + *s = '\0'; + len = (int)(s - str); + for(i = 0; i < len/2; ++i) + { + char ch = str[i]; + str[i] = str[len-1-i]; + str[len-1-i] = ch; + } + + return (int)(s - begin); +} + +static int decimal(char *dst, int len, int decimal_point, flags flags) { + int leading_zero = flags & LEADING_ZERO; + int d_exp = decimal_point - len, written = 0, shift = 0, j; + if (decimal_point <= 0) { + if (leading_zero) shift++; + shift++; // decimal point + // shift digits over, then fill in zeros and decimal point + memmove(dst + shift - decimal_point, dst, len); + if (leading_zero) dst[0] = '0'; + dst[leading_zero ? 1 : 0] = '.'; + for (j = 0; j < -decimal_point; ++j) { + dst[shift + j] = '0'; + } + written += shift + j; + } else if (decimal_point > len) { + // right pad with zeros + while (d_exp-- > 0) { + dst[len + written++] = '0'; + } + } else if (len > 1 && d_exp < 0) { + // stick decimal in the middle + memmove(dst + decimal_point + 1, dst + decimal_point, -d_exp); + dst[decimal_point] = '.'; + written++; + } + return written; +} + +static int scientific(char *dst, int len, int decimal_point, flags flags) { + int written = 0; + int exponent = decimal_point - 1; + + // if only 1 digit, don't need a decimal, e.g. 0.0000001 -> 1e-7 or + // 100000000000000000000.0 -> 1e21 + if (len > 1) { + memmove(dst + 1, dst, len); + dst[1] = '.'; + written++; + } + dst[len + written++] = 'e'; + + if (flags & PLUS_IN_EXPONENT && exponent > 0) { + dst[len + written++] = '+'; + } + if (flags & PAD_EXPONENT && exponent > 0 && exponent < 10) { + dst[len + written++] = '0'; + } + + written += i_to_str(exponent, dst+len+written); + + return written; +} + +static int exponential(char *dst, int len, int d_exp) { + int written = 0; + dst[len + written++] = 'e'; + written += i_to_str(d_exp, dst+len+written); + return written; +} + +static int shortest_dtoa(double v, char *dst, int low_exp, int high_exp) { + int d_exp, len, success, i, decimal_point, exp_len; + char *s2 = dst; + assert(dst); + + i = handle_special(v, dst); + if (i >= 0) return i; + + // Grisu3 doesn't handle negative values, so emit the - and negate `v`. + if (v < 0) { + *s2++ = '-'; + v = -v; + } + + success = fast_dtoa(v, FAST_DTOA_SHORTEST, 0, s2, &len, &decimal_point); + // If grisu3 was not able to convert the number to a string, then use old + // sprintf (suboptimal). + if (!success) { + return snprintf(s2, 20, "%.17g", v) + (int)(s2 - dst); + } + + // We now have an integer string of form "151324135" in s2 and a base-10 + // exponent for that number in d_exp. + + d_exp = decimal_point - len; + + // calculate length of printing the exponent. + // cheaper than `(int)log10(d_exp) + 1` (?) and handles negatives + exp_len = + d_exp < -99 ? 4 : + d_exp < -9 ? 3 : + d_exp < 0 ? 2 : + d_exp < 10 ? 1 : + d_exp < 100 ? 2 : + 3; + + if (decimal_point < 0) { + // decimal length = len + 1 (for .) + abs(len+d_exp) (for 0's) + // exponential length = len + 1 (for e) + exp_len + if (-decimal_point <= exp_len) { + len += decimal(s2, len, decimal_point, NO_FLAGS); + } else { + len += exponential(s2, len, d_exp); + } + } else { + // decimal length = len + d_exp (for 0's) + // exponential length = len + 1 (for e) + exp_len + if (d_exp <= 1 + exp_len) { + len += decimal(s2, len, decimal_point, NO_FLAGS); + } else { + len += exponential(s2, len, d_exp); + } + } + + s2[len] = '\0'; + return (int)(s2+len-dst); +} + +static int ecma_dtoa(double v, char *dst) { + int decimal_point, d_exp, len, success, i, exponent; + char *s2 = dst; + assert(dst); + + i = handle_special(v, dst); + if (i >= 0) return i; + + // Grisu3 doesn't handle negative values, so emit the - and negate `v`. + if (v < 0) { + *s2++ = '-'; + v = -v; + } + + success = fast_dtoa(v, FAST_DTOA_SHORTEST, 0, s2, &len, &decimal_point); + // If grisu3 was not able to convert the number to a string, then use old + // sprintf (suboptimal). + if (!success) { + return snprintf(s2, 20, "%.17g", v) + (int)(s2 - dst); + } + + // We now have an integer string of form "151324135" in s2 and the offset of + // the decimal point. + // + // if printed in scientific notation, what would the exponent be? e.g. + // 0.999999 is s2 == "999999", len == 6, decimal_point == 0, so in scientific + // notation it would be 9.99999e-1. + // + // note: we don't actually print 0.999999 in scientific notation given the + // decimal_in_shortest_low rule mentioned above. this is just an example. + exponent = decimal_point - 1; + + if (-6 <= exponent && exponent < 21) { + len += decimal(s2, len, decimal_point, LEADING_ZERO); + } else { + len += scientific(s2, len, decimal_point, PLUS_IN_EXPONENT); + } + + s2[len] = '\0'; + return (int)(s2+len-dst); +} + +// Like David M. Gay's g_fmt, but using grisu3. +// http://www.netlib.org/fp/g_fmt.c +static int grisu3_g_fmt(double v, char *dst) { + int d_exp, len, success, i, decimal_point; + char *s2 = dst; + assert(dst); + + i = handle_special(v, dst); + if (i >= 0) return i; + + // Grisu3 doesn't handle negative values, so emit the - and negate `v`. + if (v < 0) { + *s2++ = '-'; + v = -v; + } + + success = fast_dtoa(v, FAST_DTOA_SHORTEST, 0, s2, &len, &decimal_point); + // If grisu3 was not able to convert the number to a string, then use old + // sprintf (suboptimal). + if (!success) { + return snprintf(s2, 20, "%.17g", v) + (int)(s2 - dst); + } + + if (-3 <= decimal_point && decimal_point < len + 6) { + len += decimal(s2, len, decimal_point, NO_FLAGS); + } else { + len += scientific(s2, len, decimal_point, PLUS_IN_EXPONENT | PAD_EXPONENT); + } + + s2[len] = '\0'; + return (int)(s2+len-dst); +} + +CAMLprim value flow_shortest_string_of_float(value num) { + CAMLparam1(num); + char str[32]; // the max length is probably 18, but better safe than sorry + int len = shortest_dtoa(Double_val(num), str, -3, 3); + assert(len > 0 && len < 25); + CAMLreturn(caml_copy_string(str)); +} + +CAMLprim value flow_ecma_string_of_float(value num) { + CAMLparam1(num); + char str[32]; // the max length is probably 18, but better safe than sorry + int len = ecma_dtoa(Double_val(num), str); + assert(len > 0 && len < 25); + CAMLreturn(caml_copy_string(str)); +} + +CAMLprim value flow_g_fmt(value num) { + CAMLparam1(num); + char str[32]; // the max length is probably 18, but better safe than sorry + int len = grisu3_g_fmt(Double_val(num), str); + assert(len > 0 && len < 25); + CAMLreturn(caml_copy_string(str)); +} diff --git a/src/fast_dtoa.c b/src/fast_dtoa.c new file mode 100644 index 0000000..a5015f7 --- /dev/null +++ b/src/fast_dtoa.c @@ -0,0 +1,425 @@ +/** + * Copyright (c) 2017-present, Facebook, Inc. + * Copyright (c) 2010, the V8 project authors. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + */ + +#include + +#include "cached_powers.h" +#include "fast_dtoa.h" +#include "ieee.h" + +// The minimal and maximal target exponent define the range of w's binary +// exponent, where 'w' is the result of multiplying the input by a cached power +// of ten. +// +// A different range might be chosen on a different platform, to optimize digit +// generation, but a smaller range requires more powers of ten to be cached. +static const int kMinimalTargetExponent = -60; +static const int kMaximalTargetExponent = -32; + + +// Adjusts the last digit of the generated number, and screens out generated +// solutions that may be inaccurate. A solution may be inaccurate if it is +// outside the safe interval, or if we cannot prove that it is closer to the +// input than a neighboring representation of the same length. +// +// Input: * buffer containing the digits of too_high / 10^kappa +// * the buffer's length +// * distance_too_high_w == (too_high - w).f() * unit +// * unsafe_interval == (too_high - too_low).f() * unit +// * rest = (too_high - buffer * 10^kappa).f() * unit +// * ten_kappa = 10^kappa * unit +// * unit = the common multiplier +// Output: returns true if the buffer is guaranteed to contain the closest +// representable number to the input. +// Modifies the generated digits in the buffer to approach (round towards) w. +static bool round_weed(char* buffer, + int length, + uint64_t distance_too_high_w, + uint64_t unsafe_interval, + uint64_t rest, + uint64_t ten_kappa, + uint64_t unit) { + uint64_t small_distance = distance_too_high_w - unit; + uint64_t big_distance = distance_too_high_w + unit; + // Let w_low = too_high - big_distance, and + // w_high = too_high - small_distance. + // Note: w_low < w < w_high + // + // The real w (* unit) must lie somewhere inside the interval + // ]w_low; w_high[ (often written as "(w_low; w_high)") + + // Basically the buffer currently contains a number in the unsafe interval + // ]too_low; too_high[ with too_low < w < too_high + // + // too_high - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + // ^v 1 unit ^ ^ ^ ^ + // boundary_high --------------------- . . . . + // ^v 1 unit . . . . + // - - - - - - - - - - - - - - - - - - - + - - + - - - - - - . . + // . . ^ . . + // . big_distance . . . + // . . . . rest + // small_distance . . . . + // v . . . . + // w_high - - - - - - - - - - - - - - - - - - . . . . + // ^v 1 unit . . . . + // w ---------------------------------------- . . . . + // ^v 1 unit v . . . + // w_low - - - - - - - - - - - - - - - - - - - - - . . . + // . . v + // buffer --------------------------------------------------+-------+-------- + // . . + // safe_interval . + // v . + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - . + // ^v 1 unit . + // boundary_low ------------------------- unsafe_interval + // ^v 1 unit v + // too_low - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + // + // + // Note that the value of buffer could lie anywhere inside the range too_low + // to too_high. + // + // boundary_low, boundary_high and w are approximations of the real boundaries + // and v (the input number). They are guaranteed to be precise up to one unit. + // In fact the error is guaranteed to be strictly less than one unit. + // + // Anything that lies outside the unsafe interval is guaranteed not to round + // to v when read again. + // Anything that lies inside the safe interval is guaranteed to round to v + // when read again. + // If the number inside the buffer lies inside the unsafe interval but not + // inside the safe interval then we simply do not know and bail out (returning + // false). + // + // Similarly we have to take into account the imprecision of 'w' when finding + // the closest representation of 'w'. If we have two potential + // representations, and one is closer to both w_low and w_high, then we know + // it is closer to the actual value v. + // + // By generating the digits of too_high we got the largest (closest to + // too_high) buffer that is still in the unsafe interval. In the case where + // w_high < buffer < too_high we try to decrement the buffer. + // This way the buffer approaches (rounds towards) w. + // There are 3 conditions that stop the decrementation process: + // 1) the buffer is already below w_high + // 2) decrementing the buffer would make it leave the unsafe interval + // 3) decrementing the buffer would yield a number below w_high and farther + // away than the current number. In other words: + // (buffer{-1} < w_high) && w_high - buffer{-1} > buffer - w_high + // Instead of using the buffer directly we use its distance to too_high. + // Conceptually rest ~= too_high - buffer + // We need to do the following tests in this order to avoid over- and + // underflows. + ASSERT(rest <= unsafe_interval); + while (rest < small_distance && // Negated condition 1 + unsafe_interval - rest >= ten_kappa && // Negated condition 2 + (rest + ten_kappa < small_distance || // buffer{-1} > w_high + small_distance - rest >= rest + ten_kappa - small_distance)) { + buffer[length - 1]--; + rest += ten_kappa; + } + + // We have approached w+ as much as possible. We now test if approaching w- + // would require changing the buffer. If yes, then we have two possible + // representations close to w, but we cannot decide which one is closer. + if (rest < big_distance && + unsafe_interval - rest >= ten_kappa && + (rest + ten_kappa < big_distance || + big_distance - rest > rest + ten_kappa - big_distance)) { + return false; + } + + // Weeding test. + // The safe interval is [too_low + 2 ulp; too_high - 2 ulp] + // Since too_low = too_high - unsafe_interval this is equivalent to + // [too_high - unsafe_interval + 4 ulp; too_high - 2 ulp] + // Conceptually we have: rest ~= too_high - buffer + return (2 * unit <= rest) && (rest <= unsafe_interval - 4 * unit); +} + + +// Returns the biggest power of ten that is less than or equal to the given +// number. We furthermore receive the maximum number of bits 'number' has. +// +// Returns power == 10^(exponent_plus_one-1) such that +// power <= number < power * 10. +// If number_bits == 0 then 0^(0-1) is returned. +// The number of bits must be <= 32. +// Precondition: number < (1 << (number_bits + 1)). + +// Inspired by the method for finding an integer log base 10 from here: +// http://graphics.stanford.edu/~seander/bithacks.html#IntegerLog10 +static unsigned int const kSmallPowersOfTen[] = + {0, 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, + 1000000000}; + +static void biggest_power_ten(uint32_t number, + int number_bits, + uint32_t* power, + int* exponent_plus_one) { + ASSERT(number < (1u << (number_bits + 1))); + // 1233/4096 is approximately 1/lg(10). + int exponent_plus_one_guess = ((number_bits + 1) * 1233 >> 12); + // We increment to skip over the first entry in the kPowersOf10 table. + // Note: kPowersOf10[i] == 10^(i-1). + exponent_plus_one_guess++; + // We don't have any guarantees that 2^number_bits <= number. + if (number < kSmallPowersOfTen[exponent_plus_one_guess]) { + exponent_plus_one_guess--; + } + *power = kSmallPowersOfTen[exponent_plus_one_guess]; + *exponent_plus_one = exponent_plus_one_guess; +} + +// Generates the digits of input number w. +// w is a floating-point number (DiyFp), consisting of a significand and an +// exponent. Its exponent is bounded by kMinimalTargetExponent and +// kMaximalTargetExponent. +// Hence -60 <= w.e <= -32. +// +// Returns false if it fails, in which case the generated digits in the buffer +// should not be used. +// Preconditions: +// * low, w and high are correct up to 1 ulp (unit in the last place). That +// is, their error must be less than a unit of their last digits. +// * low.e == w.e == high.e +// * low < w < high, and taking into account their error: low~ <= high~ +// * kMinimalTargetExponent <= w.e <= kMaximalTargetExponent +// Postconditions: returns false if procedure fails. +// otherwise: +// * buffer is not null-terminated, but len contains the number of digits. +// * buffer contains the shortest possible decimal digit-sequence +// such that LOW < buffer * 10^kappa < HIGH, where LOW and HIGH are the +// correct values of low and high (without their error). +// * if more than one decimal representation gives the minimal number of +// decimal digits then the one closest to W (where W is the correct value +// of w) is chosen. +// Remark: this procedure takes into account the imprecision of its input +// numbers. If the precision is not enough to guarantee all the postconditions +// then false is returned. This usually happens rarely (~0.5%). +// +// Say, for the sake of example, that +// w.e == -48, and w.f == 0x1234567890abcdef +// w's value can be computed by w.f * 2^w.e +// We can obtain w's integral digits by simply shifting w.f by -w.e. +// -> w's integral part is 0x1234 +// w's fractional part is therefore 0x567890abcdef. +// Printing w's integral part is easy (simply print 0x1234 in decimal). +// In order to print its fraction we repeatedly multiply the fraction by 10 and +// get each digit. Example the first digit after the point would be computed by +// (0x567890abcdef * 10) >> 48. -> 3 +// The whole thing becomes slightly more complicated because we want to stop +// once we have enough digits. That is, once the digits inside the buffer +// represent 'w' we can stop. Everything inside the interval low - high +// represents w. However we have to pay attention to low, high and w's +// imprecision. +static bool digit_gen(diy_fp low, + diy_fp w, + diy_fp high, + char* buffer, + int* length, + int* kappa) { + ASSERT(low.e == w.e && w.e == high.e); + ASSERT(low.f + 1 <= high.f - 1); + ASSERT(kMinimalTargetExponent <= w.e && w.e <= kMaximalTargetExponent); + // low, w and high are imprecise, but by less than one ulp (unit in the last + // place). + // If we remove (resp. add) 1 ulp from low (resp. high) we are certain that + // the new numbers are outside of the interval we want the final + // representation to lie in. + // Inversely adding (resp. removing) 1 ulp from low (resp. high) would yield + // numbers that are certain to lie in the interval. We will use this fact + // later on. + // We will now start by generating the digits within the uncertain + // interval. Later we will weed out representations that lie outside the safe + // interval and thus _might_ lie outside the correct interval. + uint64_t unit = 1; + diy_fp too_low = { low.f - unit, low.e }; + diy_fp too_high = { high.f + unit, high.e }; + // too_low and too_high are guaranteed to lie outside the interval we want the + // generated number in. + diy_fp unsafe_interval = diy_fp_minus(too_high, too_low); + // We now cut the input number into two parts: the integral digits and the + // fractionals. We will not write any decimal separator though, but adapt + // kappa instead. + // Reminder: we are currently computing the digits (stored inside the buffer) + // such that: too_low < buffer * 10^kappa < too_high + // We use too_high for the digit_generation and stop as soon as possible. + // If we stop early we effectively round down. + diy_fp one = { (uint64_t)(1) << -w.e, w.e }; + // Division by one is a shift. + uint32_t integrals = (uint32_t)(too_high.f >> -one.e); + // Modulo by one is an and. + uint64_t fractionals = too_high.f & (one.f - 1); + uint32_t divisor; + int divisor_exponent_plus_one; + biggest_power_ten(integrals, kDiyFpSignificandSize - (-one.e), + &divisor, &divisor_exponent_plus_one); + *kappa = divisor_exponent_plus_one; + *length = 0; + // Loop invariant: buffer = too_high / 10^kappa (integer division) + // The invariant holds for the first iteration: kappa has been initialized + // with the divisor exponent + 1. And the divisor is the biggest power of ten + // that is smaller than integrals. + while (*kappa > 0) { + int digit = integrals / divisor; + ASSERT(digit <= 9); + buffer[*length] = (char)('0' + digit); + (*length)++; + integrals %= divisor; + (*kappa)--; + // Note that kappa now equals the exponent of the divisor and that the + // invariant thus holds again. + uint64_t rest = ((uint64_t)(integrals) << -one.e) + fractionals; + // Invariant: too_high = buffer * 10^kappa + DiyFp(rest, one.e) + // Reminder: unsafe_interval.e == one.e + if (rest < unsafe_interval.f) { + // Rounding down (by not emitting the remaining digits) yields a number + // that lies within the unsafe interval. + return round_weed(buffer, *length, diy_fp_minus(too_high, w).f, + unsafe_interval.f, rest, + (uint64_t)(divisor) << -one.e, unit); + } + divisor /= 10; + } + + // The integrals have been generated. We are at the point of the decimal + // separator. In the following loop we simply multiply the remaining digits by + // 10 and divide by one. We just need to pay attention to multiply associated + // data (like the interval or 'unit'), too. + // Note that the multiplication by 10 does not overflow, because w.e >= -60 + // and thus one.e >= -60. + ASSERT(one.e >= -60); + ASSERT(fractionals < one.f); + ASSERT(UINT64_2PART_C(0xFFFFFFFF, FFFFFFFF) / 10 >= one.f); + for (;;) { + fractionals *= 10; + unit *= 10; + unsafe_interval.f = unsafe_interval.f * 10; + // Integer division by one. + int digit = (int)(fractionals >> -one.e); + ASSERT(digit <= 9); + buffer[*length] = (char)('0' + digit); + (*length)++; + fractionals &= one.f - 1; // Modulo by one. + (*kappa)--; + if (fractionals < unsafe_interval.f) { + return round_weed(buffer, *length, diy_fp_minus(too_high, w).f * unit, + unsafe_interval.f, fractionals, one.f, unit); + } + } +} + +// Provides a decimal representation of v. +// Returns true if it succeeds, otherwise the result cannot be trusted. +// There will be *length digits inside the buffer (not null-terminated). +// If the function returns true then +// v == (double) (buffer * 10^decimal_exponent). +// The digits in the buffer are the shortest representation possible: no +// 0.09999999999999999 instead of 0.1. The shorter representation will even be +// chosen even if the longer one would be closer to v. +// The last digit will be closest to the actual v. That is, even if several +// digits might correctly yield 'v' when read again, the closest will be +// computed. +static bool grisu3(double v, + FastDtoaMode mode, + char* buffer, + int* length, + int* decimal_exponent) { + diy_fp w = double_as_normalized_diy_fp(v); + + // boundary_minus and boundary_plus are the boundaries between v and its + // closest floating-point neighbors. Any number strictly between + // boundary_minus and boundary_plus will round to v when convert to a double. + // Grisu3 will never output representations that lie exactly on a boundary. + diy_fp boundary_minus, boundary_plus; + if (mode == FAST_DTOA_SHORTEST) { + double_normalized_boundaries(v, &boundary_minus, &boundary_plus); + } else { + ASSERT(mode == FAST_DTOA_SHORTEST_SINGLE); + float single_v = (float)(v); + // Single(single_v).NormalizedBoundaries(&boundary_minus, &boundary_plus); + } + ASSERT(boundary_plus.e == w.e); + diy_fp ten_mk; // Cached power of ten: 10^-k + int mk; // -k + int ten_mk_minimal_binary_exponent = + kMinimalTargetExponent - (w.e + kDiyFpSignificandSize); + int ten_mk_maximal_binary_exponent = + kMaximalTargetExponent - (w.e + kDiyFpSignificandSize); + cached_power_for_binary_exponent_range( + ten_mk_minimal_binary_exponent, + ten_mk_maximal_binary_exponent, + &ten_mk, &mk); + ASSERT((kMinimalTargetExponent <= w.e + ten_mk.e + kDiyFpSignificandSize) && + (kMaximalTargetExponent >= w.e + ten_mk.e + kDiyFpSignificandSize)); + // Note that ten_mk is only an approximation of 10^-k. A diy_fp only contains + // a 64 bit significand and ten_mk is thus only precise up to 64 bits. + + // The diy_fp_multiply procedure rounds its result, and ten_mk is approximated + // too. The variable scaled_w (as well as scaled_boundary_minus/plus) are now + // off by a small amount. + // In fact: scaled_w - w*10^k < 1ulp (unit in the last place) of scaled_w. + // In other words: let f = scaled_w.f and e = scaled_w.e, then + // (f-1) * 2^e < w*10^k < (f+1) * 2^e + diy_fp scaled_w = diy_fp_multiply(w, ten_mk); + ASSERT(scaled_w.e == boundary_plus.e + ten_mk.e + kDiyFpSignificandSize); + // In theory it would be possible to avoid some recomputations by computing + // the difference between w and boundary_minus/plus (a power of 2) and to + // compute scaled_boundary_minus/plus by subtracting/adding from + // scaled_w. However the code becomes much less readable and the speed + // enhancements are not terriffic. + diy_fp scaled_boundary_minus = diy_fp_multiply(boundary_minus, ten_mk); + diy_fp scaled_boundary_plus = diy_fp_multiply(boundary_plus, ten_mk); + + // digit_gen will generate the digits of scaled_w. Therefore we have + // v == (double) (scaled_w * 10^-mk). + // Set decimal_exponent == -mk and pass it to digit_gen. If scaled_w is not an + // integer than it will be updated. For instance if scaled_w == 1.23 then + // the buffer will be filled with "123" und the decimal_exponent will be + // decreased by 2. + int kappa; + bool result = digit_gen(scaled_boundary_minus, scaled_w, scaled_boundary_plus, + buffer, length, &kappa); + *decimal_exponent = -mk + kappa; + return result; +} + +bool fast_dtoa(double v, + FastDtoaMode mode, + int requested_digits, + char* buffer, + int* length, + int* decimal_point) { + ASSERT(v > 0); + ASSERT(!double_is_special(v)); + + bool result = false; + int decimal_exponent = 0; + switch (mode) { + case FAST_DTOA_SHORTEST: + case FAST_DTOA_SHORTEST_SINGLE: + result = grisu3(v, mode, buffer, length, &decimal_exponent); + break; + case FAST_DTOA_PRECISION: + // result = Grisu3Counted(v, requested_digits, + // buffer, length, &decimal_exponent); + break; + // default: + // UNREACHABLE(); + } + if (result) { + *decimal_point = *length + decimal_exponent; + buffer[*length] = '\0'; + } + return result; +} diff --git a/src/fast_dtoa.h b/src/fast_dtoa.h new file mode 100644 index 0000000..926d046 --- /dev/null +++ b/src/fast_dtoa.h @@ -0,0 +1,61 @@ +/** + * Copyright (c) 2017-present, Facebook, Inc. + * Copyright (c) 2010, the V8 project authors. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + */ + +#ifndef OCAML_DTOA_FAST_DTOA_H_ +#define OCAML_DTOA_FAST_DTOA_H_ + +#include + +typedef enum { + // Computes the shortest representation of the given input. The returned + // result will be the most accurate number of this length. Longer + // representations might be more accurate. + FAST_DTOA_SHORTEST, + // Same as FAST_DTOA_SHORTEST but for single-precision floats. + FAST_DTOA_SHORTEST_SINGLE, + // Computes a representation where the precision (number of digits) is + // given as input. The precision is independent of the decimal point. + FAST_DTOA_PRECISION +} FastDtoaMode; + +// Provides a decimal representation of v. +// The result should be interpreted as buffer * 10^(point - length). +// +// Precondition: +// * v must be a strictly positive finite double. +// +// Returns true if it succeeds, otherwise the result can not be trusted. +// There will be *length digits inside the buffer followed by a null terminator. +// If the function returns true and mode equals +// - FAST_DTOA_SHORTEST, then +// the parameter requested_digits is ignored. +// The result satisfies +// v == (double) (buffer * 10^(point - length)). +// The digits in the buffer are the shortest representation possible. E.g. +// if 0.099999999999 and 0.1 represent the same double then "1" is returned +// with point = 0. +// The last digit will be closest to the actual v. That is, even if several +// digits might correctly yield 'v' when read again, the buffer will contain +// the one closest to v. +// - FAST_DTOA_PRECISION, then +// the buffer contains requested_digits digits. +// the difference v - (buffer * 10^(point-length)) is closest to zero for +// all possible representations of requested_digits digits. +// If there are two values that are equally close, then fast_dtoa returns +// false. +// For both modes the buffer must be large enough to hold the result. +bool fast_dtoa(double d, + FastDtoaMode mode, + int requested_digits, + char* buffer, + int* length, + int* decimal_point); + +#endif // OCAML_DTOA_FAST_DTOA_H_ diff --git a/src/ieee.c b/src/ieee.c new file mode 100644 index 0000000..228e53e --- /dev/null +++ b/src/ieee.c @@ -0,0 +1,143 @@ +/** + * Copyright (c) 2017-present, Facebook, Inc. + * Copyright (c) 2010, the V8 project authors. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + */ + +#include +#include + +#include "diy_fp.h" +#include "ieee.h" +#include "utils.h" + +const uint64_t kSignMask = UINT64_2PART_C(0x80000000, 00000000); +const uint64_t kExponentMask = UINT64_2PART_C(0x7FF00000, 00000000); +const uint64_t kSignificandMask = UINT64_2PART_C(0x000FFFFF, FFFFFFFF); +const uint64_t kHiddenBit = UINT64_2PART_C(0x00100000, 00000000); +const int kDoublePhysicalSignificandSize = 52; // Excludes the hidden bit. +const int kDoubleSignificandSize = 53; + +const int kDoubleExponentBias = 0x433; // 0x3FF + kDoublePhysicalSignificandSize +const int kDoubleDenormalExponent = -0x432; // -kDoubleExponentBias + 1; + +static uint64_t as_uint64(double d) { + union { double d; uint64_t u64; } u; + u.d = d; + return u.u64; +} + +static bool d64_is_denormal(uint64_t d64) { + return (d64 & kExponentMask) == 0; +} + +// Returns true if the double is a denormal. +bool double_is_denormal(double d) { + return d64_is_denormal(as_uint64(d)); +} + +static bool d64_is_special(uint64_t d64) { + return (d64 & kExponentMask) == kExponentMask; +} + +bool double_is_special(double d) { + return d64_is_special(as_uint64(d)); +} + +static uint64_t d64_significand(uint64_t d64) { + uint64_t significand = d64 & kSignificandMask; + if (!d64_is_denormal(d64)) { + return significand + kHiddenBit; + } else { + return significand; + } +} + +static int d64_exponent(uint64_t d64) { + if (d64_is_denormal(d64)) return kDoubleDenormalExponent; + + int biased_e = (int)((d64 & kExponentMask) >> kDoublePhysicalSignificandSize); + return biased_e - kDoubleExponentBias; +} + +static int d64_sign(uint64_t d64) { + return (d64 & kSignMask) == 0? 1: -1; +} + +int double_sign(double d) { + return d64_sign(as_uint64(d)); +} + +// The value encoded by this Double must be greater or equal to +0.0. +// It must not be special (infinity, or NaN). +diy_fp double_as_diy_fp(double d) { + diy_fp r; + uint64_t d64 = as_uint64(d); + ASSERT(d64_sign(d64) > 0); + ASSERT(!d64_is_special(d64)); + r.f = d64_significand(d64); + r.e = d64_exponent(d64); + return r; +} + +diy_fp double_as_normalized_diy_fp(double d) { + diy_fp r; + ASSERT(d > 0.0); + uint64_t d64 = as_uint64(d); + uint64_t f = d64_significand(d64); + int e = d64_exponent(d64); + + // The current double could be a denormal. + while ((f & kHiddenBit) == 0) { + f <<= 1; + e--; + } + // Do the final shifts in one go. + f <<= kDiyFpSignificandSize - kDoubleSignificandSize; + e -= kDiyFpSignificandSize - kDoubleSignificandSize; + r.f = f; + r.e = e; + return r; +} + +static bool double_lower_boundary_is_closer(double d) { + // The boundary is closer if the significand is of the form f == 2^p-1 then + // the lower boundary is closer. + // Think of v = 1000e10 and v- = 9999e9. + // Then the boundary (== (v - v-)/2) is not just at a distance of 1e9 but + // at a distance of 1e8. + // The only exception is for the smallest normal: the largest denormal is + // at the same distance as its successor. + // Note: denormals have the same exponent as the smallest normals. + uint64_t d64 = as_uint64(d); + bool physical_significand_is_zero = ((d64 & kSignificandMask) == 0); + return physical_significand_is_zero && + (d64_exponent(d64) != kDoubleDenormalExponent); +} + +void double_normalized_boundaries( + double d, + diy_fp* out_m_minus, + diy_fp* out_m_plus +) { + ASSERT(d > 0.0); + diy_fp v = double_as_diy_fp(d); + diy_fp tmp = { (v.f << 1) + 1, v.e - 1 }; + diy_fp m_plus = diy_fp_normalize(tmp); + diy_fp m_minus; + if (double_lower_boundary_is_closer(d)) { + m_minus.f = (v.f << 2) - 1; + m_minus.e = v.e - 2; + } else { + m_minus.f = (v.f << 1) - 1; + m_minus.e = v.e - 1; + } + m_minus.f = m_minus.f << (m_minus.e - m_plus.e); + m_minus.e = m_plus.e; + *out_m_plus = m_plus; + *out_m_minus = m_minus; +} diff --git a/src/ieee.h b/src/ieee.h new file mode 100644 index 0000000..2bb2ff8 --- /dev/null +++ b/src/ieee.h @@ -0,0 +1,43 @@ +/** + * Copyright (c) 2017-present, Facebook, Inc. + * Copyright (c) 2010, the V8 project authors. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + */ + +#ifndef OCAML_DTOA_IEEE_H_ +#define OCAML_DTOA_IEEE_H_ + +#include +#include "diy_fp.h" +#include "utils.h" + +extern const uint64_t kSignMask; +extern const uint64_t kExponentMask; +extern const uint64_t kSignificandMask; +extern const uint64_t kHiddenBit; +extern const int kPhysicalSignificandSize; +extern const int kDoubleSignificandSize; + +// Returns true if the double is a denormal. +bool double_is_denormal(double d); + +// We consider denormals not to be special. +// Hence only Infinity and NaN are special. +bool double_is_special(double d); + +// The value encoded by this double must be strictly greater than 0. +diy_fp double_as_normalized_diy_fp(double d); + +// Computes the two boundaries of this. +// The bigger boundary (m_plus) is normalized. The lower boundary has the same +// exponent as m_plus. +// Precondition: the value encoded by this Single must be greater than 0. +void double_normalized_boundaries( + double d, diy_fp* out_m_minus, diy_fp* out_m_plus +); + +#endif // OCAML_DTOA_IEEE_H_ diff --git a/src/utils.h b/src/utils.h new file mode 100644 index 0000000..00e99b7 --- /dev/null +++ b/src/utils.h @@ -0,0 +1,37 @@ +/** + * Copyright (c) 2017-present, Facebook, Inc. + * Copyright (c) 2010, the V8 project authors. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + */ + +#ifndef OCAML_DTOA_UTILS_H_ +#define OCAML_DTOA_UTILS_H_ + +#include + +#include +#ifndef ASSERT +#define ASSERT(condition) \ + assert(condition); +#endif + +// The following macro works on both 32 and 64-bit platforms. +// Usage: instead of writing 0x1234567890123456 +// write UINT64_2PART_C(0x12345678,90123456); +#define UINT64_2PART_C(a, b) (((uint64_t)(a) << 32) + 0x##b##u) + +// The expression ARRAY_SIZE(a) is a compile-time constant of type +// size_t which represents the number of elements of the given +// array. You should only use ARRAY_SIZE on statically allocated +// arrays. +#ifndef ARRAY_SIZE +#define ARRAY_SIZE(a) \ + ((sizeof(a) / sizeof(*(a))) / \ + (size_t)(!(sizeof(a) % sizeof(*(a))))) +#endif + +#endif // OCAML_DTOA_UTILS_H_ diff --git a/test/test.ml b/test/test.ml new file mode 100644 index 0000000..ce93af2 --- /dev/null +++ b/test/test.ml @@ -0,0 +1,21 @@ +(*--------------------------------------------------------------------------- + Copyright (c) 2017 Marshall Roch. All rights reserved. + Distributed under the ISC license, see terms at the end of the file. + %%NAME%% %%VERSION%% + ---------------------------------------------------------------------------*) + +(*--------------------------------------------------------------------------- + Copyright (c) 2017 Marshall Roch + + Permission to use, copy, modify, and/or distribute this software for any + purpose with or without fee is hereby granted, provided that the above + copyright notice and this permission notice appear in all copies. + + THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES + WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF + MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR + ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES + WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN + ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF + OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + ---------------------------------------------------------------------------*) diff --git a/tests/ecma_test.ml b/tests/ecma_test.ml new file mode 100644 index 0000000..4696c96 --- /dev/null +++ b/tests/ecma_test.ml @@ -0,0 +1,63 @@ +(** + * Copyright (c) 2017-present, Facebook, Inc. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + *) + +open OUnit2 +open Dtoa + +let eq ctxt = assert_equal ~ctxt ~printer:(fun x -> x) + +let tests = "ecma_string_of_float" >::: [ + "sanity" >:: begin fun ctxt -> + let eq = eq ctxt in + eq "1" (ecma_string_of_float 1.); + eq "0" (ecma_string_of_float 0.); + eq "-1" (ecma_string_of_float (-.1.)); + end; + + "big_and_small" >:: begin fun ctxt -> + let eq = eq ctxt in + + eq "0.30000000000000004" (ecma_string_of_float (0.1 +. 0.2)); + eq "2592000000" (ecma_string_of_float 2.592e+09); + eq "2592000000" (ecma_string_of_float 2592000000.); + eq "0.9999999999999999" (ecma_string_of_float 0.99999999999999994); + eq "1" (ecma_string_of_float 0.99999999999999995); + eq "100" (ecma_string_of_float 100.); + eq "1000" (ecma_string_of_float 1000.); + + eq "0.000001" (ecma_string_of_float 0.000001); + eq "0.999999" (ecma_string_of_float 0.999999); + eq "1e-7" (ecma_string_of_float 0.0000001); + eq "0.9999999" (ecma_string_of_float 0.9999999); + eq "0.000001" (ecma_string_of_float 1.0e-6); + eq "1e-7" (ecma_string_of_float 0.1e-6); + eq "1e-7" (ecma_string_of_float 1.0e-7); + eq "111111111111111110000" (ecma_string_of_float 111111111111111111111.0); + eq "100000000000000000000" (ecma_string_of_float 100000000000000000000.0); + eq "1.1111111111111111e+21" (ecma_string_of_float 1111111111111111111111.0); + eq "1e+21" (ecma_string_of_float 1000000000000000000000.0); + eq "4294967272" (ecma_string_of_float 4294967272.0); + eq "4.185580496821357e+298" (ecma_string_of_float 4.1855804968213567e298); + end; + + "edge_cases" >:: begin fun ctxt -> + let eq = eq ctxt in + + eq "5e-324" (ecma_string_of_float 5e-324); + eq "1.7976931348623157e+308" (ecma_string_of_float 1.7976931348623157e308); + eq "NaN" (ecma_string_of_float Pervasives.nan); + eq "Infinity" (ecma_string_of_float Pervasives.infinity); + eq "-Infinity" (ecma_string_of_float Pervasives.neg_infinity); + end; + + "round_trip" >:: begin fun ctxt -> + let x = ecma_string_of_float epsilon_float |> float_of_string in + assert_equal ~ctxt epsilon_float x; + end; +] diff --git a/tests/g_fmt_test.ml b/tests/g_fmt_test.ml new file mode 100644 index 0000000..34c76a5 --- /dev/null +++ b/tests/g_fmt_test.ml @@ -0,0 +1,63 @@ +(** + * Copyright (c) 2017-present, Facebook, Inc. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + *) + +open OUnit2 +open Dtoa + +let eq ctxt = assert_equal ~ctxt ~printer:(fun x -> x) + +let tests = "g_fmt" >::: [ + "sanity" >:: begin fun ctxt -> + let eq = eq ctxt in + eq "1" (g_fmt 1.); + eq "0" (g_fmt 0.); + eq "-1" (g_fmt (-.1.)); + end; + + "big_and_small" >:: begin fun ctxt -> + let eq = eq ctxt in + + eq ".30000000000000004" (g_fmt (0.1 +. 0.2)); + eq "2.592e+09" (g_fmt 2.592e+09); + eq "2.592e+09" (g_fmt 2592000000.); + eq ".9999999999999999" (g_fmt 0.99999999999999994); + eq "1" (g_fmt 0.99999999999999995); + eq "100" (g_fmt 100.); + eq "1000" (g_fmt 1000.); + + eq "1e-6" (g_fmt 0.000001); + eq ".999999" (g_fmt 0.999999); + eq "1e-7" (g_fmt 0.0000001); + eq ".9999999" (g_fmt 0.9999999); + eq "1e-6" (g_fmt 1.0e-6); + eq "1e-7" (g_fmt 0.1e-6); + eq "1e-7" (g_fmt 1.0e-7); + eq "111111111111111110000" (g_fmt 111111111111111111111.0); + eq "1e+20" (g_fmt 100000000000000000000.0); + eq "1111111111111111100000" (g_fmt 1111111111111111111111.0); + eq "1e+21" (g_fmt 1000000000000000000000.0); + eq "4294967272" (g_fmt 4294967272.0); + eq "4.185580496821357e+298" (g_fmt 4.1855804968213567e298); + end; + + "edge_cases" >:: begin fun ctxt -> + let eq = eq ctxt in + + eq "5e-324" (g_fmt 5e-324); + eq "1.7976931348623157e+308" (g_fmt 1.7976931348623157e308); + eq "NaN" (g_fmt Pervasives.nan); + eq "Infinity" (g_fmt Pervasives.infinity); + eq "-Infinity" (g_fmt Pervasives.neg_infinity); + end; + + "round_trip" >:: begin fun ctxt -> + let x = g_fmt epsilon_float |> float_of_string in + assert_equal ~ctxt epsilon_float x; + end; +] diff --git a/tests/shortest_test.ml b/tests/shortest_test.ml new file mode 100644 index 0000000..e9868ae --- /dev/null +++ b/tests/shortest_test.ml @@ -0,0 +1,67 @@ +(** + * Copyright (c) 2017-present, Facebook, Inc. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + *) + +open OUnit2 +open Dtoa + +let eq ctxt = assert_equal ~ctxt ~printer:(fun x -> x) + +let tests = "shortest_string_of_float" >::: [ + "sanity" >:: begin fun ctxt -> + let eq = eq ctxt in + eq "1" (shortest_string_of_float 1.); + eq "0" (shortest_string_of_float 0.); + eq "-1" (shortest_string_of_float (-.1.)); + end; + + "big_and_small" >:: begin fun ctxt -> + let eq = eq ctxt in + + eq ".30000000000000004" (shortest_string_of_float (0.1 +. 0.2)); + eq "2592e6" (shortest_string_of_float 2.592e+09); + eq "2592e6" (shortest_string_of_float 2592000000.); + eq ".025" (shortest_string_of_float 0.025); + eq ".9999999999999999" (shortest_string_of_float 0.99999999999999994); + eq "1" (shortest_string_of_float 0.99999999999999995); + eq "100" (shortest_string_of_float 100.); + eq "1e3" (shortest_string_of_float 1000.); + eq ".11" (shortest_string_of_float 1.1e-1); + eq ".1" (shortest_string_of_float 1.0e-1); + eq ".01" (shortest_string_of_float 1.0e-2); + eq ".001" (shortest_string_of_float 1.0e-3); + eq "1e-6" (shortest_string_of_float 0.000001); + eq ".999999" (shortest_string_of_float 0.999999); + eq "1e-7" (shortest_string_of_float 0.0000001); + eq ".9999999" (shortest_string_of_float 0.9999999); + eq "1e-6" (shortest_string_of_float 1.0e-6); + eq "1e-7" (shortest_string_of_float 0.1e-6); + eq "1e-7" (shortest_string_of_float 1.0e-7); + eq "11111111111111111e4" (shortest_string_of_float 111111111111111111111.0); + eq "1e20" (shortest_string_of_float 100000000000000000000.0); + eq "11111111111111111e5" (shortest_string_of_float 1111111111111111111111.0); + eq "1e21" (shortest_string_of_float 1000000000000000000000.0); + eq "4294967272" (shortest_string_of_float 4294967272.0); + eq "4185580496821357e283" (shortest_string_of_float 4.1855804968213567e298); + end; + + "edge_cases" >:: begin fun ctxt -> + let eq = eq ctxt in + + eq "5e-324" (shortest_string_of_float 5e-324); + eq "17976931348623157e292" (shortest_string_of_float 1.7976931348623157e308); + eq "NaN" (shortest_string_of_float Pervasives.nan); + eq "Infinity" (shortest_string_of_float Pervasives.infinity); + eq "-Infinity" (shortest_string_of_float Pervasives.neg_infinity); + end; + + "round_trip" >:: begin fun ctxt -> + let x = shortest_string_of_float epsilon_float |> float_of_string in + assert_equal ~ctxt ~printer:string_of_float epsilon_float x; + end; +] diff --git a/tests/test.ml b/tests/test.ml new file mode 100644 index 0000000..d32cae5 --- /dev/null +++ b/tests/test.ml @@ -0,0 +1,18 @@ +(** + * Copyright (c) 2017-present, Facebook, Inc. + * All rights reserved. + * + * This source code is licensed under the BSD-style license found in the + * LICENSE file in the root directory of this source tree. An additional grant + * of patent rights can be found in the PATENTS file in the same directory. + *) + +open OUnit2 + +let tests = "dtoa" >::: [ + Ecma_test.tests; + Shortest_test.tests; + G_fmt_test.tests; +] + +let () = run_test_tt_main tests