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