Quelle  floor.c

/* ceil()
 *       floor()
 *       frexp()
 *       ldexp()
 *       signbit()
 *       isnan()
 *       isfinite()
 *
 * Floating point numeric utilities
 *
 *
 *
 * SYNOPSIS:
 *
 * double ceil(), floor(), frexp(), ldexp();
 * int signbit(), isnan(), isfinite();
 * double x, y;
 * int expnt, n;
 *
 * y = floor(x);
 * y = ceil(x);
 * y = frexp( x, &expnt );
 * y = ldexp( x, n );
 * n = signbit(x);
 * n = isnan(x);
 * n = isfinite(x);
 *
 *
 *
 * DESCRIPTION:
 *
 * All four routines return a double precision floating point
 * result.
 *
 * floor() returns the largest integer less than or equal to x.
 * It truncates toward minus infinity.
 *
 * ceil() returns the smallest integer greater than or equal
 * to x.  It truncates toward plus infinity.
 *
 * frexp() extracts the exponent from x.  It returns an integer
 * power of two to expnt and the significand between 0.5 and 1
 * to y.  Thus  x = y * 2**expn.
 *
 * ldexp() multiplies x by 2**n.
 *
 * signbit(x) returns 1 if the sign bit of x is 1, else 0.
 *
 * These functions are part of the standard C run time library
 * for many but not all C compilers.  The ones supplied are
 * written in C for either DEC or IEEE arithmetic.  They should
 * be used only if your compiler library does not already have
 * them.
 *
 * The IEEE versions assume that denormal numbers are implemented
 * in the arithmetic.  Some modifications will be required if
 * the arithmetic has abrupt rather than gradual underflow.
 */

/*
Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/

#include "mconf.h"

#ifdef UNK
/* ceil(), floor(), frexp(), ldexp() may need to be rewritten. */
#undef UNK
#if BIGENDIAN
#define MIEEE 1
#else
#define IBMPC 1
#endif
#endif

#ifdef DEC
#define EXPMSK 0x807f
#define MEXP 255
#define NBITS 56
#endif

#ifdef IBMPC
#define EXPMSK 0x800f
#define MEXP 0x7ff
#define NBITS 53
#endif

#ifdef MIEEE
#define EXPMSK 0x800f
#define MEXP 0x7ff
#define NBITS 53
#endif

extern double MAXNUM, NEGZERO;
#ifdef ANSIPROT
double ignore_floor(double);
int isnan(double);
int isfinite(double);
double ldexp(double, int);
#else
double floor();
int isnan(), isfinite();
double ldexp();
#endif

double ignore_ceil(x) double x;
{
  double y;

#ifdef UNK
  mtherr("ceil", DOMAIN);
  return (0.0);
#endif
#ifdef NANS
  if (isnan(x))
    return (x);
#endif
#ifdef INFINITIES
  if (!isfinite(x))
    return (x);
#endif

  y = ignore_floor(x);
  if (y < x)
    y += 1.0;
#ifdef MINUSZERO
  if (y == 0.0 && x < 0.0)
    return (NEGZERO);
#endif
  return (y);
}

/* Bit clearing masks: */

static unsigned short bmask[] = {
    0xffff, 0xfffe, 0xfffc, 0xfff8, 0xfff0, 0xffe0, 0xffc0, 0xff80, 0xff00,
    0xfe00, 0xfc00, 0xf800, 0xf000, 0xe000, 0xc000, 0x8000, 0x0000,
};

double ignore_floor(x) double x;
{
  union {
    double y;
    unsigned short sh[4];
  } u;
  unsigned short *p;
  int e;

#ifdef UNK
  mtherr("floor", DOMAIN);
  return (0.0);
#endif
#ifdef NANS
  if (isnan(x))
    return (x);
#endif
#ifdef INFINITIES
  if (!isfinite(x))
    return (x);
#endif
#ifdef MINUSZERO
  if (x == 0.0L)
    return (x);
#endif
  u.y = x;
/* find the exponent (power of 2) */
#ifdef DEC
  p = (unsigned short *)&u.sh[0];
  e = ((*p >> 7) & 0377) - 0201;
  p += 3;
#endif

#ifdef IBMPC
  p = (unsigned short *)&u.sh[3];
  e = ((*p >> 4) & 0x7ff) - 0x3ff;
  p -= 3;
#endif

#ifdef MIEEE
  p = (unsigned short *)&u.sh[0];
  e = ((*p >> 4) & 0x7ff) - 0x3ff;
  p += 3;
#endif

  if (e < 0) {
    if (u.y < 0.0)
      return (-1.0);
    else
      return (0.0);
  }

  e = (NBITS - 1) - e;
  /* clean out 16 bits at a time */
  while (e >= 16) {
#ifdef IBMPC
    *p++ = 0;
#endif

#ifdef DEC
    *p-- = 0;
#endif

#ifdef MIEEE
    *p-- = 0;
#endif
    e -= 16;
  }

  /* clear the remaining bits */
  if (e > 0)
    *p &= bmask[e];

  if ((x < 0) && (u.y != x))
    u.y -= 1.0;

  return (u.y);
}

