Read and write numbers exactly without strtod or printf in the common cases

5 files changed, 172 insertions, 2 deletions

diff --git a/Makefile b/Makefile
index 8848d7b..2ad79d1 100644
--- a/Makefile
+++ b/Makefile
@@ -16,7 +16,8 @@ BUILD ?= build
 #   cli/    command-line hosts.
 # libpure.a is the app model with the tree store: what a host embeds (Vega's
 # surfaces link it). prototypes/ keeps the earlier document runtime (000-077).
-# Unit tests sit by their part (core/test-module.c, store/test-tree-store.c);
+# Unit tests sit by their part (core/test-module.c, store/test-tree-store.c,
+# serve/test-message-json.c);
 # tests/check.py checks the hosts end to end.
 INCLUDES = -Icore
 MODEL = core/module.c core/module-list.c core/module-save.c core/module-text.c core/module-store.c core/flow.c
@@ -59,6 +60,11 @@ $(BUILD)/test-module: core/test-module.c $(MODEL) $(STORE) $(HEADERS)
 	@mkdir -p $(BUILD)
 	$(CC) $(CPPFLAGS) $(CFLAGS) $(INCLUDES) -DPURE_STORE_TEST core/test-module.c $(MODEL) $(STORE) $(LDFLAGS) -pthread -lm $(LDLIBS) -o $@
 
+$(BUILD)/test-message-json: serve/test-message-json.c serve/message-json.c $(BUILD)/libpure.a $(HEADERS)
+	@mkdir -p $(BUILD)
+	$(CC) $(CPPFLAGS) $(CFLAGS) $(INCLUDES) -Iserve serve/test-message-json.c serve/message-json.c $(BUILD)/libpure.a \
+		$(LDFLAGS) -lm $(LDLIBS) -o $@
+
 $(BUILD)/test-tree-store: store/test-tree-store.c $(STORE) $(HEADERS)
 	@mkdir -p $(BUILD)
 	$(CC) $(CPPFLAGS) $(CFLAGS) $(INCLUDES) -DPURE_STORE_TEST store/test-tree-store.c $(STORE) $(LDFLAGS) $(LDLIBS) -o $@
@@ -66,9 +72,10 @@ $(BUILD)/test-tree-store: store/test-tree-store.c $(STORE) $(HEADERS)
 # The HTTP host's checks build their apps with aup (the sibling checkout) and
 # compare documents with the earlier document server.
 AUP ?= ../aup
-check: all $(BUILD)/test-module $(BUILD)/test-tree-store
+check: all $(BUILD)/test-module $(BUILD)/test-tree-store $(BUILD)/test-message-json
 	$(BUILD)/test-tree-store
 	$(BUILD)/test-module
+	$(BUILD)/test-message-json
 	$(MAKE) -C $(AUP) build/aup
 	$(MAKE) -C prototypes/connections BUILD=build/tree STORE=store-tree.c build/tree/pure-connections
 	python3 tests/check.py --aup $(AUP)/build/aup --aup-root $(AUP)
diff --git a/core/module-text.c b/core/module-text.c
index 68c03b0..4ce1769 100644
--- a/core/module-text.c
+++ b/core/module-text.c
@@ -72,6 +72,35 @@ size_t pure_format_number(double value, char out[32]) {
         out[length] = 0;
         return length;
     }
