#!/usr/bin/perl

# root($func, $lo_guess, $hi_guess, $epsilon) uses Newton's method
# to find the root of the function $func (provided as a code reference).
# It searches between $lo_guess and $hi_guess, and returns as soon as
# the root is known within $epsilon.
#
sub root {
    my ($func, $lo_guess, $hi_guess, $epsilon) = @_;
    my ($lo, $hi, $root, $step, $step_old, $value, $i);
    my ($lo_val, $hi_val) =
           ( &{$func}( $lo_guess ), &{$func}( $hi_guess ) );
    my $deriv = deriv($func, $hi_guess);

    use constant ITERATIONS => 128;

    return undef if $lo_val > 0 && $hi_val > 0 or
                        $lo_val < 0 && $hi_val < 0;

    # Are we already at a root?
    #
    return $lo_guess if abs($lo_val) < $epsilon;
    return $hi_guess if abs($hi_val) < $epsilon;

    if ($lo_val < 0) { ($lo, $hi) = ($lo_guess, $hi_guess) }
    else             { ($lo, $hi) = ($hi_guess, $lo_guess) }

    $root  = ($lo_guess + $hi_guess) / 2;
    $step  = $step_old = abs($hi_guess - $lo_guess);
    $value = &{$func}( $root );
    return $root if abs($value) < $epsilon;
    $deriv = deriv($func, $root);

    for ($i = 0; $i < ITERATIONS; $i++) {

        # Is Newton's method applicable?  If so, use it.
        #
        if ( ( $deriv * ($root - $hi) - $value ) *
             ( $deriv * ($root - $lo) - $value ) < 0  and
            abs($value * 2) <= abs($step_old * $deriv) ) {

            ($step_old, $step) = ($step, $value / $deriv);
            return $root if $step == 0 and abs($value) < $epsilon;
            $root -= $step;

        # Otherwise, bisect the current high and low guesses.
        #
        } else {
            ($step_old, $step) = ($step, ($hi - $lo) / 2);
            $root = $lo + $step;
            return $root if $lo == $root and abs($value) < $epsilon;
        }

        return $root if abs($step) < $epsilon and abs($value) < $epsilon;

        $value = &{$func}( $root );
        $deriv = deriv($func, $root);

        if ($value < 0) { $lo = $root } else { $hi = $root }
    }
    return;   # Maximum number of iterations reached.
}

sub cashflow {
    my $x = shift;
    return unless $x > 0;
    return ($x ** 3) - (100 * log($x));
}

print root(\&cashflow, 2, 100, .001);


# deriv($func, $x, $delta) approximates the derivative of $func (a
# code reference) at $x. If provided, $delta is used as dx; otherwise,
# it begins at $delta = 1e-31 and increments by an order of magnitude
# until a just noticeable difference is reached.
#
# If the function $func is discontinuous, all bets are off.
#
sub deriv {
    my ($func, $x, $delta) = @_;

    # Choose a delta if one wasn't provided.
    #
    unless ($delta) {
        $delta = 1e-31;
        while ($x == $x + $delta) { $delta *= 10 }
    }

    # Compute and return an approximation to the derivative.
    #
    return ( &{$func}( $x + $delta ) - &{$func}( $x ) ) / $delta;
}