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