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