Listing 6. C code to determine optimal scale and offset for some approximation

const double Epsilon = 0.0001;
const double ThresholdError = 0.0001;

// Find values for A and B minimizing the maximal error given n and m
void FindAB( double n, double m, double &A, double &B )
{
   double k0 = 0;
   double k1 = 1;
   double k;
   double ErrorLeft, ErrorRight;

   // Binary search for optimal crossing point
   do
   {
      k = (k0 + k1)/2.0;
      A = ( pow( k, n/m ) - 1.0 ) / ( k - 1.0 );
      B = 1.0 - A;

      EvaluateMaxError( n, m, A, B, k, ErrorLeft, ErrorRight );

      if( ErrorLeft < ErrorRight )
         k0 = k;
      else
         k1 = k;

   } while( fabs( ErrorLeft - ErrorRight ) > ThresholdError );
}

// Evaluate the maximal error to the left and right of the crossing point
// for given values of A and B
// The crossing point k is given to reduce computation
void EvaluateMaxError( double n, double m, double A, double B, double k,
double &ErrorLeft, double &ErrorRight )
{
   double x;

   double DerErr;

   // Find maximal point of the error function on the left of the crossing point
   double x0 = 0;
   double x1 = k;
   do
   {
      x = (x0 + x1)/2.0;

      DerErr = DerError( x, n, m, A, B );
      if( DerErr < 0 )
         x0 = x;
      else
         x1 = x;
   } while( fabs( x0-x1 ) > Epsilon );

   // Evaluate the error at that point
   ErrorLeft = fabs( Error( x, n, m, A, B ) );

   x0 = k;
   x1 = 1;

   // Find maximal point of the error function on the right of the crossing point
   do
   {
      x = (x0 + x1)/2.0;

      DerErr = DerError( x, n, m, A, B );
      if( DerErr > 0 )
         x0 = x;
      else
         x1 = x;
   } while( fabs( x0-x1 ) > Epsilon );

   // Evaluate the error at that point
   ErrorRight = fabs( Error( x, n, m, A, B ) );

}

// Evaluate the approximation function
double Approx( double x, double m, double A, double B )
{
   if( x < -B/A )
      return 0;
   return pow( A * x + B, m );
}

// Evaluate the derivative of the approximation function
double DerApprox( double x, double m, double A, double B )
{
   if( x < -B/A )
      return 0;
   return A*m*pow( A * x + B, m-1 );
}

// Evaluate the error function
double Error( double x, double n, double m, double A, double B )
{
   return Approx( x, m, A, B ) - pow( x, n );
}

// Evaluate the derivative of the error function
double DerError( double x, double n, double m, double A, double B )
{
   return DerApprox( x, m, A, B ) - n*pow( x, n-1 );
}