12116Sjkh/* @(#)s_rint.c 5.1 93/09/24 */ 22116Sjkh/* 32116Sjkh * ==================================================== 42116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 52116Sjkh * 62116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business. 72116Sjkh * Permission to use, copy, modify, and distribute this 88870Srgrimes * software is freely granted, provided that this notice 92116Sjkh * is preserved. 102116Sjkh * ==================================================== 112116Sjkh */ 122116Sjkh 13176451Sdas#include <sys/cdefs.h> 14176451Sdas__FBSDID("$FreeBSD$"); 152116Sjkh 162116Sjkh/* 172116Sjkh * rint(x) 182116Sjkh * Return x rounded to integral value according to the prevailing 192116Sjkh * rounding mode. 202116Sjkh * Method: 212116Sjkh * Using floating addition. 222116Sjkh * Exception: 232116Sjkh * Inexact flag raised if x not equal to rint(x). 242116Sjkh */ 252116Sjkh 26175309Sdas#include <float.h> 27175309Sdas 282116Sjkh#include "math.h" 292116Sjkh#include "math_private.h" 302116Sjkh 31153043Sbdestatic const double 322116SjkhTWO52[2]={ 332116Sjkh 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ 342116Sjkh -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ 352116Sjkh}; 362116Sjkh 3797413Salfreddouble 38117912Speterrint(double x) 392116Sjkh{ 402116Sjkh int32_t i0,j0,sx; 412116Sjkh u_int32_t i,i1; 422116Sjkh double w,t; 432116Sjkh EXTRACT_WORDS(i0,i1,x); 442116Sjkh sx = (i0>>31)&1; 452116Sjkh j0 = ((i0>>20)&0x7ff)-0x3ff; 462116Sjkh if(j0<20) { 478870Srgrimes if(j0<0) { 482116Sjkh if(((i0&0x7fffffff)|i1)==0) return x; 492116Sjkh i1 |= (i0&0x0fffff); 502116Sjkh i0 &= 0xfffe0000; 512116Sjkh i0 |= ((i1|-i1)>>12)&0x80000; 522116Sjkh SET_HIGH_WORD(x,i0); 53175480Sbde STRICT_ASSIGN(double,w,TWO52[sx]+x); 542116Sjkh t = w-TWO52[sx]; 552116Sjkh GET_HIGH_WORD(i0,t); 562116Sjkh SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31)); 572116Sjkh return t; 582116Sjkh } else { 592116Sjkh i = (0x000fffff)>>j0; 602116Sjkh if(((i0&i)|i1)==0) return x; /* x is integral */ 612116Sjkh i>>=1; 622116Sjkh if(((i0&i)|i1)!=0) { 63153042Sbde /* 64153042Sbde * Some bit is set after the 0.5 bit. To avoid the 65153042Sbde * possibility of errors from double rounding in 66153042Sbde * w = TWO52[sx]+x, adjust the 0.25 bit to a lower 67153042Sbde * guard bit. We do this for all j0<=51. The 68153042Sbde * adjustment is trickiest for j0==18 and j0==19 69153042Sbde * since then it spans the word boundary. 70153042Sbde */ 712116Sjkh if(j0==19) i1 = 0x40000000; else 72153042Sbde if(j0==18) i1 = 0x80000000; else 732116Sjkh i0 = (i0&(~i))|((0x20000)>>j0); 742116Sjkh } 752116Sjkh } 762116Sjkh } else if (j0>51) { 772116Sjkh if(j0==0x400) return x+x; /* inf or NaN */ 782116Sjkh else return x; /* x is integral */ 792116Sjkh } else { 802116Sjkh i = ((u_int32_t)(0xffffffff))>>(j0-20); 812116Sjkh if((i1&i)==0) return x; /* x is integral */ 822116Sjkh i>>=1; 832116Sjkh if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20)); 842116Sjkh } 852116Sjkh INSERT_WORDS(x,i0,i1); 86175480Sbde STRICT_ASSIGN(double,w,TWO52[sx]+x); 872116Sjkh return w-TWO52[sx]; 882116Sjkh} 89175309Sdas 90175309Sdas#if (LDBL_MANT_DIG == 53) 91175309Sdas__weak_reference(rint, rintl); 92175309Sdas#endif 93