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 #include <math_ldbl_opt.h>
00033 
00034 #ifdef __STDC__
00035 static const long double one = 1.0;
00036 #else
00037 static long double one = 1.0;
00038 #endif
00039 
00040 #ifdef __STDC__
00041        long double __modfl(long double x, long double *iptr)
00042 #else
00043        long double __modfl(x, iptr)
00044        long double x,*iptr;
00045 #endif
00046 {
00047        int64_t i0,i1,j0;
00048        u_int64_t i;
00049        GET_LDOUBLE_WORDS64(i0,i1,x);
00050        i1 &= 0x000fffffffffffffLL;
00051        j0 = ((i0>>52)&0x7ff)-0x3ff;       /* exponent of x */
00052        if(j0<52) {                 /* integer part in high x */
00053            if(j0<0) {                     /* |x|<1 */
00054               /* *iptr = +-0 */
00055                SET_LDOUBLE_WORDS64(*iptr,i0&0x8000000000000000ULL,0);
00056               return x;
00057            } else {
00058               i = (0x000fffffffffffffLL)>>j0;
00059               if(((i0&i)|(i1&0x7fffffffffffffffLL))==0) {             /* x is integral */
00060                   *iptr = x;
00061                   /* return +-0 */
00062                   SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0);
00063                   return x;
00064               } else {
00065                   SET_LDOUBLE_WORDS64(*iptr,i0&(~i),0);
00066                   return x - *iptr;
00067               }
00068            }
00069        } else if (j0>103) {        /* no fraction part */
00070            *iptr = x*one;
00071            /* We must handle NaNs separately.  */
00072            if (j0 == 0x400 && ((i0 & 0x000fffffffffffffLL) | i1))
00073              return x*one;
00074            /* return +-0 */
00075            SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0);
00076            return x;
00077        } else {                    /* fraction part in low x */
00078            i = -1ULL>>(j0-52);
00079            if((i1&i)==0) {         /* x is integral */
00080               *iptr = x;
00081               /* return +-0 */
00082               SET_LDOUBLE_WORDS64(x,i0&0x8000000000000000ULL,0);
00083               return x;
00084            } else {
00085               SET_LDOUBLE_WORDS64(*iptr,i0,i1&(~i));
00086               return x - *iptr;
00087            }
00088        }
00089 }
00090 #ifdef IS_IN_libm
00091 long_double_symbol (libm, __modfl, modfl);
00092 #else
00093 long_double_symbol (libc, __modfl, modfl);
00094 #endif