编写自己的数学 floor 函数实现，C答案

【问题标题】：Write your own implementation of math's floor function, C编写自己的数学 floor 函数实现，C
【发布时间】：2017-01-25 16:31:44
【问题描述】：

我在考虑math.h 中提供的floor 函数。使用起来非常简单：

#include <stdio.h>
#include <math.h>
 
int main(void)
{
  for (double a = 12.5; a < 13.4; a += 0.1)
    printf("floor of  %.1lf is  %.1lf\n", a, floor(a));
  return 0;
}

如果我想编写自己的实现怎么办？会不会是这样的：

#include <stdio.h>
#include <math.h>

double my_floor(double num)
{
    return (int)num;
}

int main(void)
{
    double a;

    for (a = 12.5; a < 13.4; a += 0.1)
        printf("floor of  %.1lf is  %.1lf\n", a, floor(a));

    printf("\n\n");

    for (a = 12.5; a < 13.4; a += 0.1)
        printf("floor of  %.1lf is  %.1lf\n", a, my_floor(a));

    return 0;
}

?

它似乎不适用于负数（my_floor），但第二个似乎没问题（my_floor_2）：

#include <stdio.h>
#include <math.h>

double my_floor(double num)
{
    return (int)num;
}

double my_floor_2(double num)
{
    if(num < 0)
        return (int)num - 1;
    else
        return (int)num;
}

int main(void)
{
    double a1 = -12.5;

    printf("%lf\n", floor(a1));
    printf("%lf\n", my_floor(a1));
    printf("%lf\n", my_floor_2(a1));

    return 0;
}

程序输出：

-13.000000
-12.000000
-13.000000

其中一个最终是否正确？

【问题讨论】：

floor(x) 产生the largest integer value less than or equal to x。如果您这样做，那么现在是正确的。
你可以试试my_floor_2(1e100)，它会以某种方式告诉你区别，

标签： c math floor

【解决方案1】：

您的两种尝试都有局限性：

如果double 值超出int 类型的范围，则转换为int 是实现定义的。
如果double 值为负但整数，则返回(int)num - 1 不正确。

这是一个（几乎）可移植的版本，它试图处理所有情况：

double my_floor_2(double num) {
    if (num >= LLONG_MAX || num <= LLONG_MIN || num != num) {
        /* handle large values, infinities and nan */
        return num;
    }
    long long n = (long long)num;
    double d = (double)n;
    if (d == num || num >= 0)
        return d;
    else
        return d - 1;
}

如果类型 long long 的值位比类型 double 多，则应该是正确的，大多数现代系统都是这种情况。

【讨论】：

