Home | History | Annotate | Line # | Download | only in math 
      1 /* s_modfl.c -- long double version of s_modf.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  * 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 }
     78