在哪里可以找到世界上最快的atof实施方案?

时间:2020-03-06 14:24:15  来源:igfitidea点击:

我正在寻找针对美国语言环境,ASCII和非科学符号进行了优化的IA32上非常快速的atof()实现。 Windows多线程CRT在这里很痛苦地崩溃,因为它在每次调用isdigit()时检查区域设置的更改。我们目前的最佳表现来自perl + tcl的atof实现的最佳表现,并且比msvcrt.dll的atof表现高出一个数量级。我想做得更好,但是没有主意。与BCD相关的x86指令似乎很有希望,但是我无法胜过perl / tcl C代码。任何SO'ers都可以挖掘出与那里最好的人的联系吗?也欢迎非基于x86程序集的解决方案。

根据初步答案作出的澄清:

对于该应用,大约2 ulp的误差是可以的。
要转换的数字将以小批量的形式通过网络以ascii消息的形式到达,我们的应用程序需要以尽可能低的延迟来转换它们。

解决方案

我记得我们有一个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

当然,测试一个数字的正确性以及负数将使其变得更加复杂。但是,如果我们可以"假设"输入正确,则可以使代码更加简单和快捷。

准确度要求是什么?如果我们真正需要它"正确"(总是获得最接近指定小数点的浮点值),则可能很难击败标准库版本(除了删除已经支持的语言环境之外),因为这需要进行任意精度的算术运算。如果我们愿意容忍一个或者两个以上的错误(并且比次标准错误更多),那么cruzer提出的这种方法可以行得通,并且可能更快,但绝对不会产生<0.5ulp的输出。我们将在精度方面做得更好,分别计算整数部分和小数部分,最后计算小数(例如,对于12345.6789,将其计算为12345 + 6789 / 10000.0,而不是6 * .1 + 7 * .01 + 8 * .001 + 9 * 0.0001),因为0.1是不合理的二进制分数,并且当我们计算0.1 ^ n时,误差会迅速累积。这也使我们可以使用整数而不是浮点数来进行大多数数学运算。

自从(IIRC)286起,BCD指令就没有在硬件中实现,如今已被微码化。它们不太可能具有很高的性能。

我们是否考虑过让GPU来完成这项工作?如果我们可以将字符串加载到GPU内存中并对其进行处理,那么我们可能会发现一种好的算法,其运行速度将比处理器快得多。

另外,也可以在FPGA中进行操作。有一些FPGA PCI-E板可用于制造任意协处理器。使用DMA将FPGA指向包含我们要转换的字符串数组的内存部分,并使其遍历它们,将转换后的值保留下来。

我们是否看过四核处理器?在大多数情况下,真正的瓶颈是内存访问。

-亚当

我已经实施了一些我们可能会觉得有用的东西。
与atof相比,它的速度快了大约x5,如果与__forceinline一起使用,则快了大约x10.
另一件好事是,它接缝的算法与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;
}