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

/* Internal prototypes */
extern long double PerfectRootGuess( long double x );
extern long double PerfectRootImprove( long double X, long double K );
extern long double PerfectRoot( long double K );
extern long double RootSearchBig( long double K );
extern long double RootSearchSmall( long double K );

/*****************************************************************************/
/* 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 PerfectRootGuess( long double x )
{
   long double Guess;
   long double Test;
   long double Min = 0.69220062755534635386542199718278976149L;
   long double Max = 1546.0e0L; /* Choose so that Max^Max = LDBL_MAX */
   int n;
   for ( n = 0; n < 16; n++ )
   {
      Guess = ( Min + Max ) * 0.5e0L;
      Test = powl( Guess, Guess );
      if ( Test > x )
         Max = Guess;
      else
         Min = Guess;
   }
   return Guess;
}
/********************************************************************/
/* 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 PerfectRootImprove( long double X, long double K )
{
   long double Num; /* The numerator */
   long double Den; /* The denominator */
   const long double Term = powl( X, X );
   const long double Lamda = 0.5e0L;
   Num = X * ( Term - K );
   Den = X * ( Term * ( logl( X ) + 1.0e0L ) ) - Lamda * Num;
   X -= ( Num / Den );
   return X;
}

/********************************************************************/
/* 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 PerfectRoot( long double K )
{
   int i;
   /* Primitive approximation */
   long double Guess = PerfectRootGuess( K );
   /* Number .6922... is (1/e)^(1/e).  This is the minimum possible. */
   if ( K < 1 )
      return RootSearchSmall( K );

   if ( K > 1e+614L )
      return RootSearchBig( K );

   /* Make seven improvement passes */
   for ( i = 0; i < 7; i++ )
      Guess = PerfectRootImprove( Guess, K );

   /* Make (up to) seven improvement passes, if needed */
   for ( i = 0; i < 7; i++ )
      if ( fabsl( K - powl( Guess, Guess ) )/K >= LDBL_EPSILON )
         Guess = PerfectRootImprove( Guess, K );
      else
         break;
   /* Return estimate of perfect root. */
   return Guess;
}
/********************************************************************/
/* 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 RootSearchBig( long double K )
{

   long double Max = 1500.0e0L;
   long double Min =  255.0e0L;
   long double Guess = ( Max + Min ) * 0.5e0L;
   int iBailout = 0;
   if ( K > 1.37053011714763899283e4764L )
   {
      struct _exceptionl Exception;
      Exception.type = DOMAIN;
      Exception.name = "RootSearchBig";
      Exception.arg1 = K;
      Exception.arg2 = 0.0e0L;
      Exception.retval = -1.0e0L;
      _matherrl( &Exception );
      return -1;
   }
   while ( ( fabsl( K - powl( Guess, Guess ) ) / K >= LDBL_EPSILON ) && ( iBailout < 1000 )  && ( Max > Min ) )
   {
      if ( powl( Guess, Guess ) > K )
         Max = Guess;
      else
         Min = Guess;
      Guess = ( Max + Min ) * 0.5e0L;
      iBailout++;
   }
   return Guess;
}
/********************************************************************/
/* 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 small for other methods.  */
/* ---------------------------------------------------------------- */
/* Copyright (C) 1996 by Dann Corbit                                */
/********************************************************************/
long double RootSearchSmall( long double K )
{

   long double Max = 1.0e0L;
   /* Min is 1/e, which is the minimum value for x^x */
   long double Min = 0.3678794411714423215955e0L;
   long double Guess = ( Max + Min ) * 0.5e0L;
   int iBailout = 0;
   if ( K < 0.692200627555346353865422e0L )
   {
      struct _exceptionl Exception;
      Exception.type = DOMAIN;
      Exception.name = "RootSearchSmall";
      Exception.arg1 = K;
      Exception.arg2 = 0.0e0L;
      Exception.retval = -1.0e0L;
      _matherrl( &Exception );
      return -1;
   }
   while ( ( fabsl( K - powl( Guess, Guess ) ) / K >= LDBL_EPSILON ) && ( iBailout < 1000 )  && ( Max > Min ) )
   {
      if ( powl( Guess, Guess ) > K )
         Max = Guess;
      else
         Min = Guess;
      Guess = ( Max + Min ) * 0.5e0L;
      iBailout++;
   }
   return Guess;
}
#ifdef TEST
#include <stdio.h>
int main()
{

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

   for ( K = 1.; K < 1000; K += 3. )
   {
      Guess = PerfectRoot( K );
      printf( "N=%.0Lf, Guess=%.20LeL, Error=%.20LeL\n",
          K, Guess, fabsl( K - powl( Guess, Guess ) )/K );
   }
   for ( K = 1.; K < 1.0e1000L; K *= 2. )
   {
      Guess = PerfectRoot( K );
      printf( "N=%Le, Guess=%.20LeL, Error=%.20LeL\n",
          K, Guess, fabsl( K - powl( Guess, Guess ) )/K );
   }
   /* Huge: */
   K = 1.0e1300L;
   Guess = PerfectRoot( K );
   printf( "N=%Le, Guess=%.20LeL, Error=%.20LeL\n",
       K, Guess, fabsl( K - powl( Guess, Guess ) )/K );
   /* Tiny: */
   K = 0.75L;
   Guess = PerfectRoot( K );
   printf( "N=%Le, Guess=%.20LeL, Error=%.20LeL\n",
       K, Guess, fabsl( K - powl( Guess, Guess ) )/K );

   /* Huge Error: */
   K = LDBL_MAX;
   Guess = PerfectRoot( K );
   printf( "N=%Le, Guess=%.20LeL, Error=%.20LeL\n",
       K, Guess, fabsl( K - powl( Guess, Guess ) )/K );
   /* Tiny: */
   K = 0.2L;
   Guess = PerfectRoot( K );
   printf( "N=%Le, Guess=%.20LeL, Error=%.20LeL\n",
       K, Guess, fabsl( K - powl( Guess, Guess ) )/K );

   return 0;
}
#endif