I just (more fully) realized how easy it is to do calculations with huge numbers in Perl so long as you don't have some unreasonable desire for a huge number of digits of precision (I settled for 10, which I think is usually too many anyway).

I threw this module prototype together in just a few minutes. Support for negative numbers and even zero is left as an exercise for the next round, however. It works surprisingly well (and fast). It is great for computing things like how many combinations / permutations some unreasonably large set of things has. 

package Math::BigPositiveOkayPrecision;

use overload
    '+' =>  \&add,
    #'-' =>  \&sub,
    '*' =>  \&mul,
    '/' =>  \&div,
    '**' => \&pow,
    '^' =>  \&prd,  # $x^$y = factorial($y)/factorial($x)
    '<<' => \&shl,
    '>>' => \&shr,
    '<=>'=> \&cmp,
    '!' =>  \&fct,  # !$x = factorial($x)
    '0+' => \&num,
    '""' => \&str,
    #'neg'=> \&neg,
;

sub new
{
    my $us= shift @_;
    if(  ! @_  ) {
        die "new() needs an object or a value"
            if  ! ref $us;
        return bless \( 0+$$us );
    }
    my $val= shift @_;
    die "Non-positive values not supported yet ($val)"
        if  $val <= 0;
    return bless \log( $val );
}
*Okay::new= \&new;  # Until we get a better module name

sub mul
{
    my( $x, $y )= @_;
    $y= $x->new( $y )
        if  ! ref $y;
    return bless \( $$x + $$y );
}

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

sub pow
{
    my( $x, $y, $rev )= @_;
    $y= $x->new( $y )
        if  ! ref $y;
    ( $x, $y )= ( $y, $x )
        if  $rev;
    return bless \( $$x * exp($$y) );
}

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

sub shl
{
    my( $x, $y, $rev )= @_;
    $y= $x->new( $y )
        if  ! ref $y;
    ( $x, $y )= ( $y, $x )
        if  $rev;
    return bless \( $$x + exp($$y)*log(2) );
}

sub shr
{
    my( $x, $y, $rev )= @_;
    $y= $x->new( $y )
        if  ! ref $y;
    ( $x, $y )= ( $y, $x )
        if  $rev;
    return bless \( $$x - exp($$y)*log(2) );
}

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

sub prd
{
    my( $x, $y, $rev )= @_;
    $y= $x->new( $y )
        if  ! ref $y;
    ( $x, $y )= ( $y, $x )
        if  $rev;
    return bless \0
        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 exp( $$x );
}

sub str
{
    my( $x )= @_;
    my $exp= int( $$x / log(10) );
    my $mant= exp( $$x - log(10)*$exp );
    $mant= sprintf "%.10f", $mant;
    $mant .=  $exp  ?  "e" . $exp  :  "";
    return $mant
        if  9 < $exp;
    return 0 + $mant;
}

__PACKAGE__;

[download]

See Re^6: One Zero variants_without_repetition (logue) for a sample use and the most recent part of my inspiration for this.

- tye        

Replies are listed 'Best First'.
Re: Math::BigPositiveOkayPrecision prototype
by dogz007 (Scribe) on Aug 08, 2007 at 19:29 UTC
    Ahh, I understand now. Quite simple, yet very effective. I'm a mechanical engineer, so I find the your problem interesting. One suggestion that I would make (although it may not help that much) would be to make the individual operation methods simpler by extracting their common code to another method. For instance, the beginning portions of mul, div, pow, add, shl, shr, cmp, prd are all practically identical. You could write:

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

    [download]

    This would allow you to simplify the operation methods in a manner similar to that below.

    sub mul {
    my $mulsub = sub { bless \( $$_[0] + $$_[1] ) };
    return &compute($mulsub,@_);
}

    [download]

    This may have adverse effects on performance since it increases the number of function calls, but I'm not experienced enough to know how much. However, it might make your code easier to understand and maintain.

[reply]
[d/l]
[select]

      Thanks. I've added: 

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

      [download]

      So that boiler-plate reduces to: my( $x, $y )= __args( @_ ); in 6 places.

      - tye        

[reply]
[d/l]
[select]
Re: Math::BigPositiveOkayPrecision prototype
by dogz007 (Scribe) on Aug 08, 2007 at 17:04 UTC
    I'm so confused. What kind of math are you working with? I'm guessing that you're storing large numbers in exponential format (meaning converting the number to some power of e), and doing math on the powers instead of the actual base-10 numbers. Correct me if I'm wrong. I just thought your module was interesting and wanted to understand it better.
[reply]

      Yes, that's a good description of what it does.

      Quoting the documentation for my up-coming Math::BigApprox (improved name and now with support for negative numbers and zero):

      Math::BigApprox stores numbers as the logarithm of their (absolute) value (along with separately tracking their sign) and so can keep track of numbers having absolute value less than about 1e10000000000 with several digits of accuracy (and quite a bit of accuracy for not quite so ridiculously large numbers).

      Therefore, Math::BigApprox numbers never require a large amount of memory and calculations with them are rather fast (about as fast as you'll get with overloaded operations in Perl).

      It can even do calculations with numbers as large as 10**(1e300) (which is too large for this space even when written as "1e..."). But for such obscenely huge numbers it only knows their magnitude, having zero digits of precision in the mantissa (or "significand") but having reasonable precision in the base-10 exponent (the number after the "e"). It would actually display such a number as something like "1e1.0217e300" since the usual XeN scientific notation is too cumbersome for such obscenely large numbers.

      The simplest calculation is multiplication because we store $x= log($X) and $y= log($Y) so we want to store $z= log( $X * $Y ) and log($X*$Y) == log($X)+log($Y) so we just store $z= $x + $y (to multiply). As Noah said, even adders can multiply on a log table.

      The most interesting calculation is addition. 

          Store $x= log($X)
    Store $y= log($Y)
    Want  $z= log( $X + $Y )
    But we don't know $X and $Y, just $x and $y so
    Want  $z= log( exp($x) + exp($y) )

      [download]

      But exp($x) and/or exp($y) are likely too big to compute so we do the following derivation to get a safer formula: 

          $z= log( exp($x) + exp($y) )
      = log( exp($x)*( 1 + exp($y)/exp($x) ) )
      = log( exp($x) * ( 1 + exp($y-$x) ) )
      = log(exp($x)) + log( 1 + exp($y-$x) )
      = $x + log( 1 + exp($y-$x) )

      [download]

      Now, the only problem is ensuring that exp($y-$x) is not a problem to compute. But if we ensure that $x <= $y (by swapping $x and $y if necessary), then the worst case is something like exp(-1e300), which just gets approximated as 0. And that case means that $y is obscenely larger than $x so $y+$x just gives us back $y, as it should, since we don't have infinite precision.

      - tye        

[reply]
[d/l]
[select]
        Math::Sliderule?
[reply]
Re: Math::BigPositiveOkayPrecision prototype (Math::BigApprox)
by tye (Sage) on Aug 10, 2007 at 05:14 UTC
[reply]
A reply falls below the community's threshold of quality. You may see it by logging in.

Back to Cool Uses for Perl