Yet more math...
authorSam Lantinga <slouken@libsdl.org>
Mon, 15 Sep 2008 06:48:41 +0000
changeset 2758045d9976f285
parent 2757 0581f49c9294
child 2759 95fccd9bf262
Yet more math...
src/libm/s_floor.c
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/libm/s_floor.c	Mon Sep 15 06:48:41 2008 +0000
     1.3 @@ -0,0 +1,96 @@
     1.4 +/* @(#)s_floor.c 5.1 93/09/24 */
     1.5 +/*
     1.6 + * ====================================================
     1.7 + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     1.8 + *
     1.9 + * Developed at SunPro, a Sun Microsystems, Inc. business.
    1.10 + * Permission to use, copy, modify, and distribute this
    1.11 + * software is freely granted, provided that this notice
    1.12 + * is preserved.
    1.13 + * ====================================================
    1.14 + */
    1.15 +
    1.16 +#if defined(LIBM_SCCS) && !defined(lint)
    1.17 +static char rcsid[] =
    1.18 +    "$NetBSD: s_floor.c,v 1.8 1995/05/10 20:47:20 jtc Exp $";
    1.19 +#endif
    1.20 +
    1.21 +/*
    1.22 + * floor(x)
    1.23 + * Return x rounded toward -inf to integral value
    1.24 + * Method:
    1.25 + *	Bit twiddling.
    1.26 + * Exception:
    1.27 + *	Inexact flag raised if x not equal to floor(x).
    1.28 + */
    1.29 +
    1.30 +#include "math.h"
    1.31 +#include "math_private.h"
    1.32 +
    1.33 +#ifdef __STDC__
    1.34 +static const double huge = 1.0e300;
    1.35 +#else
    1.36 +static double huge = 1.0e300;
    1.37 +#endif
    1.38 +
    1.39 +libm_hidden_proto(floor)
    1.40 +#ifdef __STDC__
    1.41 +     double floor(double x)
    1.42 +#else
    1.43 +     double floor(x)
    1.44 +     double x;
    1.45 +#endif
    1.46 +{
    1.47 +    int32_t i0, i1, j0;
    1.48 +    u_int32_t i, j;
    1.49 +    EXTRACT_WORDS(i0, i1, x);
    1.50 +    j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
    1.51 +    if (j0 < 20) {
    1.52 +        if (j0 < 0) {           /* raise inexact if x != 0 */
    1.53 +            if (huge + x > 0.0) {       /* return 0*sign(x) if |x|<1 */
    1.54 +                if (i0 >= 0) {
    1.55 +                    i0 = i1 = 0;
    1.56 +                } else if (((i0 & 0x7fffffff) | i1) != 0) {
    1.57 +                    i0 = 0xbff00000;
    1.58 +                    i1 = 0;
    1.59 +                }
    1.60 +            }
    1.61 +        } else {
    1.62 +            i = (0x000fffff) >> j0;
    1.63 +            if (((i0 & i) | i1) == 0)
    1.64 +                return x;       /* x is integral */
    1.65 +            if (huge + x > 0.0) {       /* raise inexact flag */
    1.66 +                if (i0 < 0)
    1.67 +                    i0 += (0x00100000) >> j0;
    1.68 +                i0 &= (~i);
    1.69 +                i1 = 0;
    1.70 +            }
    1.71 +        }
    1.72 +    } else if (j0 > 51) {
    1.73 +        if (j0 == 0x400)
    1.74 +            return x + x;       /* inf or NaN */
    1.75 +        else
    1.76 +            return x;           /* x is integral */
    1.77 +    } else {
    1.78 +        i = ((u_int32_t) (0xffffffff)) >> (j0 - 20);
    1.79 +        if ((i1 & i) == 0)
    1.80 +            return x;           /* x is integral */
    1.81 +        if (huge + x > 0.0) {   /* raise inexact flag */
    1.82 +            if (i0 < 0) {
    1.83 +                if (j0 == 20)
    1.84 +                    i0 += 1;
    1.85 +                else {
    1.86 +                    j = i1 + (1 << (52 - j0));
    1.87 +                    if (j < i1)
    1.88 +                        i0 += 1;        /* got a carry */
    1.89 +                    i1 = j;
    1.90 +                }
    1.91 +            }
    1.92 +            i1 &= (~i);
    1.93 +        }
    1.94 +    }
    1.95 +    INSERT_WORDS(x, i0, i1);
    1.96 +    return x;
    1.97 +}
    1.98 +
    1.99 +libm_hidden_def(floor)