我可以帮你翻译成中文:我在哪里可以找到世界上最快的atof实现?

20
我正在寻找一个针对IA32优化、面向美国英语环境、ASCII字符集及非科学计数法的极速atof()实现。Windows多线程CRT在此处效率极低,因为每次调用isdigit()时都会检查区域设置更改。我们目前最好的实现是基于perl和tcl的atof实现,并且比msvcrt.dll的atof快一个数量级。我想做得更好,但已经没有更好的方法了。BCD相关的x86指令似乎很有希望,但我无法使它的性能超过perl/tcl C代码。有没有SO用户能找到最好的链接?非基于x86汇编的解决方案也可以考虑。需要澄清的是,对于这个应用程序,误差约为2 ulp是可以接受的。要转换的数字将以ASCII消息的形式通过网络以小批量到达,我们的应用程序需要在尽可能低的延迟下进行转换。

1
检查 isdigit 的区域设置更改?也许他们应该查看 ISO C 标准。isdigit 没有区域设置相关的行为;它必须检查字符是否是集合 09 中的元素,仅此而已。 - Kaz
你能给我们一个问题域的概念吗?我猜它不是金融领域,否则你会使用定点算术。这是用于控制系统,例如定位吗?你有(硬件或软件)实时要求吗? - Toby Speight
1
如果您可以修改消息格式,显然发送二进制浮点数(或二进制的简单文本编码)将节省在另一端昂贵的解析。例如,如果二进制不行,可以将浮点数转储为十六进制整数。 - Peter Cordes
https://github.com/lemire/fast_double_parser(来自链接式答案)。Daniel Lemire是优化领域的知名人士,因此很可能具有高质量。另请参阅他相关的博客文章:https://lemire.me/blog/2020/09/10/parsing-floats-in-c-benchmarking-strtod-vs-from_chars/。 - Peter Cordes
6个回答

11
你的精度要求是多少?如果你真的需要它是“正确的”(总是得到最接近指定小数的浮点数),那么很难超越标准库版本(除了删除语言环境支持,这已经做过了),因为这需要进行任意精度算术运算。如果你可以容忍一个或两个ulp的误差(对于子规格而言可能更多),cruzer提出的方法可能有效并且速度更快,但它肯定不会产生<0.5ulp的输出。从精度的角度来看,你最好将整数和小数部分分别计算,并在最后计算分数(例如,对于12345.6789,将其计算为12345 + 6789 / 10000.0,而不是6* .1 + 7*.01 + 8*.001 + 9 * 0.0001),因为0.1是一个无理二进制分数,并且当你计算0.1^n时误差会迅速累积。这样还可以使用整数而不是浮点数进行大部分计算。

自 (我IRC) 286以来,BCD指令就没有在硬件中实现了,现在只是微码化的。它们不太可能特别高效.


你提出的方法实际上是我们从perl/tcl中获取的想法所做的。正如你所提到的,它更快、更精确。谢谢你关于BCD历史的提示 - 我之前没有听说过,但循环计数似乎非常高。 - Dave Moore
BCD指令(如AAM等)甚至在64位模式下都不存在,这是不使用它们的另一个原因。 - Peter Cordes

6
我刚完成的这个实现在我的桌面电脑上比内置的'atof'快两倍。它可以在2秒内将1024 * 1024 * 39个数字输入转换,而系统标准的gnu 'atof'需要4秒。(包括设置时间、获取内存等等)。
更新: 抱歉,我必须撤回我的两倍速度的说法。如果你要转换的东西已经在一个字符串中,那么它会更快,但是如果你传递的是硬编码的字符串文字,它就和atof差不多了。然而,我还是会把它留在这里,因为可能通过对ragel文件和状态机进行一些调整,你可以为特定目的生成更快的代码。

https://github.com/matiu2/yajp

您感兴趣的文件为:

https://github.com/matiu2/yajp/blob/master/tests/test_number.cpp

https://github.com/matiu2/yajp/blob/master/number.hpp

此外,您可能会对执行转换的状态机感兴趣:

Number parsing state machine diagram


