Add minimal soft floating point implementation

6 files changed, 408 insertions, 1 deletions

diff --git a/src/ieee_float.h b/src/ieee_float.h
index adf2b6cf..a7f7cdf9 100644
--- a/src/ieee_float.h
+++ b/src/ieee_float.h
@@ -34,7 +34,7 @@ static const int      dbl_subnormal_expt    = -0x3FF - DBL_MANT_DIG + 2;
 static inline uint64_t
 double_to_rep(const double d)
 {
-	uint64_t rep;
+	uint64_t rep = 0;
 	memcpy(&rep, &d, sizeof(rep));
 	return rep;
 }
diff --git a/src/soft_float.c b/src/soft_float.c
new file mode 100644
index 00000000..6c84f964
--- /dev/null
+++ b/src/soft_float.c
@@ -0,0 +1,161 @@
+/*
+  Copyright 2019-2020 David Robillard <http://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 "soft_float.h"
+
+#include "ieee_float.h"
+#include "int_math.h"
+
+#include <assert.h>
+#include <math.h>
+#include <stdint.h>
+
+/// 10^k for k = min_dec_expt, min_dec_expt + dec_expt_step, ..., max_dec_expt
+static const SerdSoftFloat soft_pow10[] = {
+        {0xFA8FD5A0081C0288, -1220}, {0xBAAEE17FA23EBF76, -1193},
+        {0x8B16FB203055AC76, -1166}, {0xCF42894A5DCE35EA, -1140},
+        {0x9A6BB0AA55653B2D, -1113}, {0xE61ACF033D1A45DF, -1087},
+        {0xAB70FE17C79AC6CA, -1060}, {0xFF77B1FCBEBCDC4F, -1034},
+        {0xBE5691EF416BD60C, -1007}, {0x8DD01FAD907FFC3C, -980},
+        {0xD3515C2831559A83, -954},  {0x9D71AC8FADA6C9B5, -927},
+        {0xEA9C227723EE8BCB, -901},  {0xAECC49914078536D, -874},
+        {0x823C12795DB6CE57, -847},  {0xC21094364DFB5637, -821},
+        {0x9096EA6F3848984F, -794},  {0xD77485CB25823AC7, -768},
+        {0xA086CFCD97BF97F4, -741},  {0xEF340A98172AACE5, -715},
+        {0xB23867FB2A35B28E, -688},  {0x84C8D4DFD2C63F3B, -661},
+        {0xC5DD44271AD3CDBA, -635},  {0x936B9FCEBB25C996, -608},
+        {0xDBAC6C247D62A584, -582},  {0xA3AB66580D5FDAF6, -555},
+        {0xF3E2F893DEC3F126, -529},  {0xB5B5ADA8AAFF80B8, -502},
+        {0x87625F056C7C4A8B, -475},  {0xC9BCFF6034C13053, -449},
+        {0x964E858C91BA2655, -422},  {0xDFF9772470297EBD, -396},
+        {0xA6DFBD9FB8E5B88F, -369},  {0xF8A95FCF88747D94, -343},
+        {0xB94470938FA89BCF, -316},  {0x8A08F0F8BF0F156B, -289},
+        {0xCDB02555653131B6, -263},  {0x993FE2C6D07B7FAC, -236},
+        {0xE45C10C42A2B3B06, -210},  {0xAA242499697392D3, -183},
+        {0xFD87B5F28300CA0E, -157},  {0xBCE5086492111AEB, -130},
+        {0x8CBCCC096F5088CC, -103},  {0xD1B71758E219652C, -77},
+        {0x9C40000000000000, -50},   {0xE8D4A51000000000, -24},
+        {0xAD78EBC5AC620000, 3},     {0x813F3978F8940984, 30},
+        {0xC097CE7BC90715B3, 56},    {0x8F7E32CE7BEA5C70, 83},
+        {0xD5D238A4ABE98068, 109},   {0x9F4F2726179A2245, 136},
+        {0xED63A231D4C4FB27, 162},   {0xB0DE65388CC8ADA8, 189},
+        {0x83C7088E1AAB65DB, 216},   {0xC45D1DF942711D9A, 242},
+        {0x924D692CA61BE758, 269},   {0xDA01EE641A708DEA, 295},
+        {0xA26DA3999AEF774A, 322},   {0xF209787BB47D6B85, 348},
+        {0xB454E4A179DD1877, 375},   {0x865B86925B9BC5C2, 402},
+        {0xC83553C5C8965D3D, 428},   {0x952AB45CFA97A0B3, 455},
+        {0xDE469FBD99A05FE3, 481},   {0xA59BC234DB398C25, 508},
+        {0xF6C69A72A3989F5C, 534},   {0xB7DCBF5354E9BECE, 561},
+        {0x88FCF317F22241E2, 588},   {0xCC20CE9BD35C78A5, 614},
+        {0x98165AF37B2153DF, 641},   {0xE2A0B5DC971F303A, 667},
+        {0xA8D9D1535CE3B396, 694},   {0xFB9B7CD9A4A7443C, 720},
+        {0xBB764C4CA7A44410, 747},   {0x8BAB8EEFB6409C1A, 774},
+        {0xD01FEF10A657842C, 800},   {0x9B10A4E5E9913129, 827},
+        {0xE7109BFBA19C0C9D, 853},   {0xAC2820D9623BF429, 880},
+        {0x80444B5E7AA7CF85, 907},   {0xBF21E44003ACDD2D, 933},
+        {0x8E679C2F5E44FF8F, 960},   {0xD433179D9C8CB841, 986},
+        {0x9E19DB92B4E31BA9, 1013},  {0xEB96BF6EBADF77D9, 1039},
+        {0xAF87023B9BF0EE6B, 1066}};
+
+SerdSoftFloat
+soft_float_from_double(const double d)
+{
+	assert(d >= 0.0);
+
+	const uint64_t rep  = double_to_rep(d);
+	const uint64_t frac = rep & dbl_mant_mask;
+	const int      expt = (int)((rep & dbl_expt_mask) >> dbl_physical_mant_dig);
+
+	if (expt == 0) { // Subnormal
+		SerdSoftFloat v = {frac, dbl_subnormal_expt};
+		return v;
+	}
+
+	const SerdSoftFloat v = {frac + dbl_hidden_bit, expt - dbl_expt_bias};
+	return v;
+}
+
+double
+soft_float_to_double(const SerdSoftFloat v)
+{
+	return ldexp((double)v.f, v.e);
+}
+
+SerdSoftFloat
+soft_float_normalize(SerdSoftFloat value)
+{
+	const unsigned amount = serd_clz64(value.f);
+
+	value.f <<= amount;
+	value.e -= (int)amount;
+
+	return value;
+}
+
+SerdSoftFloat
+soft_float_multiply(const SerdSoftFloat lhs, const SerdSoftFloat rhs)
+{
+	static const uint64_t mask  = 0xFFFFFFFF;
+	static const uint64_t round = 1u << 31;
+
+	const uint64_t l0   = lhs.f >> 32;
+	const uint64_t l1   = lhs.f & mask;
+	const uint64_t r0   = rhs.f >> 32;
+	const uint64_t r1   = rhs.f & mask;
+	const uint64_t l0r0 = l0 * r0;
+	const uint64_t l1r0 = l1 * r0;
+	const uint64_t l0r1 = l0 * r1;
+	const uint64_t l1r1 = l1 * r1;
+	const uint64_t mid  = (l1r1 >> 32) + (l0r1 & mask) + (l1r0 & mask) + round;
+
+	const SerdSoftFloat r = {l0r0 + (l0r1 >> 32) + (l1r0 >> 32) + (mid >> 32),
+	                         lhs.e + rhs.e + 64};
+
+	return r;
+}
+
+SerdSoftFloat
+soft_float_exact_pow10(const int expt)
+{
+	static const SerdSoftFloat table[8] = {{0xA000000000000000, -60},
+	                                       {0xC800000000000000, -57},
+	                                       {0xFA00000000000000, -54},
+	                                       {0x9C40000000000000, -50},
+	                                       {0xC350000000000000, -47},
+	                                       {0xF424000000000000, -44},
+	                                       {0x9896800000000000, -40}};
+
+	assert(expt > 0);
+	assert(expt < dec_expt_step);
+	return table[expt - 1];
+}
+
+SerdSoftFloat
+soft_float_pow10_under(const int exponent, int* pow10_exponent)
+{
+	assert(exponent >= min_dec_expt);
+	assert(exponent < max_dec_expt + dec_expt_step);
+
+	const int cache_offset = -min_dec_expt;
+	const int index        = (exponent + cache_offset) / dec_expt_step;
+
+	*pow10_exponent = min_dec_expt + index * dec_expt_step;
+
+	assert(*pow10_exponent <= exponent);
+	assert(exponent < *pow10_exponent + dec_expt_step);
+
+	return soft_pow10[index];
+}
diff --git a/src/soft_float.h b/src/soft_float.h
new file mode 100644
index 00000000..d13ec69b
--- /dev/null
+++ b/src/soft_float.h
@@ -0,0 +1,66 @@
+/*
+  Copyright 2019-2020 David Robillard <http://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.
+*/
+
+#ifndef SERD_SOFT_FLOAT_H
+#define SERD_SOFT_FLOAT_H
+
+#include <stdint.h>
+
+typedef struct
+{
+	uint64_t f; ///< Significand
+	int      e; ///< Exponent
+} SerdSoftFloat;
+
+static const int min_dec_expt  = -348;
+static const int max_dec_expt  = 340;
+static const int dec_expt_step = 8;
+
+/// Convert `d` to a soft float
+SerdSoftFloat
+soft_float_from_double(double d);
+
+/// Convert `v` to a double
+double
+soft_float_to_double(SerdSoftFloat v);
+
+/// Normalize `value` so the MSb of its significand is 1
+SerdSoftFloat
+soft_float_normalize(SerdSoftFloat value);
+
+/// Multiply `lhs` by `rhs` and return the result
+SerdSoftFloat
+soft_float_multiply(SerdSoftFloat lhs, SerdSoftFloat rhs);
+
+/// Return exactly 10^e for e in [0...dec_expt_step]
+SerdSoftFloat
+soft_float_exact_pow10(int expt);
+
+/**
+   Return a cached power of 10 with exponent not greater than `max_exponent`.
+
+   Valid only for `max_exponent` values from min_dec_expt to max_dec_expt +
+   dec_expt_step.  The returned power's exponent is a multiple of
+   dec_expt_step.
+
+   @param max_exponent Maximum decimal exponent of the result.
+   @param[out] pow10_exponent Set to the decimal exponent of the result.
+   @return A cached power of 10 as a soft float.
+*/
+SerdSoftFloat
+soft_float_pow10_under(int max_exponent, int* pow10_exponent);
+
+#endif // SERD_SOFT_FLOAT_H
diff --git a/tests/soft_float_test.c b/tests/soft_float_test.c
new file mode 100644
index 00000000..85e9e4b6
--- /dev/null
+++ b/tests/soft_float_test.c
@@ -0,0 +1,119 @@
+/*
+  Copyright 2011-2020 David Robillard <http://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.
+*/
+
+#undef NDEBUG
+
+#include "test_data.h"
+
+#include "../src/ieee_float.h"
+#include "../src/soft_float.h"
+
+#include <assert.h>
+#include <math.h>
+#include <stdbool.h>
+#include <stdint.h>
+
+static uint64_t
+ulp_distance(const double a, const double b)
+{
+	assert(a >= 0.0);
+	assert(b >= 0.0);
+
+	if (a == b) {
+		return 0;
+	} else if (isnan(a) || isnan(b)) {
+		return UINT64_MAX;
+	} else if (isinf(a) || isinf(b)) {
+		return UINT64_MAX;
+	}
+
+	const uint64_t ia = double_to_rep(a);
+	const uint64_t ib = double_to_rep(b);
+
+	return ia > ib ? ia - ib : ib - ia;
+}
+
+static bool
+check_multiply(const double lhs, const double rhs)
+{
+	assert(lhs >= 0.0);
+	assert(rhs >= 0.0);
+
+	const SerdSoftFloat sl = soft_float_normalize(soft_float_from_double(lhs));
+	const SerdSoftFloat sr = soft_float_normalize(soft_float_from_double(rhs));
+	const SerdSoftFloat sp = soft_float_multiply(sl, sr);
+	const double        ep = lhs * rhs;
+	const double        dp = soft_float_to_double(sp);
+
+	return ulp_distance(dp, ep) <= 1;
+}
+
+static void
+test_multiply(void)
+{
+	assert(check_multiply(1.0, 1.0));
+	assert(check_multiply(1.0, 8.0));
+	assert(check_multiply(8.0, 1.0));
+	assert(check_multiply(2.0, 4.0));
+	assert(check_multiply(1e100, 1e-100));
+
+	uint64_t seed = 1;
+	for (int i = 0; i < 1000000; ++i) {
+		const double l = fabs(double_from_rep(seed = lcg64(seed)));
+		const double r = fabs(double_from_rep(seed = lcg64(seed)));
+		if (isfinite(l) && isfinite(r)) {
+			assert(check_multiply(l, r));
+		}
+	}
+}
+
+static void
+test_exact_pow10(void)
+{
+	for (int i = 1; i < dec_expt_step; ++i) {
+		const SerdSoftFloat power = soft_float_exact_pow10(i);
+
+		const double d = soft_float_to_double(power);
+		const double p = pow(10, i);
+		assert(ulp_distance(d, p) <= 1);
+	}
+}
+
+static void
+test_pow10_under(void)
+{
+	for (int i = min_dec_expt; i < max_dec_expt + dec_expt_step; ++i) {
+		int                 expt10 = 0;
+		const SerdSoftFloat power  = soft_float_pow10_under(i, &expt10);
+
+		assert(expt10 <= i);
+		assert(i - expt10 < dec_expt_step);
+
+		const double d = soft_float_to_double(power);
+		const double p = pow(10, expt10);
+		assert(ulp_distance(d, p) <= 1);
+	}
+}
+
+int
+main(void)
+{
+	test_multiply();
+	test_exact_pow10();
+	test_pow10_under();
+
+	return 0;
+}
diff --git a/tests/test_data.h b/tests/test_data.h
new file mode 100644
index 00000000..fa32f529
--- /dev/null
+++ b/tests/test_data.h
@@ -0,0 +1,56 @@
+/*
+  Copyright 2019 David Robillard <http://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 <stdint.h>
+#include <string.h>
+
+/// Linear Congruential Generator for making random floats
+static inline uint32_t
+lcg32(const uint32_t i)
+{
+	static const uint32_t a = 134775813u;
+	static const uint32_t c = 1u;
+
+	return (a * i) + c;
+}
+
+/// Linear Congruential Generator for making random doubles
+static inline uint64_t
+lcg64(const uint64_t i)
+{
+	static const uint64_t a = 6364136223846793005ull;
+	static const uint64_t c = 1ull;
+
+	return (a * i) + c;
+}
+
+/// Return the float with representation `rep`
+static inline float
+float_from_rep(const uint32_t rep)
+{
+	float f = 0.0f;
+	memcpy(&f, &rep, sizeof(f));
+	return f;
+}
+
+/// Return the double with representation `rep`
+static inline double
+double_from_rep(const uint64_t rep)
+{
+	double d = 0.0;
+	memcpy(&d, &rep, sizeof(d));
+	return d;
+}
diff --git a/wscript b/wscript
index b43626d0..2943731d 100644
--- a/wscript
+++ b/wscript
@@ -61,6 +61,7 @@ def configure(conf):
                 '-Wno-covered-switch-default',
                 '-Wno-disabled-macro-expansion',
                 '-Wno-double-promotion',
+                '-Wno-float-equal',
                 '-Wno-format-nonliteral',
                 '-Wno-implicit-fallthrough',
                 '-Wno-padded',
@@ -68,6 +69,7 @@ def configure(conf):
             ],
             'gcc': [
                 '-Wno-cast-align',
+                '-Wno-float-equal',
                 '-Wno-padded',
             ],
             'msvc': [
@@ -181,6 +183,7 @@ lib_source = ['src/base64.c',
               'src/range.c',
               'src/reader.c',
               'src/sink.c',
+              'src/soft_float.c',
               'src/statement.c',
               'src/string.c',
               'src/syntax.c',
@@ -267,6 +270,7 @@ def build(bld):
                      ('read_chunk_test', 'tests/read_chunk_test.c'),
                      ('serd_test', 'tests/serd_test.c'),
                      ('sink_test', 'tests/sink_test.c'),
+                     ('soft_float_test', 'tests/soft_float_test.c'),
                      ('statement_test', 'tests/statement_test.c'),
                      ('terse_write_test', 'tests/terse_write_test.c')]:
             bld(features     = 'c cprogram',
@@ -699,6 +703,7 @@ def test(tst):
         check(['./read_chunk_test'])
         check(['./serd_test'])
         check(['./sink_test'])
+        check(['./soft_float_test'])
         check(['./statement_test'])
         check(['./terse_write_test'])