1  /* Round to nearest integer value, rounding halfway cases to even.
       2     ldbl-128 version.
       3     Copyright (C) 2016-2023 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     <https://www.gnu.org/licenses/>.  */
      19  
      20  #define NO_MATH_REDIRECT
      21  #include <math.h>
      22  #include <math_private.h>
      23  #include <libm-alias-ldouble.h>
      24  #include <math-use-builtins.h>
      25  #include <stdint.h>
      26  
      27  #define BIAS 0x3fff
      28  #define MANT_DIG 113
      29  #define MAX_EXP (2 * BIAS + 1)
      30  
      31  _Float128
      32  __roundevenl (_Float128 x)
      33  {
      34  #if USE_ROUNDEVENL_BUILTIN
      35    return __builtin_roundevenl (x);
      36  #else
      37    uint64_t hx, lx, uhx;
      38    GET_LDOUBLE_WORDS64 (hx, lx, x);
      39    uhx = hx & 0x7fffffffffffffffULL;
      40    int exponent = uhx >> (MANT_DIG - 1 - 64);
      41    if (exponent >= BIAS + MANT_DIG - 1)
      42      {
      43        /* Integer, infinity or NaN.  */
      44        if (exponent == MAX_EXP)
      45  	/* Infinity or NaN; quiet signaling NaNs.  */
      46  	return x + x;
      47        else
      48  	return x;
      49      }
      50    else if (exponent >= BIAS + MANT_DIG - 64)
      51      {
      52        /* Not necessarily an integer; integer bit is in low word.
      53  	 Locate the bits with exponents 0 and -1.  */
      54        int int_pos = (BIAS + MANT_DIG - 1) - exponent;
      55        int half_pos = int_pos - 1;
      56        uint64_t half_bit = 1ULL << half_pos;
      57        uint64_t int_bit = 1ULL << int_pos;
      58        if ((lx & (int_bit | (half_bit - 1))) != 0)
      59  	{
      60  	  /* Carry into the exponent works correctly.  No need to test
      61  	     whether HALF_BIT is set.  */
      62  	  lx += half_bit;
      63  	  hx += lx < half_bit;
      64  	}
      65        lx &= ~(int_bit - 1);
      66      }
      67    else if (exponent == BIAS + MANT_DIG - 65)
      68      {
      69        /* Not necessarily an integer; integer bit is bottom of high
      70  	 word, half bit is top of low word.  */
      71        if (((hx & 1) | (lx & 0x7fffffffffffffffULL)) != 0)
      72  	{
      73  	  lx += 0x8000000000000000ULL;
      74  	  hx += lx < 0x8000000000000000ULL;
      75  	}
      76        lx = 0;
      77      }
      78    else if (exponent >= BIAS)
      79      {
      80        /* At least 1; not necessarily an integer, integer bit and half
      81  	 bit are in the high word.  Locate the bits with exponents 0
      82  	 and -1 (when the unbiased exponent is 0, the bit with
      83  	 exponent 0 is implicit, but as the bias is odd it is OK to
      84  	 take it from the low bit of the exponent).  */
      85        int int_pos = (BIAS + MANT_DIG - 65) - exponent;
      86        int half_pos = int_pos - 1;
      87        uint64_t half_bit = 1ULL << half_pos;
      88        uint64_t int_bit = 1ULL << int_pos;
      89        if (((hx & (int_bit | (half_bit - 1))) | lx) != 0)
      90  	hx += half_bit;
      91        hx &= ~(int_bit - 1);
      92        lx = 0;
      93      }
      94    else if (exponent == BIAS - 1 && (uhx > 0x3ffe000000000000ULL || lx != 0))
      95      {
      96        /* Interval (0.5, 1).  */
      97        hx = (hx & 0x8000000000000000ULL) | 0x3fff000000000000ULL;
      98        lx = 0;
      99      }
     100    else
     101      {
     102        /* Rounds to 0.  */
     103        hx &= 0x8000000000000000ULL;
     104        lx = 0;
     105      }
     106    SET_LDOUBLE_WORDS64 (x, hx, lx);
     107    return x;
     108  #endif /* ! USE_ROUNDEVENL_BUILTIN  */
     109  }
     110  libm_alias_ldouble (__roundeven, roundeven)