[p5sagit/p5-mst-13.2.git] / lib / Math / BigFloat.pm

package Math::BigFloat;

use Math::BigInt;

use Exporter;  # just for use to be happy
@ISA = (Exporter);
$VERSION = '0.03';

use overload
'+'	=>	sub {new Math::BigFloat &fadd},
'-'	=>	sub {new Math::BigFloat
		       $_[2]? fsub($_[1],${$_[0]}) : fsub(${$_[0]},$_[1])},
'<=>'	=>	sub {$_[2]? fcmp($_[1],${$_[0]}) : fcmp(${$_[0]},$_[1])},
'cmp'	=>	sub {$_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])},
'*'	=>	sub {new Math::BigFloat &fmul},
'/'	=>	sub {new Math::BigFloat
		       $_[2]? scalar fdiv($_[1],${$_[0]}) :
			 scalar fdiv(${$_[0]},$_[1])},
'%'	=>	sub {new Math::BigFloat
		       $_[2]? scalar fmod($_[1],${$_[0]}) :
			 scalar fmod(${$_[0]},$_[1])},
'neg'	=>	sub {new Math::BigFloat &fneg},
'abs'	=>	sub {new Math::BigFloat &fabs},
'int'	=>	sub {new Math::BigInt &f2int},

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

sub new {
  my ($class) = shift;
  my ($foo) = fnorm(shift);
  bless \$foo, $class;
}

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

    my $minus = ($n =~ s/^([+-])// && $1 eq '-');
    $n =~ s/E//;

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

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

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

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

    return $n;
}

sub import {
  shift;
  return unless @_;
  if (@_ == 1 && $_[0] ne ':constant') {
    if ($_[0] > 0) {
      if ($VERSION < $_[0]) {
        die __PACKAGE__.": $_[0] required--this is only version $VERSION";
      }
    } else {
      die __PACKAGE__.": unknown import: $_[0]";
    }
  }
  overload::constant float => sub {Math::BigFloat->new(shift)};
}

$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;

# 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
    local($_) = @_;
    s/\s+//g;                               # strip white space
    no warnings;	# $4 and $5 below might legitimately be undefined
    if (/^([+-]?)(\d*)(\.(\d*))?([Ee]([+-]?\d+))?$/ && "$2$4" ne '') {
	&norm(($1 ? "$1$2$4" : "+$2$4"),(($4 ne '') ? $6-length($4) : $6));
    } else {
	'NaN';
    }
}

# normalize number -- for internal use
sub norm { #(mantissa, exponent) return fnum_str
    local($_, $exp) = @_;
	$exp = 0 unless defined $exp;
    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);
	}
    }
}

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

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

# 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(Math::BigInt::bmul($xm,$ym),$xe+$ye);
    }
}

# 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(Math::BigInt::badd($ym,$xm.('0' x ($xe-$ye))),$ye);
    }
}

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

# 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(Math::BigInt::bdiv($xm.('0' x $scale),$ym),
		    Math::BigInt::babs($ym)),
	    $xe-$ye-$scale);
    }
}

# modular division
#   args are dividend, divisor
sub fmod #(fnum_str, fnum_str) return fnum_str
{
    local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1]));
    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);
	if ( $xe < $ye )
	{
		$ym .= ('0' x ($ye-$xe));
	}
	else
	{
		$xm .= ('0' x ($xe-$ye));
	}
	&norm(Math::BigInt::bmod($xm,$ym));
    }
}
# round int $q based on fraction $r/$base using $rnd_mode
sub round { #(int_str, int_str, int_str) return int_str
    local($q,$r,$base) = @_;
    if ($q eq 'NaN' || $r eq 'NaN') {
	'NaN';
    } elsif ($rnd_mode eq 'trunc') {
	$q;                         # just truncate
    } else {
	local($cmp) = Math::BigInt::bcmp(Math::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 {
	    Math::BigInt::badd($q, ((substr($q,$[,1) eq '-') ? '-1' : '+1'));
				    # round up
	}
    }
}

# 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"."0" x length(substr($xm,$[+$scale+1))),
		  $xe+length($xm)-$scale-1);
	}
    }
}

