Add support for f128 values in std.fmt.parseFloat.

[danielchasehooper]

Add f128 support to std.fmt.parseFloat by transitioning it to sqlite's float parser, adapted to support 128 floats. The old implementation worked in f64 internally, and then upcast to f128 at the last moment.

Additionally, this fixes a typo in parse_f128.c - my initial approach for this PR was to port that parser to zig as recommended in #11169 (comment), but found it had issues parsing large f64 values like 144115188075855870, even after fixing the typo. I suspect there may be other issues with parse_f128.c, but due to its style, it's hard to know what else might be wrong. I also explored adapting the existing parser to truly support f128 (and not just upcast from f64 to f128), but in early tests it exhibited incorrect behavior for parsing large f64 values, and it relied on some custom 96bit-integer code that made it difficult to adapt to f128. Sqlite's f64 parser behaved correctly out of the box and is much easier to read, so was easy to adapt to support f128.

I'm not sure what the implications are for relying on f128 math in parseFloat now - vs working in f64 and then casting to f16/f32/f64/f128 at the last moment like the old implementation did. Is f128 math emulated in software for targets that don't natively support it? Should we make parseFloat work in f64 when f128 isn't the requested return type?

Fixes #11169

[danielchasehooper]

Not sure what's going on with the test failure:
992/1116 test.comptime union field access... FAIL (TestUnexpectedResult)

zig test test/behavior/union.zig works locally for me (arm64 macOS). Any insight would be appreciated.

[matu3ba]

but found it had issues parsing large f64 values like 144115188075855870, even after fixing the typo. I suspect there may be other issues with parse_f128.c, but due to its style, it's hard to know what else might be wrong.

1. Do you have the test examples anywhere, which break/are wrong with parse_f128.c?
2. How did you ensure that the deviation from f64 was not due to better accuracy?
3. Did you cross-compare with other implementations like these https://nigeltao.github.io/blog/2020/eisel-lemire.html ?
4. Any information on test data would help on this.

