#include <math.h>
#include <float.h>

/* Internal prototypes */
extern long double ldPerfectRootGuess( long double x );
extern long double ldPerfectRootImprove( long double ldX, long double ldK );
extern long double ldPerfectRoot( long double ldK );
extern long double ldRootSearch( long double ldK );

/*****************************************************************************/
/* Crude approximation for the perfect root of a number.  We want to be able */
/* to calculate what number, when raised to its own power, will be equal to  */
/* the input number.  For instance, if given 27, the answer will be 3, since */
/* 3 to the 3rd power is 27.  This function will provide an estimate for a   */
/* root solving function to begin iteration.                                 */
/*****************************************************************************/
long double ldPerfectRootGuess( long double x )
{
   long double f;
   /* First crude approximation */
   f = .136 * logl( x ) + 1.;

   /* Subtract (roughly) parabolic error term */
   x = f - ( 0.0002347645197348 * ( f - 574. ) * ( f - 574. ) - 77.088 );

   /* Subtract (roughly) sinusoidal error term */
   x += 21.73 * sinl( 0.01030030378226 * x );
   return x;
}

/********************************************************************/
/* This is Schroder's Method(A, Lamda=0.5, 1) for x^x - k = 0.      */
/* Reference: "On Infinitely Many Algorithms For Solving Equations" */
/*            by Ernst Schroder, Translated by G. W. Stewart.       */
/*            Pages 41 and 42 are the key passages for this work.   */
/* ---------------------------------------------------------------- */
/* Copyright (C) 1996 by Dann Corbit                                */
/********************************************************************/
long double ldPerfectRootImprove( long double ldX, long double ldK )
{
   long double ldNum; /* The numerator */
   long double ldDen; /* The denominator */
   const long double ldTerm = powl( ldX, ldX );
   const long double ldLamda = 0.5e0L;
   ldNum = ldX * ( ldTerm - ldK );
   ldDen = ldX * ( ldTerm * ( logl( ldX ) + 1.0e0L ) ) - ldLamda * ldNum;
   ldX -= ( ldNum / ldDen );
   return ldX;
}

/********************************************************************/
/* Calculate the solution in x for the equation:  x^x - k = 0.      */
/* Use whatever means are necessary. Report errors via _matherrl(). */
/* ---------------------------------------------------------------- */
/* Copyright (C) 1996 by Dann Corbit                                */
/********************************************************************/
long double ldPerfectRoot( long double ldK )
{
   int i;
   /* Primitive approximation */
   long double ldGuess = ldPerfectRootGuess( ldK );

   if ( ldK < 1e0L )
   {
      struct _exceptionl ldException;
      ldException.type = DOMAIN;
      ldException.name = "ldPerfectRoot";
      ldException.arg1 = ldK;
      ldException.arg2 = 0.0e0L;
      ldException.retval = -1.0e0L;
      _matherrl( &ldException );
   }
   if ( ldK > 1e+614L )
      return ldRootSearch( ldK );

   /* Make seven improvement passes */
   for ( i = 0; i < 7; i++ )
      ldGuess = ldPerfectRootImprove( ldGuess, ldK );

   /* Make (up to) seven improvement passes, if needed */
   for ( i = 0; i < 7; i++ )
      if ( fabsl( ldK - powl( ldGuess, ldGuess ) )/ldK >= LDBL_EPSILON*2.0e0L )
         ldGuess = ldPerfectRootImprove( ldGuess, ldK );
      else
         break;
   /* Return estimate of perfect root. */
   return ldGuess;
}
/********************************************************************/
/* Calculate the solution in x for the equation:  x^x - k = 0.      */
/* Use a binary search of the solution space to find the answer.    */
/* Use this routine when the input is too large for other methods.  */
/* ---------------------------------------------------------------- */
/* Copyright (C) 1996 by Dann Corbit                                */
/********************************************************************/
long double ldRootSearch( long double ldK )
{

   long double ldMax = 1500.0e0L;
   long double ldMin =  255.0e0L;
   long double ldGuess = ( ldMax + ldMin ) * 0.5e0L;
   int iBailout = 0;
   if ( ldK > 1.37053011714763899283e4764L )
   {
      struct _exceptionl ldException;
      ldException.type = DOMAIN;
      ldException.name = "ldRootSearch";
      ldException.arg1 = ldK;
      ldException.arg2 = 0.0e0L;
      ldException.retval = -1.0e0L;
      _matherrl( &ldException );
   }
   while ( ( fabsl( ldK - powl( ldGuess, ldGuess ) ) / ldK >= LDBL_EPSILON*2.0e0L ) && ( iBailout < 1000 )  && ( ldMax > ldMin ) )
   {
      if ( powl( ldGuess, ldGuess ) > ldK )
         ldMax = ldGuess;
      else
         ldMin = ldGuess;
      ldGuess = ( ldMax + ldMin ) * 0.5e0L;
      iBailout++;
   }
   return ldGuess;
}
#ifdef TEST
#include <stdio.h>
int main()
{

   long double ldK = 50.0e0L;
   long double ldGuess = ldPerfectRootGuess( ldK );
   printf( "0th guess is %30.20LeL for 50.\n", ldGuess );
   printf( "1st guess is %30.20LeL for 50.\n", ldGuess = ldPerfectRootImprove( ldGuess, ldK ) );
   printf( "2nd guess is %30.20LeL for 50.\n", ldGuess = ldPerfectRootImprove( ldGuess, ldK ) );
   printf( "3rd guess is %30.20LeL for 50.\n", ldGuess = ldPerfectRootImprove( ldGuess, ldK ) );
   printf( "4th guess is %30.20LeL for 50.\n", ldGuess = ldPerfectRootImprove( ldGuess, ldK ) );
   printf( "5th guess is %30.20LeL for 50.\n", ldGuess = ldPerfectRootImprove( ldGuess, ldK ) );
   printf( "6th guess is %30.20LeL for 50.\n", ldGuess = ldPerfectRootImprove( ldGuess, ldK ) );
   printf( "7th guess is %30.20LeL for 50.\n", ldGuess = ldPerfectRootImprove( ldGuess, ldK ) );

   for ( ldK = 1.; ldK < 1000; ldK += 3. )
   {
      ldGuess = ldPerfectRoot( ldK );
      printf( "N=%.0Lf, Guess=%.20LeL, Error=%.20LeL\n",
          ldK, ldGuess, fabsl( ldK - powl( ldGuess, ldGuess ) )/ldK );
   }
   for ( ldK = 1.; ldK < 1.0e1000L; ldK *= 2. )
   {
      ldGuess = ldPerfectRoot( ldK );
      printf( "N=%Le, Guess=%.20LeL, Error=%.20LeL\n",
          ldK, ldGuess, fabsl( ldK - powl( ldGuess, ldGuess ) )/ldK );
   }
   return 0;
}
#endif