lib/libm/src/s_floorl.c

   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_floor.c 5.1 93/09/24
  12  */
  13
  14 #include <sys/cdefs.h>
  15 __RCSID("$NetBSD: s_floorl.c,v 1.1 2013/11/11 23:57:34 joerg Exp $");
  16 #if 0
  17 __FBSDID("$FreeBSD: head/lib/msun/src/s_floorl.c 176280 2008-02-14 15:10:34Z bde $");
  18 #endif
  19
  20 /*
  21  * floorl(x)
  22  * Return x rounded toward -inf to integral value
  23  * Method:
  24  *      Bit twiddling.
  25  * Exception:
  26  *      Inexact flag raised if x not equal to floorl(x).
  27  */
  28 #include "namespace.h"
  29
  30 #include <float.h>
  31 #include <math.h>
  32 #include <stdint.h>
  33 #include <machine/ieee.h>
  34
  35 #ifdef __HAVE_LONG_DOUBLE
  36
  37 #ifdef __weak_alias
  38 __weak_alias(floorl, _floorl)
  39 #endif
  40
  41 #ifdef LDBL_IMPLICIT_NBIT
  42 #define MANH_SIZE       (EXT_FRACHBITS + 1)
  43 #define INC_MANH(ux, c) do {                                    \
  44         uint64_t oi = ux.extu_frach;                            \
  45         ux.extu_frach += (c);                                   \
  46         if (ux.extu_frach < oi)                                 \
  47                 ux.extu_exp++;                                  \
  48 } while (0)
  49 #else
  50 #define MANH_SIZE       EXT_FRACHBITS
  51 #define INC_MANH(ux, c) do {                                    \
  52         uint64_t oi = ux.extu_frach;                            \
  53         ux.extu_frach += (c);                                   \
  54         if (ux.extu_frach < oi) {                               \
  55                 ux.extu_exp++;                                  \
  56                 ux.extu_frach |= 1llu << (EXT_FRACHBITS - 1);   \
  57         }                                                       \
  58 } while (0)
  59 #endif
  60
  61 static const long double huge = 1.0e300;
  62
  63 long double
  64 floorl(long double x)
  65 {
  66         union ieee_ext_u ux = { .extu_ld = x, };
  67         int e = ux.extu_exp - LDBL_MAX_EXP + 1;
  68
  69         if (e < MANH_SIZE - 1) {
  70                 if (e < 0) {                    /* raise inexact if x != 0 */
  71                         if (huge + x > 0.0)
  72                                 if (ux.extu_exp > 0 ||
  73                                     (ux.extu_frach | ux.extu_fracl) != 0)
  74                                         ux.extu_ld = ux.extu_sign ? -1.0 : 0.0;
  75                 } else {
  76                         uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
  77                         if (((ux.extu_frach & m) | ux.extu_fracl) == 0)
  78                                 return (x);     /* x is integral */
  79                         if (ux.extu_sign) {
  80 #ifdef LDBL_IMPLICIT_NBIT
  81                                 if (e == 0)
  82                                         ux.extu_exp++;
  83                                 else
  84 #endif
  85                                 INC_MANH(ux, 1llu << (MANH_SIZE - e - 1));
  86                         }
  87                         if (huge + x > 0.0) {   /* raise inexact flag */
  88                                 ux.extu_frach &= ~m;
  89                                 ux.extu_fracl = 0;
  90                         }
  91                 }
  92         } else if (e < LDBL_MANT_DIG - 1) {
  93                 uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
  94                 if ((ux.extu_fracl & m) == 0)
  95                         return (x);     /* x is integral */
  96                 if (ux.extu_sign) {
  97                         if (e == MANH_SIZE - 1)
  98                                 INC_MANH(ux, 1);
  99                         else {
 100                                 uint64_t o = ux.extu_fracl;
 101                                 ux.extu_fracl += 1llu << (LDBL_MANT_DIG - e - 1);
 102                                 if (ux.extu_fracl < o)  /* got a carry */
 103                                         INC_MANH(ux, 1);
 104                         }
 105                 }
 106                 if (huge + x > 0.0)             /* raise inexact flag */
 107                         ux.extu_fracl &= ~m;
 108         }
 109         return (ux.extu_ld);
 110 }
 111
 112 #endif /* __HAVE_LONG_DOUBLE */