Back to index

glibc  2.9
s_rintl.c
Go to the documentation of this file.
00001 /* s_rintl.c -- long double version of s_rint.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 #if defined(LIBM_SCCS) && !defined(lint)
00017 static char rcsid[] = "$NetBSD: $";
00018 #endif
00019 
00020 /*
00021  * rintl(x)
00022  * Return x rounded to integral value according to the prevailing
00023  * rounding mode.
00024  * Method:
00025  *     Using floating addition.
00026  * Exception:
00027  *     Inexact flag raised if x not equal to rintl(x).
00028  */
00029 
00030 #include "math.h"
00031 #include "math_private.h"
00032 
00033 #ifdef __STDC__
00034 static const long double
00035 #else
00036 static long double
00037 #endif
00038 TWO112[2]={
00039   5.19229685853482762853049632922009600E+33L, /* 0x406F000000000000, 0 */
00040  -5.19229685853482762853049632922009600E+33L  /* 0xC06F000000000000, 0 */
00041 };
00042 
00043 #ifdef __STDC__
00044        long double __rintl(long double x)
00045 #else
00046        long double __rintl(x)
00047        long double x;
00048 #endif
00049 {
00050        int64_t i0,j0,sx;
00051        u_int64_t i,i1;
00052        long double w,t;
00053        GET_LDOUBLE_WORDS64(i0,i1,x);
00054        sx = (((u_int64_t)i0)>>63);
00055        j0 = ((i0>>48)&0x7fff)-0x3fff;
00056        if(j0<48) {
00057            if(j0<0) {
00058               if(((i0&0x7fffffffffffffffLL)|i1)==0) return x;
00059               i1 |= (i0&0x0000ffffffffffffLL);
00060               i0 &= 0xffffe00000000000ULL;
00061               i0 |= ((i1|-i1)>>16)&0x0000800000000000LL;
00062               SET_LDOUBLE_MSW64(x,i0);
00063                w = TWO112[sx]+x;
00064                t = w-TWO112[sx];
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        w = TWO112[sx]+x;
00088        return w-TWO112[sx];
00089 }
00090 weak_alias (__rintl, rintl)