diff --git a/lib/std/core/floatparse.c3 b/lib/std/core/floatparse.c3 new file mode 100644 index 000000000..da0bb3427 --- /dev/null +++ b/lib/std/core/floatparse.c3 @@ -0,0 +1,488 @@ +module std::core::str; +import std::math; + +// Float parsing based on code in Musl floatscan.c by Rich Felker. +// Musl uses the MIT license, copied below: +// ---------------------------------------------------------------------- +// Copyright © 2005-2014 Rich Felker, et al. +// +// Permission is hereby granted, free of charge, to any person obtaining +// a copy of this software and associated documentation files (the +// "Software"), to deal in the Software without restriction, including +// without limitation the rights to use, copy, modify, merge, publish, +// distribute, sublicense, and/or sell copies of the Software, and to +// permit persons to whom the Software is furnished to do so, subject to +// the following conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +// IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +// CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +// TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +// SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +// ---------------------------------------------------------------------- + +const KMAX = 128; +const MASK = KMAX - 1; +const B1B_DIG = 2; +const uint[2] B1B_MAX = { 9007199, 254740991 }; + +/** + * @require chars.len > 0 + **/ +macro double! decfloat(char[] chars, int $bits, int $emin, int sign) +{ + uint[KMAX] x; + const uint[2] TH = B1B_MAX; + int emax = - $emin - $bits + 3; + + const int[*] P10S = { 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000 }; + usz index = 0; + bool got_digit = chars[0] == '0'; + bool got_rad; + usz len = chars.len; + usz last_char = len - 1; + long lrp = 0; + + assert(len); + + char c = void; + // Skip past first characters + while ((c = chars[index]) == '0') + { + if (index == last_char) return sign * 0.0; + index++; + } + + if (c == '.') + { + got_rad = true; + if (index == last_char) + { + if (!got_digit) return NumberConversion.MALFORMED_FLOAT!; + return sign * 0.0; + } + if (index != last_char && (c = chars[++index]) == '0') + { + lrp--; + got_digit = true; + while (last_char != index && (c = chars[++index]) == '0') + { + lrp--; + } + } + } + + long dc = 0; + int k; + int j; + int lnz; + while (c - '0' < 10u || c == '.') + { + switch + { + case c == '.': + if (got_rad) return NumberConversion.MALFORMED_FLOAT!; + got_rad = true; + lrp = dc; + case k < KMAX - 3: + dc++; + if (c != '0') lnz = (int)dc; + if (j) + { + x[k] = x[k] * 10 + c - '0'; + } + else + { + x[k] = c - '0'; + } + if (++j == 9) + { + k++; + j = 0; + } + got_digit = true; + default: + dc++; + if (c != '0') x[KMAX - 4] |= 1; + + } + if (index == last_char) break; + assert(index < last_char); + c = chars[++index]; + } + if (!got_rad) lrp = dc; + if (!got_digit) return NumberConversion.MALFORMED_FLOAT!; + if ((c | 32) == 'e') + { + if (last_char == index) return NumberConversion.MALFORMED_FLOAT!; + long! e10 = str::to_long((String)chars[index + 1..]); + if (catch e10) return NumberConversion.MALFORMED_FLOAT!; + lrp += e10; + } + // Handle zero specially to avoid nasty special cases later + if (!x[0]) return sign * 0.0; + + // Optimize small integers (w/no exponent) and over/under-flow + if (lrp == dc && dc < 10 && ($bits > 30 || (ulong)x[0] >> $bits == 0)) return sign * (double)x[0]; + if (lrp > - $emin / 2) return NumberConversion.FLOAT_OUT_OF_RANGE!; + if (lrp < $emin - 2 * math::DOUBLE_MANT_DIG) return NumberConversion.FLOAT_OUT_OF_RANGE!; + + // Align incomplete final B1B digit + if (j) + { + for (; j < 9; j++) x[k] *= 10; + k++; + j = 0; + } + + int a; + int z = k; + int e2; + long rp = lrp; + + // Optimize small to mid-size integers (even in exp. notation) + if (lnz < 9 && lnz <= rp && rp < 18) + { + if (rp == 9) return sign * (double)x[0]; + if (rp < 9) return sign * (double)x[0] / P10S[8 - rp]; + int bitlim = $bits - 3 * (int)(rp - 9); + if (bitlim > 30 || x[0] >> bitlim == 0) return sign * (double)x[0] * P10S[rp - 10]; + } + + // Align radix point to B1B digit boundary + if (rp % 9) + { + long rpm9 = rp >= 0 ? rp % 9 : rp % 9 + 9; + int p10 = P10S[8 - rpm9]; + uint carry = 0; + for (k = a; k != z; k++) + { + uint tmp = x[k] % p10; + x[k] = x[k] / p10 + carry; + carry = 1000000000 / p10 * tmp; + if (k == a && !x[k]) + { + a = (a + 1) & MASK; + rp -= 9; + } + } + if (carry) x[z++] = carry; + rp += 9 - rpm9; + } + + // Upscale until desired number of bits are left of radix point + while (rp < 9 * B1B_DIG || (rp == 9 * B1B_DIG && x[a] < TH[0])) + { + uint carry = 0; + e2 -= 29; + for (k = (z - 1) & MASK; ; k = (k - 1) & MASK) + { + ulong tmp = (ulong)x[k] << 29 + carry; + if (tmp > 1000000000) + { + carry = (uint)(tmp / 1000000000); + x[k] = (uint)(tmp % 1000000000); + } + else + { + carry = 0; + x[k] = (uint)tmp; + } + if (k == (z - 1) & MASK && k != a && !x[k]) z = k; + if (k == a) break; + } + if (carry) + { + rp += 9; + a = (a - 1) & MASK; + if (a == z) + { + z = (z - 1) & MASK; + x[(z - 1) & MASK] |= x[z]; + } + x[a] = carry; + } + } + + // Downscale until exactly number of bits are left of radix point + while (true) + { + uint carry = 0; + int sh = 1; + int i; + for (i = 0; i < B1B_DIG; i++) + { + k = (a + i) & MASK; + if (k == z || x[k] < TH[i]) + { + i = B1B_DIG; + break; + } + if (x[(a + i) & MASK] > TH[i]) break; + } + if (i == B1B_DIG && rp == 9 * B1B_DIG) break; + if (rp > 9 + 9 * B1B_DIG) sh = 9; + e2 += sh; + for (k = a; k != z; k = (k+1) & MASK) + { + uint tmp = x[k] & (1 << sh - 1); + x[k] = x[k] >> sh + carry; + carry = (1000000000 >> sh) * tmp; + if (k == a && !x[k]) + { + a = (a + 1) & MASK; + i--; + rp -= 9; + } + } + if (carry) + { + if ((z + 1) & MASK != a) + { + x[z] = carry; + z = (z + 1) & MASK; + } + else + { + x[(z - 1) & MASK] |= 1; + } + } + } + + // Assemble desired bits into floating point variable + double y; + int i; + for (i = 0; i < B1B_DIG; i++) + { + if ((a + i) & MASK == z) x[(z = (z + 1) & MASK) - 1] = 0; + y = 1000000000.0 * y + x[(a + i) & MASK]; + } + + y *= sign; + + bool denormal; + /* Limit precision for denormal results */ + uint bits = $bits; + if (bits > math::DOUBLE_MANT_DIG + e2 - $emin) + { + bits = math::DOUBLE_MANT_DIG + e2 - $emin; + if (bits < 0) bits = 0; + denormal = true; + } + + // Calculate bias term to force rounding, move out lower bits + double bias; + double frac; + if (bits < math::DOUBLE_MANT_DIG) + { + bias = math::copysign(math::scalbn(1, 2 * math::DOUBLE_MANT_DIG - bits - 1), y); + frac = y % math::scalbn(1, math::DOUBLE_MANT_DIG - bits); + y -= frac; + y += bias; + } + + // Process tail of decimal input so it can affect rounding + if ((a + i) & MASK != z) + { + uint t = x[(a + i) & MASK]; + switch + { + case t < 500000000 && (t || (a + i + 1) & MASK != z): + frac += 0.25 * sign; + case t > 500000000: + frac += 0.75 * sign; + case t == 500000000: + if ((a + i + 1) & MASK == z) + { + frac += 0.5 * sign; + } + else + { + frac += 0.75 * sign; + } + } + if (math::DOUBLE_MANT_DIG - bits >= 2 && !(frac % 1)) frac++; + } + + y += frac; + y -= bias; + + if (((e2 + math::DOUBLE_MANT_DIG) & int.max) > emax - 5) + { + if (math::abs(y) >= 0x1p53) + { + if (denormal && bits == math::DOUBLE_MANT_DIG + e2 - $emin) denormal = false; + y *= 0.5; + e2++; + } + if (e2 + math::DOUBLE_MANT_DIG > emax || (denormal && frac)) return NumberConversion.MALFORMED_FLOAT!; + } + return math::scalbn(y, e2); +} + +macro double! hexfloat(char[] chars, int $bits, int $emin, int sign) +{ + double scale = 1; + uint x; + long rp; + long dc; + char c = void; + bool got_rad; + bool got_digit; + bool got_tail; + usz len = chars.len; + usz last_char = len - 1; + usz index; + double y; + + // Skip past first characters + while ((c = chars[index]) == '0') + { + if (index == last_char) return 0.0; + index++; + } + if (c == '.') + { + got_rad = true; + if (index == last_char) + { + if (!got_digit) return NumberConversion.MALFORMED_FLOAT!; + return sign * 0.0; + } + if (index != last_char && (c = chars[++index]) == '0') + { + rp--; + got_digit = true; + while (last_char != index && (c = chars[++index]) == '0') + { + rp--; + } + } + } + + while ((c - '0') < 10u || ((c | 32) - 'a') < 6u || c == '.') + { + if (c == '.') + { + if (got_rad) return NumberConversion.MALFORMED_FLOAT!; + got_rad = true; + rp = dc; + } + else + { + got_digit = true; + int d = {| + if (c > '9') return (c | 32) + 10 - 'a'; + return c - '0'; + |}; + switch + { + case dc < 8: + x = x * 16 + d; + case dc < math::DOUBLE_MANT_DIG / 4 + 1: + y += d * (scale /= 16); + got_tail = true; + case d && !got_tail: + y += 0.5 * scale; + got_tail = true; + } + dc++; + } + if (index == last_char) break; + c = chars[++index]; + } + if (!got_digit) return NumberConversion.MALFORMED_FLOAT!; + if (!got_rad) rp = dc; + for (; dc < 8; dc++) x *= 16; + + long e2; + if ((c | 32) == 'p') + { + long! e2val = str::to_long((String)chars[index + 1..]); + if (catch e2val) return NumberConversion.MALFORMED_FLOAT!; + e2 = e2val; + } + e2 += 4 * rp - 32; + if (!x) return sign * 0.0; + if (e2 > -$emin) return NumberConversion.FLOAT_OUT_OF_RANGE!; + if (e2 < $emin - 2 * math::DOUBLE_MANT_DIG) return NumberConversion.FLOAT_OUT_OF_RANGE!; + + while (x < 0x80000000) + { + if (y >= 0.5) + { + x += x + 1; + y += y - 1; + } + else + { + x += x; + y += y; + } + e2--; + } + int bits = $bits; + if ($bits > 32 + e2 - $emin) + { + bits = (int)(32 + e2 - $emin); + if (bits < 0) bits = 0; + } + double bias; + if (bits < math::DOUBLE_MANT_DIG) + { + bias = math::copysign(math::scalbn(1, 32 + math::DOUBLE_MANT_DIG - bits - 1), (double)sign); + } + + if (bits < 32 && y && !(x & 1)) + { + x++; + y = 0; + } + y = bias + sign * (double)x + sign * y; + y -= bias; + if (!y) return NumberConversion.FLOAT_OUT_OF_RANGE!; + + return math::scalbn(y, (int)e2); +} + +private macro floatparse(String chars, $Type) +{ + int sign = 1; + $switch ($Type): + $case float: + const int BITS = math::FLOAT_MANT_DIG; + const int EMIN = math::FLOAT_MIN_EXP - BITS; + $case double: + const int BITS = math::DOUBLE_MANT_DIG; + const int EMIN = math::DOUBLE_MIN_EXP - BITS; + $case float128: + $assert(false, "Not yet supported"); + $default: + $assert(false, "Unexpected type"); + $endswitch; + + while (chars.len && chars[0] == ' ') chars = chars[1..]; + if (!chars.len) return NumberConversion.MALFORMED_FLOAT!; + switch (chars[0]) + { + case '-': + sign = -1; + nextcase; + case '+': + chars = chars[1..]; + } + if (chars == "infinity" || chars == "INFINITY") return sign * $Type.inf; + if (chars == "NAN" || chars == "nan") return $Type.nan; + + if (chars.len > 2 && chars[0] == '0' && (chars[1] | 32) == 'x') + { + return ($Type)hexfloat((char[])chars[2..], BITS, EMIN, sign); + } + return ($Type)decfloat((char[])chars, BITS, EMIN, sign); +} + diff --git a/lib/std/core/str.c3 b/lib/std/core/str.c3 index 7d09211c3..6c6632826 100644 --- a/lib/std/core/str.c3 +++ b/lib/std/core/str.c3 @@ -19,6 +19,8 @@ fault NumberConversion NEGATIVE_VALUE, MALFORMED_INTEGER, INTEGER_OVERFLOW, + MALFORMED_FLOAT, + FLOAT_OUT_OF_RANGE, } fn VarString join(String[] s, String joiner) { @@ -123,6 +125,8 @@ private macro to_integer($Type, String string) return value; } +fn float! to_float(String string) => floatparse(string, float); +fn double! to_double(String string) => floatparse(string, double); fn int128! to_int128(String string) => to_integer(int128, string); fn long! to_long(String string) => to_integer(long, string); fn int! to_int(String string) => to_integer(int, string); diff --git a/lib/std/io/io_formatter_private.c3 b/lib/std/io/io_formatter_private.c3 index 70961813c..52fee86c2 100644 --- a/lib/std/io/io_formatter_private.c3 +++ b/lib/std/io/io_formatter_private.c3 @@ -151,17 +151,10 @@ private fn void! Formatter.out_substr(Formatter *this, String str) return this.left_adjust(l); } - -union ConvUnion -{ - ulong u; - double f; -} - private fn void! Formatter.etoa(Formatter* this, FloatType value) { // check for NaN and special values - if (value != value || value < -FloatType.max || value > FloatType.max) + if (!value || value != value || value < -FloatType.max || value > FloatType.max) { return this.ftoa(value); } @@ -178,25 +171,24 @@ private fn void! Formatter.etoa(Formatter* this, FloatType value) // determine the decimal exponent // based on the algorithm by David Gay (https://www.ampl.com/netlib/fp/dtoa.c) - ConvUnion conv; - - conv.f = (double)value; - int exp2 = (int)(conv.u >> 52 & 0x7FF) - 1023; // effectively log2 - conv.u = (conv.u & (1u64 << 52 - 1)) | (1023u64 << 52); // drop the exponent so conv.F is now in [1,2) + ulong u = bitcast((double)value, ulong); + int exp2 = (int)(u >> 52 & 0x7FFu) - 1023; // effectively log2 + u = (u & (1u64 << 52 - 1u)) | (1023u64 << 52u); // drop the exponent so conv.F is now in [1,2) // now approximate log10 from the log2 integer part and an expansion of ln around 1.5 - int expval = (int)(0.1760912590558 + exp2 * 0.301029995663981 + (conv.f - 1.5) * 0.289529654602168); + double f = bitcast(u, double); + int expval = (int)(0.1760912590558 + exp2 * 0.301029995663981 + (f - 1.5) * 0.289529654602168); // now we want to compute 10^expval but we want to be sure it won't overflow exp2 = (int)(expval * 3.321928094887362 + 0.5); double z = expval * 2.302585092994046 - exp2 * 0.6931471805599453; double z2 = z * z; - conv.u = (ulong)(exp2 + 1023) << 52; + f = bitcast((ulong)(exp2 + 1023) << 52, double); // compute exp(z) using continued fractions, see https://en.wikipedia.org/wiki/Exponential_function#Continued_fractions_for_ex - conv.f *= 1 + 2 * z / (2 - z + (z2 / (6 + (z2 / (10 + z2 / 14))))); + f *= 1 + 2 * z / (2 - z + (z2 / (6 + (z2 / (10 + z2 / 14))))); // correct for rounding errors - if (value < conv.f) + if (value < f) { expval--; - conv.f /= 10; + f /= 10; } // the exponent format is "%+03d" and largest value is "307", so set aside 4-5 characters @@ -220,7 +212,6 @@ private fn void! Formatter.etoa(Formatter* this, FloatType value) if (this.prec > 0 && this.flags.precision) this.prec--; } } - // Adjust width uint fwidth = this.width > minwidth ? this.width - minwidth : 0; @@ -228,8 +219,7 @@ private fn void! Formatter.etoa(Formatter* this, FloatType value) if (this.flags.left && minwidth) fwidth = 0; // rescale the float value - if (expval) value /= conv.f; - + if (expval) value /= f; // output the floating part usz start_idx = this.idx; PrintFlags old = this.flags; @@ -303,14 +293,12 @@ private fn void! Formatter.ftoa(Formatter* this, FloatType value) buf[len++] = '0'; this.prec--; } - // Safe due to 1e9 limit. int whole = (int)value; FloatType tmp = (value - whole) * POW10[this.prec]; ulong frac = (ulong)tmp; diff = tmp - frac; - - switch (true) + switch { case diff > 0.5: ++frac; diff --git a/lib/std/math/math.c3 b/lib/std/math/math.c3 index 5996a29e5..1bda8f50b 100644 --- a/lib/std/math/math.c3 +++ b/lib/std/math/math.c3 @@ -688,6 +688,7 @@ macro is_nan(x) $endswitch; } +macro double scalbn(double x, int n) => _scalbn(x, n); extern fn double _atan(double x) @extern("atan"); extern fn float _atanf(float x) @extern("atanf"); @@ -696,4 +697,5 @@ extern fn float _atan2f(float, float) @extern("atan2f"); extern fn void _sincos(double, double*) @extern("sincos"); extern fn void _sincosf(float, float*) @extern("sincosf"); extern fn double _tan(double x) @extern("tan"); -extern fn float _tanf(float x) @extern("tanf"); \ No newline at end of file +extern fn float _tanf(float x) @extern("tanf"); +extern fn double _scalbn(double x, int n) @extern("scalbn"); \ No newline at end of file diff --git a/src/compiler/llvm_codegen_expr.c b/src/compiler/llvm_codegen_expr.c index cdcb1be78..dd78d6689 100644 --- a/src/compiler/llvm_codegen_expr.c +++ b/src/compiler/llvm_codegen_expr.c @@ -2330,7 +2330,7 @@ static void llvm_emit_unary_expr(GenContext *c, BEValue *value, Expr *expr) Type *vec_type = type_vector_type(type); if (type_is_float(vec_type)) { - llvm_value = LLVMBuildFCmp(c->builder, LLVMRealUNE, value->value, llvm_get_zero(c, type), "not"); + llvm_value = LLVMBuildFCmp(c->builder, LLVMRealUEQ, value->value, llvm_get_zero(c, type), "not"); } else { @@ -2345,7 +2345,7 @@ static void llvm_emit_unary_expr(GenContext *c, BEValue *value, Expr *expr) { case ALL_FLOATS: llvm_value_rvalue(c, value); - llvm_value = LLVMBuildFCmp(c->builder, LLVMRealUNE, value->value, llvm_get_zero(c, type), "not"); + llvm_value = LLVMBuildFCmp(c->builder, LLVMRealUEQ, value->value, llvm_get_zero(c, type), "not"); break; case TYPE_BOOL: llvm_value_rvalue(c, value); diff --git a/src/version.h b/src/version.h index f711c57e8..446cb6239 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define COMPILER_VERSION "0.4.54" \ No newline at end of file +#define COMPILER_VERSION "0.4.55" \ No newline at end of file diff --git a/test/test_suite/vector/vector_bit.c3t b/test/test_suite/vector/vector_bit.c3t index 6727d6eda..cfc851607 100644 --- a/test/test_suite/vector/vector_bit.c3t +++ b/test/test_suite/vector/vector_bit.c3t @@ -44,7 +44,7 @@ entry: %w = alloca <4 x i32>, align 16 store <4 x float> , ptr %y, align 16 %0 = load <4 x float>, ptr %y, align 16 - %not = fcmp une <4 x float> %0, zeroinitializer + %not = fcmp ueq <4 x float> %0, zeroinitializer %1 = sext <4 x i1> %not to <4 x i32> store <4 x i32> %1, ptr %w, align 16 %2 = load <4 x i32>, ptr %w, align 16 diff --git a/test/unit/stdlib/string_to_float.c3 b/test/unit/stdlib/string_to_float.c3 new file mode 100644 index 000000000..71486fe54 --- /dev/null +++ b/test/unit/stdlib/string_to_float.c3 @@ -0,0 +1,29 @@ +module string_to_float_tests; + +fn void! test_float() @test +{ + assert(str::to_float("1.2")? == 1.2f); + assert(str::to_float("10")? == 10f); + assert(str::to_float(".7647834")? == 0.7647834f); + assert(str::to_float("0.213232")? == 0.213232f); + assert(str::to_float("000001.487348")? == 000001.487348f); + assert(str::to_float("3.54500000")? == 3.54500000f); + assert(str::to_float("4.0")? == 4.0f); + assert(str::to_float("-23.545")? == -23.545f); + assert(str::to_float("1.5555555555555")? == 1.5555555555555f); + assert(str::to_float("1.5555555555556666")? == 1.5555555555556666f); +} + +fn void! test_double() @test +{ + assert(str::to_double("1.2")? == 1.2); + assert(str::to_double("10")? == 10); + assert(str::to_double(".7647834")? == 0.7647834); + assert(str::to_double("0.213232")? == 0.213232); + assert(str::to_double("000001.487348")? == 000001.487348); + assert(str::to_double("3.54500000")? == 3.54500000); + assert(str::to_double("4.0")? == 4.0); + assert(str::to_double("-23.545")? == -23.545); + assert(str::to_double("1.5555555555555")? == 1.5555555555555); + assert(str::to_double("1.5555555555556666")? == 1.5555555555556666); +} \ No newline at end of file