Back to index

glibc  2.9
s_roundl.c
Go to the documentation of this file.
00001 /* Round to int long double floating-point values.
00002    IBM extended format long double version.
00003    Copyright (C) 2006, 2007 Free Software Foundation, Inc.
00004    This file is part of the GNU C Library.
00005 
00006    The GNU C Library is free software; you can redistribute it and/or
00007    modify it under the terms of the GNU Lesser General Public
00008    License as published by the Free Software Foundation; either
00009    version 2.1 of the License, or (at your option) any later version.
00010 
00011    The GNU C Library is distributed in the hope that it will be useful,
00012    but WITHOUT ANY WARRANTY; without even the implied warranty of
00013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014    Lesser General Public License for more details.
00015 
00016    You should have received a copy of the GNU Lesser General Public
00017    License along with the GNU C Library; if not, write to the Free
00018    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00019    02111-1307 USA.  */
00020 
00021 /* This has been coded in assembler because GCC makes such a mess of it
00022    when it's coded in C.  */
00023 
00024 #include <math.h>
00025 #include <math_ldbl_opt.h>
00026 #include <float.h>
00027 #include <ieee754.h>
00028 
00029 
00030 #ifdef __STDC__
00031 long double
00032 __roundl (long double x)
00033 #else
00034 long double
00035 __roundl (x)
00036      long double x;
00037 #endif
00038 {
00039   double xh, xl, hi, lo;
00040 
00041   ldbl_unpack (x, &xh, &xl);
00042 
00043   /* Return Inf, Nan, +/-0 unchanged.  */
00044   if (__builtin_expect (xh != 0.0
00045                      && __builtin_isless (__builtin_fabs (xh),
00046                                         __builtin_inf ()), 1))
00047     {
00048       double orig_xh;
00049 
00050       /* Long double arithmetic, including the canonicalisation below,
00051         only works in round-to-nearest mode.  */
00052 
00053       /* Convert the high double to integer.  */
00054       orig_xh = xh;
00055       hi = ldbl_nearbyint (xh);
00056 
00057       /* Subtract integral high part from the value.  */
00058       xh -= hi;
00059       ldbl_canonicalize (&xh, &xl);
00060 
00061       /* Now convert the low double, adjusted for any remainder from the
00062         high double.  */
00063       lo = ldbl_nearbyint (xh);
00064 
00065       /* Adjust the result when the remainder is exactly 0.5.  nearbyint
00066         rounds values halfway between integers to the nearest even
00067         integer.  roundl must round away from zero.
00068         Also correct cases where nearbyint returns an incorrect value
00069         for LO.  */
00070       xh -= lo;
00071       ldbl_canonicalize (&xh, &xl);
00072       if (xh == 0.5)
00073        {
00074          if (xl > 0.0 || (xl == 0.0 && orig_xh > 0.0))
00075            lo += 1.0;
00076        }
00077       else if (-xh == 0.5)
00078        {
00079          if (xl < 0.0 || (xl == 0.0 && orig_xh < 0.0))
00080            lo -= 1.0;
00081        }
00082 
00083       /* Ensure the final value is canonical.  In certain cases,
00084         rounding causes hi,lo calculated so far to be non-canonical.  */
00085       xh = hi;
00086       xl = lo;
00087       ldbl_canonicalize (&xh, &xl);
00088     }
00089 
00090   return ldbl_pack (xh, xl);
00091 }
00092 
00093 long_double_symbol (libm, __roundl, roundl);