1N/Apackage bigfloat;
1N/Arequire "bigint.pl";
1N/A#
1N/A# This library is no longer being maintained, and is included for backward
1N/A# compatibility with Perl 4 programs which may require it.
1N/A#
1N/A# In particular, this should not be used as an example of modern Perl
1N/A# programming techniques.
1N/A#
1N/A# Suggested alternative: Math::BigFloat
1N/A#
1N/A# Arbitrary length float math package
1N/A#
1N/A# by Mark Biggar
1N/A#
1N/A# number format
1N/A# canonical strings have the form /[+-]\d+E[+-]\d+/
1N/A# Input values can have embedded whitespace
1N/A# Error returns
1N/A# 'NaN' An input parameter was "Not a Number" or
1N/A# divide by zero or sqrt of negative number
1N/A# Division is computed to
1N/A# max($div_scale,length(dividend)+length(divisor))
1N/A# digits by default.
1N/A# Also used for default sqrt scale
1N/A
1N/A$div_scale = 40;
1N/A
1N/A# Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
1N/A
1N/A$rnd_mode = 'even';
1N/A
1N/A# bigfloat routines
1N/A#
1N/A# fadd(NSTR, NSTR) return NSTR addition
1N/A# fsub(NSTR, NSTR) return NSTR subtraction
1N/A# fmul(NSTR, NSTR) return NSTR multiplication
1N/A# fdiv(NSTR, NSTR[,SCALE]) returns NSTR division to SCALE places
1N/A# fneg(NSTR) return NSTR negation
1N/A# fabs(NSTR) return NSTR absolute value
1N/A# fcmp(NSTR,NSTR) return CODE compare undef,<0,=0,>0
1N/A# fround(NSTR, SCALE) return NSTR round to SCALE digits
1N/A# ffround(NSTR, SCALE) return NSTR round at SCALEth place
1N/A# fnorm(NSTR) return (NSTR) normalize
1N/A# fsqrt(NSTR[, SCALE]) return NSTR sqrt to SCALE places
1N/A
1N/A# Convert a number to canonical string form.
1N/A# Takes something that looks like a number and converts it to
1N/A# the form /^[+-]\d+E[+-]\d+$/.
1N/Asub main'fnorm { #(string) return fnum_str
1N/A local($_) = @_;
1N/A s/\s+//g; # strip white space
1N/A if (/^([+-]?)(\d*)(\.(\d*))?([Ee]([+-]?\d+))?$/
1N/A && ($2 ne '' || defined($4))) {
1N/A my $x = defined($4) ? $4 : '';
1N/A &norm(($1 ? "$1$2$x" : "+$2$x"), (($x ne '') ? $6-length($x) : $6));
1N/A } else {
1N/A 'NaN';
1N/A }
1N/A}
1N/A
1N/A# normalize number -- for internal use
1N/Asub norm { #(mantissa, exponent) return fnum_str
1N/A local($_, $exp) = @_;
1N/A if ($_ eq 'NaN') {
1N/A 'NaN';
1N/A } else {
1N/A s/^([+-])0+/$1/; # strip leading zeros
1N/A if (length($_) == 1) {
1N/A '+0E+0';
1N/A } else {
1N/A $exp += length($1) if (s/(0+)$//); # strip trailing zeros
1N/A sprintf("%sE%+ld", $_, $exp);
1N/A }
1N/A }
1N/A}
1N/A
1N/A# negation
1N/Asub main'fneg { #(fnum_str) return fnum_str
1N/A local($_) = &'fnorm($_[$[]);
1N/A vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0E+0'; # flip sign
1N/A if ( ord("\t") == 9 ) { # ascii
1N/A s/^H/N/;
1N/A }
1N/A else { # ebcdic character set
1N/A s/\373/N/;
1N/A }
1N/A $_;
1N/A}
1N/A
1N/A# absolute value
1N/Asub main'fabs { #(fnum_str) return fnum_str
1N/A local($_) = &'fnorm($_[$[]);
1N/A s/^-/+/; # mash sign
1N/A $_;
1N/A}
1N/A
1N/A# multiplication
1N/Asub main'fmul { #(fnum_str, fnum_str) return fnum_str
1N/A local($x,$y) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]));
1N/A if ($x eq 'NaN' || $y eq 'NaN') {
1N/A 'NaN';
1N/A } else {
1N/A local($xm,$xe) = split('E',$x);
1N/A local($ym,$ye) = split('E',$y);
1N/A &norm(&'bmul($xm,$ym),$xe+$ye);
1N/A }
1N/A}
1N/A
1N/A# addition
1N/Asub main'fadd { #(fnum_str, fnum_str) return fnum_str
1N/A local($x,$y) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]));
1N/A if ($x eq 'NaN' || $y eq 'NaN') {
1N/A 'NaN';
1N/A } else {
1N/A local($xm,$xe) = split('E',$x);
1N/A local($ym,$ye) = split('E',$y);
1N/A ($xm,$xe,$ym,$ye) = ($ym,$ye,$xm,$xe) if ($xe < $ye);
1N/A &norm(&'badd($ym,$xm.('0' x ($xe-$ye))),$ye);
1N/A }
1N/A}
1N/A
1N/A# subtraction
1N/Asub main'fsub { #(fnum_str, fnum_str) return fnum_str
1N/A &'fadd($_[$[],&'fneg($_[$[+1]));
1N/A}
1N/A
1N/A# division
1N/A# args are dividend, divisor, scale (optional)
1N/A# result has at most max(scale, length(dividend), length(divisor)) digits
1N/Asub main'fdiv #(fnum_str, fnum_str[,scale]) return fnum_str
1N/A{
1N/A local($x,$y,$scale) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]),$_[$[+2]);
1N/A if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0E+0') {
1N/A 'NaN';
1N/A } else {
1N/A local($xm,$xe) = split('E',$x);
1N/A local($ym,$ye) = split('E',$y);
1N/A $scale = $div_scale if (!$scale);
1N/A $scale = length($xm)-1 if (length($xm)-1 > $scale);
1N/A $scale = length($ym)-1 if (length($ym)-1 > $scale);
1N/A $scale = $scale + length($ym) - length($xm);
1N/A &norm(&round(&'bdiv($xm.('0' x $scale),$ym),&'babs($ym)),
1N/A $xe-$ye-$scale);
1N/A }
1N/A}
1N/A
1N/A# round int $q based on fraction $r/$base using $rnd_mode
1N/Asub round { #(int_str, int_str, int_str) return int_str
1N/A local($q,$r,$base) = @_;
1N/A if ($q eq 'NaN' || $r eq 'NaN') {
1N/A 'NaN';
1N/A } elsif ($rnd_mode eq 'trunc') {
1N/A $q; # just truncate
1N/A } else {
1N/A local($cmp) = &'bcmp(&'bmul($r,'+2'),$base);
1N/A if ( $cmp < 0 ||
1N/A ($cmp == 0 &&
1N/A ( $rnd_mode eq 'zero' ||
1N/A ($rnd_mode eq '-inf' && (substr($q,$[,1) eq '+')) ||
1N/A ($rnd_mode eq '+inf' && (substr($q,$[,1) eq '-')) ||
1N/A ($rnd_mode eq 'even' && $q =~ /[24680]$/) ||
1N/A ($rnd_mode eq 'odd' && $q =~ /[13579]$/) )) ) {
1N/A $q; # round down
1N/A } else {
1N/A &'badd($q, ((substr($q,$[,1) eq '-') ? '-1' : '+1'));
1N/A # round up
1N/A }
1N/A }
1N/A}
1N/A
1N/A# round the mantissa of $x to $scale digits
1N/Asub main'fround { #(fnum_str, scale) return fnum_str
1N/A local($x,$scale) = (&'fnorm($_[$[]),$_[$[+1]);
1N/A if ($x eq 'NaN' || $scale <= 0) {
1N/A $x;
1N/A } else {
1N/A local($xm,$xe) = split('E',$x);
1N/A if (length($xm)-1 <= $scale) {
1N/A $x;
1N/A } else {
1N/A &norm(&round(substr($xm,$[,$scale+1),
1N/A "+0".substr($xm,$[+$scale+1,1),"+10"),
1N/A $xe+length($xm)-$scale-1);
1N/A }
1N/A }
1N/A}
1N/A
1N/A# round $x at the 10 to the $scale digit place
1N/Asub main'ffround { #(fnum_str, scale) return fnum_str
1N/A local($x,$scale) = (&'fnorm($_[$[]),$_[$[+1]);
1N/A if ($x eq 'NaN') {
1N/A 'NaN';
1N/A } else {
1N/A local($xm,$xe) = split('E',$x);
1N/A if ($xe >= $scale) {
1N/A $x;
1N/A } else {
1N/A $xe = length($xm)+$xe-$scale;
1N/A if ($xe < 1) {
1N/A '+0E+0';
1N/A } elsif ($xe == 1) {
1N/A # The first substr preserves the sign, which means that
1N/A # we'll pass a non-normalized "-0" to &round when rounding
1N/A # -0.006 (for example), purely so that &round won't lose
1N/A # the sign.
1N/A &norm(&round(substr($xm,$[,1).'0',
1N/A "+0".substr($xm,$[+1,1),"+10"), $scale);
1N/A } else {
1N/A &norm(&round(substr($xm,$[,$xe),
1N/A "+0".substr($xm,$[+$xe,1),"+10"), $scale);
1N/A }
1N/A }
1N/A }
1N/A}
1N/A
1N/A# compare 2 values returns one of undef, <0, =0, >0
1N/A# returns undef if either or both input value are not numbers
1N/Asub main'fcmp #(fnum_str, fnum_str) return cond_code
1N/A{
1N/A local($x, $y) = (&'fnorm($_[$[]),&'fnorm($_[$[+1]));
1N/A if ($x eq "NaN" || $y eq "NaN") {
1N/A undef;
1N/A } else {
1N/A ord($y) <=> ord($x)
1N/A ||
1N/A ( local($xm,$xe,$ym,$ye) = split('E', $x."E$y"),
1N/A (($xe <=> $ye) * (substr($x,$[,1).'1')
1N/A || &bigint'cmp($xm,$ym))
1N/A );
1N/A }
1N/A}
1N/A
1N/A# square root by Newtons method.
1N/Asub main'fsqrt { #(fnum_str[, scale]) return fnum_str
1N/A local($x, $scale) = (&'fnorm($_[$[]), $_[$[+1]);
1N/A if ($x eq 'NaN' || $x =~ /^-/) {
1N/A 'NaN';
1N/A } elsif ($x eq '+0E+0') {
1N/A '+0E+0';
1N/A } else {
1N/A local($xm, $xe) = split('E',$x);
1N/A $scale = $div_scale if (!$scale);
1N/A $scale = length($xm)-1 if ($scale < length($xm)-1);
1N/A local($gs, $guess) = (1, sprintf("1E%+d", (length($xm)+$xe-1)/2));
1N/A while ($gs < 2*$scale) {
1N/A $guess = &'fmul(&'fadd($guess,&'fdiv($x,$guess,$gs*2)),".5");
1N/A $gs *= 2;
1N/A }
1N/A &'fround($guess, $scale);
1N/A }
1N/A}
1N/A
1N/A1;