sysdeps/ieee754/ldbl-96/s_roundl.c

   1 /* Round long double to integer away from zero.
   2    Copyright (C) 1997-2014 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    <http://www.gnu.org/licenses/>.  */
  19
  20 #include <math.h>
  21
  22 #include <math_private.h>
  23
  24
  25 static const long double huge = 1.0e4930L;
  26
  27
  28 long double
  29 __roundl (long double x)
  30 {
  31   int32_t j0;
  32   u_int32_t se, i1, i0;
  33
  34   GET_LDOUBLE_WORDS (se, i0, i1, x);
  35   j0 = (se & 0x7fff) - 0x3fff;
  36   if (j0 < 31)
  37     {
  38       if (j0 < 0)
  39         {
  40           math_force_eval (huge + x);
  41           se &= 0x8000;
  42           i0 = i1 = 0;
  43           if (j0 == -1)
  44             {
  45               se |= 0x3fff;
  46               i0 = 0x80000000;
  47             }
  48         }
  49       else
  50         {
  51           u_int32_t i = 0x7fffffff >> j0;
  52           if (((i0 & i) | i1) == 0)
  53             /* X is integral.  */
  54             return x;
  55
  56           /* Raise inexact if x != 0.  */
  57           math_force_eval (huge + x);
  58           u_int32_t j = i0 + (0x40000000 >> j0);
  59           if (j < i0)
  60             se += 1;
  61           i0 = (j & ~i) | 0x80000000;
  62           i1 = 0;
  63         }
  64     }
  65   else if (j0 > 62)
  66     {
  67       if (j0 == 0x4000)
  68         /* Inf or NaN.  */
  69         return x + x;
  70       else
  71         return x;
  72     }
  73   else
  74     {
  75       u_int32_t i = 0xffffffff >> (j0 - 31);
  76       if ((i1 & i) == 0)
  77         /* X is integral.  */
  78         return x;
  79
  80       math_force_eval (huge + x);
  81       /* Raise inexact if x != 0.  */
  82       u_int32_t j = i1 + (1 << (62 - j0));
  83       if (j < i1)
  84         {
  85           u_int32_t k = i0 + 1;
  86           if (k < i0)
  87             {
  88               se += 1;
  89               k |= 0x80000000;
  90             }
  91           i0 = k;
  92         }
  93       i1 = j;
  94       i1 &= ~i;
  95     }
  96
  97   SET_LDOUBLE_WORDS (x, se, i0, i1);
  98   return x;
  99 }
 100 weak_alias (__roundl, roundl)