sysdeps/ieee754/ldbl-96/s_roundl.c

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