Also: I found this library, which is MIT licensed: http://www.kerbaya.com/ieee754lib/index.html
Also: There appears to be only proprietary licensed test programs for 128bit floats. :(
"IEEE 754 has no official tests for floating-point but there are well-known third party tools to check such as John Hauser's TestFloat." https://libre-soc.org/resources/
see "IeeeCC754++" and its predecessor IEEECC754 https://www.uantwerpen.be/en/research-groups/cma/research/software/
"With the current release, the following binary formats can be tested: 16-bit half-precision, 32-bit single-precision, 64-bit double-precision, 80-bit double-extended-precision, and/or 128-bit quadruple-precision. TestFloat cannot test decimal floating-point."
http://www.jhauser.us/arithmetic/TestFloat.html
Do you think testfloat can be hacked up to generate useful test cases for parsing?

casting to f16/f32/f64/f128 at the last moment like the old implementation did

Accuracy is lower with less bits, so one wants to use all of them and not upcast things.

[danielchasehooper]

@matu3ba You probably know this, but for clarity I'll mention it: parse_f128.c and the previous implementation of parseFloat that this PR replaces are two separate parsers. This PR is not directly about parse_f128.c - it's about making parseFloat pass tests that were previously failing. Those tests now pass.

1. Do you have the test examples anywhere, which break/are wrong with parse_f128.c?

Mentioning parse_f128.c may have been a distraction. I only mentioned it because at first my approach for this PR was to port it to zig in order to make parseFloat support f128. Once I noticed that parse_f128.c had typos and my port had strange behavior for some numbers, I abandoned that path because the port was hard enough even when I had faith in the source material. It's possible that parse_f128.c is perfect (besides that one typo this PR corrects) and the issue was in my zig port of it. But the code is so unreadable I abandoned it. Regardless, any potential issues with it are separate from this PR. I'd be happy to write a new GitHub issue about testing parse_f128.c if you think it's worth it. Probably best to forget I mentioned parse_f128.c at all.

2. How did you ensure that the deviation from f64 was not due to better accuracy?

144115188075855870 can be represented exactly by both f64 and f128. The previous implementation of parseFloat didn't get it right when parsing to f64 OR f128. Even if it was able to get that specific number correct, the previous parseFloat used f64s internally, and so you always got f64's range and precision, even if asking for a f128.

3. Did you cross-compare with other implementations like these https://nigeltao.github.io/blog/2020/eisel-lemire.html ?

I compared the behavior of this PR against c's atof and found that they matched for 64 bit floats. atof and that implementation you linked to do not support f128, so I wasn't able to verify any f128 numbers. One of the tests for f128 I manually figured out the 128 bits the number should be parsed as - the test has that written as hex.

Do you think testfloat can be hacked up to generate useful test cases for parsing?

I'll look into this... more rigorous testing beyond what we had for the previous implementation of parseFloat may be a separate task.

Accuracy is lower with less bits, so one wants to use all of them and not upcast things.

I agree, and the parser working in a higher bit depth is the point of this PR. I was suggesting that maybe it would be worthwhile to make parseFloat be written generically so that its calculations happen in f64 when an f64 number is requested. This would have benefits for targets which f128 is software emulated/not supported. If a f128 type was requested, then yes, it would work in the higher bit depth for accuracy.

[topolarity]

Really clean and readable implementation - Very nice work 👍

Can you clarify the rounding behavior, and any support for subnormals? Based on the code, those probably depend on the configuration of the FPU/softfloat library at the time of parsing?

[topolarity]

This looks like it's off by 1 bit to me - should be 0x7ff8136c69ce8adff4397b050cae44c6

Based on a rudimentary calculation in Python, intended to get the error in ulp:

>>> exp = 0x7ff8 - 0x3fff
>>> (10**4930 - 0x1136c69ce8adff4397b050cae44c7 * 2**(exp-112)) / 2**(exp-112)
-1.2945749668158448
>>> (10**4930 - 0x1136c69ce8adff4397b050cae44c6 * 2**(exp-112)) / 2**(exp-112)
-0.2945749668158449

I might have made a mistake in the calculation, of course

[andrewrk]

Not sure what's going on with the test failure: 992/1116 test.comptime union field access... FAIL (TestUnexpectedResult)

zig test test/behavior/union.zig works locally for me (arm64 macOS). Any insight would be appreciated.

Disable the non-LLVM backends in the newly passing test case:

    if (builtin.zig_backend == .stage2_wasm) return error.SkipZigTest; // TODO
    if (builtin.zig_backend == .stage2_c) return error.SkipZigTest; // TODO
    if (builtin.zig_backend == .stage2_x86_64) return error.SkipZigTest; // TODO
    if (builtin.zig_backend == .stage2_arm) return error.SkipZigTest; // TODO
    if (builtin.zig_backend == .stage2_aarch64) return error.SkipZigTest; // TODO

These other backends do not handle f128 yet.

[danielchasehooper]

@andrewrk

Disable the non-LLVM backends in the newly passing test case:

I'm confused, The failing test cases on the ci server seem unrelated to the behavior/floatop.zig "float literal at compile time not lossy" test - How does disabling this test make other tests start succeeding in the ci (test.comptime union field access)?

Additionally, parseFloat now requires f128 internally - is that going to be an issue for the backends that don't support f128? My musings from my previous comment are related, and I'm curious to get your take:

maybe it would be worthwhile to make parseFloat be written generically so that its calculations happen in f64 when an f64 number is requested. This would have benefits for targets which f128 is software emulated/not supported. If a f128 type was requested, then yes, it would work in the higher bit depth for accuracy.

[andrewrk]

I'm confused, The failing test cases on the ci server seem unrelated to the behavior/floatop.zig "float literal at compile time not lossy" test - How does disabling this test make other tests start succeeding in the ci (test.comptime union field access)?

I didn't look closely when I made this suggestion. You will hit the problem that I mentioned next after you solve this one.

The command that is failing, on both x86_64-macos and x86_64-linux, is:

+ stage2/bin/zig test test/behavior.zig -I test -fLLVM
553/1116 test.anytype params... FAIL (TestUnexpectedResult)
838/1116 test.fully anonymous struct... FAIL (TestUnexpectedResult)
839/1116 test.fully anonymous list literal... FAIL (TestUnexpectedResult)
992/1116 test.comptime union field access... FAIL (TestUnexpectedResult)
1036/1116 test.anonymous union literal syntax... FAIL (TestUnexpectedResult)

All of these tests have float literals in them.

Regarding reliance on f128 - it's OK to rely on f128. On targets where such type is not supported by the hardware, zig provides compiler-rt functions for softfloat implementations. However, it would likely be more efficient to limit use to f64 when parsing an f64, and limit use to f32 when parsing an f32, and I would like to see that improvement made eventually. That improvement could come as a follow-up to this PR, or it could be part of this one.

[topolarity]

However, it would likely be more efficient to limit use to f64 when parsing an f64, and limit use to f32 when parsing an f32, and I would like to see that improvement made eventually. That improvement could come as a follow-up to this PR, or it could be part of this one.

Another note for future optimization: A lot of the fastest float parsers seem to favor a "fast path" where the math can be done easily in smaller types, and then fall back to slower options only for extremely long/precise literals. This tries to give the "best of both worlds": exact precision/rounding in the difficult case, and fast performance in the common case

The fast path can usually be chosen based on the number of digits in the literal, and the decimal exponent. fast_double_parser takes this approach, for example, and I think the stage1 float parser has a few similar tricks

[andrewrk]

I would like to repeat my suggestion to port musl's float parsing code (which has already been ported to self-contained .c code in parse_f128.c). If you want to port a different implementation instead (such as the one in this PR), can you provide some benchmark results to compare perf against it?

[danielchasehooper]

@andrewrk

can you provide some benchmark results to compare perf against it?

Sure! Since I don't have a zig port of musl's parser, I put both the sqlite and musl parsers in a .c file and had them parse argv[1] a million times. sqlite's implementation is 3 to 150 times faster than musl's, depending on the number being parsed. I believe sqlite is faster because it uses the FPU to construct the float (including denormals), rather than doing a bunch of bit twiddling. I included the source for my test below for anyone that would like to check my work.

If you're satisfied with this, I'll get to making the tests pass and look into making it work in f64/f32 when f128 isn't requested.

@topolarity I noticed in this test that this unaltered sqlite implementation differs strtod/musl by 4 ULP for numbers with large exponents on arm64. If Andrew is ok with pursuing this sqlite implementation, I can look into the source of the difference.

x64 in rosetta mode

~/Developer/float_parse_perf $ clang -O3 main.c
~/Developer/float_parse_perf $ file a.out
a.out: Mach-O 64-bit executable x86_64
~/Developer/float_parse_perf $ ./a.out 3.141592654
sql 1.0271s
musl 2.59674s
~/Developer/float_parse_perf $ ./a.out 3141592654
sql 0.020399s
musl 1.78777s
~/Developer/float_parse_perf $ ./a.out +1.53e238
sql 1.88123s
musl 14.0651s
~/Developer/float_parse_perf $ ./a.out -1.53e-238
sql 1.99978s
musl 5.0289s

arm64

~/Developer/float_parse_perf $ clang -O3 main.c
~/Developer/float_parse_perf $ file a.out
a.out: Mach-O 64-bit executable arm64
~/Developer/float_parse_perf $ ./a.out 3.141592654
parsing '3.141592654' (strtod f64 bits 0x400921fb54524550)
sqlite 0.020455s
musl 0.146159s
~/Developer/float_parse_perf $ ./a.out 3141592654
parsing '3141592654' (strtod f64 bits 0x41e7681cc9c00000)
sqlite 0.016496s
musl 0.121425s
~/Developer/float_parse_perf $ ./a.out 31e184
parsing '31e184' (strtod f64 bits 0x66723d3840ca92e3)
sqlite didn't match: 0x66723d3840ca92e5
sqlite 0.034005s
musl 5.10116s
~/Developer/float_parse_perf $ ./a.out 31e-184
parsing '31e-184' (strtod f64 bits 0x1a0a582d82d507cb)
sqlite didn't match: 0x1a0a582d82d507c8
sqlite 0.025871s
musl 1.95019s

main.c source

#include <stdint.h>
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <limits.h>
#include <errno.h>
#include <ctype.h>
#include <stdbool.h>
#include <string.h>
#include <assert.h>
#include <time.h>
#include <stdlib.h>

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024

#define LD_B1B_DIG 2
#define LD_B1B_MAX 9007199, 254740991
#define KMAX 128

#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384

#define LD_B1B_DIG 3
#define LD_B1B_MAX 18, 446744073, 709551615
#define KMAX 2048

#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384

#define LD_B1B_DIG 4
#define LD_B1B_MAX 10384593, 717069655, 257060992, 658440191
#define KMAX 2048

#else
#error Unsupported long double representation
#endif
#define MASK (KMAX-1)

typedef struct {
    const char *str;
    uint64_t off;
    uint64_t len;
} StringFile;

inline void shunget(StringFile* f) {
    f->off--;
    if (f->off < 0) f->off = 0;
}

inline int shgetc(StringFile* f) {
    if (f->off >= f->len) return 0;
    int res = f->str[f->off];
    f->off++;
    return res;
}

static long long scanexp(StringFile *f)
{
    int c;
    int x;
    long long y;
    int neg = 0;
    
    c = shgetc(f);
    if (c=='+' || c=='-') {
        neg = (c=='-');
        c = shgetc(f);
        if (c-'0'>=10U) shunget(f);
    }
    if (c-'0'>=10U) {
        shunget(f);
        return LLONG_MIN;
    }
    for (x=0; c-'0'<10U && x<INT_MAX/10; c = shgetc(f))
        x = 10*x + c-'0';
    for (y=x; c-'0'<10U && y<LLONG_MAX/100; c = shgetc(f))
        y = 10*y + c-'0';
    for (; c-'0'<10U; c = shgetc(f));
    shunget(f);
    return neg ? -y : y;
}

static long double decfloat(StringFile *f, int c, int bits, int emin, int sign)
{
    uint32_t x[KMAX];
    static const uint32_t th[] = { LD_B1B_MAX };
    int i, j, k, a, z;
    long long lrp=0, dc=0;
    long long e10=0;
    int lnz = 0;
    int gotdig = 0, gotrad = 0;
    int rp;
    int e2;
    int emax = -emin-bits+3;
    int denormal = 0;
    long double y;
    long double frac=0;
    long double bias=0;
    static const int p10s[] = { 10, 100, 1000, 10000,
        100000, 1000000, 10000000, 100000000 };

    j=0;
    k=0;


    /* Don't let leading zeros consume buffer space */
    for (; c=='0'; c = shgetc(f)) gotdig=1;
    if (c=='.') {
        gotrad = 1;
        for (c = shgetc(f); c=='0'; c = shgetc(f)) gotdig=1, lrp--;
    }

    x[0] = 0;
    for (; c-'0'<10U || c=='.'; c = shgetc(f)) {
        if (c == '.') {
            if (gotrad) break;
            gotrad = 1;
            lrp = dc;
        } else if (k < KMAX-3) {
            dc++;
            if (c!='0') lnz = dc;
            if (j) x[k] = x[k]*10 + c-'0';
            else x[k] = c-'0';
            if (++j==9) {
                k++;
                j=0;
            }
            gotdig=1;
        } else {
            dc++;
            if (c!='0') {
                lnz = (KMAX-4)*9;
                x[KMAX-4] |= 1;
            }
        }
    }
    if (!gotrad) lrp=dc;

    if (gotdig && (c|32)=='e') {
        e10 = scanexp(f);
        if (e10 == LLONG_MIN) {
            shunget(f);
            e10 = 0;
        }
        lrp += e10;
    } else if (c>=0) {
        shunget(f);
    }
    if (!gotdig) {
        errno = EINVAL;
        // shlim(f, 0);
        return 0;
    }

    /* 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 || x[0]>>bits==0))
        return sign * (long double)x[0];
    if (lrp > -emin/2) {
        errno = ERANGE;
        return sign * LDBL_MAX * LDBL_MAX;
    }
    if (lrp < emin-2*LDBL_MANT_DIG) {
        errno = ERANGE;
        return sign * LDBL_MIN * LDBL_MIN;
    }

    /* Align incomplete final B1B digit */
    if (j) {
        for (; j<9; j++) x[k]*=10;
        k++;
        j=0;
    }

    a = 0;
    z = k;
    e2 = 0;
    rp = lrp;

    /* Optimize small to mid-size integers (even in exp. notation) */
    if (lnz<9 && lnz<=rp && rp < 18) {
        if (rp == 9) return sign * (long double)x[0];
        if (rp < 9) return sign * (long double)x[0] / p10s[8-rp];
        int bitlim = bits-3*(int)(rp-9);
        if (bitlim>30 || x[0]>>bitlim==0)
            return sign * (long double)x[0] * p10s[rp-10];
    }

    /* Drop trailing zeros */
    for (; !x[z-1]; z--);

    /* Align radix point to B1B digit boundary */
    if (rp % 9) {
        int rpm9 = rp>=0 ? rp%9 : rp%9+9;
        int p10 = p10s[8-rpm9];
        uint32_t carry = 0;
        for (k=a; k!=z; k++) {
            uint32_t 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*LD_B1B_DIG || (rp == 9*LD_B1B_DIG && x[a]<th[0])) {
        uint32_t carry = 0;
        e2 -= 29;
        for (k=(z-1 & MASK); ; k=(k-1 & MASK)) {
            uint64_t tmp = ((uint64_t)x[k] << 29) + carry;
            if (tmp > 1000000000) {
                carry = tmp / 1000000000;
                x[k] = tmp % 1000000000;
            } else {
                carry = 0;
                x[k] = 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 */
    for (;;) {
        uint32_t carry = 0;
        int sh = 1;
        for (i=0; i<LD_B1B_DIG; i++) {
            k = (a+i & MASK);
            if (k == z || x[k] < th[i]) {
                i=LD_B1B_DIG;
                break;
            }
            if (x[a+i & MASK] > th[i]) break;
        }
        if (i==LD_B1B_DIG && rp==9*LD_B1B_DIG) break;
        /* FIXME: find a way to compute optimal sh */
        if (rp > 9+9*LD_B1B_DIG) sh = 9;
        e2 += sh;
        for (k=a; k!=z; k=(k+1 & MASK)) {
            uint32_t 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 */
    for (y=i=0; i<LD_B1B_DIG; i++) {
        if ((a+i & MASK)==z) x[(z=(z+1 & MASK))-1] = 0;
        y = 1000000000.0L * y + x[a+i & MASK];
    }

    y *= sign;

    /* Limit precision for denormal results */
    if (bits > LDBL_MANT_DIG+e2-emin) {
        bits = LDBL_MANT_DIG+e2-emin;
        if (bits<0) bits=0;
        denormal = 1;
    }

    /* Calculate bias term to force rounding, move out lower bits */
    if (bits < LDBL_MANT_DIG) {
        bias = copysignl(scalbn(1, 2*LDBL_MANT_DIG-bits-1), y);
        frac = fmodl(y, scalbn(1, LDBL_MANT_DIG-bits));
        y -= frac;
        y += bias;
    }

    /* Process tail of decimal input so it can affect rounding */
    if ((a+i & MASK) != z) {
        uint32_t t = x[a+i & MASK];
        if (t < 500000000 && (t || (a+i+1 & MASK) != z))
            frac += 0.25*sign;
        else if (t > 500000000)
            frac += 0.75*sign;
        else if (t == 500000000) {
            if ((a+i+1 & MASK) == z)
                frac += 0.5*sign;
            else
                frac += 0.75*sign;
        }
        if (LDBL_MANT_DIG-bits >= 2 && !fmodl(frac, 1))
            frac++;
    }

    y += frac;
    y -= bias;

    if ((e2+LDBL_MANT_DIG & INT_MAX) > emax-5) {
        if (fabsl(y) >= 2/LDBL_EPSILON) {
            if (denormal && bits==LDBL_MANT_DIG+e2-emin)
                denormal = 0;
            y *= 0.5;
            e2++;
        }
        if (e2+LDBL_MANT_DIG>emax || (denormal && frac))
            errno = ERANGE;
    }

    return scalbnl(y, e2);
}

long double __floatscan(StringFile *f, int prec)
{
    int sign = 1;
    size_t i;
    int bits;
    int emin;
    int c;

    switch (prec) {
    case 0:
        bits = FLT_MANT_DIG;
        emin = FLT_MIN_EXP-bits;
        break;
    case 1:
        bits = DBL_MANT_DIG;
        emin = DBL_MIN_EXP-bits;
        break;
    case 2:
        bits = LDBL_MANT_DIG;
        emin = LDBL_MIN_EXP-bits;
        break;
    default:
        return 0;
    }

    while (isspace((c=shgetc(f))));

    if (c=='+' || c=='-') {
        sign -= 2*(c=='-');
        c = shgetc(f);
    }

    for (i=0; i<8 && (c|32)=="infinity"[i]; i++)
        if (i<7) c = shgetc(f);
    if (i==3 || i==8 || (i>3 && true)) {
        if (i!=8) {
            shunget(f);
            if (true) for (; i>3; i--) shunget(f);
        }
        return sign * INFINITY;
    }
    if (!i) for (i=0; i<3 && (c|32)=="nan"[i]; i++)
        if (i<2) c = shgetc(f);
    if (i==3) {
        if (shgetc(f) != '(') {
            shunget(f);
            return NAN;
        }
        for (i=1; ; i++) {
            c = shgetc(f);
            if (c-'0'<10U || c-'A'<26U || c-'a'<26U || c=='_')
                continue;
            if (c==')') return NAN;
            shunget(f);
            while (i--) shunget(f);
            return NAN;
        }
        return NAN;
    }

    if (i) {
        shunget(f);
        errno = EINVAL;
        // shlim(f, 0);
        return 0;
    }

    if (c=='0') {
        c = shgetc(f);
        if ((c|32) == 'x') {
            // hex floats aren't supported
            return 0;
        }
        shunget(f);
        c = '0';
    }

    return decfloat(f, c, bits, emin, sign);
}

#define LARGEST_INT64 INT64_MAX
#define LONGDOUBLE_TYPE long double

/*
** Compute 10 to the E-th power.  Examples:  E==1 results in 10.
** E==2 results in 100.  E==50 results in 1.0e50.
**
** This routine only works for values of E between 1 and 341.
*/
static LONGDOUBLE_TYPE sqlite3Pow10(int E){
#if defined(_MSC_VER)
  static const LONGDOUBLE_TYPE x[] = {
    1.0e+001,
    1.0e+002,
    1.0e+004,
    1.0e+008,
    1.0e+016,
    1.0e+032,
    1.0e+064,
    1.0e+128,
    1.0e+256
  };
  LONGDOUBLE_TYPE r = 1.0;
  int i;
  assert( E>=0 && E<=307 );
  for(i=0; E!=0; i++, E >>=1){
    if( E & 1 ) r *= x[i];
  }
  return r;
#else
  LONGDOUBLE_TYPE x = 10.0;
  LONGDOUBLE_TYPE r = 1.0;
  while(1){
    if( E & 1 ) r *= x;
    E >>= 1;
    if( E==0 ) break;
    x *= x;
  }
  return r; 
#endif
}

int sqlite3AtoF(const char *z, double *pResult, int length){
  int incr = 1;
  const char *zEnd = z + length;
  /* sign * significand * (10 ^ (esign * exponent)) */
  int sign = 1;    /* sign of significand */
  int64_t s = 0;       /* significand */
  int d = 0;       /* adjust exponent for shifting decimal point */
  int esign = 1;   /* sign of exponent */
  int e = 0;       /* exponent */
  int eValid = 1;  /* True exponent is either not used or is well-formed */
  double result;
  int nDigits = 0;
  int nonNum = 0;  /* True if input contains UTF16 with high byte non-zero */

  *pResult = 0.0;   /* Default return value, in case of an error */

  /* skip leading spaces */
  while( z<zEnd && isspace(*z) ) z+=incr;
  if( z>=zEnd ) return 0;

  /* get sign of significand */
  if( *z=='-' ){
    sign = -1;
    z+=incr;
  }else if( *z=='+' ){
    z+=incr;
  }

  /* copy max significant digits to significand */
  while( z<zEnd && isdigit(*z) && s<((LARGEST_INT64-9)/10) ){
    s = s*10 + (*z - '0');
    z+=incr; nDigits++;
  }

  /* skip non-significant significand digits
  ** (increase exponent by d to shift decimal left) */
  while( z<zEnd && isdigit(*z) ){ z+=incr; nDigits++; d++; }
  if( z>=zEnd ) goto do_atof_calc;

  /* if decimal point is present */
  if( *z=='.' ){
    z+=incr;
    /* copy digits from after decimal to significand
    ** (decrease exponent by d to shift decimal right) */
    while( z<zEnd && isdigit(*z) ){
      if( s<((LARGEST_INT64-9)/10) ){
        s = s*10 + (*z - '0');
        d--;
      }
      z+=incr; nDigits++;
    }
  }
  if( z>=zEnd ) goto do_atof_calc;

  /* if exponent is present */
  if( *z=='e' || *z=='E' ){
    z+=incr;
    eValid = 0;

    /* This branch is needed to avoid a (harmless) buffer overread.  The 
    ** special comment alerts the mutation tester that the correct answer
    ** is obtained even if the branch is omitted */
    if( z>=zEnd ) goto do_atof_calc;              /*PREVENTS-HARMLESS-OVERREAD*/

    /* get sign of exponent */
    if( *z=='-' ){
      esign = -1;
      z+=incr;
    }else if( *z=='+' ){
      z+=incr;
    }
    /* copy digits to exponent */
    while( z<zEnd && isdigit(*z) ){
      e = e<10000 ? (e*10 + (*z - '0')) : 10000;
      z+=incr;
      eValid = 1;
    }
  }

  /* skip trailing spaces */
  while( z<zEnd && isspace(*z) ) z+=incr;

do_atof_calc:
  /* adjust exponent by d, and update sign */
  e = (e*esign) + d;
  if( e<0 ) {
    esign = -1;
    e *= -1;
  } else {
    esign = 1;
  }

  if( s==0 ) {
    /* In the IEEE 754 standard, zero is signed. */
    result = sign<0 ? -(double)0 : (double)0;
  } else {
    /* Attempt to reduce exponent.
    **
    ** Branches that are not required for the correct answer but which only
    ** help to obtain the correct answer faster are marked with special
    ** comments, as a hint to the mutation tester.
    */
    while( e>0 ){                                       /*OPTIMIZATION-IF-TRUE*/
      if( esign>0 ){
        if( s>=(LARGEST_INT64/10) ) break;             /*OPTIMIZATION-IF-FALSE*/
        s *= 10;
      }else{
        if( s%10!=0 ) break;                           /*OPTIMIZATION-IF-FALSE*/
        s /= 10;
      }
      e--;
    }

    /* adjust the sign of significand */
    s = sign<0 ? -s : s;

    if( e==0 ){                                         /*OPTIMIZATION-IF-TRUE*/
      result = (double)s;
    }else{
      /* attempt to handle extremely small/large numbers better */
      if( e>307 ){                                      /*OPTIMIZATION-IF-TRUE*/
        if( e<342 ){                                    /*OPTIMIZATION-IF-TRUE*/
          LONGDOUBLE_TYPE scale = sqlite3Pow10(e-308);
          if( esign<0 ){
            result = s / scale;
            result /= 1.0e+308;
          }else{
            result = s * scale;
            result *= 1.0e+308;
          }
        }else{ assert( e>=342 );
          if( esign<0 ){
            result = 0.0*s;
          }else{
#ifdef INFINITY
            result = INFINITY*s;
#else
            result = 1e308*1e308*s;  /* Infinity */
#endif
          }
        }
      }else{
        LONGDOUBLE_TYPE scale = sqlite3Pow10(e);
        if( esign<0 ){
          result = s / scale;
        }else{
          result = s * scale;
        }
      }
    }
  }

  /* store the result */
  *pResult = result;

  /* return true if number and no extra non-whitespace chracters after */
  return z==zEnd && nDigits>0 && eValid && nonNum==0;
}

int main(int argc, char *argv[]) {
    assert(argc==2);
    const char *num = argv[1];
    int len = strlen(num);
    const double ref_answer = strtod(num, NULL);

    printf("parsing '%s' (strtod f64 bits 0x%llx)\n", num, *(uint64_t *)&ref_answer);


    const int iterations = 1000000;

    double musl_time, sql_time;
    {
        clock_t start = clock();

        double result;
        for (int i=0; i < iterations; i++) {
            sqlite3AtoF(num, &result, len);
            
        }
        if(ref_answer != result) {
            printf("sqlite didn't match: 0x%llx\n", *(uint64_t *)&result);
        }

        sql_time = ((double) (clock() - start)) / CLOCKS_PER_SEC;
    }
    
    {
        clock_t start = clock();

        StringFile f = {.str = num, .len = len};
        double result;
        for (int i=0; i < iterations; i++) {
            f.off=0;
            result = __floatscan(&f, 1);
        }
        if(ref_answer != result) {
            printf("musl didn't match: 0x%llx\n", *(uint64_t *)&result);
        }
        musl_time = ((double) (clock() - start)) / CLOCKS_PER_SEC;
    }

    printf("sqlite %gs\n", sql_time);
    printf("musl %gs\n", musl_time);
}

Test Notes
This source is the original f64 versions of the parsers, but I think sqlite will still be faster for f128 zig adaptations.
Performed on a 10 Core M1 Pro, (arm64, macos).

[andrewrk]

If it's faster and still correct that is a clear win. Nice!

[andrewrk]

x64 in rosetta mode

I'm taking this opportunity to make a general statement: I don't care about perf of x64 in rosetta mode. It's completely irrelevant to the near future of computing, in which Apple deprecates x86 entirely and we are left with only native ARM software on macOS.

Zig is young enough that we can aim a little bit for the future and not the present, and this is one of those times. arm64 perf matters. Native x86_64 perf matters. Rosetta-emulated x86_64 perf does not matter.

Same deal with effort spent on compatibility, QA, testing, etc. It's all a waste of time. Given that we have limited efforts to spend, we should spend our efforts on native x86_64 and native arm64 and forget that Rosetta exists. I personally don't even have Rosetta installed.

[andrewrk]

Closing stale PR. Feel free to open a new one if you decide to complete this.