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