5
我觉得你想手动构建一个状态机,其中每个状态处理第N个输入数字或指数数字;这个状态机将形状像一棵树(没有循环!)。目标是尽可能地进行整数运算,并且(显然)在状态中隐含地记住状态变量(“前导负号”,“小数点位于位置3”),以避免赋值、存储和稍后获取/测试此类值。仅使用输入字符的普通“if”语句实现状态机(因此您的树将成为一组嵌套的ifs)。内联访问缓冲区字符;您不希望调用getchar函数使速度变慢。
可以简单地抑制前导零;您可能需要一个循环来处理非常长的前导零序列。第一个非零数字可以收集而无需将累加器清零或乘以十。前4-9个非零数字(对于16位或32位整数)可以通过与常量值十相乘的整数乘法来收集(由大多数编译器转换为几次移位和加法)。【过度:直到找到非零数字并且需要连续N个零的乘以10^N,所有这些都可以连接到状态机中】。后面的数字可以使用32位或64位乘法收集,具体取决于您的机器字长。由于您不关心准确性,因此在收集了32位或64位的数字后,可以简单地忽略后面的数字;我猜您实际上可以根据应用程序对这些数字的实际处理情况停止某个固定数量的非零数字。在数字字符串中找到小数点只需在状态机树中引起分支。该分支知道小数点的隐含位置,因此稍后如何适当地按10的幂进行缩放。如果您不喜欢此代码的大小,可以努力将一些状态机子树组合在一起。
【过度:将整数和小数部分保持为单独的(小)整数。这将需要在最后进行额外的浮点运算以组合整数和小数部分,可能不值得】。
【过度:将2个字符收集为数字对的16位值,查找16位值。这避免了寄存器中的乘法,但要换取内存访问,可能在现代机器上没有优势】。
遇到“E”时,像上面那样将指数收集为整数;在预先计算/缩放的十次幂表中查找准确的乘数(如果指数中存在“-”号,则为倒数),并将收集的尾数乘以该乘数。(永远不要执行浮点除法)。由于每个指数收集例程都在树的不同分支(叶子)中,因此必须通过偏移十次幂索引来调整小数点的表面或实际位置。

[越过顶端:如果您知道数字字符线性存储在缓冲区中并且不跨越缓冲区边界,则可以避免ptr ++的成本。 在树枝的第k个状态下,您可以将第k个字符访问为*(start+k)。 一个好的编译器通常可以隐藏寻址模式中的“... + k”索引偏移量。]

如果正确完成,此方案每个非零数字大约需要一次廉价的乘加法操作,一个浮点数强制转换操作和一个浮点数乘法操作以按指数和小数点位置进行比例缩放结果。

我没有实现上述内容。 我已经实现了循环版本,它们非常快速。


3

我实现了一些你可能会发现有用的东西。

与atof相比,它大约快了5倍,如果与__forceinline一起使用,则大约快了10倍。

另一个好处是它似乎具有与crt实现完全相同的算术。

当然,它也有一些缺点:

  • 它仅支持单精度浮点数,
  • 并且不扫描任何特殊值,如#INF等...
__forceinline bool float_scan(const wchar_t* wcs, float* val)
{
int hdr=0;
while (wcs[hdr]==L' ')
    hdr++;

int cur=hdr;

bool negative=false;
bool has_sign=false;

if (wcs[cur]==L'+' || wcs[cur]==L'-')
{
    if (wcs[cur]==L'-')
        negative=true;
    has_sign=true;
    cur++;
}
else
    has_sign=false;

int quot_digs=0;
int frac_digs=0;

bool full=false;

wchar_t period=0;
int binexp=0;
int decexp=0;
unsigned long value=0;

while (wcs[cur]>=L'0' && wcs[cur]<=L'9')
{
    if (!full)
    {
        if (value>=0x19999999 && wcs[cur]-L'0'>5 || value>0x19999999)
        {
            full=true;
            decexp++;
        }
        else
            value=value*10+wcs[cur]-L'0';
    }
    else
        decexp++;

    quot_digs++;
    cur++;
}

if (wcs[cur]==L'.' || wcs[cur]==L',')
{
    period=wcs[cur];
    cur++;

    while (wcs[cur]>=L'0' && wcs[cur]<=L'9')
    {
        if (!full)
        {
            if (value>=0x19999999 && wcs[cur]-L'0'>5 || value>0x19999999)
                full=true;
            else
            {
                decexp--;
                value=value*10+wcs[cur]-L'0';
            }
        }

        frac_digs++;
        cur++;
    }
}

if (!quot_digs && !frac_digs)
    return false;

wchar_t exp_char=0;

int decexp2=0; // explicit exponent
bool exp_negative=false;
bool has_expsign=false;
int exp_digs=0;

// even if value is 0, we still need to eat exponent chars
if (wcs[cur]==L'e' || wcs[cur]==L'E')
{
    exp_char=wcs[cur];
    cur++;

    if (wcs[cur]==L'+' || wcs[cur]==L'-')
    {
        has_expsign=true;
        if (wcs[cur]=='-')
            exp_negative=true;
        cur++;
    }

    while (wcs[cur]>=L'0' && wcs[cur]<=L'9')
    {
        if (decexp2>=0x19999999)
            return false;
        decexp2=10*decexp2+wcs[cur]-L'0';
        exp_digs++;
        cur++;
    }

    if (exp_negative)
        decexp-=decexp2;
    else
        decexp+=decexp2;
}

// end of wcs scan, cur contains value's tail

if (value)
{
    while (value<=0x19999999)
    {
        decexp--;
        value=value*10;
    }

    if (decexp)
    {
        // ensure 1bit space for mul by something lower than 2.0
        if (value&0x80000000)
        {
            value>>=1;
            binexp++;
        }

        if (decexp>308 || decexp<-307)
            return false;

        // convert exp from 10 to 2 (using FPU)
        int E;
        double v=pow(10.0,decexp);
        double m=frexp(v,&E);
        m=2.0*m;
        E--;
        value=(unsigned long)floor(value*m);

        binexp+=E;
    }

    binexp+=23; // rebase exponent to 23bits of mantisa


    // so the value is: +/- VALUE * pow(2,BINEXP);
    // (normalize manthisa to 24bits, update exponent)
    while (value&0xFE000000)
    {
        value>>=1;
        binexp++;
    }
    if (value&0x01000000)
    {
        if (value&1)
            value++;
        value>>=1;
        binexp++;
        if (value&0x01000000)
        {
            value>>=1;
            binexp++;
        }
    }

    while (!(value&0x00800000))
    {
        value<<=1;
        binexp--;
    }

    if (binexp<-127)
    {
        // underflow
        value=0;
        binexp=-127;
    }
    else
    if (binexp>128)
        return false;

    //exclude "implicit 1"
    value&=0x007FFFFF;

    // encode exponent
    unsigned long exponent=(binexp+127)<<23;
    value |= exponent;
}

// encode sign
unsigned long sign=negative<<31;
value |= sign;

if (val)
{
    *(unsigned long*)val=value;
}

return true;
}

