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 used internally: 8# sign : "+", "-", "+inf", "-inf", or "NaN" if not a number 9# _m : mantissa ($LIB thingy) 10# _es : sign of _e 11# _e : exponent ($LIB thingy) 12# _a : accuracy 13# _p : precision 14 15use 5.006001; 16use strict; 17use warnings; 18 19use Carp qw< carp croak >; 20use Math::BigInt (); 21 22our $VERSION = '1.999816'; 23 24our @ISA = qw/Math::BigInt/; 25our @EXPORT_OK = qw/bpi/; 26 27# $_trap_inf/$_trap_nan are internal and should never be accessed from outside 28our ($AUTOLOAD, $accuracy, $precision, $div_scale, $round_mode, $rnd_mode, 29 $upgrade, $downgrade, $_trap_nan, $_trap_inf); 30 31my $class = "Math::BigFloat"; 32 33use overload 34 35 # overload key: with_assign 36 37 '+' => sub { $_[0] -> copy() -> badd($_[1]); }, 38 39 '-' => sub { my $c = $_[0] -> copy(); 40 $_[2] ? $c -> bneg() -> badd($_[1]) 41 : $c -> bsub($_[1]); }, 42 43 '*' => sub { $_[0] -> copy() -> bmul($_[1]); }, 44 45 '/' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bdiv($_[0]) 46 : $_[0] -> copy() -> bdiv($_[1]); }, 47 48 '%' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bmod($_[0]) 49 : $_[0] -> copy() -> bmod($_[1]); }, 50 51 '**' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bpow($_[0]) 52 : $_[0] -> copy() -> bpow($_[1]); }, 53 54 '<<' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> blsft($_[0]) 55 : $_[0] -> copy() -> blsft($_[1]); }, 56 57 '>>' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> brsft($_[0]) 58 : $_[0] -> copy() -> brsft($_[1]); }, 59 60 # overload key: assign 61 62 '+=' => sub { $_[0]->badd($_[1]); }, 63 64 '-=' => sub { $_[0]->bsub($_[1]); }, 65 66 '*=' => sub { $_[0]->bmul($_[1]); }, 67 68 '/=' => sub { scalar $_[0]->bdiv($_[1]); }, 69 70 '%=' => sub { $_[0]->bmod($_[1]); }, 71 72 '**=' => sub { $_[0]->bpow($_[1]); }, 73 74 '<<=' => sub { $_[0]->blsft($_[1]); }, 75 76 '>>=' => sub { $_[0]->brsft($_[1]); }, 77 78# 'x=' => sub { }, 79 80# '.=' => sub { }, 81 82 # overload key: num_comparison 83 84 '<' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> blt($_[0]) 85 : $_[0] -> blt($_[1]); }, 86 87 '<=' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> ble($_[0]) 88 : $_[0] -> ble($_[1]); }, 89 90 '>' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bgt($_[0]) 91 : $_[0] -> bgt($_[1]); }, 92 93 '>=' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bge($_[0]) 94 : $_[0] -> bge($_[1]); }, 95 96 '==' => sub { $_[0] -> beq($_[1]); }, 97 98 '!=' => sub { $_[0] -> bne($_[1]); }, 99 100 # overload key: 3way_comparison 101 102 '<=>' => sub { my $cmp = $_[0] -> bcmp($_[1]); 103 defined($cmp) && $_[2] ? -$cmp : $cmp; }, 104 105 'cmp' => sub { $_[2] ? "$_[1]" cmp $_[0] -> bstr() 106 : $_[0] -> bstr() cmp "$_[1]"; }, 107 108 # overload key: str_comparison 109 110# 'lt' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrlt($_[0]) 111# : $_[0] -> bstrlt($_[1]); }, 112# 113# 'le' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrle($_[0]) 114# : $_[0] -> bstrle($_[1]); }, 115# 116# 'gt' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrgt($_[0]) 117# : $_[0] -> bstrgt($_[1]); }, 118# 119# 'ge' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrge($_[0]) 120# : $_[0] -> bstrge($_[1]); }, 121# 122# 'eq' => sub { $_[0] -> bstreq($_[1]); }, 123# 124# 'ne' => sub { $_[0] -> bstrne($_[1]); }, 125 126 # overload key: binary 127 128 '&' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> band($_[0]) 129 : $_[0] -> copy() -> band($_[1]); }, 130 131 '&=' => sub { $_[0] -> band($_[1]); }, 132 133 '|' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bior($_[0]) 134 : $_[0] -> copy() -> bior($_[1]); }, 135 136 '|=' => sub { $_[0] -> bior($_[1]); }, 137 138 '^' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bxor($_[0]) 139 : $_[0] -> copy() -> bxor($_[1]); }, 140 141 '^=' => sub { $_[0] -> bxor($_[1]); }, 142 143# '&.' => sub { }, 144 145# '&.=' => sub { }, 146 147# '|.' => sub { }, 148 149# '|.=' => sub { }, 150 151# '^.' => sub { }, 152 153# '^.=' => sub { }, 154 155 # overload key: unary 156 157 'neg' => sub { $_[0] -> copy() -> bneg(); }, 158 159# '!' => sub { }, 160 161 '~' => sub { $_[0] -> copy() -> bnot(); }, 162 163# '~.' => sub { }, 164 165 # overload key: mutators 166 167 '++' => sub { $_[0] -> binc() }, 168 169 '--' => sub { $_[0] -> bdec() }, 170 171 # overload key: func 172 173 'atan2' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> batan2($_[0]) 174 : $_[0] -> copy() -> batan2($_[1]); }, 175 176 'cos' => sub { $_[0] -> copy() -> bcos(); }, 177 178 'sin' => sub { $_[0] -> copy() -> bsin(); }, 179 180 'exp' => sub { $_[0] -> copy() -> bexp($_[1]); }, 181 182 'abs' => sub { $_[0] -> copy() -> babs(); }, 183 184 'log' => sub { $_[0] -> copy() -> blog(); }, 185 186 'sqrt' => sub { $_[0] -> copy() -> bsqrt(); }, 187 188 'int' => sub { $_[0] -> copy() -> bint(); }, 189 190 # overload key: conversion 191 192 'bool' => sub { $_[0] -> is_zero() ? '' : 1; }, 193 194 '""' => sub { $_[0] -> bstr(); }, 195 196 '0+' => sub { $_[0] -> numify(); }, 197 198 '=' => sub { $_[0]->copy(); }, 199 200 ; 201 202############################################################################## 203# global constants, flags and assorted stuff 204 205# the following are public, but their usage is not recommended. Use the 206# accessor methods instead. 207 208# class constants, use Class->constant_name() to access 209# one of 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common' 210$round_mode = 'even'; 211$accuracy = undef; 212$precision = undef; 213$div_scale = 40; 214 215$upgrade = undef; 216$downgrade = undef; 217# the package we are using for our private parts, defaults to: 218# Math::BigInt->config('lib') 219my $LIB = 'Math::BigInt::Calc'; 220 221# are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config() 222$_trap_nan = 0; 223# the same for infinity 224$_trap_inf = 0; 225 226# constant for easier life 227my $nan = 'NaN'; 228 229my $IMPORT = 0; # was import() called yet? used to make require work 230 231# some digits of accuracy for blog(undef, 10); which we use in blog() for speed 232my $LOG_10 = 233 '2.3025850929940456840179914546843642076011014886287729760333279009675726097'; 234my $LOG_10_A = length($LOG_10)-1; 235# ditto for log(2) 236my $LOG_2 = 237 '0.6931471805599453094172321214581765680755001343602552541206800094933936220'; 238my $LOG_2_A = length($LOG_2)-1; 239my $HALF = '0.5'; # made into an object if nec. 240 241############################################################################## 242# the old code had $rnd_mode, so we need to support it, too 243 244sub TIESCALAR { 245 my ($class) = @_; 246 bless \$round_mode, $class; 247} 248 249sub FETCH { 250 return $round_mode; 251} 252 253sub STORE { 254 $rnd_mode = $_[0]->round_mode($_[1]); 255} 256 257BEGIN { 258 # when someone sets $rnd_mode, we catch this and check the value to see 259 # whether it is valid or not. 260 $rnd_mode = 'even'; 261 tie $rnd_mode, 'Math::BigFloat'; 262 263 # we need both of them in this package: 264 *as_int = \&as_number; 265} 266 267sub DESTROY { 268 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub 269} 270 271sub AUTOLOAD { 272 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx() 273 my $name = $AUTOLOAD; 274 275 $name =~ s/(.*):://; # split package 276 my $c = $1 || $class; 277 no strict 'refs'; 278 $c->import() if $IMPORT == 0; 279 if (!_method_alias($name)) { 280 if (!defined $name) { 281 # delayed load of Carp and avoid recursion 282 croak("$c: Can't call a method without name"); 283 } 284 if (!_method_hand_up($name)) { 285 # delayed load of Carp and avoid recursion 286 croak("Can't call $c\-\>$name, not a valid method"); 287 } 288 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx() 289 $name =~ s/^f/b/; 290 return &{"Math::BigInt"."::$name"}(@_); 291 } 292 my $bname = $name; 293 $bname =~ s/^f/b/; 294 $c .= "::$name"; 295 *{$c} = \&{$bname}; 296 &{$c}; # uses @_ 297} 298 299############################################################################## 300 301{ 302 # valid method aliases for AUTOLOAD 303 my %methods = map { $_ => 1 } 304 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm 305 fint facmp fcmp fzero fnan finf finc fdec ffac fneg 306 fceil ffloor frsft flsft fone flog froot fexp 307 /; 308 # valid methods that can be handed up (for AUTOLOAD) 309 my %hand_ups = map { $_ => 1 } 310 qw / is_nan is_inf is_negative is_positive is_pos is_neg 311 accuracy precision div_scale round_mode fabs fnot 312 objectify upgrade downgrade 313 bone binf bnan bzero 314 bsub 315 /; 316 317 sub _method_alias { exists $methods{$_[0]||''}; } 318 sub _method_hand_up { exists $hand_ups{$_[0]||''}; } 319} 320 321sub isa { 322 my ($self, $class) = @_; 323 return if $class =~ /^Math::BigInt/; # we aren't one of these 324 UNIVERSAL::isa($self, $class); 325} 326 327sub config { 328 # return (later set?) configuration data as hash ref 329 my $class = shift || 'Math::BigFloat'; 330 331 # Getter/accessor. 332 333 if (@_ == 1 && ref($_[0]) ne 'HASH') { 334 my $param = shift; 335 return $class if $param eq 'class'; 336 return $LIB if $param eq 'with'; 337 return $class->SUPER::config($param); 338 } 339 340 # Setter. 341 342 my $cfg = $class->SUPER::config(@_); 343 344 # now we need only to override the ones that are different from our parent 345 $cfg->{class} = $class; 346 $cfg->{with} = $LIB; 347 $cfg; 348} 349 350############################################################################### 351# Constructor methods 352############################################################################### 353 354sub new { 355 # Create a new Math::BigFloat object from a string or another bigfloat object. 356 # _e: exponent 357 # _m: mantissa 358 # sign => ("+", "-", "+inf", "-inf", or "NaN") 359 360 my $self = shift; 361 my $selfref = ref $self; 362 my $class = $selfref || $self; 363 364 my ($wanted, @r) = @_; 365 366 # avoid numify-calls by not using || on $wanted! 367 368 unless (defined $wanted) { 369 #carp("Use of uninitialized value in new"); 370 return $self->bzero(@r); 371 } 372 373 # Using $wanted->isa("Math::BigFloat") here causes a 'Deep recursion on 374 # subroutine "Math::BigFloat::as_number"' in some tests. Fixme! 375 376 if (UNIVERSAL::isa($wanted, 'Math::BigFloat')) { 377 my $copy = $wanted -> copy(); 378 if ($selfref) { # if new() called as instance method 379 %$self = %$copy; 380 } else { # if new() called as class method 381 $self = $copy; 382 } 383 return $copy; 384 } 385 386 $class->import() if $IMPORT == 0; # make require work 387 388 # If called as a class method, initialize a new object. 389 390 $self = bless {}, $class unless $selfref; 391 392 # shortcut for bigints and its subclasses 393 if ((ref($wanted)) && $wanted -> can("as_number")) { 394 $self->{_m} = $wanted->as_number()->{value}; # get us a bigint copy 395 $self->{_e} = $LIB->_zero(); 396 $self->{_es} = '+'; 397 $self->{sign} = $wanted->sign(); 398 return $self->bnorm(); 399 } 400 401 # else: got a string or something masquerading as number (with overload) 402 403 # Handle Infs. 404 405 if ($wanted =~ /^\s*([+-]?)inf(inity)?\s*\z/i) { 406 return $downgrade->new($wanted) if $downgrade; 407 my $sgn = $1 || '+'; 408 $self->{sign} = $sgn . 'inf'; # set a default sign for bstr() 409 return $self->binf($sgn); 410 } 411 412 # Handle explicit NaNs (not the ones returned due to invalid input). 413 414 if ($wanted =~ /^\s*([+-]?)nan\s*\z/i) { 415 return $downgrade->new($wanted) if $downgrade; 416 $self = $class -> bnan(); 417 $self->round(@r) unless @r >= 2 && !defined $r[0] && !defined $r[1]; 418 return $self; 419 } 420 421 # Handle hexadecimal numbers. 422 423 if ($wanted =~ /^\s*[+-]?0[Xx]/) { 424 $self = $class -> from_hex($wanted); 425 $self->round(@r) unless @r >= 2 && !defined $r[0] && !defined $r[1]; 426 return $self; 427 } 428 429 # Handle binary numbers. 430 431 if ($wanted =~ /^\s*[+-]?0[Bb]/) { 432 $self = $class -> from_bin($wanted); 433 $self->round(@r) unless @r >= 2 && !defined $r[0] && !defined $r[1]; 434 return $self; 435 } 436 437 # Shortcut for simple forms like '12' that have no trailing zeros. 438 if ($wanted =~ /^([+-]?)0*([1-9][0-9]*[1-9])$/) { 439 $self->{_e} = $LIB -> _zero(); 440 $self->{_es} = '+'; 441 $self->{sign} = $1 || '+'; 442 $self->{_m} = $LIB -> _new($2); 443 if (!$downgrade) { 444 $self->round(@r) unless @r >= 2 && !defined $r[0] && !defined $r[1]; 445 return $self; 446 } 447 } 448 449 my ($mis, $miv, $mfv, $es, $ev) = Math::BigInt::_split($wanted); 450 if (!ref $mis) { 451 if ($_trap_nan) { 452 croak("$wanted is not a number initialized to $class"); 453 } 454 455 return $downgrade->bnan() if $downgrade; 456 457 $self->{_e} = $LIB->_zero(); 458 $self->{_es} = '+'; 459 $self->{_m} = $LIB->_zero(); 460 $self->{sign} = $nan; 461 } else { 462 # make integer from mantissa by adjusting exp, then convert to int 463 $self->{_e} = $LIB->_new($$ev); # exponent 464 $self->{_es} = $$es || '+'; 465 my $mantissa = "$$miv$$mfv"; # create mant. 466 $mantissa =~ s/^0+(\d)/$1/; # strip leading zeros 467 $self->{_m} = $LIB->_new($mantissa); # create mant. 468 469 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5 470 if (CORE::length($$mfv) != 0) { 471 my $len = $LIB->_new(CORE::length($$mfv)); 472 ($self->{_e}, $self->{_es}) = 473 _e_sub($self->{_e}, $len, $self->{_es}, '+'); 474 } 475 # we can only have trailing zeros on the mantissa if $$mfv eq '' 476 else { 477 # Use a regexp to count the trailing zeros in $$miv instead of 478 # _zeros() because that is faster, especially when _m is not stored 479 # in base 10. 480 my $zeros = 0; 481 $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/; 482 if ($zeros != 0) { 483 my $z = $LIB->_new($zeros); 484 # turn '120e2' into '12e3' 485 $self->{_m} = $LIB->_rsft($self->{_m}, $z, 10); 486 ($self->{_e}, $self->{_es}) = 487 _e_add($self->{_e}, $z, $self->{_es}, '+'); 488 } 489 } 490 $self->{sign} = $$mis; 491 492 # for something like 0Ey, set y to 0, and -0 => +0 493 # Check $$miv for being '0' and $$mfv eq '', because otherwise _m could not 494 # have become 0. That's faster than to call $LIB->_is_zero(). 495 $self->{sign} = '+', $self->{_e} = $LIB->_zero() 496 if $$miv eq '0' and $$mfv eq ''; 497 498 if (!$downgrade) { 499 $self->round(@r) unless @r >= 2 && !defined $r[0] && !defined $r[1]; 500 return $self; 501 } 502 } 503 504 # if downgrade, inf, NaN or integers go down 505 506 if ($downgrade && $self->{_es} eq '+') { 507 if ($LIB->_is_zero($self->{_e})) { 508 return $downgrade->new($$mis . $LIB->_str($self->{_m})); 509 } 510 return $downgrade->new($self->bsstr()); 511 } 512 $self->bnorm(); 513 $self->round(@r) unless @r >= 2 && !defined $r[0] && !defined $r[1]; 514 return $self; 515} 516 517sub from_hex { 518 my $self = shift; 519 my $selfref = ref $self; 520 my $class = $selfref || $self; 521 522 # Don't modify constant (read-only) objects. 523 524 return if $selfref && $self->modify('from_hex'); 525 526 my $str = shift; 527 528 # If called as a class method, initialize a new object. 529 530 $self = $class -> bzero() unless $selfref; 531 532 if ($str =~ s/ 533 ^ 534 \s* 535 536 # sign 537 ( [+-]? ) 538 539 # optional "hex marker" 540 (?: 0? x )? 541 542 # significand using the hex digits 0..9 and a..f 543 ( 544 [0-9a-fA-F]+ (?: _ [0-9a-fA-F]+ )* 545 (?: 546 \. 547 (?: [0-9a-fA-F]+ (?: _ [0-9a-fA-F]+ )* )? 548 )? 549 | 550 \. 551 [0-9a-fA-F]+ (?: _ [0-9a-fA-F]+ )* 552 ) 553 554 # exponent (power of 2) using decimal digits 555 (?: 556 [Pp] 557 ( [+-]? ) 558 ( \d+ (?: _ \d+ )* ) 559 )? 560 561 \s* 562 $ 563 //x) 564 { 565 my $s_sign = $1 || '+'; 566 my $s_value = $2; 567 my $e_sign = $3 || '+'; 568 my $e_value = $4 || '0'; 569 $s_value =~ tr/_//d; 570 $e_value =~ tr/_//d; 571 572 # The significand must be multiplied by 2 raised to this exponent. 573 574 my $two_expon = $class -> new($e_value); 575 $two_expon -> bneg() if $e_sign eq '-'; 576 577 # If there is a dot in the significand, remove it and adjust the 578 # exponent according to the number of digits in the fraction part of 579 # the significand. Since the digits in the significand are in base 16, 580 # but the exponent is only in base 2, multiply the exponent adjustment 581 # value by log(16) / log(2) = 4. 582 583 my $idx = index($s_value, '.'); 584 if ($idx >= 0) { 585 substr($s_value, $idx, 1) = ''; 586 $two_expon -= $class -> new(CORE::length($s_value)) 587 -> bsub($idx) 588 -> bmul("4"); 589 } 590 591 $self -> {sign} = $s_sign; 592 $self -> {_m} = $LIB -> _from_hex('0x' . $s_value); 593 594 if ($two_expon > 0) { 595 my $factor = $class -> new("2") -> bpow($two_expon); 596 $self -> bmul($factor); 597 } elsif ($two_expon < 0) { 598 my $factor = $class -> new("0.5") -> bpow(-$two_expon); 599 $self -> bmul($factor); 600 } 601 602 return $self; 603 } 604 605 return $self->bnan(); 606} 607 608sub from_oct { 609 my $self = shift; 610 my $selfref = ref $self; 611 my $class = $selfref || $self; 612 613 # Don't modify constant (read-only) objects. 614 615 return if $selfref && $self->modify('from_oct'); 616 617 my $str = shift; 618 619 # If called as a class method, initialize a new object. 620 621 $self = $class -> bzero() unless $selfref; 622 623 if ($str =~ s/ 624 ^ 625 \s* 626 627 # sign 628 ( [+-]? ) 629 630 # significand using the octal digits 0..7 631 ( 632 [0-7]+ (?: _ [0-7]+ )* 633 (?: 634 \. 635 (?: [0-7]+ (?: _ [0-7]+ )* )? 636 )? 637 | 638 \. 639 [0-7]+ (?: _ [0-7]+ )* 640 ) 641 642 # exponent (power of 2) using decimal digits 643 (?: 644 [Pp] 645 ( [+-]? ) 646 ( \d+ (?: _ \d+ )* ) 647 )? 648 649 \s* 650 $ 651 //x) 652 { 653 my $s_sign = $1 || '+'; 654 my $s_value = $2; 655 my $e_sign = $3 || '+'; 656 my $e_value = $4 || '0'; 657 $s_value =~ tr/_//d; 658 $e_value =~ tr/_//d; 659 660 # The significand must be multiplied by 2 raised to this exponent. 661 662 my $two_expon = $class -> new($e_value); 663 $two_expon -> bneg() if $e_sign eq '-'; 664 665 # If there is a dot in the significand, remove it and adjust the 666 # exponent according to the number of digits in the fraction part of 667 # the significand. Since the digits in the significand are in base 8, 668 # but the exponent is only in base 2, multiply the exponent adjustment 669 # value by log(8) / log(2) = 3. 670 671 my $idx = index($s_value, '.'); 672 if ($idx >= 0) { 673 substr($s_value, $idx, 1) = ''; 674 $two_expon -= $class -> new(CORE::length($s_value)) 675 -> bsub($idx) 676 -> bmul("3"); 677 } 678 679 $self -> {sign} = $s_sign; 680 $self -> {_m} = $LIB -> _from_oct($s_value); 681 682 if ($two_expon > 0) { 683 my $factor = $class -> new("2") -> bpow($two_expon); 684 $self -> bmul($factor); 685 } elsif ($two_expon < 0) { 686 my $factor = $class -> new("0.5") -> bpow(-$two_expon); 687 $self -> bmul($factor); 688 } 689 690 return $self; 691 } 692 693 return $self->bnan(); 694} 695 696sub from_bin { 697 my $self = shift; 698 my $selfref = ref $self; 699 my $class = $selfref || $self; 700 701 # Don't modify constant (read-only) objects. 702 703 return if $selfref && $self->modify('from_bin'); 704 705 my $str = shift; 706 707 # If called as a class method, initialize a new object. 708 709 $self = $class -> bzero() unless $selfref; 710 711 if ($str =~ s/ 712 ^ 713 \s* 714 715 # sign 716 ( [+-]? ) 717 718 # optional "bin marker" 719 (?: 0? b )? 720 721 # significand using the binary digits 0 and 1 722 ( 723 [01]+ (?: _ [01]+ )* 724 (?: 725 \. 726 (?: [01]+ (?: _ [01]+ )* )? 727 )? 728 | 729 \. 730 [01]+ (?: _ [01]+ )* 731 ) 732 733 # exponent (power of 2) using decimal digits 734 (?: 735 [Pp] 736 ( [+-]? ) 737 ( \d+ (?: _ \d+ )* ) 738 )? 739 740 \s* 741 $ 742 //x) 743 { 744 my $s_sign = $1 || '+'; 745 my $s_value = $2; 746 my $e_sign = $3 || '+'; 747 my $e_value = $4 || '0'; 748 $s_value =~ tr/_//d; 749 $e_value =~ tr/_//d; 750 751 # The significand must be multiplied by 2 raised to this exponent. 752 753 my $two_expon = $class -> new($e_value); 754 $two_expon -> bneg() if $e_sign eq '-'; 755 756 # If there is a dot in the significand, remove it and adjust the 757 # exponent according to the number of digits in the fraction part of 758 # the significand. 759 760 my $idx = index($s_value, '.'); 761 if ($idx >= 0) { 762 substr($s_value, $idx, 1) = ''; 763 $two_expon -= $class -> new(CORE::length($s_value)) 764 -> bsub($idx); 765 } 766 767 $self -> {sign} = $s_sign; 768 $self -> {_m} = $LIB -> _from_bin('0b' . $s_value); 769 770 if ($two_expon > 0) { 771 my $factor = $class -> new("2") -> bpow($two_expon); 772 $self -> bmul($factor); 773 } elsif ($two_expon < 0) { 774 my $factor = $class -> new("0.5") -> bpow(-$two_expon); 775 $self -> bmul($factor); 776 } 777 778 return $self; 779 } 780 781 return $self->bnan(); 782} 783 784sub bzero { 785 # create/assign '+0' 786 787 if (@_ == 0) { 788 #carp("Using bone() as a function is deprecated;", 789 # " use bone() as a method instead"); 790 unshift @_, __PACKAGE__; 791 } 792 793 my $self = shift; 794 my $selfref = ref $self; 795 my $class = $selfref || $self; 796 797 $self->import() if $IMPORT == 0; # make require work 798 return if $selfref && $self->modify('bzero'); 799 800 $self = bless {}, $class unless $selfref; 801 802 $self -> {sign} = '+'; 803 $self -> {_m} = $LIB -> _zero(); 804 $self -> {_es} = '+'; 805 $self -> {_e} = $LIB -> _zero(); 806 807 # If rounding parameters are given as arguments, use them. If no rounding 808 # parameters are given, and if called as a class method initialize the new 809 # instance with the class variables. 810 811 if (@_) { 812 croak "can't specify both accuracy and precision" 813 if @_ >= 2 && defined $_[0] && defined $_[1]; 814 $self->{_a} = $_[0]; 815 $self->{_p} = $_[1]; 816 } else { 817 unless($selfref) { 818 $self->{_a} = $class -> accuracy(); 819 $self->{_p} = $class -> precision(); 820 } 821 } 822 823 return $self; 824} 825 826sub bone { 827 # Create or assign '+1' (or -1 if given sign '-'). 828 829 if (@_ == 0 || (defined($_[0]) && ($_[0] eq '+' || $_[0] eq '-'))) { 830 #carp("Using bone() as a function is deprecated;", 831 # " use bone() as a method instead"); 832 unshift @_, __PACKAGE__; 833 } 834 835 my $self = shift; 836 my $selfref = ref $self; 837 my $class = $selfref || $self; 838 839 $self->import() if $IMPORT == 0; # make require work 840 return if $selfref && $self->modify('bone'); 841 842 my $sign = shift; 843 $sign = defined $sign && $sign =~ /^\s*-/ ? "-" : "+"; 844 845 $self = bless {}, $class unless $selfref; 846 847 $self -> {sign} = $sign; 848 $self -> {_m} = $LIB -> _one(); 849 $self -> {_es} = '+'; 850 $self -> {_e} = $LIB -> _zero(); 851 852 # If rounding parameters are given as arguments, use them. If no rounding 853 # parameters are given, and if called as a class method initialize the new 854 # instance with the class variables. 855 856 if (@_) { 857 croak "can't specify both accuracy and precision" 858 if @_ >= 2 && defined $_[0] && defined $_[1]; 859 $self->{_a} = $_[0]; 860 $self->{_p} = $_[1]; 861 } else { 862 unless($selfref) { 863 $self->{_a} = $class -> accuracy(); 864 $self->{_p} = $class -> precision(); 865 } 866 } 867 868 return $self; 869} 870 871sub binf { 872 # create/assign a '+inf' or '-inf' 873 874 if (@_ == 0 || (defined($_[0]) && !ref($_[0]) && 875 $_[0] =~ /^\s*[+-](inf(inity)?)?\s*$/)) 876 { 877 #carp("Using binf() as a function is deprecated;", 878 # " use binf() as a method instead"); 879 unshift @_, __PACKAGE__; 880 } 881 882 my $self = shift; 883 my $selfref = ref $self; 884 my $class = $selfref || $self; 885 886 { 887 no strict 'refs'; 888 if (${"${class}::_trap_inf"}) { 889 croak("Tried to create +-inf in $class->binf()"); 890 } 891 } 892 893 $self->import() if $IMPORT == 0; # make require work 894 return if $selfref && $self->modify('binf'); 895 896 my $sign = shift; 897 $sign = defined $sign && $sign =~ /^\s*-/ ? "-" : "+"; 898 899 $self = bless {}, $class unless $selfref; 900 901 $self -> {sign} = $sign . 'inf'; 902 $self -> {_m} = $LIB -> _zero(); 903 $self -> {_es} = '+'; 904 $self -> {_e} = $LIB -> _zero(); 905 906 # If rounding parameters are given as arguments, use them. If no rounding 907 # parameters are given, and if called as a class method initialize the new 908 # instance with the class variables. 909 910 if (@_) { 911 croak "can't specify both accuracy and precision" 912 if @_ >= 2 && defined $_[0] && defined $_[1]; 913 $self->{_a} = $_[0]; 914 $self->{_p} = $_[1]; 915 } else { 916 unless($selfref) { 917 $self->{_a} = $class -> accuracy(); 918 $self->{_p} = $class -> precision(); 919 } 920 } 921 922 return $self; 923} 924 925sub bnan { 926 # create/assign a 'NaN' 927 928 if (@_ == 0) { 929 #carp("Using bnan() as a function is deprecated;", 930 # " use bnan() as a method instead"); 931 unshift @_, __PACKAGE__; 932 } 933 934 my $self = shift; 935 my $selfref = ref $self; 936 my $class = $selfref || $self; 937 938 { 939 no strict 'refs'; 940 if (${"${class}::_trap_nan"}) { 941 croak("Tried to create NaN in $class->bnan()"); 942 } 943 } 944 945 $self->import() if $IMPORT == 0; # make require work 946 return if $selfref && $self->modify('bnan'); 947 948 $self = bless {}, $class unless $selfref; 949 950 $self -> {sign} = $nan; 951 $self -> {_m} = $LIB -> _zero(); 952 $self -> {_es} = '+'; 953 $self -> {_e} = $LIB -> _zero(); 954 955 # If rounding parameters are given as arguments, use them. If no rounding 956 # parameters are given, and if called as a class method initialize the new 957 # instance with the class variables. 958 959 if (@_) { 960 croak "can't specify both accuracy and precision" 961 if @_ >= 2 && defined $_[0] && defined $_[1]; 962 $self->{_a} = $_[0]; 963 $self->{_p} = $_[1]; 964 } else { 965 unless($selfref) { 966 $self->{_a} = $class -> accuracy(); 967 $self->{_p} = $class -> precision(); 968 } 969 } 970 971 return $self; 972} 973 974sub bpi { 975 976 # Called as Argument list 977 # --------- ------------- 978 # Math::BigFloat->bpi() ("Math::BigFloat") 979 # Math::BigFloat->bpi(10) ("Math::BigFloat", 10) 980 # $x->bpi() ($x) 981 # $x->bpi(10) ($x, 10) 982 # Math::BigFloat::bpi() () 983 # Math::BigFloat::bpi(10) (10) 984 # 985 # In ambiguous cases, we favour the OO-style, so the following case 986 # 987 # $n = Math::BigFloat->new("10"); 988 # $x = Math::BigFloat->bpi($n); 989 # 990 # which gives an argument list with the single element $n, is resolved as 991 # 992 # $n->bpi(); 993 994 my $self = shift; 995 my $selfref = ref $self; 996 my $class = $selfref || $self; 997 998 my @r; # rounding paramters 999 1000 # If bpi() is called as a function ... 1001 # 1002 # This cludge is necessary because we still support bpi() as a function. If 1003 # bpi() is called with either no argument or one argument, and that one 1004 # argument is either undefined or a scalar that looks like a number, then 1005 # we assume bpi() is called as a function. 1006 1007 if (@_ == 0 && 1008 (defined($self) && !ref($self) && $self =~ /^\s*[+-]?\d/i) 1009 || 1010 !defined($self)) 1011 { 1012 $r[0] = $self; 1013 $class = __PACKAGE__; 1014 $self = $class -> bzero(@r); # initialize 1015 } 1016 1017 # ... or if bpi() is called as a method ... 1018 1019 else { 1020 @r = @_; 1021 if ($selfref) { # bpi() called as instance method 1022 return $self if $self -> modify('bpi'); 1023 } else { # bpi() called as class method 1024 $self = $class -> bzero(@r); # initialize 1025 } 1026 } 1027 1028 ($self, @r) = $self -> _find_round_parameters(@r); 1029 1030 # The accuracy, i.e., the number of digits. Pi has one digit before the 1031 # dot, so a precision of 4 digits is equivalent to an accuracy of 5 digits. 1032 1033 my $n = defined $r[0] ? $r[0] 1034 : defined $r[1] ? 1 - $r[1] 1035 : $self -> div_scale(); 1036 1037 my $rmode = defined $r[2] ? $r[2] : $self -> round_mode(); 1038 1039 my $pi; 1040 1041 if ($n <= 1000) { 1042 1043 # 75 x 14 = 1050 digits 1044 1045 my $all_digits = <<EOF; 1046314159265358979323846264338327950288419716939937510582097494459230781640628 1047620899862803482534211706798214808651328230664709384460955058223172535940812 1048848111745028410270193852110555964462294895493038196442881097566593344612847 1049564823378678316527120190914564856692346034861045432664821339360726024914127 1050372458700660631558817488152092096282925409171536436789259036001133053054882 1051046652138414695194151160943305727036575959195309218611738193261179310511854 1052807446237996274956735188575272489122793818301194912983367336244065664308602 1053139494639522473719070217986094370277053921717629317675238467481846766940513 1054200056812714526356082778577134275778960917363717872146844090122495343014654 1055958537105079227968925892354201995611212902196086403441815981362977477130996 1056051870721134999999837297804995105973173281609631859502445945534690830264252 1057230825334468503526193118817101000313783875288658753320838142061717766914730 1058359825349042875546873115956286388235378759375195778185778053217122680661300 1059192787661119590921642019893809525720106548586327886593615338182796823030195 1060EOF 1061 1062 # Should we round up? 1063 1064 my $round_up; 1065 1066 # From the string above, we need to extract the number of digits we 1067 # want plus extra characters for the newlines. 1068 1069 my $nchrs = $n + int($n / 75); 1070 1071 # Extract the digits we want. 1072 1073 my $digits = substr($all_digits, 0, $nchrs); 1074 1075 # Find out whether we should round up or down. Rounding is easy, since 1076 # pi is trancendental. With directed rounding, it doesn't matter what 1077 # the following digits are. With rounding to nearest, we only have to 1078 # look at one extra digit. 1079 1080 if ($rmode eq 'trunc') { 1081 $round_up = 0; 1082 } else { 1083 my $next_digit = substr($all_digits, $nchrs, 1); 1084 $round_up = $next_digit lt '5' ? 0 : 1; 1085 } 1086 1087 # Remove the newlines. 1088 1089 $digits =~ tr/0-9//cd; 1090 1091 # Now do the rounding. We could easily make the regex substitution 1092 # handle all cases, but we avoid using the regex engine when it is 1093 # simple to avoid it. 1094 1095 if ($round_up) { 1096 my $last_digit = substr($digits, -1, 1); 1097 if ($last_digit lt '9') { 1098 substr($digits, -1, 1) = ++$last_digit; 1099 } else { 1100 $digits =~ s/([0-8])(9+)$/ ($1 + 1) . ("0" x CORE::length($2)) /e; 1101 } 1102 } 1103 1104 # Append the exponent and convert to an object. 1105 1106 $pi = Math::BigFloat -> new($digits . 'e-' . ($n - 1)); 1107 1108 } else { 1109 1110 # For large accuracy, the arctan formulas become very inefficient with 1111 # Math::BigFloat, so use Brent-Salamin (aka AGM or Gauss-Legendre). 1112 1113 # Use a few more digits in the intermediate computations. 1114 $n += 8; 1115 1116 $HALF = $class -> new($HALF) unless ref($HALF); 1117 my ($an, $bn, $tn, $pn) = ($class -> bone, $HALF -> copy() -> bsqrt($n), 1118 $HALF -> copy() -> bmul($HALF), $class -> bone); 1119 while ($pn < $n) { 1120 my $prev_an = $an -> copy(); 1121 $an -> badd($bn) -> bmul($HALF, $n); 1122 $bn -> bmul($prev_an) -> bsqrt($n); 1123 $prev_an -> bsub($an); 1124 $tn -> bsub($pn * $prev_an * $prev_an); 1125 $pn -> badd($pn); 1126 } 1127 $an -> badd($bn); 1128 $an -> bmul($an, $n) -> bdiv(4 * $tn, $n); 1129 1130 $an -> round(@r); 1131 $pi = $an; 1132 } 1133 1134 if (defined $r[0]) { 1135 $pi -> accuracy($r[0]); 1136 } elsif (defined $r[1]) { 1137 $pi -> precision($r[1]); 1138 } 1139 1140 for my $key (qw/ sign _m _es _e _a _p /) { 1141 $self -> {$key} = $pi -> {$key}; 1142 } 1143 1144 return $self; 1145} 1146 1147sub copy { 1148 my $self = shift; 1149 my $selfref = ref $self; 1150 my $class = $selfref || $self; 1151 1152 # If called as a class method, the object to copy is the next argument. 1153 1154 $self = shift() unless $selfref; 1155 1156 my $copy = bless {}, $class; 1157 1158 $copy->{sign} = $self->{sign}; 1159 $copy->{_es} = $self->{_es}; 1160 $copy->{_m} = $LIB->_copy($self->{_m}); 1161 $copy->{_e} = $LIB->_copy($self->{_e}); 1162 $copy->{_a} = $self->{_a} if exists $self->{_a}; 1163 $copy->{_p} = $self->{_p} if exists $self->{_p}; 1164 1165 return $copy; 1166} 1167 1168sub as_number { 1169 # return copy as a bigint representation of this Math::BigFloat number 1170 my ($class, $x) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 1171 1172 return $x if $x->modify('as_number'); 1173 1174 if (!$x->isa('Math::BigFloat')) { 1175 # if the object can as_number(), use it 1176 return $x->as_number() if $x->can('as_number'); 1177 # otherwise, get us a float and then a number 1178 $x = $x->can('as_float') ? $x->as_float() : $class->new(0+"$x"); 1179 } 1180 1181 return Math::BigInt->binf($x->sign()) if $x->is_inf(); 1182 return Math::BigInt->bnan() if $x->is_nan(); 1183 1184 my $z = $LIB->_copy($x->{_m}); 1185 if ($x->{_es} eq '-') { # < 0 1186 $z = $LIB->_rsft($z, $x->{_e}, 10); 1187 } elsif (! $LIB->_is_zero($x->{_e})) { # > 0 1188 $z = $LIB->_lsft($z, $x->{_e}, 10); 1189 } 1190 $z = Math::BigInt->new($x->{sign} . $LIB->_str($z)); 1191 $z; 1192} 1193 1194############################################################################### 1195# Boolean methods 1196############################################################################### 1197 1198sub is_zero { 1199 # return true if arg (BFLOAT or num_str) is zero 1200 my ($class, $x) = ref($_[0]) ? (undef, $_[0]) : objectify(1, @_); 1201 1202 ($x->{sign} eq '+' && $LIB->_is_zero($x->{_m})) ? 1 : 0; 1203} 1204 1205sub is_one { 1206 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given 1207 my ($class, $x, $sign) = ref($_[0]) ? (undef, @_) : objectify(1, @_); 1208 1209 $sign = '+' if !defined $sign || $sign ne '-'; 1210 1211 ($x->{sign} eq $sign && 1212 $LIB->_is_zero($x->{_e}) && 1213 $LIB->_is_one($x->{_m})) ? 1 : 0; 1214} 1215 1216sub is_odd { 1217 # return true if arg (BFLOAT or num_str) is odd or false if even 1218 my ($class, $x) = ref($_[0]) ? (undef, $_[0]) : objectify(1, @_); 1219 1220 (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't 1221 ($LIB->_is_zero($x->{_e})) && 1222 ($LIB->_is_odd($x->{_m}))) ? 1 : 0; 1223} 1224 1225sub is_even { 1226 # return true if arg (BINT or num_str) is even or false if odd 1227 my ($class, $x) = ref($_[0]) ? (undef, $_[0]) : objectify(1, @_); 1228 1229 (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't 1230 ($x->{_es} eq '+') && # 123.45 isn't 1231 ($LIB->_is_even($x->{_m}))) ? 1 : 0; # but 1200 is 1232} 1233 1234sub is_int { 1235 # return true if arg (BFLOAT or num_str) is an integer 1236 my ($class, $x) = ref($_[0]) ? (undef, $_[0]) : objectify(1, @_); 1237 1238 (($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't 1239 ($x->{_es} eq '+')) ? 1 : 0; # 1e-1 => no integer 1240} 1241 1242############################################################################### 1243# Comparison methods 1244############################################################################### 1245 1246sub bcmp { 1247 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort) 1248 1249 # set up parameters 1250 my ($class, $x, $y) = (ref($_[0]), @_); 1251 1252 # objectify is costly, so avoid it 1253 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 1254 ($class, $x, $y) = objectify(2, @_); 1255 } 1256 1257 return $upgrade->bcmp($x, $y) if defined $upgrade && 1258 ((!$x->isa($class)) || (!$y->isa($class))); 1259 1260 # Handle all 'nan' cases. 1261 1262 return undef if ($x->{sign} eq $nan) || ($y->{sign} eq $nan); 1263 1264 # Handle all '+inf' and '-inf' cases. 1265 1266 return 0 if ($x->{sign} eq '+inf' && $y->{sign} eq '+inf' || 1267 $x->{sign} eq '-inf' && $y->{sign} eq '-inf'); 1268 return +1 if $x->{sign} eq '+inf'; # x = +inf and y < +inf 1269 return -1 if $x->{sign} eq '-inf'; # x = -inf and y > -inf 1270 return -1 if $y->{sign} eq '+inf'; # x < +inf and y = +inf 1271 return +1 if $y->{sign} eq '-inf'; # x > -inf and y = -inf 1272 1273 # Handle all cases with opposite signs. 1274 1275 return +1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # also does 0 <=> -y 1276 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # also does -x <=> 0 1277 1278 # Handle all remaining zero cases. 1279 1280 my $xz = $x->is_zero(); 1281 my $yz = $y->is_zero(); 1282 return 0 if $xz && $yz; # 0 <=> 0 1283 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y 1284 return +1 if $yz && $x->{sign} eq '+'; # +x <=> 0 1285 1286 # Both arguments are now finite, non-zero numbers with the same sign. 1287 1288 my $cmp; 1289 1290 # The next step is to compare the exponents, but since each mantissa is an 1291 # integer of arbitrary value, the exponents must be normalized by the length 1292 # of the mantissas before we can compare them. 1293 1294 my $mxl = $LIB->_len($x->{_m}); 1295 my $myl = $LIB->_len($y->{_m}); 1296 1297 # If the mantissas have the same length, there is no point in normalizing the 1298 # exponents by the length of the mantissas, so treat that as a special case. 1299 1300 if ($mxl == $myl) { 1301 1302 # First handle the two cases where the exponents have different signs. 1303 1304 if ($x->{_es} eq '+' && $y->{_es} eq '-') { 1305 $cmp = +1; 1306 } elsif ($x->{_es} eq '-' && $y->{_es} eq '+') { 1307 $cmp = -1; 1308 } 1309 1310 # Then handle the case where the exponents have the same sign. 1311 1312 else { 1313 $cmp = $LIB->_acmp($x->{_e}, $y->{_e}); 1314 $cmp = -$cmp if $x->{_es} eq '-'; 1315 } 1316 1317 # Adjust for the sign, which is the same for x and y, and bail out if 1318 # we're done. 1319 1320 $cmp = -$cmp if $x->{sign} eq '-'; # 124 > 123, but -124 < -123 1321 return $cmp if $cmp; 1322 1323 } 1324 1325 # We must normalize each exponent by the length of the corresponding 1326 # mantissa. Life is a lot easier if we first make both exponents 1327 # non-negative. We do this by adding the same positive value to both 1328 # exponent. This is safe, because when comparing the exponents, only the 1329 # relative difference is important. 1330 1331 my $ex; 1332 my $ey; 1333 1334 if ($x->{_es} eq '+') { 1335 1336 # If the exponent of x is >= 0 and the exponent of y is >= 0, there is no 1337 # need to do anything special. 1338 1339 if ($y->{_es} eq '+') { 1340 $ex = $LIB->_copy($x->{_e}); 1341 $ey = $LIB->_copy($y->{_e}); 1342 } 1343 1344 # If the exponent of x is >= 0 and the exponent of y is < 0, add the 1345 # absolute value of the exponent of y to both. 1346 1347 else { 1348 $ex = $LIB->_copy($x->{_e}); 1349 $ex = $LIB->_add($ex, $y->{_e}); # ex + |ey| 1350 $ey = $LIB->_zero(); # -ex + |ey| = 0 1351 } 1352 1353 } else { 1354 1355 # If the exponent of x is < 0 and the exponent of y is >= 0, add the 1356 # absolute value of the exponent of x to both. 1357 1358 if ($y->{_es} eq '+') { 1359 $ex = $LIB->_zero(); # -ex + |ex| = 0 1360 $ey = $LIB->_copy($y->{_e}); 1361 $ey = $LIB->_add($ey, $x->{_e}); # ey + |ex| 1362 } 1363 1364 # If the exponent of x is < 0 and the exponent of y is < 0, add the 1365 # absolute values of both exponents to both exponents. 1366 1367 else { 1368 $ex = $LIB->_copy($y->{_e}); # -ex + |ey| + |ex| = |ey| 1369 $ey = $LIB->_copy($x->{_e}); # -ey + |ex| + |ey| = |ex| 1370 } 1371 1372 } 1373 1374 # Now we can normalize the exponents by adding lengths of the mantissas. 1375 1376 $ex = $LIB->_add($ex, $LIB->_new($mxl)); 1377 $ey = $LIB->_add($ey, $LIB->_new($myl)); 1378 1379 # We're done if the exponents are different. 1380 1381 $cmp = $LIB->_acmp($ex, $ey); 1382 $cmp = -$cmp if $x->{sign} eq '-'; # 124 > 123, but -124 < -123 1383 return $cmp if $cmp; 1384 1385 # Compare the mantissas, but first normalize them by padding the shorter 1386 # mantissa with zeros (shift left) until it has the same length as the longer 1387 # mantissa. 1388 1389 my $mx = $x->{_m}; 1390 my $my = $y->{_m}; 1391 1392 if ($mxl > $myl) { 1393 $my = $LIB->_lsft($LIB->_copy($my), $LIB->_new($mxl - $myl), 10); 1394 } elsif ($mxl < $myl) { 1395 $mx = $LIB->_lsft($LIB->_copy($mx), $LIB->_new($myl - $mxl), 10); 1396 } 1397 1398 $cmp = $LIB->_acmp($mx, $my); 1399 $cmp = -$cmp if $x->{sign} eq '-'; # 124 > 123, but -124 < -123 1400 return $cmp; 1401 1402} 1403 1404sub bacmp { 1405 # Compares 2 values, ignoring their signs. 1406 # Returns one of undef, <0, =0, >0. (suitable for sort) 1407 1408 # set up parameters 1409 my ($class, $x, $y) = (ref($_[0]), @_); 1410 # objectify is costly, so avoid it 1411 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 1412 ($class, $x, $y) = objectify(2, @_); 1413 } 1414 1415 return $upgrade->bacmp($x, $y) if defined $upgrade && 1416 ((!$x->isa($class)) || (!$y->isa($class))); 1417 1418 # handle +-inf and NaN's 1419 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/) { 1420 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan)); 1421 return 0 if ($x->is_inf() && $y->is_inf()); 1422 return 1 if ($x->is_inf() && !$y->is_inf()); 1423 return -1; 1424 } 1425 1426 # shortcut 1427 my $xz = $x->is_zero(); 1428 my $yz = $y->is_zero(); 1429 return 0 if $xz && $yz; # 0 <=> 0 1430 return -1 if $xz && !$yz; # 0 <=> +y 1431 return 1 if $yz && !$xz; # +x <=> 0 1432 1433 # adjust so that exponents are equal 1434 my $lxm = $LIB->_len($x->{_m}); 1435 my $lym = $LIB->_len($y->{_m}); 1436 my ($xes, $yes) = (1, 1); 1437 $xes = -1 if $x->{_es} ne '+'; 1438 $yes = -1 if $y->{_es} ne '+'; 1439 # the numify somewhat limits our length, but makes it much faster 1440 my $lx = $lxm + $xes * $LIB->_num($x->{_e}); 1441 my $ly = $lym + $yes * $LIB->_num($y->{_e}); 1442 my $l = $lx - $ly; 1443 return $l <=> 0 if $l != 0; 1444 1445 # lengths (corrected by exponent) are equal 1446 # so make mantissa equal-length by padding with zero (shift left) 1447 my $diff = $lxm - $lym; 1448 my $xm = $x->{_m}; # not yet copy it 1449 my $ym = $y->{_m}; 1450 if ($diff > 0) { 1451 $ym = $LIB->_copy($y->{_m}); 1452 $ym = $LIB->_lsft($ym, $LIB->_new($diff), 10); 1453 } elsif ($diff < 0) { 1454 $xm = $LIB->_copy($x->{_m}); 1455 $xm = $LIB->_lsft($xm, $LIB->_new(-$diff), 10); 1456 } 1457 $LIB->_acmp($xm, $ym); 1458} 1459 1460############################################################################### 1461# Arithmetic methods 1462############################################################################### 1463 1464sub bneg { 1465 # (BINT or num_str) return BINT 1466 # negate number or make a negated number from string 1467 my ($class, $x) = ref($_[0]) ? (undef, $_[0]) : objectify(1, @_); 1468 1469 return $x if $x->modify('bneg'); 1470 1471 # for +0 do not negate (to have always normalized +0). Does nothing for 'NaN' 1472 $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $LIB->_is_zero($x->{_m})); 1473 $x; 1474} 1475 1476sub bnorm { 1477 # adjust m and e so that m is smallest possible 1478 my ($class, $x) = ref($_[0]) ? (undef, $_[0]) : objectify(1, @_); 1479 1480 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc 1481 1482 my $zeros = $LIB->_zeros($x->{_m}); # correct for trailing zeros 1483 if ($zeros != 0) { 1484 my $z = $LIB->_new($zeros); 1485 $x->{_m} = $LIB->_rsft($x->{_m}, $z, 10); 1486 if ($x->{_es} eq '-') { 1487 if ($LIB->_acmp($x->{_e}, $z) >= 0) { 1488 $x->{_e} = $LIB->_sub($x->{_e}, $z); 1489 $x->{_es} = '+' if $LIB->_is_zero($x->{_e}); 1490 } else { 1491 $x->{_e} = $LIB->_sub($LIB->_copy($z), $x->{_e}); 1492 $x->{_es} = '+'; 1493 } 1494 } else { 1495 $x->{_e} = $LIB->_add($x->{_e}, $z); 1496 } 1497 } else { 1498 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing 1499 # zeros). So, for something like 0Ey, set y to 1, and -0 => +0 1500 $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $LIB->_one() 1501 if $LIB->_is_zero($x->{_m}); 1502 } 1503 1504 $x; 1505} 1506 1507sub binc { 1508 # increment arg by one 1509 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 1510 1511 return $x if $x->modify('binc'); 1512 1513 if ($x->{_es} eq '-') { 1514 return $x->badd($class->bone(), @r); # digits after dot 1515 } 1516 1517 if (!$LIB->_is_zero($x->{_e})) # _e == 0 for NaN, inf, -inf 1518 { 1519 # 1e2 => 100, so after the shift below _m has a '0' as last digit 1520 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # 1e2 => 100 1521 $x->{_e} = $LIB->_zero(); # normalize 1522 $x->{_es} = '+'; 1523 # we know that the last digit of $x will be '1' or '9', depending on the 1524 # sign 1525 } 1526 # now $x->{_e} == 0 1527 if ($x->{sign} eq '+') { 1528 $x->{_m} = $LIB->_inc($x->{_m}); 1529 return $x->bnorm()->bround(@r); 1530 } elsif ($x->{sign} eq '-') { 1531 $x->{_m} = $LIB->_dec($x->{_m}); 1532 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # -1 +1 => -0 => +0 1533 return $x->bnorm()->bround(@r); 1534 } 1535 # inf, nan handling etc 1536 $x->badd($class->bone(), @r); # badd() does round 1537} 1538 1539sub bdec { 1540 # decrement arg by one 1541 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 1542 1543 return $x if $x->modify('bdec'); 1544 1545 if ($x->{_es} eq '-') { 1546 return $x->badd($class->bone('-'), @r); # digits after dot 1547 } 1548 1549 if (!$LIB->_is_zero($x->{_e})) { 1550 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # 1e2 => 100 1551 $x->{_e} = $LIB->_zero(); # normalize 1552 $x->{_es} = '+'; 1553 } 1554 # now $x->{_e} == 0 1555 my $zero = $x->is_zero(); 1556 # <= 0 1557 if (($x->{sign} eq '-') || $zero) { 1558 $x->{_m} = $LIB->_inc($x->{_m}); 1559 $x->{sign} = '-' if $zero; # 0 => 1 => -1 1560 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # -1 +1 => -0 => +0 1561 return $x->bnorm()->round(@r); 1562 } 1563 # > 0 1564 elsif ($x->{sign} eq '+') { 1565 $x->{_m} = $LIB->_dec($x->{_m}); 1566 return $x->bnorm()->round(@r); 1567 } 1568 # inf, nan handling etc 1569 $x->badd($class->bone('-'), @r); # does round 1570} 1571 1572sub badd { 1573 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first) 1574 # return result as BFLOAT 1575 1576 # set up parameters 1577 my ($class, $x, $y, @r) = (ref($_[0]), @_); 1578 # objectify is costly, so avoid it 1579 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 1580 ($class, $x, $y, @r) = objectify(2, @_); 1581 } 1582 1583 return $x if $x->modify('badd'); 1584 1585 # inf and NaN handling 1586 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/)) { 1587 # NaN first 1588 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan)); 1589 # inf handling 1590 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/)) { 1591 # +inf++inf or -inf+-inf => same, rest is NaN 1592 return $x if $x->{sign} eq $y->{sign}; 1593 return $x->bnan(); 1594 } 1595 # +-inf + something => +inf; something +-inf => +-inf 1596 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/; 1597 return $x; 1598 } 1599 1600 return $upgrade->badd($x, $y, @r) if defined $upgrade && 1601 ((!$x->isa($class)) || (!$y->isa($class))); 1602 1603 $r[3] = $y; # no push! 1604 1605 # speed: no add for 0+y or x+0 1606 return $x->bround(@r) if $y->is_zero(); # x+0 1607 if ($x->is_zero()) # 0+y 1608 { 1609 # make copy, clobbering up x (modify in place!) 1610 $x->{_e} = $LIB->_copy($y->{_e}); 1611 $x->{_es} = $y->{_es}; 1612 $x->{_m} = $LIB->_copy($y->{_m}); 1613 $x->{sign} = $y->{sign} || $nan; 1614 return $x->round(@r); 1615 } 1616 1617 # take lower of the two e's and adapt m1 to it to match m2 1618 my $e = $y->{_e}; 1619 $e = $LIB->_zero() if !defined $e; # if no BFLOAT? 1620 $e = $LIB->_copy($e); # make copy (didn't do it yet) 1621 1622 my $es; 1623 1624 ($e, $es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es}); 1625 1626 my $add = $LIB->_copy($y->{_m}); 1627 1628 if ($es eq '-') # < 0 1629 { 1630 $x->{_m} = $LIB->_lsft($x->{_m}, $e, 10); 1631 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es); 1632 } elsif (!$LIB->_is_zero($e)) # > 0 1633 { 1634 $add = $LIB->_lsft($add, $e, 10); 1635 } 1636 # else: both e are the same, so just leave them 1637 1638 if ($x->{sign} eq $y->{sign}) { 1639 # add 1640 $x->{_m} = $LIB->_add($x->{_m}, $add); 1641 } else { 1642 ($x->{_m}, $x->{sign}) = 1643 _e_add($x->{_m}, $add, $x->{sign}, $y->{sign}); 1644 } 1645 1646 # delete trailing zeros, then round 1647 $x->bnorm()->round(@r); 1648} 1649 1650sub bsub { 1651 # (BINT or num_str, BINT or num_str) return BINT 1652 # subtract second arg from first, modify first 1653 1654 # set up parameters 1655 my ($class, $x, $y, @r) = (ref($_[0]), @_); 1656 1657 # objectify is costly, so avoid it 1658 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 1659 ($class, $x, $y, @r) = objectify(2, @_); 1660 } 1661 1662 return $x if $x -> modify('bsub'); 1663 1664 return $upgrade -> new($x) -> bsub($upgrade -> new($y), @r) 1665 if defined $upgrade && (!$x -> isa($class) || !$y -> isa($class)); 1666 1667 return $x -> round(@r) if $y -> is_zero(); 1668 1669 # To correctly handle the lone special case $x -> bsub($x), we note the 1670 # sign of $x, then flip the sign from $y, and if the sign of $x did change, 1671 # too, then we caught the special case: 1672 1673 my $xsign = $x -> {sign}; 1674 $y -> {sign} =~ tr/+-/-+/; # does nothing for NaN 1675 if ($xsign ne $x -> {sign}) { 1676 # special case of $x -> bsub($x) results in 0 1677 return $x -> bzero(@r) if $xsign =~ /^[+-]$/; 1678 return $x -> bnan(); # NaN, -inf, +inf 1679 } 1680 $x -> badd($y, @r); # badd does not leave internal zeros 1681 $y -> {sign} =~ tr/+-/-+/; # refix $y (does nothing for NaN) 1682 $x; # already rounded by badd() or no rounding 1683} 1684 1685sub bmul { 1686 # multiply two numbers 1687 1688 # set up parameters 1689 my ($class, $x, $y, @r) = (ref($_[0]), @_); 1690 # objectify is costly, so avoid it 1691 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 1692 ($class, $x, $y, @r) = objectify(2, @_); 1693 } 1694 1695 return $x if $x->modify('bmul'); 1696 1697 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan)); 1698 1699 # inf handling 1700 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/)) { 1701 return $x->bnan() if $x->is_zero() || $y->is_zero(); 1702 # result will always be +-inf: 1703 # +inf * +/+inf => +inf, -inf * -/-inf => +inf 1704 # +inf * -/-inf => -inf, -inf * +/+inf => -inf 1705 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/); 1706 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/); 1707 return $x->binf('-'); 1708 } 1709 1710 return $upgrade->bmul($x, $y, @r) if defined $upgrade && 1711 ((!$x->isa($class)) || (!$y->isa($class))); 1712 1713 # aEb * cEd = (a*c)E(b+d) 1714 $x->{_m} = $LIB->_mul($x->{_m}, $y->{_m}); 1715 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es}); 1716 1717 $r[3] = $y; # no push! 1718 1719 # adjust sign: 1720 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+'; 1721 $x->bnorm->round(@r); 1722} 1723 1724sub bmuladd { 1725 # multiply two numbers and add the third to the result 1726 1727 # set up parameters 1728 my ($class, $x, $y, $z, @r) = objectify(3, @_); 1729 1730 return $x if $x->modify('bmuladd'); 1731 1732 return $x->bnan() if (($x->{sign} eq $nan) || 1733 ($y->{sign} eq $nan) || 1734 ($z->{sign} eq $nan)); 1735 1736 # inf handling 1737 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/)) { 1738 return $x->bnan() if $x->is_zero() || $y->is_zero(); 1739 # result will always be +-inf: 1740 # +inf * +/+inf => +inf, -inf * -/-inf => +inf 1741 # +inf * -/-inf => -inf, -inf * +/+inf => -inf 1742 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/); 1743 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/); 1744 return $x->binf('-'); 1745 } 1746 1747 return $upgrade->bmul($x, $y, @r) if defined $upgrade && 1748 ((!$x->isa($class)) || (!$y->isa($class))); 1749 1750 # aEb * cEd = (a*c)E(b+d) 1751 $x->{_m} = $LIB->_mul($x->{_m}, $y->{_m}); 1752 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es}); 1753 1754 $r[3] = $y; # no push! 1755 1756 # adjust sign: 1757 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+'; 1758 1759 # z=inf handling (z=NaN handled above) 1760 $x->{sign} = $z->{sign}, return $x if $z->{sign} =~ /^[+-]inf$/; 1761 1762 # take lower of the two e's and adapt m1 to it to match m2 1763 my $e = $z->{_e}; 1764 $e = $LIB->_zero() if !defined $e; # if no BFLOAT? 1765 $e = $LIB->_copy($e); # make copy (didn't do it yet) 1766 1767 my $es; 1768 1769 ($e, $es) = _e_sub($e, $x->{_e}, $z->{_es} || '+', $x->{_es}); 1770 1771 my $add = $LIB->_copy($z->{_m}); 1772 1773 if ($es eq '-') # < 0 1774 { 1775 $x->{_m} = $LIB->_lsft($x->{_m}, $e, 10); 1776 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es); 1777 } elsif (!$LIB->_is_zero($e)) # > 0 1778 { 1779 $add = $LIB->_lsft($add, $e, 10); 1780 } 1781 # else: both e are the same, so just leave them 1782 1783 if ($x->{sign} eq $z->{sign}) { 1784 # add 1785 $x->{_m} = $LIB->_add($x->{_m}, $add); 1786 } else { 1787 ($x->{_m}, $x->{sign}) = 1788 _e_add($x->{_m}, $add, $x->{sign}, $z->{sign}); 1789 } 1790 1791 # delete trailing zeros, then round 1792 $x->bnorm()->round(@r); 1793} 1794 1795sub bdiv { 1796 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return 1797 # (BFLOAT, BFLOAT) (quo, rem) or BFLOAT (only quo) 1798 1799 # set up parameters 1800 my ($class, $x, $y, $a, $p, $r) = (ref($_[0]), @_); 1801 # objectify is costly, so avoid it 1802 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 1803 ($class, $x, $y, $a, $p, $r) = objectify(2, @_); 1804 } 1805 1806 return $x if $x->modify('bdiv'); 1807 1808 my $wantarray = wantarray; # call only once 1809 1810 # At least one argument is NaN. This is handled the same way as in 1811 # Math::BigInt -> bdiv(). 1812 1813 if ($x -> is_nan() || $y -> is_nan()) { 1814 return $wantarray ? ($x -> bnan(), $class -> bnan()) : $x -> bnan(); 1815 } 1816 1817 # Divide by zero and modulo zero. This is handled the same way as in 1818 # Math::BigInt -> bdiv(). See the comment in the code for Math::BigInt -> 1819 # bdiv() for further details. 1820 1821 if ($y -> is_zero()) { 1822 my ($quo, $rem); 1823 if ($wantarray) { 1824 $rem = $x -> copy(); 1825 } 1826 if ($x -> is_zero()) { 1827 $quo = $x -> bnan(); 1828 } else { 1829 $quo = $x -> binf($x -> {sign}); 1830 } 1831 return $wantarray ? ($quo, $rem) : $quo; 1832 } 1833 1834 # Numerator (dividend) is +/-inf. This is handled the same way as in 1835 # Math::BigInt -> bdiv(). See the comment in the code for Math::BigInt -> 1836 # bdiv() for further details. 1837 1838 if ($x -> is_inf()) { 1839 my ($quo, $rem); 1840 $rem = $class -> bnan() if $wantarray; 1841 if ($y -> is_inf()) { 1842 $quo = $x -> bnan(); 1843 } else { 1844 my $sign = $x -> bcmp(0) == $y -> bcmp(0) ? '+' : '-'; 1845 $quo = $x -> binf($sign); 1846 } 1847 return $wantarray ? ($quo, $rem) : $quo; 1848 } 1849 1850 # Denominator (divisor) is +/-inf. This is handled the same way as in 1851 # Math::BigInt -> bdiv(), with one exception: In scalar context, 1852 # Math::BigFloat does true division (although rounded), not floored division 1853 # (F-division), so a finite number divided by +/-inf is always zero. See the 1854 # comment in the code for Math::BigInt -> bdiv() for further details. 1855 1856 if ($y -> is_inf()) { 1857 my ($quo, $rem); 1858 if ($wantarray) { 1859 if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) { 1860 $rem = $x -> copy(); 1861 $quo = $x -> bzero(); 1862 } else { 1863 $rem = $class -> binf($y -> {sign}); 1864 $quo = $x -> bone('-'); 1865 } 1866 return ($quo, $rem); 1867 } else { 1868 if ($y -> is_inf()) { 1869 if ($x -> is_nan() || $x -> is_inf()) { 1870 return $x -> bnan(); 1871 } else { 1872 return $x -> bzero(); 1873 } 1874 } 1875 } 1876 } 1877 1878 # At this point, both the numerator and denominator are finite numbers, and 1879 # the denominator (divisor) is non-zero. 1880 1881 # x == 0? 1882 return wantarray ? ($x, $class->bzero()) : $x if $x->is_zero(); 1883 1884 # upgrade ? 1885 return $upgrade->bdiv($upgrade->new($x), $y, $a, $p, $r) if defined $upgrade; 1886 1887 # we need to limit the accuracy to protect against overflow 1888 my $fallback = 0; 1889 my (@params, $scale); 1890 ($x, @params) = $x->_find_round_parameters($a, $p, $r, $y); 1891 1892 return $x if $x->is_nan(); # error in _find_round_parameters? 1893 1894 # no rounding at all, so must use fallback 1895 if (scalar @params == 0) { 1896 # simulate old behaviour 1897 $params[0] = $class->div_scale(); # and round to it as accuracy 1898 $scale = $params[0]+4; # at least four more for proper round 1899 $params[2] = $r; # round mode by caller or undef 1900 $fallback = 1; # to clear a/p afterwards 1901 } else { 1902 # the 4 below is empirical, and there might be cases where it is not 1903 # enough... 1904 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 1905 } 1906 1907 my $rem; 1908 $rem = $class -> bzero() if wantarray; 1909 1910 $y = $class->new($y) unless $y->isa('Math::BigFloat'); 1911 1912 my $lx = $LIB -> _len($x->{_m}); my $ly = $LIB -> _len($y->{_m}); 1913 $scale = $lx if $lx > $scale; 1914 $scale = $ly if $ly > $scale; 1915 my $diff = $ly - $lx; 1916 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx! 1917 1918 # check that $y is not 1 nor -1 and cache the result: 1919 my $y_not_one = !($LIB->_is_zero($y->{_e}) && $LIB->_is_one($y->{_m})); 1920 1921 # flipping the sign of $y will also flip the sign of $x for the special 1922 # case of $x->bsub($x); so we can catch it below: 1923 my $xsign = $x->{sign}; 1924 $y->{sign} =~ tr/+-/-+/; 1925 1926 if ($xsign ne $x->{sign}) { 1927 # special case of $x /= $x results in 1 1928 $x->bone(); # "fixes" also sign of $y, since $x is $y 1929 } else { 1930 # correct $y's sign again 1931 $y->{sign} =~ tr/+-/-+/; 1932 # continue with normal div code: 1933 1934 # make copy of $x in case of list context for later remainder calculation 1935 if (wantarray && $y_not_one) { 1936 $rem = $x->copy(); 1937 } 1938 1939 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+'; 1940 1941 # check for / +-1 (+/- 1E0) 1942 if ($y_not_one) { 1943 # promote BigInts and it's subclasses (except when already a Math::BigFloat) 1944 $y = $class->new($y) unless $y->isa('Math::BigFloat'); 1945 1946 # calculate the result to $scale digits and then round it 1947 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d) 1948 $x->{_m} = $LIB->_lsft($x->{_m}, $LIB->_new($scale), 10); 1949 $x->{_m} = $LIB->_div($x->{_m}, $y->{_m}); # a/c 1950 1951 # correct exponent of $x 1952 ($x->{_e}, $x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es}); 1953 # correct for 10**scale 1954 ($x->{_e}, $x->{_es}) = _e_sub($x->{_e}, $LIB->_new($scale), $x->{_es}, '+'); 1955 $x->bnorm(); # remove trailing 0's 1956 } 1957 } # end else $x != $y 1958 1959 # shortcut to not run through _find_round_parameters again 1960 if (defined $params[0]) { 1961 delete $x->{_a}; # clear before round 1962 $x->bround($params[0], $params[2]); # then round accordingly 1963 } else { 1964 delete $x->{_p}; # clear before round 1965 $x->bfround($params[1], $params[2]); # then round accordingly 1966 } 1967 if ($fallback) { 1968 # clear a/p after round, since user did not request it 1969 delete $x->{_a}; delete $x->{_p}; 1970 } 1971 1972 if (wantarray) { 1973 if ($y_not_one) { 1974 $x -> bfloor(); 1975 $rem->bmod($y, @params); # copy already done 1976 } 1977 if ($fallback) { 1978 # clear a/p after round, since user did not request it 1979 delete $rem->{_a}; delete $rem->{_p}; 1980 } 1981 return ($x, $rem); 1982 } 1983 $x; 1984} 1985 1986sub bmod { 1987 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return remainder 1988 1989 # set up parameters 1990 my ($class, $x, $y, $a, $p, $r) = (ref($_[0]), @_); 1991 # objectify is costly, so avoid it 1992 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 1993 ($class, $x, $y, $a, $p, $r) = objectify(2, @_); 1994 } 1995 1996 return $x if $x->modify('bmod'); 1997 1998 # At least one argument is NaN. This is handled the same way as in 1999 # Math::BigInt -> bmod(). 2000 2001 if ($x -> is_nan() || $y -> is_nan()) { 2002 return $x -> bnan(); 2003 } 2004 2005 # Modulo zero. This is handled the same way as in Math::BigInt -> bmod(). 2006 2007 if ($y -> is_zero()) { 2008 return $x; 2009 } 2010 2011 # Numerator (dividend) is +/-inf. This is handled the same way as in 2012 # Math::BigInt -> bmod(). 2013 2014 if ($x -> is_inf()) { 2015 return $x -> bnan(); 2016 } 2017 2018 # Denominator (divisor) is +/-inf. This is handled the same way as in 2019 # Math::BigInt -> bmod(). 2020 2021 if ($y -> is_inf()) { 2022 if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) { 2023 return $x; 2024 } else { 2025 return $x -> binf($y -> sign()); 2026 } 2027 } 2028 2029 return $x->bzero() if $x->is_zero() 2030 || ($x->is_int() && 2031 # check that $y == +1 or $y == -1: 2032 ($LIB->_is_zero($y->{_e}) && $LIB->_is_one($y->{_m}))); 2033 2034 my $cmp = $x->bacmp($y); # equal or $x < $y? 2035 if ($cmp == 0) { # $x == $y => result 0 2036 return $x -> bzero($a, $p); 2037 } 2038 2039 # only $y of the operands negative? 2040 my $neg = $x->{sign} ne $y->{sign} ? 1 : 0; 2041 2042 $x->{sign} = $y->{sign}; # calc sign first 2043 if ($cmp < 0 && $neg == 0) { # $x < $y => result $x 2044 return $x -> round($a, $p, $r); 2045 } 2046 2047 my $ym = $LIB->_copy($y->{_m}); 2048 2049 # 2e1 => 20 2050 $ym = $LIB->_lsft($ym, $y->{_e}, 10) 2051 if $y->{_es} eq '+' && !$LIB->_is_zero($y->{_e}); 2052 2053 # if $y has digits after dot 2054 my $shifty = 0; # correct _e of $x by this 2055 if ($y->{_es} eq '-') # has digits after dot 2056 { 2057 # 123 % 2.5 => 1230 % 25 => 5 => 0.5 2058 $shifty = $LIB->_num($y->{_e}); # no more digits after dot 2059 $x->{_m} = $LIB->_lsft($x->{_m}, $y->{_e}, 10); # 123 => 1230, $y->{_m} is already 25 2060 } 2061 # $ym is now mantissa of $y based on exponent 0 2062 2063 my $shiftx = 0; # correct _e of $x by this 2064 if ($x->{_es} eq '-') # has digits after dot 2065 { 2066 # 123.4 % 20 => 1234 % 200 2067 $shiftx = $LIB->_num($x->{_e}); # no more digits after dot 2068 $ym = $LIB->_lsft($ym, $x->{_e}, 10); # 123 => 1230 2069 } 2070 # 123e1 % 20 => 1230 % 20 2071 if ($x->{_es} eq '+' && !$LIB->_is_zero($x->{_e})) { 2072 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # es => '+' here 2073 } 2074 2075 $x->{_e} = $LIB->_new($shiftx); 2076 $x->{_es} = '+'; 2077 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0; 2078 $x->{_e} = $LIB->_add($x->{_e}, $LIB->_new($shifty)) if $shifty != 0; 2079 2080 # now mantissas are equalized, exponent of $x is adjusted, so calc result 2081 2082 $x->{_m} = $LIB->_mod($x->{_m}, $ym); 2083 2084 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # fix sign for -0 2085 $x->bnorm(); 2086 2087 if ($neg != 0 && ! $x -> is_zero()) # one of them negative => correct in place 2088 { 2089 my $r = $y - $x; 2090 $x->{_m} = $r->{_m}; 2091 $x->{_e} = $r->{_e}; 2092 $x->{_es} = $r->{_es}; 2093 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # fix sign for -0 2094 $x->bnorm(); 2095 } 2096 2097 $x->round($a, $p, $r, $y); # round and return 2098} 2099 2100sub bmodpow { 2101 # takes a very large number to a very large exponent in a given very 2102 # large modulus, quickly, thanks to binary exponentiation. Supports 2103 # negative exponents. 2104 my ($class, $num, $exp, $mod) = objectify(3, @_); 2105 2106 return $num if $num->modify('bmodpow'); 2107 2108 # check modulus for valid values 2109 return $num->bnan() if ($mod->{sign} ne '+' # NaN, -, -inf, +inf 2110 || $mod->is_zero()); 2111 2112 # check exponent for valid values 2113 if ($exp->{sign} =~ /\w/) { 2114 # i.e., if it's NaN, +inf, or -inf... 2115 return $num->bnan(); 2116 } 2117 2118 $num->bmodinv ($mod) if ($exp->{sign} eq '-'); 2119 2120 # check num for valid values (also NaN if there was no inverse but $exp < 0) 2121 return $num->bnan() if $num->{sign} !~ /^[+-]$/; 2122 2123 # $mod is positive, sign on $exp is ignored, result also positive 2124 2125 # XXX TODO: speed it up when all three numbers are integers 2126 $num->bpow($exp)->bmod($mod); 2127} 2128 2129sub bpow { 2130 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT 2131 # compute power of two numbers, second arg is used as integer 2132 # modifies first argument 2133 2134 # set up parameters 2135 my ($class, $x, $y, $a, $p, $r) = (ref($_[0]), @_); 2136 # objectify is costly, so avoid it 2137 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 2138 ($class, $x, $y, $a, $p, $r) = objectify(2, @_); 2139 } 2140 2141 return $x if $x->modify('bpow'); 2142 2143 # $x and/or $y is a NaN 2144 return $x->bnan() if $x->is_nan() || $y->is_nan(); 2145 2146 # $x and/or $y is a +/-Inf 2147 if ($x->is_inf("-")) { 2148 return $x->bzero() if $y->is_negative(); 2149 return $x->bnan() if $y->is_zero(); 2150 return $x if $y->is_odd(); 2151 return $x->bneg(); 2152 } elsif ($x->is_inf("+")) { 2153 return $x->bzero() if $y->is_negative(); 2154 return $x->bnan() if $y->is_zero(); 2155 return $x; 2156 } elsif ($y->is_inf("-")) { 2157 return $x->bnan() if $x -> is_one("-"); 2158 return $x->binf("+") if $x > -1 && $x < 1; 2159 return $x->bone() if $x -> is_one("+"); 2160 return $x->bzero(); 2161 } elsif ($y->is_inf("+")) { 2162 return $x->bnan() if $x -> is_one("-"); 2163 return $x->bzero() if $x > -1 && $x < 1; 2164 return $x->bone() if $x -> is_one("+"); 2165 return $x->binf("+"); 2166 } 2167 2168 # we don't support complex numbers, so return NaN 2169 return $x->bnan() if $x->is_negative() && !$y->is_int(); 2170 2171 # cache the result of is_zero 2172 my $y_is_zero = $y->is_zero(); 2173 return $x->bone() if $y_is_zero; 2174 return $x if $x->is_one() || $y->is_one(); 2175 2176 my $x_is_zero = $x->is_zero(); 2177 return $x->_pow($y, $a, $p, $r) if !$x_is_zero && !$y->is_int(); 2178 2179 my $y1 = $y->as_number()->{value}; # make MBI part 2180 2181 if ($x->is_one("-")) { 2182 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1 2183 return $LIB->_is_odd($y1) ? $x : $x->babs(1); 2184 } 2185 if ($x_is_zero) { 2186 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0) 2187 # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf) 2188 return $x->binf(); 2189 } 2190 2191 my $new_sign = '+'; 2192 $new_sign = $LIB->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+'; 2193 2194 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster) 2195 $x->{_m} = $LIB->_pow($x->{_m}, $y1); 2196 $x->{_e} = $LIB->_mul ($x->{_e}, $y1); 2197 2198 $x->{sign} = $new_sign; 2199 $x->bnorm(); 2200 if ($y->{sign} eq '-') { 2201 # modify $x in place! 2202 my $z = $x->copy(); $x->bone(); 2203 return scalar $x->bdiv($z, $a, $p, $r); # round in one go (might ignore y's A!) 2204 } 2205 $x->round($a, $p, $r, $y); 2206} 2207 2208sub blog { 2209 # Return the logarithm of the operand. If a second operand is defined, that 2210 # value is used as the base, otherwise the base is assumed to be Euler's 2211 # constant. 2212 2213 my ($class, $x, $base, $a, $p, $r); 2214 2215 # Don't objectify the base, since an undefined base, as in $x->blog() or 2216 # $x->blog(undef) signals that the base is Euler's number. 2217 2218 if (!ref($_[0]) && $_[0] =~ /^[A-Za-z]|::/) { 2219 # E.g., Math::BigFloat->blog(256, 2) 2220 ($class, $x, $base, $a, $p, $r) = 2221 defined $_[2] ? objectify(2, @_) : objectify(1, @_); 2222 } else { 2223 # E.g., Math::BigFloat::blog(256, 2) or $x->blog(2) 2224 ($class, $x, $base, $a, $p, $r) = 2225 defined $_[1] ? objectify(2, @_) : objectify(1, @_); 2226 } 2227 2228 return $x if $x->modify('blog'); 2229 2230 return $x -> bnan() if $x -> is_nan(); 2231 2232 # we need to limit the accuracy to protect against overflow 2233 my $fallback = 0; 2234 my ($scale, @params); 2235 ($x, @params) = $x->_find_round_parameters($a, $p, $r); 2236 2237 # no rounding at all, so must use fallback 2238 if (scalar @params == 0) { 2239 # simulate old behaviour 2240 $params[0] = $class->div_scale(); # and round to it as accuracy 2241 $params[1] = undef; # P = undef 2242 $scale = $params[0]+4; # at least four more for proper round 2243 $params[2] = $r; # round mode by caller or undef 2244 $fallback = 1; # to clear a/p afterwards 2245 } else { 2246 # the 4 below is empirical, and there might be cases where it is not 2247 # enough... 2248 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 2249 } 2250 2251 my $done = 0; 2252 if (defined $base) { 2253 $base = $class -> new($base) unless ref $base; 2254 if ($base -> is_nan() || $base -> is_one()) { 2255 $x -> bnan(); 2256 $done = 1; 2257 } elsif ($base -> is_inf() || $base -> is_zero()) { 2258 if ($x -> is_inf() || $x -> is_zero()) { 2259 $x -> bnan(); 2260 } else { 2261 $x -> bzero(@params); 2262 } 2263 $done = 1; 2264 } elsif ($base -> is_negative()) { # -inf < base < 0 2265 if ($x -> is_one()) { # x = 1 2266 $x -> bzero(@params); 2267 } elsif ($x == $base) { 2268 $x -> bone('+', @params); # x = base 2269 } else { 2270 $x -> bnan(); # otherwise 2271 } 2272 $done = 1; 2273 } elsif ($x == $base) { 2274 $x -> bone('+', @params); # 0 < base && 0 < x < inf 2275 $done = 1; 2276 } 2277 } 2278 2279 # We now know that the base is either undefined or positive and finite. 2280 2281 unless ($done) { 2282 if ($x -> is_inf()) { # x = +/-inf 2283 my $sign = defined $base && $base < 1 ? '-' : '+'; 2284 $x -> binf($sign); 2285 $done = 1; 2286 } elsif ($x -> is_neg()) { # -inf < x < 0 2287 $x -> bnan(); 2288 $done = 1; 2289 } elsif ($x -> is_one()) { # x = 1 2290 $x -> bzero(@params); 2291 $done = 1; 2292 } elsif ($x -> is_zero()) { # x = 0 2293 my $sign = defined $base && $base < 1 ? '+' : '-'; 2294 $x -> binf($sign); 2295 $done = 1; 2296 } 2297 } 2298 2299 if ($done) { 2300 if ($fallback) { 2301 # clear a/p after round, since user did not request it 2302 delete $x->{_a}; 2303 delete $x->{_p}; 2304 } 2305 return $x; 2306 } 2307 2308 # when user set globals, they would interfere with our calculation, so 2309 # disable them and later re-enable them 2310 no strict 'refs'; 2311 my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef; 2312 my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef; 2313 # we also need to disable any set A or P on $x (_find_round_parameters took 2314 # them already into account), since these would interfere, too 2315 delete $x->{_a}; delete $x->{_p}; 2316 # need to disable $upgrade in BigInt, to avoid deep recursion 2317 local $Math::BigInt::upgrade = undef; 2318 local $Math::BigFloat::downgrade = undef; 2319 2320 # upgrade $x if $x is not a Math::BigFloat (handle BigInt input) 2321 # XXX TODO: rebless! 2322 if (!$x->isa('Math::BigFloat')) { 2323 $x = Math::BigFloat->new($x); 2324 $class = ref($x); 2325 } 2326 2327 $done = 0; 2328 2329 # If the base is defined and an integer, try to calculate integer result 2330 # first. This is very fast, and in case the real result was found, we can 2331 # stop right here. 2332 if (defined $base && $base->is_int() && $x->is_int()) { 2333 my $xint = Math::BigInt -> new($x -> bdstr()); 2334 my $bint = Math::BigInt -> new($base -> bdstr()); 2335 $xint->blog($bint); 2336 2337 # if we found the exact result, we're done 2338 if ($bint -> bpow($xint) == $x) { 2339 my $xflt = Math::BigFloat -> new($xint -> bdstr()); 2340 $x->{sign} = $xflt->{sign}; 2341 $x->{_m} = $xflt->{_m}; 2342 $x->{_es} = $xflt->{_es}; 2343 $x->{_e} = $xflt->{_e}; 2344 $done = 1; 2345 } 2346 } 2347 2348 if ($done == 0) { 2349 # First calculate the log to base e (using reduction by 10 and possibly 2350 # also by 2): 2351 $x->_log_10($scale); 2352 2353 # and if a different base was requested, convert it 2354 if (defined $base) { 2355 $base = Math::BigFloat->new($base) 2356 unless $base->isa('Math::BigFloat'); 2357 # log_b(x) = ln(x) / ln(b), so compute ln(b) 2358 my $base_log_e = $base->copy()->_log_10($scale); 2359 $x->bdiv($base_log_e, $scale); 2360 } 2361 } 2362 2363 # shortcut to not run through _find_round_parameters again 2364 if (defined $params[0]) { 2365 $x->bround($params[0], $params[2]); # then round accordingly 2366 } else { 2367 $x->bfround($params[1], $params[2]); # then round accordingly 2368 } 2369 if ($fallback) { 2370 # clear a/p after round, since user did not request it 2371 delete $x->{_a}; 2372 delete $x->{_p}; 2373 } 2374 # restore globals 2375 $$abr = $ab; 2376 $$pbr = $pb; 2377 2378 $x; 2379} 2380 2381sub bexp { 2382 # Calculate e ** X (Euler's number to the power of X) 2383 my ($class, $x, $a, $p, $r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 2384 2385 return $x if $x->modify('bexp'); 2386 2387 return $x->binf() if $x->{sign} eq '+inf'; 2388 return $x->bzero() if $x->{sign} eq '-inf'; 2389 2390 # we need to limit the accuracy to protect against overflow 2391 my $fallback = 0; 2392 my ($scale, @params); 2393 ($x, @params) = $x->_find_round_parameters($a, $p, $r); 2394 2395 # also takes care of the "error in _find_round_parameters?" case 2396 return $x if $x->{sign} eq 'NaN'; 2397 2398 # no rounding at all, so must use fallback 2399 if (scalar @params == 0) { 2400 # simulate old behaviour 2401 $params[0] = $class->div_scale(); # and round to it as accuracy 2402 $params[1] = undef; # P = undef 2403 $scale = $params[0]+4; # at least four more for proper round 2404 $params[2] = $r; # round mode by caller or undef 2405 $fallback = 1; # to clear a/p afterwards 2406 } else { 2407 # the 4 below is empirical, and there might be cases where it's not 2408 # enough ... 2409 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 2410 } 2411 2412 return $x->bone(@params) if $x->is_zero(); 2413 2414 if (!$x->isa('Math::BigFloat')) { 2415 $x = Math::BigFloat->new($x); 2416 $class = ref($x); 2417 } 2418 2419 # when user set globals, they would interfere with our calculation, so 2420 # disable them and later re-enable them 2421 no strict 'refs'; 2422 my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef; 2423 my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef; 2424 # we also need to disable any set A or P on $x (_find_round_parameters took 2425 # them already into account), since these would interfere, too 2426 delete $x->{_a}; 2427 delete $x->{_p}; 2428 # need to disable $upgrade in BigInt, to avoid deep recursion 2429 local $Math::BigInt::upgrade = undef; 2430 local $Math::BigFloat::downgrade = undef; 2431 2432 my $x_org = $x->copy(); 2433 2434 # We use the following Taylor series: 2435 2436 # x x^2 x^3 x^4 2437 # e = 1 + --- + --- + --- + --- ... 2438 # 1! 2! 3! 4! 2439 2440 # The difference for each term is X and N, which would result in: 2441 # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term 2442 2443 # But it is faster to compute exp(1) and then raising it to the 2444 # given power, esp. if $x is really big and an integer because: 2445 2446 # * The numerator is always 1, making the computation faster 2447 # * the series converges faster in the case of x == 1 2448 # * We can also easily check when we have reached our limit: when the 2449 # term to be added is smaller than "1E$scale", we can stop - f.i. 2450 # scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5. 2451 # * we can compute the *exact* result by simulating bigrat math: 2452 2453 # 1 1 gcd(3, 4) = 1 1*24 + 1*6 5 2454 # - + - = ---------- = -- 2455 # 6 24 6*24 24 2456 2457 # We do not compute the gcd() here, but simple do: 2458 # 1 1 1*24 + 1*6 30 2459 # - + - = --------- = -- 2460 # 6 24 6*24 144 2461 2462 # In general: 2463 # a c a*d + c*b and note that c is always 1 and d = (b*f) 2464 # - + - = --------- 2465 # b d b*d 2466 2467 # This leads to: which can be reduced by b to: 2468 # a 1 a*b*f + b a*f + 1 2469 # - + - = --------- = ------- 2470 # b b*f b*b*f b*f 2471 2472 # The first terms in the series are: 2473 2474 # 1 1 1 1 1 1 1 1 13700 2475 # -- + -- + -- + -- + -- + --- + --- + ---- = ----- 2476 # 1 1 2 6 24 120 720 5040 5040 2477 2478 # Note that we cannot simply reduce 13700/5040 to 685/252, but must keep 2479 # the numerator and the denominator! 2480 2481 if ($scale <= 75) { 2482 # set $x directly from a cached string form 2483 $x->{_m} = $LIB->_new("2718281828459045235360287471352662497757" . 2484 "2470936999595749669676277240766303535476"); 2485 $x->{sign} = '+'; 2486 $x->{_es} = '-'; 2487 $x->{_e} = $LIB->_new(79); 2488 } else { 2489 # compute A and B so that e = A / B. 2490 2491 # After some terms we end up with this, so we use it as a starting point: 2492 my $A = $LIB->_new("9093339520860578540197197" . 2493 "0164779391644753259799242"); 2494 my $F = $LIB->_new(42); 2495 my $step = 42; 2496 2497 # Compute how many steps we need to take to get $A and $B sufficiently big 2498 my $steps = _len_to_steps($scale - 4); 2499 # print STDERR "# Doing $steps steps for ", $scale-4, " digits\n"; 2500 while ($step++ <= $steps) { 2501 # calculate $a * $f + 1 2502 $A = $LIB->_mul($A, $F); 2503 $A = $LIB->_inc($A); 2504 # increment f 2505 $F = $LIB->_inc($F); 2506 } 2507 # compute $B as factorial of $steps (this is faster than doing it manually) 2508 my $B = $LIB->_fac($LIB->_new($steps)); 2509 2510 # print "A ", $LIB->_str($A), "\nB ", $LIB->_str($B), "\n"; 2511 2512 # compute A/B with $scale digits in the result (truncate, not round) 2513 $A = $LIB->_lsft($A, $LIB->_new($scale), 10); 2514 $A = $LIB->_div($A, $B); 2515 2516 $x->{_m} = $A; 2517 $x->{sign} = '+'; 2518 $x->{_es} = '-'; 2519 $x->{_e} = $LIB->_new($scale); 2520 } 2521 2522 # $x contains now an estimate of e, with some surplus digits, so we can round 2523 if (!$x_org->is_one()) { 2524 # Reduce size of fractional part, followup with integer power of two. 2525 my $lshift = 0; 2526 while ($lshift < 30 && $x_org->bacmp(2 << $lshift) > 0) { 2527 $lshift++; 2528 } 2529 # Raise $x to the wanted power and round it. 2530 if ($lshift == 0) { 2531 $x->bpow($x_org, @params); 2532 } else { 2533 my($mul, $rescale) = (1 << $lshift, $scale+1+$lshift); 2534 $x->bpow(scalar $x_org->bdiv($mul, $rescale), $rescale)->bpow($mul, @params); 2535 } 2536 } else { 2537 # else just round the already computed result 2538 delete $x->{_a}; 2539 delete $x->{_p}; 2540 # shortcut to not run through _find_round_parameters again 2541 if (defined $params[0]) { 2542 $x->bround($params[0], $params[2]); # then round accordingly 2543 } else { 2544 $x->bfround($params[1], $params[2]); # then round accordingly 2545 } 2546 } 2547 if ($fallback) { 2548 # clear a/p after round, since user did not request it 2549 delete $x->{_a}; 2550 delete $x->{_p}; 2551 } 2552 # restore globals 2553 $$abr = $ab; 2554 $$pbr = $pb; 2555 2556 $x; # return modified $x 2557} 2558 2559sub bnok { 2560 # Calculate n over k (binomial coefficient or "choose" function) as integer. 2561 # set up parameters 2562 my ($class, $x, $y, @r) = (ref($_[0]), @_); 2563 2564 # objectify is costly, so avoid it 2565 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 2566 ($class, $x, $y, @r) = objectify(2, @_); 2567 } 2568 2569 return $x if $x->modify('bnok'); 2570 2571 return $x->bnan() if $x->is_nan() || $y->is_nan(); 2572 return $x->bnan() if (($x->is_finite() && !$x->is_int()) || 2573 ($y->is_finite() && !$y->is_int())); 2574 2575 my $xint = Math::BigInt -> new($x -> bsstr()); 2576 my $yint = Math::BigInt -> new($y -> bsstr()); 2577 $xint -> bnok($yint); 2578 my $xflt = Math::BigFloat -> new($xint); 2579 2580 $x->{_m} = $xflt->{_m}; 2581 $x->{_e} = $xflt->{_e}; 2582 $x->{_es} = $xflt->{_es}; 2583 $x->{sign} = $xflt->{sign}; 2584 2585 return $x; 2586} 2587 2588sub bsin { 2589 # Calculate a sinus of x. 2590 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 2591 2592 # taylor: x^3 x^5 x^7 x^9 2593 # sin = x - --- + --- - --- + --- ... 2594 # 3! 5! 7! 9! 2595 2596 # we need to limit the accuracy to protect against overflow 2597 my $fallback = 0; 2598 my ($scale, @params); 2599 ($x, @params) = $x->_find_round_parameters(@r); 2600 2601 # constant object or error in _find_round_parameters? 2602 return $x if $x->modify('bsin') || $x->is_nan(); 2603 2604 return $x->bzero(@r) if $x->is_zero(); 2605 2606 # no rounding at all, so must use fallback 2607 if (scalar @params == 0) { 2608 # simulate old behaviour 2609 $params[0] = $class->div_scale(); # and round to it as accuracy 2610 $params[1] = undef; # disable P 2611 $scale = $params[0]+4; # at least four more for proper round 2612 $params[2] = $r[2]; # round mode by caller or undef 2613 $fallback = 1; # to clear a/p afterwards 2614 } else { 2615 # the 4 below is empirical, and there might be cases where it is not 2616 # enough... 2617 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 2618 } 2619 2620 # when user set globals, they would interfere with our calculation, so 2621 # disable them and later re-enable them 2622 no strict 'refs'; 2623 my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef; 2624 my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef; 2625 # we also need to disable any set A or P on $x (_find_round_parameters took 2626 # them already into account), since these would interfere, too 2627 delete $x->{_a}; 2628 delete $x->{_p}; 2629 # need to disable $upgrade in BigInt, to avoid deep recursion 2630 local $Math::BigInt::upgrade = undef; 2631 2632 my $over = $x * $x; # X ^ 2 2633 my $x2 = $over->copy(); # X ^ 2; difference between terms 2634 $over->bmul($x); # X ^ 3 as starting value 2635 my $sign = 1; # start with -= 2636 my $below = $class->new(6); my $factorial = $class->new(4); 2637 delete $x->{_a}; 2638 delete $x->{_p}; 2639 2640 my $limit = $class->new("1E-". ($scale-1)); 2641 #my $steps = 0; 2642 while (3 < 5) { 2643 # we calculate the next term, and add it to the last 2644 # when the next term is below our limit, it won't affect the outcome 2645 # anymore, so we stop: 2646 my $next = $over->copy()->bdiv($below, $scale); 2647 last if $next->bacmp($limit) <= 0; 2648 2649 if ($sign == 0) { 2650 $x->badd($next); 2651 } else { 2652 $x->bsub($next); 2653 } 2654 $sign = 1-$sign; # alternate 2655 # calculate things for the next term 2656 $over->bmul($x2); # $x*$x 2657 $below->bmul($factorial); $factorial->binc(); # n*(n+1) 2658 $below->bmul($factorial); $factorial->binc(); # n*(n+1) 2659 } 2660 2661 # shortcut to not run through _find_round_parameters again 2662 if (defined $params[0]) { 2663 $x->bround($params[0], $params[2]); # then round accordingly 2664 } else { 2665 $x->bfround($params[1], $params[2]); # then round accordingly 2666 } 2667 if ($fallback) { 2668 # clear a/p after round, since user did not request it 2669 delete $x->{_a}; 2670 delete $x->{_p}; 2671 } 2672 # restore globals 2673 $$abr = $ab; 2674 $$pbr = $pb; 2675 $x; 2676} 2677 2678sub bcos { 2679 # Calculate a cosinus of x. 2680 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 2681 2682 # Taylor: x^2 x^4 x^6 x^8 2683 # cos = 1 - --- + --- - --- + --- ... 2684 # 2! 4! 6! 8! 2685 2686 # we need to limit the accuracy to protect against overflow 2687 my $fallback = 0; 2688 my ($scale, @params); 2689 ($x, @params) = $x->_find_round_parameters(@r); 2690 2691 # constant object or error in _find_round_parameters? 2692 return $x if $x->modify('bcos') || $x->is_nan(); 2693 2694 return $x->bone(@r) if $x->is_zero(); 2695 2696 # no rounding at all, so must use fallback 2697 if (scalar @params == 0) { 2698 # simulate old behaviour 2699 $params[0] = $class->div_scale(); # and round to it as accuracy 2700 $params[1] = undef; # disable P 2701 $scale = $params[0]+4; # at least four more for proper round 2702 $params[2] = $r[2]; # round mode by caller or undef 2703 $fallback = 1; # to clear a/p afterwards 2704 } else { 2705 # the 4 below is empirical, and there might be cases where it is not 2706 # enough... 2707 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 2708 } 2709 2710 # when user set globals, they would interfere with our calculation, so 2711 # disable them and later re-enable them 2712 no strict 'refs'; 2713 my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef; 2714 my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef; 2715 # we also need to disable any set A or P on $x (_find_round_parameters took 2716 # them already into account), since these would interfere, too 2717 delete $x->{_a}; delete $x->{_p}; 2718 # need to disable $upgrade in BigInt, to avoid deep recursion 2719 local $Math::BigInt::upgrade = undef; 2720 2721 my $over = $x * $x; # X ^ 2 2722 my $x2 = $over->copy(); # X ^ 2; difference between terms 2723 my $sign = 1; # start with -= 2724 my $below = $class->new(2); 2725 my $factorial = $class->new(3); 2726 $x->bone(); 2727 delete $x->{_a}; 2728 delete $x->{_p}; 2729 2730 my $limit = $class->new("1E-". ($scale-1)); 2731 #my $steps = 0; 2732 while (3 < 5) { 2733 # we calculate the next term, and add it to the last 2734 # when the next term is below our limit, it won't affect the outcome 2735 # anymore, so we stop: 2736 my $next = $over->copy()->bdiv($below, $scale); 2737 last if $next->bacmp($limit) <= 0; 2738 2739 if ($sign == 0) { 2740 $x->badd($next); 2741 } else { 2742 $x->bsub($next); 2743 } 2744 $sign = 1-$sign; # alternate 2745 # calculate things for the next term 2746 $over->bmul($x2); # $x*$x 2747 $below->bmul($factorial); $factorial->binc(); # n*(n+1) 2748 $below->bmul($factorial); $factorial->binc(); # n*(n+1) 2749 } 2750 2751 # shortcut to not run through _find_round_parameters again 2752 if (defined $params[0]) { 2753 $x->bround($params[0], $params[2]); # then round accordingly 2754 } else { 2755 $x->bfround($params[1], $params[2]); # then round accordingly 2756 } 2757 if ($fallback) { 2758 # clear a/p after round, since user did not request it 2759 delete $x->{_a}; 2760 delete $x->{_p}; 2761 } 2762 # restore globals 2763 $$abr = $ab; 2764 $$pbr = $pb; 2765 $x; 2766} 2767 2768sub batan { 2769 # Calculate a arcus tangens of x. 2770 2771 my $self = shift; 2772 my $selfref = ref $self; 2773 my $class = $selfref || $self; 2774 2775 my (@r) = @_; 2776 2777 # taylor: x^3 x^5 x^7 x^9 2778 # atan = x - --- + --- - --- + --- ... 2779 # 3 5 7 9 2780 2781 # We need to limit the accuracy to protect against overflow. 2782 2783 my $fallback = 0; 2784 my ($scale, @params); 2785 ($self, @params) = $self->_find_round_parameters(@r); 2786 2787 # Constant object or error in _find_round_parameters? 2788 2789 return $self if $self->modify('batan') || $self->is_nan(); 2790 2791 if ($self->{sign} =~ /^[+-]inf\z/) { 2792 # +inf result is PI/2 2793 # -inf result is -PI/2 2794 # calculate PI/2 2795 my $pi = $class->bpi(@r); 2796 # modify $self in place 2797 $self->{_m} = $pi->{_m}; 2798 $self->{_e} = $pi->{_e}; 2799 $self->{_es} = $pi->{_es}; 2800 # -y => -PI/2, +y => PI/2 2801 $self->{sign} = substr($self->{sign}, 0, 1); # "+inf" => "+" 2802 $self -> {_m} = $LIB->_div($self->{_m}, $LIB->_new(2)); 2803 return $self; 2804 } 2805 2806 return $self->bzero(@r) if $self->is_zero(); 2807 2808 # no rounding at all, so must use fallback 2809 if (scalar @params == 0) { 2810 # simulate old behaviour 2811 $params[0] = $class->div_scale(); # and round to it as accuracy 2812 $params[1] = undef; # disable P 2813 $scale = $params[0]+4; # at least four more for proper round 2814 $params[2] = $r[2]; # round mode by caller or undef 2815 $fallback = 1; # to clear a/p afterwards 2816 } else { 2817 # the 4 below is empirical, and there might be cases where it is not 2818 # enough... 2819 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 2820 } 2821 2822 # 1 or -1 => PI/4 2823 # inlined is_one() && is_one('-') 2824 if ($LIB->_is_one($self->{_m}) && $LIB->_is_zero($self->{_e})) { 2825 my $pi = $class->bpi($scale - 3); 2826 # modify $self in place 2827 $self->{_m} = $pi->{_m}; 2828 $self->{_e} = $pi->{_e}; 2829 $self->{_es} = $pi->{_es}; 2830 # leave the sign of $self alone (+1 => +PI/4, -1 => -PI/4) 2831 $self->{_m} = $LIB->_div($self->{_m}, $LIB->_new(4)); 2832 return $self; 2833 } 2834 2835 # This series is only valid if -1 < x < 1, so for other x we need to 2836 # calculate PI/2 - atan(1/x): 2837 my $pi = undef; 2838 if ($self->bacmp($self->copy()->bone) >= 0) { 2839 # calculate PI/2 2840 $pi = $class->bpi($scale - 3); 2841 $pi->{_m} = $LIB->_div($pi->{_m}, $LIB->_new(2)); 2842 # calculate 1/$self: 2843 my $self_copy = $self->copy(); 2844 # modify $self in place 2845 $self->bone(); 2846 $self->bdiv($self_copy, $scale); 2847 } 2848 2849 my $fmul = 1; 2850 foreach (0 .. int($scale / 20)) { 2851 $fmul *= 2; 2852 $self->bdiv($self->copy()->bmul($self)->binc->bsqrt($scale + 4)->binc, $scale + 4); 2853 } 2854 2855 # When user set globals, they would interfere with our calculation, so 2856 # disable them and later re-enable them. 2857 no strict 'refs'; 2858 my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef; 2859 my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef; 2860 # We also need to disable any set A or P on $self (_find_round_parameters 2861 # took them already into account), since these would interfere, too 2862 delete $self->{_a}; 2863 delete $self->{_p}; 2864 # Need to disable $upgrade in BigInt, to avoid deep recursion. 2865 local $Math::BigInt::upgrade = undef; 2866 2867 my $over = $self * $self; # X ^ 2 2868 my $self2 = $over->copy(); # X ^ 2; difference between terms 2869 $over->bmul($self); # X ^ 3 as starting value 2870 my $sign = 1; # start with -= 2871 my $below = $class->new(3); 2872 my $two = $class->new(2); 2873 delete $self->{_a}; 2874 delete $self->{_p}; 2875 2876 my $limit = $class->new("1E-". ($scale-1)); 2877 #my $steps = 0; 2878 while (1) { 2879 # We calculate the next term, and add it to the last. When the next 2880 # term is below our limit, it won't affect the outcome anymore, so we 2881 # stop: 2882 my $next = $over->copy()->bdiv($below, $scale); 2883 last if $next->bacmp($limit) <= 0; 2884 2885 if ($sign == 0) { 2886 $self->badd($next); 2887 } else { 2888 $self->bsub($next); 2889 } 2890 $sign = 1-$sign; # alternatex 2891 # calculate things for the next term 2892 $over->bmul($self2); # $self*$self 2893 $below->badd($two); # n += 2 2894 } 2895 $self->bmul($fmul); 2896 2897 if (defined $pi) { 2898 my $self_copy = $self->copy(); 2899 # modify $self in place 2900 $self->{_m} = $pi->{_m}; 2901 $self->{_e} = $pi->{_e}; 2902 $self->{_es} = $pi->{_es}; 2903 # PI/2 - $self 2904 $self->bsub($self_copy); 2905 } 2906 2907 # Shortcut to not run through _find_round_parameters again. 2908 if (defined $params[0]) { 2909 $self->bround($params[0], $params[2]); # then round accordingly 2910 } else { 2911 $self->bfround($params[1], $params[2]); # then round accordingly 2912 } 2913 if ($fallback) { 2914 # Clear a/p after round, since user did not request it. 2915 delete $self->{_a}; 2916 delete $self->{_p}; 2917 } 2918 2919 # restore globals 2920 $$abr = $ab; 2921 $$pbr = $pb; 2922 $self; 2923} 2924 2925sub batan2 { 2926 # $y -> batan2($x) returns the arcus tangens of $y / $x. 2927 2928 # Set up parameters. 2929 my ($class, $y, $x, @r) = (ref($_[0]), @_); 2930 2931 # Objectify is costly, so avoid it if we can. 2932 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 2933 ($class, $y, $x, @r) = objectify(2, @_); 2934 } 2935 2936 # Quick exit if $y is read-only. 2937 return $y if $y -> modify('batan2'); 2938 2939 # Handle all NaN cases. 2940 return $y -> bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan; 2941 2942 # We need to limit the accuracy to protect against overflow. 2943 my $fallback = 0; 2944 my ($scale, @params); 2945 ($y, @params) = $y -> _find_round_parameters(@r); 2946 2947 # Error in _find_round_parameters? 2948 return $y if $y->is_nan(); 2949 2950 # No rounding at all, so must use fallback. 2951 if (scalar @params == 0) { 2952 # Simulate old behaviour 2953 $params[0] = $class -> div_scale(); # and round to it as accuracy 2954 $params[1] = undef; # disable P 2955 $scale = $params[0] + 4; # at least four more for proper round 2956 $params[2] = $r[2]; # round mode by caller or undef 2957 $fallback = 1; # to clear a/p afterwards 2958 } else { 2959 # The 4 below is empirical, and there might be cases where it is not 2960 # enough ... 2961 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 2962 } 2963 2964 if ($x -> is_inf("+")) { # x = inf 2965 if ($y -> is_inf("+")) { # y = inf 2966 $y -> bpi($scale) -> bmul("0.25"); # pi/4 2967 } elsif ($y -> is_inf("-")) { # y = -inf 2968 $y -> bpi($scale) -> bmul("-0.25"); # -pi/4 2969 } else { # -inf < y < inf 2970 return $y -> bzero(@r); # 0 2971 } 2972 } elsif ($x -> is_inf("-")) { # x = -inf 2973 if ($y -> is_inf("+")) { # y = inf 2974 $y -> bpi($scale) -> bmul("0.75"); # 3/4 pi 2975 } elsif ($y -> is_inf("-")) { # y = -inf 2976 $y -> bpi($scale) -> bmul("-0.75"); # -3/4 pi 2977 } elsif ($y >= 0) { # y >= 0 2978 $y -> bpi($scale); # pi 2979 } else { # y < 0 2980 $y -> bpi($scale) -> bneg(); # -pi 2981 } 2982 } elsif ($x > 0) { # 0 < x < inf 2983 if ($y -> is_inf("+")) { # y = inf 2984 $y -> bpi($scale) -> bmul("0.5"); # pi/2 2985 } elsif ($y -> is_inf("-")) { # y = -inf 2986 $y -> bpi($scale) -> bmul("-0.5"); # -pi/2 2987 } else { # -inf < y < inf 2988 $y -> bdiv($x, $scale) -> batan($scale); # atan(y/x) 2989 } 2990 } elsif ($x < 0) { # -inf < x < 0 2991 my $pi = $class -> bpi($scale); 2992 if ($y >= 0) { # y >= 0 2993 $y -> bdiv($x, $scale) -> batan() # atan(y/x) + pi 2994 -> badd($pi); 2995 } else { # y < 0 2996 $y -> bdiv($x, $scale) -> batan() # atan(y/x) - pi 2997 -> bsub($pi); 2998 } 2999 } else { # x = 0 3000 if ($y > 0) { # y > 0 3001 $y -> bpi($scale) -> bmul("0.5"); # pi/2 3002 } elsif ($y < 0) { # y < 0 3003 $y -> bpi($scale) -> bmul("-0.5"); # -pi/2 3004 } else { # y = 0 3005 return $y -> bzero(@r); # 0 3006 } 3007 } 3008 3009 $y -> round(@r); 3010 3011 if ($fallback) { 3012 delete $y->{_a}; 3013 delete $y->{_p}; 3014 } 3015 3016 return $y; 3017} 3018############################################################################## 3019 3020sub bsqrt { 3021 # calculate square root 3022 my ($class, $x, $a, $p, $r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 3023 3024 return $x if $x->modify('bsqrt'); 3025 3026 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0 3027 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf 3028 return $x->round($a, $p, $r) if $x->is_zero() || $x->is_one(); 3029 3030 # we need to limit the accuracy to protect against overflow 3031 my $fallback = 0; 3032 my (@params, $scale); 3033 ($x, @params) = $x->_find_round_parameters($a, $p, $r); 3034 3035 return $x if $x->is_nan(); # error in _find_round_parameters? 3036 3037 # no rounding at all, so must use fallback 3038 if (scalar @params == 0) { 3039 # simulate old behaviour 3040 $params[0] = $class->div_scale(); # and round to it as accuracy 3041 $scale = $params[0]+4; # at least four more for proper round 3042 $params[2] = $r; # round mode by caller or undef 3043 $fallback = 1; # to clear a/p afterwards 3044 } else { 3045 # the 4 below is empirical, and there might be cases where it is not 3046 # enough... 3047 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 3048 } 3049 3050 # when user set globals, they would interfere with our calculation, so 3051 # disable them and later re-enable them 3052 no strict 'refs'; 3053 my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef; 3054 my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef; 3055 # we also need to disable any set A or P on $x (_find_round_parameters took 3056 # them already into account), since these would interfere, too 3057 delete $x->{_a}; 3058 delete $x->{_p}; 3059 # need to disable $upgrade in BigInt, to avoid deep recursion 3060 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI 3061 3062 my $i = $LIB->_copy($x->{_m}); 3063 $i = $LIB->_lsft($i, $x->{_e}, 10) unless $LIB->_is_zero($x->{_e}); 3064 my $xas = Math::BigInt->bzero(); 3065 $xas->{value} = $i; 3066 3067 my $gs = $xas->copy()->bsqrt(); # some guess 3068 3069 if (($x->{_es} ne '-') # guess can't be accurate if there are 3070 # digits after the dot 3071 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head? 3072 { 3073 # exact result, copy result over to keep $x 3074 $x->{_m} = $gs->{value}; 3075 $x->{_e} = $LIB->_zero(); 3076 $x->{_es} = '+'; 3077 $x->bnorm(); 3078 # shortcut to not run through _find_round_parameters again 3079 if (defined $params[0]) { 3080 $x->bround($params[0], $params[2]); # then round accordingly 3081 } else { 3082 $x->bfround($params[1], $params[2]); # then round accordingly 3083 } 3084 if ($fallback) { 3085 # clear a/p after round, since user did not request it 3086 delete $x->{_a}; 3087 delete $x->{_p}; 3088 } 3089 # re-enable A and P, upgrade is taken care of by "local" 3090 ${"$class\::accuracy"} = $ab; 3091 ${"$class\::precision"} = $pb; 3092 return $x; 3093 } 3094 3095 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy 3096 # of the result by multiplying the input by 100 and then divide the integer 3097 # result of sqrt(input) by 10. Rounding afterwards returns the real result. 3098 3099 # The following steps will transform 123.456 (in $x) into 123456 (in $y1) 3100 my $y1 = $LIB->_copy($x->{_m}); 3101 3102 my $length = $LIB->_len($y1); 3103 3104 # Now calculate how many digits the result of sqrt(y1) would have 3105 my $digits = int($length / 2); 3106 3107 # But we need at least $scale digits, so calculate how many are missing 3108 my $shift = $scale - $digits; 3109 3110 # This happens if the input had enough digits 3111 # (we take care of integer guesses above) 3112 $shift = 0 if $shift < 0; 3113 3114 # Multiply in steps of 100, by shifting left two times the "missing" digits 3115 my $s2 = $shift * 2; 3116 3117 # We now make sure that $y1 has the same odd or even number of digits than 3118 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left, 3119 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not 3120 # steps of 10. The length of $x does not count, since an even or odd number 3121 # of digits before the dot is not changed by adding an even number of digits 3122 # after the dot (the result is still odd or even digits long). 3123 $s2++ if $LIB->_is_odd($x->{_e}); 3124 3125 $y1 = $LIB->_lsft($y1, $LIB->_new($s2), 10); 3126 3127 # now take the square root and truncate to integer 3128 $y1 = $LIB->_sqrt($y1); 3129 3130 # By "shifting" $y1 right (by creating a negative _e) we calculate the final 3131 # result, which is than later rounded to the desired scale. 3132 3133 # calculate how many zeros $x had after the '.' (or before it, depending 3134 # on sign of $dat, the result should have half as many: 3135 my $dat = $LIB->_num($x->{_e}); 3136 $dat = -$dat if $x->{_es} eq '-'; 3137 $dat += $length; 3138 3139 if ($dat > 0) { 3140 # no zeros after the dot (e.g. 1.23, 0.49 etc) 3141 # preserve half as many digits before the dot than the input had 3142 # (but round this "up") 3143 $dat = int(($dat+1)/2); 3144 } else { 3145 $dat = int(($dat)/2); 3146 } 3147 $dat -= $LIB->_len($y1); 3148 if ($dat < 0) { 3149 $dat = abs($dat); 3150 $x->{_e} = $LIB->_new($dat); 3151 $x->{_es} = '-'; 3152 } else { 3153 $x->{_e} = $LIB->_new($dat); 3154 $x->{_es} = '+'; 3155 } 3156 $x->{_m} = $y1; 3157 $x->bnorm(); 3158 3159 # shortcut to not run through _find_round_parameters again 3160 if (defined $params[0]) { 3161 $x->bround($params[0], $params[2]); # then round accordingly 3162 } else { 3163 $x->bfround($params[1], $params[2]); # then round accordingly 3164 } 3165 if ($fallback) { 3166 # clear a/p after round, since user did not request it 3167 delete $x->{_a}; 3168 delete $x->{_p}; 3169 } 3170 # restore globals 3171 $$abr = $ab; 3172 $$pbr = $pb; 3173 $x; 3174} 3175 3176sub broot { 3177 # calculate $y'th root of $x 3178 3179 # set up parameters 3180 my ($class, $x, $y, $a, $p, $r) = (ref($_[0]), @_); 3181 # objectify is costly, so avoid it 3182 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 3183 ($class, $x, $y, $a, $p, $r) = objectify(2, @_); 3184 } 3185 3186 return $x if $x->modify('broot'); 3187 3188 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0 3189 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() || 3190 $y->{sign} !~ /^\+$/; 3191 3192 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one(); 3193 3194 # we need to limit the accuracy to protect against overflow 3195 my $fallback = 0; 3196 my (@params, $scale); 3197 ($x, @params) = $x->_find_round_parameters($a, $p, $r); 3198 3199 return $x if $x->is_nan(); # error in _find_round_parameters? 3200 3201 # no rounding at all, so must use fallback 3202 if (scalar @params == 0) { 3203 # simulate old behaviour 3204 $params[0] = $class->div_scale(); # and round to it as accuracy 3205 $scale = $params[0]+4; # at least four more for proper round 3206 $params[2] = $r; # round mode by caller or undef 3207 $fallback = 1; # to clear a/p afterwards 3208 } else { 3209 # the 4 below is empirical, and there might be cases where it is not 3210 # enough... 3211 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 3212 } 3213 3214 # when user set globals, they would interfere with our calculation, so 3215 # disable them and later re-enable them 3216 no strict 'refs'; 3217 my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef; 3218 my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef; 3219 # we also need to disable any set A or P on $x (_find_round_parameters took 3220 # them already into account), since these would interfere, too 3221 delete $x->{_a}; 3222 delete $x->{_p}; 3223 # need to disable $upgrade in BigInt, to avoid deep recursion 3224 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI 3225 3226 # remember sign and make $x positive, since -4 ** (1/2) => -2 3227 my $sign = 0; 3228 $sign = 1 if $x->{sign} eq '-'; 3229 $x->{sign} = '+'; 3230 3231 my $is_two = 0; 3232 if ($y->isa('Math::BigFloat')) { 3233 $is_two = ($y->{sign} eq '+' && $LIB->_is_two($y->{_m}) && $LIB->_is_zero($y->{_e})); 3234 } else { 3235 $is_two = ($y == 2); 3236 } 3237 3238 # normal square root if $y == 2: 3239 if ($is_two) { 3240 $x->bsqrt($scale+4); 3241 } elsif ($y->is_one('-')) { 3242 # $x ** -1 => 1/$x 3243 my $u = $class->bone()->bdiv($x, $scale); 3244 # copy private parts over 3245 $x->{_m} = $u->{_m}; 3246 $x->{_e} = $u->{_e}; 3247 $x->{_es} = $u->{_es}; 3248 } else { 3249 # calculate the broot() as integer result first, and if it fits, return 3250 # it rightaway (but only if $x and $y are integer): 3251 3252 my $done = 0; # not yet 3253 if ($y->is_int() && $x->is_int()) { 3254 my $i = $LIB->_copy($x->{_m}); 3255 $i = $LIB->_lsft($i, $x->{_e}, 10) unless $LIB->_is_zero($x->{_e}); 3256 my $int = Math::BigInt->bzero(); 3257 $int->{value} = $i; 3258 $int->broot($y->as_number()); 3259 # if ($exact) 3260 if ($int->copy()->bpow($y) == $x) { 3261 # found result, return it 3262 $x->{_m} = $int->{value}; 3263 $x->{_e} = $LIB->_zero(); 3264 $x->{_es} = '+'; 3265 $x->bnorm(); 3266 $done = 1; 3267 } 3268 } 3269 if ($done == 0) { 3270 my $u = $class->bone()->bdiv($y, $scale+4); 3271 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts 3272 $x->bpow($u, $scale+4); # el cheapo 3273 } 3274 } 3275 $x->bneg() if $sign == 1; 3276 3277 # shortcut to not run through _find_round_parameters again 3278 if (defined $params[0]) { 3279 $x->bround($params[0], $params[2]); # then round accordingly 3280 } else { 3281 $x->bfround($params[1], $params[2]); # then round accordingly 3282 } 3283 if ($fallback) { 3284 # clear a/p after round, since user did not request it 3285 delete $x->{_a}; 3286 delete $x->{_p}; 3287 } 3288 # restore globals 3289 $$abr = $ab; 3290 $$pbr = $pb; 3291 $x; 3292} 3293 3294sub bfac { 3295 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT 3296 # compute factorial number, modifies first argument 3297 3298 # set up parameters 3299 my ($class, $x, @r) = (ref($_[0]), @_); 3300 # objectify is costly, so avoid it 3301 ($class, $x, @r) = objectify(1, @_) if !ref($x); 3302 3303 # inf => inf 3304 return $x if $x->modify('bfac') || $x->{sign} eq '+inf'; 3305 3306 return $x->bnan() 3307 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN 3308 ($x->{_es} ne '+')); # digits after dot? 3309 3310 if (! $LIB->_is_zero($x->{_e})) { 3311 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # change 12e1 to 120e0 3312 $x->{_e} = $LIB->_zero(); # normalize 3313 $x->{_es} = '+'; 3314 } 3315 $x->{_m} = $LIB->_fac($x->{_m}); # calculate factorial 3316 $x->bnorm()->round(@r); # norm again and round result 3317} 3318 3319sub bdfac { 3320 # compute double factorial 3321 3322 # set up parameters 3323 my ($class, $x, @r) = (ref($_[0]), @_); 3324 # objectify is costly, so avoid it 3325 ($class, $x, @r) = objectify(1, @_) if !ref($x); 3326 3327 # inf => inf 3328 return $x if $x->modify('bfac') || $x->{sign} eq '+inf'; 3329 3330 return $x->bnan() 3331 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN 3332 ($x->{_es} ne '+')); # digits after dot? 3333 3334 croak("bdfac() requires a newer version of the $LIB library.") 3335 unless $LIB->can('_dfac'); 3336 3337 if (! $LIB->_is_zero($x->{_e})) { 3338 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # change 12e1 to 120e0 3339 $x->{_e} = $LIB->_zero(); # normalize 3340 $x->{_es} = '+'; 3341 } 3342 $x->{_m} = $LIB->_dfac($x->{_m}); # calculate factorial 3343 $x->bnorm()->round(@r); # norm again and round result 3344} 3345 3346sub blsft { 3347 # shift left by $y (multiply by $b ** $y) 3348 3349 # set up parameters 3350 my ($class, $x, $y, $b, $a, $p, $r) = (ref($_[0]), @_); 3351 3352 # objectify is costly, so avoid it 3353 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 3354 ($class, $x, $y, $b, $a, $p, $r) = objectify(2, @_); 3355 } 3356 3357 return $x if $x -> modify('blsft'); 3358 return $x if $x -> {sign} !~ /^[+-]$/; # nan, +inf, -inf 3359 3360 $b = 2 if !defined $b; 3361 $b = $class -> new($b) unless ref($b) && $b -> isa($class); 3362 3363 return $x -> bnan() if $x -> is_nan() || $y -> is_nan() || $b -> is_nan(); 3364 3365 # shift by a negative amount? 3366 return $x -> brsft($y -> copy() -> babs(), $b) if $y -> {sign} =~ /^-/; 3367 3368 $x -> bmul($b -> bpow($y), $a, $p, $r, $y); 3369} 3370 3371sub brsft { 3372 # shift right by $y (divide $b ** $y) 3373 3374 # set up parameters 3375 my ($class, $x, $y, $b, $a, $p, $r) = (ref($_[0]), @_); 3376 3377 # objectify is costly, so avoid it 3378 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 3379 ($class, $x, $y, $b, $a, $p, $r) = objectify(2, @_); 3380 } 3381 3382 return $x if $x -> modify('brsft'); 3383 return $x if $x -> {sign} !~ /^[+-]$/; # nan, +inf, -inf 3384 3385 $b = 2 if !defined $b; 3386 $b = $class -> new($b) unless ref($b) && $b -> isa($class); 3387 3388 return $x -> bnan() if $x -> is_nan() || $y -> is_nan() || $b -> is_nan(); 3389 3390 # shift by a negative amount? 3391 return $x -> blsft($y -> copy() -> babs(), $b) if $y -> {sign} =~ /^-/; 3392 3393 # the following call to bdiv() will return either quotient (scalar context) 3394 # or quotient and remainder (list context). 3395 $x -> bdiv($b -> bpow($y), $a, $p, $r, $y); 3396} 3397 3398############################################################################### 3399# Bitwise methods 3400############################################################################### 3401 3402sub band { 3403 my $x = shift; 3404 my $xref = ref($x); 3405 my $class = $xref || $x; 3406 3407 croak 'band() is an instance method, not a class method' unless $xref; 3408 croak 'Not enough arguments for band()' if @_ < 1; 3409 3410 return if $x -> modify('band'); 3411 3412 my $y = shift; 3413 $y = $class -> new($y) unless ref($y); 3414 3415 my @r = @_; 3416 3417 my $xtmp = Math::BigInt -> new($x -> bint()); # to Math::BigInt 3418 $xtmp -> band($y); 3419 $xtmp = $class -> new($xtmp); # back to Math::BigFloat 3420 3421 $x -> {sign} = $xtmp -> {sign}; 3422 $x -> {_m} = $xtmp -> {_m}; 3423 $x -> {_es} = $xtmp -> {_es}; 3424 $x -> {_e} = $xtmp -> {_e}; 3425 3426 return $x -> round(@r); 3427} 3428 3429sub bior { 3430 my $x = shift; 3431 my $xref = ref($x); 3432 my $class = $xref || $x; 3433 3434 croak 'bior() is an instance method, not a class method' unless $xref; 3435 croak 'Not enough arguments for bior()' if @_ < 1; 3436 3437 return if $x -> modify('bior'); 3438 3439 my $y = shift; 3440 $y = $class -> new($y) unless ref($y); 3441 3442 my @r = @_; 3443 3444 my $xtmp = Math::BigInt -> new($x -> bint()); # to Math::BigInt 3445 $xtmp -> bior($y); 3446 $xtmp = $class -> new($xtmp); # back to Math::BigFloat 3447 3448 $x -> {sign} = $xtmp -> {sign}; 3449 $x -> {_m} = $xtmp -> {_m}; 3450 $x -> {_es} = $xtmp -> {_es}; 3451 $x -> {_e} = $xtmp -> {_e}; 3452 3453 return $x -> round(@r); 3454} 3455 3456sub bxor { 3457 my $x = shift; 3458 my $xref = ref($x); 3459 my $class = $xref || $x; 3460 3461 croak 'bxor() is an instance method, not a class method' unless $xref; 3462 croak 'Not enough arguments for bxor()' if @_ < 1; 3463 3464 return if $x -> modify('bxor'); 3465 3466 my $y = shift; 3467 $y = $class -> new($y) unless ref($y); 3468 3469 my @r = @_; 3470 3471 my $xtmp = Math::BigInt -> new($x -> bint()); # to Math::BigInt 3472 $xtmp -> bxor($y); 3473 $xtmp = $class -> new($xtmp); # back to Math::BigFloat 3474 3475 $x -> {sign} = $xtmp -> {sign}; 3476 $x -> {_m} = $xtmp -> {_m}; 3477 $x -> {_es} = $xtmp -> {_es}; 3478 $x -> {_e} = $xtmp -> {_e}; 3479 3480 return $x -> round(@r); 3481} 3482 3483sub bnot { 3484 my $x = shift; 3485 my $xref = ref($x); 3486 my $class = $xref || $x; 3487 3488 croak 'bnot() is an instance method, not a class method' unless $xref; 3489 3490 return if $x -> modify('bnot'); 3491 3492 my @r = @_; 3493 3494 my $xtmp = Math::BigInt -> new($x -> bint()); # to Math::BigInt 3495 $xtmp -> bnot(); 3496 $xtmp = $class -> new($xtmp); # back to Math::BigFloat 3497 3498 $x -> {sign} = $xtmp -> {sign}; 3499 $x -> {_m} = $xtmp -> {_m}; 3500 $x -> {_es} = $xtmp -> {_es}; 3501 $x -> {_e} = $xtmp -> {_e}; 3502 3503 return $x -> round(@r); 3504} 3505 3506############################################################################### 3507# Rounding methods 3508############################################################################### 3509 3510sub bround { 3511 # accuracy: preserve $N digits, and overwrite the rest with 0's 3512 my $x = shift; 3513 my $class = ref($x) || $x; 3514 $x = $class->new(shift) if !ref($x); 3515 3516 if (($_[0] || 0) < 0) { 3517 croak('bround() needs positive accuracy'); 3518 } 3519 3520 my ($scale, $mode) = $x->_scale_a(@_); 3521 return $x if !defined $scale || $x->modify('bround'); # no-op 3522 3523 # scale is now either $x->{_a}, $accuracy, or the user parameter 3524 # test whether $x already has lower accuracy, do nothing in this case 3525 # but do round if the accuracy is the same, since a math operation might 3526 # want to round a number with A=5 to 5 digits afterwards again 3527 return $x if defined $x->{_a} && $x->{_a} < $scale; 3528 3529 # scale < 0 makes no sense 3530 # scale == 0 => keep all digits 3531 # never round a +-inf, NaN 3532 return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/; 3533 3534 # 1: never round a 0 3535 # 2: if we should keep more digits than the mantissa has, do nothing 3536 if ($x->is_zero() || $LIB->_len($x->{_m}) <= $scale) { 3537 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale; 3538 return $x; 3539 } 3540 3541 # pass sign to bround for '+inf' and '-inf' rounding modes 3542 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt'; 3543 3544 $m->bround($scale, $mode); # round mantissa 3545 $x->{_m} = $m->{value}; # get our mantissa back 3546 $x->{_a} = $scale; # remember rounding 3547 delete $x->{_p}; # and clear P 3548 $x->bnorm(); # del trailing zeros gen. by bround() 3549} 3550 3551sub bfround { 3552 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.' 3553 # $n == 0 means round to integer 3554 # expects and returns normalized numbers! 3555 my $x = shift; 3556 my $class = ref($x) || $x; 3557 $x = $class->new(shift) if !ref($x); 3558 3559 my ($scale, $mode) = $x->_scale_p(@_); 3560 return $x if !defined $scale || $x->modify('bfround'); # no-op 3561 3562 # never round a 0, +-inf, NaN 3563 if ($x->is_zero()) { 3564 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2 3565 return $x; 3566 } 3567 return $x if $x->{sign} !~ /^[+-]$/; 3568 3569 # don't round if x already has lower precision 3570 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p}); 3571 3572 $x->{_p} = $scale; # remember round in any case 3573 delete $x->{_a}; # and clear A 3574 if ($scale < 0) { 3575 # round right from the '.' 3576 3577 return $x if $x->{_es} eq '+'; # e >= 0 => nothing to round 3578 3579 $scale = -$scale; # positive for simplicity 3580 my $len = $LIB->_len($x->{_m}); # length of mantissa 3581 3582 # the following poses a restriction on _e, but if _e is bigger than a 3583 # scalar, you got other problems (memory etc) anyway 3584 my $dad = -(0+ ($x->{_es}.$LIB->_num($x->{_e}))); # digits after dot 3585 my $zad = 0; # zeros after dot 3586 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style 3587 3588 # print "scale $scale dad $dad zad $zad len $len\n"; 3589 # number bsstr len zad dad 3590 # 0.123 123e-3 3 0 3 3591 # 0.0123 123e-4 3 1 4 3592 # 0.001 1e-3 1 2 3 3593 # 1.23 123e-2 3 0 2 3594 # 1.2345 12345e-4 5 0 4 3595 3596 # do not round after/right of the $dad 3597 return $x if $scale > $dad; # 0.123, scale >= 3 => exit 3598 3599 # round to zero if rounding inside the $zad, but not for last zero like: 3600 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case) 3601 return $x->bzero() if $scale < $zad; 3602 if ($scale == $zad) # for 0.006, scale -3 and trunc 3603 { 3604 $scale = -$len; 3605 } else { 3606 # adjust round-point to be inside mantissa 3607 if ($zad != 0) { 3608 $scale = $scale-$zad; 3609 } else { 3610 my $dbd = $len - $dad; 3611 $dbd = 0 if $dbd < 0; # digits before dot 3612 $scale = $dbd+$scale; 3613 } 3614 } 3615 } else { 3616 # round left from the '.' 3617 3618 # 123 => 100 means length(123) = 3 - $scale (2) => 1 3619 3620 my $dbt = $LIB->_len($x->{_m}); 3621 # digits before dot 3622 my $dbd = $dbt + ($x->{_es} . $LIB->_num($x->{_e})); 3623 # should be the same, so treat it as this 3624 $scale = 1 if $scale == 0; 3625 # shortcut if already integer 3626 return $x if $scale == 1 && $dbt <= $dbd; 3627 # maximum digits before dot 3628 ++$dbd; 3629 3630 if ($scale > $dbd) { 3631 # not enough digits before dot, so round to zero 3632 return $x->bzero; 3633 } elsif ($scale == $dbd) { 3634 # maximum 3635 $scale = -$dbt; 3636 } else { 3637 $scale = $dbd - $scale; 3638 } 3639 } 3640 # pass sign to bround for rounding modes '+inf' and '-inf' 3641 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt'; 3642 $m->bround($scale, $mode); 3643 $x->{_m} = $m->{value}; # get our mantissa back 3644 $x->bnorm(); 3645} 3646 3647sub bfloor { 3648 # round towards minus infinity 3649 my ($class, $x, $a, $p, $r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 3650 3651 return $x if $x->modify('bfloor'); 3652 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf 3653 3654 # if $x has digits after dot 3655 if ($x->{_es} eq '-') { 3656 $x->{_m} = $LIB->_rsft($x->{_m}, $x->{_e}, 10); # cut off digits after dot 3657 $x->{_e} = $LIB->_zero(); # trunc/norm 3658 $x->{_es} = '+'; # abs e 3659 $x->{_m} = $LIB->_inc($x->{_m}) if $x->{sign} eq '-'; # increment if negative 3660 } 3661 $x->round($a, $p, $r); 3662} 3663 3664sub bceil { 3665 # round towards plus infinity 3666 my ($class, $x, $a, $p, $r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 3667 3668 return $x if $x->modify('bceil'); 3669 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf 3670 3671 # if $x has digits after dot 3672 if ($x->{_es} eq '-') { 3673 $x->{_m} = $LIB->_rsft($x->{_m}, $x->{_e}, 10); # cut off digits after dot 3674 $x->{_e} = $LIB->_zero(); # trunc/norm 3675 $x->{_es} = '+'; # abs e 3676 if ($x->{sign} eq '+') { 3677 $x->{_m} = $LIB->_inc($x->{_m}); # increment if positive 3678 } else { 3679 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # avoid -0 3680 } 3681 } 3682 $x->round($a, $p, $r); 3683} 3684 3685sub bint { 3686 # round towards zero 3687 my ($class, $x, $a, $p, $r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 3688 3689 return $x if $x->modify('bint'); 3690 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf 3691 3692 # if $x has digits after the decimal point 3693 if ($x->{_es} eq '-') { 3694 $x->{_m} = $LIB->_rsft($x->{_m}, $x->{_e}, 10); # cut off digits after dot 3695 $x->{_e} = $LIB->_zero(); # truncate/normalize 3696 $x->{_es} = '+'; # abs e 3697 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # avoid -0 3698 } 3699 $x->round($a, $p, $r); 3700} 3701 3702############################################################################### 3703# Other mathematical methods 3704############################################################################### 3705 3706sub bgcd { 3707 # (BINT or num_str, BINT or num_str) return BINT 3708 # does not modify arguments, but returns new object 3709 3710 unshift @_, __PACKAGE__ 3711 unless ref($_[0]) || $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i; 3712 3713 my ($class, @args) = objectify(0, @_); 3714 3715 my $x = shift @args; 3716 $x = ref($x) && $x -> isa($class) ? $x -> copy() : $class -> new($x); 3717 return $class->bnan() unless $x -> is_int(); 3718 3719 while (@args) { 3720 my $y = shift @args; 3721 $y = $class->new($y) unless ref($y) && $y -> isa($class); 3722 return $class->bnan() unless $y -> is_int(); 3723 3724 # greatest common divisor 3725 while (! $y->is_zero()) { 3726 ($x, $y) = ($y->copy(), $x->copy()->bmod($y)); 3727 } 3728 3729 last if $x -> is_one(); 3730 } 3731 return $x -> babs(); 3732} 3733 3734sub blcm { 3735 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT 3736 # does not modify arguments, but returns new object 3737 # Least Common Multiple 3738 3739 unshift @_, __PACKAGE__ 3740 unless ref($_[0]) || $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i; 3741 3742 my ($class, @args) = objectify(0, @_); 3743 3744 my $x = shift @args; 3745 $x = ref($x) && $x -> isa($class) ? $x -> copy() : $class -> new($x); 3746 return $class->bnan() if $x->{sign} !~ /^[+-]$/; # x NaN? 3747 3748 while (@args) { 3749 my $y = shift @args; 3750 $y = $class -> new($y) unless ref($y) && $y -> isa($class); 3751 return $x->bnan() unless $y -> is_int(); 3752 my $gcd = $x -> bgcd($y); 3753 $x -> bdiv($gcd) -> bmul($y); 3754 } 3755 3756 return $x -> babs(); 3757} 3758 3759############################################################################### 3760# Object property methods 3761############################################################################### 3762 3763sub length { 3764 my $x = shift; 3765 my $class = ref($x) || $x; 3766 $x = $class->new(shift) unless ref($x); 3767 3768 return 1 if $LIB->_is_zero($x->{_m}); 3769 3770 my $len = $LIB->_len($x->{_m}); 3771 $len += $LIB->_num($x->{_e}) if $x->{_es} eq '+'; 3772 if (wantarray()) { 3773 my $t = 0; 3774 $t = $LIB->_num($x->{_e}) if $x->{_es} eq '-'; 3775 return ($len, $t); 3776 } 3777 $len; 3778} 3779 3780sub mantissa { 3781 # return a copy of the mantissa 3782 my ($class, $x) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 3783 3784 if ($x->{sign} !~ /^[+-]$/) { 3785 my $s = $x->{sign}; 3786 $s =~ s/^[+]//; 3787 return Math::BigInt->new($s, undef, undef); # -inf, +inf => +inf 3788 } 3789 my $m = Math::BigInt->new($LIB->_str($x->{_m}), undef, undef); 3790 $m->bneg() if $x->{sign} eq '-'; 3791 3792 $m; 3793} 3794 3795sub exponent { 3796 # return a copy of the exponent 3797 my ($class, $x) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 3798 3799 if ($x->{sign} !~ /^[+-]$/) { 3800 my $s = $x->{sign}; 3801$s =~ s/^[+-]//; 3802 return Math::BigInt->new($s, undef, undef); # -inf, +inf => +inf 3803 } 3804 Math::BigInt->new($x->{_es} . $LIB->_str($x->{_e}), undef, undef); 3805} 3806 3807sub parts { 3808 # return a copy of both the exponent and the mantissa 3809 my ($class, $x) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 3810 3811 if ($x->{sign} !~ /^[+-]$/) { 3812 my $s = $x->{sign}; 3813$s =~ s/^[+]//; 3814my $se = $s; 3815$se =~ s/^[-]//; 3816 return ($class->new($s), $class->new($se)); # +inf => inf and -inf, +inf => inf 3817 } 3818 my $m = Math::BigInt->bzero(); 3819 $m->{value} = $LIB->_copy($x->{_m}); 3820 $m->bneg() if $x->{sign} eq '-'; 3821 ($m, Math::BigInt->new($x->{_es} . $LIB->_num($x->{_e}))); 3822} 3823 3824sub sparts { 3825 my $self = shift; 3826 my $class = ref $self; 3827 3828 croak("sparts() is an instance method, not a class method") 3829 unless $class; 3830 3831 # Not-a-number. 3832 3833 if ($self -> is_nan()) { 3834 my $mant = $self -> copy(); # mantissa 3835 return $mant unless wantarray; # scalar context 3836 my $expo = $class -> bnan(); # exponent 3837 return ($mant, $expo); # list context 3838 } 3839 3840 # Infinity. 3841 3842 if ($self -> is_inf()) { 3843 my $mant = $self -> copy(); # mantissa 3844 return $mant unless wantarray; # scalar context 3845 my $expo = $class -> binf('+'); # exponent 3846 return ($mant, $expo); # list context 3847 } 3848 3849 # Finite number. 3850 3851 my $mant = $class -> bzero(); 3852 $mant -> {sign} = $self -> {sign}; 3853 $mant -> {_m} = $LIB->_copy($self -> {_m}); 3854 return $mant unless wantarray; 3855 3856 my $expo = $class -> bzero(); 3857 $expo -> {sign} = $self -> {_es}; 3858 $expo -> {_m} = $LIB->_copy($self -> {_e}); 3859 3860 return ($mant, $expo); 3861} 3862 3863sub nparts { 3864 my $self = shift; 3865 my $class = ref $self; 3866 3867 croak("nparts() is an instance method, not a class method") 3868 unless $class; 3869 3870 # Not-a-number. 3871 3872 if ($self -> is_nan()) { 3873 my $mant = $self -> copy(); # mantissa 3874 return $mant unless wantarray; # scalar context 3875 my $expo = $class -> bnan(); # exponent 3876 return ($mant, $expo); # list context 3877 } 3878 3879 # Infinity. 3880 3881 if ($self -> is_inf()) { 3882 my $mant = $self -> copy(); # mantissa 3883 return $mant unless wantarray; # scalar context 3884 my $expo = $class -> binf('+'); # exponent 3885 return ($mant, $expo); # list context 3886 } 3887 3888 # Finite number. 3889 3890 my ($mant, $expo) = $self -> sparts(); 3891 3892 if ($mant -> bcmp(0)) { 3893 my ($ndigtot, $ndigfrac) = $mant -> length(); 3894 my $expo10adj = $ndigtot - $ndigfrac - 1; 3895 3896 if ($expo10adj != 0) { 3897 my $factor = "1e" . -$expo10adj; 3898 $mant -> bmul($factor); 3899 return $mant unless wantarray; 3900 $expo -> badd($expo10adj); 3901 return ($mant, $expo); 3902 } 3903 } 3904 3905 return $mant unless wantarray; 3906 return ($mant, $expo); 3907} 3908 3909sub eparts { 3910 my $self = shift; 3911 my $class = ref $self; 3912 3913 croak("eparts() is an instance method, not a class method") 3914 unless $class; 3915 3916 # Not-a-number and Infinity. 3917 3918 return $self -> sparts() if $self -> is_nan() || $self -> is_inf(); 3919 3920 # Finite number. 3921 3922 my ($mant, $expo) = $self -> nparts(); 3923 3924 my $c = $expo -> copy() -> bmod(3); 3925 $mant -> blsft($c, 10); 3926 return $mant unless wantarray; 3927 3928 $expo -> bsub($c); 3929 return ($mant, $expo); 3930} 3931 3932sub dparts { 3933 my $self = shift; 3934 my $class = ref $self; 3935 3936 croak("dparts() is an instance method, not a class method") 3937 unless $class; 3938 3939 # Not-a-number and Infinity. 3940 3941 if ($self -> is_nan() || $self -> is_inf()) { 3942 my $int = $self -> copy(); 3943 return $int unless wantarray; 3944 my $frc = $class -> bzero(); 3945 return ($int, $frc); 3946 } 3947 3948 my $int = $self -> copy(); 3949 my $frc = $class -> bzero(); 3950 3951 # If the input has a fraction part. 3952 3953 if ($int->{_es} eq '-') { 3954 $int->{_m} = $LIB -> _rsft($int->{_m}, $int->{_e}, 10); 3955 $int->{_e} = $LIB -> _zero(); 3956 $int->{_es} = '+'; 3957 $int->{sign} = '+' if $LIB->_is_zero($int->{_m}); # avoid -0 3958 3959 return $int unless wantarray; 3960 $frc = $self -> copy() -> bsub($int); 3961 return ($int, $frc); 3962 } 3963 3964 return $int unless wantarray; 3965 return ($int, $frc); 3966} 3967 3968############################################################################### 3969# String conversion methods 3970############################################################################### 3971 3972sub bstr { 3973 # (ref to BFLOAT or num_str) return num_str 3974 # Convert number from internal format to (non-scientific) string format. 3975 # internal format is always normalized (no leading zeros, "-0" => "+0") 3976 my ($class, $x) = ref($_[0]) ? (undef, $_[0]) : objectify(1, @_); 3977 3978 if ($x->{sign} !~ /^[+-]$/) { 3979 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 3980 return 'inf'; # +inf 3981 } 3982 3983 my $es = '0'; 3984my $len = 1; 3985my $cad = 0; 3986my $dot = '.'; 3987 3988 # $x is zero? 3989 my $not_zero = !($x->{sign} eq '+' && $LIB->_is_zero($x->{_m})); 3990 if ($not_zero) { 3991 $es = $LIB->_str($x->{_m}); 3992 $len = CORE::length($es); 3993 my $e = $LIB->_num($x->{_e}); 3994 $e = -$e if $x->{_es} eq '-'; 3995 if ($e < 0) { 3996 $dot = ''; 3997 # if _e is bigger than a scalar, the following will blow your memory 3998 if ($e <= -$len) { 3999 my $r = abs($e) - $len; 4000 $es = '0.'. ('0' x $r) . $es; 4001 $cad = -($len+$r); 4002 } else { 4003 substr($es, $e, 0) = '.'; 4004 $cad = $LIB->_num($x->{_e}); 4005 $cad = -$cad if $x->{_es} eq '-'; 4006 } 4007 } elsif ($e > 0) { 4008 # expand with zeros 4009 $es .= '0' x $e; 4010$len += $e; 4011$cad = 0; 4012 } 4013 } # if not zero 4014 4015 $es = '-'.$es if $x->{sign} eq '-'; 4016 # if set accuracy or precision, pad with zeros on the right side 4017 if ((defined $x->{_a}) && ($not_zero)) { 4018 # 123400 => 6, 0.1234 => 4, 0.001234 => 4 4019 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340 4020 $zeros = $x->{_a} - $len if $cad != $len; 4021 $es .= $dot.'0' x $zeros if $zeros > 0; 4022 } elsif ((($x->{_p} || 0) < 0)) { 4023 # 123400 => 6, 0.1234 => 4, 0.001234 => 6 4024 my $zeros = -$x->{_p} + $cad; 4025 $es .= $dot.'0' x $zeros if $zeros > 0; 4026 } 4027 $es; 4028} 4029 4030# Decimal notation, e.g., "12345.6789". 4031 4032sub bdstr { 4033 my $x = shift; 4034 4035 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 4036 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 4037 return 'inf'; # +inf 4038 } 4039 4040 my $mant = $LIB->_str($x->{_m}); 4041 my $expo = $x -> exponent(); 4042 4043 my $str = $mant; 4044 if ($expo >= 0) { 4045 $str .= "0" x $expo; 4046 } else { 4047 my $mantlen = CORE::length($mant); 4048 my $c = $mantlen + $expo; 4049 $str = "0" x (1 - $c) . $str if $c <= 0; 4050 substr($str, $expo, 0) = '.'; 4051 } 4052 4053 return $x->{sign} eq '-' ? "-$str" : $str; 4054} 4055 4056# Scientific notation with significand/mantissa as an integer, e.g., "12345.6789" 4057# is written as "123456789e-4". 4058 4059sub bsstr { 4060 my ($class, $x) = ref($_[0]) ? (undef, $_[0]) : objectify(1, @_); 4061 4062 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 4063 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 4064 return 'inf'; # +inf 4065 } 4066 4067 my $str = $LIB->_str($x->{_m}) . 'e' . $x->{_es}. $LIB->_str($x->{_e}); 4068 return $x->{sign} eq '-' ? "-$str" : $str; 4069} 4070 4071# Normalized notation, e.g., "12345.6789" is written as "1.23456789e+4". 4072 4073sub bnstr { 4074 my $x = shift; 4075 4076 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 4077 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 4078 return 'inf'; # +inf 4079 } 4080 4081 my ($mant, $expo) = $x -> nparts(); 4082 4083 my $esgn = $expo < 0 ? '-' : '+'; 4084 my $eabs = $expo -> babs() -> bfround(0) -> bstr(); 4085 #$eabs = '0' . $eabs if length($eabs) < 2; 4086 4087 return $mant . 'e' . $esgn . $eabs; 4088} 4089 4090# Engineering notation, e.g., "12345.6789" is written as "12.3456789e+3". 4091 4092sub bestr { 4093 my $x = shift; 4094 4095 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 4096 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 4097 return 'inf'; # +inf 4098 } 4099 4100 my ($mant, $expo) = $x -> eparts(); 4101 4102 my $esgn = $expo < 0 ? '-' : '+'; 4103 my $eabs = $expo -> babs() -> bfround(0) -> bstr(); 4104 #$eabs = '0' . $eabs if length($eabs) < 2; 4105 4106 return $mant . 'e' . $esgn . $eabs; 4107} 4108 4109sub to_hex { 4110 # return number as hexadecimal string (only for integers defined) 4111 4112 my ($class, $x) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 4113 4114 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc 4115 return '0' if $x->is_zero(); 4116 4117 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex? 4118 4119 my $z = $LIB->_copy($x->{_m}); 4120 if (! $LIB->_is_zero($x->{_e})) { # > 0 4121 $z = $LIB->_lsft($z, $x->{_e}, 10); 4122 } 4123 my $str = $LIB->_to_hex($z); 4124 return $x->{sign} eq '-' ? "-$str" : $str; 4125} 4126 4127sub to_oct { 4128 # return number as octal digit string (only for integers defined) 4129 4130 my ($class, $x) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 4131 4132 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc 4133 return '0' if $x->is_zero(); 4134 4135 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in octal? 4136 4137 my $z = $LIB->_copy($x->{_m}); 4138 if (! $LIB->_is_zero($x->{_e})) { # > 0 4139 $z = $LIB->_lsft($z, $x->{_e}, 10); 4140 } 4141 my $str = $LIB->_to_oct($z); 4142 return $x->{sign} eq '-' ? "-$str" : $str; 4143} 4144 4145sub to_bin { 4146 # return number as binary digit string (only for integers defined) 4147 4148 my ($class, $x) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 4149 4150 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc 4151 return '0' if $x->is_zero(); 4152 4153 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in binary? 4154 4155 my $z = $LIB->_copy($x->{_m}); 4156 if (! $LIB->_is_zero($x->{_e})) { # > 0 4157 $z = $LIB->_lsft($z, $x->{_e}, 10); 4158 } 4159 my $str = $LIB->_to_bin($z); 4160 return $x->{sign} eq '-' ? "-$str" : $str; 4161} 4162 4163sub as_hex { 4164 # return number as hexadecimal string (only for integers defined) 4165 4166 my ($class, $x) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 4167 4168 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc 4169 return '0x0' if $x->is_zero(); 4170 4171 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex? 4172 4173 my $z = $LIB->_copy($x->{_m}); 4174 if (! $LIB->_is_zero($x->{_e})) { # > 0 4175 $z = $LIB->_lsft($z, $x->{_e}, 10); 4176 } 4177 my $str = $LIB->_as_hex($z); 4178 return $x->{sign} eq '-' ? "-$str" : $str; 4179} 4180 4181sub as_oct { 4182 # return number as octal digit string (only for integers defined) 4183 4184 my ($class, $x) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 4185 4186 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc 4187 return '00' if $x->is_zero(); 4188 4189 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in octal? 4190 4191 my $z = $LIB->_copy($x->{_m}); 4192 if (! $LIB->_is_zero($x->{_e})) { # > 0 4193 $z = $LIB->_lsft($z, $x->{_e}, 10); 4194 } 4195 my $str = $LIB->_as_oct($z); 4196 return $x->{sign} eq '-' ? "-$str" : $str; 4197} 4198 4199sub as_bin { 4200 # return number as binary digit string (only for integers defined) 4201 4202 my ($class, $x) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 4203 4204 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc 4205 return '0b0' if $x->is_zero(); 4206 4207 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in binary? 4208 4209 my $z = $LIB->_copy($x->{_m}); 4210 if (! $LIB->_is_zero($x->{_e})) { # > 0 4211 $z = $LIB->_lsft($z, $x->{_e}, 10); 4212 } 4213 my $str = $LIB->_as_bin($z); 4214 return $x->{sign} eq '-' ? "-$str" : $str; 4215} 4216 4217sub numify { 4218 # Make a Perl scalar number from a Math::BigFloat object. 4219 my ($class, $x) = ref($_[0]) ? (undef, $_[0]) : objectify(1, @_); 4220 4221 if ($x -> is_nan()) { 4222 require Math::Complex; 4223 my $inf = Math::Complex::Inf(); 4224 return $inf - $inf; 4225 } 4226 4227 if ($x -> is_inf()) { 4228 require Math::Complex; 4229 my $inf = Math::Complex::Inf(); 4230 return $x -> is_negative() ? -$inf : $inf; 4231 } 4232 4233 # Create a string and let Perl's atoi()/atof() handle the rest. 4234 return 0 + $x -> bsstr(); 4235} 4236 4237############################################################################### 4238# Private methods and functions. 4239############################################################################### 4240 4241sub import { 4242 my $class = shift; 4243 my $l = scalar @_; 4244 my $lib = ''; 4245my @a; 4246 my $lib_kind = 'try'; 4247 $IMPORT=1; 4248 for (my $i = 0; $i < $l ; $i++) { 4249 if ($_[$i] eq ':constant') { 4250 # This causes overlord er load to step in. 'binary' and 'integer' 4251 # are handled by BigInt. 4252 overload::constant float => sub { $class->new(shift); }; 4253 } elsif ($_[$i] eq 'upgrade') { 4254 # this causes upgrading 4255 $upgrade = $_[$i+1]; # or undef to disable 4256 $i++; 4257 } elsif ($_[$i] eq 'downgrade') { 4258 # this causes downgrading 4259 $downgrade = $_[$i+1]; # or undef to disable 4260 $i++; 4261 } elsif ($_[$i] =~ /^(lib|try|only)\z/) { 4262 # alternative library 4263 $lib = $_[$i+1] || ''; # default Calc 4264 $lib_kind = $1; # lib, try or only 4265 $i++; 4266 } elsif ($_[$i] eq 'with') { 4267 # alternative class for our private parts() 4268 # XXX: no longer supported 4269 # $LIB = $_[$i+1] || 'Math::BigInt'; 4270 $i++; 4271 } else { 4272 push @a, $_[$i]; 4273 } 4274 } 4275 4276 $lib =~ tr/a-zA-Z0-9,://cd; # restrict to sane characters 4277 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work 4278 my $mbilib = eval { Math::BigInt->config('lib') }; 4279 if ((defined $mbilib) && ($LIB eq 'Math::BigInt::Calc')) { 4280 # $LIB already loaded 4281 Math::BigInt->import($lib_kind, "$lib, $mbilib", 'objectify'); 4282 } else { 4283 # $LIB not loaded, or with ne "Math::BigInt::Calc" 4284 $lib .= ",$mbilib" if defined $mbilib; 4285 $lib =~ s/^,//; # don't leave empty 4286 4287 # replacement library can handle lib statement, but also could ignore it 4288 4289 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is 4290 # used in the same script, or eval inside import(). So we require MBI: 4291 require Math::BigInt; 4292 Math::BigInt->import($lib_kind => $lib, 'objectify'); 4293 } 4294 if ($@) { 4295 croak("Couldn't load $lib: $! $@"); 4296 } 4297 # find out which one was actually loaded 4298 $LIB = Math::BigInt->config('lib'); 4299 4300 # register us with MBI to get notified of future lib changes 4301 Math::BigInt::_register_callback($class, sub { $LIB = $_[0]; }); 4302 4303 $class->export_to_level(1, $class, @a); # export wanted functions 4304} 4305 4306sub _len_to_steps { 4307 # Given D (digits in decimal), compute N so that N! (N factorial) is 4308 # at least D digits long. D should be at least 50. 4309 my $d = shift; 4310 4311 # two constants for the Ramanujan estimate of ln(N!) 4312 my $lg2 = log(2 * 3.14159265) / 2; 4313 my $lg10 = log(10); 4314 4315 # D = 50 => N => 42, so L = 40 and R = 50 4316 my $l = 40; 4317my $r = $d; 4318 4319 # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :( 4320 $l = $l->numify if ref($l); 4321 $r = $r->numify if ref($r); 4322 $lg2 = $lg2->numify if ref($lg2); 4323 $lg10 = $lg10->numify if ref($lg10); 4324 4325 # binary search for the right value (could this be written as the reverse of lg(n!)?) 4326 while ($r - $l > 1) { 4327 my $n = int(($r - $l) / 2) + $l; 4328 my $ramanujan = 4329 int(($n * log($n) - $n + log($n * (1 + 4*$n*(1+2*$n))) / 6 + $lg2) / $lg10); 4330 $ramanujan > $d ? $r = $n : $l = $n; 4331 } 4332 $l; 4333} 4334 4335sub _log { 4336 # internal log function to calculate ln() based on Taylor series. 4337 # Modifies $x in place. 4338 my ($x, $scale) = @_; 4339 my $class = ref $x; 4340 4341 # in case of $x == 1, result is 0 4342 return $x->bzero() if $x->is_one(); 4343 4344 # XXX TODO: rewrite this in a similar manner to bexp() 4345 4346 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log 4347 4348 # u = x-1, v = x+1 4349 # _ _ 4350 # Taylor: | u 1 u^3 1 u^5 | 4351 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0 4352 # |_ v 3 v^3 5 v^5 _| 4353 4354 # This takes much more steps to calculate the result and is thus not used 4355 # u = x-1 4356 # _ _ 4357 # Taylor: | u 1 u^2 1 u^3 | 4358 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2 4359 # |_ x 2 x^2 3 x^3 _| 4360 4361 my ($limit, $v, $u, $below, $factor, $next, $over, $f); 4362 4363 $v = $x->copy(); $v->binc(); # v = x+1 4364 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1 4365 $x->bdiv($v, $scale); # first term: u/v 4366 $below = $v->copy(); 4367 $over = $u->copy(); 4368 $u *= $u; $v *= $v; # u^2, v^2 4369 $below->bmul($v); # u^3, v^3 4370 $over->bmul($u); 4371 $factor = $class->new(3); $f = $class->new(2); 4372 4373 my $steps = 0; 4374 $limit = $class->new("1E-". ($scale-1)); 4375 4376 while (3 < 5) { 4377 # we calculate the next term, and add it to the last 4378 # when the next term is below our limit, it won't affect the outcome 4379 # anymore, so we stop 4380 4381 # calculating the next term simple from over/below will result in quite 4382 # a time hog if the input has many digits, since over and below will 4383 # accumulate more and more digits, and the result will also have many 4384 # digits, but in the end it is rounded to $scale digits anyway. So if we 4385 # round $over and $below first, we save a lot of time for the division 4386 # (not with log(1.2345), but try log (123**123) to see what I mean. This 4387 # can introduce a rounding error if the division result would be f.i. 4388 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but 4389 # if we truncated $over and $below we might get 0.12345. Does this matter 4390 # for the end result? So we give $over and $below 4 more digits to be 4391 # on the safe side (unscientific error handling as usual... :+D 4392 4393 $next = $over->copy()->bround($scale+4) 4394 ->bdiv($below->copy()->bmul($factor)->bround($scale+4), 4395 $scale); 4396 4397 ## old version: 4398 ## $next = $over->copy()->bdiv($below->copy()->bmul($factor), $scale); 4399 4400 last if $next->bacmp($limit) <= 0; 4401 4402 delete $next->{_a}; 4403 delete $next->{_p}; 4404 $x->badd($next); 4405 # calculate things for the next term 4406 $over *= $u; 4407 $below *= $v; 4408 $factor->badd($f); 4409 } 4410 $x->bmul($f); # $x *= 2 4411} 4412 4413sub _log_10 { 4414 # Internal log function based on reducing input to the range of 0.1 .. 9.99 4415 # and then "correcting" the result to the proper one. Modifies $x in place. 4416 my ($x, $scale) = @_; 4417 my $class = ref $x; 4418 4419 # Taking blog() from numbers greater than 10 takes a *very long* time, so we 4420 # break the computation down into parts based on the observation that: 4421 # blog(X*Y) = blog(X) + blog(Y) 4422 # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller 4423 # $x is the faster it gets. Since 2*$x takes about 10 times as 4424 # long, we make it faster by about a factor of 100 by dividing $x by 10. 4425 4426 # The same observation is valid for numbers smaller than 0.1, e.g. computing 4427 # log(1) is fastest, and the further away we get from 1, the longer it takes. 4428 # So we also 'break' this down by multiplying $x with 10 and subtract the 4429 # log(10) afterwards to get the correct result. 4430 4431 # To get $x even closer to 1, we also divide by 2 and then use log(2) to 4432 # correct for this. For instance if $x is 2.4, we use the formula: 4433 # blog(2.4 * 2) == blog(1.2) + blog(2) 4434 # and thus calculate only blog(1.2) and blog(2), which is faster in total 4435 # than calculating blog(2.4). 4436 4437 # In addition, the values for blog(2) and blog(10) are cached. 4438 4439 # Calculate nr of digits before dot. x = 123, dbd = 3; x = 1.23, dbd = 1; 4440 # x = 0.0123, dbd = -1; x = 0.000123, dbd = -3, etc. 4441 4442 my $dbd = $LIB->_num($x->{_e}); 4443 $dbd = -$dbd if $x->{_es} eq '-'; 4444 $dbd += $LIB->_len($x->{_m}); 4445 4446 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid 4447 # infinite recursion 4448 4449 my $calc = 1; # do some calculation? 4450 4451 # disable the shortcut for 10, since we need log(10) and this would recurse 4452 # infinitely deep 4453 if ($x->{_es} eq '+' && # $x == 10 4454 ($LIB->_is_one($x->{_e}) && 4455 $LIB->_is_one($x->{_m}))) 4456 { 4457 $dbd = 0; # disable shortcut 4458 # we can use the cached value in these cases 4459 if ($scale <= $LOG_10_A) { 4460 $x->bzero(); 4461 $x->badd($LOG_10); # modify $x in place 4462 $calc = 0; # no need to calc, but round 4463 } 4464 # if we can't use the shortcut, we continue normally 4465 } else { 4466 # disable the shortcut for 2, since we maybe have it cached 4467 if (($LIB->_is_zero($x->{_e}) && # $x == 2 4468 $LIB->_is_two($x->{_m}))) 4469 { 4470 $dbd = 0; # disable shortcut 4471 # we can use the cached value in these cases 4472 if ($scale <= $LOG_2_A) { 4473 $x->bzero(); 4474 $x->badd($LOG_2); # modify $x in place 4475 $calc = 0; # no need to calc, but round 4476 } 4477 # if we can't use the shortcut, we continue normally 4478 } 4479 } 4480 4481 # if $x = 0.1, we know the result must be 0-log(10) 4482 if ($calc != 0 && 4483 ($x->{_es} eq '-' && # $x == 0.1 4484 ($LIB->_is_one($x->{_e}) && 4485 $LIB->_is_one($x->{_m})))) 4486 { 4487 $dbd = 0; # disable shortcut 4488 # we can use the cached value in these cases 4489 if ($scale <= $LOG_10_A) { 4490 $x->bzero(); 4491 $x->bsub($LOG_10); 4492 $calc = 0; # no need to calc, but round 4493 } 4494 } 4495 4496 return $x if $calc == 0; # already have the result 4497 4498 # default: these correction factors are undef and thus not used 4499 my $l_10; # value of ln(10) to A of $scale 4500 my $l_2; # value of ln(2) to A of $scale 4501 4502 my $two = $class->new(2); 4503 4504 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1 4505 # so don't do this shortcut for 1 or 0 4506 if (($dbd > 1) || ($dbd < 0)) { 4507 # convert our cached value to an object if not already (avoid doing this 4508 # at import() time, since not everybody needs this) 4509 $LOG_10 = $class->new($LOG_10, undef, undef) unless ref $LOG_10; 4510 4511 #print "x = $x, dbd = $dbd, calc = $calc\n"; 4512 # got more than one digit before the dot, or more than one zero after the 4513 # dot, so do: 4514 # log(123) == log(1.23) + log(10) * 2 4515 # log(0.0123) == log(1.23) - log(10) * 2 4516 4517 if ($scale <= $LOG_10_A) { 4518 # use cached value 4519 $l_10 = $LOG_10->copy(); # copy for mul 4520 } else { 4521 # else: slower, compute and cache result 4522 # also disable downgrade for this code path 4523 local $Math::BigFloat::downgrade = undef; 4524 4525 # shorten the time to calculate log(10) based on the following: 4526 # log(1.25 * 8) = log(1.25) + log(8) 4527 # = log(1.25) + log(2) + log(2) + log(2) 4528 4529 # first get $l_2 (and possible compute and cache log(2)) 4530 $LOG_2 = $class->new($LOG_2, undef, undef) unless ref $LOG_2; 4531 if ($scale <= $LOG_2_A) { 4532 # use cached value 4533 $l_2 = $LOG_2->copy(); # copy() for the mul below 4534 } else { 4535 # else: slower, compute and cache result 4536 $l_2 = $two->copy(); 4537 $l_2->_log($scale); # scale+4, actually 4538 $LOG_2 = $l_2->copy(); # cache the result for later 4539 # the copy() is for mul below 4540 $LOG_2_A = $scale; 4541 } 4542 4543 # now calculate log(1.25): 4544 $l_10 = $class->new('1.25'); 4545 $l_10->_log($scale); # scale+4, actually 4546 4547 # log(1.25) + log(2) + log(2) + log(2): 4548 $l_10->badd($l_2); 4549 $l_10->badd($l_2); 4550 $l_10->badd($l_2); 4551 $LOG_10 = $l_10->copy(); # cache the result for later 4552 # the copy() is for mul below 4553 $LOG_10_A = $scale; 4554 } 4555 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1 4556 $l_10->bmul($class->new($dbd)); # log(10) * (digits_before_dot-1) 4557 my $dbd_sign = '+'; 4558 if ($dbd < 0) { 4559 $dbd = -$dbd; 4560 $dbd_sign = '-'; 4561 } 4562 ($x->{_e}, $x->{_es}) = 4563 _e_sub($x->{_e}, $LIB->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23 4564 4565 } 4566 4567 # Now: 0.1 <= $x < 10 (and possible correction in l_10) 4568 4569 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div 4570 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1) 4571 4572 $HALF = $class->new($HALF) unless ref($HALF); 4573 4574 my $twos = 0; # default: none (0 times) 4575 while ($x->bacmp($HALF) <= 0) { # X <= 0.5 4576 $twos--; 4577 $x->bmul($two); 4578 } 4579 while ($x->bacmp($two) >= 0) { # X >= 2 4580 $twos++; 4581 $x->bdiv($two, $scale+4); # keep all digits 4582 } 4583 $x->bround($scale+4); 4584 # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both) 4585 # So calculate correction factor based on ln(2): 4586 if ($twos != 0) { 4587 $LOG_2 = $class->new($LOG_2, undef, undef) unless ref $LOG_2; 4588 if ($scale <= $LOG_2_A) { 4589 # use cached value 4590 $l_2 = $LOG_2->copy(); # copy() for the mul below 4591 } else { 4592 # else: slower, compute and cache result 4593 # also disable downgrade for this code path 4594 local $Math::BigFloat::downgrade = undef; 4595 $l_2 = $two->copy(); 4596 $l_2->_log($scale); # scale+4, actually 4597 $LOG_2 = $l_2->copy(); # cache the result for later 4598 # the copy() is for mul below 4599 $LOG_2_A = $scale; 4600 } 4601 $l_2->bmul($twos); # * -2 => subtract, * 2 => add 4602 } else { 4603 undef $l_2; 4604 } 4605 4606 $x->_log($scale); # need to do the "normal" way 4607 $x->badd($l_10) if defined $l_10; # correct it by ln(10) 4608 $x->badd($l_2) if defined $l_2; # and maybe by ln(2) 4609 4610 # all done, $x contains now the result 4611 $x; 4612} 4613 4614sub _e_add { 4615 # Internal helper sub to take two positive integers and their signs and 4616 # then add them. Input ($LIB, $LIB, ('+'|'-'), ('+'|'-')), output 4617 # ($LIB, ('+'|'-')). 4618 4619 my ($x, $y, $xs, $ys) = @_; 4620 4621 # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8) 4622 if ($xs eq $ys) { 4623 $x = $LIB->_add($x, $y); # +a + +b or -a + -b 4624 } else { 4625 my $a = $LIB->_acmp($x, $y); 4626 if ($a == 0) { 4627 # This does NOT modify $x in-place. TODO: Fix this? 4628 $x = $LIB->_zero(); # result is 0 4629 $xs = '+'; 4630 return ($x, $xs); 4631 } 4632 if ($a > 0) { 4633 $x = $LIB->_sub($x, $y); # abs sub 4634 } else { # a < 0 4635 $x = $LIB->_sub ($y, $x, 1); # abs sub 4636 $xs = $ys; 4637 } 4638 } 4639 4640 $xs = '+' if $xs eq '-' && $LIB->_is_zero($x); # no "-0" 4641 4642 return ($x, $xs); 4643} 4644 4645sub _e_sub { 4646 # Internal helper sub to take two positive integers and their signs and 4647 # then subtract them. Input ($LIB, $LIB, ('+'|'-'), ('+'|'-')), 4648 # output ($LIB, ('+'|'-')) 4649 my ($x, $y, $xs, $ys) = @_; 4650 4651 # flip sign 4652 $ys = $ys eq '+' ? '-' : '+'; # swap sign of second operand ... 4653 _e_add($x, $y, $xs, $ys); # ... and let _e_add() do the job 4654} 4655 4656sub _pow { 4657 # Calculate a power where $y is a non-integer, like 2 ** 0.3 4658 my ($x, $y, @r) = @_; 4659 my $class = ref($x); 4660 4661 # if $y == 0.5, it is sqrt($x) 4662 $HALF = $class->new($HALF) unless ref($HALF); 4663 return $x->bsqrt(@r, $y) if $y->bcmp($HALF) == 0; 4664 4665 # Using: 4666 # a ** x == e ** (x * ln a) 4667 4668 # u = y * ln x 4669 # _ _ 4670 # Taylor: | u u^2 u^3 | 4671 # x ** y = 1 + | --- + --- + ----- + ... | 4672 # |_ 1 1*2 1*2*3 _| 4673 4674 # we need to limit the accuracy to protect against overflow 4675 my $fallback = 0; 4676 my ($scale, @params); 4677 ($x, @params) = $x->_find_round_parameters(@r); 4678 4679 return $x if $x->is_nan(); # error in _find_round_parameters? 4680 4681 # no rounding at all, so must use fallback 4682 if (scalar @params == 0) { 4683 # simulate old behaviour 4684 $params[0] = $class->div_scale(); # and round to it as accuracy 4685 $params[1] = undef; # disable P 4686 $scale = $params[0]+4; # at least four more for proper round 4687 $params[2] = $r[2]; # round mode by caller or undef 4688 $fallback = 1; # to clear a/p afterwards 4689 } else { 4690 # the 4 below is empirical, and there might be cases where it is not 4691 # enough... 4692 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 4693 } 4694 4695 # when user set globals, they would interfere with our calculation, so 4696 # disable them and later re-enable them 4697 no strict 'refs'; 4698 my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef; 4699 my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef; 4700 # we also need to disable any set A or P on $x (_find_round_parameters took 4701 # them already into account), since these would interfere, too 4702 delete $x->{_a}; 4703 delete $x->{_p}; 4704 # need to disable $upgrade in BigInt, to avoid deep recursion 4705 local $Math::BigInt::upgrade = undef; 4706 4707 my ($limit, $v, $u, $below, $factor, $next, $over); 4708 4709 $u = $x->copy()->blog(undef, $scale)->bmul($y); 4710 my $do_invert = ($u->{sign} eq '-'); 4711 $u->bneg() if $do_invert; 4712 $v = $class->bone(); # 1 4713 $factor = $class->new(2); # 2 4714 $x->bone(); # first term: 1 4715 4716 $below = $v->copy(); 4717 $over = $u->copy(); 4718 4719 $limit = $class->new("1E-". ($scale-1)); 4720 #my $steps = 0; 4721 while (3 < 5) { 4722 # we calculate the next term, and add it to the last 4723 # when the next term is below our limit, it won't affect the outcome 4724 # anymore, so we stop: 4725 $next = $over->copy()->bdiv($below, $scale); 4726 last if $next->bacmp($limit) <= 0; 4727 $x->badd($next); 4728 # calculate things for the next term 4729 $over *= $u; 4730 $below *= $factor; 4731 $factor->binc(); 4732 4733 last if $x->{sign} !~ /^[-+]$/; 4734 4735 #$steps++; 4736 } 4737 4738 if ($do_invert) { 4739 my $x_copy = $x->copy(); 4740 $x->bone->bdiv($x_copy, $scale); 4741 } 4742 4743 # shortcut to not run through _find_round_parameters again 4744 if (defined $params[0]) { 4745 $x->bround($params[0], $params[2]); # then round accordingly 4746 } else { 4747 $x->bfround($params[1], $params[2]); # then round accordingly 4748 } 4749 if ($fallback) { 4750 # clear a/p after round, since user did not request it 4751 delete $x->{_a}; 4752 delete $x->{_p}; 4753 } 4754 # restore globals 4755 $$abr = $ab; 4756 $$pbr = $pb; 4757 $x; 4758} 4759 47601; 4761 4762__END__ 4763 4764=pod 4765 4766=head1 NAME 4767 4768Math::BigFloat - Arbitrary size floating point math package 4769 4770=head1 SYNOPSIS 4771 4772 use Math::BigFloat; 4773 4774 # Configuration methods (may be used as class methods and instance methods) 4775 4776 Math::BigFloat->accuracy(); # get class accuracy 4777 Math::BigFloat->accuracy($n); # set class accuracy 4778 Math::BigFloat->precision(); # get class precision 4779 Math::BigFloat->precision($n); # set class precision 4780 Math::BigFloat->round_mode(); # get class rounding mode 4781 Math::BigFloat->round_mode($m); # set global round mode, must be one of 4782 # 'even', 'odd', '+inf', '-inf', 'zero', 4783 # 'trunc', or 'common' 4784 Math::BigFloat->config("lib"); # name of backend math library 4785 4786 # Constructor methods (when the class methods below are used as instance 4787 # methods, the value is assigned the invocand) 4788 4789 $x = Math::BigFloat->new($str); # defaults to 0 4790 $x = Math::BigFloat->new('0x123'); # from hexadecimal 4791 $x = Math::BigFloat->new('0b101'); # from binary 4792 $x = Math::BigFloat->from_hex('0xc.afep+3'); # from hex 4793 $x = Math::BigFloat->from_hex('cafe'); # ditto 4794 $x = Math::BigFloat->from_oct('1.3267p-4'); # from octal 4795 $x = Math::BigFloat->from_oct('0377'); # ditto 4796 $x = Math::BigFloat->from_bin('0b1.1001p-4'); # from binary 4797 $x = Math::BigFloat->from_bin('0101'); # ditto 4798 $x = Math::BigFloat->bzero(); # create a +0 4799 $x = Math::BigFloat->bone(); # create a +1 4800 $x = Math::BigFloat->bone('-'); # create a -1 4801 $x = Math::BigFloat->binf(); # create a +inf 4802 $x = Math::BigFloat->binf('-'); # create a -inf 4803 $x = Math::BigFloat->bnan(); # create a Not-A-Number 4804 $x = Math::BigFloat->bpi(); # returns pi 4805 4806 $y = $x->copy(); # make a copy (unlike $y = $x) 4807 $y = $x->as_int(); # return as BigInt 4808 4809 # Boolean methods (these don't modify the invocand) 4810 4811 $x->is_zero(); # if $x is 0 4812 $x->is_one(); # if $x is +1 4813 $x->is_one("+"); # ditto 4814 $x->is_one("-"); # if $x is -1 4815 $x->is_inf(); # if $x is +inf or -inf 4816 $x->is_inf("+"); # if $x is +inf 4817 $x->is_inf("-"); # if $x is -inf 4818 $x->is_nan(); # if $x is NaN 4819 4820 $x->is_positive(); # if $x > 0 4821 $x->is_pos(); # ditto 4822 $x->is_negative(); # if $x < 0 4823 $x->is_neg(); # ditto 4824 4825 $x->is_odd(); # if $x is odd 4826 $x->is_even(); # if $x is even 4827 $x->is_int(); # if $x is an integer 4828 4829 # Comparison methods 4830 4831 $x->bcmp($y); # compare numbers (undef, < 0, == 0, > 0) 4832 $x->bacmp($y); # compare absolutely (undef, < 0, == 0, > 0) 4833 $x->beq($y); # true if and only if $x == $y 4834 $x->bne($y); # true if and only if $x != $y 4835 $x->blt($y); # true if and only if $x < $y 4836 $x->ble($y); # true if and only if $x <= $y 4837 $x->bgt($y); # true if and only if $x > $y 4838 $x->bge($y); # true if and only if $x >= $y 4839 4840 # Arithmetic methods 4841 4842 $x->bneg(); # negation 4843 $x->babs(); # absolute value 4844 $x->bsgn(); # sign function (-1, 0, 1, or NaN) 4845 $x->bnorm(); # normalize (no-op) 4846 $x->binc(); # increment $x by 1 4847 $x->bdec(); # decrement $x by 1 4848 $x->badd($y); # addition (add $y to $x) 4849 $x->bsub($y); # subtraction (subtract $y from $x) 4850 $x->bmul($y); # multiplication (multiply $x by $y) 4851 $x->bmuladd($y,$z); # $x = $x * $y + $z 4852 $x->bdiv($y); # division (floored), set $x to quotient 4853 # return (quo,rem) or quo if scalar 4854 $x->btdiv($y); # division (truncated), set $x to quotient 4855 # return (quo,rem) or quo if scalar 4856 $x->bmod($y); # modulus (x % y) 4857 $x->btmod($y); # modulus (truncated) 4858 $x->bmodinv($mod); # modular multiplicative inverse 4859 $x->bmodpow($y,$mod); # modular exponentiation (($x ** $y) % $mod) 4860 $x->bpow($y); # power of arguments (x ** y) 4861 $x->blog(); # logarithm of $x to base e (Euler's number) 4862 $x->blog($base); # logarithm of $x to base $base (e.g., base 2) 4863 $x->bexp(); # calculate e ** $x where e is Euler's number 4864 $x->bnok($y); # x over y (binomial coefficient n over k) 4865 $x->bsin(); # sine 4866 $x->bcos(); # cosine 4867 $x->batan(); # inverse tangent 4868 $x->batan2($y); # two-argument inverse tangent 4869 $x->bsqrt(); # calculate square root 4870 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root) 4871 $x->bfac(); # factorial of $x (1*2*3*4*..$x) 4872 4873 $x->blsft($n); # left shift $n places in base 2 4874 $x->blsft($n,$b); # left shift $n places in base $b 4875 # returns (quo,rem) or quo (scalar context) 4876 $x->brsft($n); # right shift $n places in base 2 4877 $x->brsft($n,$b); # right shift $n places in base $b 4878 # returns (quo,rem) or quo (scalar context) 4879 4880 # Bitwise methods 4881 4882 $x->band($y); # bitwise and 4883 $x->bior($y); # bitwise inclusive or 4884 $x->bxor($y); # bitwise exclusive or 4885 $x->bnot(); # bitwise not (two's complement) 4886 4887 # Rounding methods 4888 $x->round($A,$P,$mode); # round to accuracy or precision using 4889 # rounding mode $mode 4890 $x->bround($n); # accuracy: preserve $n digits 4891 $x->bfround($n); # $n > 0: round to $nth digit left of dec. point 4892 # $n < 0: round to $nth digit right of dec. point 4893 $x->bfloor(); # round towards minus infinity 4894 $x->bceil(); # round towards plus infinity 4895 $x->bint(); # round towards zero 4896 4897 # Other mathematical methods 4898 4899 $x->bgcd($y); # greatest common divisor 4900 $x->blcm($y); # least common multiple 4901 4902 # Object property methods (do not modify the invocand) 4903 4904 $x->sign(); # the sign, either +, - or NaN 4905 $x->digit($n); # the nth digit, counting from the right 4906 $x->digit(-$n); # the nth digit, counting from the left 4907 $x->length(); # return number of digits in number 4908 ($xl,$f) = $x->length(); # length of number and length of fraction 4909 # part, latter is always 0 digits long 4910 # for Math::BigInt objects 4911 $x->mantissa(); # return (signed) mantissa as BigInt 4912 $x->exponent(); # return exponent as BigInt 4913 $x->parts(); # return (mantissa,exponent) as BigInt 4914 $x->sparts(); # mantissa and exponent (as integers) 4915 $x->nparts(); # mantissa and exponent (normalised) 4916 $x->eparts(); # mantissa and exponent (engineering notation) 4917 $x->dparts(); # integer and fraction part 4918 4919 # Conversion methods (do not modify the invocand) 4920 4921 $x->bstr(); # decimal notation, possibly zero padded 4922 $x->bsstr(); # string in scientific notation with integers 4923 $x->bnstr(); # string in normalized notation 4924 $x->bestr(); # string in engineering notation 4925 $x->bdstr(); # string in decimal notation 4926 $x->as_hex(); # as signed hexadecimal string with prefixed 0x 4927 $x->as_bin(); # as signed binary string with prefixed 0b 4928 $x->as_oct(); # as signed octal string with prefixed 0 4929 4930 # Other conversion methods 4931 4932 $x->numify(); # return as scalar (might overflow or underflow) 4933 4934=head1 DESCRIPTION 4935 4936Math::BigFloat provides support for arbitrary precision floating point. 4937Overloading is also provided for Perl operators. 4938 4939All operators (including basic math operations) are overloaded if you 4940declare your big floating point numbers as 4941 4942 $x = Math::BigFloat -> new('12_3.456_789_123_456_789E-2'); 4943 4944Operations with overloaded operators preserve the arguments, which is 4945exactly what you expect. 4946 4947=head2 Input 4948 4949Input values to these routines may be any scalar number or string that looks 4950like a number and represents a floating point number. 4951 4952=over 4953 4954=item * 4955 4956Leading and trailing whitespace is ignored. 4957 4958=item * 4959 4960Leading and trailing zeros are ignored. 4961 4962=item * 4963 4964If the string has a "0x" prefix, it is interpreted as a hexadecimal number. 4965 4966=item * 4967 4968If the string has a "0b" prefix, it is interpreted as a binary number. 4969 4970=item * 4971 4972For hexadecimal and binary numbers, the exponent must be separated from the 4973significand (mantissa) by the letter "p" or "P", not "e" or "E" as with decimal 4974numbers. 4975 4976=item * 4977 4978One underline is allowed between any two digits, including hexadecimal and 4979binary digits. 4980 4981=item * 4982 4983If the string can not be interpreted, NaN is returned. 4984 4985=back 4986 4987Octal numbers are typically prefixed by "0", but since leading zeros are 4988stripped, these methods can not automatically recognize octal numbers, so use 4989the constructor from_oct() to interpret octal strings. 4990 4991Some examples of valid string input 4992 4993 Input string Resulting value 4994 123 123 4995 1.23e2 123 4996 12300e-2 123 4997 0xcafe 51966 4998 0b1101 13 4999 67_538_754 67538754 5000 -4_5_6.7_8_9e+0_1_0 -4567890000000 5001 0x1.921fb5p+1 3.14159262180328369140625e+0 5002 0b1.1001p-4 9.765625e-2 5003 5004=head2 Output 5005 5006Output values are usually Math::BigFloat objects. 5007 5008Boolean operators C<is_zero()>, C<is_one()>, C<is_inf()>, etc. return true or 5009false. 5010 5011Comparison operators C<bcmp()> and C<bacmp()>) return -1, 0, 1, or 5012undef. 5013 5014=head1 METHODS 5015 5016Math::BigFloat supports all methods that Math::BigInt supports, except it 5017calculates non-integer results when possible. Please see L<Math::BigInt> for a 5018full description of each method. Below are just the most important differences: 5019 5020=head2 Configuration methods 5021 5022=over 5023 5024=item accuracy() 5025 5026 $x->accuracy(5); # local for $x 5027 CLASS->accuracy(5); # global for all members of CLASS 5028 # Note: This also applies to new()! 5029 5030 $A = $x->accuracy(); # read out accuracy that affects $x 5031 $A = CLASS->accuracy(); # read out global accuracy 5032 5033Set or get the global or local accuracy, aka how many significant digits the 5034results have. If you set a global accuracy, then this also applies to new()! 5035 5036Warning! The accuracy I<sticks>, e.g. once you created a number under the 5037influence of C<< CLASS->accuracy($A) >>, all results from math operations with 5038that number will also be rounded. 5039 5040In most cases, you should probably round the results explicitly using one of 5041L<Math::BigInt/round()>, L<Math::BigInt/bround()> or L<Math::BigInt/bfround()> 5042or by passing the desired accuracy to the math operation as additional 5043parameter: 5044 5045 my $x = Math::BigInt->new(30000); 5046 my $y = Math::BigInt->new(7); 5047 print scalar $x->copy()->bdiv($y, 2); # print 4300 5048 print scalar $x->copy()->bdiv($y)->bround(2); # print 4300 5049 5050=item precision() 5051 5052 $x->precision(-2); # local for $x, round at the second 5053 # digit right of the dot 5054 $x->precision(2); # ditto, round at the second digit 5055 # left of the dot 5056 5057 CLASS->precision(5); # Global for all members of CLASS 5058 # This also applies to new()! 5059 CLASS->precision(-5); # ditto 5060 5061 $P = CLASS->precision(); # read out global precision 5062 $P = $x->precision(); # read out precision that affects $x 5063 5064Note: You probably want to use L</accuracy()> instead. With L</accuracy()> you 5065set the number of digits each result should have, with L</precision()> you 5066set the place where to round! 5067 5068=back 5069 5070=head2 Constructor methods 5071 5072=over 5073 5074=item from_hex() 5075 5076 $x -> from_hex("0x1.921fb54442d18p+1"); 5077 $x = Math::BigFloat -> from_hex("0x1.921fb54442d18p+1"); 5078 5079Interpret input as a hexadecimal string.A prefix ("0x", "x", ignoring case) is 5080optional. A single underscore character ("_") may be placed between any two 5081digits. If the input is invalid, a NaN is returned. The exponent is in base 2 5082using decimal digits. 5083 5084If called as an instance method, the value is assigned to the invocand. 5085 5086=item from_oct() 5087 5088 $x -> from_oct("1.3267p-4"); 5089 $x = Math::BigFloat -> from_oct("1.3267p-4"); 5090 5091Interpret input as an octal string. A single underscore character ("_") may be 5092placed between any two digits. If the input is invalid, a NaN is returned. The 5093exponent is in base 2 using decimal digits. 5094 5095If called as an instance method, the value is assigned to the invocand. 5096 5097=item from_bin() 5098 5099 $x -> from_bin("0b1.1001p-4"); 5100 $x = Math::BigFloat -> from_bin("0b1.1001p-4"); 5101 5102Interpret input as a hexadecimal string. A prefix ("0b" or "b", ignoring case) 5103is optional. A single underscore character ("_") may be placed between any two 5104digits. If the input is invalid, a NaN is returned. The exponent is in base 2 5105using decimal digits. 5106 5107If called as an instance method, the value is assigned to the invocand. 5108 5109=item bpi() 5110 5111 print Math::BigFloat->bpi(100), "\n"; 5112 5113Calculate PI to N digits (including the 3 before the dot). The result is 5114rounded according to the current rounding mode, which defaults to "even". 5115 5116This method was added in v1.87 of Math::BigInt (June 2007). 5117 5118=back 5119 5120=head2 Arithmetic methods 5121 5122=over 5123 5124=item bmuladd() 5125 5126 $x->bmuladd($y,$z); 5127 5128Multiply $x by $y, and then add $z to the result. 5129 5130This method was added in v1.87 of Math::BigInt (June 2007). 5131 5132=item bdiv() 5133 5134 $q = $x->bdiv($y); 5135 ($q, $r) = $x->bdiv($y); 5136 5137In scalar context, divides $x by $y and returns the result to the given or 5138default accuracy/precision. In list context, does floored division 5139(F-division), returning an integer $q and a remainder $r so that $x = $q * $y + 5140$r. The remainer (modulo) is equal to what is returned by C<< $x->bmod($y) >>. 5141 5142=item bmod() 5143 5144 $x->bmod($y); 5145 5146Returns $x modulo $y. When $x is finite, and $y is finite and non-zero, the 5147result is identical to the remainder after floored division (F-division). If, 5148in addition, both $x and $y are integers, the result is identical to the result 5149from Perl's % operator. 5150 5151=item bexp() 5152 5153 $x->bexp($accuracy); # calculate e ** X 5154 5155Calculates the expression C<e ** $x> where C<e> is Euler's number. 5156 5157This method was added in v1.82 of Math::BigInt (April 2007). 5158 5159=item bnok() 5160 5161 $x->bnok($y); # x over y (binomial coefficient n over k) 5162 5163Calculates the binomial coefficient n over k, also called the "choose" 5164function. The result is equivalent to: 5165 5166 ( n ) n! 5167 | - | = ------- 5168 ( k ) k!(n-k)! 5169 5170This method was added in v1.84 of Math::BigInt (April 2007). 5171 5172=item bsin() 5173 5174 my $x = Math::BigFloat->new(1); 5175 print $x->bsin(100), "\n"; 5176 5177Calculate the sinus of $x, modifying $x in place. 5178 5179This method was added in v1.87 of Math::BigInt (June 2007). 5180 5181=item bcos() 5182 5183 my $x = Math::BigFloat->new(1); 5184 print $x->bcos(100), "\n"; 5185 5186Calculate the cosinus of $x, modifying $x in place. 5187 5188This method was added in v1.87 of Math::BigInt (June 2007). 5189 5190=item batan() 5191 5192 my $x = Math::BigFloat->new(1); 5193 print $x->batan(100), "\n"; 5194 5195Calculate the arcus tanges of $x, modifying $x in place. See also L</batan2()>. 5196 5197This method was added in v1.87 of Math::BigInt (June 2007). 5198 5199=item batan2() 5200 5201 my $y = Math::BigFloat->new(2); 5202 my $x = Math::BigFloat->new(3); 5203 print $y->batan2($x), "\n"; 5204 5205Calculate the arcus tanges of C<$y> divided by C<$x>, modifying $y in place. 5206See also L</batan()>. 5207 5208This method was added in v1.87 of Math::BigInt (June 2007). 5209 5210=item as_float() 5211 5212This method is called when Math::BigFloat encounters an object it doesn't know 5213how to handle. For instance, assume $x is a Math::BigFloat, or subclass 5214thereof, and $y is defined, but not a Math::BigFloat, or subclass thereof. If 5215you do 5216 5217 $x -> badd($y); 5218 5219$y needs to be converted into an object that $x can deal with. This is done by 5220first checking if $y is something that $x might be upgraded to. If that is the 5221case, no further attempts are made. The next is to see if $y supports the 5222method C<as_float()>. The method C<as_float()> is expected to return either an 5223object that has the same class as $x, a subclass thereof, or a string that 5224C<ref($x)-E<gt>new()> can parse to create an object. 5225 5226In Math::BigFloat, C<as_float()> has the same effect as C<copy()>. 5227 5228=back 5229 5230=head2 ACCURACY AND PRECISION 5231 5232See also: L<Rounding|/Rounding>. 5233 5234Math::BigFloat supports both precision (rounding to a certain place before or 5235after the dot) and accuracy (rounding to a certain number of digits). For a 5236full documentation, examples and tips on these topics please see the large 5237section about rounding in L<Math::BigInt>. 5238 5239Since things like C<sqrt(2)> or C<1 / 3> must presented with a limited 5240accuracy lest a operation consumes all resources, each operation produces 5241no more than the requested number of digits. 5242 5243If there is no global precision or accuracy set, B<and> the operation in 5244question was not called with a requested precision or accuracy, B<and> the 5245input $x has no accuracy or precision set, then a fallback parameter will 5246be used. For historical reasons, it is called C<div_scale> and can be accessed 5247via: 5248 5249 $d = Math::BigFloat->div_scale(); # query 5250 Math::BigFloat->div_scale($n); # set to $n digits 5251 5252The default value for C<div_scale> is 40. 5253 5254In case the result of one operation has more digits than specified, 5255it is rounded. The rounding mode taken is either the default mode, or the one 5256supplied to the operation after the I<scale>: 5257 5258 $x = Math::BigFloat->new(2); 5259 Math::BigFloat->accuracy(5); # 5 digits max 5260 $y = $x->copy()->bdiv(3); # gives 0.66667 5261 $y = $x->copy()->bdiv(3,6); # gives 0.666667 5262 $y = $x->copy()->bdiv(3,6,undef,'odd'); # gives 0.666667 5263 Math::BigFloat->round_mode('zero'); 5264 $y = $x->copy()->bdiv(3,6); # will also give 0.666667 5265 5266Note that C<< Math::BigFloat->accuracy() >> and C<< Math::BigFloat->precision() >> 5267set the global variables, and thus B<any> newly created number will be subject 5268to the global rounding B<immediately>. This means that in the examples above, the 5269C<3> as argument to C<bdiv()> will also get an accuracy of B<5>. 5270 5271It is less confusing to either calculate the result fully, and afterwards 5272round it explicitly, or use the additional parameters to the math 5273functions like so: 5274 5275 use Math::BigFloat; 5276 $x = Math::BigFloat->new(2); 5277 $y = $x->copy()->bdiv(3); 5278 print $y->bround(5),"\n"; # gives 0.66667 5279 5280 or 5281 5282 use Math::BigFloat; 5283 $x = Math::BigFloat->new(2); 5284 $y = $x->copy()->bdiv(3,5); # gives 0.66667 5285 print "$y\n"; 5286 5287=head2 Rounding 5288 5289=over 5290 5291=item bfround ( +$scale ) 5292 5293Rounds to the $scale'th place left from the '.', counting from the dot. 5294The first digit is numbered 1. 5295 5296=item bfround ( -$scale ) 5297 5298Rounds to the $scale'th place right from the '.', counting from the dot. 5299 5300=item bfround ( 0 ) 5301 5302Rounds to an integer. 5303 5304=item bround ( +$scale ) 5305 5306Preserves accuracy to $scale digits from the left (aka significant digits) and 5307pads the rest with zeros. If the number is between 1 and -1, the significant 5308digits count from the first non-zero after the '.' 5309 5310=item bround ( -$scale ) and bround ( 0 ) 5311 5312These are effectively no-ops. 5313 5314=back 5315 5316All rounding functions take as a second parameter a rounding mode from one of 5317the following: 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'. 5318 5319The default rounding mode is 'even'. By using 5320C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default 5321mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is 5322no longer supported. 5323The second parameter to the round functions then overrides the default 5324temporarily. 5325 5326The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses 5327'trunc' as rounding mode to make it equivalent to: 5328 5329 $x = 2.5; 5330 $y = int($x) + 2; 5331 5332You can override this by passing the desired rounding mode as parameter to 5333C<as_number()>: 5334 5335 $x = Math::BigFloat->new(2.5); 5336 $y = $x->as_number('odd'); # $y = 3 5337 5338=head1 Autocreating constants 5339 5340After C<use Math::BigFloat ':constant'> all the floating point constants 5341in the given scope are converted to C<Math::BigFloat>. This conversion 5342happens at compile time. 5343 5344In particular 5345 5346 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"' 5347 5348prints the value of C<2E-100>. Note that without conversion of 5349constants the expression 2E-100 will be calculated as normal floating point 5350number. 5351 5352Please note that ':constant' does not affect integer constants, nor binary 5353nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to 5354work. 5355 5356=head2 Math library 5357 5358Math with the numbers is done (by default) by a module called 5359Math::BigInt::Calc. This is equivalent to saying: 5360 5361 use Math::BigFloat lib => 'Calc'; 5362 5363You can change this by using: 5364 5365 use Math::BigFloat lib => 'GMP'; 5366 5367B<Note>: General purpose packages should not be explicit about the library 5368to use; let the script author decide which is best. 5369 5370Note: The keyword 'lib' will warn when the requested library could not be 5371loaded. To suppress the warning use 'try' instead: 5372 5373 use Math::BigFloat try => 'GMP'; 5374 5375If your script works with huge numbers and Calc is too slow for them, 5376you can also for the loading of one of these libraries and if none 5377of them can be used, the code will die: 5378 5379 use Math::BigFloat only => 'GMP,Pari'; 5380 5381The following would first try to find Math::BigInt::Foo, then 5382Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc: 5383 5384 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar'; 5385 5386See the respective low-level library documentation for further details. 5387 5388Please note that Math::BigFloat does B<not> use the denoted library itself, 5389but it merely passes the lib argument to Math::BigInt. So, instead of the need 5390to do: 5391 5392 use Math::BigInt lib => 'GMP'; 5393 use Math::BigFloat; 5394 5395you can roll it all into one line: 5396 5397 use Math::BigFloat lib => 'GMP'; 5398 5399It is also possible to just require Math::BigFloat: 5400 5401 require Math::BigFloat; 5402 5403This will load the necessary things (like BigInt) when they are needed, and 5404automatically. 5405 5406See L<Math::BigInt> for more details than you ever wanted to know about using 5407a different low-level library. 5408 5409=head2 Using Math::BigInt::Lite 5410 5411For backwards compatibility reasons it is still possible to 5412request a different storage class for use with Math::BigFloat: 5413 5414 use Math::BigFloat with => 'Math::BigInt::Lite'; 5415 5416However, this request is ignored, as the current code now uses the low-level 5417math library for directly storing the number parts. 5418 5419=head1 EXPORTS 5420 5421C<Math::BigFloat> exports nothing by default, but can export the C<bpi()> method: 5422 5423 use Math::BigFloat qw/bpi/; 5424 5425 print bpi(10), "\n"; 5426 5427=head1 CAVEATS 5428 5429Do not try to be clever to insert some operations in between switching 5430libraries: 5431 5432 require Math::BigFloat; 5433 my $matter = Math::BigFloat->bone() + 4; # load BigInt and Calc 5434 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too 5435 my $anti_matter = Math::BigFloat->bone()+4; # now use Pari 5436 5437This will create objects with numbers stored in two different backend libraries, 5438and B<VERY BAD THINGS> will happen when you use these together: 5439 5440 my $flash_and_bang = $matter + $anti_matter; # Don't do this! 5441 5442=over 5443 5444=item stringify, bstr() 5445 5446Both stringify and bstr() now drop the leading '+'. The old code would return 5447'+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for 5448reasoning and details. 5449 5450=item brsft() 5451 5452The following will probably not print what you expect: 5453 5454 my $c = Math::BigFloat->new('3.14159'); 5455 print $c->brsft(3,10),"\n"; # prints 0.00314153.1415 5456 5457It prints both quotient and remainder, since print calls C<brsft()> in list 5458context. Also, C<< $c->brsft() >> will modify $c, so be careful. 5459You probably want to use 5460 5461 print scalar $c->copy()->brsft(3,10),"\n"; 5462 # or if you really want to modify $c 5463 print scalar $c->brsft(3,10),"\n"; 5464 5465instead. 5466 5467=item Modifying and = 5468 5469Beware of: 5470 5471 $x = Math::BigFloat->new(5); 5472 $y = $x; 5473 5474It will not do what you think, e.g. making a copy of $x. Instead it just makes 5475a second reference to the B<same> object and stores it in $y. Thus anything 5476that modifies $x will modify $y (except overloaded math operators), and vice 5477versa. See L<Math::BigInt> for details and how to avoid that. 5478 5479=item precision() vs. accuracy() 5480 5481A common pitfall is to use L</precision()> when you want to round a result to 5482a certain number of digits: 5483 5484 use Math::BigFloat; 5485 5486 Math::BigFloat->precision(4); # does not do what you 5487 # think it does 5488 my $x = Math::BigFloat->new(12345); # rounds $x to "12000"! 5489 print "$x\n"; # print "12000" 5490 my $y = Math::BigFloat->new(3); # rounds $y to "0"! 5491 print "$y\n"; # print "0" 5492 $z = $x / $y; # 12000 / 0 => NaN! 5493 print "$z\n"; 5494 print $z->precision(),"\n"; # 4 5495 5496Replacing L</precision()> with L</accuracy()> is probably not what you want, either: 5497 5498 use Math::BigFloat; 5499 5500 Math::BigFloat->accuracy(4); # enables global rounding: 5501 my $x = Math::BigFloat->new(123456); # rounded immediately 5502 # to "12350" 5503 print "$x\n"; # print "123500" 5504 my $y = Math::BigFloat->new(3); # rounded to "3 5505 print "$y\n"; # print "3" 5506 print $z = $x->copy()->bdiv($y),"\n"; # 41170 5507 print $z->accuracy(),"\n"; # 4 5508 5509What you want to use instead is: 5510 5511 use Math::BigFloat; 5512 5513 my $x = Math::BigFloat->new(123456); # no rounding 5514 print "$x\n"; # print "123456" 5515 my $y = Math::BigFloat->new(3); # no rounding 5516 print "$y\n"; # print "3" 5517 print $z = $x->copy()->bdiv($y,4),"\n"; # 41150 5518 print $z->accuracy(),"\n"; # undef 5519 5520In addition to computing what you expected, the last example also does B<not> 5521"taint" the result with an accuracy or precision setting, which would 5522influence any further operation. 5523 5524=back 5525 5526=head1 BUGS 5527 5528Please report any bugs or feature requests to 5529C<bug-math-bigint at rt.cpan.org>, or through the web interface at 5530L<https://rt.cpan.org/Ticket/Create.html?Queue=Math-BigInt> 5531(requires login). 5532We will be notified, and then you'll automatically be notified of progress on 5533your bug as I make changes. 5534 5535=head1 SUPPORT 5536 5537You can find documentation for this module with the perldoc command. 5538 5539 perldoc Math::BigFloat 5540 5541You can also look for information at: 5542 5543=over 4 5544 5545=item * RT: CPAN's request tracker 5546 5547L<https://rt.cpan.org/Public/Dist/Display.html?Name=Math-BigInt> 5548 5549=item * AnnoCPAN: Annotated CPAN documentation 5550 5551L<http://annocpan.org/dist/Math-BigInt> 5552 5553=item * CPAN Ratings 5554 5555L<http://cpanratings.perl.org/dist/Math-BigInt> 5556 5557=item * Search CPAN 5558 5559L<http://search.cpan.org/dist/Math-BigInt/> 5560 5561=item * CPAN Testers Matrix 5562 5563L<http://matrix.cpantesters.org/?dist=Math-BigInt> 5564 5565=item * The Bignum mailing list 5566 5567=over 4 5568 5569=item * Post to mailing list 5570 5571C<bignum at lists.scsys.co.uk> 5572 5573=item * View mailing list 5574 5575L<http://lists.scsys.co.uk/pipermail/bignum/> 5576 5577=item * Subscribe/Unsubscribe 5578 5579L<http://lists.scsys.co.uk/cgi-bin/mailman/listinfo/bignum> 5580 5581=back 5582 5583=back 5584 5585=head1 LICENSE 5586 5587This program is free software; you may redistribute it and/or modify it under 5588the same terms as Perl itself. 5589 5590=head1 SEE ALSO 5591 5592L<Math::BigFloat> and L<Math::BigInt> as well as the backends 5593L<Math::BigInt::FastCalc>, L<Math::BigInt::GMP>, and L<Math::BigInt::Pari>. 5594 5595The pragmas L<bignum>, L<bigint> and L<bigrat> also might be of interest 5596because they solve the autoupgrading/downgrading issue, at least partly. 5597 5598=head1 AUTHORS 5599 5600=over 4 5601 5602=item * 5603 5604Mark Biggar, overloaded interface by Ilya Zakharevich, 1996-2001. 5605 5606=item * 5607 5608Completely rewritten by Tels L<http://bloodgate.com> in 2001-2008. 5609 5610=item * 5611 5612Florian Ragwitz E<lt>flora@cpan.orgE<gt>, 2010. 5613 5614=item * 5615 5616Peter John Acklam E<lt>pjacklam@online.noE<gt>, 2011-. 5617 5618=back 5619 5620=cut 5621