Fix of bug in !floatval codegen. Added float parsing.

This commit is contained in:
Christoffer Lerno
2023-02-10 22:17:55 +01:00
parent b3f15a867c
commit 74d868d113
8 changed files with 540 additions and 29 deletions

488
lib/std/core/floatparse.c3 Normal file
View File

@@ -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);
}

View File

@@ -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);

View File

@@ -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;

View File

@@ -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");
extern fn float _tanf(float x) @extern("tanf");
extern fn double _scalbn(double x, int n) @extern("scalbn");

View File

@@ -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);

View File

@@ -1 +1 @@
#define COMPILER_VERSION "0.4.54"
#define COMPILER_VERSION "0.4.55"

View File

@@ -44,7 +44,7 @@ entry:
%w = alloca <4 x i32>, align 16
store <4 x float> <float 0.000000e+00, float 2.000000e+00, float 3.000000e+00, float 0.000000e+00>, 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

View File

@@ -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);
}