2
这段代码看起来很费力地追求精确,但最终失败了。对于非常小的decexp值,pow(10.0,decexp)是不精确的。而且你好像把次标准结果发送到0,而不是它们的值。如果我错了,请纠正我。 - R.. GitHub STOP HELPING ICE

2

我记得我们有一个Winforms应用程序,在解析一些数据交换文件时表现非常缓慢,我们都认为是数据库服务器出现了问题,但我们聪明的老板实际上发现瓶颈在于将解析的字符串转换为十进制数的调用!

最简单的方法是循环遍历字符串中的每个数字(字符),保持一个运行总数,将总数乘以10,然后加上下一个数字的值。继续这样做,直到您到达字符串的末尾或遇到一个点。如果遇到一个点,请将整数部分与小数部分分开,然后有一个乘法器,每个数字都会将自己除以10。随着你的前进不断地添加它们。

例如:123.456

累加总数=0,加1(现在是1) 累加总数=1 * 10 = 10,加2(现在是12) 累加总数=12 * 10 = 120,加3(现在是123) 遇到小数点,准备处理小数部分 乘数=0.1,乘以4,得到0.4,加到累加总数中,变成了123.4 乘数=0.1 / 10 = 0.01,乘以5,得到0.05,加到累加总数中,变成了123.45 乘数=0.01 / 10 = 0.001,乘以6,得到0.006,加到累加总数中,变成了123.456

当然,测试数字的正确性以及负数会使它更加复杂。但是,如果您可以“假设”输入是正确的,那么您可以使代码更简单、更快速。


1
我认为你应该只处理小数部分一次。收集小数点之前的所有内容,当你到达小数点时,开始一个新数字和数字1。将两个数字都乘以10。最后,你有3个浮点数,整数部分,小数部分作为整数和使其成为小数所需的除数。现在进行除法运算并相加:intpart + decpart/divisor。 - jmucchiello

0

你有没有考虑让GPU来处理这个任务?如果你可以将字符串加载到GPU内存中并让它们进行处理,你可能会发现一个好的算法可以比你的处理器运行得更快。

或者,你可以使用FPGA来完成这项工作——有一些FPGA PCI-E板可以用来制作任意协处理器。使用DMA将FPGA指向包含要转换的字符串数组的内存部分,让它快速地完成转换并留下转换后的值。

你有考虑过四核处理器吗?在大多数情况下,真正的瓶颈都是内存访问...

-Adam


2
这些浮点数经常通过网络到达,但时间是随机的。我们强调延迟而不是吞吐量,否则我认为您的GPU或FPGA建议会很可靠。我们使用8和16核CPU,并且像上面提到的循环展开+转换一样,内存/缓存很重要。 - Dave Moore

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接