1  /* Round to nearest integer value, rounding halfway cases to even.
       2     flt-32 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-float.h>
      24  #include <math-use-builtins.h>
      25  #include <stdint.h>
      26  
      27  #define BIAS 0x7f
      28  #define MANT_DIG 24
      29  #define MAX_EXP (2 * BIAS + 1)
      30  
      31  float
      32  __roundevenf (float x)
      33  {
      34  #if USE_ROUNDEVENF_BUILTIN
      35    return __builtin_roundevenf (x);
      36  #else
      37    uint32_t ix, ux;
      38    GET_FLOAT_WORD (ix, x);
      39    ux = ix & 0x7fffffff;
      40    int exponent = ux >> (MANT_DIG - 1);
      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)
      51      {
      52        /* At least 1; not necessarily an integer.  Locate the bits with
      53  	 exponents 0 and -1 (when the unbiased exponent is 0, the bit
      54  	 with exponent 0 is implicit, but as the bias is odd it is OK
      55  	 to take it from the low bit of the exponent).  */
      56        int int_pos = (BIAS + MANT_DIG - 1) - exponent;
      57        int half_pos = int_pos - 1;
      58        uint32_t half_bit = 1U << half_pos;
      59        uint32_t int_bit = 1U << int_pos;
      60        if ((ix & (int_bit | (half_bit - 1))) != 0)
      61  	/* Carry into the exponent works correctly.  No need to test
      62  	   whether HALF_BIT is set.  */
      63  	ix += half_bit;
      64        ix &= ~(int_bit - 1);
      65      }
      66    else if (exponent == BIAS - 1 && ux > 0x3f000000)
      67      /* Interval (0.5, 1).  */
      68      ix = (ix & 0x80000000) | 0x3f800000;
      69    else
      70      /* Rounds to 0.  */
      71      ix &= 0x80000000;
      72    SET_FLOAT_WORD (x, ix);
      73    return x;
      74  #endif /* ! USE_ROUNDEVENF_BUILTIN  */
      75  }
      76  #ifndef __roundevenf
      77  libm_alias_float (__roundeven, roundeven)
      78  #endif