我怎样才能在不将它们转换为 double/long double 的情况下以不错的精度实现 fmod 函数？-解网

问：

我正在尝试在汇编中实现fmodf函数。我使用了如下公式。

fmodf（浮点数 x，浮点数 y） = x - 截断（x / y） * y

问题是，当 x>>y 时，函数的精度会越来越差。事实上，当 abs（x/y）>2^23 时，trunc 不会去掉小数部分，因为它没有小数部分（我知道 IEEE-754 浮点标准的尾数是 23 位）

我已经研究了 cmath 中 fmod 的实现，它只是使用了 __built_in_remainderl基本上使用 long double 而不是 float，我无法在我的汇编中使用它（它不支持 double 或 long double）

算法浮点精度计算

fmod()如果性能不重要，可以通过简单的二进制长手除法以直接的方式实现。长手除法是我们在小学学到的除法过程，每步产生一个商数。当使用基数二时，它的工作方式相同，每步生成一位商。与任何长手除法一样，它提供了 ISO-C 和 ISO-C++ 规范要求的确切余数。fmod()

二进制除法过程可以采用整数算术或浮点运算法执行。Asker 没有说明其平台的细节，但提到对单精度计算的支持表明它可能是一个具有低延迟硬件支持的处理器，支持使用该格式的单精度 IEEE-754 算法，例如 ARM Cortex-M4。因此，在下面的 ISO C99 代码中，我将演示基于浮点的计算。binary32

代码首先处理一些涉及 NaN、无穷大和零的特殊情况。如果除数的大小大于被除数的大小，则结果为被除数。否则，我们必须通过使除数左对齐来准备除法过程，即将除数与除数达到相同的数量级。在每个除法步骤中，下一个商位基于部分余数与除数的比较。如果它大于或等于除数，则商位为 1，我们需要从部分余数中减去除数。由于商本身从来不需要，所以它没有被记录下来。my_fmodf()fmod()

该代码假定映射到 IEEE-754 格式，并且用于 32 位整数数据和 32 位浮点数据的存储约定相同。它是在支持次正态操作数的平台上开发和“冒烟”测试的，因此它可能在非IEEE754刷新归零（FTZ）状态下正常工作，也可能不正常工作。我使用英特尔 C/C++ 编译器构建，用于最严格地遵守 IEEE-754 浮点语义。floatbinary32/fp:strict

下面的 C 代码几乎可以一对一地翻译成汇编语言。并且用于操作浮点数的指数字段，并且对于此处使用的有限功能，可以用一些汇编指令替换，每个指令在常见的处理器体系结构上。类似于对浮点操作数底层位模式的简单检查。类似地，如果硬件指令不直接支持，则和是使用整数逻辑指令对浮点操作数的符号位的简单操作。frexpf()ldexpf()isnan()isinf()fabsf()copysignf()

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>

float uint32_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof r);
    return r;
}

uint32_t float_as_uint32 (float a)
{
    uint32_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

/* returns the floating-point remainder of a/b (rounded towards zero) */
float my_fmodf (float a, float b)
{
    const float NAN_INDEFINITE = uint32_as_float (0xffc00000);
    float r;
    if (isnan (a) || isnan (b)) {
        r = a + b;
    } else if (isinf (a) || (b == 0.0f)) {
        r = NAN_INDEFINITE;
    } else {
        float fa, fb, dividend, divisor;
        int expo_a, expo_b;
        fa = fabsf (a);
        fb = fabsf (b);
        if (fa >= fb) {
            dividend = fa;
            /* normalize divisor */
            (void)frexpf (fa, &expo_a);
            (void)frexpf (fb, &expo_b);
            divisor = ldexpf (fb, expo_a - expo_b);
            if (divisor <= 0.5f * dividend) {
                divisor += divisor;
            }
            /* compute quotient one bit at a time */
            while (divisor >= fb) {
                if (dividend >= divisor) {
                    dividend -= divisor;
                }
                divisor *= 0.5f;
            }
            /* dividend now represents remainder */
            r = copysignf (dividend, a);
        } else {
            r = a;
        }
    }
    return r;
}

// George Marsaglia's KISS PRNG, period 2**123. Newsgroup sci.math, 21 Jan 1999
// Bug fix: Greg Rose, "KISS: A Bit Too Simple" http://eprint.iacr.org/2011/007
static uint32_t kiss_z=362436069, kiss_w=521288629;
static uint32_t kiss_jsr=123456789, kiss_jcong=380116160;
#define znew (kiss_z=36969*(kiss_z&65535)+(kiss_z>>16))
#define wnew (kiss_w=18000*(kiss_w&65535)+(kiss_w>>16))
#define MWC  ((znew<<16)+wnew )
#define SHR3 (kiss_jsr^=(kiss_jsr<<13),kiss_jsr^=(kiss_jsr>>17), \
              kiss_jsr^=(kiss_jsr<<5))
#define CONG (kiss_jcong=69069*kiss_jcong+1234567)
#define KISS ((MWC^CONG)+SHR3)

int main (void)
{
    double a, b, res, ref;
    uint64_t i = 0;
    do {
        a = uint32_as_float (KISS);
        b = uint32_as_float (KISS);
        res = my_fmodf (a, b);
        ref = fmodf (a, b);

        if (float_as_uint32 (res) != float_as_uint32 (ref)) {
            printf ("error: a=% 15.8e b=% 15.8e res=% 15.8e ref=% 15.8e\n", 
                    a, b, res, ref);
            return EXIT_FAILURE;
        }
        i++;
        if (!(i & 0xfffff)) printf ("\r%llu", i);
    } while (i);
    return EXIT_SUCCESS;
}

上一个：Javascript 精度，小数点后 15 位以上

下一个：Python，浮点数可以有多长？

我怎样才能在不将它们转换为 double/long double 的情况下以不错的精度实现 fmod 函数？

How can I implement a fmod function with decent precision without casting them into double/long double?

评论