diff --git a/docs/system.md b/docs/system.md index ed5dddb8..cafd51bd 100644 --- a/docs/system.md +++ b/docs/system.md @@ -22,6 +22,7 @@ See [the BQN specification](https://mlochbaum.github.io/BQN/spec/system.html) fo | `•Out` | | | `•Show` | | | `•Repr` | | +| `•ParseFloat` | Should exactly round floats with up to 17 significant digits, but won't necessarily round correctly with more | | `•Fmt` | | | `•term` | Fields: `Flush`, `RawMode`, `CharB`, `CharN`; has extensions | | `•SH` | See [•SH](#sh) | diff --git a/src/builtins.h b/src/builtins.h index c1691095..a4fd7d7a 100644 --- a/src/builtins.h +++ b/src/builtins.h @@ -8,9 +8,10 @@ /* sfns.c*/A(shifta,"«") A(take,"↑") A(drop,"↓") A(group,"⊔") A(reverse,"⌽") A(transp,"⍉") \ /* sort.c*/A(gradeUp,"⍋") A(gradeDown,"⍒") \ /* everything before the definition of •Type is defined to be pure, and everything after is not */ \ -/* sysfn.c*/M(type,"•Type") M(decp,"•Decompose") M(primInd,"•PrimInd") M(glyph,"•Glyph") A(fill,"•FillFn") M(sys,"•getsys") A(grLen,"•GroupLen") D(grOrd,"•GroupOrd") \ -/* sysfn.c*/M(repr,"•Repr") M(fmt,"•Fmt") A(asrt,"!") A(casrt,"!") M(out,"•Out") M(show,"•Show") A(bqn,"•BQN") A(sh,"•SH") M(fromUtf8,"•FromUTF8") M(toUtf8,"•ToUTF8") M(currentError,"•CurrentError") \ -/* sysfn.c*/D(cmp,"•Cmp") A(hash,"•Hash") M(unixTime,"•UnixTime") M(monoTime,"•MonoTime") M(delay,"•Delay") M(makeRand,"•MakeRand") M(reBQN,"•ReBQN") M(exit,"•Exit") M(getLine,"•GetLine") \ +/* sysfn.c*/M(type,"•Type") M(decp,"•Decompose") M(primInd,"•PrimInd") M(glyph,"•Glyph") A(fill,"•FillFn") M(sys,"•getsys") A(grLen,"•GroupLen") \ +/* sysfn.c*/D(grOrd,"•GroupOrd") M(repr,"•Repr") M(parseFloat,"•ParseFloat") M(fmt,"•Fmt") A(asrt,"!") A(casrt,"!") M(out,"•Out") M(show,"•Show") A(bqn,"•BQN") \ +/* sysfn.c*/A(sh,"•SH") M(fromUtf8,"•FromUTF8") M(toUtf8,"•ToUTF8") M(currentError,"•CurrentError") D(cmp,"•Cmp") A(hash,"•Hash") M(unixTime,"•UnixTime")\ +/* sysfn.c*/M(monoTime,"•MonoTime") M(delay,"•Delay") M(makeRand,"•MakeRand") M(reBQN,"•ReBQN") M(exit,"•Exit") M(getLine,"•GetLine") \ /* sysfn.c*/M(fName,"•file.Name") M(fParent,"•file.Parent") \ /* sysfn.c*/M(tRawMode,"•term.RawMode") M(tFlush,"•term.Flush") M(tCharB,"•term.CharB") M(tCharN,"•term.CharN") M(tOutRaw,"•term.OutRaw") M(tErrRaw,"•term.ErrRaw") \ /* inverse.c*/M(setInvReg,"(SetInvReg)") M(setInvSwap,"(SetInvSwap)") M(nativeInvReg,"(NativeInvReg)") M(nativeInvSwap,"(NativeInvSwap)") \ diff --git a/src/builtins/sysfn.c b/src/builtins/sysfn.c index c6fc148a..114fdb31 100644 --- a/src/builtins/sysfn.c +++ b/src/builtins/sysfn.c @@ -83,6 +83,7 @@ B glyph_c1(B t, B x) { #if !NO_RYU B ryu_d2s(double f); +bool ryu_s2d_n(u8* buffer, int len, f64* result); #endif B repr_c1(B t, B x) { @@ -110,6 +111,27 @@ B fmt_c1(B t, B x) { #endif } +#if NO_RYU +B parseFloat_c1(B t, B x) { thrM("•ParseFloat: Not supported with Ryu disabled"); } +#else +B parseFloat_c1(B t, B x) { + if (isAtm(x)) thrM("•ParseFloat: Expected a character list argument"); + if (TI(x,elType)!=el_c8) { + x = chr_squeeze(x); + if (TI(x,elType)!=el_c8) thrM("•ParseFloat: Expected a character list argument"); + } + usz ia = IA(x); + if (RNK(x)!=1) thrM("•ParseFloat: Input must have rank 1"); + if (ia==0) thrM("•ParseFloat: Input was empty"); + if (ia >= (1<<20)) thrM("•ParseFloat: Input too long"); // otherwise things like + u8* data = c8any_ptr(x); + f64 res; + if (!ryu_s2d_n(data, ia, &res)) thrM("•ParseFloat: Malformed input"); + decG(x); + return m_f64(res); +} +#endif + B fill_c1(B t, B x) { B r = getFillE(x); dec(x); @@ -1416,6 +1438,7 @@ static Body* file_nsGen; F("delay", U"•Delay", bi_delay) \ F("hash", U"•Hash", bi_hash) \ F("repr", U"•Repr", bi_repr) \ + F("parsefloat", U"•ParseFloat", bi_parseFloat) \ F("fmt", U"•Fmt", bi_fmt) \ F("glyph", U"•Glyph", bi_glyph) \ F("makerand", U"•MakeRand", bi_makeRand) \ diff --git a/src/utils/ryu.c b/src/utils/ryu.c index 8b913b93..0911ef76 100644 --- a/src/utils/ryu.c +++ b/src/utils/ryu.c @@ -23,6 +23,7 @@ void ryu_init() { } #else +// original file: d2s.h #include "../core.h" #include "ryu/ryu_common.h" @@ -270,6 +271,7 @@ static inline floating_decimal_64 d2d(const uint64_t mantissa, const uint32_t ex + #define W2D(RES, NUM) memcpy((RES), DIGIT_TABLE + (NUM)*2, 2) #define W4D(RES, NUM) ({ char* r4_ = (RES); AUTO n4_ = (NUM); W2D(r4_, n4_/100 ); W2D(r4_+2, n4_%100 ); }) @@ -298,7 +300,6 @@ static inline bool d2d_small_int(const uint64_t mantissa, const uint32_t exponen return true; } - static inline B to_chars(const floating_decimal_64 v, const bool sign, const bool forcePositional) { char buf[25]; char* dec = buf+25; @@ -380,7 +381,6 @@ static inline B to_chars(const floating_decimal_64 v, const bool sign, const boo return r; } - static B fmt_nan, fmt_inf[2], fmt_zero[2]; B ryu_d2s(double f) { const uint64_t bits = double_to_bits(f); // decode the floating-point number, and unify normalized and subnormal cases. @@ -429,6 +429,212 @@ B ryu_d2s(double f) { return to_chars(v, sign, forcePositional); } + + + + + +// original file: s2d.h +#if defined(_MSC_VER) + #include + static inline uint32_t floor_log2(const uint64_t value) { + long index; + return _BitScanReverse64(&index, value) ? index : 64; + } +#else + static inline uint32_t floor_log2(const uint64_t value) { return 63 - CLZ(value); } +#endif + +// The max function is already defined on Windows. +static inline int32_t max32(int32_t a, int32_t b) { return a < b ? b : a; } + +static inline double int64Bits2Double(uint64_t bits) { + double f; + memcpy(&f, &bits, sizeof(double)); + return f; +} + +bool ryu_s2d_n(u8* buffer, int len, f64* result) { + assert(len>0 && len<(1<<20)); // max length so that '0.0000[a billion zeroes]0000e1000000000' doesn't have to be handled + int m10digits = 0; + int dotIndex = len; + int eIndex = len; + uint64_t m10 = 0; + int32_t e10 = 0; + int32_t offset = 0; + bool signedM = false; + bool signedE = false; + + int i = 0; + if (buffer[i] == '-') { + signedM = true; + i++; + } + int i1 = i; + for (; i < len; i++) { + char c = buffer[i]; + if (c == '.') { + if (dotIndex!=len) return false; + dotIndex = i; + continue; + } + if ((c<'0') || (c>'9')) break; + if (m10digits < 18) { + m10 = 10*m10 + (c-'0'); + if (m10!=0) m10digits++; + } else offset++; + } + if (i-i1<=1) { + if (i==i1) return false; // "-", "e2", "-e2" + if (i1==dotIndex) return false; // ".", ".e2" + } + if (i= len) return false; // "123e", "123e+" + for (; i < len; i++) { + char c = buffer[i]; + if ((c<'0') || (c>'9')) return false; + + // 1<<25 >> 2⋆20, 10×2⋆20 < 2⋆32; some leeway because why not, and also to make sure we're really out of the range where the exponent could in any way compensate a long input + // else, just leave e10 as-is; continue loop to finish checking input validity, but an infinity or ±0 + // from e.g. ±"0e100000000000000" or "123e-100000000000000" will get output regardless of the precise e10 value + if (e10 < (1<<25)) { + e10 = 10*e10 + (c-'0'); + } + } + } + if (i < len) return false; + if (signedE) e10 = -e10; + + e10+= offset; + e10-= dotIndex= 310) { + // Number is larger than 1e+309, which should be rounded to +/-Infinity. + uint64_t ieee = (((uint64_t) signedM) << (DOUBLE_EXPONENT_BITS + DOUBLE_MANTISSA_BITS)) | (0x7ffull << DOUBLE_MANTISSA_BITS); + *result = int64Bits2Double(ieee); + return true; + } + + // Convert to binary float m2 * 2^e2, while retaining information about whether the conversion + // was exact (trailingZeros). + int32_t e2; + uint64_t m2; + bool trailingZeros; + if (e10 >= 0) { + // The length of m * 10^e in bits is: + // log2(m10 * 10^e10) = log2(m10) + e10 log2(10) = log2(m10) + e10 + e10 * log2(5) + // + // We want to compute the DOUBLE_MANTISSA_BITS + 1 top-most bits (+1 for the implicit leading + // one in IEEE format). We therefore choose a binary output exponent of + // log2(m10 * 10^e10) - (DOUBLE_MANTISSA_BITS + 1). + // + // We use floor(log2(5^e10)) so that we get at least this many bits; better to + // have an additional bit than to not have enough bits. + e2 = floor_log2(m10) + e10 + log2pow5(e10) - (DOUBLE_MANTISSA_BITS + 1); + + // We now compute [m10 * 10^e10 / 2^e2] = [m10 * 5^e10 / 2^(e2-e10)]. + // To that end, we use the DOUBLE_POW5_SPLIT table. + int j = e2 - e10 - ceil_log2pow5(e10) + DOUBLE_POW5_BITCOUNT; + assert(j >= 0); +#if defined(RYU_OPTIMIZE_SIZE) + uint64_t pow5[2]; + double_computePow5(e10, pow5); + m2 = mulShift64(m10, pow5, j); +#else + assert(e10 < DOUBLE_POW5_TABLE_SIZE); + m2 = mulShift64(m10, DOUBLE_POW5_SPLIT[e10], j); +#endif + // We also compute if the result is exact, i.e., + // [m10 * 10^e10 / 2^e2] == m10 * 10^e10 / 2^e2. + // This can only be the case if 2^e2 divides m10 * 10^e10, which in turn requires that the + // largest power of 2 that divides m10 + e10 is greater than e2. If e2 is less than e10, then + // the result must be exact. Otherwise we use the existing multipleOfPowerOf2 function. + trailingZeros = e2 < e10 || (e2 - e10 < 64 && multipleOfPowerOf2(m10, e2 - e10)); + } else { + e2 = floor_log2(m10) + e10 - ceil_log2pow5(-e10) - (DOUBLE_MANTISSA_BITS + 1); + int j = e2 - e10 + ceil_log2pow5(-e10) - 1 + DOUBLE_POW5_INV_BITCOUNT; +#if defined(RYU_OPTIMIZE_SIZE) + uint64_t pow5[2]; + double_computeInvPow5(-e10, pow5); + m2 = mulShift64(m10, pow5, j); +#else + assert(-e10 < DOUBLE_POW5_INV_TABLE_SIZE); + m2 = mulShift64(m10, DOUBLE_POW5_INV_SPLIT[-e10], j); +#endif + trailingZeros = multipleOfPowerOf5(m10, -e10); + } + +#ifdef RYU_DEBUG + printf("m2 * 2^e2 = %" PRIu64 " * 2^%d\n", m2, e2); +#endif + + // Compute the final IEEE exponent. + uint32_t ieee_e2 = (uint32_t) max32(0, e2 + DOUBLE_BIAS + floor_log2(m2)); + + if (ieee_e2 > 0x7fe) { + // Final IEEE exponent is larger than the maximum representable; return +/-Infinity. + uint64_t ieee = (((uint64_t) signedM) << (DOUBLE_EXPONENT_BITS + DOUBLE_MANTISSA_BITS)) | (0x7ffull << DOUBLE_MANTISSA_BITS); + *result = int64Bits2Double(ieee); + return true; + } + + // We need to figure out how much we need to shift m2. The tricky part is that we need to take + // the final IEEE exponent into account, so we need to reverse the bias and also special-case + // the value 0. + int32_t shift = (ieee_e2 == 0 ? 1 : ieee_e2) - e2 - DOUBLE_BIAS - DOUBLE_MANTISSA_BITS; + assert(shift >= 0); +#ifdef RYU_DEBUG + printf("ieee_e2 = %d\n", ieee_e2); + printf("shift = %d\n", shift); +#endif + + // We need to round up if the exact value is more than 0.5 above the value we computed. That's + // equivalent to checking if the last removed bit was 1 and either the value was not just + // trailing zeros or the result would otherwise be odd. + // + // We need to update trailingZeros given that we have the exact output exponent ieee_e2 now. + trailingZeros &= (m2 & ((1ull << (shift - 1)) - 1)) == 0; + uint64_t lastRemovedBit = (m2 >> (shift - 1)) & 1; + bool roundUp = (lastRemovedBit != 0) && (!trailingZeros || (((m2 >> shift) & 1) != 0)); + +#ifdef RYU_DEBUG + printf("roundUp = %d\n", roundUp); + printf("ieee_m2 = %" PRIu64 "\n", (m2 >> shift) + roundUp); +#endif + uint64_t ieee_m2 = (m2 >> shift) + roundUp; + assert(ieee_m2 <= (1ull << (DOUBLE_MANTISSA_BITS + 1))); + ieee_m2 &= (1ull << DOUBLE_MANTISSA_BITS) - 1; + if (ieee_m2 == 0 && roundUp) { + // Due to how the IEEE represents +/-Infinity, we don't need to check for overflow here. + ieee_e2++; + } + + uint64_t ieee = (((((uint64_t) signedM) << DOUBLE_EXPONENT_BITS) | (uint64_t)ieee_e2) << DOUBLE_MANTISSA_BITS) | ieee_m2; + *result = int64Bits2Double(ieee); + return true; +} + void ryu_init() { gc_add(fmt_nan = m_c8vec("NaN", 3)); gc_add(fmt_zero[0] = m_c8vec("0", 1)); diff --git a/src/utils/ryu/d2s_small_table.h b/src/utils/ryu/d2s_small_table.h index f3c260c0..a5fbcc05 100644 --- a/src/utils/ryu/d2s_small_table.h +++ b/src/utils/ryu/d2s_small_table.h @@ -17,9 +17,6 @@ #ifndef RYU_D2S_SMALL_TABLE_H #define RYU_D2S_SMALL_TABLE_H -// Defines HAS_UINT128 and uint128_t if applicable. -#include "ryu_common.h" - // These tables are generated by PrintDoubleLookupTable. #define DOUBLE_POW5_INV_BITCOUNT 125 #define DOUBLE_POW5_BITCOUNT 125 diff --git a/src/utils/ryu/ryu_common.h b/src/utils/ryu/ryu_common.h index 46be17ee..81842515 100644 --- a/src/utils/ryu/ryu_common.h +++ b/src/utils/ryu/ryu_common.h @@ -36,6 +36,16 @@ static const char DIGIT_TABLE[200] = { #define RYU_32_BIT_PLATFORM #endif +// Returns e == 0 ? 1 : [log_2(5^e)]; requires 0 <= e <= 3528. +static inline int32_t log2pow5(const int32_t e) { + // This approximation works up to the point that the multiplication overflows at e = 3529. + // If the multiplication were done in 64 bits, it would fail at 5^4004 which is just greater + // than 2^9297. + assert(e >= 0); + assert(e <= 3528); + return (int32_t) ((((uint32_t) e) * 1217359) >> 19); +} + // Returns e == 0 ? 1 : ceil(log_2(5^e)); requires 0 <= e <= 3528. static inline int32_t pow5bits(const int32_t e) { // This approximation works up to the point that the multiplication overflows at e = 3529. @@ -46,6 +56,11 @@ static inline int32_t pow5bits(const int32_t e) { return (int32_t) (((((uint32_t) e) * 1217359) >> 19) + 1); } +// Returns e == 0 ? 1 : ceil(log_2(5^e)); requires 0 <= e <= 3528. +static inline int32_t ceil_log2pow5(const int32_t e) { + return log2pow5(e) + 1; +} + // Returns floor(log_10(2^e)); requires 0 <= e <= 1650. static inline uint32_t log10Pow2(const int32_t e) { // The first value this approximation fails for is 2^1651 which is just greater than 10^297. diff --git a/test/ffi/test.expected b/test/ffi/test.expected index a7535790..8207b92c 100644 --- a/test/ffi/test.expected +++ b/test/ffi/test.expected @@ -39,7 +39,7 @@ args: -123 -12323 -212312312 250 50000 3123456789 3.141592741012573242 0.3333333 args: -123 -12323 -212312312 250 50000 3123456789 3.141592741012573242 0.333333333333333315 @ args: -123 -12323 -212312312 250 50000 3123456789 3.141592741012573242 0.333333333333333315 -⟨ ⟨ ¯123 ⟩ ⟨ ¯12323 ⟩ ⟨ ¯212312312 ⟩ ⟨ 250 ⟩ ⟨ 50000 ⟩ ⟨ 3123456789 ⟩ ⟨ 3.141592741012573 ⟩ ⟨ 0.3333333333333333 ⟩ ⟩ +⟨ ⟨ ¯123 ⟩ ⟨ ¯12323 ⟩ ⟨ ¯212312312 ⟩ ⟨ 250 ⟩ ⟨ 50000 ⟩ ⟨ 3123456789 ⟩ ⟨ 3.1415927410125732 ⟩ ⟨ 0.3333333333333333 ⟩ ⟩ 726f776f6c6c6568 6161616161616161 3837363534333231 7478657474786574 12345678 ff7fdfefefdf7bb4 ff7fdfefefdf7bb4 fefffdfff7ffbffb bffff7fffdfffeff @@ -74,7 +74,7 @@ ff7fdfefefdf7bb4 ff7fdfefefdf7bb4 fefffdfff7ffbffb bffff7fffdfffeff # u64 tests 51539608786 -4503599628419072 +4.503599628419072e15 ⟨ ¯2045800064 28744 ⟩ ⟨ 0 0 0 1 0 1 1 0 1 0 1 0 0 1 1 0 0 0 1 1 0 1 1 0 0 0 1 1 0 1 1 0 1 1 1 1 0 1 1 0 1 1 1 0 1 1 1 0 1 1 1 1 0 1 1 0 0 1 0 0 1 1 1 0 ⟩