From ebd1a21af1d0260c4b7e0494c7c882f6d9cbc6c6 Mon Sep 17 00:00:00 2001 From: David Robillard Date: Sun, 29 Sep 2019 22:39:51 +0200 Subject: Add minimal soft floating point implementation --- src/soft_float.c | 161 ++++++++++++++++++++++++++++++++++++++++++++++++ src/soft_float.h | 66 ++++++++++++++++++++ tests/soft_float_test.c | 119 +++++++++++++++++++++++++++++++++++ tests/test_data.h | 56 +++++++++++++++++ wscript | 3 + 5 files changed, 405 insertions(+) create mode 100644 src/soft_float.c create mode 100644 src/soft_float.h create mode 100644 tests/soft_float_test.c create mode 100644 tests/test_data.h diff --git a/src/soft_float.c b/src/soft_float.c new file mode 100644 index 00000000..3738490f --- /dev/null +++ b/src/soft_float.c @@ -0,0 +1,161 @@ +/* + Copyright 2019 David Robillard + + 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 +#include +#include + +/// 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..58627bbe --- /dev/null +++ b/src/soft_float.h @@ -0,0 +1,66 @@ +/* + Copyright 2019 David Robillard + + 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 + +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..58b4dac0 --- /dev/null +++ b/tests/soft_float_test.c @@ -0,0 +1,119 @@ +/* + Copyright 2011-2019 David Robillard + + 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 +#include +#include +#include + +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..ba8b860b --- /dev/null +++ b/tests/test_data.h @@ -0,0 +1,56 @@ +/* + Copyright 2019 David Robillard + + 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 +#include + +/// 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; + 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; + memcpy(&d, &rep, sizeof(d)); + return d; +} diff --git a/wscript b/wscript index f17f2547..db0c76a6 100644 --- a/wscript +++ b/wscript @@ -122,6 +122,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', @@ -202,6 +203,7 @@ def build(bld): ('sink_test', 'tests/sink_test.c'), ('serd_test', 'tests/serd_test.c'), ('read_chunk_test', 'tests/read_chunk_test.c'), + ('soft_float_test', 'tests/soft_float_test.c'), ('terse_write_test', 'tests/terse_write_test.c'), ('nodes_test', 'tests/nodes_test.c'), ('overflow_test', 'tests/overflow_test.c'), @@ -559,6 +561,7 @@ def test(tst): check(['./nodes_test']) check(['./overflow_test']) check(['./serd_test']) + check(['./soft_float_test']) check(['./terse_write_test']) check(['./read_chunk_test']) -- cgit v1.2.1