1  /* Round to nearest integer value, rounding halfway cases to even.
       2     ldbl-96 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 <stdint.h>
      25  
      26  #define BIAS 0x3fff
      27  #define MANT_DIG 64
      28  #define MAX_EXP (2 * BIAS + 1)
      29  
      30  long double
      31  __roundevenl (long double x)
      32  {
      33    uint16_t se;
      34    uint32_t hx, lx;
      35    GET_LDOUBLE_WORDS (se, hx, lx, x);
      36    int exponent = se & 0x7fff;
      37    if (exponent >= BIAS + MANT_DIG - 1)
      38      {
      39        /* Integer, infinity or NaN.  */
      40        if (exponent == MAX_EXP)
      41  	/* Infinity or NaN; quiet signaling NaNs.  */
      42  	return x + x;
      43        else
      44  	return x;
      45      }
      46    else if (exponent >= BIAS + MANT_DIG - 32)
      47      {
      48        /* Not necessarily an integer; integer bit is in low word.
      49  	 Locate the bits with exponents 0 and -1.  */
      50        int int_pos = (BIAS + MANT_DIG - 1) - exponent;
      51        int half_pos = int_pos - 1;
      52        uint32_t half_bit = 1U << half_pos;
      53        uint32_t int_bit = 1U << int_pos;
      54        if ((lx & (int_bit | (half_bit - 1))) != 0)
      55  	{
      56  	  /* No need to test whether HALF_BIT is set.  */
      57  	  lx += half_bit;
      58  	  if (lx < half_bit)
      59  	    {
      60  	      hx++;
      61  	      if (hx == 0)
      62  		{
      63  		  hx = 0x80000000;
      64  		  se++;
      65  		}
      66  	    }
      67  	}
      68        lx &= ~(int_bit - 1);
      69      }
      70    else if (exponent == BIAS + MANT_DIG - 33)
      71      {
      72        /* Not necessarily an integer; integer bit is bottom of high
      73  	 word, half bit is top of low word.  */
      74        if (((hx & 1) | (lx & 0x7fffffff)) != 0)
      75  	{
      76  	  lx += 0x80000000;
      77  	  if (lx < 0x80000000)
      78  	    {
      79  	      hx++;
      80  	      if (hx == 0)
      81  		{
      82  		  hx = 0x80000000;
      83  		  se++;
      84  		}
      85  	    }
      86  	}
      87        lx = 0;
      88      }
      89    else if (exponent >= BIAS)
      90      {
      91        /* At least 1; not necessarily an integer, integer bit and half
      92  	 bit are in the high word.  Locate the bits with exponents 0
      93  	 and -1.  */
      94        int int_pos = (BIAS + MANT_DIG - 33) - exponent;
      95        int half_pos = int_pos - 1;
      96        uint32_t half_bit = 1U << half_pos;
      97        uint32_t int_bit = 1U << int_pos;
      98        if (((hx & (int_bit | (half_bit - 1))) | lx) != 0)
      99  	{
     100  	  hx += half_bit;
     101  	  if (hx < half_bit)
     102  	    {
     103  	      hx = 0x80000000;
     104  	      se++;
     105  	    }
     106  	}
     107        hx &= ~(int_bit - 1);
     108        lx = 0;
     109      }
     110    else if (exponent == BIAS - 1 && (hx > 0x80000000 || lx != 0))
     111      {
     112        /* Interval (0.5, 1).  */
     113        se = (se & 0x8000) | 0x3fff;
     114        hx = 0x80000000;
     115        lx = 0;
     116      }
     117    else
     118      {
     119        /* Rounds to 0.  */
     120        se &= 0x8000;
     121        hx = 0;
     122        lx = 0;
     123      }
     124    SET_LDOUBLE_WORDS (x, se, hx, lx);
     125    return x;
     126  }
     127  libm_alias_ldouble (__roundeven, roundeven)