sysdeps/ieee754/ldbl-128ibm/s_frexpl.c

   1 /* s_frexpl.c -- long double version of s_frexp.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  * for non-zero x
  22  *      x = frexpl(arg,&exp);
  23  * return a long double fp quantity x such that 0.5 <= |x| <1.0
  24  * and the corresponding binary exponent "exp". That is
  25  *      arg = x*2^exp.
  26  * If arg is inf, 0.0, or NaN, then frexpl(arg,&exp) returns arg
  27  * with *exp=0.
  28  */
  29
  30 #include <math.h>
  31 #include <math_private.h>
  32 #include <math_ldbl_opt.h>
  33
  34 long double __frexpl(long double x, int *eptr)
  35 {
  36   uint64_t hx, lx, ix, ixl;
  37   int64_t explo, expon;
  38   double xhi, xlo;
  39
  40   ldbl_unpack (x, &xhi, &xlo);
  41   EXTRACT_WORDS64 (hx, xhi);
  42   EXTRACT_WORDS64 (lx, xlo);
  43   ixl = 0x7fffffffffffffffULL & lx;
  44   ix  = 0x7fffffffffffffffULL & hx;
  45   expon = 0;
  46   if (ix >= 0x7ff0000000000000ULL || ix == 0)
  47     {
  48       /* 0,inf,nan.  */
  49       *eptr = expon;
  50       return x + x;
  51     }
  52   expon = ix >> 52;
  53   if (expon == 0)
  54     {
  55       /* Denormal high double, the low double must be 0.0.  */
  56       int cnt;
  57
  58       /* Normalize.  */
  59       if (sizeof (ix) == sizeof (long))
  60         cnt = __builtin_clzl (ix);
  61       else if ((ix >> 32) != 0)
  62         cnt = __builtin_clzl ((long) (ix >> 32));
  63       else
  64         cnt = __builtin_clzl ((long) ix) + 32;
  65       cnt = cnt - 12;
  66       expon -= cnt;
  67       ix <<= cnt + 1;
  68     }
  69   expon -= 1022;
  70   ix &= 0x000fffffffffffffULL;
  71   hx &= 0x8000000000000000ULL;
  72   hx |= (1022LL << 52) | ix;
  73
  74   if (ixl != 0)
  75     {
  76       /* If the high double is an exact power of two and the low
  77          double has the opposite sign, then the exponent calculated
  78          from the high double is one too big.  */
  79       if (ix == 0
  80           && (int64_t) (hx ^ lx) < 0)
  81         {
  82           hx += 1LL << 52;
  83           expon -= 1;
  84         }
  85
  86       explo = ixl >> 52;
  87       if (explo == 0)
  88         {
  89           /* The low double started out as a denormal.  Normalize its
  90              mantissa and adjust the exponent.  */
  91           int cnt;
  92
  93           if (sizeof (ixl) == sizeof (long))
  94             cnt = __builtin_clzl (ixl);
  95           else if ((ixl >> 32) != 0)
  96             cnt = __builtin_clzl ((long) (ixl >> 32));
  97           else
  98             cnt = __builtin_clzl ((long) ixl) + 32;
  99           cnt = cnt - 12;
 100           explo -= cnt;
 101           ixl <<= cnt + 1;
 102         }
 103
 104       /* With variable precision we can't assume much about the
 105          magnitude of the returned low double.  It may even be a
 106          denormal.  */
 107       explo -= expon;
 108       ixl &= 0x000fffffffffffffULL;
 109       lx  &= 0x8000000000000000ULL;
 110       if (explo <= 0)
 111         {
 112           /* Handle denormal low double.  */
 113           if (explo > -52)
 114             {
 115               ixl |= 1LL << 52;
 116               ixl >>= 1 - explo;
 117             }
 118           else
 119             {
 120               ixl = 0;
 121               lx = 0;
 122               if ((hx & 0x7ff0000000000000ULL) == (1023LL << 52))
 123                 {
 124                   /* Oops, the adjustment we made above for values a
 125                      little smaller than powers of two turned out to
 126                      be wrong since the returned low double will be
 127                      zero.  This can happen if the input was
 128                      something weird like 0x1p1000 - 0x1p-1000.  */
 129                   hx -= 1LL << 52;
 130                   expon += 1;
 131                 }
 132             }
 133           explo = 0;
 134         }
 135       lx |= (explo << 52) | ixl;
 136     }
 137
 138   INSERT_WORDS64 (xhi, hx);
 139   INSERT_WORDS64 (xlo, lx);
 140   x = ldbl_pack (xhi, xlo);
 141   *eptr = expon;
 142   return x;
 143 }
 144 #if IS_IN (libm)
 145 long_double_symbol (libm, __frexpl, frexpl);
 146 #else
 147 long_double_symbol (libc, __frexpl, frexpl);
 148 #endif