Home | History | Annotate | Line # | Download | only in math 
      1 /* s_rintl.c -- long double version of s_rint.c.
      2  * Conversion to IEEE quad long double by Jakub Jelinek, jj (at) ultra.linux.cz.
      3  */
      4 
      5 /*
      6  * ====================================================
      7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
      8  *
      9  * Developed at SunPro, a Sun Microsystems, Inc. business.
     10  * Permission to use, copy, modify, and distribute this
     11  * software is freely granted, provided that this notice
     12  * is preserved.
     13  * ====================================================
     14  */
     15 
     16 #if defined(LIBM_SCCS) && !defined(lint)
     17 static char rcsid[] = "NetBSD: ";
     18 #endif
     19 
     20 /*
     21  * rintq(x)
     22  * Return x rounded to integral value according to the prevailing
     23  * rounding mode.
     24  * Method:
     25  *	Using floating addition.
     26  * Exception:
     27  *	Inexact flag raised if x not equal to rintq(x).
     28  */
     29 
     30 #define NO_MATH_REDIRECT
     31 
     32 #include "quadmath-imp.h"
     33 
     34 static const __float128
     35 TWO112[2]={
     36   5.19229685853482762853049632922009600E+33L, /* 0x406F000000000000, 0 */
     37  -5.19229685853482762853049632922009600E+33L  /* 0xC06F000000000000, 0 */
     38 };
     39 
     40 __float128 rintq(__float128 x)
     41 {
     42 	int64_t i0,j0,sx;
     43 	uint64_t i1 __attribute__ ((unused));
     44 	__float128 w,t;
     45 	GET_FLT128_WORDS64(i0,i1,x);
     46 	sx = (((uint64_t)i0)>>63);
     47 	j0 = ((i0>>48)&0x7fff)-0x3fff;
     48 	if(j0<112) {
     49 	    if(j0<0) {
     50 	        w = TWO112[sx]+x;
     51 	        t = w-TWO112[sx];
     52 		GET_FLT128_MSW64(i0,t);
     53 		SET_FLT128_MSW64(t,(i0&0x7fffffffffffffffLL)|(sx<<63));
     54 	        return t;
     55 	    }
     56 	} else {
     57 	    if(j0==0x4000) return x+x;	/* inf or NaN */
     58 	    else return x;		/* x is integral */
     59 	}
     60 	w = TWO112[sx]+x;
     61 	return w-TWO112[sx];
     62 }
     63