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