June 19 2019 Multisection: When You Need to Bisect More Than Once Back

Polynomial Root Problem Solved in Perl

        # 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";
        }


Home Last TOC Copyright © 2019 James E Keenan Back Next