# 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) {
		# The first substr preserves the sign, passing a non-
		# normalized "-0" to &round when rounding -0.006 (for
		# example), purely so &round won't lose the sign.
		&norm(&round(substr($xm,$[,1).'0',
		      "+0".substr($xm,$[+1),
		      "+1"."0" x length(substr($xm,$[+1))), $scale);
	    } else {
		&norm(&round(substr($xm,$[,$xe),
		      "+0".substr($xm,$[+$xe),
		      "+1"."0" x length(substr($xm,$[+$xe))), $scale);
	    }
	}
    }
}

# Calculate the integer part of $x
sub f2int { #(fnum_str) return inum_str
    local($x) = ${$_[$[]};
    if ($x eq 'NaN') {
	die "Attempt to take int(NaN)";
    } else {
	local($xm,$xe) = split('E',$x);
	if ($xe >= 0) {
	    $xm . '0' x $xe;
	} else {
	    $xe = length($xm)+$xe;
	    if ($xe <= 1) {
		'+0';
	    } else {
	        substr($xm,$[,$xe);
	    }
	}
    }
}

# 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 {
	local($xm,$xe,$ym,$ye) = split('E', $x."E$y");
	if ($xm eq '+0' || $ym eq '+0') {
	    return $xm <=> $ym;
	}
	if ( $xe < $ye )	# adjust the exponents to be equal
	{
		$ym .= '0' x ($ye - $xe);
		$ye = $xe;
	}
	elsif ( $ye < $xe )	# same here
	{
		$xm .= '0' x ($xe - $ye);
		$xe = $ye;
	}
	return Math::BigInt::cmp($xm,$ym);
    }
}

# 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 Math::BigFloat &fround($guess, $scale);
    }
}

1;
__END__

=head1 NAME

Math::BigFloat - Arbitrary length float math package

=head1 SYNOPSIS

  use Math::BigFloat;
  $f = Math::BigFloat->new($string);

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

=head1 DESCRIPTION

All basic math operations are overloaded if you declare your big
floats as

    $float = new Math::BigFloat "2.123123123123123123123123123123123";

=over 2

=item number format

canonical strings have the form /[+-]\d+E[+-]\d+/ .  Input values can
have embedded whitespace.

=item Error returns 'NaN'

An input parameter was "Not a Number" or divide by zero or sqrt of
negative number.

=item Division is computed to

C<max($Math::BigFloat::div_scale,length(dividend)+length(divisor))>
digits by default.
Also used for default sqrt scale.

=item Rounding is performed

according to the value of
C<$Math::BigFloat::rnd_mode>:

  trunc     truncate the value
  zero      round towards 0
  +inf      round towards +infinity (round up)
  -inf      round towards -infinity (round down)
  even      round to the nearest, .5 to the even digit
  odd       round to the nearest, .5 to the odd digit

The default is C<even> rounding.

=back

=head1 BUGS

The current version of this module is a preliminary version of the
real thing that is currently (as of perl5.002) under development.

The printf subroutine does not use the value of
C<$Math::BigFloat::rnd_mode> when rounding values for printing.
Consequently, the way to print rounded values is
to specify the number of digits both as an
argument to C<ffround> and in the C<%f> printf string,
as follows:

  printf "%.3f\n", $bigfloat->ffround(-3);

=head1 AUTHOR

Mark Biggar
Patches by John Peacock Apr 2001
=cut
Commit	Line	Data
a0d0e21e	1	package Math::BigFloat;
	2
	3	use Math::BigInt;
	4
	5	use Exporter; # just for use to be happy
	6	@ISA = (Exporter);
288d023a	7	$VERSION = '0.03';
a0d0e21e	8
a5f75d66	9	use overload
	10	'+' => sub {new Math::BigFloat &fadd},
	11	'-' => sub {new Math::BigFloat
a0d0e21e	12	$_[2]? fsub($_[1],${$_[0]}) : fsub(${$_[0]},$_[1])},
5d7098d5	13	'<=>' => sub {$_[2]? fcmp($_[1],${$_[0]}) : fcmp(${$_[0]},$_[1])},
5d7098d5	14	'cmp' => sub {$_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])},
a5f75d66	15	'*' => sub {new Math::BigFloat &fmul},
288d023a	16	'/' => sub {new Math::BigFloat
a0d0e21e	17	$_[2]? scalar fdiv($_[1],${$_[0]}) :
a0d0e21e	18	scalar fdiv(${$_[0]},$_[1])},
288d023a	19	'%' => sub {new Math::BigFloat
	20	$_[2]? scalar fmod($_[1],${$_[0]}) :
	21	scalar fmod(${$_[0]},$_[1])},
