Perl's brain-dead floating point math 
Author Message
 Perl's brain-dead floating point math


:There is now a BigInt package in perl.  A FixedDecimal package could
:be fairly easily implemented as a subclass of this package.  Or
:someone could do a true FixedDecimal package from scratch.  Several
:years ago I wrote an Aribitrary Precision math package in C for use
:within a financial-based system I was developing, and so perhaps if I
:have some time (next decade? :) I might be able to use my experience
:in this area to write the FixedDecimal package in perl.

I created a BigFloat.pm -- I haven't done the BigRat.pm yet, but
it should be straightforward.  Would that make you happy?

I've been thinking that a good printing function is needed, and wondered
whether $BigFloat::OVERLOAD{'"'} might be better than BigFloat->fmt().
That way it would interpolate into strings in a nicer way.  Or I could
tie it with a you a TIE_SCLAR and a FETCH method.

--tom



Fri, 14 Mar 1997 14:25:14 GMT  
 Perl's brain-dead floating point math

:I've been thinking that a good printing function is needed, and wondered
:whether $BigFloat::OVERLOAD{'"'} might be better than BigFloat->fmt().
:That way it would interpolate into strings in a nicer way.  Or I could
:tie it with a you a TIE_SCLAR and a FETCH method.

Here's my test version of BigFloat.pm.  You may use $n->fsqrt().  There is no
pow() or mod() arithmagical method.  It has a useful (?) stringify() method,
so you should be able to just print out numbers.  Her'es my test.

    #!/usr/bin/perl
    use English; $ORS = "\n";
    use BigFloat;
    $num = new BigFloat 43545.0;
    print $num;
    for (1..10) { $num /= 10.0; print $num; }
    for (1..3) { $num *= $num; print $num; }
    print $num=$num->fsqrt();
    for (1..50) { $num /= 10.0; print $num; }

--tom

package BigFloat;

use Assert;
use BigInt;

%OVERLOAD = (
                                # Anonymous subroutines:
'+'     =>   sub {new BigFloat &fadd},
'-'     =>   sub {new BigFloat
                       $_[2]? fsub($_[1],${$_[0]}) : fsub(${$_[0]},$_[1])},
'<=>'     =>   sub {new BigFloat
                       $_[2]? fcmp($_[1],${$_[0]}) : fcmp(${$_[0]},$_[1])},
'cmp'   =>   sub {new BigFloat
                       $_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])},
'*'     =>   sub {new BigFloat &fmul},
'/'     =>   sub {new BigFloat
                       $_[2]? scalar fdiv($_[1],${$_[0]}) :
                         scalar fdiv(${$_[0]},$_[1])},
'neg'   =>   sub {new BigFloat &fneg},
'abs'   =>   sub {new BigFloat &fabs},

qw(
""    stringify
0+      numify)                 # Order of arguments unsignificant
);

sub new {
  my $foo = fnorm($_[1]);
  panic("Not a number initialized to BigFloat") if $foo eq "NaN";
  bless \$foo;

Quote:
}

sub numify { 0 + "${$_[0]}" } # Not needed, additional overhead
                                # comparing to direct compilation based on
                                # stringify
sub stringify {
    my $n = ${$_[0]};

    $n =~ s/^\+//;
    $n =~ s/E//;

    $n =~ s/([-+]\d+)$//;

    my $e = $1;
    my $ln = length($n);

    if ($e > 0) {
        $n .= "0" x $e . '.';
    } elsif (abs($e) < $ln) {
        substr($n, $ln + $e, 0) = '.';
    } else {
        $n = '.' . ("0" x (abs($e) - $ln)) . $n;
    }

    # 1 while $n =~ s/(.*\d)(\d\d\d)/$1,$2/;

    return $n;

Quote:
}

# Arbitrary length float math package
#
# by Mark Biggar
#
# number format
#   canonical strings have the form /[+-]\d+E[+-]\d+/
#   Input values can have inbedded whitespace
# Error returns
#   'NaN'           An input parameter was "Not a Number" or
#                       divide by zero or sqrt of negative number
# Division is computed to
#   max($div_scale,length(dividend)+length(divisor))
#   digits by default.
# Also used for default sqrt scale

$div_scale = 40;

# Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.

$rnd_mode = 'even';

sub fadd; sub fsub; sub fmul; sub fdiv;
sub fneg; sub fabs; sub fcmp;
sub fround; sub ffround;
sub fnorm; sub fsqrt;

#   bigfloat routines
#
#   fadd(NSTR, NSTR) return NSTR            addition
#   fsub(NSTR, NSTR) return NSTR            subtraction
#   fmul(NSTR, NSTR) return NSTR            multiplication
#   fdiv(NSTR, NSTR[,SCALE]) returns NSTR   division to SCALE places
#   fneg(NSTR) return NSTR                  negation
#   fabs(NSTR) return NSTR                  absolute value
#   fcmp(NSTR,NSTR) return CODE             compare undef,<0,=0,>0
#   fround(NSTR, SCALE) return NSTR         round to SCALE digits
#   ffround(NSTR, SCALE) return NSTR        round at SCALEth place
#   fnorm(NSTR) return (NSTR)               normalize
#   fsqrt(NSTR[, SCALE]) return NSTR        sqrt to SCALE places

# Convert a number to canonical string form.
#   Takes something that looks like a number and converts it to
#   the form /^[+-]\d+E[+-]\d+$/.
sub fnorm { #(string) return fnum_str

    s/\s+//g;                               # strip white space
    if (/^([+-]?)(\d*)(\.(\d*))?([Ee]([+-]?\d+))?$/ && "$2$4" ne '') {
        &norm(($1 ? "$1$2$4" : "+$2$4"),(($4 ne '') ? $6-length($4) : $6));
    } else {
        'NaN';
    }

Quote:
}

# normalize number -- for internal use
sub norm { #(mantissa, exponent) return fnum_str

    if ($_ eq 'NaN') {
        'NaN';
    } else {
        s/^([+-])0+/$1/;                        # strip leading zeros
        if (length($_) == 1) {
            '+0E+0';
        } else {
            $exp += length($1) if (s/(0+)$//);  # strip trailing zeros
            sprintf("%sE%+ld", $_, $exp);
        }
    }

Quote:
}

# negation
sub fneg { #(fnum_str) return fnum_str
    local($_) = fnorm($_[$[]);
    vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0E+0'; # flip sign
    s/^H/N/;
    $_;

Quote:
}

# absolute value
sub fabs { #(fnum_str) return fnum_str
    local($_) = fnorm($_[$[]);
    s/^-/+/;                                   # mash sign
    $_;

Quote:
}

# multiplication
sub fmul { #(fnum_str, fnum_str) return fnum_str
    local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1]));
    if ($x eq 'NaN' || $y eq 'NaN') {
        'NaN';
    } else {
        local($xm,$xe) = split('E',$x);
        local($ym,$ye) = split('E',$y);
        &norm(BigInt::bmul($xm,$ym),$xe+$ye);
    }

Quote:
}

# addition
sub fadd { #(fnum_str, fnum_str) return fnum_str
    local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1]));
    if ($x eq 'NaN' || $y eq 'NaN') {
        'NaN';
    } else {
        local($xm,$xe) = split('E',$x);
        local($ym,$ye) = split('E',$y);
        ($xm,$xe,$ym,$ye) = ($ym,$ye,$xm,$xe) if ($xe < $ye);
        &norm(BigInt::badd($ym,$xm.('0' x ($xe-$ye))),$ye);
    }

Quote:
}

