#include <math.h>
#include <float.h>
/*****************************************************************************/
/* 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 (Two or three correct digits).   */
/*****************************************************************************/
long double PerfectRootGuess( long double x )
{
   long double Guess=0;
   long double Test;

   /* Min is (1/e)^(1/e).  For smaller x, there are two distinct answers. */
   long double Min = 0.69220062755534635386542199718278976149L;

   /* Choose Max such that Max^Max <= LDBL_MAX */
   long double Max = 1546.0e0L;
   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, Omega=4 ) for x^x - c = 0. */
/* Reference: "On Infinitely Many Algorithms For Solving Equations" */
/*            by Ernst Schroder, Translated by G. W. Stewart.       */
/* ---------------------------------------------------------------- */
/* Copyright (C) 1996 by Dann Corbit                                */
/********************************************************************/
long double PerfectRootImprove( long double c, long double x )
{
   long double t1 = x * x;
   long double t2 = t1 * x;
   long double t3 = logl( x );
   long double t4 = t3 * t3;
   long double t5 = t4 * t3;
   long double t6 = t2 * t5;
   long double t7 = t1 * t1;
   long double t8 = powl( t7, x );
   long double t10 = t8 * t1;
   long double t11 = t10 * t3;
   long double t12 = powl( t1, 2.0 * x );
   long double t13 = t12 * t7;
   long double t14 = t4 * t4;
   long double t15 = t13 * t14;
   long double t16 = t13 * t4;
   long double t18 = powl( x, x + 4.0 );
   long double t20 = c * c;
   long double t21 = t20 * c;
   long double t22 = t18 * t4 * t21;
   long double t23 = t12 * t2;
   long double t24 = t23 * t4;
   long double t25 = t7 * t8;
   long double t26 = powl( t1, x );
   long double t27 = t7 * t26;
   long double t28 = t3 * t20;
   long double t29 = t27 * t28;
   long double t30 = t4 * t20;
   long double t31 = t27 * t30;
   long double t32 = t7 * t4;
   long double t33 = powl( t2, x );
   long double t34 = t33 * c;
   long double t35 = t32 * t34;
   long double t36 = t32 * t8;
   long double t37 = t7 * t14;
   long double t38 = t37 * t8;
   long double t40 = t27 * t14 * t20;
   long double t41 =  28.0 * t6 * t8 + 8.0 * t11 - 14.0 * t15 - 84.0 * t16
                     - 6.0 * t22 - 108.0 * t24 + 13.0 * t25 - 44.0 * t29
                    - 66.0 * t31 - 66.0 * t35 + 78.0 * t36 + 13.0 * t38
                    - 11.0 * t40 - 14.0 * t13;
   long double t42 = t7 * t5;
   long double t43 = t42 * t8;
   long double t44 = t12 * t1;
   long double t45 = t7 * t3;
   long double t46 = t45 * t34;
   long double t47 = t2 * t26;
   long double t48 = t47 * t30;
   long double t50 = powl( x, x + 3.0 );
   long double t52 = t50 * t4 * t21;
   long double t54 = t50 * t3 * t21;
   long double t56 = t2 * t3 * t34;
   long double t58 = t27 * t5 * t20;
   long double t60 = powl( x, x + 1.0 );
   long double t61 = t60 * t21;
   long double t64 = t26 * t3 * t1 * t20;
   long double t65 = t42 * t34;
   long double t66 = t27 * t20;
   long double t68 = t26 * t1 * t20;
   long double t69 =  - 60.0 * t23 + 52.0 * t43 - 22.0 * t44 - 44.0 * t46
                      - 54.0 * t48 + 11.0 * t10 - 18.0 * t52 - 24.0 * t54
                      + 72.0 * t56 - 44.0 * t58 + 2.0 * t61 + 8.0 * t64
                      - 44.0 * t65 - 11.0 * t66 + 11.0 * t68;
   long double t71 = t8 * t2;
   long double t72 = t71 * t3;
   long double t74 = t2 * t4 * t34;
   long double t75 = t47 * t28;
   long double t76 = t37 * t34;
   long double t78 = t7 * t33 * c;
   long double t79 = t8 * x;
   long double t80 = t33 * t1;
   long double t82 = t80 * t3 * c;
   long double t83 = t50 * t21;
   long double t84 = t33 * t2;
   long double t88 = powl( x, x + 2.0 );
   long double t89 = t88 * t21;
   long double t93 = t18 * t3 * t21;
   long double t94 = t18 * t21;
   long double t96 = t18 * t14 * t21;
   long double t97 = 168.0 * t72 + 54.0 * t74 - 72.0 * t75 - 11.0 * t76
                    - 11.0 * t78 - 2.0 * t79 + 8.0 * t82 - 10.0 * t83
                    + 12.0 * t84 * t5 * c - 11.0 * t89 + 70.0 * t71
                    - 12.0 * t6 * t26 * t20 - 4.0 * t93 - t94 - t96;
   long double t99 = t88 * t3 * t21;
   long double t100 = t45 * t8;
   long double t105 = t18 * t5 * t21;
   long double t106 = t13 * t3;
   long double t107 = t13 * t5;
   long double t108 = t23 * t3;
   long double t109 = t47 * t20;
   long double t110 = t84 * c;
   long double t112 = t12 * t3 * t1;
   long double t114 = t33 * x * c;
   long double t116 = t26 * x * t20;
   long double t117 = t71 * t4;
   long double t118 = t80 * c;
   long double t119 =  - 8.0 * t99 + 52.0 * t100 - 24.0 * t23 * t5
                       - 4.0 * t50 * t5 * t21 - 4.0 * t105 - 56.0 * t106
                       - 56.0 * t107 - 144.0 * t108 - 30.0 * t109
                       + 30.0 * t110 - 16.0 * t112 + 6.0 * t114 - 6.0 * t116
                       + 126.0 * t117 + 11.0 * t118;
   long double t123 =  - 4.0 * t11 - 14.0 * t15 - 84.0 * t16 - 6.0 * t22
                       - 36.0 * t24 + 13.0 * t25 - 44.0 * t29 - 66.0 * t31
                       - 66.0 * t35 + 78.0 * t36 + 13.0 * t38 - 11.0 * t40
                       - 14.0 * t13;
   long double t124 =  - 36.0 * t23 + 52.0 * t43 + 2.0 * t44 - 44.0 * t46
                       - 18.0 * t48 - t10 - 6.0 * t52 - 12.0 * t54
                       + 36.0 * t56 - 44.0 * t58 - 2.0 * t61 - 4.0 * t64
                       - 44.0 * t65 - 11.0 * t66;
   long double t126 =  - t68 + 84.0 * t72 + 18.0 * t74 - 36.0 * t75
                       - 11.0 * t76 - 11.0 * t78 + 2.0 * t79 - 4.0 * t82
                       - 6.0 * t83 + t89 + 42.0 * t71 - 4.0 * t93 - t94;
   long double t127 =  - t96 + 4.0 * t99 + 52.0 * t100 - 4.0 * t105
                       - 56.0 * t106 - 56.0 * t107 - 72.0 * t108 - 18.0 * t109
                       + 18.0 * t110 + 8.0 * t112 - 6.0 * t114 + 6.0 * t116
                       + 42.0 * t117 - t118;
   return x * ( t41 + t69 + t97 + t119 ) / ( t123 + t124 + t126 + t127 );
}
/********************************************************************/
/* This is Schroder's Method A(Lamda=1/2, Omega=1) for x^x - c = 0. */
/* Reference: "On Infinitely Many Algorithms For Solving Equations" */
/*            by Ernst Schroder, Translated by G. W. Stewart.       */
/* ---------------------------------------------------------------- */
/* Copyright (C) 1996 by Dann Corbit                                */
/********************************************************************/
long double PerfectRootImproveSlow( long double c, long double x )
{
   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 - c );
   Den = x * ( Term * ( logl( x ) + 1.0e0L ) ) - Lamda * Num;
   x -= ( Num / Den );
   return x;
}
long double PerfectRoot( long double c )
{
   int i;
   long double x = PerfectRootImproveSlow( c, PerfectRootGuess( c ) );
   for ( i = 0; i < 20; i++ )
      if ( fabsl( c - powl( x, x ) ) / c >= LDBL_EPSILON*2.0e0L )
      {
         if ( c < 2.0e37 ) /* Then we can use the faster method... */
            x = PerfectRootImprove( c, x );
         else /* Much slower, but less prone to overflow. */
            x = PerfectRootImproveSlow( c, x );
      }
      else
         break;
   return x;
}
#ifdef TEST
#include <stdio.h>
int main()
{

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

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

   return 0;
}
#endif