lib/libm/floor.c

   1 /*                                                      ceil()
   2  *                                                      floor()
   3  *                                                      frexp()
   4  *                                                      ldexp()
   5  *                                                      signbit()
   6  *                                                      isnan()
   7  *                                                      isfinite()
   8  *
   9  *      Floating point numeric utilities
  10  *
  11  *
  12  *
  13  * SYNOPSIS:
  14  *
  15  * double ceil(), floor(), frexp(), ldexp();
  16  * int signbit(), isnan(), isfinite();
  17  * double x, y;
  18  * int expnt, n;
  19  *
  20  * y = floor(x);
  21  * y = ceil(x);
  22  * y = frexp( x, &expnt );
  23  * y = ldexp( x, n );
  24  * n = signbit(x);
  25  * n = isnan(x);
  26  * n = isfinite(x);
  27  *
  28  *
  29  *
  30  * DESCRIPTION:
  31  *
  32  * All four routines return a double precision floating point
  33  * result.
  34  *
  35  * floor() returns the largest integer less than or equal to x.
  36  * It truncates toward minus infinity.
  37  *
  38  * ceil() returns the smallest integer greater than or equal
  39  * to x.  It truncates toward plus infinity.
  40  *
  41  * frexp() extracts the exponent from x.  It returns an integer
  42  * power of two to expnt and the significand between 0.5 and 1
  43  * to y.  Thus  x = y * 2**expn.
  44  *
  45  * ldexp() multiplies x by 2**n.
  46  *
  47  * signbit(x) returns 1 if the sign bit of x is 1, else 0.
  48  *
  49  * These functions are part of the standard C run time library
  50  * for many but not all C compilers.  The ones supplied are
  51  * written in C for either DEC or IEEE arithmetic.  They should
  52  * be used only if your compiler library does not already have
  53  * them.
  54  *
  55  * The IEEE versions assume that denormal numbers are implemented
  56  * in the arithmetic.  Some modifications will be required if
  57  * the arithmetic has abrupt rather than gradual underflow.
  58  */
  59 \f
  60
  61 /*
  62 Cephes Math Library Release 2.8:  June, 2000
  63 Copyright 1984, 1995, 2000 by Stephen L. Moshier
  64 */
  65
  66
  67 #include "mconf.h"
  68
  69 #ifdef UNK
  70 /* ceil(), floor(), frexp(), ldexp() may need to be rewritten. */
  71 #undef UNK
  72 #if BIGENDIAN
  73 #define MIEEE 1
  74 #else
  75 #define IBMPC 1
  76 #endif
  77 #endif
  78
  79 #ifdef DEC
  80 #define EXPMSK 0x807f
  81 #define MEXP 255
  82 #define NBITS 56
  83 #endif
  84
  85 #ifdef IBMPC
  86 #define EXPMSK 0x800f
  87 #define MEXP 0x7ff
  88 #define NBITS 53
  89 #endif
  90
  91 #ifdef MIEEE
  92 #define EXPMSK 0x800f
  93 #define MEXP 0x7ff
  94 #define NBITS 53
  95 #endif
  96
  97 extern double MAXNUM, NEGZERO;
  98 #ifdef ANSIPROT
  99 double floor ( double );
 100 int isnan ( double );
 101 int isfinite ( double );
 102 double ldexp ( double, int );
 103 #else
 104 double floor();
 105 int isnan(), isfinite();
 106 double ldexp();
 107 #endif
 108
 109 double ceil(x)
 110 double x;
 111 {
 112 double y;
 113
 114 #ifdef UNK
 115 mtherr( "ceil", DOMAIN );
 116 return(0.0);
 117 #endif
 118 #ifdef NANS
 119 if( isnan(x) )
 120         return( x );
 121 #endif
 122 #ifdef INFINITIES
 123 if(!isfinite(x))
 124         return(x);
 125 #endif
 126
 127 y = floor(x);
 128 if( y < x )
 129         y += 1.0;
 130 #ifdef MINUSZERO
 131 if( y == 0.0 && x < 0.0 )
 132         return( NEGZERO );
 133 #endif
 134 return(y);
 135 }
 136
 137
 138
 139
 140 /* Bit clearing masks: */
 141
 142 static unsigned short bmask[] = {
 143 0xffff,
 144 0xfffe,
 145 0xfffc,
 146 0xfff8,
 147 0xfff0,
 148 0xffe0,
 149 0xffc0,
 150 0xff80,
 151 0xff00,
 152 0xfe00,
 153 0xfc00,
 154 0xf800,
 155 0xf000,
 156 0xe000,
 157 0xc000,
 158 0x8000,
 159 0x0000,
 160 };
 161
 162
 163
 164
 165
 166 double floor(x)
 167 double x;
 168 {
 169 union
 170         {
 171         double y;
 172         unsigned short sh[4];
 173         } u;
 174 unsigned short *p;
 175 int e;
 176
 177 #ifdef UNK
 178 mtherr( "floor", DOMAIN );
 179 return(0.0);
 180 #endif
 181 #ifdef NANS
 182 if( isnan(x) )
 183         return( x );
 184 #endif
 185 #ifdef INFINITIES
 186 if(!isfinite(x))
 187         return(x);
 188 #endif
 189 #ifdef MINUSZERO
 190 if(x == 0.0)
 191         return(x);
 192 #endif
 193 u.y = x;
 194 /* find the exponent (power of 2) */
 195 #ifdef DEC
 196 p = (unsigned short *)&u.sh[0];
 197 e = (( *p  >> 7) & 0377) - 0201;
 198 p += 3;
 199 #endif
 200
 201 #ifdef IBMPC
 202 p = (unsigned short *)&u.sh[3];
 203 e = (( *p >> 4) & 0x7ff) - 0x3ff;
 204 p -= 3;
 205 #endif
 206
 207 #ifdef MIEEE
 208 p = (unsigned short *)&u.sh[0];
 209 e = (( *p >> 4) & 0x7ff) - 0x3ff;
 210 p += 3;
 211 #endif
 212
 213 if( e < 0 )
 214         {
 215         if( u.y < 0.0 )
 216                 return( -1.0 );
 217         else
 218                 return( 0.0 );
 219         }
 220
 221 e = (NBITS -1) - e;
 222 /* clean out 16 bits at a time */
 223 while( e >= 16 )
 224         {
 225 #ifdef IBMPC
 226         *p++ = 0;
 227 #endif
 228
 229 #ifdef DEC
 230         *p-- = 0;
 231 #endif
 232
 233 #ifdef MIEEE
 234         *p-- = 0;
 235 #endif
 236         e -= 16;
 237         }
 238
 239 /* clear the remaining bits */
 240 if( e > 0 )
 241         *p &= bmask[e];
 242
 243 if( (x < 0) && (u.y != x) )
 244         u.y -= 1.0;
 245
 246 return(u.y);
 247 }
 248
 249
 250
 251
 252 double frexp( x, pw2 )
 253 double x;
 254 int *pw2;
 255 {
 256 union
 257         {
 258         double y;
 259         unsigned short sh[4];
 260         } u;
 261 int i;
 262 #ifdef DENORMAL
 263 int k;
 264 #endif
 265 short *q;
 266
 267 u.y = x;
 268
 269 #ifdef UNK
 270 mtherr( "frexp", DOMAIN );
 271 return(0.0);
 272 #endif
 273
 274 #ifdef IBMPC
 275 q = (short *)&u.sh[3];
 276 #endif
 277
 278 #ifdef DEC
 279 q = (short *)&u.sh[0];
 280 #endif
 281
 282 #ifdef MIEEE
 283 q = (short *)&u.sh[0];
 284 #endif
 285
 286 /* find the exponent (power of 2) */
 287 #ifdef DEC
 288 i  = ( *q >> 7) & 0377;
 289 if( i == 0 )
 290         {
 291         *pw2 = 0;
 292         return(0.0);
 293         }
 294 i -= 0200;
 295 *pw2 = i;
 296 *q &= 0x807f;   /* strip all exponent bits */
 297 *q |= 040000;   /* mantissa between 0.5 and 1 */
 298 return(u.y);
 299 #endif
 300
 301 #ifdef IBMPC
 302 i  = ( *q >> 4) & 0x7ff;
 303 if( i != 0 )
 304         goto ieeedon;
 305 #endif
 306
 307 #ifdef MIEEE
 308 i  =  *q >> 4;
 309 i &= 0x7ff;
 310 if( i != 0 )
 311         goto ieeedon;
 312 #ifdef DENORMAL
 313
 314 #else
 315 *pw2 = 0;
 316 return(0.0);
 317 #endif
 318
 319 #endif
 320
 321
 322 #ifndef DEC
 323 /* Number is denormal or zero */
 324 #ifdef DENORMAL
 325 if( u.y == 0.0 )
 326         {
 327         *pw2 = 0;
 328         return( 0.0 );
 329         }
 330
 331
 332 /* Handle denormal number. */
 333 do
 334         {
 335         u.y *= 2.0;
 336         i -= 1;
 337         k  = ( *q >> 4) & 0x7ff;
 338         }
 339 while( k == 0 );
 340 i = i + k;
 341 #endif /* DENORMAL */
 342
 343 ieeedon:
 344
 345 i -= 0x3fe;
 346 *pw2 = i;
 347 *q &= 0x800f;
 348 *q |= 0x3fe0;
 349 return( u.y );
 350 #endif
 351 }
 352
 353
 354
 355
 356
 357
 358
 359 double ldexp( x, pw2 )
 360 double x;
 361 int pw2;
 362 {
 363 union
 364         {
 365         double y;
 366         unsigned short sh[4];
 367         } u;
 368 short *q;
 369 int e;
 370
 371 #ifdef UNK
 372 mtherr( "ldexp", DOMAIN );
 373 return(0.0);
 374 #endif
 375
 376 u.y = x;
 377 #ifdef DEC
 378 q = (short *)&u.sh[0];
 379 e  = ( *q >> 7) & 0377;
 380 if( e == 0 )
 381         return(0.0);
 382 #else
 383
 384 #ifdef IBMPC
 385 q = (short *)&u.sh[3];
 386 #endif
 387 #ifdef MIEEE
 388 q = (short *)&u.sh[0];
 389 #endif
 390 while( (e = (*q & 0x7ff0) >> 4) == 0 )
 391         {
 392         if( u.y == 0.0 )
 393                 {
 394                 return( 0.0 );
 395                 }
 396 /* Input is denormal. */
 397         if( pw2 > 0 )
 398                 {
 399                 u.y *= 2.0;
 400                 pw2 -= 1;
 401                 }
 402         if( pw2 < 0 )
 403                 {
 404                 if( pw2 < -53 )
 405                         return(0.0);
 406                 u.y /= 2.0;
 407                 pw2 += 1;
 408                 }
 409         if( pw2 == 0 )
 410                 return(u.y);
 411         }
 412 #endif /* not DEC */
 413
 414 e += pw2;
 415
 416 /* Handle overflow */
 417 #ifdef DEC
 418 if( e > MEXP )
 419         return( MAXNUM );
 420 #else
 421 if( e >= MEXP )
 422         return( 2.0*MAXNUM );
 423 #endif
 424
 425 /* Handle denormalized results */
 426 if( e < 1 )
 427         {
 428 #ifdef DENORMAL
 429         if( e < -53 )
 430                 return(0.0);
 431         *q &= 0x800f;
 432         *q |= 0x10;
 433         /* For denormals, significant bits may be lost even
 434            when dividing by 2.  Construct 2^-(1-e) so the result
 435            is obtained with only one multiplication.  */
 436         u.y *= ldexp(1.0, e-1);
 437         return(u.y);
 438 #else
 439         return(0.0);
 440 #endif
 441         }
 442 else
 443         {
 444 #ifdef DEC
 445         *q &= 0x807f;   /* strip all exponent bits */
 446         *q |= (e & 0xff) << 7;
 447 #else
 448         *q &= 0x800f;
 449         *q |= (e & 0x7ff) << 4;
 450 #endif
 451         return(u.y);
 452         }
 453 }