# subtraction
sub fsub { #(fnum_str, fnum_str) return fnum_str
    fadd($_[$[],fneg($_[$[+1]));    

Quote:
}

# division
#   args are dividend, divisor, scale (optional)
#   result has at most max(scale, length(dividend), length(divisor)) digits
sub fdiv #(fnum_str, fnum_str[,scale]) return fnum_str
{
    local($x,$y,$scale) = (fnorm($_[$[]),fnorm($_[$[+1]),$_[$[+2]);
    if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0E+0') {
        'NaN';
    } else {
        local($xm,$xe) = split('E',$x);
        local($ym,$ye) = split('E',$y);
        $scale = $div_scale if (!$scale);
        $scale = length($xm)-1 if (length($xm)-1 > $scale);
        $scale = length($ym)-1 if (length($ym)-1 > $scale);
        $scale = $scale + length($ym) - length($xm);
        &norm(&round(BigInt::bdiv($xm.('0' x $scale),$ym),$ym),
            $xe-$ye-$scale);
    }

Quote:
}

# round int $q based on fraction $r/$base using $rnd_mode
sub round { #(int_str, int_str, int_str) return int_str

    if ($q eq 'NaN' || $r eq 'NaN') {
        'NaN';
    } elsif ($rnd_mode eq 'trunc') {
        $q;                         # just truncate
    } else {
        local($cmp) = BigInt::bcmp(BigInt::bmul($r,'+2'),$base);
        if ( $cmp < 0 ||
                 ($cmp == 0 &&
                  ( $rnd_mode eq 'zero'                             ||
                   ($rnd_mode eq '-inf' && (substr($q,$[,1) eq '+')) ||
                   ($rnd_mode eq '+inf' && (substr($q,$[,1) eq '-')) ||
                   ($rnd_mode eq 'even' && $q =~ /[24680]$/)        ||
                   ($rnd_mode eq 'odd'  && $q =~ /[13579]$/)        )) ) {
            $q;                     # round down
        } else {
            BigInt::badd($q, ((substr($q,$[,1) eq '-') ? '-1' : '+1'));
                                    # round up
        }
    }

Quote:
}

# round the mantissa of $x to $scale digits
sub fround { #(fnum_str, scale) return fnum_str
    local($x,$scale) = (fnorm($_[$[]),$_[$[+1]);
    if ($x eq 'NaN' || $scale <= 0) {
        $x;
    } else {
        local($xm,$xe) = split('E',$x);
        if (length($xm)-1 <= $scale) {
            $x;
        } else {
            &norm(&round(substr($xm,$[,$scale+1),
                         "+0".substr($xm,$[+$scale+1,1),"+10"),
                  $xe+length($xm)-$scale-1);
        }
    }

Quote:
}

# round $x at the 10 to the $scale digit place
sub ffround { #(fnum_str, scale) return fnum_str
    local($x,$scale) = (fnorm($_[$[]),$_[$[+1]);
    if ($x eq 'NaN') {
        'NaN';
    } else {
        local($xm,$xe) = split('E',$x);
        if ($xe >= $scale) {
            $x;
        } else {
            $xe = length($xm)+$xe-$scale;
            if ($xe < 1) {
                '+0E+0';
            } elsif ($xe == 1) {
                &norm(&round('+0',"+0".substr($xm,$[+1,1),"+10"), $scale);
            } else {
                &norm(&round(substr($xm,$[,$xe),
                      "+0".substr($xm,$[+$xe,1),"+10"), $scale);
            }
        }
    }

Quote:
}

# compare 2 values returns one of undef, <0, =0, >0
#   returns undef if either or both input value are not numbers
sub fcmp #(fnum_str, fnum_str) return cond_code
{
    local($x, $y) = (fnorm($_[$[]),fnorm($_[$[+1]));
    if ($x eq "NaN" || $y eq "NaN") {
        undef;
    } else {
        ord($y) <=> ord($x)
        ||
        (  local($xm,$xe,$ym,$ye) = split('E', $x."E$y"),
             (($xe <=> $ye) * (substr($x,$[,1).'1')
             || BigInt::cmp($xm,$ym))
        );
    }

Quote:
}

# square root by Newtons method.
sub fsqrt { #(fnum_str[, scale]) return fnum_str
    local($x, $scale) = (fnorm($_[$[]), $_[$[+1]);
    if ($x eq 'NaN' || $x =~ /^-/) {
        'NaN';
    } elsif ($x eq '+0E+0') {
        '+0E+0';
    } else {
        local($xm, $xe) = split('E',$x);
        $scale = $div_scale if (!$scale);
        $scale = length($xm)-1 if ($scale < length($xm)-1);
        local($gs, $guess) = (1, sprintf("1E%+d", (length($xm)+$xe-1)/2));
        while ($gs < 2*$scale) {
            $guess = fmul(fadd($guess,fdiv($x,$guess,$gs*2)),".5");
            $gs *= 2;
        }
        new BigFloat fround($guess, $scale);
    }

Quote:
}

1;


Fri, 14 Mar 1997 15:56:34 GMT  
 Perl's brain-dead floating point math
Here, you'll need Assert.pm as well.

--tom

# Assert.pm

#
# Usage:


#
# That is, if the first expression evals false, we blow up.  The
# rest of the args, if any, are nice to know because they will
# be printed out by &panic, which is just the stack-backtrace
# routine shamelessly borrowed from the perl de{*filter*}.

package Assert;


sub assert {
    my $code = $_[$[];
    my $ok;
    if (ref($code) eq 'CODE') {
        $ok = &$code;
    } else {
        $ok = eval $_[$[];
    }

Quote:
}

sub panic {
    select(STDERR);

    # stack traceback gratefully borrowed from perl5 de{*filter*}
    local($i,$_);

    for ($i = 1; ($p,$f,$l,$s,$h,$w) = caller($i); $i++) {


            s/'/\\'/g;
            s/([^\0]*)/'$1'/ unless /^-?[\d.]+$/;
            s/([\200-\377])/sprintf("M-%c",ord($1)&0177)/eg;
            s/([\0-\37\177])/sprintf("^%c",ord($1)^64)/eg;
        }



        last if $signal;
    }
    for ($i=0; $i <= $#sub; $i++) {
        last if $signal;
        print STDERR $sub[$i];
    }
    die "bailing out from panic";

Quote:
}

1;


Fri, 14 Mar 1997 15:59:00 GMT  
 Perl's brain-dead floating point math

Quote:

>Why is Perl so stupid when it comes to floating point math?  

>Example:

>    $num = 43545.0;

>    $num = $num / 10.0;
>    print "$num\n";

>    $num = $num / 10.0;
>    print "$num\n";

>    $num = $num / 10.0;
>    print "$num\n";

>Results:

>4354.5
>435.44999999999998863
>43.545000000000001705

>This example program has worked incorrectly on every machine I have
>tried it on, therefore I have concluded that it isn't a problem with
>the installation on my machine.

>Any thoughts?

Perl 4 defaults to using significant 20 digits when printing floats.
Most floating point implementations only have 15 significant digits.
Perl 5 fixes this:

$ ./perl5
$num = 43545.0;
$num = $num / 10.0;
print "$num\n";
$num = $num / 10.0;
print "$num\n";
$num = $num / 10.0;
print "$num\n";
^D
4354.5
435.45
43.545
$

Quote:
>Rob

Regards,
Tim Bunce.

(Actually the perl5-porters snuck this in when Larry took a break for a
week or so a few months back. Larry blessed it on his return. At the
same time the old sprintf() method was changed to the much faster
gconvert() on platforms which support it.)



Fri, 14 Mar 1997 20:33:43 GMT  
 Perl's brain-dead floating point math

Quote:


> :There is now a BigInt package in perl.  A FixedDecimal package could
> :be fairly easily implemented as a subclass of this package.  ...
> I created a BigFloat.pm -- I haven't done the BigRat.pm yet, but
> it should be straightforward.  Would that make you happy?

Is BigFloat floating point math or fixed decimal math?  If it's just
floating point with more precision than standard floats and doubles,
it still may not solve the general problem of precision that fixed
decimals solve (see my previous post in this thread about fixed
decimals and financial applications).

If BigFloat is really fixed decimal, I recommend a name change to remove
the word "Float" from the package name.

A BigRat package would be really nice!  With this package, fixed
decimal math would be trivial ... just use denominators of 10, 100,
1000, etc.  and change the way the results are printed in "human
readable" form (i.e., display "123456 / 100" as "1234.56").

Quote:
> I've been thinking that a good printing function is needed, and wondered
> whether $BigFloat::OVERLOAD{'"'} might be better than BigFloat->fmt().
> That way it would interpolate into strings in a nicer way.  Or I could
> tie it with a you a TIE_SCLAR and a FETCH method.

Can '"' really be overloaded?  If so, that's really cool!

--
Lloyd Zusman            01234567 <-- The world famous Indent-o-Meter.

   To get my PGP public key automatically mailed to you, please
   send me email with the following string as the subject:
                    mail-request public-key



Sat, 15 Mar 1997 00:24:08 GMT  
 Perl's brain-dead floating point math

Quote:

> Here's my test version of BigFloat.pm.  You may use $n->fsqrt().  There is no
> pow() or mod() arithmagical method.  It has a useful (?) stringify() method,
> so you should be able to just print out numbers.  Her'es my test.

I will check this out and report back here ... upon cursory perusal of
the code, it looks like this is more of an arbitrary-precision float
package than a fixed decimal package.  It probably could be subclassed
into some sort of fixed decimal thing, but I think the BigRat package
you mentioned in your previous post would be better suited for this.

This looks really good, though!

--
Lloyd Zusman            01234567 <-- The world famous Indent-o-Meter.

   To get my PGP public key automatically mailed to you, please
   send me email with the following string as the subject:
                    mail-request public-key



Sat, 15 Mar 1997 00:30:36 GMT  
 Perl's brain-dead floating point math

Quote:


>> Here's my test version of BigFloat.pm.  You may use $n->fsqrt().  There is no
>> pow() or mod() arithmagical method.  It has a useful (?) stringify() method,
>> so you should be able to just print out numbers.  Her'es my test.
>I will check this out and report back here ... upon cursory perusal of
>the code, it looks like this is more of an arbitrary-precision float
>package than a fixed decimal package.  It probably could be subclassed
>into some sort of fixed decimal thing, but I think the BigRat package
>you mentioned in your previous post would be better suited for this.

BigFloat.pm and the original bigfloat.pl packages are arbitrary precision
decimal floating point math built on top of the BigInt package.
Values are all of the form /[+-]\d+E[+-]\d+/ (i.e. the assumed decimal
point is always at the right end of the mantissa.  It would be very easy
to subclass out a real fixed point type where blessed variables internally
contained a scale (and maybe a precision) and did real fixed point math.
As you said you could also base it on bigrat, but I thing that bigrat's
least common demoninator format would get in the way, so I'd stick with
BigFloat.  The main thing missing from Bigfloat is pretty printing routines
to do %f, %e and %g formating.

Hey Larry, now that we have operator overloading could we get hooks
into printf and friends so a class can define % formaters for the class?
This would be easy, if printf or sprintf see a blessed reference as an argument
it would call the printf method for the object passing the next % formatter
as an argument and then put the return string value into the output.

--
Perl's Maternal Uncle
Mark Biggar



Sat, 15 Mar 1997 05:26:43 GMT  
 Perl's brain-dead floating point math
Floating point is a crutch.  Real programmers use integer arithmetic and
scale as necessary. (I'm a Perl fanatic in spite of its lack of integers.)
    ;-)
               ___    Dick Chiles
              /___\
      ____/---(___)     AT&T Global Information Solutions
    /  _ __   \___/ \    17095 Via del Campo
   |_|__/  \      .  \    San Diego, CA 92127-1711
    ---/\__/\\  _/   /

   _/// `/    .'  |     Voice: (619) 485-3858
   \// //_-\_ |   |    Fax:   (619) 485-3788
   // //__   \_\  |
   / //    `.    /
  `_/        ` -'   My employer has nothing to do with this.


Sun, 16 Mar 1997 22:37:42 GMT  
 Perl's brain-dead floating point math

Quote:
>Floating point is a crutch.  Real programmers use integer arithmetic and
>scale as necessary. (I'm a Perl fanatic in spite of its lack of integers.)

Get Perl5, inside a package you can say

use integer;

and all math will then be done in 32 ints instead of float.

--
Mark Biggar



Mon, 17 Mar 1997 01:18:09 GMT  
 Perl's brain-dead floating point math

Quote:

>Hey Larry, now that we have operator overloading could we get hooks
>into printf and friends so a class can define % formaters for the class?
>This would be easy, if printf or sprintf see a blessed reference as an argument
>it would call the printf method for the object passing the next % formatter
>as an argument and then put the return string value into the output.

An *excellent* idea!

Quote:
>Perl's Maternal Uncle
>Mark Biggar


Regards,
Tim Bunce.


Mon, 17 Mar 1997 16:29:58 GMT  
 
 [ 36 post ]  Go to page: [1] [2] [3]

 Relevant Pages 
 

 
Powered by phpBB® Forum Software