sysdeps/ieee754/ldbl-96/s_roundl.c

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