Back to index

glibc  2.9
s_modfl.c
Go to the documentation of this file.
00001 /* s_modfl.c -- long double version of s_modf.c.
00002  * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
00003  */
00004 
00005 /*
00006  * ====================================================
00007  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00008  *
00009  * Developed at SunPro, a Sun Microsystems, Inc. business.
00010  * Permission to use, copy, modify, and distribute this
00011  * software is freely granted, provided that this notice
00012  * is preserved.
00013  * ====================================================
00014  */
00015 
00016 #if defined(LIBM_SCCS) && !defined(lint)
00017 static char rcsid[] = "$NetBSD: $";
00018 #endif
00019 
00020 /*
00021  * modfl(long double x, long double *iptr)
00022  * return fraction part of x, and return x's integral part in *iptr.
00023  * Method:
00024  *     Bit twiddling.
00025  *
00026  * Exception:
00027  *     No exception.
00028  */
00029 
00030 #include "math.h"
00031 #include "math_private.h"
00032 
00033 #ifdef __STDC__
00034 static const long double one = 1.0;
00035 #else
00036 static long double one = 1.0;
00037 #endif
00038 
00039 #ifdef __STDC__
00040        long double __modfl(long double x, long double *iptr)
00041 #else
00042        long double __modfl(x, iptr)
00043        long double x,*iptr;
00044 #endif
00045 {
00046        int64_t i0,i1,j0;
00047        u_int64_t i;
00048        GET_LDOUBLE_WORDS64(i0,i1,x);
00049        j0 = ((i0>>48)&0x7fff)-0x3fff;     /* exponent of x */
00050        if(j0<48) {                 /* integer part in high x */
00051            if(j0<0) {                     /* |x|<1 */
00052               /* *iptr = +-0 */
00053                SET_LDOUBLE_WORDS64(*iptr,i0&0x8000000000000000ULL,0);
00054               return x;
00055            } else {
00056               i = (0x0000ffffffffffffLL)>>j0;
00057               if(((i0&i)|i1)==0) {        /* x is integral */
00058                   *iptr = x;
00059                   /* return +-0 */
00060                   SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0);
00061                   return x;
00062               } else {
00063                   SET_LDOUBLE_WORDS64(*iptr,i0&(~i),0);
00064                   return x - *iptr;
00065               }
00066            }
00067        } else if (j0>111) {        /* no fraction part */
00068            *iptr = x*one;
00069            /* We must handle NaNs separately.  */
00070            if (j0 == 0x4000 && ((i0 & 0x0000ffffffffffffLL) | i1))
00071              return x*one;
00072            /* return +-0 */
00073            SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0);
00074            return x;
00075        } else {                    /* fraction part in low x */
00076            i = -1ULL>>(j0-48);
00077            if((i1&i)==0) {         /* x is integral */
00078               *iptr = x;
00079               /* return +-0 */
00080               SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0);
00081               return x;
00082            } else {
00083               SET_LDOUBLE_WORDS64(*iptr,i0,i1&(~i));
00084               return x - *iptr;
00085            }
00086        }
00087 }
00088 weak_alias (__modfl, modfl)