File Coverage

blib/lib/Math/BigApprox.pm

Criterion	Covered	Total	%
statement	130	133	97.7
branch	83	88	94.3
condition	14	15	93.3
subroutine	30	31	96.7
pod	6	6	100.0
total	263	273	96.3

package Math::BigApprox;
use strict;

# use POSIX qw( floor );    # Maybe, maybe not.

use vars qw( $VERSION @EXPORT_OK );
BEGIN {
    $VERSION= 0.001_005;

    require Exporter;
    *import= \&Exporter::import;
    @EXPORT_OK= qw( c Fact Prod $SigDigs );

    my $useFloor= 0;
    if(  eval { require POSIX; 1 }  ) {
        # On "long double" Perls, POSIX::floor converts to "double"
        $useFloor=  "4e+102" eq POSIX::floor(4e102);
    }
    if(  $useFloor  ) {
        *_floor= \&POSIX::floor;
    } else {
        *_floor= sub {
            return  $_[0] < 0  ?  -int(-$_[0])  :  int($_[0]);
        };
    }
}

use overload(
    '+' =>  \&_add,
    '-' =>  \&_sub,
    neg =>  \&_neg,
    '*' =>  \&_mul,
    '/' =>  \&_div,
    '**' => \&_pow,
    '^' =>  \&_prd,  # $x^$y = factorial($y)/factorial($x)
    '<<' => \&_shl,
    '>>' => \&_shr,
    '<=>'=> \&_cmp,
    '!' =>  \&_fct,  # !$x = factorial($x)
    '0+' => \&_num,
    '""' => \&_str,
    '=' =>  \&_clone,
    abs =>  \&_abs,
    log =>  \&_log,
    sqrt => \&_half,
    bool => \&_chastise,
    fallback => 1,  # autogenerate if possible
);

# Don't use these as a methods
sub c
{
    return __PACKAGE__->new( @_ );
}

sub Prod
{
    my( $x, $y )= @_;
    my $p= c( 1 );
    while(  $x <= $y  ) {
        $p *= $x++;
    }
    return $p;
}

sub Fact
{
    return Prod( 1, @_ );
}

# Public methods:
sub new
{
    my $us= shift @_;
    if(  ! @_  ) {
        die "new() needs an object or a value"
            if  ! ref $us;
        return bless [ @$us ];
    }
    my $val= shift @_;
    return bless [ 0, 0 ]
        if  0 == $val;
    return bless [ log( -$val ), -1 ]
        if  $val < 0;
    return bless [ log( $val ), 1 ]
}

sub Get
{
    my( $x )= @_;
    return wantarray ? (0,0) : 0
        if  ! $x->[1];
    return wantarray ? @$x : $x->[0];
}

sub Sign
{
    my $x= shift @_;
    my $old= $x->[1];
    if(  @_  ) {
        my $new= shift @_;
        $x->[1]= 0<$new ? 1 : 0==$new ? 0 : -1;
    }
    return $old;
}

# Private methods (please don't call these):

# This just ignores extra arguments so overload.pm can use it
sub _clone
{
    return $_[0]->new();
}

sub _chastise
{
    die "You can't use a Math::BigApprox as a Boolean.\n";
}

sub _mul
{
    my( $x, $y )= @_;
    $y= $x->new( $y )
        if  ! ref $y;
    return bless [  $x->[0] + $y->[0],  $x->[1] * $y->[1]  ];
}

sub __args
{
    my( $x, $y, $rev )= @_;
    $y= $x->new( $y )
        if  ! ref $y;
    ( $x, $y )= ( $y, $x )
        if  $rev;
    return( $x, $y );
}

sub _div
{
    my( $x, $y )= __args( @_ );
    die "Division by 0"
        if  ! $y->[1];
    return bless [  $x->[0] - $y->[0],  $x->[1] * $y->[1]  ];
}

sub _pow
{
    my( $x, $y )= __args( @_ );
    return bless [ 0, 1 ]
        if  ! $y->[1];
    return bless [ 0, 0 ]
        if  ! $x->[1];
    # We rather "punt" on the sign.
    my $sign= 1;
    if(     $x->[1] < 0
        &&  0 < $y->[1]
        &&  $y->[0] < log(1e9)
        &&  1 == int( 0.5 + exp($y->[0]) ) % 2
    ) {
        $sign= -1;
    }
    return bless [  $x->[0] * $y->_num(),  $sign  ];
}

