1  /* s_modfl.c -- long double version of s_modf.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   * modfq(long double x, long double *iptr)
      22   * return fraction part of x, and return x's integral part in *iptr.
      23   * Method:
      24   *	Bit twiddling.
      25   *
      26   * Exception:
      27   *	No exception.
      28   */
      29  
      30  #include "quadmath-imp.h"
      31  
      32  static const __float128 one = 1.0;
      33  
      34  __float128 modfq(__float128 x, __float128 *iptr)
      35  {
      36  	int64_t i0,i1,j0;
      37  	uint64_t i;
      38  	GET_FLT128_WORDS64(i0,i1,x);
      39  	j0 = ((i0>>48)&0x7fff)-0x3fff;	/* exponent of x */
      40  	if(j0<48) {			/* integer part in high x */
      41  	    if(j0<0) {			/* |x|<1 */
      42  		/* *iptr = +-0 */
      43  	        SET_FLT128_WORDS64(*iptr,i0&0x8000000000000000ULL,0);
      44  		return x;
      45  	    } else {
      46  		i = (0x0000ffffffffffffLL)>>j0;
      47  		if(((i0&i)|i1)==0) {		/* x is integral */
      48  		    *iptr = x;
      49  		    /* return +-0 */
      50  		    SET_FLT128_WORDS64(x,i0&0x8000000000000000ULL,0);
      51  		    return x;
      52  		} else {
      53  		    SET_FLT128_WORDS64(*iptr,i0&(~i),0);
      54  		    return x - *iptr;
      55  		}
      56  	    }
      57  	} else if (j0>111) {		/* no fraction part */
      58  	    *iptr = x*one;
      59  	    /* We must handle NaNs separately.  */
      60  	    if (j0 == 0x4000 && ((i0 & 0x0000ffffffffffffLL) | i1))
      61  	      return x*one;
      62  	    /* return +-0 */
      63  	    SET_FLT128_WORDS64(x,i0&0x8000000000000000ULL,0);
      64  	    return x;
      65  	} else {			/* fraction part in low x */
      66  	    i = -1ULL>>(j0-48);
      67  	    if((i1&i)==0) { 		/* x is integral */
      68  		*iptr = x;
      69  		/* return +-0 */
      70  		SET_FLT128_WORDS64(x,i0&0x8000000000000000ULL,0);
      71  		return x;
      72  	    } else {
      73  		SET_FLT128_WORDS64(*iptr,i0,i1&(~i));
      74  		return x - *iptr;
      75  	    }
      76  	}
      77  }