1  /* s_rintl.c -- long double version of s_rint.c.
       2   * Conversion to IEEE quad long double by Jakub Jelinek, jj@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  }