File Coverage

blib/lib/Math/Rational/Approx.pm

Criterion	Covered	Total	%
statement	72	72	100.0
branch	18	22	81.8
condition	9	12	75.0
subroutine	16	16	100.0
pod	3	3	100.0
total	118	125	94.4

# --8<--8<--8<--8<--
#
# Copyright (C) 2012 Smithsonian Astrophysical Observatory
#
# This file is part of Math::Rational::Approx
#
# Math::Rational::Approx is free software: you can redistribute it
# and/or modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see .
#
# -->8-->8-->8-->8--

package Math::Rational::Approx;

use strict;
use warnings;
use Carp;


our $VERSION = '0.02';

use base 'Exporter';
our %EXPORT_TAGS = ( all => [ qw( maxD contfrac contfrac_nd ) ],
                   );
our @EXPORT_OK = map { @$_} values %EXPORT_TAGS;

use POSIX 'floor';
use Math::BigFloat;

use MooX::Types::MooseLike::Numeric ':all';
use Params::Validate qw[ validate_pos ARRAYREF ];

sub _check_bounds {

        my ( $x, $bounds ) = @_;

        croak( "incorrect number of elements in bounds\n" )
          unless 4 == @$bounds;

        my ( $a, $b, $c, $d ) = @$bounds;
        # ensure that a/b < c/d

        croak( "a/b is not less than c/d\n" )
          unless $a / $b < $c / $d;


        croak( "a/b and c/d do not bound x\n" )
          unless $a / $b < $x && $x < $c / $d;

}

sub maxD {

        my ($x, $maxD, $bounds )
          = validate_pos( @_,
                          { callbacks => { 'positive number' => sub { is_PositiveNum($_[0]) } },
                          },
                          { callbacks => { 'positive integer' => sub { is_PositiveInt($_[0]) } },
                          },
                         {  type => ARRAYREF,
                            default => [],
                         },
                        ) ;

        if ( defined $bounds && @$bounds ) {

                _check_bounds( $x, $bounds );
        }

        else {

                my $base = floor( $x );

                @{$bounds} = ( $base, 1, $base+1, 1 );

        }

        my ( $a, $b, $c, $d ) = @{$bounds};

        my ( $nom, $denom );

        while ( $b <= $maxD && $d <= $maxD ) {

                my $mediant = ( $a + $c ) / ( $b + $d );

                if ( $x == $mediant ) {

                        ( $nom, $denom ) =
                          $b + $d <= $maxD ? ( $a + $c, $b + $d )
                        : $d > $b          ? ( $c, $d )
                        :                    ( $a, $b );

                        last;
                }

                elsif ( $x > $mediant ) {

                        ($a, $b) = ( $a + $c, $b + $d );

                }

                else {

                        ($c, $d) = ( $a + $c, $b + $d );

                }

        }

        if ( ! defined $nom && ! defined $denom ) {
                ( $nom, $denom ) = $b > $maxD ? ( $c, $d )
                                              : ( $a, $b );
        }

        @{$bounds} = ( $a, $b, $c, $d );

        return ( $nom, $denom, $bounds );
}

sub contfrac {

        my ( $x, $max, $p ) =
          validate_pos( @_,
                      { callbacks => {
                                      'positive number' => sub { is_PositiveNum($_[0]) },
                                     },
                      },
                      { callbacks => {
                                      'positive integer' => sub { is_PositiveInt($_[0]) },
                                     },
                      },
                      { type => ARRAYREF, default => [] },
                      );

        $x = Math::BigFloat->new( $x );

        my $one = Math::BigFloat->bone;

        for my $n ( 0 .. $max-1 ) {

                my $int = $x->copy->bfloor;
                push @$p, $int;

                # if $x is actually rational, round off error in ( $x - $int )
                # can drive the iteration beyond its true end point, causing
                # bogus results.  Let's hope that 10 digits is enough.
                last if ( $x - $int)->bfround(-10)->is_zero;
                $x = $one->copy->bdiv( $x - $int );
        }

        return $p, $x;
}

sub contfrac_nd {

	# ignore extra parameter to ease use of contfrac_nd( contfrac ( ... ) )
    my ( $terms ) = validate_pos( @_, { type => ARRAYREF }, 0 );
    my @p = reverse @$terms;

    my $n = Math::BigInt->bone;
    my $d = Math::BigInt->new( (shift @p) );

    for my $p ( @p ) {

            $n += $d * $p;

            ( $n, $d ) = ( $d, $n );

    }

    ( $n, $d ) = ( $d, $n );
    return ( $n, $d );
}


1;



__END__