nsz Git - musl/blob - src/math/floor.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/s_floor.c */
   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  * floor(x)
  14  * Return x rounded toward -inf to integral value
  15  * Method:
  16  *      Bit twiddling.
  17  * Exception:
  18  *      Inexact flag raised if x not equal to floor(x).
  19  */
  20
  21 #include "libm.h"
  22
  23 static const double huge = 1.0e300;
  24
  25 double floor(double x)
  26 {
  27         int32_t i0,i1,j0;
  28         uint32_t i,j;
  29
  30         EXTRACT_WORDS(i0, i1, x);
  31         // FIXME: signed shift
  32         j0 = ((i0>>20)&0x7ff) - 0x3ff;
  33         if (j0 < 20) {
  34                 if (j0 < 0) {  /* |x| < 1 */
  35                         /* raise inexact if x != 0 */
  36                         if (huge+x > 0.0) {
  37                                 if (i0 >= 0) {  /* x >= 0 */
  38                                         i0 = i1 = 0;
  39                                 } else if (((i0&0x7fffffff)|i1) != 0) {
  40                                         i0 = 0xbff00000;
  41                                         i1 = 0;
  42                                 }
  43                         }
  44                 } else {
  45                         i = 0x000fffff>>j0;
  46                         if (((i0&i)|i1) == 0)
  47                                 return x; /* x is integral */
  48                          /* raise inexact flag */
  49                         if (huge+x > 0.0) {
  50                                 if (i0 < 0)
  51                                         i0 += 0x00100000>>j0;
  52                                 i0 &= ~i;
  53                                 i1 = 0;
  54                         }
  55                 }
  56         } else if (j0 > 51) {
  57                 if (j0 == 0x400)
  58                         return x+x; /* inf or NaN */
  59                 else
  60                         return x;   /* x is integral */
  61         } else {
  62                 i = (uint32_t)0xffffffff>>(j0-20);
  63                 if ((i1&i) == 0)
  64                         return x;   /* x is integral */
  65                 /* raise inexact flag */
  66                 if (huge+x > 0.0) {
  67                         if (i0 < 0) {
  68                                 if (j0 == 20)
  69                                         i0++;
  70                                 else {
  71                                         j = i1+(1<<(52-j0));
  72                                         if (j < i1)
  73                                                 i0++; /* got a carry */
  74                                         i1 = j;
  75                                 }
  76                         }
  77                         i1 &= ~i;
  78                 }
  79         }
  80         INSERT_WORDS(x, i0, i1);
  81         return x;
  82 }