Corner: num >= LLONG_MAX 存在舍入模式问题，因为 LLONG_MAX 通常不会完全转换为 double。备选：(num >= (LLONG_MAX/2 + 1)*2.0)
代码可以通过 2 个测试捕获 NAN 案例：if (!(num < (LLONG_MAX/2 + 1)*2.0 && num >= LLONG_MIN)) { ...
@chux-ReinstateMonica：在这种特殊情况下，num >= LLONG_MAX、num 中的舍入问题是无关紧要的：如果long long 的值位比double 类型多，如果测试成功， num 是一个整数，如果失败，num 将正确转换为long long。要使转换失败，num 的值必须大于 LLONG_MAX，而 LLONG_MAX 必然大于或等于 (double)LLONG_MAX。
在这种情况下使用足够公平。我的普遍担忧源于隐含的num >= (double) LLONG_MAX 和(double) LLONG_MAX 成为不同于LLONG_MAX 的值。

【解决方案2】：

不，你不能这样处理。编写自己的实现的最佳方法是从您平台上的 C 标准库中窃取一个。但请注意，它可能包含特定于平台的细微差别，因此可能无法移植。

C 标准库floor 函数通常很聪明，因为它不能通过转换为整数类型来工作。如果是这样，那么您将面临signed 整数溢出的风险，其行为未定义。（请注意，int 的最小可能范围是 -32767 到 +32767）。

精确的实现还取决于您平台上使用的浮点方案。

对于使用 IEEE754 浮点的平台和 long long 类型，您可以采用这种方案：

如果数字的大小大于 2 的 53 次方，请将其返回（因为它已经是整数）。
否则，转换为 64 位类型（long long），然后返回。

【讨论】：

在num < 0 时，投射到long long 不会返回地板。
该方法在isnan(num) 时也会失败。案例 1. 将受益于细化。
" 编写您自己的实现的最佳方式是从您平台上的 C 标准库中窃取一个。"

【解决方案3】：

在 C++ 和 32 位算术中，例如可以这样做：

//---------------------------------------------------------------------------
// IEEE 754 double MSW masks
const DWORD _f64_sig    =0x80000000;    // sign
const DWORD _f64_exp    =0x7FF00000;    // exponent
const DWORD _f64_exp_sig=0x40000000;    // exponent sign
const DWORD _f64_exp_bia=0x3FF00000;    // exponent bias
const DWORD _f64_exp_lsb=0x00100000;    // exponent LSB
const DWORD _f64_exp_pos=        20;    // exponent LSB bit position
const DWORD _f64_man    =0x000FFFFF;    // mantisa
const DWORD _f64_man_msb=0x00080000;    // mantisa MSB
const DWORD _f64_man_bits=       52;    // mantisa bits
// IEEE 754 single masks
const DWORD _f32_sig    =0x80000000;    // sign
const DWORD _f32_exp    =0x7F800000;    // exponent
const DWORD _f32_exp_sig=0x40000000;    // exponent sign
const DWORD _f32_exp_bia=0x3F800000;    // exponent bias
const DWORD _f32_exp_lsb=0x00800000;    // exponent LSB
const DWORD _f32_exp_pos=        23;    // exponent LSB bit position
const DWORD _f32_man    =0x007FFFFF;    // mantisa
const DWORD _f32_man_msb=0x00400000;    // mantisa MSB
const DWORD _f32_man_bits=       23;    // mantisa bits
//---------------------------------------------------------------------------
double f64_floor(double x)
    {
    const int h=1;      // may be platform dependent MSB/LSB order
    const int l=0;
    union _f64          // semi result
        {
        double f;       // 64bit floating point
        DWORD u[2];     // 2x32 bit uint
        } y;
    DWORD m,a;
    int sig,exp,sh;
    y.f=x;
    // extract sign
    sig =y.u[h]&_f64_sig;
    // extract exponent
    exp =((y.u[h]&_f64_exp)>>_f64_exp_pos)-(_f64_exp_bia>>_f64_exp_pos);
    // floor bit shift
    sh=_f64_man_bits-exp; a=0;
    if (exp<0)
        {
        a=y.u[l]|(y.u[h]&_f64_man);
        if (sig) return -1.0;
        return 0.0;
        }
    // LSW
    if (sh>0)
        {
        if (sh<32) m=(0xFFFFFFFF>>sh)<<sh; else m=0;
        a=y.u[l]&(m^0xFFFFFFFF); y.u[l]&=m;
        }
    // MSW
    sh-=32;
    if (sh>0)
        {
        if (sh<_f64_exp_pos) m=(0xFFFFFFFF>>sh)<<sh; else m=_f64_sig|_f64_exp;
        a|=y.u[h]&(m^0xFFFFFFFF); y.u[h]&=m;
        }
    if ((sig)&&(a)) y.f--;
    return y.f;
    }
//---------------------------------------------------------------------------
float f32_floor(float x)
    {
    union               // semi result
        {
        float f;        // 32bit floating point
        DWORD u;        // 32 bit uint
        } y;
    DWORD m,a;
    int sig,exp,sh;
    y.f=x;
    // extract sign
    sig =y.u&_f32_sig;
    // extract exponent
    exp =((y.u&_f32_exp)>>_f32_exp_pos)-(_f32_exp_bia>>_f32_exp_pos);
    // floor bit shift
    sh=_f32_man_bits-exp; a=0;
    if (exp<0)
        {
        a=y.u&_f32_man;
        if (sig) return -1.0;
        return 0.0;
        }
    if (sh>0)
        {
        if (sh<_f32_exp_pos) m=(0xFFFFFFFF>>sh)<<sh; else m=_f32_sig|_f32_exp;
        a|=y.u&(m^0xFFFFFFFF); y.u&=m;
        }
    if ((sig)&&(a)) y.f--;
    return y.f;
    }
//---------------------------------------------------------------------------

关键是制作掩码，从尾数中清除十进制位，并且在负输入和非零清除位的情况下减少结果。要访问单个位，您可以使用 union 将浮点值转换为整数表示（如示例中所示）或使用指针。

我在简单的 VCL 应用中对此进行了测试，如下所示：

float f32;
double f64;
AnsiString txt="";
// 64 bit
txt+="[double]\r\n";
for (f64=-10.0;f64<=10.0;f64+=0.1)
 if (fabs(floor(f64)-f64_floor(f64))>1e-6)
    {
    txt+=AnsiString().sprintf("%5.3lf %5.3lf %5.3lf\r\n",f64,floor(f64),f64_floor(f64));
    f64_floor(f64);
    }
for (f64=1;f64<=1e307;f64*=1.1)
    {
    if (fabs(floor( f64)-f64_floor( f64))>1e-6) { txt+=AnsiString().sprintf("%lf lf lf\r\n", f64,floor( f64),f64_floor( f64));
    f64_floor( f64); }
    if (fabs(floor(-f64)-f64_floor(-f64))>1e-6) { txt+=AnsiString().sprintf("%lf lf lf\r\n",-f64,floor(-f64),f64_floor(-f64));
    f64_floor(-f64); }
    }
// 32 bit
txt+="[float]\r\n";
for (f32=-10.0;f32<=10.0;f32+=0.1)
 if (fabs(floor(f32)-f32_floor(f32))>1e-6)
    {
    txt+=AnsiString().sprintf("%5.3lf %5.3lf %5.3lf\r\n",f32,floor(f32),f32_floor(f32));
    f32_floor(f32);
    }
for (f32=1;f32<=1e37;f32*=1.1)
    {
    if (fabs(floor( f32)-f32_floor( f32))>1e-6) { txt+=AnsiString().sprintf("%lf lf lf\r\n", f32,floor( f32),f32_floor( f32));
    f32_floor( f32); }
    if (fabs(floor(-f32)-f32_floor(-f32))>1e-6) { txt+=AnsiString().sprintf("%lf lf lf\r\n",-f32,floor(-f32),f32_floor(-f32));
    f32_floor(-f32); }
    }
mm_log->Lines->Add(txt);

没有差异结果（因此在所有测试用例中，它都匹配 math.h floor() 值。如果您想在 VCL 之外尝试一下，那么只需将 AnsiString 更改为任何字符串类型您可以将输出从 TMemo::mm_log 更改为您获得的任何内容（例如控制台 cout 或其他）

fxx_floor() 的双重调用以防出现差异是为了调试目的（您可以在错误情况下直接放置断点并单步执行）。

[备注]

注意单词的顺序 (MSW,LSW) 取决于平台，因此您应该相应地调整 h,l 常量。这段代码没有经过优化，所以很容易理解，所以不要指望它会很快。

【讨论】：

【解决方案4】：

当浮点类型的精度与宽整数类型相比足够小时，当浮点值在整数范围内时强制转换为该整数类型。

查看函数是否有超出long long 范围、NAN、无穷大和 -0.0 沙调整和所需的值。

#if DBL_MANT_DIG >= 64
#error TBD code
#endif

// LLONG_MAX is not exact as a double, yet LLONG_MAX + 1 is
#define LLONG_MAX_P1 ((LLONG_MAX/2 + 1)*2.0)

double my_floor(double x) {
  if (x >= 0.0) {
    if (x < LLONG_MAX_P1) {
      return (double)(long long)x;
    }
    return x;
  } else if (x < 0.0) {
    if (x >= LLONG_MIN) {
      long long ix = (long long) x; 
      return (ix == x) ? x : (double)(ix-1);
    }
    return x;
  }
  return x; // NAN
}

【讨论】：

【解决方案5】：

无需从 C 标准库中复制！这个floor 函数已经在Try It Online (tio.run) 和Online GDB 上进行了测试。该函数本身没有#include 任何文件。这里是：

// the return type of myFloor is supposed to be the largest size
long double myFloor(long double x)
{
    long double xcopy=x<0?x*-1:x;
    unsigned int zeros=0;
    long double n=1;
    for(n=1;xcopy>n*10;n*=10,++zeros);
    for(xcopy-=n;zeros!=-1;xcopy-=n)
        if(xcopy<0)
        {
            xcopy+=n;
            n/=10;
            --zeros;
        }
    xcopy+=n;
    return x<0?(xcopy==0?x:x-(1-xcopy)):(x-xcopy);
}

浮点数的下限是小于或等于它的最大整数。以下是一些示例：

floor(5.7)  = 5
floor(3)    = 3
floor(9.9)  = 9
floor(7.0)  = 7
floor(-7.9) = -8
floor(-5.0) = -5
floor(-3.3) = -3
floor(0)    = 0
floor(-0.0) = -0
floor(-0)   = -0

【讨论】：