Back to index

glibc  2.9
s_frexpl.c
Go to the documentation of this file.
00001 /* s_frexpl.c -- long double version of s_frexp.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  * for non-zero x
00022  *     x = frexpl(arg,&exp);
00023  * return a long double fp quantity x such that 0.5 <= |x| <1.0
00024  * and the corresponding binary exponent "exp". That is
00025  *     arg = x*2^exp.
00026  * If arg is inf, 0.0, or NaN, then frexpl(arg,&exp) returns arg
00027  * with *exp=0.
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
00036 #else
00037 static long double
00038 #endif
00039 two107 = 162259276829213363391578010288128.0; /* 0x4670000000000000, 0 */
00040 
00041 #ifdef __STDC__
00042        long double __frexpl(long double x, int *eptr)
00043 #else
00044        long double __frexpl(x, eptr)
00045        long double x; int *eptr;
00046 #endif
00047 {
00048        u_int64_t hx, lx, ix, ixl;
00049        int64_t explo;
00050        GET_LDOUBLE_WORDS64(hx,lx,x);
00051        ixl = 0x7fffffffffffffffULL&lx;
00052        ix =  0x7fffffffffffffffULL&hx;
00053        *eptr = 0;
00054        if(ix>=0x7ff0000000000000ULL||((ix|ixl)==0)) return x;  /* 0,inf,nan */
00055        if (ix<0x0010000000000000ULL) {           /* subnormal */
00056            x *= two107;
00057            GET_LDOUBLE_MSW64(hx,x);
00058            ix = hx&0x7fffffffffffffffULL;
00059            *eptr = -107;
00060        }
00061        *eptr += (ix>>52)-1022;
00062 
00063        if (ixl != 0ULL) {
00064          explo = (ixl>>52) - (ix>>52) + 0x3fe;
00065          if ((ixl&0x7ff0000000000000ULL) == 0LL) {
00066            /* the lower double is a denomal so we need to correct its
00067               mantissa and perhaps its exponent.  */
00068            int cnt;
00069 
00070            if (sizeof (ixl) == sizeof (long))
00071              cnt = __builtin_clzl (ixl);
00072            else if ((ixl >> 32) != 0)
00073              cnt = __builtin_clzl ((long) (ixl >> 32));
00074            else
00075              cnt = __builtin_clzl ((long) ixl) + 32;
00076            cnt = cnt - 12;
00077            lx = (lx&0x8000000000000000ULL) | ((explo-cnt)<<52)
00078               | ((ixl<<(cnt+1))&0x000fffffffffffffULL);
00079          } else
00080            lx = (lx&0x800fffffffffffffULL) | (explo<<52);
00081        } else
00082          lx = 0ULL;
00083 
00084        hx = (hx&0x800fffffffffffffULL) | 0x3fe0000000000000ULL;
00085        SET_LDOUBLE_WORDS64(x,hx,lx);
00086        return x;
00087 }
00088 #ifdef IS_IN_libm
00089 long_double_symbol (libm, __frexpl, frexpl);
00090 #else
00091 long_double_symbol (libc, __frexpl, frexpl);
00092 #endif