text: /*							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 floor(double);
  int isnan(double);
  int isfinite(double);
  double ldexp(double, int);
  #else
  double floor();
  int isnan(), isfinite();
  double ldexp();
  #endif
  
  double 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 = 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 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 = &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;
      }
  }
  