1 2 /* @(#)s_rint.c 1.3 95/01/18 */ 3 /* 4 * ==================================================== 5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6 * 7 * Developed at SunSoft, a Sun Microsystems, Inc. business. 8 * Permission to use, copy, modify, and distribute this 9 * software is freely granted, provided that this notice 10 * is preserved. 11 * ==================================================== 12 */ 13 14 /* 15 * ieee_rint(x) 16 * Return x rounded to integral value according to the prevailing 17 * rounding mode. 18 * Method: 19 * Using floating addition. 20 * Exception: 21 * Inexact flag raised if x not equal to ieee_rint(x). 22 */ 23 24 #include "fdlibm.h" 25 26 #ifdef __STDC__ 27 static const double 28 #else 29 static double 30 #endif 31 TWO52[2]={ 32 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ 33 -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ 34 }; 35 36 #ifdef __STDC__ ieee_rint(double x)37 double ieee_rint(double x) 38 #else 39 double ieee_rint(x) 40 double x; 41 #endif 42 { 43 int i0,j0,sx; 44 unsigned i,i1; 45 double w,t; 46 i0 = __HI(x); 47 sx = (i0>>31)&1; 48 i1 = __LO(x); 49 j0 = ((i0>>20)&0x7ff)-0x3ff; 50 if(j0<20) { 51 if(j0<0) { 52 if(((i0&0x7fffffff)|i1)==0) return x; 53 i1 |= (i0&0x0fffff); 54 i0 &= 0xfffe0000; 55 i0 |= ((i1|-i1)>>12)&0x80000; 56 __HI(x)=i0; 57 w = TWO52[sx]+x; 58 t = w-TWO52[sx]; 59 i0 = __HI(t); 60 __HI(t) = (i0&0x7fffffff)|(sx<<31); 61 return t; 62 } else { 63 i = (0x000fffff)>>j0; 64 if(((i0&i)|i1)==0) return x; /* x is integral */ 65 i>>=1; 66 if(((i0&i)|i1)!=0) { 67 if(j0==19) i1 = 0x40000000; else 68 i0 = (i0&(~i))|((0x20000)>>j0); 69 } 70 } 71 } else if (j0>51) { 72 if(j0==0x400) return x+x; /* inf or NaN */ 73 else return x; /* x is integral */ 74 } else { 75 i = ((unsigned)(0xffffffff))>>(j0-20); 76 if((i1&i)==0) return x; /* x is integral */ 77 i>>=1; 78 if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20)); 79 } 80 __HI(x) = i0; 81 __LO(x) = i1; 82 w = TWO52[sx]+x; 83 return w-TWO52[sx]; 84 } 85