4

我想ceil()在C. 在库中搜索源代码并在此处找到，但似乎很难理解。我想要干净优雅的代码。

我也在 SO 上进行了搜索，在这里找到了一些答案。似乎没有一个答案是正确的。答案之一是：

#define CEILING_POS(X) ((X-(int)(X)) > 0 ? (int)(X+1) : (int)(X))
#define CEILING_NEG(X) ((X-(int)(X)) < 0 ? (int)(X-1) : (int)(X))
#define CEILING(X) ( ((X) > 0) ? CEILING_POS(X) : CEILING_NEG(X) )

AFAIK，返回类型ceil()不是 int。宏在这里是类型安全的吗？此外，上述实现是否适用于负数？

实施它的最佳方式是什么？

你能提供干净的代码吗？

4

4 回答 4

2

您引用的宏对于大于INT_MAX但仍可以完全表示为双精度的数字肯定不能正常工作。

正确实现的唯一方法ceil()（假设您不能使用等效的汇编指令实现它）是对浮点数的二进制表示进行位旋转，就像在s_ceil.c第一个链接后面的源文件中所做的那样。理解代码是如何工作的需要理解底层平台的浮点表示——这种表示很可能是IEEE 754——但是没有办法解决这个问题。

编辑：

其中的一些复杂性s_ceil.c源于它处理的特殊情况（NaN、无穷大）以及它需要在无法假设存在 64 位整数类型的情况下完成其工作的事实。

所有位旋转的基本思想是屏蔽尾数的小数位，如果数字大于零，则在其上加 1...但是还涉及一些额外的逻辑以确保您执行在所有情况下都是正确的。

ceil()这是我拼凑在一起的 for floats 的说明性版本。注意：这不能正确处理特殊情况，并且没有经过广泛的测试——所以不要实际使用它。然而，它确实有助于说明比特旋转所涉及的原理。我试图广泛地评论该例程，但这些评论确实假设您了解浮点数是如何以 IEEE 754 格式表示的。

union float_int
{
    float f;
    int i;
};

float myceil(float x)
{
    float_int val;
    val.f=x;

    // Extract sign, exponent and mantissa
    // Bias is removed from exponent
    int sign=val.i >> 31;
    int exponent=((val.i & 0x7fffffff) >> 23) - 127;
    int mantissa=val.i & 0x7fffff;

    // Is the exponent less than zero?
    if(exponent<0) 
    {   
        // In this case, x is in the open interval (-1, 1)
        if(x<=0.0f)
            return 0.0f;
        else
            return 1.0f;
    }
    else
    {
        // Construct a bit mask that will mask off the
        // fractional part of the mantissa
        int mask=0x7fffff >> exponent;

        // Is x already an integer (i.e. are all the
        // fractional bits zero?)
        if((mantissa & mask) == 0)
            return x;
        else
        {
            // If x is positive, we need to add 1 to it
            // before clearing the fractional bits
            if(!sign)
            {
                mantissa+=1 << (23-exponent);

                // Did the mantissa overflow?
                if(mantissa & 0x800000)
                {
                    // The mantissa can only overflow if all the
                    // integer bits were previously 1 -- so we can
                    // just clear out the mantissa and increment
                    // the exponent
                    mantissa=0;
                    exponent++;
                }
            }

            // Clear the fractional bits
            mantissa&=~mask;
        }
    }

    // Put sign, exponent and mantissa together again
    val.i=(sign << 31) | ((exponent+127) << 23) | mantissa;

    return val.f;
}
于 2012-09-05T11:04:03.793 回答
2

没有什么比使用标准库实现更优雅了。没有任何代码总是比优雅的代码更优雅。

除此之外，这种方法有两个主要缺陷：

  • 如果X大于INT_MAX + 1或小于INT_MIN - 1，则宏的行为未定义。这意味着您的实现可能对所有浮点数的近一半给出不正确的结果。与 IEEE-754 相反，您还将引发无效标志。
  • 它得到了 -0、+/-infinity 和 nan 错误的边缘情况。事实上，它唯一正确的边缘情况是+0。

