/*
	Copyright (C) 2003, 2004 by J. David Sexton.
	Some rights NOT reserved!

	This source code should compile with any C compiler that meets the
	ISO 9899:1990 ("C89") standard.

	This source code is not public domain or Open Source, but I hereby license
	you to copy it and give it to others so long as you don't profit by doing
	so. In addition, you may include these functions in your own software, so
	long as that software is not allowed to be sold or used for profit.
	If I sold this, I would guarantee it; but as I don't, I won't.
*/

#include <float.h>
#include <errno.h>

#ifndef HUGE_VALL
#define HUGE_VALL (LDBL_MAX / LDBL_MIN)
#endif
#ifndef NAN
#define NAN (0.0L / 0.0L)
#endif

long double Log2OfX (long double theX, long double *theInt);
long double TwoToTheX (long double theX);

/*
	Log2OfX returns the mantissa of the base 2 logarithm of theX.
	The integer part of the logarithm is returned in the variable
	pointed to by theInt; this will be the largest integer not
	greater than the logarithm. E.g., if the logarithm is 2.1,
	Log2OfX will return 0.1 and 2.0; if the logarithm is -2.1,
	it will return 0.9 and -3.0.
	If FLT_RADIX is 2, the mantissa returned should be accurate to
	LDBL_MANT_DIG bits. The value returned in the variable pointed
	to by theInt should be precise.
*/

long double
Log2OfX (long double theX, long double *theInt)
{
  long double theSum, lastSum, pSum, qSum;
  long double ldX, ldY, ldZ, theIndex;

  if (theX < 0.0L)
    {
      errno = EDOM;
      *theInt = NAN;
      return (0.0L);
    }
  if (theX == 0.0L)
    {
      errno = ERANGE;
      *theInt = -HUGE_VALL;
      return (0.0L);
    }
  if (theX == HUGE_VALL)
    {
      errno = ERANGE;
      *theInt = HUGE_VALL;
      return (0.0L);
    }
  ldX = theX;
  theSum = pSum = 0.0L;
  ldZ = 1.0L;
  while (ldX / ldZ < 1.0L)
    {
      theIndex = 1.0L;
      ldY = 2.0L;
      do
	{
	  theSum -= theIndex;
	  ldZ /= ldY;
	  theIndex += theIndex;
	  ldY *= ldY;
	}
      while (ldZ / ldX >= ldY);
    }
  while (ldX / ldZ >= 2.0L)
    {
      theIndex = 1.0L;
      ldY = 2.0L;
      do
	{
	  theSum += theIndex;
	  ldZ *= ldY;
	  theIndex += theIndex;
	  ldY *= ldY;
	}
      while (ldX / ldZ >= ldY);
    }
  ldX = (ldX - ldZ) / ldZ;
  if (ldX > 0.0L)
    {
      theIndex = 1.0L;
      do
	{
	  theIndex /= 2.0L;
	  ldY = ldX * ldX;
	  ldZ = ldX + ldX;
	  ldX = ldY + ldZ;
	}
      while (ldX < 1.0L);
      ldX = (ldY + --ldZ) / 2.0L;
      qSum = theIndex;
      if (ldX > 0.0L)
	{
	  do
	    {
	      lastSum = pSum;
	      do
		{
		  theIndex /= 2.0L;
		  ldY = ldX * ldX;
		  ldZ = ldX + ldX;
		  ldX = ldY + ldZ;
		}
	      while (ldX < 1.0L);
	      ldX = (ldY + --ldZ) / 2.0L;
	      pSum += theIndex;
	    }
	  while (ldX > 0.0L && pSum != lastSum);
	}
      pSum += qSum;
    }
  *theInt = theSum;
  return (pSum);
}

/*
	TwoToTheX returns 2 to the power of a long double argument.
	If FLT_RADIX is 2, the value returned should be accurate to
	LDBL_MANT_DIG bits.
*/

long double
TwoToTheX (long double theX)
{
  long double ldX, theIndex, theProd, pProd;
  long double thePower, prevPower, lastApprox;

  if (theX == -HUGE_VALL)
    return (errno = ERANGE, 0.0L);
  if (theX == HUGE_VALL)
    return (errno = ERANGE, HUGE_VALL);
  ldX = theX;
  theProd = 1.0L;
  while (ldX < 0.0L)
    {
      theIndex = 1.0L;
      thePower = 2.0L;
      do
	{
	  theProd /= thePower;
	  ldX += theIndex;
	  if (theProd == 0.0L)
	    return (errno = ERANGE, theProd);
	  theIndex += theIndex;
	  thePower *= thePower;
	}
      while (-ldX >= theIndex);
    }
  while (ldX >= 1.0L)
    {
      theIndex = 1.0L;
      thePower = 2.0L;
      do
	{
	  theProd *= thePower;
	  ldX -= theIndex;
	  if (theProd == HUGE_VALL)
	    return (errno = ERANGE, theProd);
	  theIndex += theIndex;
	  thePower *= thePower;
	}
      while (ldX >= theIndex);
    }
  if (ldX > 0.0L)
    {
      pProd = 0.0L;
      thePower = theIndex = 1.0L;
      do
	{
	  prevPower = thePower;
	  thePower = 0.0L;
	  do
	    {
	      lastApprox = thePower;
	      thePower += (prevPower - thePower) / (thePower + 1.0L);
	      thePower /= 2.0L;
	    }
	  while (thePower != lastApprox);
	  theIndex /= 2.0L;
	  if (ldX >= theIndex)
	    {
	      pProd += (pProd * thePower) + thePower;
	      ldX -= theIndex;
	    }
	}
      while (ldX > 0.0L && (thePower + 1.0L) > 1.0L);
      theProd += theProd * pProd;
      if (theProd == HUGE_VALL)
	errno = ERANGE;
    }
  return (theProd);
}