December 12 2016 |
Multisection: When You Need to Bisect More Than Once |
Back
Next
|
# perl
use strict;
use warnings;
use 5.10.1;
use Carp;
use Scalar::Util qw( looks_like_number );
my ($f, $g) = (1,2);
my $nmax = 1000;
my $pre = 6;
my $polynomial_function = sub {
my $x = shift;
my $rv = $x**3 - $x - 2;
return $rv;
};
my ($rv, $n) = bisect_for_root( {
lower => $f,
upper => $g,
nmax => $nmax,
pre => $pre,
poly => $polynomial_function,
} );
my $format = '%' . $pre . 'f';
my $rounded = sprintf($format => $rv);
say "result: $rounded in $n steps";
##### SUBROUTINES #####
sub sgn {
my $x = shift;
croak "Must be a number" unless looks_like_number($x);
return -1 if $x < 0;
return 1 if $x > 0;
return 0;
}
sub bisect_for_root {
my $href = shift;
croak "Must supply hashref to bisect_for_root()"
unless ref($href) eq 'HASH';
for my $k ( qw| lower upper pre nmax | ) {
croak "Must supply element with key '$k'"
unless exists $href->{$k};
}
my $epsilon = 10 ** (0-$href->{pre});
my $n = 1;
while ($n < $href->{nmax}) {
my $h = ($href->{lower} + $href->{upper}) / 2;
if ( ($href->{poly}($h) == 0) or ((($href->{upper} - $href->{lower}) / 2) < $epsilon) ) {
return ($h, $n);
}
(sgn($href->{poly}($h)) == sgn($href->{poly}($href->{lower})))
? $href->{lower} = $h
: $href->{upper} = $h;
$n++;
}
croak "Unable to calculate result in $n attempts";
}