Back to index

glibc  2.9
s_nearbyintl.c
Go to the documentation of this file.
00001 /* s_nearbyintl.c -- long double version of s_nearbyint.c.
00002  * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
00003  */
00004 
00005 /*
00006  * ====================================================
00007  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00008  *
00009  * Developed at SunPro, a Sun Microsystems, Inc. business.
00010  * Permission to use, copy, modify, and distribute this
00011  * software is freely granted, provided that this notice
00012  * is preserved.
00013  * ====================================================
00014  */
00015 
00016 /*
00017  * nearbyintl(x)
00018  * Return x rounded to integral value according to the prevailing
00019  * rounding mode.
00020  * Method:
00021  *     Using floating addition.
00022  * Exception:
00023  *     Inexact flag raised if x not equal to rintl(x).
00024  */
00025 
00026 #include <fenv.h>
00027 #include "math.h"
00028 #include "math_private.h"
00029 
00030 #ifdef __STDC__
00031 static const long double
00032 #else
00033 static long double
00034 #endif
00035 TWO112[2]={
00036   5.19229685853482762853049632922009600E+33L, /* 0x406F000000000000, 0 */
00037  -5.19229685853482762853049632922009600E+33L  /* 0xC06F000000000000, 0 */
00038 };
00039 
00040 #ifdef __STDC__
00041        long double __nearbyintl(long double x)
00042 #else
00043        long double __nearbyintl(x)
00044        long double x;
00045 #endif
00046 {
00047        fenv_t env;
00048        int64_t i0,j0,sx;
00049        u_int64_t i,i1;
00050        long double w,t;
00051        GET_LDOUBLE_WORDS64(i0,i1,x);
00052        sx = (((u_int64_t)i0)>>63);
00053        j0 = ((i0>>48)&0x7fff)-0x3fff;
00054        if(j0<48) {
00055            if(j0<0) {
00056               if(((i0&0x7fffffffffffffffLL)|i1)==0) return x;
00057               i1 |= (i0&0x0000ffffffffffffLL);
00058               i0 &= 0xffffe00000000000ULL;
00059               i0 |= ((i1|-i1)>>16)&0x0000800000000000LL;
00060               SET_LDOUBLE_MSW64(x,i0);
00061               feholdexcept (&env);
00062                w = TWO112[sx]+x;
00063                t = w-TWO112[sx];
00064                fesetenv (&env);
00065               GET_LDOUBLE_MSW64(i0,t);
00066               SET_LDOUBLE_MSW64(t,(i0&0x7fffffffffffffffLL)|(sx<<63));
00067                return t;
00068            } else {
00069               i = (0x0000ffffffffffffLL)>>j0;
00070               if(((i0&i)|i1)==0) return x; /* x is integral */
00071               i>>=1;
00072               if(((i0&i)|i1)!=0) {
00073                   if(j0==47) i1 = 0x4000000000000000ULL; else
00074                   i0 = (i0&(~i))|((0x0000200000000000LL)>>j0);
00075               }
00076            }
00077        } else if (j0>111) {
00078            if(j0==0x4000) return x+x;     /* inf or NaN */
00079            else return x;          /* x is integral */
00080        } else {
00081            i = -1ULL>>(j0-48);
00082            if((i1&i)==0) return x; /* x is integral */
00083            i>>=1;
00084            if((i1&i)!=0) i1 = (i1&(~i))|((0x4000000000000000LL)>>(j0-48));
00085        }
00086        SET_LDOUBLE_WORDS64(x,i0,i1);
00087        feholdexcept (&env);
00088        w = TWO112[sx]+x;
00089        t = w-TWO112[sx];
00090        fesetenv (&env);
00091        return t;
00092 }
00093 weak_alias (__nearbyintl, nearbyintl)