Update.
[glibc.git] / sysdeps / libm-ieee754 / s_nearbyintl.c
blob3b0715e5a53e87279167752db67415867014c318
1 /* s_rintl.c -- long double version of s_rint.c.
2 * Conversion to long double by Ulrich Drepper,
3 * Cygnus Support, drepper@cygnus.com.
4 */
5 /* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. */
7 /*
8 * ====================================================
9 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
11 * Developed at SunPro, a Sun Microsystems, Inc. business.
12 * Permission to use, copy, modify, and distribute this
13 * software is freely granted, provided that this notice
14 * is preserved.
15 * ====================================================
19 * rintl(x)
20 * Return x rounded to integral value according to the prevailing
21 * rounding mode.
22 * Method:
23 * Using floating addition.
24 * Exception:
25 * Inexact flag raised if x not equal to rintl(x).
28 #include <fenv.h>
29 #include "math.h"
30 #include "math_private.h"
32 #ifdef __STDC__
33 static const long double
34 #else
35 static long double
36 #endif
37 TWO63[2]={
38 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */
39 -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */
42 #ifdef __STDC__
43 long double __nearbyintl(long double x)
44 #else
45 long double __nearbyintl(x)
46 long double x;
47 #endif
49 fenv_t env;
50 int32_t se,j0,sx;
51 u_int32_t i,i0,i1;
52 long double w,t;
53 GET_LDOUBLE_WORDS(se,i0,i1,x);
54 sx = (se>>15)&1;
55 j0 = (se&0x7fff)-0x3fff;
56 if(j0<31) {
57 if(j0<0) {
58 if(((se&0x7fff)|i0|i1)==0) return x;
59 i1 |= i0;
60 i0 &= 0xe0000000;
61 i0 |= (i1|-i1)&0x80000000;
62 SET_LDOUBLE_MSW(x,i0);
63 feholdexcept (&env);
64 w = TWO63[sx]+x;
65 t = w-TWO63[sx];
66 fesetenv (&env);
67 GET_LDOUBLE_EXP(i0,t);
68 SET_LDOUBLE_EXP(t,(i0&0x7fff)|(sx<<15));
69 return t;
70 } else {
71 i = (0x7fffffff)>>j0;
72 if(((i0&i)|i1)==0) return x; /* x is integral */
73 i>>=1;
74 if(((i0&i)|i1)!=0) {
75 if(j0==31) i1 = 0x40000000; else
76 i0 = (i0&(~i))|((0x20000000)>>j0);
77 /* Shouldn't this be
78 if (j0 >= 30) i1 = 0x80000000 >> (j0 - 30);
79 i0 = (i0&(~i))|((0x20000000)>>j0);
80 If yes, this should be correct in s_rint and
81 s_rintf, too. -- drepper@cygnus.com */
84 } else if (j0>62) {
85 if(j0==0x4000) return x+x; /* inf or NaN */
86 else return x; /* x is integral */
87 } else {
88 i = ((u_int32_t)(0xffffffff))>>(j0-31);
89 if((i1&i)==0) return x; /* x is integral */
90 i>>=1;
91 if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-31));
93 SET_LDOUBLE_WORDS(x,se,i0,i1);
94 feholdexcept (&env);
95 w = TWO63[sx]+x;
96 t = w-TWO63[sx];
97 fesetenv (&env);
98 return t;
100 weak_alias (__nearbyintl, nearbyintl)