src/libm/s_floor.c

   1 /* @(#)s_floor.c 5.1 93/09/24 */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunPro, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12
  13 #if defined(LIBM_SCCS) && !defined(lint)
  14 static const char rcsid[] =
  15     "$NetBSD: s_floor.c,v 1.8 1995/05/10 20:47:20 jtc Exp $";
  16 #endif
  17
  18 /*
  19  * floor(x)
  20  * Return x rounded toward -inf to integral value
  21  * Method:
  22  *      Bit twiddling.
  23  * Exception:
  24  *      Inexact flag raised if x not equal to floor(x).
  25  */
  26
  27 #include "math_libm.h"
  28 #include "math_private.h"
  29
  30 #ifdef __STDC__
  31 static const double huge_val = 1.0e300;
  32 #else
  33 static double huge_val = 1.0e300;
  34 #endif
  35
  36 libm_hidden_proto(floor)
  37 #ifdef __STDC__
  38      double floor(double x)
  39 #else
  40      double floor(x)
  41      double x;
  42 #endif
  43 {
  44     int32_t i0, i1, j0;
  45     u_int32_t i, j;
  46     EXTRACT_WORDS(i0, i1, x);
  47     j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
  48     if (j0 < 20) {
  49         if (j0 < 0) {           /* raise inexact if x != 0 */
  50             if (huge_val + x > 0.0) {       /* return 0*sign(x) if |x|<1 */
  51                 if (i0 >= 0) {
  52                     i0 = i1 = 0;
  53                 } else if (((i0 & 0x7fffffff) | i1) != 0) {
  54                     i0 = 0xbff00000;
  55                     i1 = 0;
  56                 }
  57             }
  58         } else {
  59             i = (0x000fffff) >> j0;
  60             if (((i0 & i) | i1) == 0)
  61                 return x;       /* x is integral */
  62             if (huge_val + x > 0.0) {       /* raise inexact flag */
  63                 if (i0 < 0)
  64                     i0 += (0x00100000) >> j0;
  65                 i0 &= (~i);
  66                 i1 = 0;
  67             }
  68         }
  69     } else if (j0 > 51) {
  70         if (j0 == 0x400)
  71             return x + x;       /* inf or NaN */
  72         else
  73             return x;           /* x is integral */
  74     } else {
  75         i = ((u_int32_t) (0xffffffff)) >> (j0 - 20);
  76         if ((i1 & i) == 0)
  77             return x;           /* x is integral */
  78         if (huge_val + x > 0.0) {   /* raise inexact flag */
  79             if (i0 < 0) {
  80                 if (j0 == 20)
  81                     i0 += 1;
  82                 else {
  83                     j = i1 + (1 << (52 - j0));
  84                     if (j < (u_int32_t) i1)
  85                         i0 += 1;        /* got a carry */
  86                     i1 = j;
  87                 }
  88             }
  89             i1 &= (~i);
  90         }
  91     }
  92     INSERT_WORDS(x, i0, i1);
  93     return x;
  94 }
  95
  96 libm_hidden_def(floor)