+    /* The fewest decimals k whose digits, the value * 10^k rounded to the
+     * nearest integer m, read back as the value: %g's digits. m and 10^k are
+     * exact (m < 2^53, k <= 22), so m / 10^k is the correctly rounded value of
+     * those digits, as reading them is; the product's rounding error (fma)
+     * tells the nearest m, and a half between two goes to %g. Else, and for
+     * what %g would write with an exponent, the fewest %g digits that do. */
+    static const double tens[] = {1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11,
+                                  1e12, 1e13, 1e14, 1e15, 1e16};
+    double magnitude = fabs(value);
+    if (magnitude >= 1e-4 && magnitude < 1e15)
+        for (int k = 1; k <= 16; ++k) {
+            double scaled = magnitude * tens[k];
+            if (scaled >= 9007199254740992.0) break;
+            uint64_t m = (uint64_t)scaled;
+            double beyond = (scaled - (double)m) + fma(magnitude, tens[k], -scaled);
+            if (beyond == 0.5) break;
+            m += beyond > 0.5;
+            if ((double)m / tens[k] != magnitude) continue;
+            char digits[24];
+            size_t count = 0, length = 0;
+            do { digits[count++] = (char)('0' + m % 10); m /= 10; } while (m);
+            while (count <= (size_t)k) digits[count++] = '0'; /* a leading zero before the point */
+            if (value < 0) out[length++] = '-';
+            while (count > (size_t)k) out[length++] = digits[--count];
+            out[length++] = '.';
+            while (count) out[length++] = digits[--count];
+            out[length] = 0;
+            return length;
+        }
     for (int precision = 1; precision <= 17; ++precision) {
         snprintf(out, 32, "%.*g", precision, value);
         if (strtod(out, NULL) == value) break;
diff --git a/core/test-module.c b/core/test-module.c
index d5b7eee..8ca031c 100644
--- a/core/test-module.c
+++ b/core/test-module.c
@@ -3987,7 +3987,43 @@ static void timers(void) {
     CHECK(load(bad, &module_b) == PURE_INVALID && strstr(error, "view with an effect"));
 }
 
+/* A number is written with the fewest digits that read back as it, the
+ * digits %g gives (with its exponent's leading zeros dropped): the same
+ * text as that search over precisions, for values of every kind. */
+static void number_texts(void) {
+    uint32_t seed = 0x6e756d73u;
+    size_t checked = 0;
+    for (unsigned i = 0; i < 300000; ++i) {
+        double value;
+        uint64_t bits = (uint64_t)next_random(&seed) << 32 | next_random(&seed);
+        switch (i % 5) {
+        case 0: value = (double)((int64_t)(bits % 2000001) - 1000000) / 100.0; break;         /* money */
+        case 1: value = (double)(bits % 20000001) / pow(10, (double)(bits >> 40 & 7)); break; /* decimals */
+        case 2: memcpy(&value, &bits, sizeof(value)); if (!isfinite(value)) value = 0.75; break;
+        case 3: value = (double)(bits % 100000) / 7.0; break;
+        default: value = ldexp((double)(bits >> 11), (int)(bits % 80) - 60); break;
+        }
+        char got[32], want[32];
+        pure_format_number(value, got);
+        if (value == 0) strcpy(want, "0");
+        else if (value == floor(value) && fabs(value) < 1e15) snprintf(want, sizeof(want), "%.0f", value);
+        else {
+            for (int precision = 1; precision <= 17; ++precision) {
+                snprintf(want, sizeof(want), "%.*g", precision, value);
+                if (strtod(want, NULL) == value) break;
+            }
+            char *exponent = strchr(want, 'e');
+            if (exponent)
+                while (exponent[2] == '0' && exponent[3]) memmove(exponent + 2, exponent + 3, strlen(exponent + 2));
+        }
+        CHECK(!strcmp(got, want) && strtod(got, NULL) == value);
+        ++checked;
+    }
+    printf("{\"number_texts\":%zu}\n", checked);
+}
+
 int main(void) {
+    number_texts();
     behaviour();
     faults();
     rejections();
diff --git a/serve/message-json.c b/serve/message-json.c
index 8c288ce..643cc28 100644
--- a/serve/message-json.c
+++ b/serve/message-json.c
@@ -46,6 +46,36 @@ static size_t json_number(const unsigned char *p, size_t size, bool *integer) {
 static bool json_double(const unsigned char *text, size_t size, double *value) {
     bool integer;
     if (!size || json_number(text, size, &integer) != size || size > 40) return false;
+    /* Digits m < 2^53 and a power of ten 10^e, |e| <= 22, are exact, so one
+     * multiplication or division rounds m * 10^e correctly (Clinger's fast
+     * path); other numbers go to strtod. */
+    static const double tens[] = {1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11,
+                                  1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22};
+    const unsigned char *p = text, *end = text + size;
+    bool negative = *p == '-', exact = true;
+    uint64_t m = 0;
+    int exponent = 0;
+    p += negative;
+    for (; p < end && *p >= '0' && *p <= '9'; ++p) {
+        exact = exact && m < UINT64_C(900719925474099);
+        m = m * 10 + (uint64_t)(*p - '0');
+    }
+    if (p < end && *p == '.')
+        for (++p; p < end && *p >= '0' && *p <= '9'; ++p, --exponent) {
+            exact = exact && m < UINT64_C(900719925474099);
+            m = m * 10 + (uint64_t)(*p - '0');
+        }
+    if (p < end) { /* an exponent */
+        bool below = *++p == '-';
+        int e = 0;
+        for (p += *p == '-' || *p == '+'; p < end; ++p) e = e < 1000 ? e * 10 + (*p - '0') : e;
+        exponent += below ? -e : e;
+    }
+    if (exact && exponent >= -22 && exponent <= 22) {
+        double v = exponent < 0 ? (double)m / tens[-exponent] : (double)m * tens[exponent];
+        *value = negative ? -v : v;
+        return true;
+    }
     char digits[48];
     memcpy(digits, text, size);
     digits[size] = 0;
diff --git a/serve/test-message-json.c b/serve/test-message-json.c
new file mode 100644
index 0000000..618a363
--- /dev/null
+++ b/serve/test-message-json.c
@@ -0,0 +1,68 @@
+/* JSON numbers are read as strtod reads them, bit for bit: in the forms a
+ * client writes (a program's shortest digits, fixed decimals, exponents), at
+ * the edges of the exact path (2^53, 10^22) and past them, and refused when
+ * not finite or not JSON. */
+#include "message-json.h"
+#include <math.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#define CHECK(x) do { if (!(x)) { fprintf(stderr, "%s:%d: %s\n", __FILE__, __LINE__, #x); exit(1); } } while (0)
+
+static uint64_t next_random(uint64_t *state) {
+    uint64_t x = *state;
+    x ^= x << 13, x ^= x >> 7, x ^= x << 17;
+    return *state = x;
+}
+
+/* The text read as a number field: true with its value, or false. */
+static bool read_number(const char *text, double *value) {
+    static const PureField number = {.type = PURE_FIELD_NUMBER};
+    unsigned char prefix[8];
+    if (!pure_json_fixed(&number, (const unsigned char *)text, strlen(text), prefix)) return false;
+    memcpy(value, prefix, sizeof(*value));
+    return true;
+}
+
+static void same_as_strtod(const char *text) {
+    double got, want = strtod(text, NULL);
+    CHECK(read_number(text, &got) && !memcmp(&got, &want, sizeof(got)));
+}
+
+int main(void) {
+    static const char *const edges[] = {
+        "0", "-0", "0.0", "-0.0", "1", "0.1", "0.2", "0.30000000000000004", "123456.789", "1e22", "1e23", "1E-22",
+        "1e-23", "9007199254740991", "9007199254740992", "9007199254740993", "900719925474099.3", "12345678901234567890",
+        "0.000000000000000000000001", "1.7976931348623157e308", "4.9406564584124654e-324", "2.2250738585072014e-308",
+        "-123.456e-7", "5e+0", "0.5e1", "1234567890.0987654321"};
+    for (size_t i = 0; i < sizeof(edges) / sizeof(*edges); ++i) same_as_strtod(edges[i]);
+    static const char *const refused[] = {"1e309", "-1e400", "01", "1.", ".5", "+1", "1e", "--1", "0x10", "1.5e+", "NaN"};
+    for (size_t i = 0; i < sizeof(refused) / sizeof(*refused); ++i) {
+        double value;
+        CHECK(!read_number(refused[i], &value));
+    }
+    uint64_t seed = 0x6a736f6eu;
+    size_t read = 0;
+    for (unsigned i = 0; i < 400000; ++i) {
+        uint64_t bits = next_random(&seed);
+        double value;
+        switch (i % 4) {
+        case 0: value = (double)(int64_t)(bits % 2000001 - 1000000) / 100.0; break;
+        case 1: value = (double)(bits % 100000000) / pow(10, (double)(bits >> 50 & 15)); break;
+        case 2: memcpy(&value, &bits, sizeof(value)); if (!isfinite(value)) value = -2.5; break;
+        default: value = ldexp((double)(bits >> 11), (int)(bits % 120) - 90); break;
+        }
+        char text[64];
+        snprintf(text, sizeof(text), "%.17g", value);
+        same_as_strtod(text);
+        snprintf(text, sizeof(text), "%.*f", (int)(bits >> 56) % 12, value);
+        if (strlen(text) < 40) same_as_strtod(text);
+        snprintf(text, sizeof(text), "%.*e", (int)(bits >> 58) % 17, value);
+        same_as_strtod(text);
+        read += 3;
+    }
+    printf("{\"json_numbers\":%zu}\n", read);
+    return 0;
+}