Branch data     Line data    Source code

       1                 :             : /*
       2                 :             :  * ====================================================
       3                 :             :  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
       4                 :             :  *
       5                 :             :  * Developed at SunPro, a Sun Microsystems, Inc. business.
       6                 :             :  * Permission to use, copy, modify, and distribute this
       7                 :             :  * software is freely granted, provided that this notice
       8                 :             :  * is preserved.
       9                 :             :  * ====================================================
      10                 :             :  *
      11                 :             :  * From: @(#)s_ceil.c 5.1 93/09/24
      12                 :             :  */
      13                 :             : 
      14                 :             : #include "cdefs-compat.h"
      15                 :             : //__FBSDID("$FreeBSD: src/lib/msun/src/s_ceill.c,v 1.9 2008/02/14 15:10:33 bde Exp $");
      16                 :             : 
      17                 :             : /*
      18                 :             :  * ceill(x)
      19                 :             :  * Return x rounded toward -inf to integral value
      20                 :             :  * Method:
      21                 :             :  *      Bit twiddling.
      22                 :             :  * Exception:
      23                 :             :  *      Inexact flag raised if x not equal to ceill(x).
      24                 :             :  */
      25                 :             : 
      26                 :             : #include <float.h>
      27                 :             : #include <openlibm_math.h>
      28                 :             : #include <stdint.h>
      29                 :             : 
      30                 :             : #include "fpmath.h"
      31                 :             : #include "math_private.h"
      32                 :             : 
      33                 :             : #ifdef LDBL_IMPLICIT_NBIT
      34                 :             : #define MANH_SIZE       (LDBL_MANH_SIZE + 1)
      35                 :             : #define INC_MANH(u, c)  do {                                    \
      36                 :             :         u_int64_t o = u.bits.manh;                              \
      37                 :             :         u.bits.manh += (c);                                     \
      38                 :             :         if (u.bits.manh < o)                                 \
      39                 :             :                 u.bits.exp++;                                   \
      40                 :             : } while (0)
      41                 :             : #else
      42                 :             : #define MANH_SIZE       LDBL_MANH_SIZE
      43                 :             : #define INC_MANH(u, c)  do {                                    \
      44                 :             :         u_int64_t o = u.bits.manh;                              \
      45                 :             :         u.bits.manh += (c);                                     \
      46                 :             :         if (u.bits.manh < o) {                                       \
      47                 :             :                 u.bits.exp++;                                   \
      48                 :             :                 u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1);      \
      49                 :             :         }                                                       \
      50                 :             : } while (0)
      51                 :             : #endif
      52                 :             : 
      53                 :             : static const long double huge = 1.0e300;
      54                 :             : 
      55                 :             : OLM_DLLEXPORT long double
      56                 :           0 : ceill(long double x)
      57                 :             : {
      58                 :           0 :         union IEEEl2bits u = { .e = x };
      59                 :           0 :         int e = u.bits.exp - LDBL_MAX_EXP + 1;
      60                 :             : 
      61         [ #  # ]:           0 :         if (e < MANH_SIZE - 1) {
      62         [ #  # ]:           0 :                 if (e < 0) {                 /* raise inexact if x != 0 */
      63         [ #  # ]:           0 :                         if (huge + x > 0.0)
      64         [ #  # ]:           0 :                                 if (u.bits.exp > 0 ||
      65         [ #  # ]:           0 :                                     (u.bits.manh | u.bits.manl) != 0)
      66         [ #  # ]:           0 :                                         u.e = u.bits.sign ? -0.0 : 1.0;
      67                 :             :                 } else {
      68                 :           0 :                         u_int64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
      69         [ #  # ]:           0 :                         if (((u.bits.manh & m) | u.bits.manl) == 0)
      70                 :           0 :                                 return (x);     /* x is integral */
      71         [ #  # ]:           0 :                         if (!u.bits.sign) {
      72                 :             : #ifdef LDBL_IMPLICIT_NBIT
      73                 :             :                                 if (e == 0)
      74                 :             :                                         u.bits.exp++;
      75                 :             :                                 else
      76                 :             : #endif
      77         [ #  # ]:           0 :                                 INC_MANH(u, 1llu << (MANH_SIZE - e - 1));
      78                 :             :                         }
      79         [ #  # ]:           0 :                         if (huge + x > 0.0) {        /* raise inexact flag */
      80                 :           0 :                                 u.bits.manh &= ~m;
      81                 :           0 :                                 u.bits.manl = 0;
      82                 :             :                         }
      83                 :             :                 }
      84         [ #  # ]:           0 :         } else if (e < LDBL_MANT_DIG - 1) {
      85                 :           0 :                 u_int64_t m = (u_int64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
      86         [ #  # ]:           0 :                 if ((u.bits.manl & m) == 0)
      87                 :           0 :                         return (x);     /* x is integral */
      88         [ #  # ]:           0 :                 if (!u.bits.sign) {
      89         [ #  # ]:           0 :                         if (e == MANH_SIZE - 1)
      90         [ #  # ]:           0 :                                 INC_MANH(u, 1);
      91                 :             :                         else {
      92                 :           0 :                                 u_int64_t o = u.bits.manl;
      93                 :           0 :                                 u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1);
      94         [ #  # ]:           0 :                                 if (u.bits.manl < o) /* got a carry */
      95         [ #  # ]:           0 :                                         INC_MANH(u, 1);
      96                 :             :                         }
      97                 :             :                 }
      98         [ #  # ]:           0 :                 if (huge + x > 0.0)          /* raise inexact flag */
      99                 :           0 :                         u.bits.manl &= ~m;
     100                 :             :         }
     101                 :           0 :         return (u.e);
     102                 :             : }