sysdeps/ieee754/ldbl-96/s_roundl.c

   1 /* Round long double to integer away from zero.
   2    Copyright (C) 1997-2024 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
  22 #include <math_private.h>
  23 #include <libm-alias-ldouble.h>
  24
  25
  26 long double
  27 __roundl (long double x)
  28 {
  29   int32_t j0;
  30   uint32_t se, i1, i0;
  31
  32   GET_LDOUBLE_WORDS (se, i0, i1, x);
  33   j0 = (se & 0x7fff) - 0x3fff;
  34   if (j0 < 31)
  35     {
  36       if (j0 < 0)
  37         {
  38           se &= 0x8000;
  39           i0 = i1 = 0;
  40           if (j0 == -1)
  41             {
  42               se |= 0x3fff;
  43               i0 = 0x80000000;
  44             }
  45         }
  46       else
  47         {
  48           uint32_t i = 0x7fffffff >> j0;
  49           if (((i0 & i) | i1) == 0)
  50             /* X is integral.  */
  51             return x;
  52
  53           uint32_t j = i0 + (0x40000000 >> j0);
  54           if (j < i0)
  55             se += 1;
  56           i0 = (j & ~i) | 0x80000000;
  57           i1 = 0;
  58         }
  59     }
  60   else if (j0 > 62)
  61     {
  62       if (j0 == 0x4000)
  63         /* Inf or NaN.  */
  64         return x + x;
  65       else
  66         return x;
  67     }
  68   else
  69     {
  70       uint32_t i = 0xffffffff >> (j0 - 31);
  71       if ((i1 & i) == 0)
  72         /* X is integral.  */
  73         return x;
  74
  75       uint32_t j = i1 + (1 << (62 - j0));
  76       if (j < i1)
  77         {
  78           uint32_t k = i0 + 1;
  79           if (k < i0)
  80             {
  81               se += 1;
  82               k |= 0x80000000;
  83             }
  84           i0 = k;
  85         }
  86       i1 = j;
  87       i1 &= ~i;
  88     }
  89
  90   SET_LDOUBLE_WORDS (x, se, i0, i1);
  91   return x;
  92 }
  93 libm_alias_ldouble (__round, round)