Update.
[glibc.git] / sysdeps / libm-ieee754 / s_nearbyint.c
blob32f5bf9447285165548793c8a67c2d1920d8067d
1 /* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. */
2 /*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
13 #if defined(LIBM_SCCS) && !defined(lint)
14 static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $";
15 #endif
18 * rint(x)
19 * Return x rounded to integral value according to the prevailing
20 * rounding mode.
21 * Method:
22 * Using floating addition.
23 * Exception:
24 * Inexact flag raised if x not equal to rint(x).
27 #include <fenv.h>
28 #include "math.h"
29 #include "math_private.h"
31 #ifdef __STDC__
32 static const double
33 #else
34 static double
35 #endif
36 TWO52[2]={
37 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
38 -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
41 #ifdef __STDC__
42 double __nearbyint(double x)
43 #else
44 double __nearbyint(x)
45 double x;
46 #endif
48 fenv_t env;
49 int32_t i0,j0,sx;
50 u_int32_t i,i1;
51 double w,t;
52 EXTRACT_WORDS(i0,i1,x);
53 sx = (i0>>31)&1;
54 j0 = ((i0>>20)&0x7ff)-0x3ff;
55 if(j0<20) {
56 if(j0<0) {
57 if(((i0&0x7fffffff)|i1)==0) return x;
58 i1 |= (i0&0x0fffff);
59 i0 &= 0xfffe0000;
60 i0 |= ((i1|-i1)>>12)&0x80000;
61 SET_HIGH_WORD(x,i0);
62 feholdexcept (&env);
63 w = TWO52[sx]+x;
64 t = w-TWO52[sx];
65 fesetenv (&env);
66 GET_HIGH_WORD(i0,t);
67 SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
68 return t;
69 } else {
70 i = (0x000fffff)>>j0;
71 if(((i0&i)|i1)==0) return x; /* x is integral */
72 i>>=1;
73 if(((i0&i)|i1)!=0) {
74 if(j0==19) i1 = 0x40000000; else
75 i0 = (i0&(~i))|((0x20000)>>j0);
78 } else if (j0>51) {
79 if(j0==0x400) return x+x; /* inf or NaN */
80 else return x; /* x is integral */
81 } else {
82 i = ((u_int32_t)(0xffffffff))>>(j0-20);
83 if((i1&i)==0) return x; /* x is integral */
84 i>>=1;
85 if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
87 INSERT_WORDS(x,i0,i1);
88 feholdexcept (&env);
89 w = TWO52[sx]+x;
90 t = w-TWO52[sx];
91 fesetenv (&env);
92 return t;
94 weak_alias (__nearbyint, nearbyint)
95 #ifdef NO_LONG_DOUBLE
96 strong_alias (__nearbyint, __nearbyintl)
97 weak_alias (__nearbyint, nearbyintl)
98 #endif