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