double frexp(x, pw2) double x;
int *pw2;
{
  union {
    double y;
    unsigned short sh[4];
  } u;
  int i;
#ifdef DENORMAL
  int k;
#endif
  short *q;

  u.y = x;

#ifdef UNK
  mtherr("frexp", DOMAIN);
  return (0.0);
#endif

#ifdef IBMPC
  q = (short *)&u.sh[3];
#endif

#ifdef DEC
  q = (short *)&u.sh[0];
#endif

#ifdef MIEEE
  q = (short *)&u.sh[0];
#endif

/* find the exponent (power of 2) */
#ifdef DEC
  i = (*q >> 7) & 0377;
  if (i == 0) {
    *pw2 = 0;
    return (0.0);
  }
  i -= 0200;
  *pw2 = i;
  *q &= 0x807f; /* strip all exponent bits */
  *q |= 040000; /* mantissa between 0.5 and 1 */
  return (u.y);
#endif

#ifdef IBMPC
  i = (*q >> 4) & 0x7ff;
  if (i != 0)
    goto ieeedon;
#endif

#ifdef MIEEE
  i = *q >> 4;
  i &= 0x7ff;
  if (i != 0)
    goto ieeedon;
#ifdef DENORMAL

#else
  *pw2 = 0;
  return (0.0);
#endif

#endif

#ifndef DEC
/* Number is denormal or zero */
#ifdef DENORMAL
  if (u.y == 0.0) {
    *pw2 = 0;
    return (0.0);
  }

  /* Handle denormal number. */
  do {
    u.y *= 2.0;
    i -= 1;
    k = (*q >> 4) & 0x7ff;
  } while (k == 0);
  i = i + k;
#endif /* DENORMAL */

ieeedon:

  i -= 0x3fe;
  *pw2 = i;
  *q &= 0x800f;
  *q |= 0x3fe0;
  return (u.y);
#endif
}

double ldexp(x, pw2) double x;
int pw2;
{
  union {
    double y;
    unsigned short sh[4];
  } u;
  short *q;
  int e;

#ifdef UNK
  mtherr("ldexp", DOMAIN);
  return (0.0);
#endif

  u.y = x;
#ifdef DEC
  q = (short *)&u.sh[0];
  e = (*q >> 7) & 0377;
  if (e == 0)
    return (0.0);
#else

#ifdef IBMPC
  q = (short *)&u.sh[3];
#endif
#ifdef MIEEE
  q = (short *)&u.sh[0];
#endif
  while ((e = (*q & 0x7ff0) >> 4) == 0) {
    if (u.y == 0.0) {
      return (0.0);
    }
    /* Input is denormal. */
    if (pw2 > 0) {
      u.y *= 2.0;
      pw2 -= 1;
    }
    if (pw2 < 0) {
      if (pw2 < -53)
        return (0.0);
      u.y /= 2.0;
      pw2 += 1;
    }
    if (pw2 == 0)
      return (u.y);
  }
#endif /* not DEC */

  e += pw2;

/* Handle overflow */
#ifdef DEC
  if (e > MEXP)
    return (MAXNUM);
#else
  if (e >= MEXP)
    return (2.0 * MAXNUM);
#endif

  /* Handle denormalized results */
  if (e < 1) {
#ifdef DENORMAL
    if (e < -53)
      return (0.0);
    *q &= 0x800f;
    *q |= 0x10;
    /* For denormals, significant bits may be lost even
       when dividing by 2.  Construct 2^-(1-e) so the result
       is obtained with only one multiplication.  */
    u.y *= ldexp(1.0, e - 1);
    return (u.y);
#else
    return (0.0);
#endif
  } else {
#ifdef DEC
    *q &= 0x807f; /* strip all exponent bits */
    *q |= (e & 0xff) << 7;
#else
    *q &= 0x800f;
    *q |= (e & 0x7ff) << 4;
#endif
    return (u.y);
  }
}