aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorDavid Robillard <d@drobilla.net>2019-09-29 22:39:51 +0200
committerDavid Robillard <d@drobilla.net>2020-10-27 13:13:58 +0100
commitfd48cffcc1838060bcc451300b38982776779c13 (patch)
tree6df1f08c1748ad255953a1b6a25a2366d8ba4bca
parent30d845380249d1de3a28e12407c5263629424c6d (diff)
downloadserd-fd48cffcc1838060bcc451300b38982776779c13.tar.gz
serd-fd48cffcc1838060bcc451300b38982776779c13.tar.bz2
serd-fd48cffcc1838060bcc451300b38982776779c13.zip
Add minimal soft floating point implementation
-rw-r--r--src/ieee_float.h2
-rw-r--r--src/soft_float.c161
-rw-r--r--src/soft_float.h66
-rw-r--r--tests/soft_float_test.c119
-rw-r--r--tests/test_data.h56
-rw-r--r--wscript5
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'])