a5f75d66	22	'neg' => sub {new Math::BigFloat &fneg},
a5f75d66	23	'abs' => sub {new Math::BigFloat &fabs},
f216259d	24	'int' => sub {new Math::BigInt &f2int},
a0d0e21e	25
	26	qw(
	27	"" stringify
	28	0+ numify) # Order of arguments unsignificant
a5f75d66	29	;
a0d0e21e	30
a0d0e21e	31	sub new {
a5f75d66	32	my ($class) = shift;
a5f75d66	33	my ($foo) = fnorm(shift);
a5f75d66	34	bless \$foo, $class;
a0d0e21e	35	}
0e8b9368	36
a0d0e21e	37	sub numify { 0 + "${$_[0]}" } # Not needed, additional overhead
	38	# comparing to direct compilation based on
	39	# stringify
	40	sub stringify {
	41	my $n = ${$_[0]};
	42
9b599b2a	43	my $minus = ($n =~ s/^([+-])// && $1 eq '-');
a0d0e21e	44	$n =~ s/E//;
	45
	46	$n =~ s/([-+]\d+)$//;
	47
	48	my $e = $1;
	49	my $ln = length($n);
	50
288d023a	51	if ( defined $e )
	52	{
	53	if ($e > 0) {
	54	$n .= "0" x $e . '.';
	55	} elsif (abs($e) < $ln) {
	56	substr($n, $ln + $e, 0) = '.';
	57	} else {
	58	$n = '.' . ("0" x (abs($e) - $ln)) . $n;
	59	}
a0d0e21e	60	}
9b599b2a	61	$n = "-$n" if $minus;
a0d0e21e	62
	63	# 1 while $n =~ s/(.*\d)(\d\d\d)/$1,$2/;
	64
	65	return $n;
	66	}
	67
f216259d	68	sub import {
	69	shift;
	70	return unless @_;
82cf049f	71	if (@_ == 1 && $_[0] ne ':constant') {
	72	if ($_[0] > 0) {
	73	if ($VERSION < $_[0]) {
	74	die __PACKAGE__.": $_[0] required--this is only version $VERSION";
	75	}
	76	} else {
	77	die __PACKAGE__.": unknown import: $_[0]";
	78	}
	79	}
f216259d	80	overload::constant float => sub {Math::BigFloat->new(shift)};
	81	}
	82
a0d0e21e	83	$div_scale = 40;
	84
	85	# Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
	86
	87	$rnd_mode = 'even';
	88
	89	sub fadd; sub fsub; sub fmul; sub fdiv;
	90	sub fneg; sub fabs; sub fcmp;
	91	sub fround; sub ffround;
	92	sub fnorm; sub fsqrt;
	93
a0d0e21e	94	# Convert a number to canonical string form.
	95	# Takes something that looks like a number and converts it to
	96	# the form /^[+-]\d+E[+-]\d+$/.
	97	sub fnorm { #(string) return fnum_str
	98	local($_) = @_;
	99	s/\s+//g; # strip white space
db376a24	100	no warnings; # $4 and $5 below might legitimately be undefined
a0d0e21e	101	if (/^([+-]?)(\d)(\.(\d))?([Ee]([+-]?\d+))?$/ && "$2$4" ne '') {
	102	&norm(($1 ? "$1$2$4" : "+$2$4"),(($4 ne '') ? $6-length($4) : $6));
	103	} else {
	104	'NaN';
	105	}
	106	}
	107
	108	# normalize number -- for internal use
	109	sub norm { #(mantissa, exponent) return fnum_str
	110	local($_, $exp) = @_;
3deb277d	111	$exp = 0 unless defined $exp;
a0d0e21e	112	if ($_ eq 'NaN') {
	113	'NaN';
	114	} else {
	115	s/^([+-])0+/$1/; # strip leading zeros
	116	if (length($_) == 1) {
	117	'+0E+0';
	118	} else {
	119	$exp += length($1) if (s/(0+)$//); # strip trailing zeros
	120	sprintf("%sE%+ld", $_, $exp);
	121	}
	122	}
	123	}
	124
	125	# negation
	126	sub fneg { #(fnum_str) return fnum_str
	127	local($_) = fnorm($_[$[]);
	128	vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0E+0'; # flip sign
	129	s/^H/N/;
	130	$_;
	131	}
	132
	133	# absolute value
	134	sub fabs { #(fnum_str) return fnum_str
	135	local($_) = fnorm($_[$[]);
	136	s/^-/+/; # mash sign
	137	$_;
	138	}
	139
	140	# multiplication
	141	sub fmul { #(fnum_str, fnum_str) return fnum_str
	142	local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1]));
	143	if ($x eq 'NaN' \|\| $y eq 'NaN') {
	144	'NaN';
	145	} else {
	146	local($xm,$xe) = split('E',$x);
	147	local($ym,$ye) = split('E',$y);
	148	&norm(Math::BigInt::bmul($xm,$ym),$xe+$ye);
	149	}
	150	}
a5f75d66	151
a0d0e21e	152	# addition
	153	sub fadd { #(fnum_str, fnum_str) return fnum_str
	154	local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1]));
	155	if ($x eq 'NaN' \|\| $y eq 'NaN') {
	156	'NaN';
	157	} else {
	158	local($xm,$xe) = split('E',$x);
	159	local($ym,$ye) = split('E',$y);
	160	($xm,$xe,$ym,$ye) = ($ym,$ye,$xm,$xe) if ($xe < $ye);
	161	&norm(Math::BigInt::badd($ym,$xm.('0' x ($xe-$ye))),$ye);
	162	}
	163	}
	164
	165	# subtraction
	166	sub fsub { #(fnum_str, fnum_str) return fnum_str
288d023a	167	fadd($_[$[],fneg($_[$[+1]));
a0d0e21e	168	}
	169
	170	# division
	171	# args are dividend, divisor, scale (optional)
	172	# result has at most max(scale, length(dividend), length(divisor)) digits
	173	sub fdiv #(fnum_str, fnum_str[,scale]) return fnum_str
	174	{
	175	local($x,$y,$scale) = (fnorm($_[$[]),fnorm($_[$[+1]),$_[$[+2]);
	176	if ($x eq 'NaN' \|\| $y eq 'NaN' \|\| $y eq '+0E+0') {
	177	'NaN';
	178	} else {
	179	local($xm,$xe) = split('E',$x);
	180	local($ym,$ye) = split('E',$y);
	181	$scale = $div_scale if (!$scale);
	182	$scale = length($xm)-1 if (length($xm)-1 > $scale);
	183	$scale = length($ym)-1 if (length($ym)-1 > $scale);
	184	$scale = $scale + length($ym) - length($xm);
5d7098d5	185	&norm(&round(Math::BigInt::bdiv($xm.('0' x $scale),$ym),
5d7098d5	186	Math::BigInt::babs($ym)),
a0d0e21e	187	$xe-$ye-$scale);
	188	}
	189	}
a5f75d66	190
288d023a	191	# modular division
	192	# args are dividend, divisor
	193	sub fmod #(fnum_str, fnum_str) return fnum_str
	194	{
	195	local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1]));
	196	if ($x eq 'NaN' \|\| $y eq 'NaN' \|\| $y eq '+0E+0') {
	197	'NaN';
	198	} else {
	199	local($xm,$xe) = split('E',$x);
	200	local($ym,$ye) = split('E',$y);
	201	if ( $xe < $ye )
	202	{
	203	$ym .= ('0' x ($ye-$xe));
	204	}
	205	else
	206	{
	207	$xm .= ('0' x ($xe-$ye));
	208	}
	209	&norm(Math::BigInt::bmod($xm,$ym));
	210	}
	211	}
a0d0e21e	212	# round int $q based on fraction $r/$base using $rnd_mode
	213	sub round { #(int_str, int_str, int_str) return int_str
	214	local($q,$r,$base) = @_;
	215	if ($q eq 'NaN' \|\| $r eq 'NaN') {
	216	'NaN';
	217	} elsif ($rnd_mode eq 'trunc') {
	218	$q; # just truncate
	219	} else {
	220	local($cmp) = Math::BigInt::bcmp(Math::BigInt::bmul($r,'+2'),$base);
	221	if ( $cmp < 0 \|\|
3deb277d	222	($cmp == 0 && (
3deb277d	223	($rnd_mode eq 'zero' ) \|\|
a0d0e21e	224	($rnd_mode eq '-inf' && (substr($q,$[,1) eq '+')) \|\|
a0d0e21e	225	($rnd_mode eq '+inf' && (substr($q,$[,1) eq '-')) \|\|
288d023a	226	($rnd_mode eq 'even' && $q =~ /[24680]$/ ) \|\|
	227	($rnd_mode eq 'odd' && $q =~ /[13579]$/ ) )
	228	)
3deb277d	229	) {
a0d0e21e	230	$q; # round down
	231	} else {
	232	Math::BigInt::badd($q, ((substr($q,$[,1) eq '-') ? '-1' : '+1'));
	233	# round up
	234	}
	235	}
	236	}
	237
	238	# round the mantissa of $x to $scale digits
	239	sub fround { #(fnum_str, scale) return fnum_str
	240	local($x,$scale) = (fnorm($_[$[]),$_[$[+1]);
	241	if ($x eq 'NaN' \|\| $scale <= 0) {
	242	$x;
	243	} else {
	244	local($xm,$xe) = split('E',$x);
	245	if (length($xm)-1 <= $scale) {
	246	$x;
	247	} else {
288d023a	248	&norm(&round(substr($xm,$[,$scale+1),
288d023a	249	"+0".substr($xm,$[+$scale+1),"+1"."0" x length(substr($xm,$[+$scale+1))),
a0d0e21e	250	$xe+length($xm)-$scale-1);
	251	}
	252	}
	253	}
a5f75d66	254
a0d0e21e	255	# round $x at the 10 to the $scale digit place
	256	sub ffround { #(fnum_str, scale) return fnum_str
	257	local($x,$scale) = (fnorm($_[$[]),$_[$[+1]);
	258	if ($x eq 'NaN') {
	259	'NaN';
	260	} else {
	261	local($xm,$xe) = split('E',$x);
	262	if ($xe >= $scale) {
	263	$x;
	264	} else {
	265	$xe = length($xm)+$xe-$scale;
	266	if ($xe < 1) {
	267	'+0E+0';
	268	} elsif ($xe == 1) {
5d7098d5	269	# The first substr preserves the sign, passing a non-
	270	# normalized "-0" to &round when rounding -0.006 (for
	271	# example), purely so &round won't lose the sign.
	272	&norm(&round(substr($xm,$[,1).'0',
288d023a	273	"+0".substr($xm,$[+1),
288d023a	274	"+1"."0" x length(substr($xm,$[+1))), $scale);
a0d0e21e	275	} else {
a0d0e21e	276	&norm(&round(substr($xm,$[,$xe),
288d023a	277	"+0".substr($xm,$[+$xe),
288d023a	278	"+1"."0" x length(substr($xm,$[+$xe))), $scale);
a0d0e21e	279	}
	280	}
	281	}
	282	}
f216259d	283
	284	# Calculate the integer part of $x
	285	sub f2int { #(fnum_str) return inum_str
	286	local($x) = ${$_[$[]};
	287	if ($x eq 'NaN') {
	288	die "Attempt to take int(NaN)";
	289	} else {
	290	local($xm,$xe) = split('E',$x);
	291	if ($xe >= 0) {
	292	$xm . '0' x $xe;
	293	} else {
	294	$xe = length($xm)+$xe;
	295	if ($xe <= 1) {
	296	'+0';
	297	} else {
	298	substr($xm,$[,$xe);
	299	}
	300	}
	301	}
	302	}
288d023a	303
a0d0e21e	304	# compare 2 values returns one of undef, <0, =0, >0
	305	# returns undef if either or both input value are not numbers
	306	sub fcmp #(fnum_str, fnum_str) return cond_code
	307	{
	308	local($x, $y) = (fnorm($_[$[]),fnorm($_[$[+1]));
	309	if ($x eq "NaN" \|\| $y eq "NaN") {
	310	undef;
	311	} else {
5dce09b1	312	local($xm,$xe,$ym,$ye) = split('E', $x."E$y");
	313	if ($xm eq '+0' \|\| $ym eq '+0') {
	314	return $xm <=> $ym;
	315	}
3deb277d	316	if ( $xe < $ye ) # adjust the exponents to be equal
	317	{
	318	$ym .= '0' x ($ye - $xe);
	319	$ye = $xe;
	320	}
	321	elsif ( $ye < $xe ) # same here
	322	{
	323	$xm .= '0' x ($xe - $ye);
	324	$xe = $ye;
	325	}
	326	return Math::BigInt::cmp($xm,$ym);
a0d0e21e	327	}
a0d0e21e	328	}
a5f75d66	329
a0d0e21e	330	# square root by Newtons method.
	331	sub fsqrt { #(fnum_str[, scale]) return fnum_str
	332	local($x, $scale) = (fnorm($_[$[]), $_[$[+1]);
	333	if ($x eq 'NaN' \|\| $x =~ /^-/) {
	334	'NaN';
	335	} elsif ($x eq '+0E+0') {
	336	'+0E+0';
	337	} else {
	338	local($xm, $xe) = split('E',$x);
	339	$scale = $div_scale if (!$scale);
	340	$scale = length($xm)-1 if ($scale < length($xm)-1);
	341	local($gs, $guess) = (1, sprintf("1E%+d", (length($xm)+$xe-1)/2));
	342	while ($gs < 2*$scale) {
	343	$guess = fmul(fadd($guess,fdiv($x,$guess,$gs*2)),".5");
	344	$gs *= 2;
	345	}
a5f75d66	346	new Math::BigFloat &fround($guess, $scale);
a0d0e21e	347	}
	348	}
	349
	350	1;
a5f75d66	351	__END__
	352
	353	=head1 NAME
	354
	355	Math::BigFloat - Arbitrary length float math package
	356
	357	=head1 SYNOPSIS
	358
a2008d6d	359	use Math::BigFloat;
a5f75d66	360	$f = Math::BigFloat->new($string);
	361
	362	$f->fadd(NSTR) return NSTR addition
	363	$f->fsub(NSTR) return NSTR subtraction
	364	$f->fmul(NSTR) return NSTR multiplication
	365	$f->fdiv(NSTR[,SCALE]) returns NSTR division to SCALE places
288d023a	366	$f->fmod(NSTR) returns NSTR modular remainder
a5f75d66	367	$f->fneg() return NSTR negation
	368	$f->fabs() return NSTR absolute value
	369	$f->fcmp(NSTR) return CODE compare undef,<0,=0,>0
	370	$f->fround(SCALE) return NSTR round to SCALE digits
	371	$f->ffround(SCALE) return NSTR round at SCALEth place
	372	$f->fnorm() return (NSTR) normalize
	373	$f->fsqrt([SCALE]) return NSTR sqrt to SCALE places
	374
	375	=head1 DESCRIPTION
	376
	377	All basic math operations are overloaded if you declare your big
	378	floats as
	379
	380	$float = new Math::BigFloat "2.123123123123123123123123123123123";
	381
	382	=over 2
	383
	384	=item number format
	385
	386	canonical strings have the form /[+-]\d+E[+-]\d+/ . Input values can
db696efe	387	have embedded whitespace.
a5f75d66	388
	389	=item Error returns 'NaN'
	390
	391	An input parameter was "Not a Number" or divide by zero or sqrt of
	392	negative number.
	393
288d023a	394	=item Division is computed to
a5f75d66	395
5d7098d5	396	C<max($Math::BigFloat::div_scale,length(dividend)+length(divisor))>
5d7098d5	397	digits by default.
a5f75d66	398	Also used for default sqrt scale.
a5f75d66	399
5d7098d5	400	=item Rounding is performed
	401
	402	according to the value of
	403	C<$Math::BigFloat::rnd_mode>:
	404
	405	trunc truncate the value
	406	zero round towards 0
	407	+inf round towards +infinity (round up)
	408	-inf round towards -infinity (round down)
	409	even round to the nearest, .5 to the even digit
	410	odd round to the nearest, .5 to the odd digit
	411
	412	The default is C<even> rounding.
	413
a5f75d66	414	=back
	415
	416	=head1 BUGS
	417
	418	The current version of this module is a preliminary version of the
	419	real thing that is currently (as of perl5.002) under development.
	420
5d7098d5	421	The printf subroutine does not use the value of
	422	C<$Math::BigFloat::rnd_mode> when rounding values for printing.
	423	Consequently, the way to print rounded values is
	424	to specify the number of digits both as an
	425	argument to C<ffround> and in the C<%f> printf string,
	426	as follows:
	427
	428	printf "%.3f\n", $bigfloat->ffround(-3);
	429
a5f75d66	430	=head1 AUTHOR
	431
	432	Mark Biggar
288d023a	433	Patches by John Peacock Apr 2001
a5f75d66	434	=cut