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                                 /* return 0*sign(x) if |x|<1 */
  38                                 if (i0 < 0) {
  39                                         i0 = 0x80000000;
  40                                         i1=0;
  41                                 } else if ((i0|i1) != 0) {
  42                                         i0=0x3ff00000;
  43                                         i1=0;
  44                                 }
  45                         }
  46                 } else {
  47                         i = (0x000fffff)>>j0;
  48                         if (((i0&i)|i1) == 0) /* x is integral */
  49                                 return x;
  50                         /* raise inexact flag */
  51                         if (huge+x > 0.0) {
  52                                 if (i0 > 0)
  53                                         i0 += 0x00100000>>j0;
  54                                 i0 &= ~i;
  55                                 i1 = 0;
  56                         }
  57                 }
  58         } else if (j0 > 51) {
  59                 if (j0 == 0x400)  /* inf or NaN */
  60                         return x+x;
  61                 return x;         /* x is integral */
  62         } else {
  63                 i = (uint32_t)0xffffffff>>(j0-20);
  64                 if ((i1&i) == 0)
  65                         return x; /* x is integral */
  66                 /* raise inexact flag */
  67                 if (huge+x > 0.0) {
  68                         if (i0 > 0) {
  69                                 if (j0 == 20)
  70                                         i0 += 1;
  71                                 else {
  72                                         j = i1 + (1<<(52-j0));
  73                                         if (j < i1)  /* got a carry */
  74                                                 i0 += 1;
  75                                         i1 = j;
  76                                 }
  77                         }
  78                         i1 &= ~i;
  79                 }
  80         }
  81         INSERT_WORDS(x, i0, i1);
  82         return x;
  83 }