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