Back to index

glibc  2.9
s_nearbyint.c
Go to the documentation of this file.
00001 /* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>.  */
00002 /*
00003  * ====================================================
00004  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00005  *
00006  * Developed at SunPro, a Sun Microsystems, Inc. business.
00007  * Permission to use, copy, modify, and distribute this
00008  * software is freely granted, provided that this notice
00009  * is preserved.
00010  * ====================================================
00011  */
00012 
00013 #if defined(LIBM_SCCS) && !defined(lint)
00014 static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $";
00015 #endif
00016 
00017 /*
00018  * rint(x)
00019  * Return x rounded to integral value according to the prevailing
00020  * rounding mode.
00021  * Method:
00022  *     Using floating addition.
00023  * Exception:
00024  *     Inexact flag raised if x not equal to rint(x).
00025  */
00026 
00027 #include <fenv.h>
00028 #include "math.h"
00029 #include "math_private.h"
00030 
00031 #ifdef __STDC__
00032 static const double
00033 #else
00034 static double
00035 #endif
00036 TWO52[2]={
00037   4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
00038  -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
00039 };
00040 
00041 #ifdef __STDC__
00042        double __nearbyint(double x)
00043 #else
00044        double __nearbyint(x)
00045        double x;
00046 #endif
00047 {
00048        fenv_t env;
00049        int32_t i0,j0,sx;
00050        u_int32_t i,i1;
00051        double w,t;
00052        EXTRACT_WORDS(i0,i1,x);
00053        sx = (i0>>31)&1;
00054        j0 = ((i0>>20)&0x7ff)-0x3ff;
00055        if(j0<20) {
00056            if(j0<0) {
00057               if(((i0&0x7fffffff)|i1)==0) return x;
00058               i1 |= (i0&0x0fffff);
00059               i0 &= 0xfffe0000;
00060               i0 |= ((i1|-i1)>>12)&0x80000;
00061               SET_HIGH_WORD(x,i0);
00062               feholdexcept (&env);
00063                w = TWO52[sx]+x;
00064                t =  w-TWO52[sx];
00065               fesetenv (&env);
00066               GET_HIGH_WORD(i0,t);
00067               SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
00068                return t;
00069            } else {
00070               i = (0x000fffff)>>j0;
00071               if(((i0&i)|i1)==0) return x; /* x is integral */
00072               i>>=1;
00073               if(((i0&i)|i1)!=0) {
00074                 if (j0==19)
00075                   i1 = 0x40000000;
00076                 else if (j0<18)
00077                   i0 = (i0&(~i))|((0x20000)>>j0);
00078                 else
00079                   {
00080                     i0 &= ~i;
00081                     i1 = 0x80000000;
00082                   }
00083               }
00084            }
00085        } else if (j0>51) {
00086            if(j0==0x400) return x+x;      /* inf or NaN */
00087            else return x;          /* x is integral */
00088        } else {
00089            i = ((u_int32_t)(0xffffffff))>>(j0-20);
00090            if((i1&i)==0) return x; /* x is integral */
00091            i>>=1;
00092            if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
00093        }
00094        INSERT_WORDS(x,i0,i1);
00095        feholdexcept (&env);
00096        w = TWO52[sx]+x;
00097        t = w-TWO52[sx];
00098        fesetenv (&env);
00099        return t;
00100 }
00101 weak_alias (__nearbyint, nearbyint)
00102 #ifdef NO_LONG_DOUBLE
00103 strong_alias (__nearbyint, __nearbyintl)
00104 weak_alias (__nearbyint, nearbyintl)
00105 #endif