您可以ceil以类似于您尝试的方式实现，就像这样（此实现假定 IEEE-754 双精度）：

#include <math.h>

double ceil(double x) {
    // All floating-point numbers larger than 2^52 are exact integers, so we
    // simply return x for those inputs.  We also handle ceil(nan) = nan here.
    if (isnan(x) || fabs(x) >= 0x1.0p52) return x;
    // Now we know that |x| < 2^52, and therefore we can use conversion to
    // long long to force truncation of x without risking undefined behavior.
    const double truncation = (long long)x;
    // If the truncation of x is smaller than x, then it is one less than the
    // desired result.  If it is greater than or equal to x, it is the result.
    // Adding one cannot produce a rounding error because `truncation` is an
    // integer smaller than 2^52.
    const double ceiling = truncation + (truncation < x);
    // Finally, we need to patch up one more thing; the standard specifies that
    // ceil(-small) be -0.0, whereas we will have 0.0 right now.  To handle this
    // correctly, we apply the sign of x to the result.
    return copysign(ceiling, x);
}

像这样的东西就像你能得到的一样优雅，而且仍然是正确的。

我对 Martin 在他的回答中提出的（通常很好！）实现提出了一些担忧。以下是我将如何实施他的方法：

#include <stdint.h>
#include <string.h>

static inline uint64_t toRep(double x) {
    uint64_t r;
    memcpy(&r, &x, sizeof x);
    return r;
}

static inline double fromRep(uint64_t r) {
    double x;
    memcpy(&x, &r, sizeof x);
    return x;
}

double ceil(double x) {

    const uint64_t signbitMask  = UINT64_C(0x8000000000000000);
    const uint64_t significandMask = UINT64_C(0x000fffffffffffff);

    const uint64_t xrep = toRep(x);
    const uint64_t xabs = xrep & signbitMask;

    // If |x| is larger than 2^52 or x is NaN, the result is just x.
    if (xabs >= toRep(0x1.0p52)) return x;

    if (xabs < toRep(1.0)) {
        // If x is in (1.0, 0.0], the result is copysign(0.0, x).
        // We can generate this value by clearing everything except the signbit.
        if (x <= 0.0) return fromRep(xrep & signbitMask);
        // Otherwise x is in (0.0, 1.0), and the result is 1.0.
        else return 1.0;
    }

    // Now we know that the exponent of x is strictly in the range [0, 51],
    // which means that x contains both integral and fractional bits.  We
    // generate a mask covering the fractional bits.
    const int exponent = xabs >> 52;
    const uint64_t fractionalBits = significandMask >> exponent;

    // If x is negative, we want to truncate, so we simply mask off the
    // fractional bits.
    if (xrep & signbitMask) return fromRep(xrep & ~fractionalBits);

    // x is positive; to force rounding to go away from zero, we first *add*
    // the fractionalBits to x, then truncate the result.  The add may
    // overflow the significand into the exponent, but this produces the
    // desired result (zero significand, incremented exponent), so we just
    // let it happen.
    return fromRep(xrep + fractionalBits & ~fractionalBits);
}

关于这种方法需要注意的一点是，它不会为非整数输入引发不精确的浮点标志。这可能会或可能不会影响您的使用。我列出的第一个实现确实提高了标志。

于 2012-09-05T11:06:22.383 回答
0

我不认为宏函数是一个好的解决方案：它不是类型安全的，并且对参数进行了多次评估（副作用）。你应该写一个干净优雅的函数。

于 2012-09-05T11:01:27.927 回答
0

正如我所期望的答案中的更多笑话，我将尝试几个

#define CEILING(X) ceil(X)

奖励：一个没有太多副作用的宏
如果你不太关心负零

#define CEILING(X) (-floor(-(X)))

如果你关心负零，那么

#define CEILING(X) (NEGATIVE_ZERO - floor(-(X)))

将 NEGATIVE_ZERO 的可移植定义留作练习……另外，它还会设置 FP 标志（OVERFLOW INVALID INEXACT）

于 2012-09-06T19:28:54.440 回答