code: public static double cubicSafe( double f0 , double g0 , double x0 ,
  									double f1 , double g1 , double x1 ,
  									double min , double max ) {
  		double theta = 3.0*(f0-f1)/(x1-x0) + g0 + g1;
  		double s = Math.max(Math.abs(theta),Math.abs(g0));
  		s= Math.max(s,Math.abs(g1));
  		double gamma = s*Math.sqrt((theta/s)*(theta/s) - (g0/s)*(g1/s));
  		if( x1 < x0 )
  			gamma = -gamma;
  		double p = (gamma-g0) + theta;
  		double q = (gamma+(g1-g0))+gamma;
  
  		double r = p/q;
  
  		// gamma == 0 only rises if the cubic does not tend to infinity in the direction of the step
  		if( r < 0 && gamma != 0 ) {
  			return x0 + r*(x1-x0);
  		} else if( x0 > x1 ) {
  			return max;
  		} else {
  			return min;
  		}
  
  	}
docs: <p>
  Use cubic interpolation only if the cubic tends to infinity in the direction of the step or if the minim of the
  cubic is beyond x1.  Otherwise the the step will be max if x0 {@code >} x1 else it will be min.
  </p>
  <p>
  [1] MINPACK-2 source code http://ftp.mcs.anl.gov/pub/MINPACK-2/dcstep.f
  </p>
  
  @param f0 Value of f(x0)
  @param g0 Derivative f'(x0)
  @param x0 First sample point
  @param f1 Value of f(x1)
  @param g1 Derivative g'(x1)
  @param x1 Second sample point
  @return Interpolated point