1  /* Round to nearest integer value, rounding halfway cases to even.
       2     Copyright (C) 2016-2023 Free Software Foundation, Inc.
       3     This file is part of the GNU C Library.
       4  
       5     The GNU C Library is free software; you can redistribute it and/or
       6     modify it under the terms of the GNU Lesser General Public
       7     License as published by the Free Software Foundation; either
       8     version 2.1 of the License, or (at your option) any later version.
       9  
      10     The GNU C Library is distributed in the hope that it will be useful,
      11     but WITHOUT ANY WARRANTY; without even the implied warranty of
      12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      13     Lesser General Public License for more details.
      14  
      15     You should have received a copy of the GNU Lesser General Public
      16     License along with the GNU C Library; if not, see
      17     <https://www.gnu.org/licenses/>.  */
      18  
      19  #define NO_MATH_REDIRECT
      20  #include <math.h>
      21  #include <math_private.h>
      22  #include <libm-alias-double.h>
      23  #include <stdint.h>
      24  #include <math-use-builtins.h>
      25  
      26  #define BIAS 0x3ff
      27  #define MANT_DIG 53
      28  #define MAX_EXP (2 * BIAS + 1)
      29  
      30  double
      31  __roundeven (double x)
      32  {
      33  #if USE_ROUNDEVEN_BUILTIN
      34    return __builtin_roundeven (x);
      35  #else
      36    uint64_t ix, ux;
      37    EXTRACT_WORDS64 (ix, x);
      38    ux = ix & 0x7fffffffffffffffULL;
      39    int exponent = ux >> (MANT_DIG - 1);
      40    if (exponent >= BIAS + MANT_DIG - 1)
      41      {
      42        /* Integer, infinity or NaN.  */
      43        if (exponent == MAX_EXP)
      44  	/* Infinity or NaN; quiet signaling NaNs.  */
      45  	return x + x;
      46        else
      47  	return x;
      48      }
      49    else if (exponent >= BIAS)
      50      {
      51        /* At least 1; not necessarily an integer.  Locate the bits with
      52  	 exponents 0 and -1 (when the unbiased exponent is 0, the bit
      53  	 with exponent 0 is implicit, but as the bias is odd it is OK
      54  	 to take it from the low bit of the exponent).  */
      55        int int_pos = (BIAS + MANT_DIG - 1) - exponent;
      56        int half_pos = int_pos - 1;
      57        uint64_t half_bit = 1ULL << half_pos;
      58        uint64_t int_bit = 1ULL << int_pos;
      59        if ((ix & (int_bit | (half_bit - 1))) != 0)
      60  	/* Carry into the exponent works correctly.  No need to test
      61  	   whether HALF_BIT is set.  */
      62  	ix += half_bit;
      63        ix &= ~(int_bit - 1);
      64      }
      65    else if (exponent == BIAS - 1 && ux > 0x3fe0000000000000ULL)
      66      /* Interval (0.5, 1).  */
      67      ix = (ix & 0x8000000000000000ULL) | 0x3ff0000000000000ULL;
      68    else
      69      /* Rounds to 0.  */
      70      ix &= 0x8000000000000000ULL;
      71    INSERT_WORDS64 (x, ix);
      72    return x;
      73  #endif /* ! USE_ROUNDEVEN_BUILTIN  */
      74  }
      75  #ifndef __roundeven
      76  libm_alias_double (__roundeven, roundeven)
      77  #endif