sysdeps/ieee754/ldbl-128/s_roundevenl.c

   1 /* Round to nearest integer value, rounding halfway cases to even.
   2    ldbl-128 version.
   3    Copyright (C) 2016-2018 Free Software Foundation, Inc.
   4    This file is part of the GNU C Library.
   5
   6    The GNU C Library is free software; you can redistribute it and/or
   7    modify it under the terms of the GNU Lesser General Public
   8    License as published by the Free Software Foundation; either
   9    version 2.1 of the License, or (at your option) any later version.
  10
  11    The GNU C Library is distributed in the hope that it will be useful,
  12    but WITHOUT ANY WARRANTY; without even the implied warranty of
  13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14    Lesser General Public License for more details.
  15
  16    You should have received a copy of the GNU Lesser General Public
  17    License along with the GNU C Library; if not, see
  18    <http://www.gnu.org/licenses/>.  */
  19
  20 #include <math.h>
  21 #include <math_private.h>
  22 #include <libm-alias-ldouble.h>
  23 #include <stdint.h>
  24
  25 #define BIAS 0x3fff
  26 #define MANT_DIG 113
  27 #define MAX_EXP (2 * BIAS + 1)
  28
  29 _Float128
  30 __roundevenl (_Float128 x)
  31 {
  32   uint64_t hx, lx, uhx;
  33   GET_LDOUBLE_WORDS64 (hx, lx, x);
  34   uhx = hx & 0x7fffffffffffffffULL;
  35   int exponent = uhx >> (MANT_DIG - 1 - 64);
  36   if (exponent >= BIAS + MANT_DIG - 1)
  37     {
  38       /* Integer, infinity or NaN.  */
  39       if (exponent == MAX_EXP)
  40         /* Infinity or NaN; quiet signaling NaNs.  */
  41         return x + x;
  42       else
  43         return x;
  44     }
  45   else if (exponent >= BIAS + MANT_DIG - 64)
  46     {
  47       /* Not necessarily an integer; integer bit is in low word.
  48          Locate the bits with exponents 0 and -1.  */
  49       int int_pos = (BIAS + MANT_DIG - 1) - exponent;
  50       int half_pos = int_pos - 1;
  51       uint64_t half_bit = 1ULL << half_pos;
  52       uint64_t int_bit = 1ULL << int_pos;
  53       if ((lx & (int_bit | (half_bit - 1))) != 0)
  54         {
  55           /* Carry into the exponent works correctly.  No need to test
  56              whether HALF_BIT is set.  */
  57           lx += half_bit;
  58           hx += lx < half_bit;
  59         }
  60       lx &= ~(int_bit - 1);
  61     }
  62   else if (exponent == BIAS + MANT_DIG - 65)
  63     {
  64       /* Not necessarily an integer; integer bit is bottom of high
  65          word, half bit is top of low word.  */
  66       if (((hx & 1) | (lx & 0x7fffffffffffffffULL)) != 0)
  67         {
  68           lx += 0x8000000000000000ULL;
  69           hx += lx < 0x8000000000000000ULL;
  70         }
  71       lx = 0;
  72     }
  73   else if (exponent >= BIAS)
  74     {
  75       /* At least 1; not necessarily an integer, integer bit and half
  76          bit are in the high word.  Locate the bits with exponents 0
  77          and -1 (when the unbiased exponent is 0, the bit with
  78          exponent 0 is implicit, but as the bias is odd it is OK to
  79          take it from the low bit of the exponent).  */
  80       int int_pos = (BIAS + MANT_DIG - 65) - exponent;
  81       int half_pos = int_pos - 1;
  82       uint64_t half_bit = 1ULL << half_pos;
  83       uint64_t int_bit = 1ULL << int_pos;
  84       if (((hx & (int_bit | (half_bit - 1))) | lx) != 0)
  85         hx += half_bit;
  86       hx &= ~(int_bit - 1);
  87       lx = 0;
  88     }
  89   else if (exponent == BIAS - 1 && (uhx > 0x3ffe000000000000ULL || lx != 0))
  90     {
  91       /* Interval (0.5, 1).  */
  92       hx = (hx & 0x8000000000000000ULL) | 0x3fff000000000000ULL;
  93       lx = 0;
  94     }
  95   else
  96     {
  97       /* Rounds to 0.  */
  98       hx &= 0x8000000000000000ULL;
  99       lx = 0;
 100     }
 101   SET_LDOUBLE_WORDS64 (x, hx, lx);
 102   return x;
 103 }
 104 libm_alias_ldouble (__roundeven, roundeven)