[perl #30609] [PATCH] BigInt v1.71 - first try
[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)
1286 if (overload::StrVal($x) eq overload::StrVal($y))
b3abae2a 1287 {
03874afe 1288 $x->bone(); # x/x => 1, rem 0
b3abae2a 1289 }
03874afe 1290 else
58cde26e 1291 {
03874afe 1292
1293 # make copy of $x in case of list context for later reminder calculation
1294 if (wantarray && !$y->is_one())
1295 {
1296 $rem = $x->copy();
1297 }
394e6ffb 1298
03874afe 1299 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
9b924220 1300
03874afe 1301 # check for / +-1 ( +/- 1E0)
1302 if (!$y->is_one())
1303 {
1304 # promote BigInts and it's subclasses (except when already a BigFloat)
1305 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1306
1307 # calculate the result to $scale digits and then round it
1308 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1309 $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1310 $MBI->_div ($x->{_m},$y->{_m}); # a/c
1311
1312 # correct exponent of $x
1313 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1314 # correct for 10**scale
1315 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1316 $x->bnorm(); # remove trailing 0's
1317 }
1318 } # ende else $x != $y
a5f75d66 1319
091c87b1 1320 # shortcut to not run through _find_round_parameters again
990fb837 1321 if (defined $params[0])
ee15d750 1322 {
ef9466ea 1323 delete $x->{_a}; # clear before round
990fb837 1324 $x->bround($params[0],$params[2]); # then round accordingly
ee15d750 1325 }
1326 else
1327 {
ef9466ea 1328 delete $x->{_p}; # clear before round
990fb837 1329 $x->bfround($params[1],$params[2]); # then round accordingly
ee15d750 1330 }
574bacfe 1331 if ($fallback)
1332 {
1333 # clear a/p after round, since user did not request it
ef9466ea 1334 delete $x->{_a}; delete $x->{_p};
574bacfe 1335 }
03874afe 1336
58cde26e 1337 if (wantarray)
1338 {
394e6ffb 1339 if (!$y->is_one())
1340 {
990fb837 1341 $rem->bmod($y,@params); # copy already done
394e6ffb 1342 }
574bacfe 1343 if ($fallback)
1344 {
1345 # clear a/p after round, since user did not request it
ef9466ea 1346 delete $rem->{_a}; delete $rem->{_p};
574bacfe 1347 }
0716bf9b 1348 return ($x,$rem);
58cde26e 1349 }
9393ace2 1350 $x;
58cde26e 1351 }
a0d0e21e 1352
58cde26e 1353sub bmod
1354 {
1355 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
f9a08e12 1356
1357 # set up parameters
1358 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1359 # objectify is costly, so avoid it
1360 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1361 {
1362 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1363 }
a0d0e21e 1364
9b924220 1365 # handle NaN, inf, -inf
61f5c3f5 1366 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1367 {
1368 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
f9a08e12 1369 $x->{sign} = $re->{sign};
1370 $x->{_e} = $re->{_e};
1371 $x->{_m} = $re->{_m};
1372 return $x->round($a,$p,$r,$y);
61f5c3f5 1373 }
9b924220 1374 if ($y->is_zero())
1375 {
1376 return $x->bnan() if $x->is_zero();
1377 return $x;
1378 }
61f5c3f5 1379 return $x->bzero() if $y->is_one() || $x->is_zero();
58cde26e 1380
61f5c3f5 1381 my $cmp = $x->bacmp($y); # equal or $x < $y?
1382 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1383
1384 # only $y of the operands negative?
1385 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1386
1387 $x->{sign} = $y->{sign}; # calc sign first
1388 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1389
9b924220 1390 my $ym = $MBI->_copy($y->{_m});
61f5c3f5 1391
1392 # 2e1 => 20
9b924220 1393 $MBI->_lsft( $ym, $y->{_e}, 10)
1394 if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
61f5c3f5 1395
1396 # if $y has digits after dot
1397 my $shifty = 0; # correct _e of $x by this
9b924220 1398 if ($y->{_es} eq '-') # has digits after dot
61f5c3f5 1399 {
1400 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
9b924220 1401 $shifty = $MBI->_num($y->{_e}); # no more digits after dot
1402 $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
61f5c3f5 1403 }
1404 # $ym is now mantissa of $y based on exponent 0
b3abae2a 1405
61f5c3f5 1406 my $shiftx = 0; # correct _e of $x by this
9b924220 1407 if ($x->{_es} eq '-') # has digits after dot
61f5c3f5 1408 {
1409 # 123.4 % 20 => 1234 % 200
9b924220 1410 $shiftx = $MBI->_num($x->{_e}); # no more digits after dot
1411 $MBI->_lsft($ym, $x->{_e}, 10); # 123 => 1230
61f5c3f5 1412 }
1413 # 123e1 % 20 => 1230 % 20
9b924220 1414 if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
61f5c3f5 1415 {
9b924220 1416 $MBI->_lsft( $x->{_m}, $x->{_e},10); # es => '+' here
61f5c3f5 1417 }
9b924220 1418
1419 $x->{_e} = $MBI->_new($shiftx);
1420 $x->{_es} = '+';
1421 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1422 $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
61f5c3f5 1423
1424 # now mantissas are equalized, exponent of $x is adjusted, so calc result
b3abae2a 1425
9b924220 1426 $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
61f5c3f5 1427
9b924220 1428 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
61f5c3f5 1429 $x->bnorm();
1430
1431 if ($neg != 0) # one of them negative => correct in place
1432 {
1433 my $r = $y - $x;
1434 $x->{_m} = $r->{_m};
1435 $x->{_e} = $r->{_e};
9b924220 1436 $x->{_es} = $r->{_es};
1437 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
61f5c3f5 1438 $x->bnorm();
1439 }
1440
1441 $x->round($a,$p,$r,$y); # round and return
58cde26e 1442 }
1443
990fb837 1444sub broot
1445 {
1446 # calculate $y'th root of $x
3a427a11 1447
1448 # set up parameters
1449 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1450 # objectify is costly, so avoid it
1451 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1452 {
1453 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1454 }
990fb837 1455
1456 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1457 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1458 $y->{sign} !~ /^\+$/;
1459
1460 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1461
1462 # we need to limit the accuracy to protect against overflow
1463 my $fallback = 0;
1464 my (@params,$scale);
1465 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1466
1467 return $x if $x->is_nan(); # error in _find_round_parameters?
1468
1469 # no rounding at all, so must use fallback
1470 if (scalar @params == 0)
1471 {
1472 # simulate old behaviour
1473 $params[0] = $self->div_scale(); # and round to it as accuracy
1474 $scale = $params[0]+4; # at least four more for proper round
9b924220 1475 $params[2] = $r; # iound mode by caller or undef
990fb837 1476 $fallback = 1; # to clear a/p afterwards
1477 }
1478 else
1479 {
1480 # the 4 below is empirical, and there might be cases where it is not
1481 # enough...
1482 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1483 }
1484
1485 # when user set globals, they would interfere with our calculation, so
1486 # disable them and later re-enable them
1487 no strict 'refs';
1488 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1489 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1490 # we also need to disable any set A or P on $x (_find_round_parameters took
1491 # them already into account), since these would interfere, too
1492 delete $x->{_a}; delete $x->{_p};
1493 # need to disable $upgrade in BigInt, to avoid deep recursion
1494 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1495
1496 # remember sign and make $x positive, since -4 ** (1/2) => -2
27e7b8bb 1497 my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
1498
1499 my $is_two = 0;
1500 if ($y->isa('Math::BigFloat'))
1501 {
1502 $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
1503 }
1504 else
1505 {
1506 $is_two = ($y == 2);
1507 }
990fb837 1508
27e7b8bb 1509 # normal square root if $y == 2:
1510 if ($is_two)
990fb837 1511 {
1512 $x->bsqrt($scale+4);
1513 }
1514 elsif ($y->is_one('-'))
1515 {
1516 # $x ** -1 => 1/$x
1517 my $u = $self->bone()->bdiv($x,$scale);
1518 # copy private parts over
1519 $x->{_m} = $u->{_m};
1520 $x->{_e} = $u->{_e};
9b924220 1521 $x->{_es} = $u->{_es};
990fb837 1522 }
1523 else
1524 {
3a427a11 1525 # calculate the broot() as integer result first, and if it fits, return
1526 # it rightaway (but only if $x and $y are integer):
1527
1528 my $done = 0; # not yet
1529 if ($y->is_int() && $x->is_int())
1530 {
9b924220 1531 my $i = $MBI->_copy( $x->{_m} );
1532 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1533 my $int = Math::BigInt->bzero();
1534 $int->{value} = $i;
3a427a11 1535 $int->broot($y->as_number());
1536 # if ($exact)
1537 if ($int->copy()->bpow($y) == $x)
1538 {
1539 # found result, return it
9b924220 1540 $x->{_m} = $int->{value};
1541 $x->{_e} = $MBI->_zero();
1542 $x->{_es} = '+';
3a427a11 1543 $x->bnorm();
1544 $done = 1;
1545 }
1546 }
1547 if ($done == 0)
1548 {
1549 my $u = $self->bone()->bdiv($y,$scale+4);
1550 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
1551 $x->bpow($u,$scale+4); # el cheapo
1552 }
990fb837 1553 }
1554 $x->bneg() if $sign == 1;
1555
091c87b1 1556 # shortcut to not run through _find_round_parameters again
990fb837 1557 if (defined $params[0])
1558 {
1559 $x->bround($params[0],$params[2]); # then round accordingly
1560 }
1561 else
1562 {
1563 $x->bfround($params[1],$params[2]); # then round accordingly
1564 }
1565 if ($fallback)
1566 {
1567 # clear a/p after round, since user did not request it
ef9466ea 1568 delete $x->{_a}; delete $x->{_p};
990fb837 1569 }
1570 # restore globals
1571 $$abr = $ab; $$pbr = $pb;
1572 $x;
1573 }
1574
58cde26e 1575sub bsqrt
1576 {
990fb837 1577 # calculate square root
ee15d750 1578 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 1579
990fb837 1580 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
1581 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
1582 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
58cde26e 1583
61f5c3f5 1584 # we need to limit the accuracy to protect against overflow
574bacfe 1585 my $fallback = 0;
990fb837 1586 my (@params,$scale);
1587 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1588
1589 return $x if $x->is_nan(); # error in _find_round_parameters?
61f5c3f5 1590
1591 # no rounding at all, so must use fallback
990fb837 1592 if (scalar @params == 0)
0716bf9b 1593 {
1594 # simulate old behaviour
990fb837 1595 $params[0] = $self->div_scale(); # and round to it as accuracy
1596 $scale = $params[0]+4; # at least four more for proper round
1597 $params[2] = $r; # round mode by caller or undef
ee15d750 1598 $fallback = 1; # to clear a/p afterwards
0716bf9b 1599 }
61f5c3f5 1600 else
1601 {
1602 # the 4 below is empirical, and there might be cases where it is not
1603 # enough...
990fb837 1604 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
61f5c3f5 1605 }
1606
1607 # when user set globals, they would interfere with our calculation, so
9393ace2 1608 # disable them and later re-enable them
61f5c3f5 1609 no strict 'refs';
1610 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
b3abae2a 1611 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
61f5c3f5 1612 # we also need to disable any set A or P on $x (_find_round_parameters took
1613 # them already into account), since these would interfere, too
1614 delete $x->{_a}; delete $x->{_p};
9393ace2 1615 # need to disable $upgrade in BigInt, to avoid deep recursion
1616 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
61f5c3f5 1617
9b924220 1618 my $i = $MBI->_copy( $x->{_m} );
1619 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1620 my $xas = Math::BigInt->bzero();
1621 $xas->{value} = $i;
1622
394e6ffb 1623 my $gs = $xas->copy()->bsqrt(); # some guess
b3abae2a 1624
9b924220 1625 if (($x->{_es} ne '-') # guess can't be accurate if there are
394e6ffb 1626 # digits after the dot
b3abae2a 1627 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
394e6ffb 1628 {
9b924220 1629 # exact result, copy result over to keep $x
1630 $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
1631 $x->bnorm();
091c87b1 1632 # shortcut to not run through _find_round_parameters again
990fb837 1633 if (defined $params[0])
61f5c3f5 1634 {
990fb837 1635 $x->bround($params[0],$params[2]); # then round accordingly
61f5c3f5 1636 }
1637 else
1638 {
990fb837 1639 $x->bfround($params[1],$params[2]); # then round accordingly
61f5c3f5 1640 }
1641 if ($fallback)
1642 {
1643 # clear a/p after round, since user did not request it
ef9466ea 1644 delete $x->{_a}; delete $x->{_p};
61f5c3f5 1645 }
9393ace2 1646 # re-enable A and P, upgrade is taken care of by "local"
b3abae2a 1647 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
61f5c3f5 1648 return $x;
394e6ffb 1649 }
2ab5f49d 1650
1651 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
1652 # of the result by multipyling the input by 100 and then divide the integer
1653 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
9b924220 1654
1655 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
1656 my $y1 = $MBI->_copy($x->{_m});
1657
1658 my $length = $MBI->_len($y1);
1659
1660 # Now calculate how many digits the result of sqrt(y1) would have
1661 my $digits = int($length / 2);
1662
1663 # But we need at least $scale digits, so calculate how many are missing
1664 my $shift = $scale - $digits;
1665
1666 # That should never happen (we take care of integer guesses above)
1667 # $shift = 0 if $shift < 0;
1668
1669 # Multiply in steps of 100, by shifting left two times the "missing" digits
1670 my $s2 = $shift * 2;
1671
2ab5f49d 1672 # We now make sure that $y1 has the same odd or even number of digits than
1673 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
1674 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
1675 # steps of 10. The length of $x does not count, since an even or odd number
1676 # of digits before the dot is not changed by adding an even number of digits
1677 # after the dot (the result is still odd or even digits long).
9b924220 1678 $s2++ if $MBI->_is_odd($x->{_e});
1679
1680 $MBI->_lsft( $y1, $MBI->_new($s2), 10);
1681
2ab5f49d 1682 # now take the square root and truncate to integer
9b924220 1683 $y1 = $MBI->_sqrt($y1);
1684
2ab5f49d 1685 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
1686 # result, which is than later rounded to the desired scale.
990fb837 1687
1688 # calculate how many zeros $x had after the '.' (or before it, depending
9b924220 1689 # on sign of $dat, the result should have half as many:
1690 my $dat = $MBI->_num($x->{_e});
1691 $dat = -$dat if $x->{_es} eq '-';
1692 $dat += $length;
990fb837 1693
1694 if ($dat > 0)
1695 {
1696 # no zeros after the dot (e.g. 1.23, 0.49 etc)
1697 # preserve half as many digits before the dot than the input had
1698 # (but round this "up")
1699 $dat = int(($dat+1)/2);
1700 }
1701 else
1702 {
1703 $dat = int(($dat)/2);
1704 }
9b924220 1705 $dat -= $MBI->_len($y1);
1706 if ($dat < 0)
1707 {
1708 $dat = abs($dat);
1709 $x->{_e} = $MBI->_new( $dat );
1710 $x->{_es} = '-';
1711 }
1712 else
1713 {
1714 $x->{_e} = $MBI->_new( $dat );
1715 $x->{_es} = '+';
1716 }
2ab5f49d 1717 $x->{_m} = $y1;
9b924220 1718 $x->bnorm();
61f5c3f5 1719
091c87b1 1720 # shortcut to not run through _find_round_parameters again
990fb837 1721 if (defined $params[0])
61f5c3f5 1722 {
990fb837 1723 $x->bround($params[0],$params[2]); # then round accordingly
61f5c3f5 1724 }
1725 else
1726 {
990fb837 1727 $x->bfround($params[1],$params[2]); # then round accordingly
61f5c3f5 1728 }
574bacfe 1729 if ($fallback)
1730 {
1731 # clear a/p after round, since user did not request it
ef9466ea 1732 delete $x->{_a}; delete $x->{_p};
574bacfe 1733 }
61f5c3f5 1734 # restore globals
b3abae2a 1735 $$abr = $ab; $$pbr = $pb;
574bacfe 1736 $x;
58cde26e 1737 }
1738
b3abae2a 1739sub bfac
1740 {
28df3e88 1741 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
091c87b1 1742 # compute factorial number, modifies first argument
b3abae2a 1743
b282a552 1744 # set up parameters
1745 my ($self,$x,@r) = (ref($_[0]),@_);
1746 # objectify is costly, so avoid it
1747 ($self,$x,@r) = objectify(1,@_) if !ref($x);
1748
1749 return $x if $x->{sign} eq '+inf'; # inf => inf
28df3e88 1750 return $x->bnan()
1751 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
9b924220 1752 ($x->{_es} ne '+')); # digits after dot?
b3abae2a 1753
b3abae2a 1754 # use BigInt's bfac() for faster calc
9b924220 1755 if (! $MBI->_is_zero($x->{_e}))
091c87b1 1756 {
9b924220 1757 $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
1758 $x->{_e} = $MBI->_zero(); # normalize
1759 $x->{_es} = '+';
091c87b1 1760 }
9b924220 1761 $MBI->_fac($x->{_m}); # calculate factorial
091c87b1 1762 $x->bnorm()->round(@r); # norm again and round result
b3abae2a 1763 }
1764
9393ace2 1765sub _pow
1766 {
1767 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1768 my ($x,$y,$a,$p,$r) = @_;
1769 my $self = ref($x);
1770
1771 # if $y == 0.5, it is sqrt($x)
27e7b8bb 1772 $HALF = $self->new($HALF) unless ref($HALF);
1773 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp($HALF) == 0;
9393ace2 1774
990fb837 1775 # Using:
1776 # a ** x == e ** (x * ln a)
1777
9393ace2 1778 # u = y * ln x
990fb837 1779 # _ _
1780 # Taylor: | u u^2 u^3 |
1781 # x ** y = 1 + | --- + --- + ----- + ... |
1782 # |_ 1 1*2 1*2*3 _|
9393ace2 1783
1784 # we need to limit the accuracy to protect against overflow
1785 my $fallback = 0;
990fb837 1786 my ($scale,@params);
1787 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1788
1789 return $x if $x->is_nan(); # error in _find_round_parameters?
9393ace2 1790
1791 # no rounding at all, so must use fallback
990fb837 1792 if (scalar @params == 0)
9393ace2 1793 {
1794 # simulate old behaviour
990fb837 1795 $params[0] = $self->div_scale(); # and round to it as accuracy
1796 $params[1] = undef; # disable P
1797 $scale = $params[0]+4; # at least four more for proper round
1798 $params[2] = $r; # round mode by caller or undef
9393ace2 1799 $fallback = 1; # to clear a/p afterwards
1800 }
1801 else
1802 {
1803 # the 4 below is empirical, and there might be cases where it is not
1804 # enough...
990fb837 1805 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
9393ace2 1806 }
1807
1808 # when user set globals, they would interfere with our calculation, so
56d9de68 1809 # disable them and later re-enable them
9393ace2 1810 no strict 'refs';
1811 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1812 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1813 # we also need to disable any set A or P on $x (_find_round_parameters took
1814 # them already into account), since these would interfere, too
1815 delete $x->{_a}; delete $x->{_p};
1816 # need to disable $upgrade in BigInt, to avoid deep recursion
1817 local $Math::BigInt::upgrade = undef;
1818
1819 my ($limit,$v,$u,$below,$factor,$next,$over);
1820
990fb837 1821 $u = $x->copy()->blog(undef,$scale)->bmul($y);
9393ace2 1822 $v = $self->bone(); # 1
1823 $factor = $self->new(2); # 2
1824 $x->bone(); # first term: 1
1825
1826 $below = $v->copy();
1827 $over = $u->copy();
1828
1829 $limit = $self->new("1E-". ($scale-1));
1830 #my $steps = 0;
1831 while (3 < 5)
1832 {
1833 # we calculate the next term, and add it to the last
1834 # when the next term is below our limit, it won't affect the outcome
1835 # anymore, so we stop
1836 $next = $over->copy()->bdiv($below,$scale);
990fb837 1837 last if $next->bacmp($limit) <= 0;
9393ace2 1838 $x->badd($next);
9393ace2 1839 # calculate things for the next term
1840 $over *= $u; $below *= $factor; $factor->binc();
9b924220 1841
1842 last if $x->{sign} !~ /^[-+]$/;
1843
9393ace2 1844 #$steps++;
1845 }
1846
091c87b1 1847 # shortcut to not run through _find_round_parameters again
990fb837 1848 if (defined $params[0])
9393ace2 1849 {
990fb837 1850 $x->bround($params[0],$params[2]); # then round accordingly
9393ace2 1851 }
1852 else
1853 {
990fb837 1854 $x->bfround($params[1],$params[2]); # then round accordingly
9393ace2 1855 }
1856 if ($fallback)
1857 {
1858 # clear a/p after round, since user did not request it
ef9466ea 1859 delete $x->{_a}; delete $x->{_p};
9393ace2 1860 }
1861 # restore globals
1862 $$abr = $ab; $$pbr = $pb;
1863 $x;
1864 }
1865
58cde26e 1866sub bpow
1867 {
1868 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1869 # compute power of two numbers, second arg is used as integer
1870 # modifies first argument
1871
f9a08e12 1872 # set up parameters
1873 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1874 # objectify is costly, so avoid it
1875 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1876 {
1877 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1878 }
58cde26e 1879
0716bf9b 1880 return $x if $x->{sign} =~ /^[+-]inf$/;
58cde26e 1881 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
574bacfe 1882 return $x->bone() if $y->is_zero();
58cde26e 1883 return $x if $x->is_one() || $y->is_one();
9393ace2 1884
d614cd8b 1885 return $x->_pow($y,$a,$p,$r) if !$y->is_int(); # non-integer power
9393ace2 1886
9b924220 1887 my $y1 = $y->as_number()->{value}; # make CALC
1888
394e6ffb 1889 # if ($x == -1)
9b924220 1890 if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
58cde26e 1891 {
1892 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
9b924220 1893 return $MBI->_is_odd($y1) ? $x : $x->babs(1);
288d023a 1894 }
28df3e88 1895 if ($x->is_zero())
1896 {
9b924220 1897 return $x->bone() if $y->is_zero();
28df3e88 1898 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
9b924220 1899 # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
1900 return $x->binf();
28df3e88 1901 }
58cde26e 1902
9b924220 1903 my $new_sign = '+';
1904 $new_sign = $y->is_odd() ? '-' : '+' if ($x->{sign} ne '+');
1905
58cde26e 1906 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
9b924220 1907 $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
1908 $MBI->_mul ($x->{_e}, $y1);
1909
1910 $x->{sign} = $new_sign;
58cde26e 1911 $x->bnorm();
1912 if ($y->{sign} eq '-')
1913 {
1914 # modify $x in place!
0716bf9b 1915 my $z = $x->copy(); $x->bzero()->binc();
58cde26e 1916 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
a0d0e21e 1917 }
28df3e88 1918 $x->round($a,$p,$r,$y);
58cde26e 1919 }
1920
1921###############################################################################
1922# rounding functions
1923
1924sub bfround
1925 {
1926 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1927 # $n == 0 means round to integer
1928 # expects and returns normalized numbers!
ee15d750 1929 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
a0d0e21e 1930
58cde26e 1931 return $x if $x->modify('bfround');
9b924220 1932
ee15d750 1933 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
58cde26e 1934 return $x if !defined $scale; # no-op
1935
574bacfe 1936 # never round a 0, +-inf, NaN
61f5c3f5 1937 if ($x->is_zero())
1938 {
1939 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1940 return $x;
1941 }
1942 return $x if $x->{sign} !~ /^[+-]$/;
58cde26e 1943
ee15d750 1944 # don't round if x already has lower precision
1945 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1946
1947 $x->{_p} = $scale; # remember round in any case
ef9466ea 1948 delete $x->{_a}; # and clear A
58cde26e 1949 if ($scale < 0)
1950 {
58cde26e 1951 # round right from the '.'
f9a08e12 1952
9b924220 1953 return $x if $x->{_es} eq '+'; # e >= 0 => nothing to round
f9a08e12 1954
58cde26e 1955 $scale = -$scale; # positive for simplicity
9b924220 1956 my $len = $MBI->_len($x->{_m}); # length of mantissa
f9a08e12 1957
1958 # the following poses a restriction on _e, but if _e is bigger than a
1959 # scalar, you got other problems (memory etc) anyway
9b924220 1960 my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e}))); # digits after dot
58cde26e 1961 my $zad = 0; # zeros after dot
f9a08e12 1962 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
9b924220 1963
1964 # p rint "scale $scale dad $dad zad $zad len $len\n";
58cde26e 1965 # number bsstr len zad dad
1966 # 0.123 123e-3 3 0 3
1967 # 0.0123 123e-4 3 1 4
1968 # 0.001 1e-3 1 2 3
1969 # 1.23 123e-2 3 0 2
1970 # 1.2345 12345e-4 5 0 4
1971
1972 # do not round after/right of the $dad
1973 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
1974
ee15d750 1975 # round to zero if rounding inside the $zad, but not for last zero like:
1976 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1977 return $x->bzero() if $scale < $zad;
1978 if ($scale == $zad) # for 0.006, scale -3 and trunc
58cde26e 1979 {
b3abae2a 1980 $scale = -$len;
58cde26e 1981 }
1982 else
1983 {
1984 # adjust round-point to be inside mantissa
1985 if ($zad != 0)
1986 {
1987 $scale = $scale-$zad;
1988 }
1989 else
1990 {
1991 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
1992 $scale = $dbd+$scale;
1993 }
1994 }
a0d0e21e 1995 }
58cde26e 1996 else
1997 {
f9a08e12 1998 # round left from the '.'
1999
58cde26e 2000 # 123 => 100 means length(123) = 3 - $scale (2) => 1
a5f75d66 2001
9b924220 2002 my $dbt = $MBI->_len($x->{_m});
b3abae2a 2003 # digits before dot
9b924220 2004 my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
b3abae2a 2005 # should be the same, so treat it as this
2006 $scale = 1 if $scale == 0;
2007 # shortcut if already integer
2008 return $x if $scale == 1 && $dbt <= $dbd;
2009 # maximum digits before dot
2010 ++$dbd;
2011
2012 if ($scale > $dbd)
2013 {
2014 # not enough digits before dot, so round to zero
2015 return $x->bzero;
2016 }
2017 elsif ( $scale == $dbd )
2018 {
2019 # maximum
2020 $scale = -$dbt;
2021 }
58cde26e 2022 else
b3abae2a 2023 {
2024 $scale = $dbd - $scale;
2025 }
a0d0e21e 2026 }
574bacfe 2027 # pass sign to bround for rounding modes '+inf' and '-inf'
9b924220 2028 my $m = Math::BigInt->new( $x->{sign} . $MBI->_str($x->{_m}));
2029 $m->bround($scale,$mode);
2030 $x->{_m} = $m->{value}; # get our mantissa back
58cde26e 2031 $x->bnorm();
2032 }
2033
2034sub bround
2035 {
2036 # accuracy: preserve $N digits, and overwrite the rest with 0's
ee15d750 2037 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
9b924220 2038
990fb837 2039 if (($_[0] || 0) < 0)
2040 {
2041 require Carp; Carp::croak ('bround() needs positive accuracy');
2042 }
58cde26e 2043
ee15d750 2044 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
2045 return $x if !defined $scale; # no-op
61f5c3f5 2046
58cde26e 2047 return $x if $x->modify('bround');
61f5c3f5 2048
ee15d750 2049 # scale is now either $x->{_a}, $accuracy, or the user parameter
2050 # test whether $x already has lower accuracy, do nothing in this case
2051 # but do round if the accuracy is the same, since a math operation might
2052 # want to round a number with A=5 to 5 digits afterwards again
2053 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
58cde26e 2054
61f5c3f5 2055 # scale < 0 makes no sense
2056 # never round a +-inf, NaN
2057 return $x if ($scale < 0) || $x->{sign} !~ /^[+-]$/;
58cde26e 2058
61f5c3f5 2059 # 1: $scale == 0 => keep all digits
2060 # 2: never round a 0
2061 # 3: if we should keep more digits than the mantissa has, do nothing
9b924220 2062 if ($scale == 0 || $x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
61f5c3f5 2063 {
2064 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
2065 return $x;
2066 }
f216259d 2067
58cde26e 2068 # pass sign to bround for '+inf' and '-inf' rounding modes
9b924220 2069 my $m = Math::BigInt->new( $x->{sign} . $MBI->_str($x->{_m}));
2070
2071 $m->bround($scale,$mode); # round mantissa
2072 $x->{_m} = $m->{value}; # get our mantissa back
ee15d750 2073 $x->{_a} = $scale; # remember rounding
ef9466ea 2074 delete $x->{_p}; # and clear P
574bacfe 2075 $x->bnorm(); # del trailing zeros gen. by bround()
58cde26e 2076 }
2077
2078sub bfloor
2079 {
2080 # return integer less or equal then $x
ee15d750 2081 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 2082
2083 return $x if $x->modify('bfloor');
2084
2085 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2086
2087 # if $x has digits after dot
9b924220 2088 if ($x->{_es} eq '-')
58cde26e 2089 {
9b924220 2090 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2091 $x->{_e} = $MBI->_zero(); # trunc/norm
2092 $x->{_es} = '+'; # abs e
2093 $MBI->_inc($x->{_m}) if $x->{sign} eq '-'; # increment if negative
f216259d 2094 }
61f5c3f5 2095 $x->round($a,$p,$r);
58cde26e 2096 }
288d023a 2097
58cde26e 2098sub bceil
2099 {
2100 # return integer greater or equal then $x
ee15d750 2101 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 2102
2103 return $x if $x->modify('bceil');
2104 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2105
2106 # if $x has digits after dot
9b924220 2107 if ($x->{_es} eq '-')
58cde26e 2108 {
9b924220 2109 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2110 $x->{_e} = $MBI->_zero(); # trunc/norm
2111 $x->{_es} = '+'; # abs e
2112 $MBI->_inc($x->{_m}) if $x->{sign} eq '+'; # increment if positive
a0d0e21e 2113 }
61f5c3f5 2114 $x->round($a,$p,$r);
58cde26e 2115 }
2116
394e6ffb 2117sub brsft
2118 {
f9a08e12 2119 # shift right by $y (divide by power of $n)
2120
2121 # set up parameters
2122 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2123 # objectify is costly, so avoid it
2124 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2125 {
2126 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2127 }
394e6ffb 2128
2129 return $x if $x->modify('brsft');
2130 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2131
f9a08e12 2132 $n = 2 if !defined $n; $n = $self->new($n);
2133 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
394e6ffb 2134 }
2135
2136sub blsft
2137 {
f9a08e12 2138 # shift left by $y (multiply by power of $n)
2139
2140 # set up parameters
2141 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2142 # objectify is costly, so avoid it
2143 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2144 {
2145 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2146 }
394e6ffb 2147
f9a08e12 2148 return $x if $x->modify('blsft');
394e6ffb 2149 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
2150
f9a08e12 2151 $n = 2 if !defined $n; $n = $self->new($n);
2152 $x->bmul($n->bpow($y),$a,$p,$r,$y);
394e6ffb 2153 }
2154
58cde26e 2155###############################################################################
a5f75d66 2156
58cde26e 2157sub DESTROY
2158 {
b282a552 2159 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
58cde26e 2160 }
2161
2162sub AUTOLOAD
2163 {
b3abae2a 2164 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
2165 # or falling back to MBI::bxxx()
58cde26e 2166 my $name = $AUTOLOAD;
2167
ef9466ea 2168 $name =~ s/(.*):://; # split package
2169 my $c = $1 || $class;
ee15d750 2170 no strict 'refs';
ef9466ea 2171 $c->import() if $IMPORT == 0;
ee15d750 2172 if (!method_alias($name))
58cde26e 2173 {
ee15d750 2174 if (!defined $name)
2175 {
2176 # delayed load of Carp and avoid recursion
2177 require Carp;
ef9466ea 2178 Carp::croak ("$c: Can't call a method without name");
ee15d750 2179 }
ee15d750 2180 if (!method_hand_up($name))
2181 {
2182 # delayed load of Carp and avoid recursion
2183 require Carp;
ef9466ea 2184 Carp::croak ("Can't call $c\-\>$name, not a valid method");
ee15d750 2185 }
2186 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
2187 $name =~ s/^f/b/;
9b924220 2188 return &{"Math::BigInt"."::$name"}(@_);
a0d0e21e 2189 }
58cde26e 2190 my $bname = $name; $bname =~ s/^f/b/;
ef9466ea 2191 $c .= "::$name";
2192 *{$c} = \&{$bname};
2193 &{$c}; # uses @_
58cde26e 2194 }
2195
2196sub exponent
2197 {
2198 # return a copy of the exponent
ee15d750 2199 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2200
ee15d750 2201 if ($x->{sign} !~ /^[+-]$/)
2202 {
2203 my $s = $x->{sign}; $s =~ s/^[+-]//;
9b924220 2204 return Math::BigInt->new($s); # -inf, +inf => +inf
ee15d750 2205 }
9b924220 2206 Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
58cde26e 2207 }
2208
2209sub mantissa
2210 {
2211 # return a copy of the mantissa
ee15d750 2212 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2213
ee15d750 2214 if ($x->{sign} !~ /^[+-]$/)
2215 {
2216 my $s = $x->{sign}; $s =~ s/^[+]//;
9b924220 2217 return Math::BigInt->new($s); # -inf, +inf => +inf
ee15d750 2218 }
9b924220 2219 my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
ee15d750 2220 $m->bneg() if $x->{sign} eq '-';
58cde26e 2221
61f5c3f5 2222 $m;
58cde26e 2223 }
2224
2225sub parts
2226 {
2227 # return a copy of both the exponent and the mantissa
ee15d750 2228 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2229
ee15d750 2230 if ($x->{sign} !~ /^[+-]$/)
2231 {
2232 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2233 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2234 }
9b924220 2235 my $m = Math::BigInt->bzero();
2236 $m->{value} = $MBI->_copy($x->{_m});
ee15d750 2237 $m->bneg() if $x->{sign} eq '-';
9b924220 2238 ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
58cde26e 2239 }
2240
2241##############################################################################
2242# private stuff (internal use only)
2243
58cde26e 2244sub import
2245 {
2246 my $self = shift;
8f675a64 2247 my $l = scalar @_;
2248 my $lib = ''; my @a;
990fb837 2249 $IMPORT=1;
8f675a64 2250 for ( my $i = 0; $i < $l ; $i++)
58cde26e 2251 {
2252 if ( $_[$i] eq ':constant' )
2253 {
091c87b1 2254 # This causes overlord er load to step in. 'binary' and 'integer'
2255 # are handled by BigInt.
58cde26e 2256 overload::constant float => sub { $self->new(shift); };
b3abae2a 2257 }
2258 elsif ($_[$i] eq 'upgrade')
2259 {
2260 # this causes upgrading
28df3e88 2261 $upgrade = $_[$i+1]; # or undef to disable
8f675a64 2262 $i++;
28df3e88 2263 }
2264 elsif ($_[$i] eq 'downgrade')
2265 {
2266 # this causes downgrading
2267 $downgrade = $_[$i+1]; # or undef to disable
8f675a64 2268 $i++;
58cde26e 2269 }
56b9c951 2270 elsif ($_[$i] eq 'lib')
2271 {
990fb837 2272 # alternative library
56b9c951 2273 $lib = $_[$i+1] || ''; # default Calc
8f675a64 2274 $i++;
56b9c951 2275 }
2276 elsif ($_[$i] eq 'with')
2277 {
990fb837 2278 # alternative class for our private parts()
9b924220 2279 # XXX: no longer supported
2280 # $MBI = $_[$i+1] || 'Math::BigInt';
8f675a64 2281 $i++;
2282 }
2283 else
2284 {
2285 push @a, $_[$i];
56b9c951 2286 }
58cde26e 2287 }
8f675a64 2288
56b9c951 2289 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2290 my $mbilib = eval { Math::BigInt->config()->{lib} };
9b924220 2291 if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
8f675a64 2292 {
2293 # MBI already loaded
9b924220 2294 Math::BigInt->import('lib',"$lib,$mbilib", 'objectify');
8f675a64 2295 }
2296 else
2297 {
9b924220 2298 # MBI not loaded, or with ne "Math::BigInt::Calc"
8f675a64 2299 $lib .= ",$mbilib" if defined $mbilib;
07d34614 2300 $lib =~ s/^,//; # don't leave empty
990fb837 2301 # replacement library can handle lib statement, but also could ignore it
8f675a64 2302 if ($] < 5.006)
2303 {
2304 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2305 # used in the same script, or eval inside import().
9b924220 2306 require Math::BigInt;
2307 Math::BigInt->import( lib => $lib, 'objectify' );
8f675a64 2308 }
2309 else
2310 {
9b924220 2311 my $rc = "use Math::BigInt lib => '$lib', 'objectify';";
8f675a64 2312 eval $rc;
2313 }
2314 }
990fb837 2315 if ($@)
2316 {
9b924220 2317 require Carp; Carp::croak ("Couldn't load $lib: $! $@");
990fb837 2318 }
9b924220 2319 $MBI = Math::BigInt->config()->{lib};
56b9c951 2320
58cde26e 2321 # any non :constant stuff is handled by our parent, Exporter
2322 # even if @_ is empty, to give it a chance
b3abae2a 2323 $self->SUPER::import(@a); # for subclasses
2324 $self->export_to_level(1,$self,@a); # need this, too
58cde26e 2325 }
2326
2327sub bnorm
2328 {
2329 # adjust m and e so that m is smallest possible
2330 # round number according to accuracy and precision settings
9b924220 2331 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
58cde26e 2332
0716bf9b 2333 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
58cde26e 2334
9b924220 2335 my $zeros = $MBI->_zeros($x->{_m}); # correct for trailing zeros
b282a552 2336 if ($zeros != 0)
2337 {
9b924220 2338 my $z = $MBI->_new($zeros);
2339 $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
2340 if ($x->{_es} eq '-')
2341 {
2342 if ($MBI->_acmp($x->{_e},$z) >= 0)
2343 {
2344 $x->{_e} = $MBI->_sub ($x->{_e}, $z);
27e7b8bb 2345 $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
9b924220 2346 }
2347 else
2348 {
2349 $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
2350 $x->{_es} = '+';
2351 }
2352 }
2353 else
2354 {
2355 $x->{_e} = $MBI->_add ($x->{_e}, $z);
2356 }
b282a552 2357 }
2358 else
2359 {
2360 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
2361 # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
9b924220 2362 $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
2363 if $MBI->_is_zero($x->{_m});
b282a552 2364 }
2365
61f5c3f5 2366 $x; # MBI bnorm is no-op, so dont call it
2367 }
58cde26e 2368
2369##############################################################################
56d9de68 2370
2371sub as_hex
2372 {
2373 # return number as hexadecimal string (only for integers defined)
2374 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2375
2376 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2377 return '0x0' if $x->is_zero();
2378
9b924220 2379 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
56d9de68 2380
9b924220 2381 my $z = $MBI->_copy($x->{_m});
2382 if (! $MBI->_is_zero($x->{_e})) # > 0
56d9de68 2383 {
9b924220 2384 $MBI->_lsft( $z, $x->{_e},10);
56d9de68 2385 }
9b924220 2386 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
56d9de68 2387 $z->as_hex();
2388 }
2389
2390sub as_bin
2391 {
2392 # return number as binary digit string (only for integers defined)
2393 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2394
2395 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2396 return '0b0' if $x->is_zero();
2397
9b924220 2398 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex!?
56d9de68 2399
9b924220 2400 my $z = $MBI->_copy($x->{_m});
2401 if (! $MBI->_is_zero($x->{_e})) # > 0
56d9de68 2402 {
9b924220 2403 $MBI->_lsft( $z, $x->{_e},10);
56d9de68 2404 }
9b924220 2405 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
56d9de68 2406 $z->as_bin();
2407 }
58cde26e 2408
2409sub as_number
2410 {
394e6ffb 2411 # return copy as a bigint representation of this BigFloat number
2412 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2413
9b924220 2414 my $z = $MBI->_copy($x->{_m});
2415 if ($x->{_es} eq '-') # < 0
58cde26e 2416 {
9b924220 2417 $MBI->_rsft( $z, $x->{_e},10);
0716bf9b 2418 }
9b924220 2419 elsif (! $MBI->_is_zero($x->{_e})) # > 0
0716bf9b 2420 {
9b924220 2421 $MBI->_lsft( $z, $x->{_e},10);
58cde26e 2422 }
9b924220 2423 $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
61f5c3f5 2424 $z;
58cde26e 2425 }
2426
2427sub length
2428 {
ee15d750 2429 my $x = shift;
2430 my $class = ref($x) || $x;
2431 $x = $class->new(shift) unless ref($x);
58cde26e 2432
9b924220 2433 return 1 if $MBI->_is_zero($x->{_m});
2434
2435 my $len = $MBI->_len($x->{_m});
2436 $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
58cde26e 2437 if (wantarray())
2438 {
9b924220 2439 my $t = 0;
2440 $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
2441 return ($len, $t);
58cde26e 2442 }
61f5c3f5 2443 $len;
58cde26e 2444 }
a0d0e21e 2445
24461;
a5f75d66 2447__END__
2448
2449=head1 NAME
2450
58cde26e 2451Math::BigFloat - Arbitrary size floating point math package
a5f75d66 2452
2453=head1 SYNOPSIS
2454
a2008d6d 2455 use Math::BigFloat;
58cde26e 2456
b3abae2a 2457 # Number creation
2458 $x = Math::BigFloat->new($str); # defaults to 0
2459 $nan = Math::BigFloat->bnan(); # create a NotANumber
2460 $zero = Math::BigFloat->bzero(); # create a +0
2461 $inf = Math::BigFloat->binf(); # create a +inf
2462 $inf = Math::BigFloat->binf('-'); # create a -inf
2463 $one = Math::BigFloat->bone(); # create a +1
2464 $one = Math::BigFloat->bone('-'); # create a -1
58cde26e 2465
2466 # Testing
b3abae2a 2467 $x->is_zero(); # true if arg is +0
2468 $x->is_nan(); # true if arg is NaN
0716bf9b 2469 $x->is_one(); # true if arg is +1
2470 $x->is_one('-'); # true if arg is -1
2471 $x->is_odd(); # true if odd, false for even
2472 $x->is_even(); # true if even, false for odd
ef9466ea 2473 $x->is_pos(); # true if >= 0
2474 $x->is_neg(); # true if < 0
b3abae2a 2475 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
2476
58cde26e 2477 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
2478 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
2479 $x->sign(); # return the sign, either +,- or NaN
b3abae2a 2480 $x->digit($n); # return the nth digit, counting from right
2481 $x->digit(-$n); # return the nth digit, counting from left
58cde26e 2482
990fb837 2483 # The following all modify their first argument. If you want to preserve
2484 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2485 # neccessary when mixing $a = $b assigments with non-overloaded math.
2486
58cde26e 2487 # set
2488 $x->bzero(); # set $i to 0
2489 $x->bnan(); # set $i to NaN
b3abae2a 2490 $x->bone(); # set $x to +1
2491 $x->bone('-'); # set $x to -1
2492 $x->binf(); # set $x to inf
2493 $x->binf('-'); # set $x to -inf
58cde26e 2494
2495 $x->bneg(); # negation
2496 $x->babs(); # absolute value
2497 $x->bnorm(); # normalize (no-op)
2498 $x->bnot(); # two's complement (bit wise not)
2499 $x->binc(); # increment x by 1
2500 $x->bdec(); # decrement x by 1
2501
2502 $x->badd($y); # addition (add $y to $x)
2503 $x->bsub($y); # subtraction (subtract $y from $x)
2504 $x->bmul($y); # multiplication (multiply $x by $y)
990fb837 2505 $x->bdiv($y); # divide, set $x to quotient
58cde26e 2506 # return (quo,rem) or quo if scalar
2507
990fb837 2508 $x->bmod($y); # modulus ($x % $y)
2509 $x->bpow($y); # power of arguments ($x ** $y)
58cde26e 2510 $x->blsft($y); # left shift
2511 $x->brsft($y); # right shift
2512 # return (quo,rem) or quo if scalar
2513
990fb837 2514 $x->blog(); # logarithm of $x to base e (Euler's number)
2515 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
61f5c3f5 2516
58cde26e 2517 $x->band($y); # bit-wise and
2518 $x->bior($y); # bit-wise inclusive or
2519 $x->bxor($y); # bit-wise exclusive or
2520 $x->bnot(); # bit-wise not (two's complement)
b3abae2a 2521
2522 $x->bsqrt(); # calculate square-root
990fb837 2523 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
b3abae2a 2524 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
2525
990fb837 2526 $x->bround($N); # accuracy: preserve $N digits
58cde26e 2527 $x->bfround($N); # precision: round to the $Nth digit
2528
990fb837 2529 $x->bfloor(); # return integer less or equal than $x
2530 $x->bceil(); # return integer greater or equal than $x
2531
58cde26e 2532 # The following do not modify their arguments:
990fb837 2533
58cde26e 2534 bgcd(@values); # greatest common divisor
2535 blcm(@values); # lowest common multiplicator
2536
2537 $x->bstr(); # return string
2538 $x->bsstr(); # return string in scientific notation
ef9466ea 2539
2540 $x->as_int(); # return $x as BigInt
58cde26e 2541 $x->exponent(); # return exponent as BigInt
2542 $x->mantissa(); # return mantissa as BigInt
2543 $x->parts(); # return (mantissa,exponent) as BigInt
2544
2545 $x->length(); # number of digits (w/o sign and '.')
2546 ($l,$f) = $x->length(); # number of digits, and length of fraction
a5f75d66 2547
f9a08e12 2548 $x->precision(); # return P of $x (or global, if P of $x undef)
2549 $x->precision($n); # set P of $x to $n
2550 $x->accuracy(); # return A of $x (or global, if A of $x undef)
723d369b 2551 $x->accuracy($n); # set A $x to $n
f9a08e12 2552
990fb837 2553 # these get/set the appropriate global value for all BigFloat objects
2554 Math::BigFloat->precision(); # Precision
2555 Math::BigFloat->accuracy(); # Accuracy
2556 Math::BigFloat->round_mode(); # rounding mode
f9a08e12 2557
a5f75d66 2558=head1 DESCRIPTION
2559
58cde26e 2560All operators (inlcuding basic math operations) are overloaded if you
2561declare your big floating point numbers as
a5f75d66 2562
58cde26e 2563 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2564
2565Operations with overloaded operators preserve the arguments, which is
2566exactly what you expect.
2567
2568=head2 Canonical notation
2569
2570Input to these routines are either BigFloat objects, or strings of the
2571following four forms:
a5f75d66 2572
2573=over 2
2574
58cde26e 2575=item *
2576
2577C</^[+-]\d+$/>
a5f75d66 2578
58cde26e 2579=item *
a5f75d66 2580
58cde26e 2581C</^[+-]\d+\.\d*$/>
a5f75d66 2582
58cde26e 2583=item *
a5f75d66 2584
58cde26e 2585C</^[+-]\d+E[+-]?\d+$/>
a5f75d66 2586
58cde26e 2587=item *
a5f75d66 2588
58cde26e 2589C</^[+-]\d*\.\d+E[+-]?\d+$/>
5d7098d5 2590
58cde26e 2591=back
2592
2593all with optional leading and trailing zeros and/or spaces. Additonally,
2594numbers are allowed to have an underscore between any two digits.
2595
2596Empty strings as well as other illegal numbers results in 'NaN'.
2597
2598bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2599are always stored in normalized form. On a string, it creates a BigFloat
2600object.
2601
2602=head2 Output
2603
2604Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2605
2606The string output will always have leading and trailing zeros stripped and drop
2607a plus sign. C<bstr()> will give you always the form with a decimal point,
990fb837 2608while C<bsstr()> (s for scientific) gives you the scientific notation.
58cde26e 2609
2610 Input bstr() bsstr()
2611 '-0' '0' '0E1'
2612 ' -123 123 123' '-123123123' '-123123123E0'
2613 '00.0123' '0.0123' '123E-4'
2614 '123.45E-2' '1.2345' '12345E-4'
2615 '10E+3' '10000' '1E4'
2616
2617Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2618C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2619return either undef, <0, 0 or >0 and are suited for sort.
2620
990fb837 2621Actual math is done by using the class defined with C<with => Class;> (which
2622defaults to BigInts) to represent the mantissa and exponent.
2623
58cde26e 2624The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2625represent the result when input arguments are not numbers, as well as
2626the result of dividing by zero.
2627
2628=head2 C<mantissa()>, C<exponent()> and C<parts()>
2629
2630C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2631as BigInts such that:
2632
2633 $m = $x->mantissa();
2634 $e = $x->exponent();
2635 $y = $m * ( 10 ** $e );
2636 print "ok\n" if $x == $y;
2637
2638C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2639
2640A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2641
2642Currently the mantissa is reduced as much as possible, favouring higher
2643exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2644This might change in the future, so do not depend on it.
2645
2646=head2 Accuracy vs. Precision
2647
2648See also: L<Rounding|Rounding>.
2649
027dc388 2650Math::BigFloat supports both precision and accuracy. For a full documentation,
2651examples and tips on these topics please see the large section in
2652L<Math::BigInt>.
5d7098d5 2653
58cde26e 2654Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2655a operation consumes all resources, each operation produces no more than
990fb837 2656the requested number of digits.
2657
2658Please refer to BigInt's documentation for the precedence rules of which
2659accuracy/precision setting will be used.
2660
2661If there is no gloabl precision set, B<and> the operation inquestion was not
2662called with a requested precision or accuracy, B<and> the input $x has no
2663accuracy or precision set, then a fallback parameter will be used. For
2664historical reasons, it is called C<div_scale> and can be accessed via:
2665
2666 $d = Math::BigFloat->div_scale(); # query
2667 Math::BigFloat->div_scale($n); # set to $n digits
2668
2669The default value is 40 digits.
58cde26e 2670
2671In case the result of one operation has more precision than specified,
2672it is rounded. The rounding mode taken is either the default mode, or the one
2673supplied to the operation after the I<scale>:
2674
2675 $x = Math::BigFloat->new(2);
990fb837 2676 Math::BigFloat->precision(5); # 5 digits max
58cde26e 2677 $y = $x->copy()->bdiv(3); # will give 0.66666
2678 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2679 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
990fb837 2680 Math::BigFloat->round_mode('zero');
58cde26e 2681 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2682
2683=head2 Rounding
2684
2685=over 2
2686
5dc6f178 2687=item ffround ( +$scale )
58cde26e 2688
0716bf9b 2689Rounds to the $scale'th place left from the '.', counting from the dot.
2690The first digit is numbered 1.
58cde26e 2691
5dc6f178 2692=item ffround ( -$scale )
58cde26e 2693
0716bf9b 2694Rounds to the $scale'th place right from the '.', counting from the dot.
58cde26e 2695
5dc6f178 2696=item ffround ( 0 )
2697
0716bf9b 2698Rounds to an integer.
5dc6f178 2699
2700=item fround ( +$scale )
2701
0716bf9b 2702Preserves accuracy to $scale digits from the left (aka significant digits)
2703and pads the rest with zeros. If the number is between 1 and -1, the
2704significant digits count from the first non-zero after the '.'
5dc6f178 2705
2706=item fround ( -$scale ) and fround ( 0 )
2707
990fb837 2708These are effectively no-ops.
5d7098d5 2709
a5f75d66 2710=back
2711
0716bf9b 2712All rounding functions take as a second parameter a rounding mode from one of
2713the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
58cde26e 2714
2715The default rounding mode is 'even'. By using
990fb837 2716C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
ee15d750 2717mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
0716bf9b 2718no longer supported.
b22b3e31 2719The second parameter to the round functions then overrides the default
0716bf9b 2720temporarily.
58cde26e 2721
990fb837 2722The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
58cde26e 2723'trunc' as rounding mode to make it equivalent to:
2724
2725 $x = 2.5;
2726 $y = int($x) + 2;
2727
2728You can override this by passing the desired rounding mode as parameter to
2729C<as_number()>:
2730
2731 $x = Math::BigFloat->new(2.5);
2732 $y = $x->as_number('odd'); # $y = 3
2733
2734=head1 EXAMPLES
2735
58cde26e 2736 # not ready yet
58cde26e 2737
2738=head1 Autocreating constants
2739
2740After C<use Math::BigFloat ':constant'> all the floating point constants
2741in the given scope are converted to C<Math::BigFloat>. This conversion
2742happens at compile time.
2743
2744In particular
2745
2746 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2747
56b9c951 2748prints the value of C<2E-100>. Note that without conversion of
58cde26e 2749constants the expression 2E-100 will be calculated as normal floating point
2750number.
2751
56b9c951 2752Please note that ':constant' does not affect integer constants, nor binary
2753nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2754work.
2755
2756=head2 Math library
2757
2758Math with the numbers is done (by default) by a module called
2759Math::BigInt::Calc. This is equivalent to saying:
2760
2761 use Math::BigFloat lib => 'Calc';
2762
2763You can change this by using:
2764
2765 use Math::BigFloat lib => 'BitVect';
2766
2767The following would first try to find Math::BigInt::Foo, then
2768Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2769
2770 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2771
2772Calc.pm uses as internal format an array of elements of some decimal base
2773(usually 1e7, but this might be differen for some systems) with the least
2774significant digit first, while BitVect.pm uses a bit vector of base 2, most
2775significant bit first. Other modules might use even different means of
2776representing the numbers. See the respective module documentation for further
2777details.
2778
2779Please note that Math::BigFloat does B<not> use the denoted library itself,
2780but it merely passes the lib argument to Math::BigInt. So, instead of the need
2781to do:
2782
2783 use Math::BigInt lib => 'GMP';
2784 use Math::BigFloat;
2785
2786you can roll it all into one line:
2787
2788 use Math::BigFloat lib => 'GMP';
2789
990fb837 2790It is also possible to just require Math::BigFloat:
2791
2792 require Math::BigFloat;
2793
2794This will load the neccessary things (like BigInt) when they are needed, and
2795automatically.
2796
2797Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
2798you ever wanted to know about loading a different library.
56b9c951 2799
2800=head2 Using Math::BigInt::Lite
2801
2802It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2803
2804 # 1
2805 use Math::BigFloat with => 'Math::BigInt::Lite';
2806
2807There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2808can combine these if you want. For instance, you may want to use
2809Math::BigInt objects in your main script, too.
2810
2811 # 2
2812 use Math::BigInt;
2813 use Math::BigFloat with => 'Math::BigInt::Lite';
2814
2815Of course, you can combine this with the C<lib> parameter.
2816
2817 # 3
2818 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2819
990fb837 2820There is no need for a "use Math::BigInt;" statement, even if you want to
2821use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
2822always loads it. But if you add it, add it B<before>:
56b9c951 2823
2824 # 4
2825 use Math::BigInt;
2826 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2827
2828Notice that the module with the last C<lib> will "win" and thus
2829it's lib will be used if the lib is available:
2830
2831 # 5
2832 use Math::BigInt lib => 'Bar,Baz';
2833 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2834
2835That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2836words, Math::BigFloat will try to retain previously loaded libs when you
990fb837 2837don't specify it onem but if you specify one, it will try to load them.
56b9c951 2838
2839Actually, the lib loading order would be "Bar,Baz,Calc", and then
2840"Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
990fb837 2841same as trying the latter load alone, except for the fact that one of Bar or
2842Baz might be loaded needlessly in an intermidiate step (and thus hang around
2843and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
2844will still be tried to be loaded, but this is not as time/memory consuming as
2845actually loading one of them. Still, this type of usage is not recommended due
2846to these issues.
56b9c951 2847
990fb837 2848The old way (loading the lib only in BigInt) still works though:
56b9c951 2849
2850 # 6
2851 use Math::BigInt lib => 'Bar,Baz';
2852 use Math::BigFloat;
2853
990fb837 2854You can even load Math::BigInt afterwards:
56b9c951 2855
990fb837 2856 # 7
2857 use Math::BigFloat;
2858 use Math::BigInt lib => 'Bar,Baz';
a5f75d66 2859
990fb837 2860But this has the same problems like #5, it will first load Calc
2861(Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
2862Baz, depending on which of them works and is usable/loadable. Since this
2863loads Calc unnecc., it is not recommended.
58cde26e 2864
990fb837 2865Since it also possible to just require Math::BigFloat, this poses the question
2866about what libary this will use:
58cde26e 2867
990fb837 2868 require Math::BigFloat;
2869 my $x = Math::BigFloat->new(123); $x += 123;
58cde26e 2870
990fb837 2871It will use Calc. Please note that the call to import() is still done, but
2872only when you use for the first time some Math::BigFloat math (it is triggered
2873via any constructor, so the first time you create a Math::BigFloat, the load
2874will happen in the background). This means:
58cde26e 2875
990fb837 2876 require Math::BigFloat;
2877 Math::BigFloat->import ( lib => 'Foo,Bar' );
58cde26e 2878
990fb837 2879would be the same as:
58cde26e 2880
990fb837 2881 use Math::BigFloat lib => 'Foo, Bar';
2882
2883But don't try to be clever to insert some operations in between:
2884
2885 require Math::BigFloat;
2886 my $x = Math::BigFloat->bone() + 4; # load BigInt and Calc
2887 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
2888 $x = Math::BigFloat->bone()+4; # now use Pari
2889
2890While this works, it loads Calc needlessly. But maybe you just wanted that?
2891
2892B<Examples #3 is highly recommended> for daily usage.
2893
2894=head1 BUGS
2895
2896Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
58cde26e 2897
990fb837 2898=head1 CAVEATS
58cde26e 2899
2900=over 1
2901
2902=item stringify, bstr()
2903
2904Both stringify and bstr() now drop the leading '+'. The old code would return
2905'+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2906reasoning and details.
2907
2908=item bdiv
2909
2910The following will probably not do what you expect:
2911
2912 print $c->bdiv(123.456),"\n";
2913
2914It prints both quotient and reminder since print works in list context. Also,
2915bdiv() will modify $c, so be carefull. You probably want to use
2916
2917 print $c / 123.456,"\n";
2918 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
2919
2920instead.
2921
2922=item Modifying and =
2923
2924Beware of:
2925
2926 $x = Math::BigFloat->new(5);
2927 $y = $x;
2928
2929It will not do what you think, e.g. making a copy of $x. Instead it just makes
2930a second reference to the B<same> object and stores it in $y. Thus anything
990fb837 2931that modifies $x will modify $y (except overloaded math operators), and vice
2932versa. See L<Math::BigInt> for details and how to avoid that.
58cde26e 2933
2934=item bpow
2935
2936C<bpow()> now modifies the first argument, unlike the old code which left
2937it alone and only returned the result. This is to be consistent with
2938C<badd()> etc. The first will modify $x, the second one won't:
2939
2940 print bpow($x,$i),"\n"; # modify $x
2941 print $x->bpow($i),"\n"; # ditto
2942 print $x ** $i,"\n"; # leave $x alone
2943
2944=back
2945
990fb837 2946=head1 SEE ALSO
2947
2948L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
2949L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
2950
2951The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
2952because they solve the autoupgrading/downgrading issue, at least partly.
2953
2954The package at
2955L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
2956more documentation including a full version history, testcases, empty
2957subclass files and benchmarks.
2958
58cde26e 2959=head1 LICENSE
a5f75d66 2960
58cde26e 2961This program is free software; you may redistribute it and/or modify it under
2962the same terms as Perl itself.
5d7098d5 2963
58cde26e 2964=head1 AUTHORS
5d7098d5 2965
58cde26e 2966Mark Biggar, overloaded interface by Ilya Zakharevich.
990fb837 2967Completely rewritten by Tels http://bloodgate.com in 2001, 2002, and still
2968at it in 2003.
a5f75d66 2969
a5f75d66 2970=cut