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-2014 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 high = u.d[0].d;
  42       feholdexcept (&env);
  43       if (high > 0.0)
  44         {
  45           high += TWO52;
  46           high -= TWO52;
  47           if (high == -0.0) high = 0.0;
  48         }
  49       else if (high < 0.0)
  50         {
  51           high -= TWO52;
  52           high += TWO52;
  53           if (high == 0.0) high = -0.0;
  54         }
  55       u.d[0].d = high;
  56       u.d[1].d = 0.0;
  57       math_force_eval (u.d[0]);
  58       math_force_eval (u.d[1]);
  59       fesetenv (&env);
  60     }
  61   else if (fabs (u.d[1].d) < TWO52 && u.d[1].d != 0.0)
  62     {
  63       double high, low, tau;
  64       /* In this case we have to round the low double and handle any
  65          adjustment to the high double that may be caused by rounding
  66          (up).  This is complicated by the fact that the high double
  67          may already be rounded and the low double may have the
  68          opposite sign to compensate.  */
  69       feholdexcept (&env);
  70       if (u.d[0].d > 0.0)
  71         {
  72           if (u.d[1].d > 0.0)
  73             {
  74               /* If the high/low doubles are the same sign then simply
  75                  round the low double.  */
  76               high = u.d[0].d;
  77               low = u.d[1].d;
  78             }
  79           else if (u.d[1].d < 0.0)
  80             {
  81               /* Else the high double is pre rounded and we need to
  82                  adjust for that.  */
  83
  84               tau = __nextafter (u.d[0].d, 0.0);
  85               tau = (u.d[0].d - tau) * 2.0;
  86               high = u.d[0].d - tau;
  87               low = u.d[1].d + tau;
  88             }
  89           low += TWO52;
  90           low -= TWO52;
  91         }
  92       else if (u.d[0].d < 0.0)
  93         {
  94           if (u.d[1].d < 0.0)
  95             {
  96               /* If the high/low doubles are the same sign then simply
  97                  round the low double.  */
  98               high = u.d[0].d;
  99               low = u.d[1].d;
 100             }
 101           else if (u.d[1].d > 0.0)
 102             {
 103               /* Else the high double is pre rounded and we need to
 104                  adjust for that.  */
 105               tau = __nextafter (u.d[0].d, 0.0);
 106               tau = (u.d[0].d - tau) * 2.0;
 107               high = u.d[0].d - tau;
 108               low = u.d[1].d + tau;
 109             }
 110           low = TWO52 - low;
 111           low = -(low - TWO52);
 112         }
 113       u.d[0].d = high + low;
 114       u.d[1].d = high - u.d[0].d + low;
 115       math_force_eval (u.d[0]);
 116       math_force_eval (u.d[1]);
 117       fesetenv (&env);
 118     }
 119
 120   return u.ld;
 121 }
 122
 123 long_double_symbol (libm, __nearbyintl, nearbyintl);