问用C语言实现ceil()

EN

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

#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);
}

#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);
}