sysdeps/ieee754/dbl-64/s_ceil.c

   1 /* @(#)s_ceil.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 /*
  14  * ceil(x)
  15  * Return x rounded toward -inf to integral value
  16  * Method:
  17  *      Bit twiddling.
  18  * Exception:
  19  *      Inexact flag raised if x not equal to ceil(x).
  20  */
  21
  22 #include <math.h>
  23 #include <math_private.h>
  24
  25 static const double huge = 1.0e300;
  26
  27 double
  28 __ceil (double x)
  29 {
  30   int32_t i0, i1, j0;
  31   u_int32_t i, j;
  32   EXTRACT_WORDS (i0, i1, x);
  33   j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
  34   if (j0 < 20)
  35     {
  36       if (j0 < 0)       /* raise inexact if x != 0 */
  37         {
  38           math_force_eval (huge + x);
  39           /* return 0*sign(x) if |x|<1 */
  40           if (i0 < 0)
  41             {
  42               i0 = 0x80000000; i1 = 0;
  43             }
  44           else if ((i0 | i1) != 0)
  45             {
  46               i0 = 0x3ff00000; i1 = 0;
  47             }
  48         }
  49       else
  50         {
  51           i = (0x000fffff) >> j0;
  52           if (((i0 & i) | i1) == 0)
  53             return x;                        /* x is integral */
  54           math_force_eval (huge + x);           /* raise inexact flag */
  55           if (i0 > 0)
  56             i0 += (0x00100000) >> j0;
  57           i0 &= (~i); i1 = 0;
  58         }
  59     }
  60   else if (j0 > 51)
  61     {
  62       if (j0 == 0x400)
  63         return x + x;                   /* inf or NaN */
  64       else
  65         return x;                       /* x is integral */
  66     }
  67   else
  68     {
  69       i = ((u_int32_t) (0xffffffff)) >> (j0 - 20);
  70       if ((i1 & i) == 0)
  71         return x;                       /* x is integral */
  72       math_force_eval (huge + x);               /* raise inexact flag */
  73       if (i0 > 0)
  74         {
  75           if (j0 == 20)
  76             i0 += 1;
  77           else
  78             {
  79               j = i1 + (1 << (52 - j0));
  80               if (j < i1)
  81                 i0 += 1;                /* got a carry */
  82               i1 = j;
  83             }
  84         }
  85       i1 &= (~i);
  86     }
  87   INSERT_WORDS (x, i0, i1);
  88   return x;
  89 }
  90 #ifndef __ceil
  91 weak_alias (__ceil, ceil)
  92 # ifdef NO_LONG_DOUBLE
  93 strong_alias (__ceil, __ceill)
  94 weak_alias (__ceil, ceill)
  95 # endif
  96 #endif