sub _add
{
    my( $x, $y )= @_;
    $y= $x->new( $y )
        if  ! ref $y;
    return $y->new()
        if  ! $x->[1];
    return $x->new()
        if  ! $y->[1];
    ( $x, $y )= ( $y, $x )
        if  $y->[0] < $x->[0];
    if(  $x->[1] == $y->[1]  ) {
        return bless [
            $y->[0] + log( 1 + exp( $x->[0] - $y->[0] ) ),
            $x->[1],
        ];
    }
    return bless [ 0, 0 ]
        if  $x->[0] == $y->[0];
    return bless [
        $y->[0] + log( 1 - exp( $x->[0] - $y->[0] ) ),
        $y->[1],
    ];
}

sub _sub
{
    my( $x, $y )= __args( @_ );
    return $x->_add( $y->_neg() );
}

sub _neg
{
    my( $x )= @_;
    return bless [ $x->[0], -$x->[1] ];
}

sub _shl
{
    my( $x, $y )= __args( @_ );
    return bless [ $x->[0] + $y->_num()*log(2), $x->[1] ];
}

sub _shr
{
    my( $x, $y )= __args( @_ );
    return bless [ $x->[0] - $y->_num()*log(2), $x->[1] ];
}

sub _log
{
    my( $x )= @_;
    die "Can't take log(" . $x . ")\n"
        if  1 != $x->[1];
    return $x->[0];
}

sub _abs
{
    my( $x )= @_;
    return bless [ $x->[0], abs( $x->[1] ) ];
}

sub _half
{
    my( $x )= @_;
    die "Can't take sqrt(" . $x . ")\n"
        if  $x->[1] < 0;
    return bless [ $x->[0]/2, $x->[1] ];
}

sub _cmp
{
    my( $x, $y, $rev )= @_;
    $y= $x->new( $y )
        if  ! ref $y;
    return 0
        if  $x eq $y;
    ( $x, $y )= ( $y, $x )
        if  $rev;
    return $x->[1] <=> $y->[1] || $x->[1]*$x->[0] <=> $y->[1]*$y->[0];
}

sub _prd
{
    my( $x, $y )= __args( @_ );
    return bless [ 0, 1 ]
        if  $y < $x;
    my $p= $x->new();
    my $m= $x + 1;
    while(  $m <= $y  ) {
        $p= $p * $m;
        $m= $m + 1;
    }
    return $p;
}

sub _fct
{
    my( $x )= @_;
    return 1^$x;
}

sub _num
{
    my( $x )= @_;
    return $x->[1] * exp( $x->[0] );
}

use vars qw( $SigDigs $FloorMag $LenMag );
# Figure out how many significant digits are in an NV on this platform:
BEGIN {

    # Cheap trick to get at how many sig digs Perl figured out it has:
    $SigDigs= length( 10 / 7 );

    # Don't call floor() on numbers larger than this, since they can't
    # have a factional part [and floor(4e102) can substract 8e80!]
    $FloorMag= "1e" . $SigDigs;

    # Minus 1 for decmical point, minus another two just
    # because our calculations lose some precision:
    $SigDigs -= 3;

    # Long doubles leave $SigDigs set too high (probably
    # not all C RTL calls are fully long-double-using):
    $SigDigs -= 3
        if  14 < $SigDigs;

    # Don't us length() to measure magnitude on numbers larger than this:
    $LenMag= "1e" . $SigDigs;
}

sub _str
{
    my( $x )= @_;
    return 0
        if  ! $x->[1];
    my $exp= $x->[0] / log(10);
    if(  $exp  &&  2*$exp == $exp  ) {
        return 0
            if  $exp < 0;
        return $x->[1] * $exp;
    }
    $exp= sprintf "%.*g", $SigDigs-1, $exp;
    $exp= _floor( $exp )
        if  abs($exp) < $FloorMag;
    my $mant= exp( $x->[0] - log(10)*$exp );
    $mant=  1
        if  2*$mant == $mant;
    $exp= "+$exp"
        if  $exp !~ /^-/;
    my $digs=  $LenMag <= abs($exp)  ?  1  :  $SigDigs - length($exp);
    $digs= 1
        if  $digs < 1;
    $mant= sprintf "%s%.*f", $x->[1] < 0 ? '-' : '', $digs-1, $mant;
    $mant =~ s/[.]?0+$//;
    $mant .=  0==$exp  ?  ""  :  "e" . $exp;
    return $mant
        if  $digs < abs($exp);
    return 0+$mant;
}

__PACKAGE__;
__END__