编写自己的math.floor函数实现，使用C语言。

Question

我在思考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

其中一个最终是正确的吗？

- darias

---

floor(x) 返回 小于或等于x的最大整数值。如果您的函数已经这样做了，那么它现在是正确的。 - Scott Hunter

你可以尝试调用my_floor_2(1e100)函数，这将以某种方式告诉你它与原来的不同。 - WhatsUp

Answer 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类型多，这在大多数现代系统上是正确的。

- chqrlie

---

角落：num >= LLONG_MAX 存在舍入模式问题，因为 LLONG_MAX 通常不能精确转换为 double。替代方案：(num >= (LLONG_MAX/2 + 1)*2.0)。 - chux - Reinstate Monica

代码可以通过两个测试来捕获NAN情况：if (!(num < (LLONG_MAX/2 + 1)*2.0 && num >= LLONG_MIN)) { ... - chux - Reinstate Monica

@chux-ReinstateMonica：在num >= LLONG_MAX中的舍入问题，num在这种情况下是无关紧要的：如果long long具有比类型double更多的值位，则如果测试成功，则num是一个整数，如果失败，则num将正确转换为long long。要使转换失败需要num具有大于LLONG_MAX的值，这必然大于或等于(double)LLONG_MAX。 - chqrlie

在这种情况下的使用是足够公平的。我的一般担忧来自于暗示的 num >= (double) LLONG_MAX 和 (double) LLONG_MAX 变成了一个与 LLONG_MAX 不同的值。 - chux - Reinstate Monica

Answer 2

不行，你不能这样去解决。编写自己的实现的最佳方法是使用您平台上的C标准库中的实现。但请注意，它可能包含特定于平台的细微差别，因此可能不可移植。

C标准库的floor函数通常很聪明，因为它不通过将浮点数转换为整型来工作。如果这样做，则有可能出现有符号整数溢出的风险，其行为未定义。（请注意，int的最小可能范围为-32767到+32767）。

精确的实现也取决于您平台所使用的浮点计算方案。

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

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

- Bathsheba

---

将数据类型转换为 long long 时，当 num < 0 时不会返回 _floor_。 - chux - Reinstate Monica

这种方法在 isnan(num) 的情况下也会失败。第一种情况需要完善。 - chux - Reinstate Monica

“编写自己的实现的最佳方法是从您所在平台的C标准库中窃取一个。”<-只有在您接受特定于平台的实现时才是正确的。 - einpoklum

Answer 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;
    }
//---------------------------------------------------------------------------

重点是制作一个掩码，以清除尾数中的小数位，并在负输入和非零清除位的情况下将结果减少。要访问单个位，可以使用联合将浮点值转换为整数表示（如示例中所示），或者改用指针。

我在简单的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常量。此代码未经过优化，因此易于理解，因此不要指望它会很快。

- Spektre

Answer 4

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

检查函数是否存在超出intmax_t范围、NAN、无穷大和-0.0的值，并根据需要进行调整。

#if DBL_MANT_DIG >= 64
#error TBD code
#endif
#include <inttypes.h>

// INTMAX_MAX is not exact as a double, yet INTMAX_MAX + 1 is an exact double
#define INTMAX_MAX_P1 ((INTMAX_MAX/2 + 1)*2.0)

double my_floor(double x) {
  if (x >= 0.0) {
    if (x < INTMAX_MAX_P1) {
      return (double)(intmax_t)x;
    }
    return x;
  } else if (x < 0.0) {
    if (x >= INTMAX_MIN) {
      intmax_t ix = (intmax_t) x; 
      return (ix == x) ? x : (double)(ix-1);
    }
    return x;
  }
  return x; // NAN
}

- chux - Reinstate Monica

Answer 5

在线尝试！。将此作为您的函数尝试：

// we need as much space as possible
typedef long double largestFloat;

largestFloat myFloor(largestFloat x)
{
    largestFloat xcopy = (x < 0) ? (x * -1) : x;
    unsigned int zeros = 0;
    largestFloat n = 1;
    
    // Count digits before the decimal
    for (n = 1; xcopy > (n * 10); n *= 10, ++zeros)
        ;
    
    // Make xcopy follow 0 <= xcopy < 1
    for (xcopy -= n; zeros != -1; xcopy -= n) {
        if (xcopy < 0) {
            xcopy += n;
            n /= 10;
            --zeros;
        }
    }
    xcopy += n;
    
    // Follow standard floor behavior
    if (x < 0)
        return (xcopy == 0) ? x : (x + xcopy - 1);
    else
        return x - xcopy;
}

这是代码的解释。

1. 创建 xcopy（x 的绝对值）
2. 使用第一个 for 循环来确定小数点前的数字数量。
3. 使用该数字不断减少 xcopy，直到满足 0 <= xcopy < 1
4. 根据 x 最初是正数还是负数，返回 x - xcopy 或 x - (1 - xcopy)。

- Lakshya Raj