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

   1 /* origin: FreeBSD /usr/src/lib/msun/src/s_ceil.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  * ceil(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 ceil(x).
  19  */
  20
  21 #include "libm.h"
  22
  23 static const double huge = 1.0e300;
  24
  25 double ceil(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) {
  35                          /* raise inexact if x != 0 */
  36                         if (huge+x > 0.0) {
  37                                 if (i0 < 0) {
  38                                         i0 = 0x80000000;
  39                                         i1=0;
  40                                 } else if ((i0|i1) != 0) {
  41                                         i0=0x3ff00000;
  42                                         i1=0;
  43                                 }
  44                         }
  45                 } else {
  46                         i = 0x000fffff>>j0;
  47                         if (((i0&i)|i1) == 0) /* x is integral */
  48                                 return x;
  49                         /* raise inexact flag */
  50                         if (huge+x > 0.0) {
  51                                 if (i0 > 0)
  52                                         i0 += 0x00100000>>j0;
  53                                 i0 &= ~i;
  54                                 i1 = 0;
  55                         }
  56                 }
  57         } else if (j0 > 51) {
  58                 if (j0 == 0x400)  /* inf or NaN */
  59                         return x+x;
  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 += 1;
  70                                 else {
  71                                         j = i1 + (1<<(52-j0));
  72                                         if (j < i1)  /* got a carry */
  73                                                 i0 += 1;
  74                                         i1 = j;
  75                                 }
  76                         }
  77                         i1 &= ~i;
  78                 }
  79         }
  80         INSERT_WORDS(x, i0, i1);
  81         return x;
  82 }