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