sysdeps/ieee754/ldbl-128ibm/s_nearbyintl.c

   1 /* Round to int long double floating-point values without raising inexact.
   2    IBM extended format long double version.
   3    Copyright (C) 2006-2015 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    <http://www.gnu.org/licenses/>.  */
  19
  20 /* This has been coded in assembler because GCC makes such a mess of it
  21    when it's coded in C.  */
  22
  23 #include <math.h>
  24 #include <math_private.h>
  25 #include <fenv.h>
  26 #include <math_ldbl_opt.h>
  27 #include <float.h>
  28 #include <ieee754.h>
  29
  30
  31 long double
  32 __nearbyintl (long double x)
  33 {
  34   fenv_t env;
  35   static const long double TWO52 = 4503599627370496.0L;
  36   union ibm_extended_long_double u;
  37   u.ld = x;
  38
  39   if (fabs (u.d[0].d) < TWO52)
  40     {
  41       double xh = u.d[0].d;
  42       double high = u.d[0].d;
  43       feholdexcept (&env);
  44       if (high > 0.0)
  45         {
  46           high += TWO52;
  47           high -= TWO52;
  48           if (high == -0.0) high = 0.0;
  49         }
  50       else if (high < 0.0)
  51         {
  52           high -= TWO52;
  53           high += TWO52;
  54           if (high == 0.0) high = -0.0;
  55         }
  56       if (u.d[1].d > 0.0 && (xh - high == 0.5))
  57         high += 1.0;
  58       else if (u.d[1].d < 0.0 && (-(xh - high) == 0.5))
  59         high -= 1.0;
  60       u.d[0].d = high;
  61       u.d[1].d = 0.0;
  62       math_force_eval (u.d[0]);
  63       math_force_eval (u.d[1]);
  64       fesetenv (&env);
  65     }
  66   else if (fabs (u.d[1].d) < TWO52 && u.d[1].d != 0.0)
  67     {
  68       double high, low, tau;
  69       /* In this case we have to round the low double and handle any
  70          adjustment to the high double that may be caused by rounding
  71          (up).  This is complicated by the fact that the high double
  72          may already be rounded and the low double may have the
  73          opposite sign to compensate.  */
  74       feholdexcept (&env);
  75       if (u.d[0].d > 0.0)
  76         {
  77           if (u.d[1].d > 0.0)
  78             {
  79               /* If the high/low doubles are the same sign then simply
  80                  round the low double.  */
  81               high = u.d[0].d;
  82               low = u.d[1].d;
  83             }
  84           else if (u.d[1].d < 0.0)
  85             {
  86               /* Else the high double is pre rounded and we need to
  87                  adjust for that.  */
  88
  89               tau = __nextafter (u.d[0].d, 0.0);
  90               tau = (u.d[0].d - tau) * 2.0;
  91               high = u.d[0].d - tau;
  92               low = u.d[1].d + tau;
  93             }
  94           low += TWO52;
  95           low -= TWO52;
  96         }
  97       else if (u.d[0].d < 0.0)
  98         {
  99           if (u.d[1].d < 0.0)
 100             {
 101               /* If the high/low doubles are the same sign then simply
 102                  round the low double.  */
 103               high = u.d[0].d;
 104               low = u.d[1].d;
 105             }
 106           else if (u.d[1].d > 0.0)
 107             {
 108               /* Else the high double is pre rounded and we need to
 109                  adjust for that.  */
 110               tau = __nextafter (u.d[0].d, 0.0);
 111               tau = (u.d[0].d - tau) * 2.0;
 112               high = u.d[0].d - tau;
 113               low = u.d[1].d + tau;
 114             }
 115           low = TWO52 - low;
 116           low = -(low - TWO52);
 117         }
 118       u.d[0].d = high + low;
 119       u.d[1].d = high - u.d[0].d + low;
 120       math_force_eval (u.d[0]);
 121       math_force_eval (u.d[1]);
 122       fesetenv (&env);
 123     }
 124
 125   return u.ld;
 126 }
 127
 128 long_double_symbol (libm, __nearbyintl, nearbyintl);