Upgrade to Locale::Maketext 1.07.
[p5sagit/p5-mst-13.2.git] / lib / Math / BigFloat.pm
CommitLineData
13a12e00 1package Math::BigFloat;
2
3#
d614cd8b 4# Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
13a12e00 5#
6
58cde26e 7# The following hash values are internally used:
8# _e: exponent (BigInt)
9# _m: mantissa (absolute BigInt)
990fb837 10# sign: +,-,+inf,-inf, or "NaN" if not a number
58cde26e 11# _a: accuracy
12# _p: precision
0716bf9b 13# _f: flags, used to signal MBI not to touch our private parts
58cde26e 14
b282a552 15$VERSION = '1.42';
58cde26e 16require 5.005;
3a427a11 17
18require Exporter;
2ab5f49d 19@ISA = qw(Exporter Math::BigInt);
394e6ffb 20
58cde26e 21use strict;
b282a552 22# $_trap_inf and $_trap_nan are internal and should never be accessed from the outside
23use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode
24 $upgrade $downgrade $_trap_nan $_trap_inf/;
58cde26e 25my $class = "Math::BigFloat";
a0d0e21e 26
a5f75d66 27use overload
bd05a461 28'<=>' => sub { $_[2] ?
29 ref($_[0])->bcmp($_[1],$_[0]) :
30 ref($_[0])->bcmp($_[0],$_[1])},
0716bf9b 31'int' => sub { $_[0]->as_number() }, # 'trunc' to bigint
a5f75d66 32;
a0d0e21e 33
0716bf9b 34##############################################################################
990fb837 35# global constants, flags and assorted stuff
0716bf9b 36
990fb837 37# the following are public, but their usage is not recommended. Use the
38# accessor methods instead.
58cde26e 39
ee15d750 40# class constants, use Class->constant_name() to access
41$round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
42$accuracy = undef;
43$precision = undef;
44$div_scale = 40;
58cde26e 45
b3abae2a 46$upgrade = undef;
47$downgrade = undef;
8f675a64 48my $MBI = 'Math::BigInt'; # the package we are using for our private parts
49 # changable by use Math::BigFloat with => 'package'
b3abae2a 50
990fb837 51# the following are private and not to be used from the outside:
52
b282a552 53sub MB_NEVER_ROUND () { 0x0001; }
990fb837 54
55# are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
56$_trap_nan = 0;
57# the same for infs
58$_trap_inf = 0;
59
60# constant for easier life
61my $nan = 'NaN';
62
63my $IMPORT = 0; # was import() called yet?
64 # used to make require work
65
66# some digits of accuracy for blog(undef,10); which we use in blog() for speed
67my $LOG_10 =
68 '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
69my $LOG_10_A = length($LOG_10)-1;
70# ditto for log(2)
71my $LOG_2 =
72 '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
73my $LOG_2_A = length($LOG_2)-1;
74
027dc388 75##############################################################################
76# the old code had $rnd_mode, so we need to support it, too
77
027dc388 78sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
79sub FETCH { return $round_mode; }
80sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
81
56b9c951 82BEGIN
990fb837 83 {
84 # when someone set's $rnd_mode, we catch this and check the value to see
85 # whether it is valid or not.
86 $rnd_mode = 'even'; tie $rnd_mode, 'Math::BigFloat';
56b9c951 87 }
027dc388 88
89##############################################################################
90
574bacfe 91# in case we call SUPER::->foo() and this wants to call modify()
92# sub modify () { 0; }
93
58cde26e 94{
ee15d750 95 # valid method aliases for AUTOLOAD
58cde26e 96 my %methods = map { $_ => 1 }
97 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
b3abae2a 98 fint facmp fcmp fzero fnan finf finc fdec flog ffac
990fb837 99 fceil ffloor frsft flsft fone flog froot
ee15d750 100 /;
61f5c3f5 101 # valid method's that can be hand-ed up (for AUTOLOAD)
ee15d750 102 my %hand_ups = map { $_ => 1 }
103 qw / is_nan is_inf is_negative is_positive
990fb837 104 accuracy precision div_scale round_mode fneg fabs fnot
28df3e88 105 objectify upgrade downgrade
13a12e00 106 bone binf bnan bzero
58cde26e 107 /;
108
ee15d750 109 sub method_alias { return exists $methods{$_[0]||''}; }
110 sub method_hand_up { return exists $hand_ups{$_[0]||''}; }
a0d0e21e 111}
0e8b9368 112
58cde26e 113##############################################################################
114# constructors
a0d0e21e 115
58cde26e 116sub new
117 {
118 # create a new BigFloat object from a string or another bigfloat object.
119 # _e: exponent
120 # _m: mantissa
121 # sign => sign (+/-), or "NaN"
a0d0e21e 122
61f5c3f5 123 my ($class,$wanted,@r) = @_;
b3abae2a 124
61f5c3f5 125 # avoid numify-calls by not using || on $wanted!
126 return $class->bzero() if !defined $wanted; # default to 0
127 return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
a0d0e21e 128
990fb837 129 $class->import() if $IMPORT == 0; # make require work
130
58cde26e 131 my $self = {}; bless $self, $class;
b22b3e31 132 # shortcut for bigints and its subclasses
0716bf9b 133 if ((ref($wanted)) && (ref($wanted) ne $class))
58cde26e 134 {
0716bf9b 135 $self->{_m} = $wanted->as_number(); # get us a bigint copy
56b9c951 136 $self->{_e} = $MBI->bzero();
58cde26e 137 $self->{_m}->babs();
138 $self->{sign} = $wanted->sign();
0716bf9b 139 return $self->bnorm();
58cde26e 140 }
141 # got string
142 # handle '+inf', '-inf' first
ee15d750 143 if ($wanted =~ /^[+-]?inf$/)
58cde26e 144 {
28df3e88 145 return $downgrade->new($wanted) if $downgrade;
146
56b9c951 147 $self->{_e} = $MBI->bzero();
148 $self->{_m} = $MBI->bzero();
58cde26e 149 $self->{sign} = $wanted;
ee15d750 150 $self->{sign} = '+inf' if $self->{sign} eq 'inf';
0716bf9b 151 return $self->bnorm();
58cde26e 152 }
153 #print "new string '$wanted'\n";
b282a552 154
58cde26e 155 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
156 if (!ref $mis)
157 {
990fb837 158 if ($_trap_nan)
159 {
160 require Carp;
161 Carp::croak ("$wanted is not a number initialized to $class");
162 }
28df3e88 163
164 return $downgrade->bnan() if $downgrade;
165
56b9c951 166 $self->{_e} = $MBI->bzero();
167 $self->{_m} = $MBI->bzero();
58cde26e 168 $self->{sign} = $nan;
169 }
170 else
171 {
172 # make integer from mantissa by adjusting exp, then convert to bigint
61f5c3f5 173 # undef,undef to signal MBI that we don't need no bloody rounding
56b9c951 174 $self->{_e} = $MBI->new("$$es$$ev",undef,undef); # exponent
175 $self->{_m} = $MBI->new("$$miv$$mfv",undef,undef); # create mant.
b282a552 176
177 # this is to prevent automatically rounding when MBI's globals are set
178 $self->{_m}->{_f} = MB_NEVER_ROUND;
179 $self->{_e}->{_f} = MB_NEVER_ROUND;
180
58cde26e 181 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
b282a552 182 $self->{_e}->bsub( $MBI->new(CORE::length($$mfv),undef,undef))
183 if CORE::length($$mfv) != 0;
027dc388 184 $self->{sign} = $$mis;
b282a552 185
186 #print "$$miv$$mfv $$es$$ev\n";
187
188 # we can only have trailing zeros on the mantissa of $$mfv eq ''
189 if (CORE::length($$mfv) == 0)
190 {
191 my $zeros = $self->{_m}->_trailing_zeros(); # correct for trailing zeros
192 if ($zeros != 0)
193 {
194 $self->{_m}->brsft($zeros,10); $self->{_e}->badd($MBI->new($zeros));
195 }
196 }
197# else
198# {
199 # for something like 0Ey, set y to 1, and -0 => +0
200 $self->{sign} = '+', $self->{_e}->bone() if $self->{_m}->is_zero();
201# }
202 return $self->round(@r) if !$downgrade;
58cde26e 203 }
28df3e88 204 # if downgrade, inf, NaN or integers go down
205
206 if ($downgrade && $self->{_e}->{sign} eq '+')
207 {
990fb837 208 #print "downgrading $$miv$$mfv"."E$$es$$ev";
28df3e88 209 if ($self->{_e}->is_zero())
210 {
211 $self->{_m}->{sign} = $$mis; # negative if wanted
212 return $downgrade->new($self->{_m});
213 }
8df1e0a2 214 return $downgrade->new($self->bsstr());
28df3e88 215 }
990fb837 216 #print "mbf new $self->{sign} $self->{_m} e $self->{_e} ",ref($self),"\n";
217 $self->bnorm()->round(@r); # first normalize, then round
58cde26e 218 }
a0d0e21e 219
13a12e00 220sub _bnan
58cde26e 221 {
990fb837 222 # used by parent class bone() to initialize number to NaN
58cde26e 223 my $self = shift;
990fb837 224
225 if ($_trap_nan)
226 {
227 require Carp;
228 my $class = ref($self);
229 Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
230 }
231
232 $IMPORT=1; # call our import only once
56b9c951 233 $self->{_m} = $MBI->bzero();
234 $self->{_e} = $MBI->bzero();
58cde26e 235 }
a0d0e21e 236
13a12e00 237sub _binf
58cde26e 238 {
990fb837 239 # used by parent class bone() to initialize number to +-inf
58cde26e 240 my $self = shift;
990fb837 241
242 if ($_trap_inf)
243 {
244 require Carp;
245 my $class = ref($self);
246 Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
247 }
248
249 $IMPORT=1; # call our import only once
56b9c951 250 $self->{_m} = $MBI->bzero();
251 $self->{_e} = $MBI->bzero();
58cde26e 252 }
a0d0e21e 253
13a12e00 254sub _bone
574bacfe 255 {
13a12e00 256 # used by parent class bone() to initialize number to 1
574bacfe 257 my $self = shift;
990fb837 258 $IMPORT=1; # call our import only once
56b9c951 259 $self->{_m} = $MBI->bone();
260 $self->{_e} = $MBI->bzero();
574bacfe 261 }
262
13a12e00 263sub _bzero
58cde26e 264 {
990fb837 265 # used by parent class bone() to initialize number to 0
58cde26e 266 my $self = shift;
990fb837 267 $IMPORT=1; # call our import only once
56b9c951 268 $self->{_m} = $MBI->bzero();
269 $self->{_e} = $MBI->bone();
58cde26e 270 }
271
9393ace2 272sub isa
273 {
274 my ($self,$class) = @_;
56b9c951 275 return if $class =~ /^Math::BigInt/; # we aren't one of these
276 UNIVERSAL::isa($self,$class);
9393ace2 277 }
278
8f675a64 279sub config
280 {
281 # return (later set?) configuration data as hash ref
282 my $class = shift || 'Math::BigFloat';
283
990fb837 284 my $cfg = $class->SUPER::config(@_);
8f675a64 285
990fb837 286 # now we need only to override the ones that are different from our parent
8f675a64 287 $cfg->{class} = $class;
288 $cfg->{with} = $MBI;
8f675a64 289 $cfg;
290 }
291
58cde26e 292##############################################################################
293# string conversation
294
295sub bstr
296 {
297 # (ref to BFLOAT or num_str ) return num_str
298 # Convert number from internal format to (non-scientific) string format.
299 # internal format is always normalized (no leading zeros, "-0" => "+0")
ee15d750 300 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 301
574bacfe 302 if ($x->{sign} !~ /^[+-]$/)
58cde26e 303 {
574bacfe 304 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
305 return 'inf'; # +inf
58cde26e 306 }
c38b2de2 307
574bacfe 308 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
309
c38b2de2 310 # $x is zero?
311 my $not_zero = !($x->{sign} eq '+' && $x->{_m}->is_zero());
574bacfe 312 if ($not_zero)
58cde26e 313 {
574bacfe 314 $es = $x->{_m}->bstr();
315 $len = CORE::length($es);
c38b2de2 316 my $e = $x->{_e}->numify();
317 if ($e < 0)
58cde26e 318 {
c38b2de2 319 $dot = '';
320 # if _e is bigger than a scalar, the following will blow your memory
321 if ($e <= -$len)
574bacfe 322 {
c38b2de2 323 #print "style: 0.xxxx\n";
324 my $r = abs($e) - $len;
325 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
574bacfe 326 }
327 else
328 {
c38b2de2 329 #print "insert '.' at $e in '$es'\n";
330 substr($es,$e,0) = '.'; $cad = $x->{_e};
574bacfe 331 }
82cf049f 332 }
c38b2de2 333 elsif ($e > 0)
334 {
335 # expand with zeros
336 $es .= '0' x $e; $len += $e; $cad = 0;
337 }
574bacfe 338 } # if not zero
c38b2de2 339 $es = '-'.$es if $x->{sign} eq '-';
340 # if set accuracy or precision, pad with zeros on the right side
574bacfe 341 if ((defined $x->{_a}) && ($not_zero))
342 {
343 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
344 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
345 $zeros = $x->{_a} - $len if $cad != $len;
574bacfe 346 $es .= $dot.'0' x $zeros if $zeros > 0;
82cf049f 347 }
c38b2de2 348 elsif ((($x->{_p} || 0) < 0))
58cde26e 349 {
574bacfe 350 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
351 my $zeros = -$x->{_p} + $cad;
574bacfe 352 $es .= $dot.'0' x $zeros if $zeros > 0;
58cde26e 353 }
56b9c951 354 $es;
82cf049f 355 }
f216259d 356
58cde26e 357sub bsstr
358 {
359 # (ref to BFLOAT or num_str ) return num_str
360 # Convert number from internal format to scientific string format.
361 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
ee15d750 362 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
a0d0e21e 363
574bacfe 364 if ($x->{sign} !~ /^[+-]$/)
365 {
366 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
367 return 'inf'; # +inf
368 }
091c87b1 369 # do $esign, because we need '1e+1', since $x->{_e}->bstr() misses the +
56d9de68 370 my $esign = $x->{_e}->{sign}; $esign = '' if $esign eq '-';
371 my $sep = 'e'.$esign;
372 my $sign = $x->{sign}; $sign = '' if $sign eq '+';
373 $sign . $x->{_m}->bstr() . $sep . $x->{_e}->bstr();
58cde26e 374 }
375
376sub numify
377 {
378 # Make a number from a BigFloat object
b282a552 379 # simple return a string and let Perl's atoi()/atof() handle the rest
380 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
56b9c951 381 $x->bsstr();
58cde26e 382 }
a0d0e21e 383
58cde26e 384##############################################################################
385# public stuff (usually prefixed with "b")
386
574bacfe 387# tels 2001-08-04
b282a552 388# XXX TODO this must be overwritten and return NaN for non-integer values
574bacfe 389# band(), bior(), bxor(), too
58cde26e 390#sub bnot
391# {
392# $class->SUPER::bnot($class,@_);
393# }
394
395sub bcmp
396 {
397 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
f9a08e12 398
399 # set up parameters
400 my ($self,$x,$y) = (ref($_[0]),@_);
401 # objectify is costly, so avoid it
402 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
403 {
404 ($self,$x,$y) = objectify(2,@_);
405 }
58cde26e 406
56d9de68 407 return $upgrade->bcmp($x,$y) if defined $upgrade &&
408 ((!$x->isa($self)) || (!$y->isa($self)));
409
0716bf9b 410 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
411 {
412 # handle +-inf and NaN
413 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
414 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
415 return +1 if $x->{sign} eq '+inf';
416 return -1 if $x->{sign} eq '-inf';
417 return -1 if $y->{sign} eq '+inf';
b3abae2a 418 return +1;
0716bf9b 419 }
420
421 # check sign for speed first
574bacfe 422 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
58cde26e 423 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
424
574bacfe 425 # shortcut
426 my $xz = $x->is_zero();
427 my $yz = $y->is_zero();
428 return 0 if $xz && $yz; # 0 <=> 0
429 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
430 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
58cde26e 431
432 # adjust so that exponents are equal
bd05a461 433 my $lxm = $x->{_m}->length();
434 my $lym = $y->{_m}->length();
28df3e88 435 # the numify somewhat limits our length, but makes it much faster
436 my $lx = $lxm + $x->{_e}->numify();
437 my $ly = $lym + $y->{_e}->numify();
438 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
bd05a461 439 return $l <=> 0 if $l != 0;
58cde26e 440
bd05a461 441 # lengths (corrected by exponent) are equal
28df3e88 442 # so make mantissa equal length by padding with zero (shift left)
bd05a461 443 my $diff = $lxm - $lym;
444 my $xm = $x->{_m}; # not yet copy it
445 my $ym = $y->{_m};
446 if ($diff > 0)
447 {
448 $ym = $y->{_m}->copy()->blsft($diff,10);
449 }
450 elsif ($diff < 0)
451 {
452 $xm = $x->{_m}->copy()->blsft(-$diff,10);
453 }
28df3e88 454 my $rc = $xm->bacmp($ym);
58cde26e 455 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
b3abae2a 456 $rc <=> 0;
58cde26e 457 }
458
459sub bacmp
460 {
461 # Compares 2 values, ignoring their signs.
462 # Returns one of undef, <0, =0, >0. (suitable for sort)
f9a08e12 463
464 # set up parameters
465 my ($self,$x,$y) = (ref($_[0]),@_);
466 # objectify is costly, so avoid it
467 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
468 {
469 ($self,$x,$y) = objectify(2,@_);
470 }
ee15d750 471
56d9de68 472 return $upgrade->bacmp($x,$y) if defined $upgrade &&
473 ((!$x->isa($self)) || (!$y->isa($self)));
474
ee15d750 475 # handle +-inf and NaN's
abcfbf51 476 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
ee15d750 477 {
478 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
479 return 0 if ($x->is_inf() && $y->is_inf());
480 return 1 if ($x->is_inf() && !$y->is_inf());
b3abae2a 481 return -1;
ee15d750 482 }
483
484 # shortcut
485 my $xz = $x->is_zero();
486 my $yz = $y->is_zero();
487 return 0 if $xz && $yz; # 0 <=> 0
488 return -1 if $xz && !$yz; # 0 <=> +y
489 return 1 if $yz && !$xz; # +x <=> 0
490
491 # adjust so that exponents are equal
492 my $lxm = $x->{_m}->length();
493 my $lym = $y->{_m}->length();
28df3e88 494 # the numify somewhat limits our length, but makes it much faster
495 my $lx = $lxm + $x->{_e}->numify();
496 my $ly = $lym + $y->{_e}->numify();
394e6ffb 497 my $l = $lx - $ly;
ee15d750 498 return $l <=> 0 if $l != 0;
58cde26e 499
ee15d750 500 # lengths (corrected by exponent) are equal
394e6ffb 501 # so make mantissa equal-length by padding with zero (shift left)
ee15d750 502 my $diff = $lxm - $lym;
503 my $xm = $x->{_m}; # not yet copy it
504 my $ym = $y->{_m};
505 if ($diff > 0)
506 {
507 $ym = $y->{_m}->copy()->blsft($diff,10);
508 }
509 elsif ($diff < 0)
510 {
511 $xm = $x->{_m}->copy()->blsft(-$diff,10);
512 }
28df3e88 513 $xm->bacmp($ym) <=> 0;
58cde26e 514 }
a0d0e21e 515
58cde26e 516sub badd
517 {
518 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
519 # return result as BFLOAT
f9a08e12 520
521 # set up parameters
522 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
523 # objectify is costly, so avoid it
524 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
525 {
526 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
527 }
58cde26e 528
574bacfe 529 # inf and NaN handling
530 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
531 {
532 # NaN first
533 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
13a12e00 534 # inf handling
574bacfe 535 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
536 {
13a12e00 537 # +inf++inf or -inf+-inf => same, rest is NaN
538 return $x if $x->{sign} eq $y->{sign};
539 return $x->bnan();
574bacfe 540 }
56b9c951 541 # +-inf + something => +inf; something +-inf => +-inf
574bacfe 542 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
543 return $x;
544 }
545
8f675a64 546 return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
547 ((!$x->isa($self)) || (!$y->isa($self)));
548
58cde26e 549 # speed: no add for 0+y or x+0
28df3e88 550 return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
58cde26e 551 if ($x->is_zero()) # 0+y
552 {
553 # make copy, clobbering up x (modify in place!)
554 $x->{_e} = $y->{_e}->copy();
555 $x->{_m} = $y->{_m}->copy();
556 $x->{sign} = $y->{sign} || $nan;
557 return $x->round($a,$p,$r,$y);
a0d0e21e 558 }
58cde26e 559
560 # take lower of the two e's and adapt m1 to it to match m2
28df3e88 561 my $e = $y->{_e};
091c87b1 562 $e = $MBI->bzero() if !defined $e; # if no BFLOAT ?
563 $e = $e->copy(); # make copy (didn't do it yet)
564 $e->bsub($x->{_e}); # Ye - Xe
58cde26e 565 my $add = $y->{_m}->copy();
091c87b1 566 if ($e->{sign} eq '-') # < 0
58cde26e 567 {
091c87b1 568 $x->{_e} += $e; # need the sign of e
569 $x->{_m}->blsft($e->babs(),10); # destroys copy of _e
58cde26e 570 }
091c87b1 571 elsif (!$e->is_zero()) # > 0
58cde26e 572 {
28df3e88 573 $add->blsft($e,10);
58cde26e 574 }
61f5c3f5 575 # else: both e are the same, so just leave them
576 $x->{_m}->{sign} = $x->{sign}; # fiddle with signs
58cde26e 577 $add->{sign} = $y->{sign};
61f5c3f5 578 $x->{_m} += $add; # finally do add/sub
579 $x->{sign} = $x->{_m}->{sign}; # re-adjust signs
580 $x->{_m}->{sign} = '+'; # mantissa always positiv
581 # delete trailing zeros, then round
091c87b1 582 $x->bnorm()->round($a,$p,$r,$y);
58cde26e 583 }
584
585sub bsub
586 {
0716bf9b 587 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
58cde26e 588 # subtract second arg from first, modify first
f9a08e12 589
590 # set up parameters
591 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
592 # objectify is costly, so avoid it
593 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
594 {
595 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
596 }
a0d0e21e 597
28df3e88 598 if ($y->is_zero()) # still round for not adding zero
e745a66c 599 {
28df3e88 600 return $x->round($a,$p,$r);
e745a66c 601 }
091c87b1 602
603 # $x - $y = -$x + $y
604 $y->{sign} =~ tr/+-/-+/; # does nothing for NaN
28df3e88 605 $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
091c87b1 606 $y->{sign} =~ tr/+-/-+/; # refix $y (does nothing for NaN)
e745a66c 607 $x; # already rounded by badd()
58cde26e 608 }
609
610sub binc
611 {
612 # increment arg by one
b282a552 613 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
e745a66c 614
615 if ($x->{_e}->sign() eq '-')
616 {
b282a552 617 return $x->badd($self->bone(),@r); # digits after dot
e745a66c 618 }
619
b282a552 620 if (!$x->{_e}->is_zero()) # _e == 0 for NaN, inf, -inf
e745a66c 621 {
b282a552 622 # 1e2 => 100, so after the shift below _m has a '0' as last digit
e745a66c 623 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
b282a552 624 $x->{_e}->bzero(); # normalize
625 # we know that the last digit of $x will be '1' or '9', depending on the
626 # sign
e745a66c 627 }
628 # now $x->{_e} == 0
629 if ($x->{sign} eq '+')
630 {
631 $x->{_m}->binc();
b282a552 632 return $x->bnorm()->bround(@r);
e745a66c 633 }
634 elsif ($x->{sign} eq '-')
635 {
636 $x->{_m}->bdec();
637 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
b282a552 638 return $x->bnorm()->bround(@r);
e745a66c 639 }
640 # inf, nan handling etc
b282a552 641 $x->badd($self->bone(),@r); # badd() does round
58cde26e 642 }
643
644sub bdec
645 {
646 # decrement arg by one
b282a552 647 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
e745a66c 648
649 if ($x->{_e}->sign() eq '-')
650 {
b282a552 651 return $x->badd($self->bone('-'),@r); # digits after dot
e745a66c 652 }
653
654 if (!$x->{_e}->is_zero())
655 {
656 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
657 $x->{_e}->bzero();
658 }
659 # now $x->{_e} == 0
660 my $zero = $x->is_zero();
661 # <= 0
662 if (($x->{sign} eq '-') || $zero)
663 {
664 $x->{_m}->binc();
665 $x->{sign} = '-' if $zero; # 0 => 1 => -1
666 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
b282a552 667 return $x->bnorm()->round(@r);
e745a66c 668 }
669 # > 0
670 elsif ($x->{sign} eq '+')
671 {
672 $x->{_m}->bdec();
b282a552 673 return $x->bnorm()->round(@r);
e745a66c 674 }
675 # inf, nan handling etc
b282a552 676 $x->badd($self->bone('-'),@r); # does round
58cde26e 677 }
678
990fb837 679sub DEBUG () { 0; }
680
61f5c3f5 681sub blog
682 {
990fb837 683 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
61f5c3f5 684
990fb837 685 # $base > 0, $base != 1; if $base == undef default to $base == e
686 # $x >= 0
9393ace2 687
b3abae2a 688 # we need to limit the accuracy to protect against overflow
689 my $fallback = 0;
990fb837 690 my ($scale,@params);
691 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
61f5c3f5 692
990fb837 693 # also takes care of the "error in _find_round_parameters?" case
694 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
091c87b1 695
b3abae2a 696 # no rounding at all, so must use fallback
990fb837 697 if (scalar @params == 0)
b3abae2a 698 {
699 # simulate old behaviour
990fb837 700 $params[0] = $self->div_scale(); # and round to it as accuracy
701 $params[1] = undef; # P = undef
702 $scale = $params[0]+4; # at least four more for proper round
703 $params[2] = $r; # round mode by caller or undef
b3abae2a 704 $fallback = 1; # to clear a/p afterwards
705 }
706 else
707 {
708 # the 4 below is empirical, and there might be cases where it is not
709 # enough...
990fb837 710 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
b3abae2a 711 }
61f5c3f5 712
b3abae2a 713 return $x->bzero(@params) if $x->is_one();
990fb837 714 # base not defined => base == Euler's constant e
715 if (defined $base)
716 {
091c87b1 717 # make object, since we don't feed it through objectify() to still get the
990fb837 718 # case of $base == undef
719 $base = $self->new($base) unless ref($base);
720 # $base > 0; $base != 1
721 return $x->bnan() if $base->is_zero() || $base->is_one() ||
722 $base->{sign} ne '+';
091c87b1 723 # if $x == $base, we know the result must be 1.0
990fb837 724 return $x->bone('+',@params) if $x->bcmp($base) == 0;
725 }
61f5c3f5 726
b3abae2a 727 # when user set globals, they would interfere with our calculation, so
56d9de68 728 # disable them and later re-enable them
b3abae2a 729 no strict 'refs';
730 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
731 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
732 # we also need to disable any set A or P on $x (_find_round_parameters took
733 # them already into account), since these would interfere, too
734 delete $x->{_a}; delete $x->{_p};
9393ace2 735 # need to disable $upgrade in BigInt, to avoid deep recursion
b3abae2a 736 local $Math::BigInt::upgrade = undef;
93c87d9d 737 local $Math::BigFloat::downgrade = undef;
990fb837 738
739 # upgrade $x if $x is not a BigFloat (handle BigInt input)
740 if (!$x->isa('Math::BigFloat'))
741 {
742 $x = Math::BigFloat->new($x);
743 $self = ref($x);
744 }
b282a552 745
746 my $done = 0;
747
748 # If the base is defined and an integer, try to calculate integer result
749 # first. This is very fast, and in case the real result was found, we can
750 # stop right here.
751 if (defined $base && $base->is_int() && $x->is_int())
752 {
753 my $int = $x->{_m}->copy();
754 $int->blsft($x->{_e},10) unless $x->{_e}->is_zero();
755 $int->blog($base->as_number());
756 # if ($exact)
757 if ($base->copy()->bpow($int) == $x)
758 {
759 # found result, return it
760 $x->{_m} = $int;
761 $x->{_e} = $MBI->bzero();
762 $x->bnorm();
763 $done = 1;
764 }
765 }
766
767 if ($done == 0)
9393ace2 768 {
b282a552 769 # first calculate the log to base e (using reduction by 10 (and probably 2))
770 $self->_log_10($x,$scale);
771
772 # and if a different base was requested, convert it
773 if (defined $base)
774 {
775 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
776 # not ln, but some other base (don't modify $base)
777 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
778 }
9393ace2 779 }
990fb837 780
091c87b1 781 # shortcut to not run through _find_round_parameters again
990fb837 782 if (defined $params[0])
b3abae2a 783 {
990fb837 784 $x->bround($params[0],$params[2]); # then round accordingly
b3abae2a 785 }
786 else
787 {
990fb837 788 $x->bfround($params[1],$params[2]); # then round accordingly
b3abae2a 789 }
790 if ($fallback)
791 {
792 # clear a/p after round, since user did not request it
793 $x->{_a} = undef; $x->{_p} = undef;
794 }
795 # restore globals
796 $$abr = $ab; $$pbr = $pb;
797
798 $x;
61f5c3f5 799 }
800
990fb837 801sub _log
802 {
091c87b1 803 # internal log function to calculate ln() based on Taylor series.
990fb837 804 # Modifies $x in place.
805 my ($self,$x,$scale) = @_;
806
091c87b1 807 # in case of $x == 1, result is 0
808 return $x->bzero() if $x->is_one();
809
990fb837 810 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
811
812 # u = x-1, v = x+1
813 # _ _
814 # Taylor: | u 1 u^3 1 u^5 |
815 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
816 # |_ v 3 v^3 5 v^5 _|
817
818 # This takes much more steps to calculate the result and is thus not used
819 # u = x-1
820 # _ _
821 # Taylor: | u 1 u^2 1 u^3 |
822 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
823 # |_ x 2 x^2 3 x^3 _|
824
990fb837 825 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
826
827 $v = $x->copy(); $v->binc(); # v = x+1
828 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
829 $x->bdiv($v,$scale); # first term: u/v
830 $below = $v->copy();
831 $over = $u->copy();
832 $u *= $u; $v *= $v; # u^2, v^2
833 $below->bmul($v); # u^3, v^3
834 $over->bmul($u);
835 $factor = $self->new(3); $f = $self->new(2);
836
837 my $steps = 0 if DEBUG;
838 $limit = $self->new("1E-". ($scale-1));
839 while (3 < 5)
840 {
841 # we calculate the next term, and add it to the last
842 # when the next term is below our limit, it won't affect the outcome
843 # anymore, so we stop
844
845 # calculating the next term simple from over/below will result in quite
846 # a time hog if the input has many digits, since over and below will
847 # accumulate more and more digits, and the result will also have many
848 # digits, but in the end it is rounded to $scale digits anyway. So if we
849 # round $over and $below first, we save a lot of time for the division
850 # (not with log(1.2345), but try log (123**123) to see what I mean. This
851 # can introduce a rounding error if the division result would be f.i.
852 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
091c87b1 853 # if we truncated $over and $below we might get 0.12345. Does this matter
854 # for the end result? So we give $over and $below 4 more digits to be
855 # on the safe side (unscientific error handling as usual... :+D
990fb837 856
857 $next = $over->copy->bround($scale+4)->bdiv(
858 $below->copy->bmul($factor)->bround($scale+4),
859 $scale);
860
861## old version:
862## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
863
864 last if $next->bacmp($limit) <= 0;
865
866 delete $next->{_a}; delete $next->{_p};
867 $x->badd($next);
868 #print "step $x\n ($next - $limit = ",$next - $limit,")\n";
869 # calculate things for the next term
870 $over *= $u; $below *= $v; $factor->badd($f);
871 if (DEBUG)
872 {
873 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
874 }
875 }
876 $x->bmul($f); # $x *= 2
877 print "took $steps steps\n" if DEBUG;
878 }
879
880sub _log_10
881 {
091c87b1 882 # Internal log function based on reducing input to the range of 0.1 .. 9.99
883 # and then "correcting" the result to the proper one. Modifies $x in place.
990fb837 884 my ($self,$x,$scale) = @_;
885
886 # taking blog() from numbers greater than 10 takes a *very long* time, so we
887 # break the computation down into parts based on the observation that:
888 # blog(x*y) = blog(x) + blog(y)
889 # We set $y here to multiples of 10 so that $x is below 1 (the smaller $x is
890 # the faster it get's, especially because 2*$x takes about 10 times as long,
891 # so by dividing $x by 10 we make it at least factor 100 faster...)
892
893 # The same observation is valid for numbers smaller than 0.1 (e.g. computing
894 # log(1) is fastest, and the farther away we get from 1, the longer it takes)
895 # so we also 'break' this down by multiplying $x with 10 and subtract the
896 # log(10) afterwards to get the correct result.
897
898 # calculate nr of digits before dot
899 my $dbd = $x->{_m}->length() + $x->{_e}->numify();
900
901 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
902 # infinite recursion
903
904 my $calc = 1; # do some calculation?
905
906 # disable the shortcut for 10, since we need log(10) and this would recurse
907 # infinitely deep
908 if ($x->{_e}->is_one() && $x->{_m}->is_one())
909 {
910 $dbd = 0; # disable shortcut
911 # we can use the cached value in these cases
912 if ($scale <= $LOG_10_A)
913 {
914 $x->bzero(); $x->badd($LOG_10);
915 $calc = 0; # no need to calc, but round
916 }
917 }
091c87b1 918 else
990fb837 919 {
091c87b1 920 # disable the shortcut for 2, since we maybe have it cached
921 if ($x->{_e}->is_zero() && $x->{_m}->bcmp(2) == 0)
990fb837 922 {
091c87b1 923 $dbd = 0; # disable shortcut
924 # we can use the cached value in these cases
925 if ($scale <= $LOG_2_A)
926 {
927 $x->bzero(); $x->badd($LOG_2);
928 $calc = 0; # no need to calc, but round
929 }
990fb837 930 }
931 }
932
933 # if $x = 0.1, we know the result must be 0-log(10)
091c87b1 934 if ($calc != 0 && $x->{_e}->is_one('-') && $x->{_m}->is_one())
990fb837 935 {
936 $dbd = 0; # disable shortcut
937 # we can use the cached value in these cases
938 if ($scale <= $LOG_10_A)
939 {
940 $x->bzero(); $x->bsub($LOG_10);
941 $calc = 0; # no need to calc, but round
942 }
943 }
944
091c87b1 945 return if $calc == 0; # already have the result
946
990fb837 947 # default: these correction factors are undef and thus not used
948 my $l_10; # value of ln(10) to A of $scale
949 my $l_2; # value of ln(2) to A of $scale
950
951 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
952 # so don't do this shortcut for 1 or 0
953 if (($dbd > 1) || ($dbd < 0))
954 {
955 # convert our cached value to an object if not already (avoid doing this
956 # at import() time, since not everybody needs this)
957 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
958
959 #print "x = $x, dbd = $dbd, calc = $calc\n";
960 # got more than one digit before the dot, or more than one zero after the
961 # dot, so do:
962 # log(123) == log(1.23) + log(10) * 2
963 # log(0.0123) == log(1.23) - log(10) * 2
964
965 if ($scale <= $LOG_10_A)
966 {
967 # use cached value
968 #print "using cached value for l_10\n";
969 $l_10 = $LOG_10->copy(); # copy for mul
970 }
971 else
972 {
973 # else: slower, compute it (but don't cache it, because it could be big)
974 # also disable downgrade for this code path
975 local $Math::BigFloat::downgrade = undef;
976 #print "l_10 = $l_10 (self = $self',
977 # ", ref(l_10) = ",ref($l_10)," scale $scale)\n";
978 #print "calculating value for l_10, scale $scale\n";
979 $l_10 = $self->new(10)->blog(undef,$scale); # scale+4, actually
980 }
981 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
982 # make object
983 $dbd = $self->new($dbd);
984 #print "dbd $dbd\n";
985 $l_10->bmul($dbd); # log(10) * (digits_before_dot-1)
986 #print "l_10 = $l_10\n";
987 #print "x = $x";
988 $x->{_e}->bsub($dbd); # 123 => 1.23
989 #print " => $x\n";
990 #print "calculating log($x) with scale=$scale\n";
991
992 }
993
994 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
995
996 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
997 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
998
091c87b1 999 my $half = $self->new('0.5');
1000 my $twos = 0; # default: none (0 times)
1001 my $two = $self->new(2);
1002 while ($x->bacmp($half) <= 0)
990fb837 1003 {
091c87b1 1004 $twos--; $x->bmul($two);
1005 }
1006 while ($x->bacmp($two) >= 0)
1007 {
1008 $twos++; $x->bdiv($two,$scale+4); # keep all digits
1009 }
1010 #print "$twos\n";
1011 # $twos > 0 => did mul 2, < 0 => did div 2 (never both)
1012 # calculate correction factor based on ln(2)
1013 if ($twos != 0)
1014 {
1015 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1016 if ($scale <= $LOG_2_A)
990fb837 1017 {
091c87b1 1018 # use cached value
1019 #print "using cached value for l_10\n";
1020 $l_2 = $LOG_2->copy(); # copy for mul
990fb837 1021 }
091c87b1 1022 else
990fb837 1023 {
091c87b1 1024 # else: slower, compute it (but don't cache it, because it could be big)
1025 # also disable downgrade for this code path
1026 local $Math::BigFloat::downgrade = undef;
1027 #print "calculating value for l_2, scale $scale\n";
1028 $l_2 = $two->blog(undef,$scale); # scale+4, actually
990fb837 1029 }
091c87b1 1030 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
990fb837 1031 }
1032
091c87b1 1033 $self->_log($x,$scale); # need to do the "normal" way
1034 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
1035 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
990fb837 1036 # all done, $x contains now the result
1037 }
1038
58cde26e 1039sub blcm
1040 {
ee15d750 1041 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
58cde26e 1042 # does not modify arguments, but returns new object
1043 # Lowest Common Multiplicator
58cde26e 1044
1045 my ($self,@arg) = objectify(0,@_);
1046 my $x = $self->new(shift @arg);
1047 while (@arg) { $x = _lcm($x,shift @arg); }
1048 $x;
1049 }
1050
1051sub bgcd
1052 {
ee15d750 1053 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
58cde26e 1054 # does not modify arguments, but returns new object
1055 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
58cde26e 1056
1057 my ($self,@arg) = objectify(0,@_);
1058 my $x = $self->new(shift @arg);
1059 while (@arg) { $x = _gcd($x,shift @arg); }
1060 $x;
1061 }
1062
b3abae2a 1063###############################################################################
1064# is_foo methods (is_negative, is_positive are inherited from BigInt)
1065
091c87b1 1066sub _is_zero_or_one
1067 {
1068 # internal, return true if BigInt arg is zero or one, saving the
1069 # two calls to is_zero() and is_one()
1070 my $x = $_[0];
1071
1072 $x->{sign} eq '+' && ($x->is_zero() || $x->is_one());
1073 }
1074
b3abae2a 1075sub is_int
1076 {
1077 # return true if arg (BFLOAT or num_str) is an integer
091c87b1 1078 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
b3abae2a 1079
1080 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1081 $x->{_e}->{sign} eq '+'; # 1e-1 => no integer
1082 0;
1083 }
1084
58cde26e 1085sub is_zero
1086 {
b3abae2a 1087 # return true if arg (BFLOAT or num_str) is zero
091c87b1 1088 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
574bacfe 1089
1090 return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
b3abae2a 1091 0;
58cde26e 1092 }
1093
1094sub is_one
1095 {
b3abae2a 1096 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
091c87b1 1097 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
ee15d750 1098
990fb837 1099 $sign = '+' if !defined $sign || $sign ne '-';
ee15d750 1100 return 1
1101 if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
b3abae2a 1102 0;
58cde26e 1103 }
1104
1105sub is_odd
1106 {
ee15d750 1107 # return true if arg (BFLOAT or num_str) is odd or false if even
091c87b1 1108 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
0716bf9b 1109
b3abae2a 1110 return 1 if ($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1111 ($x->{_e}->is_zero() && $x->{_m}->is_odd());
1112 0;
58cde26e 1113 }
1114
1115sub is_even
1116 {
b22b3e31 1117 # return true if arg (BINT or num_str) is even or false if odd
091c87b1 1118 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
0716bf9b 1119
1120 return 0 if $x->{sign} !~ /^[+-]$/; # NaN & +-inf aren't
b3abae2a 1121 return 1 if ($x->{_e}->{sign} eq '+' # 123.45 is never
1122 && $x->{_m}->is_even()); # but 1200 is
1123 0;
58cde26e 1124 }
1125
1126sub bmul
1127 {
1128 # multiply two numbers -- stolen from Knuth Vol 2 pg 233
1129 # (BINT or num_str, BINT or num_str) return BINT
f9a08e12 1130
1131 # set up parameters
1132 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1133 # objectify is costly, so avoid it
1134 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1135 {
1136 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1137 }
58cde26e 1138
58cde26e 1139 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1140
574bacfe 1141 # inf handling
1142 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1143 {
13a12e00 1144 return $x->bnan() if $x->is_zero() || $y->is_zero();
574bacfe 1145 # result will always be +-inf:
1146 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1147 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1148 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1149 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1150 return $x->binf('-');
1151 }
13a12e00 1152 # handle result = 0
1153 return $x->bzero() if $x->is_zero() || $y->is_zero();
8f675a64 1154
1155 return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
1156 ((!$x->isa($self)) || (!$y->isa($self)));
574bacfe 1157
58cde26e 1158 # aEb * cEd = (a*c)E(b+d)
394e6ffb 1159 $x->{_m}->bmul($y->{_m});
1160 $x->{_e}->badd($y->{_e});
58cde26e 1161 # adjust sign:
1162 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
394e6ffb 1163 return $x->bnorm()->round($a,$p,$r,$y);
58cde26e 1164 }
1165
1166sub bdiv
1167 {
1168 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
9393ace2 1169 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
f9a08e12 1170
1171 # set up parameters
1172 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1173 # objectify is costly, so avoid it
1174 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1175 {
1176 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1177 }
58cde26e 1178
13a12e00 1179 return $self->_div_inf($x,$y)
1180 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
574bacfe 1181
13a12e00 1182 # x== 0 # also: or y == 1 or y == -1
394e6ffb 1183 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
0716bf9b 1184
9393ace2 1185 # upgrade ?
1186 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
13a12e00 1187
58cde26e 1188 # we need to limit the accuracy to protect against overflow
574bacfe 1189 my $fallback = 0;
990fb837 1190 my (@params,$scale);
1191 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1192
1193 return $x if $x->is_nan(); # error in _find_round_parameters?
ee15d750 1194
1195 # no rounding at all, so must use fallback
990fb837 1196 if (scalar @params == 0)
58cde26e 1197 {
0716bf9b 1198 # simulate old behaviour
990fb837 1199 $params[0] = $self->div_scale(); # and round to it as accuracy
1200 $scale = $params[0]+4; # at least four more for proper round
1201 $params[2] = $r; # round mode by caller or undef
ee15d750 1202 $fallback = 1; # to clear a/p afterwards
1203 }
1204 else
1205 {
1206 # the 4 below is empirical, and there might be cases where it is not
1207 # enough...
990fb837 1208 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
a0d0e21e 1209 }
0716bf9b 1210 my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
58cde26e 1211 $scale = $lx if $lx > $scale;
58cde26e 1212 $scale = $ly if $ly > $scale;
0716bf9b 1213 my $diff = $ly - $lx;
1214 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
b3abae2a 1215
1216 # make copy of $x in case of list context for later reminder calculation
1217 my $rem;
1218 if (wantarray && !$y->is_one())
1219 {
1220 $rem = $x->copy();
1221 }
a0d0e21e 1222
58cde26e 1223 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
a0d0e21e 1224
58cde26e 1225 # check for / +-1 ( +/- 1E0)
394e6ffb 1226 if (!$y->is_one())
58cde26e 1227 {
394e6ffb 1228 # promote BigInts and it's subclasses (except when already a BigFloat)
1229 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1230
9393ace2 1231 # need to disable $upgrade in BigInt, to avoid deep recursion
1232 local $Math::BigInt::upgrade = undef; # should be parent class vs MBI
1233
394e6ffb 1234 # calculate the result to $scale digits and then round it
1235 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1236 $x->{_m}->blsft($scale,10);
1237 $x->{_m}->bdiv( $y->{_m} ); # a/c
1238 $x->{_e}->bsub( $y->{_e} ); # b-d
1239 $x->{_e}->bsub($scale); # correct for 10**scale
1240 $x->bnorm(); # remove trailing 0's
a0d0e21e 1241 }
a5f75d66 1242
091c87b1 1243 # shortcut to not run through _find_round_parameters again
990fb837 1244 if (defined $params[0])
ee15d750 1245 {
56d9de68 1246 $x->{_a} = undef; # clear before round
990fb837 1247 $x->bround($params[0],$params[2]); # then round accordingly
ee15d750 1248 }
1249 else
1250 {
56d9de68 1251 $x->{_p} = undef; # clear before round
990fb837 1252 $x->bfround($params[1],$params[2]); # then round accordingly
ee15d750 1253 }
574bacfe 1254 if ($fallback)
1255 {
1256 # clear a/p after round, since user did not request it
ee15d750 1257 $x->{_a} = undef; $x->{_p} = undef;
574bacfe 1258 }
0716bf9b 1259
58cde26e 1260 if (wantarray)
1261 {
394e6ffb 1262 if (!$y->is_one())
1263 {
990fb837 1264 $rem->bmod($y,@params); # copy already done
394e6ffb 1265 }
1266 else
1267 {
1268 $rem = $self->bzero();
1269 }
574bacfe 1270 if ($fallback)
1271 {
1272 # clear a/p after round, since user did not request it
ee15d750 1273 $rem->{_a} = undef; $rem->{_p} = undef;
574bacfe 1274 }
0716bf9b 1275 return ($x,$rem);
58cde26e 1276 }
9393ace2 1277 $x;
58cde26e 1278 }
a0d0e21e 1279
58cde26e 1280sub bmod
1281 {
1282 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
f9a08e12 1283
1284 # set up parameters
1285 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1286 # objectify is costly, so avoid it
1287 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1288 {
1289 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1290 }
a0d0e21e 1291
61f5c3f5 1292 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1293 {
1294 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
f9a08e12 1295 $x->{sign} = $re->{sign};
1296 $x->{_e} = $re->{_e};
1297 $x->{_m} = $re->{_m};
1298 return $x->round($a,$p,$r,$y);
61f5c3f5 1299 }
1300 return $x->bnan() if $x->is_zero() && $y->is_zero();
1301 return $x if $y->is_zero();
1302 return $x->bnan() if $x->is_nan() || $y->is_nan();
1303 return $x->bzero() if $y->is_one() || $x->is_zero();
58cde26e 1304
61f5c3f5 1305 # inf handling is missing here
1306
1307 my $cmp = $x->bacmp($y); # equal or $x < $y?
1308 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1309
1310 # only $y of the operands negative?
1311 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1312
1313 $x->{sign} = $y->{sign}; # calc sign first
1314 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1315
1316 my $ym = $y->{_m}->copy();
1317
1318 # 2e1 => 20
1319 $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
1320
1321 # if $y has digits after dot
1322 my $shifty = 0; # correct _e of $x by this
1323 if ($y->{_e}->{sign} eq '-') # has digits after dot
1324 {
1325 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1326 $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
1327 $x->blsft($shifty,10); # 123 => 1230, $y->{_m} is already 25
1328 }
1329 # $ym is now mantissa of $y based on exponent 0
b3abae2a 1330
61f5c3f5 1331 my $shiftx = 0; # correct _e of $x by this
1332 if ($x->{_e}->{sign} eq '-') # has digits after dot
1333 {
1334 # 123.4 % 20 => 1234 % 200
1335 $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
1336 $ym->blsft($shiftx,10);
1337 }
1338 # 123e1 % 20 => 1230 % 20
1339 if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
1340 {
1341 $x->{_m}->blsft($x->{_e},10);
1342 }
56b9c951 1343 $x->{_e} = $MBI->bzero() unless $x->{_e}->is_zero();
61f5c3f5 1344
1345 $x->{_e}->bsub($shiftx) if $shiftx != 0;
1346 $x->{_e}->bsub($shifty) if $shifty != 0;
1347
1348 # now mantissas are equalized, exponent of $x is adjusted, so calc result
b3abae2a 1349
61f5c3f5 1350 $x->{_m}->bmod($ym);
1351
1352 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
1353 $x->bnorm();
1354
1355 if ($neg != 0) # one of them negative => correct in place
1356 {
1357 my $r = $y - $x;
1358 $x->{_m} = $r->{_m};
1359 $x->{_e} = $r->{_e};
1360 $x->{sign} = '+' if $x->{_m}->is_zero(); # fix sign for -0
1361 $x->bnorm();
1362 }
1363
1364 $x->round($a,$p,$r,$y); # round and return
58cde26e 1365 }
1366
990fb837 1367sub broot
1368 {
1369 # calculate $y'th root of $x
3a427a11 1370
1371 # set up parameters
1372 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1373 # objectify is costly, so avoid it
1374 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1375 {
1376 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1377 }
990fb837 1378
1379 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1380 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1381 $y->{sign} !~ /^\+$/;
1382
1383 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1384
1385 # we need to limit the accuracy to protect against overflow
1386 my $fallback = 0;
1387 my (@params,$scale);
1388 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1389
1390 return $x if $x->is_nan(); # error in _find_round_parameters?
1391
1392 # no rounding at all, so must use fallback
1393 if (scalar @params == 0)
1394 {
1395 # simulate old behaviour
1396 $params[0] = $self->div_scale(); # and round to it as accuracy
1397 $scale = $params[0]+4; # at least four more for proper round
1398 $params[2] = $r; # round mode by caller or undef
1399 $fallback = 1; # to clear a/p afterwards
1400 }
1401 else
1402 {
1403 # the 4 below is empirical, and there might be cases where it is not
1404 # enough...
1405 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1406 }
1407
1408 # when user set globals, they would interfere with our calculation, so
1409 # disable them and later re-enable them
1410 no strict 'refs';
1411 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1412 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1413 # we also need to disable any set A or P on $x (_find_round_parameters took
1414 # them already into account), since these would interfere, too
1415 delete $x->{_a}; delete $x->{_p};
1416 # need to disable $upgrade in BigInt, to avoid deep recursion
1417 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1418
1419 # remember sign and make $x positive, since -4 ** (1/2) => -2
1420 my $sign = 0; $sign = 1 if $x->is_negative(); $x->babs();
1421
1422 if ($y->bcmp(2) == 0) # normal square root
1423 {
1424 $x->bsqrt($scale+4);
1425 }
1426 elsif ($y->is_one('-'))
1427 {
1428 # $x ** -1 => 1/$x
1429 my $u = $self->bone()->bdiv($x,$scale);
1430 # copy private parts over
1431 $x->{_m} = $u->{_m};
1432 $x->{_e} = $u->{_e};
1433 }
1434 else
1435 {
3a427a11 1436 # calculate the broot() as integer result first, and if it fits, return
1437 # it rightaway (but only if $x and $y are integer):
1438
1439 my $done = 0; # not yet
1440 if ($y->is_int() && $x->is_int())
1441 {
1442 my $int = $x->{_m}->copy();
1443 $int->blsft($x->{_e},10) unless $x->{_e}->is_zero();
1444 $int->broot($y->as_number());
1445 # if ($exact)
1446 if ($int->copy()->bpow($y) == $x)
1447 {
1448 # found result, return it
1449 $x->{_m} = $int;
1450 $x->{_e} = $MBI->bzero();
1451 $x->bnorm();
1452 $done = 1;
1453 }
1454 }
1455 if ($done == 0)
1456 {
1457 my $u = $self->bone()->bdiv($y,$scale+4);
1458 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
1459 $x->bpow($u,$scale+4); # el cheapo
1460 }
990fb837 1461 }
1462 $x->bneg() if $sign == 1;
1463
091c87b1 1464 # shortcut to not run through _find_round_parameters again
990fb837 1465 if (defined $params[0])
1466 {
1467 $x->bround($params[0],$params[2]); # then round accordingly
1468 }
1469 else
1470 {
1471 $x->bfround($params[1],$params[2]); # then round accordingly
1472 }
1473 if ($fallback)
1474 {
1475 # clear a/p after round, since user did not request it
1476 $x->{_a} = undef; $x->{_p} = undef;
1477 }
1478 # restore globals
1479 $$abr = $ab; $$pbr = $pb;
1480 $x;
1481 }
1482
58cde26e 1483sub bsqrt
1484 {
990fb837 1485 # calculate square root
ee15d750 1486 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 1487
990fb837 1488 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
1489 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
1490 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
58cde26e 1491
61f5c3f5 1492 # we need to limit the accuracy to protect against overflow
574bacfe 1493 my $fallback = 0;
990fb837 1494 my (@params,$scale);
1495 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1496
1497 return $x if $x->is_nan(); # error in _find_round_parameters?
61f5c3f5 1498
1499 # no rounding at all, so must use fallback
990fb837 1500 if (scalar @params == 0)
0716bf9b 1501 {
1502 # simulate old behaviour
990fb837 1503 $params[0] = $self->div_scale(); # and round to it as accuracy
1504 $scale = $params[0]+4; # at least four more for proper round
1505 $params[2] = $r; # round mode by caller or undef
ee15d750 1506 $fallback = 1; # to clear a/p afterwards
0716bf9b 1507 }
61f5c3f5 1508 else
1509 {
1510 # the 4 below is empirical, and there might be cases where it is not
1511 # enough...
990fb837 1512 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
61f5c3f5 1513 }
1514
1515 # when user set globals, they would interfere with our calculation, so
9393ace2 1516 # disable them and later re-enable them
61f5c3f5 1517 no strict 'refs';
1518 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
b3abae2a 1519 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
61f5c3f5 1520 # we also need to disable any set A or P on $x (_find_round_parameters took
1521 # them already into account), since these would interfere, too
1522 delete $x->{_a}; delete $x->{_p};
9393ace2 1523 # need to disable $upgrade in BigInt, to avoid deep recursion
1524 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
61f5c3f5 1525
394e6ffb 1526 my $xas = $x->as_number();
1527 my $gs = $xas->copy()->bsqrt(); # some guess
b3abae2a 1528
394e6ffb 1529 if (($x->{_e}->{sign} ne '-') # guess can't be accurate if there are
1530 # digits after the dot
b3abae2a 1531 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
394e6ffb 1532 {
1533 # exact result
56b9c951 1534 $x->{_m} = $gs; $x->{_e} = $MBI->bzero(); $x->bnorm();
091c87b1 1535 # shortcut to not run through _find_round_parameters again
990fb837 1536 if (defined $params[0])
61f5c3f5 1537 {
990fb837 1538 $x->bround($params[0],$params[2]); # then round accordingly
61f5c3f5 1539 }
1540 else
1541 {
990fb837 1542 $x->bfround($params[1],$params[2]); # then round accordingly
61f5c3f5 1543 }
1544 if ($fallback)
1545 {
1546 # clear a/p after round, since user did not request it
1547 $x->{_a} = undef; $x->{_p} = undef;
1548 }
9393ace2 1549 # re-enable A and P, upgrade is taken care of by "local"
b3abae2a 1550 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
61f5c3f5 1551 return $x;
394e6ffb 1552 }
2ab5f49d 1553
1554 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
1555 # of the result by multipyling the input by 100 and then divide the integer
1556 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
1557 # this will transform 123.456 (in $x) into 123456 (in $y1)
1558 my $y1 = $x->{_m}->copy();
1559 # We now make sure that $y1 has the same odd or even number of digits than
1560 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
1561 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
1562 # steps of 10. The length of $x does not count, since an even or odd number
1563 # of digits before the dot is not changed by adding an even number of digits
1564 # after the dot (the result is still odd or even digits long).
990fb837 1565 my $length = $y1->length();
2ab5f49d 1566 $y1->bmul(10) if $x->{_e}->is_odd();
1567 # now calculate how many digits the result of sqrt(y1) would have
990fb837 1568 my $digits = int($length / 2);
2ab5f49d 1569 # but we need at least $scale digits, so calculate how many are missing
1570 my $shift = $scale - $digits;
1571 # that should never happen (we take care of integer guesses above)
1572 # $shift = 0 if $shift < 0;
1573 # multiply in steps of 100, by shifting left two times the "missing" digits
1574 $y1->blsft($shift*2,10);
1575 # now take the square root and truncate to integer
1576 $y1->bsqrt();
1577 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
1578 # result, which is than later rounded to the desired scale.
990fb837 1579
1580 # calculate how many zeros $x had after the '.' (or before it, depending
1581 # on sign of $dat, the result should have half as many:
1582 my $dat = $length + $x->{_e}->numify();
1583
1584 if ($dat > 0)
1585 {
1586 # no zeros after the dot (e.g. 1.23, 0.49 etc)
1587 # preserve half as many digits before the dot than the input had
1588 # (but round this "up")
1589 $dat = int(($dat+1)/2);
1590 }
1591 else
1592 {
1593 $dat = int(($dat)/2);
1594 }
1595 $x->{_e}= $MBI->new( $dat - $y1->length() );
1596
2ab5f49d 1597 $x->{_m} = $y1;
61f5c3f5 1598
091c87b1 1599 # shortcut to not run through _find_round_parameters again
990fb837 1600 if (defined $params[0])
61f5c3f5 1601 {
990fb837 1602 $x->bround($params[0],$params[2]); # then round accordingly
61f5c3f5 1603 }
1604 else
1605 {
990fb837 1606 $x->bfround($params[1],$params[2]); # then round accordingly
61f5c3f5 1607 }
574bacfe 1608 if ($fallback)
1609 {
1610 # clear a/p after round, since user did not request it
ee15d750 1611 $x->{_a} = undef; $x->{_p} = undef;
574bacfe 1612 }
61f5c3f5 1613 # restore globals
b3abae2a 1614 $$abr = $ab; $$pbr = $pb;
574bacfe 1615 $x;
58cde26e 1616 }
1617
b3abae2a 1618sub bfac
1619 {
28df3e88 1620 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
091c87b1 1621 # compute factorial number, modifies first argument
b3abae2a 1622
b282a552 1623 # set up parameters
1624 my ($self,$x,@r) = (ref($_[0]),@_);
1625 # objectify is costly, so avoid it
1626 ($self,$x,@r) = objectify(1,@_) if !ref($x);
1627
1628 return $x if $x->{sign} eq '+inf'; # inf => inf
28df3e88 1629 return $x->bnan()
1630 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1631 ($x->{_e}->{sign} ne '+')); # digits after dot?
b3abae2a 1632
b3abae2a 1633 # use BigInt's bfac() for faster calc
b282a552 1634 if (! $x->{_e}->is_zero())
091c87b1 1635 {
b282a552 1636 $x->{_m}->blsft($x->{_e},10); # change 12e1 to 120e0
1637 $x->{_e}->bzero();
091c87b1 1638 }
091c87b1 1639 $x->{_m}->bfac(); # calculate factorial
1640 $x->bnorm()->round(@r); # norm again and round result
b3abae2a 1641 }
1642
9393ace2 1643sub _pow
1644 {
1645 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1646 my ($x,$y,$a,$p,$r) = @_;
1647 my $self = ref($x);
1648
1649 # if $y == 0.5, it is sqrt($x)
1650 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
1651
990fb837 1652 # Using:
1653 # a ** x == e ** (x * ln a)
1654
9393ace2 1655 # u = y * ln x
990fb837 1656 # _ _
1657 # Taylor: | u u^2 u^3 |
1658 # x ** y = 1 + | --- + --- + ----- + ... |
1659 # |_ 1 1*2 1*2*3 _|
9393ace2 1660
1661 # we need to limit the accuracy to protect against overflow
1662 my $fallback = 0;
990fb837 1663 my ($scale,@params);
1664 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1665
1666 return $x if $x->is_nan(); # error in _find_round_parameters?
9393ace2 1667
1668 # no rounding at all, so must use fallback
990fb837 1669 if (scalar @params == 0)
9393ace2 1670 {
1671 # simulate old behaviour
990fb837 1672 $params[0] = $self->div_scale(); # and round to it as accuracy
1673 $params[1] = undef; # disable P
1674 $scale = $params[0]+4; # at least four more for proper round
1675 $params[2] = $r; # round mode by caller or undef
9393ace2 1676 $fallback = 1; # to clear a/p afterwards
1677 }
1678 else
1679 {
1680 # the 4 below is empirical, and there might be cases where it is not
1681 # enough...
990fb837 1682 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
9393ace2 1683 }
1684
1685 # when user set globals, they would interfere with our calculation, so
56d9de68 1686 # disable them and later re-enable them
9393ace2 1687 no strict 'refs';
1688 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1689 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1690 # we also need to disable any set A or P on $x (_find_round_parameters took
1691 # them already into account), since these would interfere, too
1692 delete $x->{_a}; delete $x->{_p};
1693 # need to disable $upgrade in BigInt, to avoid deep recursion
1694 local $Math::BigInt::upgrade = undef;
1695
1696 my ($limit,$v,$u,$below,$factor,$next,$over);
1697
990fb837 1698 $u = $x->copy()->blog(undef,$scale)->bmul($y);
9393ace2 1699 $v = $self->bone(); # 1
1700 $factor = $self->new(2); # 2
1701 $x->bone(); # first term: 1
1702
1703 $below = $v->copy();
1704 $over = $u->copy();
1705
1706 $limit = $self->new("1E-". ($scale-1));
1707 #my $steps = 0;
1708 while (3 < 5)
1709 {
1710 # we calculate the next term, and add it to the last
1711 # when the next term is below our limit, it won't affect the outcome
1712 # anymore, so we stop
1713 $next = $over->copy()->bdiv($below,$scale);
990fb837 1714 last if $next->bacmp($limit) <= 0;
9393ace2 1715 $x->badd($next);
9393ace2 1716 # calculate things for the next term
1717 $over *= $u; $below *= $factor; $factor->binc();
1718 #$steps++;
1719 }
1720
091c87b1 1721 # shortcut to not run through _find_round_parameters again
990fb837 1722 if (defined $params[0])
9393ace2 1723 {
990fb837 1724 $x->bround($params[0],$params[2]); # then round accordingly
9393ace2 1725 }
1726 else
1727 {
990fb837 1728 $x->bfround($params[1],$params[2]); # then round accordingly
9393ace2 1729 }
1730 if ($fallback)
1731 {
1732 # clear a/p after round, since user did not request it
1733 $x->{_a} = undef; $x->{_p} = undef;
1734 }
1735 # restore globals
1736 $$abr = $ab; $$pbr = $pb;
1737 $x;
1738 }
1739
58cde26e 1740sub bpow
1741 {
1742 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1743 # compute power of two numbers, second arg is used as integer
1744 # modifies first argument
1745
f9a08e12 1746 # set up parameters
1747 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1748 # objectify is costly, so avoid it
1749 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1750 {
1751 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1752 }
58cde26e 1753
0716bf9b 1754 return $x if $x->{sign} =~ /^[+-]inf$/;
58cde26e 1755 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
574bacfe 1756 return $x->bone() if $y->is_zero();
58cde26e 1757 return $x if $x->is_one() || $y->is_one();
9393ace2 1758
d614cd8b 1759 return $x->_pow($y,$a,$p,$r) if !$y->is_int(); # non-integer power
9393ace2 1760
1761 my $y1 = $y->as_number(); # make bigint
394e6ffb 1762 # if ($x == -1)
1763 if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
58cde26e 1764 {
1765 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
0716bf9b 1766 return $y1->is_odd() ? $x : $x->babs(1);
288d023a 1767 }
28df3e88 1768 if ($x->is_zero())
1769 {
1770 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1771 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1772 $x->binf();
1773 }
58cde26e 1774
1775 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1776 $y1->babs();
1777 $x->{_m}->bpow($y1);
1778 $x->{_e}->bmul($y1);
1779 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1780 $x->bnorm();
1781 if ($y->{sign} eq '-')
1782 {
1783 # modify $x in place!
0716bf9b 1784 my $z = $x->copy(); $x->bzero()->binc();
58cde26e 1785 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
a0d0e21e 1786 }
28df3e88 1787 $x->round($a,$p,$r,$y);
58cde26e 1788 }
1789
1790###############################################################################
1791# rounding functions
1792
1793sub bfround
1794 {
1795 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1796 # $n == 0 means round to integer
1797 # expects and returns normalized numbers!
ee15d750 1798 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
a0d0e21e 1799
58cde26e 1800 return $x if $x->modify('bfround');
1801
ee15d750 1802 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
58cde26e 1803 return $x if !defined $scale; # no-op
1804
574bacfe 1805 # never round a 0, +-inf, NaN
61f5c3f5 1806 if ($x->is_zero())
1807 {
1808 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1809 return $x;
1810 }
1811 return $x if $x->{sign} !~ /^[+-]$/;
58cde26e 1812
ee15d750 1813 # don't round if x already has lower precision
1814 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1815
1816 $x->{_p} = $scale; # remember round in any case
1817 $x->{_a} = undef; # and clear A
58cde26e 1818 if ($scale < 0)
1819 {
58cde26e 1820 # round right from the '.'
f9a08e12 1821
1822 return $x if $x->{_e}->{sign} eq '+'; # e >= 0 => nothing to round
1823
58cde26e 1824 $scale = -$scale; # positive for simplicity
1825 my $len = $x->{_m}->length(); # length of mantissa
f9a08e12 1826
1827 # the following poses a restriction on _e, but if _e is bigger than a
1828 # scalar, you got other problems (memory etc) anyway
1829 my $dad = -($x->{_e}->numify()); # digits after dot
58cde26e 1830 my $zad = 0; # zeros after dot
f9a08e12 1831 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
1832
ee15d750 1833 #print "scale $scale dad $dad zad $zad len $len\n";
58cde26e 1834 # number bsstr len zad dad
1835 # 0.123 123e-3 3 0 3
1836 # 0.0123 123e-4 3 1 4
1837 # 0.001 1e-3 1 2 3
1838 # 1.23 123e-2 3 0 2
1839 # 1.2345 12345e-4 5 0 4
1840
1841 # do not round after/right of the $dad
1842 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
1843
ee15d750 1844 # round to zero if rounding inside the $zad, but not for last zero like:
1845 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1846 return $x->bzero() if $scale < $zad;
1847 if ($scale == $zad) # for 0.006, scale -3 and trunc
58cde26e 1848 {
b3abae2a 1849 $scale = -$len;
58cde26e 1850 }
1851 else
1852 {
1853 # adjust round-point to be inside mantissa
1854 if ($zad != 0)
1855 {
1856 $scale = $scale-$zad;
1857 }
1858 else
1859 {
1860 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
1861 $scale = $dbd+$scale;
1862 }
1863 }
a0d0e21e 1864 }
58cde26e 1865 else
1866 {
f9a08e12 1867 # round left from the '.'
1868
58cde26e 1869 # 123 => 100 means length(123) = 3 - $scale (2) => 1
a5f75d66 1870
b3abae2a 1871 my $dbt = $x->{_m}->length();
1872 # digits before dot
f9a08e12 1873 my $dbd = $dbt + $x->{_e}->numify();
b3abae2a 1874 # should be the same, so treat it as this
1875 $scale = 1 if $scale == 0;
1876 # shortcut if already integer
1877 return $x if $scale == 1 && $dbt <= $dbd;
1878 # maximum digits before dot
1879 ++$dbd;
1880
1881 if ($scale > $dbd)
1882 {
1883 # not enough digits before dot, so round to zero
1884 return $x->bzero;
1885 }
1886 elsif ( $scale == $dbd )
1887 {
1888 # maximum
1889 $scale = -$dbt;
1890 }
58cde26e 1891 else
b3abae2a 1892 {
1893 $scale = $dbd - $scale;
1894 }
a0d0e21e 1895 }
574bacfe 1896 # pass sign to bround for rounding modes '+inf' and '-inf'
58cde26e 1897 $x->{_m}->{sign} = $x->{sign};
1898 $x->{_m}->bround($scale,$mode);
1899 $x->{_m}->{sign} = '+'; # fix sign back
1900 $x->bnorm();
1901 }
1902
1903sub bround
1904 {
1905 # accuracy: preserve $N digits, and overwrite the rest with 0's
ee15d750 1906 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1907
990fb837 1908 if (($_[0] || 0) < 0)
1909 {
1910 require Carp; Carp::croak ('bround() needs positive accuracy');
1911 }
58cde26e 1912
ee15d750 1913 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1914 return $x if !defined $scale; # no-op
61f5c3f5 1915
58cde26e 1916 return $x if $x->modify('bround');
61f5c3f5 1917
ee15d750 1918 # scale is now either $x->{_a}, $accuracy, or the user parameter
1919 # test whether $x already has lower accuracy, do nothing in this case
1920 # but do round if the accuracy is the same, since a math operation might
1921 # want to round a number with A=5 to 5 digits afterwards again
1922 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
58cde26e 1923
61f5c3f5 1924 # scale < 0 makes no sense
1925 # never round a +-inf, NaN
1926 return $x if ($scale < 0) || $x->{sign} !~ /^[+-]$/;
58cde26e 1927
61f5c3f5 1928 # 1: $scale == 0 => keep all digits
1929 # 2: never round a 0
1930 # 3: if we should keep more digits than the mantissa has, do nothing
1931 if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1932 {
1933 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1934 return $x;
1935 }
f216259d 1936
58cde26e 1937 # pass sign to bround for '+inf' and '-inf' rounding modes
1938 $x->{_m}->{sign} = $x->{sign};
1939 $x->{_m}->bround($scale,$mode); # round mantissa
1940 $x->{_m}->{sign} = '+'; # fix sign back
61f5c3f5 1941 # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
ee15d750 1942 $x->{_a} = $scale; # remember rounding
1943 $x->{_p} = undef; # and clear P
574bacfe 1944 $x->bnorm(); # del trailing zeros gen. by bround()
58cde26e 1945 }
1946
1947sub bfloor
1948 {
1949 # return integer less or equal then $x
ee15d750 1950 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 1951
1952 return $x if $x->modify('bfloor');
1953
1954 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1955
1956 # if $x has digits after dot
1957 if ($x->{_e}->{sign} eq '-')
1958 {
28df3e88 1959 $x->{_e}->{sign} = '+'; # negate e
1960 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1961 $x->{_e}->bzero(); # trunc/norm
1962 $x->{_m}->binc() if $x->{sign} eq '-'; # decrement if negative
f216259d 1963 }
61f5c3f5 1964 $x->round($a,$p,$r);
58cde26e 1965 }
288d023a 1966
58cde26e 1967sub bceil
1968 {
1969 # return integer greater or equal then $x
ee15d750 1970 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 1971
1972 return $x if $x->modify('bceil');
1973 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1974
1975 # if $x has digits after dot
1976 if ($x->{_e}->{sign} eq '-')
1977 {
28df3e88 1978 #$x->{_m}->brsft(-$x->{_e},10);
1979 #$x->{_e}->bzero();
1980 #$x++ if $x->{sign} eq '+';
1981
1982 $x->{_e}->{sign} = '+'; # negate e
1983 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1984 $x->{_e}->bzero(); # trunc/norm
1985 $x->{_m}->binc() if $x->{sign} eq '+'; # decrement if negative
a0d0e21e 1986 }
61f5c3f5 1987 $x->round($a,$p,$r);
58cde26e 1988 }
1989
394e6ffb 1990sub brsft
1991 {
f9a08e12 1992 # shift right by $y (divide by power of $n)
1993
1994 # set up parameters
1995 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1996 # objectify is costly, so avoid it
1997 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1998 {
1999 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2000 }
394e6ffb 2001
2002 return $x if $x->modify('brsft');
2003 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2004
f9a08e12 2005 $n = 2 if !defined $n; $n = $self->new($n);
2006 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
394e6ffb 2007 }
2008
2009sub blsft
2010 {
f9a08e12 2011 # shift left by $y (multiply by power of $n)
2012
2013 # set up parameters
2014 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2015 # objectify is costly, so avoid it
2016 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2017 {
2018 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2019 }
394e6ffb 2020
f9a08e12 2021 return $x if $x->modify('blsft');
394e6ffb 2022 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2023
f9a08e12 2024 $n = 2 if !defined $n; $n = $self->new($n);
2025 $x->bmul($n->bpow($y),$a,$p,$r,$y);
394e6ffb 2026 }
2027
58cde26e 2028###############################################################################
a5f75d66 2029
58cde26e 2030sub DESTROY
2031 {
b282a552 2032 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
58cde26e 2033 }
2034
2035sub AUTOLOAD
2036 {
b3abae2a 2037 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
2038 # or falling back to MBI::bxxx()
58cde26e 2039 my $name = $AUTOLOAD;
2040
2041 $name =~ s/.*:://; # split package
ee15d750 2042 no strict 'refs';
990fb837 2043 $class->import() if $IMPORT == 0;
ee15d750 2044 if (!method_alias($name))
58cde26e 2045 {
ee15d750 2046 if (!defined $name)
2047 {
2048 # delayed load of Carp and avoid recursion
2049 require Carp;
2050 Carp::croak ("Can't call a method without name");
2051 }
ee15d750 2052 if (!method_hand_up($name))
2053 {
2054 # delayed load of Carp and avoid recursion
2055 require Carp;
2056 Carp::croak ("Can't call $class\-\>$name, not a valid method");
2057 }
2058 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
2059 $name =~ s/^f/b/;
56b9c951 2060 return &{"$MBI"."::$name"}(@_);
a0d0e21e 2061 }
58cde26e 2062 my $bname = $name; $bname =~ s/^f/b/;
b3abae2a 2063 *{$class."::$name"} = \&$bname;
58cde26e 2064 &$bname; # uses @_
2065 }
2066
2067sub exponent
2068 {
2069 # return a copy of the exponent
ee15d750 2070 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2071
ee15d750 2072 if ($x->{sign} !~ /^[+-]$/)
2073 {
2074 my $s = $x->{sign}; $s =~ s/^[+-]//;
2075 return $self->new($s); # -inf, +inf => +inf
2076 }
2077 return $x->{_e}->copy();
58cde26e 2078 }
2079
2080sub mantissa
2081 {
2082 # return a copy of the mantissa
ee15d750 2083 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2084
ee15d750 2085 if ($x->{sign} !~ /^[+-]$/)
2086 {
2087 my $s = $x->{sign}; $s =~ s/^[+]//;
2088 return $self->new($s); # -inf, +inf => +inf
2089 }
2090 my $m = $x->{_m}->copy(); # faster than going via bstr()
2091 $m->bneg() if $x->{sign} eq '-';
58cde26e 2092
61f5c3f5 2093 $m;
58cde26e 2094 }
2095
2096sub parts
2097 {
2098 # return a copy of both the exponent and the mantissa
ee15d750 2099 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2100
ee15d750 2101 if ($x->{sign} !~ /^[+-]$/)
2102 {
2103 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2104 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2105 }
2106 my $m = $x->{_m}->copy(); # faster than going via bstr()
2107 $m->bneg() if $x->{sign} eq '-';
2108 return ($m,$x->{_e}->copy());
58cde26e 2109 }
2110
2111##############################################################################
2112# private stuff (internal use only)
2113
58cde26e 2114sub import
2115 {
2116 my $self = shift;
8f675a64 2117 my $l = scalar @_;
2118 my $lib = ''; my @a;
990fb837 2119 $IMPORT=1;
8f675a64 2120 for ( my $i = 0; $i < $l ; $i++)
58cde26e 2121 {
2122 if ( $_[$i] eq ':constant' )
2123 {
091c87b1 2124 # This causes overlord er load to step in. 'binary' and 'integer'
2125 # are handled by BigInt.
58cde26e 2126 overload::constant float => sub { $self->new(shift); };
b3abae2a 2127 }
2128 elsif ($_[$i] eq 'upgrade')
2129 {
2130 # this causes upgrading
28df3e88 2131 $upgrade = $_[$i+1]; # or undef to disable
8f675a64 2132 $i++;
28df3e88 2133 }
2134 elsif ($_[$i] eq 'downgrade')
2135 {
2136 # this causes downgrading
2137 $downgrade = $_[$i+1]; # or undef to disable
8f675a64 2138 $i++;
58cde26e 2139 }
56b9c951 2140 elsif ($_[$i] eq 'lib')
2141 {
990fb837 2142 # alternative library
56b9c951 2143 $lib = $_[$i+1] || ''; # default Calc
8f675a64 2144 $i++;
56b9c951 2145 }
2146 elsif ($_[$i] eq 'with')
2147 {
990fb837 2148 # alternative class for our private parts()
56b9c951 2149 $MBI = $_[$i+1] || 'Math::BigInt'; # default Math::BigInt
8f675a64 2150 $i++;
2151 }
2152 else
2153 {
2154 push @a, $_[$i];
56b9c951 2155 }
58cde26e 2156 }
8f675a64 2157
56b9c951 2158 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2159 my $mbilib = eval { Math::BigInt->config()->{lib} };
8f675a64 2160 if ((defined $mbilib) && ($MBI eq 'Math::BigInt'))
2161 {
2162 # MBI already loaded
2163 $MBI->import('lib',"$lib,$mbilib", 'objectify');
2164 }
2165 else
2166 {
2167 # MBI not loaded, or with ne "Math::BigInt"
2168 $lib .= ",$mbilib" if defined $mbilib;
07d34614 2169 $lib =~ s/^,//; # don't leave empty
990fb837 2170 # replacement library can handle lib statement, but also could ignore it
8f675a64 2171 if ($] < 5.006)
2172 {
2173 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2174 # used in the same script, or eval inside import().
2175 my @parts = split /::/, $MBI; # Math::BigInt => Math BigInt
2176 my $file = pop @parts; $file .= '.pm'; # BigInt => BigInt.pm
07d34614 2177 require File::Spec;
8f675a64 2178 $file = File::Spec->catfile (@parts, $file);
07d34614 2179 eval { require "$file"; };
2180 $MBI->import( lib => $lib, 'objectify' );
8f675a64 2181 }
2182 else
2183 {
2184 my $rc = "use $MBI lib => '$lib', 'objectify';";
2185 eval $rc;
2186 }
2187 }
990fb837 2188 if ($@)
2189 {
2190 require Carp; Carp::croak ("Couldn't load $MBI: $! $@");
2191 }
56b9c951 2192
58cde26e 2193 # any non :constant stuff is handled by our parent, Exporter
2194 # even if @_ is empty, to give it a chance
b3abae2a 2195 $self->SUPER::import(@a); # for subclasses
2196 $self->export_to_level(1,$self,@a); # need this, too
58cde26e 2197 }
2198
2199sub bnorm
2200 {
2201 # adjust m and e so that m is smallest possible
2202 # round number according to accuracy and precision settings
ee15d750 2203 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2204
0716bf9b 2205 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
58cde26e 2206
b282a552 2207 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
2208 if ($zeros != 0)
2209 {
2210 my $z = $MBI->new($zeros,undef,undef);
2211 $x->{_m}->brsft($z,10); $x->{_e}->badd($z);
2212 }
2213 else
2214 {
2215 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
2216 # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
28df3e88 2217 $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
b282a552 2218 }
2219
ee15d750 2220 # this is to prevent automatically rounding when MBI's globals are set
0716bf9b 2221 $x->{_m}->{_f} = MB_NEVER_ROUND;
2222 $x->{_e}->{_f} = MB_NEVER_ROUND;
ee15d750 2223 # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
2224 $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
2225 $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
61f5c3f5 2226 $x; # MBI bnorm is no-op, so dont call it
2227 }
58cde26e 2228
2229##############################################################################
56d9de68 2230
2231sub as_hex
2232 {
2233 # return number as hexadecimal string (only for integers defined)
2234 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2235
2236 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2237 return '0x0' if $x->is_zero();
2238
990fb837 2239 return $nan if $x->{_e}->{sign} ne '+'; # how to do 1e-1 in hex!?
56d9de68 2240
2241 my $z = $x->{_m}->copy();
2242 if (!$x->{_e}->is_zero()) # > 0
2243 {
2244 $z->blsft($x->{_e},10);
2245 }
2246 $z->{sign} = $x->{sign};
2247 $z->as_hex();
2248 }
2249
2250sub as_bin
2251 {
2252 # return number as binary digit string (only for integers defined)
2253 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2254
2255 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2256 return '0b0' if $x->is_zero();
2257
990fb837 2258 return $nan if $x->{_e}->{sign} ne '+'; # how to do 1e-1 in hex!?
56d9de68 2259
2260 my $z = $x->{_m}->copy();
2261 if (!$x->{_e}->is_zero()) # > 0
2262 {
2263 $z->blsft($x->{_e},10);
2264 }
2265 $z->{sign} = $x->{sign};
2266 $z->as_bin();
2267 }
58cde26e 2268
2269sub as_number
2270 {
394e6ffb 2271 # return copy as a bigint representation of this BigFloat number
2272 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2273
28df3e88 2274 my $z = $x->{_m}->copy();
2275 if ($x->{_e}->{sign} eq '-') # < 0
58cde26e 2276 {
28df3e88 2277 $x->{_e}->{sign} = '+'; # flip
2278 $z->brsft($x->{_e},10);
2279 $x->{_e}->{sign} = '-'; # flip back
0716bf9b 2280 }
28df3e88 2281 elsif (!$x->{_e}->is_zero()) # > 0
0716bf9b 2282 {
2283 $z->blsft($x->{_e},10);
58cde26e 2284 }
58cde26e 2285 $z->{sign} = $x->{sign};
61f5c3f5 2286 $z;
58cde26e 2287 }
2288
2289sub length
2290 {
ee15d750 2291 my $x = shift;
2292 my $class = ref($x) || $x;
2293 $x = $class->new(shift) unless ref($x);
58cde26e 2294
ee15d750 2295 return 1 if $x->{_m}->is_zero();
58cde26e 2296 my $len = $x->{_m}->length();
2297 $len += $x->{_e} if $x->{_e}->sign() eq '+';
2298 if (wantarray())
2299 {
56b9c951 2300 my $t = $MBI->bzero();
58cde26e 2301 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
2302 return ($len,$t);
2303 }
61f5c3f5 2304 $len;
58cde26e 2305 }
a0d0e21e 2306
23071;
a5f75d66 2308__END__
2309
2310=head1 NAME
2311
58cde26e 2312Math::BigFloat - Arbitrary size floating point math package
a5f75d66 2313
2314=head1 SYNOPSIS
2315
a2008d6d 2316 use Math::BigFloat;
58cde26e 2317
b3abae2a 2318 # Number creation
2319 $x = Math::BigFloat->new($str); # defaults to 0
2320 $nan = Math::BigFloat->bnan(); # create a NotANumber
2321 $zero = Math::BigFloat->bzero(); # create a +0
2322 $inf = Math::BigFloat->binf(); # create a +inf
2323 $inf = Math::BigFloat->binf('-'); # create a -inf
2324 $one = Math::BigFloat->bone(); # create a +1
2325 $one = Math::BigFloat->bone('-'); # create a -1
58cde26e 2326
2327 # Testing
b3abae2a 2328 $x->is_zero(); # true if arg is +0
2329 $x->is_nan(); # true if arg is NaN
0716bf9b 2330 $x->is_one(); # true if arg is +1
2331 $x->is_one('-'); # true if arg is -1
2332 $x->is_odd(); # true if odd, false for even
2333 $x->is_even(); # true if even, false for odd
2334 $x->is_positive(); # true if >= 0
2335 $x->is_negative(); # true if < 0
b3abae2a 2336 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
2337
58cde26e 2338 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
2339 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
2340 $x->sign(); # return the sign, either +,- or NaN
b3abae2a 2341 $x->digit($n); # return the nth digit, counting from right
2342 $x->digit(-$n); # return the nth digit, counting from left
58cde26e 2343
990fb837 2344 # The following all modify their first argument. If you want to preserve
2345 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2346 # neccessary when mixing $a = $b assigments with non-overloaded math.
2347
58cde26e 2348 # set
2349 $x->bzero(); # set $i to 0
2350 $x->bnan(); # set $i to NaN
b3abae2a 2351 $x->bone(); # set $x to +1
2352 $x->bone('-'); # set $x to -1
2353 $x->binf(); # set $x to inf
2354 $x->binf('-'); # set $x to -inf
58cde26e 2355
2356 $x->bneg(); # negation
2357 $x->babs(); # absolute value
2358 $x->bnorm(); # normalize (no-op)
2359 $x->bnot(); # two's complement (bit wise not)
2360 $x->binc(); # increment x by 1
2361 $x->bdec(); # decrement x by 1
2362
2363 $x->badd($y); # addition (add $y to $x)
2364 $x->bsub($y); # subtraction (subtract $y from $x)
2365 $x->bmul($y); # multiplication (multiply $x by $y)
990fb837 2366 $x->bdiv($y); # divide, set $x to quotient
58cde26e 2367 # return (quo,rem) or quo if scalar
2368
990fb837 2369 $x->bmod($y); # modulus ($x % $y)
2370 $x->bpow($y); # power of arguments ($x ** $y)
58cde26e 2371 $x->blsft($y); # left shift
2372 $x->brsft($y); # right shift
2373 # return (quo,rem) or quo if scalar
2374
990fb837 2375 $x->blog(); # logarithm of $x to base e (Euler's number)
2376 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
61f5c3f5 2377
58cde26e 2378 $x->band($y); # bit-wise and
2379 $x->bior($y); # bit-wise inclusive or
2380 $x->bxor($y); # bit-wise exclusive or
2381 $x->bnot(); # bit-wise not (two's complement)
b3abae2a 2382
2383 $x->bsqrt(); # calculate square-root
990fb837 2384 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
b3abae2a 2385 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
2386
990fb837 2387 $x->bround($N); # accuracy: preserve $N digits
58cde26e 2388 $x->bfround($N); # precision: round to the $Nth digit
2389
990fb837 2390 $x->bfloor(); # return integer less or equal than $x
2391 $x->bceil(); # return integer greater or equal than $x
2392
58cde26e 2393 # The following do not modify their arguments:
990fb837 2394
58cde26e 2395 bgcd(@values); # greatest common divisor
2396 blcm(@values); # lowest common multiplicator
2397
2398 $x->bstr(); # return string
2399 $x->bsstr(); # return string in scientific notation
b3abae2a 2400
58cde26e 2401 $x->exponent(); # return exponent as BigInt
2402 $x->mantissa(); # return mantissa as BigInt
2403 $x->parts(); # return (mantissa,exponent) as BigInt
2404
2405 $x->length(); # number of digits (w/o sign and '.')
2406 ($l,$f) = $x->length(); # number of digits, and length of fraction
a5f75d66 2407
f9a08e12 2408 $x->precision(); # return P of $x (or global, if P of $x undef)
2409 $x->precision($n); # set P of $x to $n
2410 $x->accuracy(); # return A of $x (or global, if A of $x undef)
723d369b 2411 $x->accuracy($n); # set A $x to $n
f9a08e12 2412
990fb837 2413 # these get/set the appropriate global value for all BigFloat objects
2414 Math::BigFloat->precision(); # Precision
2415 Math::BigFloat->accuracy(); # Accuracy
2416 Math::BigFloat->round_mode(); # rounding mode
f9a08e12 2417
a5f75d66 2418=head1 DESCRIPTION
2419
58cde26e 2420All operators (inlcuding basic math operations) are overloaded if you
2421declare your big floating point numbers as
a5f75d66 2422
58cde26e 2423 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2424
2425Operations with overloaded operators preserve the arguments, which is
2426exactly what you expect.
2427
2428=head2 Canonical notation
2429
2430Input to these routines are either BigFloat objects, or strings of the
2431following four forms:
a5f75d66 2432
2433=over 2
2434
58cde26e 2435=item *
2436
2437C</^[+-]\d+$/>
a5f75d66 2438
58cde26e 2439=item *
a5f75d66 2440
58cde26e 2441C</^[+-]\d+\.\d*$/>
a5f75d66 2442
58cde26e 2443=item *
a5f75d66 2444
58cde26e 2445C</^[+-]\d+E[+-]?\d+$/>
a5f75d66 2446
58cde26e 2447=item *
a5f75d66 2448
58cde26e 2449C</^[+-]\d*\.\d+E[+-]?\d+$/>
5d7098d5 2450
58cde26e 2451=back
2452
2453all with optional leading and trailing zeros and/or spaces. Additonally,
2454numbers are allowed to have an underscore between any two digits.
2455
2456Empty strings as well as other illegal numbers results in 'NaN'.
2457
2458bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2459are always stored in normalized form. On a string, it creates a BigFloat
2460object.
2461
2462=head2 Output
2463
2464Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2465
2466The string output will always have leading and trailing zeros stripped and drop
2467a plus sign. C<bstr()> will give you always the form with a decimal point,
990fb837 2468while C<bsstr()> (s for scientific) gives you the scientific notation.
58cde26e 2469
2470 Input bstr() bsstr()
2471 '-0' '0' '0E1'
2472 ' -123 123 123' '-123123123' '-123123123E0'
2473 '00.0123' '0.0123' '123E-4'
2474 '123.45E-2' '1.2345' '12345E-4'
2475 '10E+3' '10000' '1E4'
2476
2477Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2478C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2479return either undef, <0, 0 or >0 and are suited for sort.
2480
990fb837 2481Actual math is done by using the class defined with C<with => Class;> (which
2482defaults to BigInts) to represent the mantissa and exponent.
2483
58cde26e 2484The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2485represent the result when input arguments are not numbers, as well as
2486the result of dividing by zero.
2487
2488=head2 C<mantissa()>, C<exponent()> and C<parts()>
2489
2490C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2491as BigInts such that:
2492
2493 $m = $x->mantissa();
2494 $e = $x->exponent();
2495 $y = $m * ( 10 ** $e );
2496 print "ok\n" if $x == $y;
2497
2498C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2499
2500A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2501
2502Currently the mantissa is reduced as much as possible, favouring higher
2503exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2504This might change in the future, so do not depend on it.
2505
2506=head2 Accuracy vs. Precision
2507
2508See also: L<Rounding|Rounding>.
2509
027dc388 2510Math::BigFloat supports both precision and accuracy. For a full documentation,
2511examples and tips on these topics please see the large section in
2512L<Math::BigInt>.
5d7098d5 2513
58cde26e 2514Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2515a operation consumes all resources, each operation produces no more than
990fb837 2516the requested number of digits.
2517
2518Please refer to BigInt's documentation for the precedence rules of which
2519accuracy/precision setting will be used.
2520
2521If there is no gloabl precision set, B<and> the operation inquestion was not
2522called with a requested precision or accuracy, B<and> the input $x has no
2523accuracy or precision set, then a fallback parameter will be used. For
2524historical reasons, it is called C<div_scale> and can be accessed via:
2525
2526 $d = Math::BigFloat->div_scale(); # query
2527 Math::BigFloat->div_scale($n); # set to $n digits
2528
2529The default value is 40 digits.
58cde26e 2530
2531In case the result of one operation has more precision than specified,
2532it is rounded. The rounding mode taken is either the default mode, or the one
2533supplied to the operation after the I<scale>:
2534
2535 $x = Math::BigFloat->new(2);
990fb837 2536 Math::BigFloat->precision(5); # 5 digits max
58cde26e 2537 $y = $x->copy()->bdiv(3); # will give 0.66666
2538 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2539 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
990fb837 2540 Math::BigFloat->round_mode('zero');
58cde26e 2541 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2542
2543=head2 Rounding
2544
2545=over 2
2546
5dc6f178 2547=item ffround ( +$scale )
58cde26e 2548
0716bf9b 2549Rounds to the $scale'th place left from the '.', counting from the dot.
2550The first digit is numbered 1.
58cde26e 2551
5dc6f178 2552=item ffround ( -$scale )
58cde26e 2553
0716bf9b 2554Rounds to the $scale'th place right from the '.', counting from the dot.
58cde26e 2555
5dc6f178 2556=item ffround ( 0 )
2557
0716bf9b 2558Rounds to an integer.
5dc6f178 2559
2560=item fround ( +$scale )
2561
0716bf9b 2562Preserves accuracy to $scale digits from the left (aka significant digits)
2563and pads the rest with zeros. If the number is between 1 and -1, the
2564significant digits count from the first non-zero after the '.'
5dc6f178 2565
2566=item fround ( -$scale ) and fround ( 0 )
2567
990fb837 2568These are effectively no-ops.
5d7098d5 2569
a5f75d66 2570=back
2571
0716bf9b 2572All rounding functions take as a second parameter a rounding mode from one of
2573the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
58cde26e 2574
2575The default rounding mode is 'even'. By using
990fb837 2576C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
ee15d750 2577mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
0716bf9b 2578no longer supported.
b22b3e31 2579The second parameter to the round functions then overrides the default
0716bf9b 2580temporarily.
58cde26e 2581
990fb837 2582The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
58cde26e 2583'trunc' as rounding mode to make it equivalent to:
2584
2585 $x = 2.5;
2586 $y = int($x) + 2;
2587
2588You can override this by passing the desired rounding mode as parameter to
2589C<as_number()>:
2590
2591 $x = Math::BigFloat->new(2.5);
2592 $y = $x->as_number('odd'); # $y = 3
2593
2594=head1 EXAMPLES
2595
58cde26e 2596 # not ready yet
58cde26e 2597
2598=head1 Autocreating constants
2599
2600After C<use Math::BigFloat ':constant'> all the floating point constants
2601in the given scope are converted to C<Math::BigFloat>. This conversion
2602happens at compile time.
2603
2604In particular
2605
2606 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2607
56b9c951 2608prints the value of C<2E-100>. Note that without conversion of
58cde26e 2609constants the expression 2E-100 will be calculated as normal floating point
2610number.
2611
56b9c951 2612Please note that ':constant' does not affect integer constants, nor binary
2613nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2614work.
2615
2616=head2 Math library
2617
2618Math with the numbers is done (by default) by a module called
2619Math::BigInt::Calc. This is equivalent to saying:
2620
2621 use Math::BigFloat lib => 'Calc';
2622
2623You can change this by using:
2624
2625 use Math::BigFloat lib => 'BitVect';
2626
2627The following would first try to find Math::BigInt::Foo, then
2628Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2629
2630 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2631
2632Calc.pm uses as internal format an array of elements of some decimal base
2633(usually 1e7, but this might be differen for some systems) with the least
2634significant digit first, while BitVect.pm uses a bit vector of base 2, most
2635significant bit first. Other modules might use even different means of
2636representing the numbers. See the respective module documentation for further
2637details.
2638
2639Please note that Math::BigFloat does B<not> use the denoted library itself,
2640but it merely passes the lib argument to Math::BigInt. So, instead of the need
2641to do:
2642
2643 use Math::BigInt lib => 'GMP';
2644 use Math::BigFloat;
2645
2646you can roll it all into one line:
2647
2648 use Math::BigFloat lib => 'GMP';
2649
990fb837 2650It is also possible to just require Math::BigFloat:
2651
2652 require Math::BigFloat;
2653
2654This will load the neccessary things (like BigInt) when they are needed, and
2655automatically.
2656
2657Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
2658you ever wanted to know about loading a different library.
56b9c951 2659
2660=head2 Using Math::BigInt::Lite
2661
2662It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2663
2664 # 1
2665 use Math::BigFloat with => 'Math::BigInt::Lite';
2666
2667There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2668can combine these if you want. For instance, you may want to use
2669Math::BigInt objects in your main script, too.
2670
2671 # 2
2672 use Math::BigInt;
2673 use Math::BigFloat with => 'Math::BigInt::Lite';
2674
2675Of course, you can combine this with the C<lib> parameter.
2676
2677 # 3
2678 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2679
990fb837 2680There is no need for a "use Math::BigInt;" statement, even if you want to
2681use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
2682always loads it. But if you add it, add it B<before>:
56b9c951 2683
2684 # 4
2685 use Math::BigInt;
2686 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2687
2688Notice that the module with the last C<lib> will "win" and thus
2689it's lib will be used if the lib is available:
2690
2691 # 5
2692 use Math::BigInt lib => 'Bar,Baz';
2693 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2694
2695That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2696words, Math::BigFloat will try to retain previously loaded libs when you
990fb837 2697don't specify it onem but if you specify one, it will try to load them.
56b9c951 2698
2699Actually, the lib loading order would be "Bar,Baz,Calc", and then
2700"Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
990fb837 2701same as trying the latter load alone, except for the fact that one of Bar or
2702Baz might be loaded needlessly in an intermidiate step (and thus hang around
2703and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
2704will still be tried to be loaded, but this is not as time/memory consuming as
2705actually loading one of them. Still, this type of usage is not recommended due
2706to these issues.
56b9c951 2707
990fb837 2708The old way (loading the lib only in BigInt) still works though:
56b9c951 2709
2710 # 6
2711 use Math::BigInt lib => 'Bar,Baz';
2712 use Math::BigFloat;
2713
990fb837 2714You can even load Math::BigInt afterwards:
56b9c951 2715
990fb837 2716 # 7
2717 use Math::BigFloat;
2718 use Math::BigInt lib => 'Bar,Baz';
a5f75d66 2719
990fb837 2720But this has the same problems like #5, it will first load Calc
2721(Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
2722Baz, depending on which of them works and is usable/loadable. Since this
2723loads Calc unnecc., it is not recommended.
58cde26e 2724
990fb837 2725Since it also possible to just require Math::BigFloat, this poses the question
2726about what libary this will use:
58cde26e 2727
990fb837 2728 require Math::BigFloat;
2729 my $x = Math::BigFloat->new(123); $x += 123;
58cde26e 2730
990fb837 2731It will use Calc. Please note that the call to import() is still done, but
2732only when you use for the first time some Math::BigFloat math (it is triggered
2733via any constructor, so the first time you create a Math::BigFloat, the load
2734will happen in the background). This means:
58cde26e 2735
990fb837 2736 require Math::BigFloat;
2737 Math::BigFloat->import ( lib => 'Foo,Bar' );
58cde26e 2738
990fb837 2739would be the same as:
58cde26e 2740
990fb837 2741 use Math::BigFloat lib => 'Foo, Bar';
2742
2743But don't try to be clever to insert some operations in between:
2744
2745 require Math::BigFloat;
2746 my $x = Math::BigFloat->bone() + 4; # load BigInt and Calc
2747 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
2748 $x = Math::BigFloat->bone()+4; # now use Pari
2749
2750While this works, it loads Calc needlessly. But maybe you just wanted that?
2751
2752B<Examples #3 is highly recommended> for daily usage.
2753
2754=head1 BUGS
2755
2756Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
58cde26e 2757
990fb837 2758=head1 CAVEATS
58cde26e 2759
2760=over 1
2761
2762=item stringify, bstr()
2763
2764Both stringify and bstr() now drop the leading '+'. The old code would return
2765'+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2766reasoning and details.
2767
2768=item bdiv
2769
2770The following will probably not do what you expect:
2771
2772 print $c->bdiv(123.456),"\n";
2773
2774It prints both quotient and reminder since print works in list context. Also,
2775bdiv() will modify $c, so be carefull. You probably want to use
2776
2777 print $c / 123.456,"\n";
2778 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
2779
2780instead.
2781
2782=item Modifying and =
2783
2784Beware of:
2785
2786 $x = Math::BigFloat->new(5);
2787 $y = $x;
2788
2789It will not do what you think, e.g. making a copy of $x. Instead it just makes
2790a second reference to the B<same> object and stores it in $y. Thus anything
990fb837 2791that modifies $x will modify $y (except overloaded math operators), and vice
2792versa. See L<Math::BigInt> for details and how to avoid that.
58cde26e 2793
2794=item bpow
2795
2796C<bpow()> now modifies the first argument, unlike the old code which left
2797it alone and only returned the result. This is to be consistent with
2798C<badd()> etc. The first will modify $x, the second one won't:
2799
2800 print bpow($x,$i),"\n"; # modify $x
2801 print $x->bpow($i),"\n"; # ditto
2802 print $x ** $i,"\n"; # leave $x alone
2803
2804=back
2805
990fb837 2806=head1 SEE ALSO
2807
2808L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
2809L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
2810
2811The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
2812because they solve the autoupgrading/downgrading issue, at least partly.
2813
2814The package at
2815L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
2816more documentation including a full version history, testcases, empty
2817subclass files and benchmarks.
2818
58cde26e 2819=head1 LICENSE
a5f75d66 2820
58cde26e 2821This program is free software; you may redistribute it and/or modify it under
2822the same terms as Perl itself.
5d7098d5 2823
58cde26e 2824=head1 AUTHORS
5d7098d5 2825
58cde26e 2826Mark Biggar, overloaded interface by Ilya Zakharevich.
990fb837 2827Completely rewritten by Tels http://bloodgate.com in 2001, 2002, and still
2828at it in 2003.
a5f75d66 2829
a5f75d66 2830=cut