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-2017 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 <stdint.h>
  23
  24 #define BIAS 0x3fff
  25 #define MANT_DIG 113
  26 #define MAX_EXP (2 * BIAS + 1)
  27
  28 _Float128
  29 roundevenl (_Float128 x)
  30 {
  31   uint64_t hx, lx, uhx;
  32   GET_LDOUBLE_WORDS64 (hx, lx, x);
  33   uhx = hx & 0x7fffffffffffffffULL;
  34   int exponent = uhx >> (MANT_DIG - 1 - 64);
  35   if (exponent >= BIAS + MANT_DIG - 1)
  36     {
  37       /* Integer, infinity or NaN.  */
  38       if (exponent == MAX_EXP)
  39         /* Infinity or NaN; quiet signaling NaNs.  */
  40         return x + x;
  41       else
  42         return x;
  43     }
  44   else if (exponent >= BIAS + MANT_DIG - 64)
  45     {
  46       /* Not necessarily an integer; integer bit is in low word.
  47          Locate the bits with exponents 0 and -1.  */
  48       int int_pos = (BIAS + MANT_DIG - 1) - exponent;
  49       int half_pos = int_pos - 1;
  50       uint64_t half_bit = 1ULL << half_pos;
  51       uint64_t int_bit = 1ULL << int_pos;
  52       if ((lx & (int_bit | (half_bit - 1))) != 0)
  53         {
  54           /* Carry into the exponent works correctly.  No need to test
  55              whether HALF_BIT is set.  */
  56           lx += half_bit;
  57           hx += lx < half_bit;
  58         }
  59       lx &= ~(int_bit - 1);
  60     }
  61   else if (exponent == BIAS + MANT_DIG - 65)
  62     {
  63       /* Not necessarily an integer; integer bit is bottom of high
  64          word, half bit is top of low word.  */
  65       if (((hx & 1) | (lx & 0x7fffffffffffffffULL)) != 0)
  66         {
  67           lx += 0x8000000000000000ULL;
  68           hx += lx < 0x8000000000000000ULL;
  69         }
  70       lx = 0;
  71     }
  72   else if (exponent >= BIAS)
  73     {
  74       /* At least 1; not necessarily an integer, integer bit and half
  75          bit are in the high word.  Locate the bits with exponents 0
  76          and -1 (when the unbiased exponent is 0, the bit with
  77          exponent 0 is implicit, but as the bias is odd it is OK to
  78          take it from the low bit of the exponent).  */
  79       int int_pos = (BIAS + MANT_DIG - 65) - exponent;
  80       int half_pos = int_pos - 1;
  81       uint64_t half_bit = 1ULL << half_pos;
  82       uint64_t int_bit = 1ULL << int_pos;
  83       if (((hx & (int_bit | (half_bit - 1))) | lx) != 0)
  84         hx += half_bit;
  85       hx &= ~(int_bit - 1);
  86       lx = 0;
  87     }
  88   else if (exponent == BIAS - 1 && (uhx > 0x3ffe000000000000ULL || lx != 0))
  89     {
  90       /* Interval (0.5, 1).  */
  91       hx = (hx & 0x8000000000000000ULL) | 0x3fff000000000000ULL;
  92       lx = 0;
  93     }
  94   else
  95     {
  96       /* Rounds to 0.  */
  97       hx &= 0x8000000000000000ULL;
  98       lx = 0;
  99     }
 100   SET_LDOUBLE_WORDS64 (x, hx, lx);
 101   return x;
 102 }