1 files changed, 243 insertions, 0 deletions
diff --git a/subprojects/exess/src/digits.c b/subprojects/exess/src/digits.c
new file mode 100644
index 00000000..dd303269
--- /dev/null
+++ b/subprojects/exess/src/digits.c
@@ -0,0 +1,243 @@
+/*
+  Copyright 2019-2021 David Robillard <d@drobilla.net>
+
+  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.
+
+  THIS 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.
+*/
+
+#include "digits.h"
+
+#include "bigint.h"
+#include "ieee_float.h"
+#include "soft_float.h"
+#include "warnings.h"
+
+#include <assert.h>
+#include <math.h>
+#include <stdbool.h>
+#include <stddef.h>
+#include <stdint.h>
+
+/*
+  This is more or less just an implementation of the classic rational number
+  based floating point print routine ("Dragon4").  See "How to Print
+  Floating-Point Numbers Accurately" by Guy L. Steele Jr. and Jon L White for
+  the canonical source.  The basic idea is to find a big rational between 1 and
+  10 where value = (numer / denom) * 10^e, then continuously divide it to
+  generate decimal digits.
+
+  Unfortunately, this algorithm requires pretty massive bigints to work
+  correctly for all doubles, and isn't particularly fast.  Something like
+  Grisu3 could be added to improve performance, but that has the annoying
+  property of needing a more precise fallback in some cases, meaning it would
+  only add more code, not replace any.  Since this is already a pretty
+  ridiculous amount of code, I'll hold off on this until it becomes a problem,
+  or somebody comes up with a better algorithm.
+*/
+
+/// Return true if the number is within the lower boundary
+static bool
+within_lower(const ExessBigint* const numer,
+             const ExessBigint* const d_lower,
+             const bool               is_even)
+{
+  return is_even ? exess_bigint_compare(numer, d_lower) <= 0
+                 : exess_bigint_compare(numer, d_lower) < 0;
+}
+
+/// Return true if the number is within the upper boundary
+static bool
+within_upper(const ExessBigint* const numer,
+             const ExessBigint* const denom,
+             const ExessBigint* const d_upper,
+             const bool               is_even)
+{
+  return is_even ? exess_bigint_plus_compare(numer, d_upper, denom) >= 0
+                 : exess_bigint_plus_compare(numer, d_upper, denom) > 0;
+}
+
+/**
+   Find values so that 0.1 <= numer/denom < 1 or 1 <= numer/denom < 10.
+
+   @param significand Double significand.
+   @param exponent Double exponent (base 2).
+   @param decimal_power Decimal exponent (log10 of the double).
+   @param[out] numer Numerator of rational number.
+   @param[out] denom Denominator of rational number.
+   @param[out] delta Distance to the lower and upper boundaries.
+*/
+static void
+calculate_initial_values(const uint64_t     significand,
+                         const int          exponent,
+                         const int          decimal_power,
+                         const bool         lower_is_closer,
+                         ExessBigint* const numer,
+                         ExessBigint* const denom,
+                         ExessBigint* const delta)
+{
+  /* Use a common denominator of 2^1 so that boundary distance is an integer.
+     If the lower boundary is closer, we need to scale everything but the
+     lower boundary to compensate, so add another factor of two here (this is
+     faster than shifting them again later as in the paper). */
+  const unsigned lg_denom = 1u + lower_is_closer;
+
+  if (exponent >= 0) {
+    // delta = 2^e
+    exess_bigint_set_u32(delta, 1);
+    exess_bigint_shift_left(delta, (unsigned)exponent);
+
+    // numer = f * 2^e
+    exess_bigint_set_u64(numer, significand);
+    exess_bigint_shift_left(numer, (unsigned)exponent + lg_denom);
+
+    // denom = 10^d
+    exess_bigint_set_pow10(denom, (unsigned)decimal_power);
+    exess_bigint_shift_left(denom, lg_denom);
+  } else if (decimal_power >= 0) {
+    // delta = 2^e, which is just 1 here since 2^-e is in the denominator
+    exess_bigint_set_u32(delta, 1);
+
+    // numer = f
+    exess_bigint_set_u64(numer, significand);
+    exess_bigint_shift_left(numer, lg_denom);
+
+    // denom = 10^d * 2^-e
+    exess_bigint_set_pow10(denom, (unsigned)decimal_power);
+    exess_bigint_shift_left(denom, (unsigned)-exponent + lg_denom);
+  } else {
+    // delta = 10^d
+    exess_bigint_set_pow10(delta, (unsigned)-decimal_power);
+
+    // numer = f * 10^-d
+    exess_bigint_set(numer, delta);
+    exess_bigint_multiply_u64(numer, significand);
+    exess_bigint_shift_left(numer, lg_denom);
+
+    // denom = 2^-exponent
+    exess_bigint_set_u32(denom, 1);
+    exess_bigint_shift_left(denom, (unsigned)-exponent + lg_denom);
+  }
+}
+
+#ifndef NDEBUG
+static bool
+check_initial_values(const ExessBigint* const numer,
+                     const ExessBigint* const denom,
+                     const ExessBigint* const d_upper)
+{
+  ExessBigint upper = *numer;
+  exess_bigint_add(&upper, d_upper);
+  assert(exess_bigint_compare(&upper, denom) >= 0);
+
+  const uint32_t div = exess_bigint_divmod(&upper, denom);
+  assert(div >= 1 && div < 10);
+  return true;
+}
+#endif
+
+static unsigned
+emit_digits(ExessBigint* const       numer,
+            const ExessBigint* const denom,
+            ExessBigint* const       d_lower,
+            ExessBigint* const       d_upper,
+            const bool               is_even,
+            char* const              buffer,
+            const size_t             max_digits)
+{
+  unsigned length = 0;
+  for (size_t i = 0; i < max_digits; ++i) {
+    // Emit the next digit
+    const uint32_t digit = exess_bigint_divmod(numer, denom);
+    assert(digit <= 9);
+    buffer[length++] = (char)('0' + digit);
+
+    // Check for termination
+    const bool within_low  = within_lower(numer, d_lower, is_even);
+    const bool within_high = within_upper(numer, denom, d_upper, is_even);
+    if (!within_low && !within_high) {
+      exess_bigint_multiply_u32(numer, 10);
+      exess_bigint_multiply_u32(d_lower, 10);
+      if (d_lower != d_upper) {
+        exess_bigint_multiply_u32(d_upper, 10);
+      }
+    } else {
+      if (!within_low || (within_high && exess_bigint_plus_compare(
+                                           numer, numer, denom) >= 0)) {
+        // In high only, or halfway and the next digit is > 5, round up
+        assert(buffer[length - 1] != '9');
+        buffer[length - 1]++;
+      }
+
+      break;
+    }
+  }
+
+  return length;
+}
+
+ExessDigitCount
+exess_digits(const double d, char* const buf, const unsigned max_digits)
+{
+  EXESS_DISABLE_CONVERSION_WARNINGS
+  assert(isfinite(d) && fpclassify(d) != FP_ZERO);
+  EXESS_RESTORE_WARNINGS
+
+  const ExessSoftFloat value           = soft_float_from_double(d);
+  const int            power           = (int)lrint(log10(d));
+  const bool           is_even         = !(value.f & 1u);
+  const bool           lower_is_closer = double_lower_boundary_is_closer(d);
+
+  // Calculate initial values so that v = (numer / denom) * 10^power
+  ExessBigint numer;
+  ExessBigint denom;
+  ExessBigint d_lower;
+  calculate_initial_values(
+    value.f, value.e, power, lower_is_closer, &numer, &denom, &d_lower);
+
+  ExessBigint  d_upper_storage;
+  ExessBigint* d_upper = NULL;
+  if (lower_is_closer) {
+    // Scale upper boundary to account for the closer lower boundary
+    // (the numerator and denominator were already scaled above)
+    d_upper_storage = d_lower;
+    d_upper         = &d_upper_storage;
+    exess_bigint_shift_left(d_upper, 1);
+  } else {
+    d_upper = &d_lower; // Boundaries are the same, reuse the lower
+  }
+
+  // Scale if necessary to make 1 <= (numer + delta) / denom < 10
+  ExessDigitCount count = {0, 0};
+  if (within_upper(&numer, &denom, d_upper, is_even)) {
+    count.expt = power;
+  } else {
+    count.expt = power - 1;
+    exess_bigint_multiply_u32(&numer, 10);
+    exess_bigint_multiply_u32(&d_lower, 10);
+    if (d_upper != &d_lower) {
+      exess_bigint_multiply_u32(d_upper, 10);
+    }
+  }
+
+  // Write digits to output
+  assert(check_initial_values(&numer, &denom, d_upper));
+  count.count =
+    emit_digits(&numer, &denom, &d_lower, d_upper, is_even, buf, max_digits);
+
+  // Trim trailing zeros
+  while (count.count > 1 && buf[count.count - 1] == '0') {
+    buf[--count.count] = 0;
+  }
+
+  buf[count.count] = '\0';
+  return count;
+}