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# 9# sign : "+", "-", "+inf", "-inf", or "NaN" 10# _m : absolute value of mantissa ($LIB thingy) 11# _es : sign of exponent ("+" or "-") 12# _e : absolute value of exponent ($LIB thingy) 13# accuracy : accuracy (scalar) 14# precision : precision (scalar) 15 16use 5.006001; 17use strict; 18use warnings; 19 20use Carp qw< carp croak >; 21use Scalar::Util qw< blessed >; 22use Math::BigInt qw< >; 23 24our $VERSION = '2.003002'; 25$VERSION =~ tr/_//d; 26 27require Exporter; 28our @ISA = qw/Math::BigInt/; 29our @EXPORT_OK = qw/bpi/; 30 31# $_trap_inf/$_trap_nan are internal and should never be accessed from outside 32our ($AUTOLOAD, $accuracy, $precision, $div_scale, $round_mode, $rnd_mode, 33 $upgrade, $downgrade, $_trap_nan, $_trap_inf); 34 35use overload 36 37 # overload key: with_assign 38 39 '+' => sub { $_[0] -> copy() -> badd($_[1]); }, 40 41 '-' => sub { my $c = $_[0] -> copy(); 42 $_[2] ? $c -> bneg() -> badd($_[1]) 43 : $c -> bsub($_[1]); }, 44 45 '*' => sub { $_[0] -> copy() -> bmul($_[1]); }, 46 47 '/' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bdiv($_[0]) 48 : $_[0] -> copy() -> bdiv($_[1]); }, 49 50 '%' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bmod($_[0]) 51 : $_[0] -> copy() -> bmod($_[1]); }, 52 53 '**' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bpow($_[0]) 54 : $_[0] -> copy() -> bpow($_[1]); }, 55 56 '<<' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bblsft($_[0]) 57 : $_[0] -> copy() -> bblsft($_[1]); }, 58 59 '>>' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bbrsft($_[0]) 60 : $_[0] -> copy() -> bbrsft($_[1]); }, 61 62 # overload key: assign 63 64 '+=' => sub { $_[0] -> badd($_[1]); }, 65 66 '-=' => sub { $_[0] -> bsub($_[1]); }, 67 68 '*=' => sub { $_[0] -> bmul($_[1]); }, 69 70 '/=' => sub { scalar $_[0] -> bdiv($_[1]); }, 71 72 '%=' => sub { $_[0] -> bmod($_[1]); }, 73 74 '**=' => sub { $_[0] -> bpow($_[1]); }, 75 76 '<<=' => sub { $_[0] -> bblsft($_[1]); }, 77 78 '>>=' => sub { $_[0] -> bbrsft($_[1]); }, 79 80# 'x=' => sub { }, 81 82# '.=' => sub { }, 83 84 # overload key: num_comparison 85 86 '<' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> blt($_[0]) 87 : $_[0] -> blt($_[1]); }, 88 89 '<=' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> ble($_[0]) 90 : $_[0] -> ble($_[1]); }, 91 92 '>' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bgt($_[0]) 93 : $_[0] -> bgt($_[1]); }, 94 95 '>=' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bge($_[0]) 96 : $_[0] -> bge($_[1]); }, 97 98 '==' => sub { $_[0] -> beq($_[1]); }, 99 100 '!=' => sub { $_[0] -> bne($_[1]); }, 101 102 # overload key: 3way_comparison 103 104 '<=>' => sub { my $cmp = $_[0] -> bcmp($_[1]); 105 defined($cmp) && $_[2] ? -$cmp : $cmp; }, 106 107 'cmp' => sub { $_[2] ? "$_[1]" cmp $_[0] -> bstr() 108 : $_[0] -> bstr() cmp "$_[1]"; }, 109 110 # overload key: str_comparison 111 112# 'lt' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrlt($_[0]) 113# : $_[0] -> bstrlt($_[1]); }, 114# 115# 'le' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrle($_[0]) 116# : $_[0] -> bstrle($_[1]); }, 117# 118# 'gt' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrgt($_[0]) 119# : $_[0] -> bstrgt($_[1]); }, 120# 121# 'ge' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrge($_[0]) 122# : $_[0] -> bstrge($_[1]); }, 123# 124# 'eq' => sub { $_[0] -> bstreq($_[1]); }, 125# 126# 'ne' => sub { $_[0] -> bstrne($_[1]); }, 127 128 # overload key: binary 129 130 '&' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> band($_[0]) 131 : $_[0] -> copy() -> band($_[1]); }, 132 133 '&=' => sub { $_[0] -> band($_[1]); }, 134 135 '|' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bior($_[0]) 136 : $_[0] -> copy() -> bior($_[1]); }, 137 138 '|=' => sub { $_[0] -> bior($_[1]); }, 139 140 '^' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bxor($_[0]) 141 : $_[0] -> copy() -> bxor($_[1]); }, 142 143 '^=' => sub { $_[0] -> bxor($_[1]); }, 144 145# '&.' => sub { }, 146 147# '&.=' => sub { }, 148 149# '|.' => sub { }, 150 151# '|.=' => sub { }, 152 153# '^.' => sub { }, 154 155# '^.=' => sub { }, 156 157 # overload key: unary 158 159 'neg' => sub { $_[0] -> copy() -> bneg(); }, 160 161# '!' => sub { }, 162 163 '~' => sub { $_[0] -> copy() -> bnot(); }, 164 165# '~.' => sub { }, 166 167 # overload key: mutators 168 169 '++' => sub { $_[0] -> binc() }, 170 171 '--' => sub { $_[0] -> bdec() }, 172 173 # overload key: func 174 175 'atan2' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> batan2($_[0]) 176 : $_[0] -> copy() -> batan2($_[1]); }, 177 178 'cos' => sub { $_[0] -> copy() -> bcos(); }, 179 180 'sin' => sub { $_[0] -> copy() -> bsin(); }, 181 182 'exp' => sub { $_[0] -> copy() -> bexp($_[1]); }, 183 184 'abs' => sub { $_[0] -> copy() -> babs(); }, 185 186 'log' => sub { $_[0] -> copy() -> blog(); }, 187 188 'sqrt' => sub { $_[0] -> copy() -> bsqrt(); }, 189 190 'int' => sub { $_[0] -> copy() -> bint(); }, 191 192 # overload key: conversion 193 194 'bool' => sub { $_[0] -> is_zero() ? '' : 1; }, 195 196 '""' => sub { $_[0] -> bstr(); }, 197 198 '0+' => sub { $_[0] -> numify(); }, 199 200 '=' => sub { $_[0] -> copy(); }, 201 202 ; 203 204############################################################################## 205# global constants, flags and assorted stuff 206 207# the following are public, but their usage is not recommended. Use the 208# accessor methods instead. 209 210# class constants, use Class->constant_name() to access 211# one of 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common' 212$round_mode = 'even'; 213$accuracy = undef; 214$precision = undef; 215$div_scale = 40; 216 217$upgrade = undef; 218$downgrade = undef; 219# the package we are using for our private parts, defaults to: 220# Math::BigInt->config('lib') 221my $LIB = 'Math::BigInt::Calc'; 222 223# are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config() 224$_trap_nan = 0; 225# the same for infinity 226$_trap_inf = 0; 227 228# constant for easier life 229my $nan = 'NaN'; 230 231# Has import() been called yet? This variable is needed to make "require" work. 232 233my $IMPORT = 0; 234 235# some digits of accuracy for blog(undef, 10); which we use in blog() for speed 236my $LOG_10 = 237 '2.3025850929940456840179914546843642076011014886287729760333279009675726097'; 238my $LOG_10_A = length($LOG_10)-1; 239# ditto for log(2) 240my $LOG_2 = 241 '0.6931471805599453094172321214581765680755001343602552541206800094933936220'; 242my $LOG_2_A = length($LOG_2)-1; 243my $HALF = '0.5'; # made into an object if nec. 244 245############################################################################## 246# the old code had $rnd_mode, so we need to support it, too 247 248sub TIESCALAR { 249 my ($class) = @_; 250 bless \$round_mode, $class; 251} 252 253sub FETCH { 254 return $round_mode; 255} 256 257sub STORE { 258 $rnd_mode = (ref $_[0]) -> round_mode($_[1]); 259} 260 261BEGIN { 262 *objectify = \&Math::BigInt::objectify; 263 264 # when someone sets $rnd_mode, we catch this and check the value to see 265 # whether it is valid or not. 266 $rnd_mode = 'even'; 267 tie $rnd_mode, 'Math::BigFloat'; 268 269 *as_number = \&as_int; 270} 271 272sub DESTROY { 273 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub 274} 275 276sub AUTOLOAD { 277 278 # Make fxxx() work by mapping fxxx() to Math::BigFloat::bxxx(). 279 280 my $name = $AUTOLOAD; 281 $name =~ s/^(.*):://; # strip package name 282 my $class = $1 || __PACKAGE__; 283 284 $class -> import() if $IMPORT == 0; 285 286 # E.g., "fabs" -> "babs", but "is_neg" -> "is_neg" 287 288 my $bname = $name; 289 $bname =~ s/^f/b/; 290 291 # Map, e.g., Math::BigFloat::fabs() to Math::BigFloat::babs() 292 293 if ($bname ne $name && Math::BigFloat -> can($bname)) { 294 no strict 'refs'; 295 return &{"Math::BigFloat::$bname"}(@_); 296 } 297 298 # Map, e.g., Math::BigFloat::babs() to Math::BigInt::babs() 299 300 elsif (Math::BigInt -> can($bname)) { 301 no strict 'refs'; 302 return &{"Math::BigInt::$bname"}(@_); 303 } 304 305 else { 306 croak("Can't call $class->$name(), not a valid method"); 307 } 308} 309 310############################################################################## 311 312# Compare the following function with @ISA above. This inheritance mess needs a 313# clean up. When doing so, also consider the BEGIN block and the AUTOLOAD code. 314# Fixme! 315 316sub isa { 317 my ($self, $class) = @_; 318 return if $class =~ /^Math::BigInt/; # we aren't one of these 319 UNIVERSAL::isa($self, $class); 320} 321 322sub config { 323 # return (later set?) configuration data as hash ref 324 my $class = shift || 'Math::BigFloat'; 325 326 # Getter/accessor. 327 328 if (@_ == 1 && ref($_[0]) ne 'HASH') { 329 my $param = shift; 330 return $class if $param eq 'class'; 331 return $LIB if $param eq 'with'; 332 return $class->SUPER::config($param); 333 } 334 335 # Setter. 336 337 my $cfg = $class->SUPER::config(@_); 338 339 # now we need only to override the ones that are different from our parent 340 $cfg->{class} = $class; 341 $cfg->{with} = $LIB; 342 $cfg; 343} 344 345############################################################################### 346# Constructor methods 347############################################################################### 348 349sub new { 350 # Create a new Math::BigFloat object from a string or another bigfloat 351 # object. 352 # _e: exponent 353 # _m: mantissa 354 # sign => ("+", "-", "+inf", "-inf", or "NaN") 355 356 my $self = shift; 357 my $selfref = ref $self; 358 my $class = $selfref || $self; 359 360 # Make "require" work. 361 362 $class -> import() if $IMPORT == 0; 363 364 # Although this use has been discouraged for more than 10 years, people 365 # apparently still use it, so we still support it. 366 367 return $class -> bzero() unless @_; 368 369 my ($wanted, @r) = @_; 370 371 if (!defined($wanted)) { 372 #if (warnings::enabled("uninitialized")) { 373 # warnings::warn("uninitialized", 374 # "Use of uninitialized value in new()"); 375 #} 376 return $class -> bzero(@r); 377 } 378 379 if (!ref($wanted) && $wanted eq "") { 380 #if (warnings::enabled("numeric")) { 381 # warnings::warn("numeric", 382 # q|Argument "" isn't numeric in new()|); 383 #} 384 #return $class -> bzero(@r); 385 return $class -> bnan(@r); 386 } 387 388 # Initialize a new object. 389 390 $self = bless {}, $class unless $selfref; 391 392 # Math::BigFloat or subclass 393 394 if (defined(blessed($wanted)) && $wanted -> isa(__PACKAGE__)) { 395 396 # Don't copy the accuracy and precision, because a new object should get 397 # them from the global configuration. 398 399 $self -> {sign} = $wanted -> {sign}; 400 $self -> {_m} = $LIB -> _copy($wanted -> {_m}); 401 $self -> {_es} = $wanted -> {_es}; 402 $self -> {_e} = $LIB -> _copy($wanted -> {_e}); 403 $self = $self->round(@r) 404 unless @r >= 2 && !defined($r[0]) && !defined($r[1]); 405 return $self; 406 } 407 408 # Shortcut for Math::BigInt and its subclasses. This should be improved. 409 410 if (defined(blessed($wanted))) { 411 if ($wanted -> isa('Math::BigInt')) { 412 $self->{sign} = $wanted -> {sign}; 413 $self->{_m} = $LIB -> _copy($wanted -> {value}); 414 $self->{_es} = '+'; 415 $self->{_e} = $LIB -> _zero(); 416 return $self -> bnorm(); 417 } 418 419 if ($wanted -> can("as_number")) { 420 $self->{sign} = $wanted -> sign(); 421 $self->{_m} = $wanted -> as_number() -> {value}; 422 $self->{_es} = '+'; 423 $self->{_e} = $LIB -> _zero(); 424 return $self -> bnorm(); 425 } 426 } 427 428 # Shortcut for simple forms like '123' that have no trailing zeros. Trailing 429 # zeros would require a non-zero exponent. 430 431 if ($wanted =~ 432 / ^ 433 \s* # optional leading whitespace 434 ( [+-]? ) # optional sign 435 0* # optional leading zeros 436 ( [1-9] (?: [0-9]* [1-9] )? ) # significand 437 \s* # optional trailing whitespace 438 $ 439 /x) 440 { 441 return $downgrade -> new($1 . $2) if defined $downgrade; 442 $self->{sign} = $1 || '+'; 443 $self->{_m} = $LIB -> _new($2); 444 $self->{_es} = '+'; 445 $self->{_e} = $LIB -> _zero(); 446 $self = $self->round(@r) 447 unless @r >= 2 && !defined $r[0] && !defined $r[1]; 448 return $self; 449 } 450 451 # Handle Infs. 452 453 if ($wanted =~ / ^ 454 \s* 455 ( [+-]? ) 456 inf (?: inity )? 457 \s* 458 \z 459 /ix) 460 { 461 my $sgn = $1 || '+'; 462 return $class -> binf($sgn, @r); 463 } 464 465 # Handle explicit NaNs (not the ones returned due to invalid input). 466 467 if ($wanted =~ / ^ 468 \s* 469 ( [+-]? ) 470 nan 471 \s* 472 \z 473 /ix) 474 { 475 return $class -> bnan(@r); 476 } 477 478 my @parts; 479 480 if ( 481 # Handle hexadecimal numbers. We auto-detect hexadecimal numbers if they 482 # have a "0x", "0X", "x", or "X" prefix, cf. CORE::oct(). 483 484 $wanted =~ /^\s*[+-]?0?[Xx]/ and 485 @parts = $class -> _hex_str_to_flt_lib_parts($wanted) 486 487 or 488 489 # Handle octal numbers. We auto-detect octal numbers if they have a 490 # "0o", "0O", "o", "O" prefix, cf. CORE::oct(). 491 492 $wanted =~ /^\s*[+-]?0?[Oo]/ and 493 @parts = $class -> _oct_str_to_flt_lib_parts($wanted) 494 495 or 496 497 # Handle binary numbers. We auto-detect binary numbers if they have a 498 # "0b", "0B", "b", or "B" prefix, cf. CORE::oct(). 499 500 $wanted =~ /^\s*[+-]?0?[Bb]/ and 501 @parts = $class -> _bin_str_to_flt_lib_parts($wanted) 502 503 or 504 505 # At this point, what is left are decimal numbers that aren't handled 506 # above and octal floating point numbers that don't have any of the 507 # "0o", "0O", "o", or "O" prefixes. First see if it is a decimal number. 508 509 @parts = $class -> _dec_str_to_flt_lib_parts($wanted) 510 or 511 512 # See if it is an octal floating point number. The extra check is 513 # included because _oct_str_to_flt_lib_parts() accepts octal numbers 514 # that don't have a prefix (this is needed to make it work with, e.g., 515 # from_oct() that don't require a prefix). However, Perl requires a 516 # prefix for octal floating point literals. For example, "1p+0" is not 517 # valid, but "01p+0" and "0__1p+0" are. 518 519 $wanted =~ /^\s*[+-]?0_*\d/ and 520 @parts = $class -> _oct_str_to_flt_lib_parts($wanted)) 521 { 522 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts; 523 524 $self = $self->round(@r) 525 unless @r >= 2 && !defined($r[0]) && !defined($r[1]); 526 527 return $downgrade -> new($self -> bdstr(), @r) 528 if defined($downgrade) && $self -> is_int(); 529 return $self; 530 } 531 532 # If we get here, the value is neither a valid decimal, binary, octal, or 533 # hexadecimal number. It is not an explicit Inf or a NaN either. 534 535 return $class -> bnan(@r); 536} 537 538sub from_dec { 539 my $self = shift; 540 my $selfref = ref $self; 541 my $class = $selfref || $self; 542 543 # Make "require" work. 544 545 $class -> import() if $IMPORT == 0; 546 547 # Don't modify constant (read-only) objects. 548 549 return $self if $selfref && $self->modify('from_dec'); 550 551 my $str = shift; 552 my @r = @_; 553 554 # If called as a class method, initialize a new object. 555 556 $self = bless {}, $class unless $selfref; 557 558 if (my @parts = $class -> _dec_str_to_flt_lib_parts($str)) { 559 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts; 560 561 $self = $self->round(@r) 562 unless @r >= 2 && !defined($r[0]) && !defined($r[1]); 563 564 return $downgrade -> new($self -> bdstr(), @r) 565 if defined($downgrade) && $self -> is_int(); 566 return $self; 567 } 568 569 return $self -> bnan(@r); 570} 571 572sub from_hex { 573 my $self = shift; 574 my $selfref = ref $self; 575 my $class = $selfref || $self; 576 577 # Make "require" work. 578 579 $class -> import() if $IMPORT == 0; 580 581 # Don't modify constant (read-only) objects. 582 583 return $self if $selfref && $self->modify('from_hex'); 584 585 my $str = shift; 586 my @r = @_; 587 588 # If called as a class method, initialize a new object. 589 590 $self = bless {}, $class unless $selfref; 591 592 if (my @parts = $class -> _hex_str_to_flt_lib_parts($str)) { 593 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts; 594 595 $self = $self->round(@r) 596 unless @r >= 2 && !defined($r[0]) && !defined($r[1]); 597 598 return $downgrade -> new($self -> bdstr(), @r) 599 if defined($downgrade) && $self -> is_int(); 600 return $self; 601 } 602 603 return $self -> bnan(@r); 604} 605 606sub from_oct { 607 my $self = shift; 608 my $selfref = ref $self; 609 my $class = $selfref || $self; 610 611 # Make "require" work. 612 613 $class -> import() if $IMPORT == 0; 614 615 # Don't modify constant (read-only) objects. 616 617 return $self if $selfref && $self->modify('from_oct'); 618 619 my $str = shift; 620 my @r = @_; 621 622 # If called as a class method, initialize a new object. 623 624 $self = bless {}, $class unless $selfref; 625 626 if (my @parts = $class -> _oct_str_to_flt_lib_parts($str)) { 627 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts; 628 629 $self = $self->round(@r) 630 unless @r >= 2 && !defined($r[0]) && !defined($r[1]); 631 632 return $downgrade -> new($self -> bdstr(), @r) 633 if defined($downgrade) && $self -> is_int(); 634 return $self; 635 } 636 637 return $self -> bnan(@r); 638} 639 640sub from_bin { 641 my $self = shift; 642 my $selfref = ref $self; 643 my $class = $selfref || $self; 644 645 # Make "require" work. 646 647 $class -> import() if $IMPORT == 0; 648 649 # Don't modify constant (read-only) objects. 650 651 return $self if $selfref && $self->modify('from_bin'); 652 653 my $str = shift; 654 my @r = @_; 655 656 # If called as a class method, initialize a new object. 657 658 $self = bless {}, $class unless $selfref; 659 660 if (my @parts = $class -> _bin_str_to_flt_lib_parts($str)) { 661 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts; 662 663 $self = $self->round(@r) 664 unless @r >= 2 && !defined($r[0]) && !defined($r[1]); 665 666 return $downgrade -> new($self -> bdstr(), @r) 667 if defined($downgrade) && $self -> is_int(); 668 return $self; 669 } 670 671 return $self -> bnan(@r); 672} 673 674sub from_ieee754 { 675 my $self = shift; 676 my $selfref = ref $self; 677 my $class = $selfref || $self; 678 679 # Make "require" work. 680 681 $class -> import() if $IMPORT == 0; 682 683 # Don't modify constant (read-only) objects. 684 685 return $self if $selfref && $self->modify('from_ieee754'); 686 687 my $in = shift; # input string (or raw bytes) 688 my $format = shift; # format ("binary32", "decimal64" etc.) 689 my $enc; # significand encoding (applies only to decimal) 690 my $k; # storage width in bits 691 my $b; # base 692 my @r = @_; # rounding parameters, if any 693 694 if ($format =~ /^binary(\d+)\z/) { 695 $k = $1; 696 $b = 2; 697 } elsif ($format =~ /^decimal(\d+)(dpd|bcd)?\z/) { 698 $k = $1; 699 $b = 10; 700 $enc = $2 || 'dpd'; # default is dencely-packed decimals (DPD) 701 } elsif ($format eq 'half') { 702 $k = 16; 703 $b = 2; 704 } elsif ($format eq 'single') { 705 $k = 32; 706 $b = 2; 707 } elsif ($format eq 'double') { 708 $k = 64; 709 $b = 2; 710 } elsif ($format eq 'quadruple') { 711 $k = 128; 712 $b = 2; 713 } elsif ($format eq 'octuple') { 714 $k = 256; 715 $b = 2; 716 } elsif ($format eq 'sexdecuple') { 717 $k = 512; 718 $b = 2; 719 } 720 721 if ($b == 2) { 722 723 # Get the parameters for this format. 724 725 my $p; # precision (in bits) 726 my $t; # number of bits in significand 727 my $w; # number of bits in exponent 728 729 if ($k == 16) { # binary16 (half-precision) 730 $p = 11; 731 $t = 10; 732 $w = 5; 733 } elsif ($k == 32) { # binary32 (single-precision) 734 $p = 24; 735 $t = 23; 736 $w = 8; 737 } elsif ($k == 64) { # binary64 (double-precision) 738 $p = 53; 739 $t = 52; 740 $w = 11; 741 } else { # binaryN (quadruple-precision and above) 742 if ($k < 128 || $k != 32 * sprintf('%.0f', $k / 32)) { 743 croak "Number of bits must be 16, 32, 64, or >= 128 and", 744 " a multiple of 32"; 745 } 746 $p = $k - sprintf('%.0f', 4 * log($k) / log(2)) + 13; 747 $t = $p - 1; 748 $w = $k - $t - 1; 749 } 750 751 # The maximum exponent, minimum exponent, and exponent bias. 752 753 my $emax = $class -> new(2) -> bpow($w - 1) -> bdec(); 754 my $emin = 1 - $emax; 755 my $bias = $emax; 756 757 # Undefined input. 758 759 unless (defined $in) { 760 carp("Input is undefined"); 761 return $self -> bzero(@r); 762 } 763 764 # Make sure input string is a string of zeros and ones. 765 766 my $len = CORE::length $in; 767 if (8 * $len == $k) { # bytes 768 $in = unpack "B*", $in; 769 } elsif (4 * $len == $k) { # hexadecimal 770 if ($in =~ /([^\da-f])/i) { 771 croak "Illegal hexadecimal digit '$1'"; 772 } 773 $in = unpack "B*", pack "H*", $in; 774 } elsif ($len == $k) { # bits 775 if ($in =~ /([^01])/) { 776 croak "Illegal binary digit '$1'"; 777 } 778 } else { 779 croak "Unknown input -- $in"; 780 } 781 782 # Split bit string into sign, exponent, and mantissa/significand. 783 784 my $sign = substr($in, 0, 1) eq '1' ? '-' : '+'; 785 my $expo = $class -> from_bin(substr($in, 1, $w)); 786 my $mant = $class -> from_bin(substr($in, $w + 1)); 787 788 my $x; 789 790 $expo = $expo -> bsub($bias); # subtract bias 791 792 if ($expo < $emin) { # zero and subnormals 793 if ($mant == 0) { # zero 794 $x = $class -> bzero(); 795 } else { # subnormals 796 # compute (1/$b)**(N) rather than ($b)**(-N) 797 $x = $class -> new("0.5"); # 1/$b 798 $x = $x -> bpow($bias + $t - 1) -> bmul($mant); 799 $x = $x -> bneg() if $sign eq '-'; 800 } 801 } 802 803 elsif ($expo > $emax) { # inf and nan 804 if ($mant == 0) { # inf 805 $x = $class -> binf($sign); 806 } else { # nan 807 $x = $class -> bnan(@r); 808 } 809 } 810 811 else { # normals 812 $mant = $class -> new(2) -> bpow($t) -> badd($mant); 813 if ($expo < $t) { 814 # compute (1/$b)**(N) rather than ($b)**(-N) 815 $x = $class -> new("0.5"); # 1/$b 816 $x = $x -> bpow($t - $expo) -> bmul($mant); 817 } else { 818 $x = $class -> new(2); 819 $x = $x -> bpow($expo - $t) -> bmul($mant); 820 } 821 $x = $x -> bneg() if $sign eq '-'; 822 } 823 824 if ($selfref) { 825 $self -> {sign} = $x -> {sign}; 826 $self -> {_m} = $x -> {_m}; 827 $self -> {_es} = $x -> {_es}; 828 $self -> {_e} = $x -> {_e}; 829 } else { 830 $self = $x; 831 } 832 833 return $downgrade -> new($self -> bdstr(), @r) 834 if defined($downgrade) && $self -> is_int(); 835 return $self -> round(@r); 836 } 837 838 croak("The format '$format' is not yet supported."); 839} 840 841sub bzero { 842 # create/assign '+0' 843 844 # Class::method(...) -> Class->method(...) 845 unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) || 846 $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i)) 847 { 848 #carp "Using ", (caller(0))[3], "() as a function is deprecated;", 849 # " use is as a method instead"; 850 unshift @_, __PACKAGE__; 851 } 852 853 my $self = shift; 854 my $selfref = ref $self; 855 my $class = $selfref || $self; 856 857 # Make "require" work. 858 859 $class -> import() if $IMPORT == 0; 860 861 # Don't modify constant (read-only) objects. 862 863 return $self if $selfref && $self->modify('bzero'); 864 865 # Get the rounding parameters, if any. 866 867 my @r = @_; 868 869 return $downgrade -> bzero(@r) if defined $downgrade; 870 871 # If called as a class method, initialize a new object. 872 873 $self = bless {}, $class unless $selfref; 874 875 $self -> {sign} = '+'; 876 $self -> {_m} = $LIB -> _zero(); 877 $self -> {_es} = '+'; 878 $self -> {_e} = $LIB -> _zero(); 879 880 # If rounding parameters are given as arguments, use them. If no rounding 881 # parameters are given, and if called as a class method initialize the new 882 # instance with the class variables. 883 884 #return $self -> round(@r); # this should work, but doesnt; fixme! 885 886 if (@r) { 887 if (@r >= 2 && defined($r[0]) && defined($r[1])) { 888 carp "can't specify both accuracy and precision"; 889 return $self -> bnan(); 890 } 891 $self->{accuracy} = $r[0]; 892 $self->{precision} = $r[1]; 893 } else { 894 unless($selfref) { 895 $self->{accuracy} = $class -> accuracy(); 896 $self->{precision} = $class -> precision(); 897 } 898 } 899 900 return $self; 901} 902 903sub bone { 904 # Create or assign '+1' (or -1 if given sign '-'). 905 906 # Class::method(...) -> Class->method(...) 907 unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) || 908 $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i)) 909 { 910 #carp "Using ", (caller(0))[3], "() as a function is deprecated;", 911 # " use is as a method instead"; 912 unshift @_, __PACKAGE__; 913 } 914 915 my $self = shift; 916 my $selfref = ref $self; 917 my $class = $selfref || $self; 918 919 # Make "require" work. 920 921 $class -> import() if $IMPORT == 0; 922 923 # Don't modify constant (read-only) objects. 924 925 return $self if $selfref && $self->modify('bone'); 926 927 return $downgrade -> bone(@_) if defined $downgrade; 928 929 # Get the sign. 930 931 my $sign = '+'; # default is to return +1 932 if (defined($_[0]) && $_[0] =~ /^\s*([+-])\s*$/) { 933 $sign = $1; 934 shift; 935 } 936 937 # Get the rounding parameters, if any. 938 939 my @r = @_; 940 941 # If called as a class method, initialize a new object. 942 943 $self = bless {}, $class unless $selfref; 944 945 $self -> {sign} = $sign; 946 $self -> {_m} = $LIB -> _one(); 947 $self -> {_es} = '+'; 948 $self -> {_e} = $LIB -> _zero(); 949 950 # If rounding parameters are given as arguments, use them. If no rounding 951 # parameters are given, and if called as a class method initialize the new 952 # instance with the class variables. 953 954 #return $self -> round(@r); # this should work, but doesnt; fixme! 955 956 if (@r) { 957 if (@r >= 2 && defined($r[0]) && defined($r[1])) { 958 carp "can't specify both accuracy and precision"; 959 return $self -> bnan(); 960 } 961 $self->{accuracy} = $_[0]; 962 $self->{precision} = $_[1]; 963 } else { 964 unless($selfref) { 965 $self->{accuracy} = $class -> accuracy(); 966 $self->{precision} = $class -> precision(); 967 } 968 } 969 970 return $self; 971} 972 973sub binf { 974 # create/assign a '+inf' or '-inf' 975 976 # Class::method(...) -> Class->method(...) 977 unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) || 978 $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i)) 979 { 980 #carp "Using ", (caller(0))[3], "() as a function is deprecated;", 981 # " use is as a method instead"; 982 unshift @_, __PACKAGE__; 983 } 984 985 my $self = shift; 986 my $selfref = ref $self; 987 my $class = $selfref || $self; 988 989 { 990 no strict 'refs'; 991 if (${"${class}::_trap_inf"}) { 992 croak("Tried to create +-inf in $class->binf()"); 993 } 994 } 995 996 # Make "require" work. 997 998 $class -> import() if $IMPORT == 0; 999 1000 # Don't modify constant (read-only) objects. 1001 1002 return $self if $selfref && $self->modify('binf'); 1003 1004 return $downgrade -> binf(@_) if $downgrade; 1005 1006 # Get the sign. 1007 1008 my $sign = '+'; # default is to return positive infinity 1009 if (defined($_[0]) && $_[0] =~ /^\s*([+-])(inf|$)/i) { 1010 $sign = $1; 1011 shift; 1012 } 1013 1014 # Get the rounding parameters, if any. 1015 1016 my @r = @_; 1017 1018 # If called as a class method, initialize a new object. 1019 1020 $self = bless {}, $class unless $selfref; 1021 1022 $self -> {sign} = $sign . 'inf'; 1023 $self -> {_m} = $LIB -> _zero(); 1024 $self -> {_es} = '+'; 1025 $self -> {_e} = $LIB -> _zero(); 1026 1027 # If rounding parameters are given as arguments, use them. If no rounding 1028 # parameters are given, and if called as a class method initialize the new 1029 # instance with the class variables. 1030 1031 #return $self -> round(@r); # this should work, but doesnt; fixme! 1032 1033 if (@r) { 1034 if (@r >= 2 && defined($r[0]) && defined($r[1])) { 1035 carp "can't specify both accuracy and precision"; 1036 return $self -> bnan(); 1037 } 1038 $self->{accuracy} = $r[0]; 1039 $self->{precision} = $r[1]; 1040 } else { 1041 unless($selfref) { 1042 $self->{accuracy} = $class -> accuracy(); 1043 $self->{precision} = $class -> precision(); 1044 } 1045 } 1046 1047 return $self; 1048} 1049 1050sub bnan { 1051 # create/assign a 'NaN' 1052 1053 # Class::method(...) -> Class->method(...) 1054 unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) || 1055 $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i)) 1056 { 1057 #carp "Using ", (caller(0))[3], "() as a function is deprecated;", 1058 # " use is as a method instead"; 1059 unshift @_, __PACKAGE__; 1060 } 1061 1062 my $self = shift; 1063 my $selfref = ref $self; 1064 my $class = $selfref || $self; 1065 1066 { 1067 no strict 'refs'; 1068 if (${"${class}::_trap_nan"}) { 1069 croak("Tried to create NaN in $class->bnan()"); 1070 } 1071 } 1072 1073 # Make "require" work. 1074 1075 $class -> import() if $IMPORT == 0; 1076 1077 # Don't modify constant (read-only) objects. 1078 1079 return $self if $selfref && $self->modify('bnan'); 1080 1081 return $downgrade -> bnan(@_) if defined $downgrade; 1082 1083 # Get the rounding parameters, if any. 1084 1085 my @r = @_; 1086 1087 # If called as a class method, initialize a new object. 1088 1089 $self = bless {}, $class unless $selfref; 1090 1091 $self -> {sign} = $nan; 1092 $self -> {_m} = $LIB -> _zero(); 1093 $self -> {_es} = '+'; 1094 $self -> {_e} = $LIB -> _zero(); 1095 1096 # If rounding parameters are given as arguments, use them. If no rounding 1097 # parameters are given, and if called as a class method initialize the new 1098 # instance with the class variables. 1099 1100 #return $self -> round(@r); # this should work, but doesnt; fixme! 1101 1102 if (@r) { 1103 if (@r >= 2 && defined($r[0]) && defined($r[1])) { 1104 carp "can't specify both accuracy and precision"; 1105 return $self -> bnan(); 1106 } 1107 $self->{accuracy} = $r[0]; 1108 $self->{precision} = $r[1]; 1109 } else { 1110 unless($selfref) { 1111 $self->{accuracy} = $class -> accuracy(); 1112 $self->{precision} = $class -> precision(); 1113 } 1114 } 1115 1116 return $self; 1117} 1118 1119sub bpi { 1120 1121 # Class::method(...) -> Class->method(...) 1122 unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) || 1123 $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i)) 1124 { 1125 #carp "Using ", (caller(0))[3], "() as a function is deprecated;", 1126 # " use is as a method instead"; 1127 unshift @_, __PACKAGE__; 1128 } 1129 1130 # Called as Argument list 1131 # --------- ------------- 1132 # Math::BigFloat->bpi() ("Math::BigFloat") 1133 # Math::BigFloat->bpi(10) ("Math::BigFloat", 10) 1134 # $x->bpi() ($x) 1135 # $x->bpi(10) ($x, 10) 1136 # Math::BigFloat::bpi() () 1137 # Math::BigFloat::bpi(10) (10) 1138 # 1139 # In ambiguous cases, we favour the OO-style, so the following case 1140 # 1141 # $n = Math::BigFloat->new("10"); 1142 # $x = Math::BigFloat->bpi($n); 1143 # 1144 # which gives an argument list with the single element $n, is resolved as 1145 # 1146 # $n->bpi(); 1147 1148 my $self = shift; 1149 my $selfref = ref $self; 1150 my $class = $selfref || $self; 1151 my @r = @_; # rounding paramters 1152 1153 # Make "require" work. 1154 1155 $class -> import() if $IMPORT == 0; 1156 1157 if ($selfref) { # bpi() called as an instance method 1158 return $self if $self -> modify('bpi'); 1159 } else { # bpi() called as a class method 1160 $self = bless {}, $class; # initialize new instance 1161 } 1162 1163 ($self, @r) = $self -> _find_round_parameters(@r); 1164 1165 # The accuracy, i.e., the number of digits. Pi has one digit before the 1166 # dot, so a precision of 4 digits is equivalent to an accuracy of 5 digits. 1167 1168 my $n = defined $r[0] ? $r[0] 1169 : defined $r[1] ? 1 - $r[1] 1170 : $self -> div_scale(); 1171 1172 my $rmode = defined $r[2] ? $r[2] : $self -> round_mode(); 1173 1174 my $pi; 1175 1176 if ($n <= 1000) { 1177 1178 # 75 x 14 = 1050 digits 1179 1180 my $all_digits = <<EOF; 1181314159265358979323846264338327950288419716939937510582097494459230781640628 1182620899862803482534211706798214808651328230664709384460955058223172535940812 1183848111745028410270193852110555964462294895493038196442881097566593344612847 1184564823378678316527120190914564856692346034861045432664821339360726024914127 1185372458700660631558817488152092096282925409171536436789259036001133053054882 1186046652138414695194151160943305727036575959195309218611738193261179310511854 1187807446237996274956735188575272489122793818301194912983367336244065664308602 1188139494639522473719070217986094370277053921717629317675238467481846766940513 1189200056812714526356082778577134275778960917363717872146844090122495343014654 1190958537105079227968925892354201995611212902196086403441815981362977477130996 1191051870721134999999837297804995105973173281609631859502445945534690830264252 1192230825334468503526193118817101000313783875288658753320838142061717766914730 1193359825349042875546873115956286388235378759375195778185778053217122680661300 1194192787661119590921642019893809525720106548586327886593615338182796823030195 1195EOF 1196 1197 # Should we round up? 1198 1199 my $round_up; 1200 1201 # From the string above, we need to extract the number of digits we 1202 # want plus extra characters for the newlines. 1203 1204 my $nchrs = $n + int($n / 75); 1205 1206 # Extract the digits we want. 1207 1208 my $digits = substr($all_digits, 0, $nchrs); 1209 1210 # Find out whether we should round up or down. Rounding is easy, since 1211 # pi is trancendental. With directed rounding, it doesn't matter what 1212 # the following digits are. With rounding to nearest, we only have to 1213 # look at one extra digit. 1214 1215 if ($rmode eq 'trunc') { 1216 $round_up = 0; 1217 } else { 1218 my $next_digit = substr($all_digits, $nchrs, 1); 1219 $round_up = $next_digit lt '5' ? 0 : 1; 1220 } 1221 1222 # Remove the newlines. 1223 1224 $digits =~ tr/0-9//cd; 1225 1226 # Now do the rounding. We could easily make the regex substitution 1227 # handle all cases, but we avoid using the regex engine when it is 1228 # simple to avoid it. 1229 1230 if ($round_up) { 1231 my $last_digit = substr($digits, -1, 1); 1232 if ($last_digit lt '9') { 1233 substr($digits, -1, 1) = ++$last_digit; 1234 } else { 1235 $digits =~ s{([0-8])(9+)$} 1236 { ($1 + 1) . ("0" x CORE::length($2)) }e; 1237 } 1238 } 1239 1240 # Convert to an object. 1241 1242 $pi = bless { 1243 sign => '+', 1244 _m => $LIB -> _new($digits), 1245 _es => '-', 1246 _e => $LIB -> _new($n - 1), 1247 }, $class; 1248 1249 } else { 1250 1251 # For large accuracy, the arctan formulas become very inefficient with 1252 # Math::BigFloat, so use Brent-Salamin (aka AGM or Gauss-Legendre). 1253 1254 # Use a few more digits in the intermediate computations. 1255 $n += 8; 1256 1257 $HALF = $class -> new($HALF) unless ref($HALF); 1258 my ($an, $bn, $tn, $pn) 1259 = ($class -> bone, $HALF -> copy() -> bsqrt($n), 1260 $HALF -> copy() -> bmul($HALF), $class -> bone); 1261 while ($pn < $n) { 1262 my $prev_an = $an -> copy(); 1263 $an = $an -> badd($bn) -> bmul($HALF, $n); 1264 $bn = $bn -> bmul($prev_an) -> bsqrt($n); 1265 $prev_an = $prev_an -> bsub($an); 1266 $tn = $tn -> bsub($pn * $prev_an * $prev_an); 1267 $pn = $pn -> badd($pn); 1268 } 1269 $an = $an -> badd($bn); 1270 $an = $an -> bmul($an, $n) -> bdiv(4 * $tn, $n); 1271 1272 $an = $an -> round(@r); 1273 $pi = $an; 1274 } 1275 1276 if (defined $r[0]) { 1277 $pi -> accuracy($r[0]); 1278 } elsif (defined $r[1]) { 1279 $pi -> precision($r[1]); 1280 } 1281 1282 for my $key (qw/ sign _m _es _e accuracy precision /) { 1283 $self -> {$key} = $pi -> {$key}; 1284 } 1285 1286 return $downgrade -> new($self -> bdstr(), @r) 1287 if defined($downgrade) && $self->is_int(); 1288 return $self; 1289} 1290 1291sub copy { 1292 my ($x, $class); 1293 if (ref($_[0])) { # $y = $x -> copy() 1294 $x = shift; 1295 $class = ref($x); 1296 } else { # $y = Math::BigInt -> copy($y) 1297 $class = shift; 1298 $x = shift; 1299 } 1300 1301 carp "Rounding is not supported for ", (caller(0))[3], "()" if @_; 1302 1303 my $copy = bless {}, $class; 1304 1305 $copy->{sign} = $x->{sign}; 1306 $copy->{_es} = $x->{_es}; 1307 $copy->{_m} = $LIB->_copy($x->{_m}); 1308 $copy->{_e} = $LIB->_copy($x->{_e}); 1309 $copy->{accuracy} = $x->{accuracy} if exists $x->{accuracy}; 1310 $copy->{precision} = $x->{precision} if exists $x->{precision}; 1311 1312 return $copy; 1313} 1314 1315sub as_int { 1316 # return copy as a bigint representation of this Math::BigFloat number 1317 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 1318 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 1319 1320 return $x -> copy() if $x -> isa("Math::BigInt"); 1321 1322 # Disable upgrading and downgrading. 1323 1324 require Math::BigInt; 1325 my $upg = Math::BigInt -> upgrade(); 1326 my $dng = Math::BigInt -> downgrade(); 1327 Math::BigInt -> upgrade(undef); 1328 Math::BigInt -> downgrade(undef); 1329 1330 # Copy the value. 1331 1332 my $y; 1333 if ($x -> is_inf()) { 1334 $y = Math::BigInt -> binf($x->sign()); 1335 } elsif ($x -> is_nan()) { 1336 $y = Math::BigInt -> bnan(); 1337 } else { 1338 $y = $LIB->_copy($x->{_m}); 1339 if ($x->{_es} eq '-') { # < 0 1340 $y = $LIB->_rsft($y, $x->{_e}, 10); 1341 } elsif (! $LIB->_is_zero($x->{_e})) { # > 0 1342 $y = $LIB->_lsft($y, $x->{_e}, 10); 1343 } 1344 $y = Math::BigInt->new($x->{sign} . $LIB->_str($y)); 1345 } 1346 1347 # Copy the remaining instance variables. 1348 1349 ($y->{accuracy}, $y->{precision}) = ($x->{accuracy}, $x->{precision}); 1350 1351 # Restore upgrading and downgrading. 1352 1353 Math::BigInt -> upgrade($upg); 1354 Math::BigInt -> downgrade($dng); 1355 1356 return $y; 1357} 1358 1359sub as_float { 1360 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 1361 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 1362 1363 return $x -> copy() if $x -> isa("Math::BigFloat"); 1364 1365 # Disable upgrading and downgrading. 1366 1367 my $upg = Math::BigFloat -> upgrade(); 1368 my $dng = Math::BigFloat -> downgrade(); 1369 Math::BigFloat -> upgrade(undef); 1370 Math::BigFloat -> downgrade(undef); 1371 1372 # Copy the value. 1373 1374 my $y = Math::BigFloat -> new($x); 1375 1376 # Copy the remaining instance variables. 1377 1378 ($y->{accuracy}, $y->{precision}) = ($x->{accuracy}, $x->{precision}); 1379 1380 # Restore upgrading and downgrading. 1381 1382 Math::BigFloat -> upgrade($upg); 1383 Math::BigFloat -> downgrade($dng); 1384 1385 return $y; 1386} 1387 1388sub as_rat { 1389 # return copy as a Math::BigRat representation of this Math::BigFloat 1390 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 1391 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 1392 1393 return $x -> copy() if $x -> isa("Math::BigRat"); 1394 1395 # Disable upgrading and downgrading. 1396 1397 require Math::BigRat; 1398 my $upg = Math::BigRat -> upgrade(); 1399 my $dng = Math::BigRat -> downgrade(); 1400 Math::BigRat -> upgrade(undef); 1401 Math::BigRat -> downgrade(undef); 1402 1403 # Copy the value. 1404 1405 my $y; 1406 if ($x -> is_inf()) { 1407 $y = Math::BigRat -> binf($x -> sign()); 1408 } elsif ($x -> is_nan()) { 1409 $y = Math::BigRat -> bnan(); 1410 } else { 1411 my @flt_parts = ($x->{sign}, $x->{_m}, $x->{_es}, $x->{_e}); 1412 my @rat_parts = $class -> _flt_lib_parts_to_rat_lib_parts(@flt_parts); 1413 $y = Math::BigRat -> new($rat_parts[0] . $LIB -> _str($rat_parts[1]) 1414 . '/' . $LIB -> _str($rat_parts[2])); 1415 } 1416 1417 # Copy the remaining instance variables. 1418 1419 ($y->{accuracy}, $y->{precision}) = ($x->{accuracy}, $x->{precision}); 1420 1421 # Restore upgrading and downgrading. 1422 1423 Math::BigRat -> upgrade($upg); 1424 Math::BigRat -> downgrade($dng); 1425 1426 return $y; 1427} 1428 1429############################################################################### 1430# Boolean methods 1431############################################################################### 1432 1433sub is_zero { 1434 # return true if arg (BFLOAT or num_str) is zero 1435 my (undef, $x) = ref($_[0]) ? (undef, @_) : objectify(1, @_); 1436 1437 ($x->{sign} eq '+' && $LIB->_is_zero($x->{_m})) ? 1 : 0; 1438} 1439 1440sub is_one { 1441 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given 1442 my (undef, $x, $sign) = ref($_[0]) ? (undef, @_) : objectify(1, @_); 1443 1444 $sign = '+' if !defined $sign || $sign ne '-'; 1445 1446 ($x->{sign} eq $sign && 1447 $LIB->_is_zero($x->{_e}) && 1448 $LIB->_is_one($x->{_m})) ? 1 : 0; 1449} 1450 1451sub is_odd { 1452 # return true if arg (BFLOAT or num_str) is odd or false if even 1453 my (undef, $x) = ref($_[0]) ? (undef, @_) : objectify(1, @_); 1454 1455 (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't 1456 ($LIB->_is_zero($x->{_e})) && 1457 ($LIB->_is_odd($x->{_m}))) ? 1 : 0; 1458} 1459 1460sub is_even { 1461 # return true if arg (BINT or num_str) is even or false if odd 1462 my (undef, $x) = ref($_[0]) ? (undef, @_) : objectify(1, @_); 1463 1464 (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't 1465 ($x->{_es} eq '+') && # 123.45 isn't 1466 ($LIB->_is_even($x->{_m}))) ? 1 : 0; # but 1200 is 1467} 1468 1469sub is_int { 1470 # return true if arg (BFLOAT or num_str) is an integer 1471 my (undef, $x) = ref($_[0]) ? (undef, @_) : objectify(1, @_); 1472 1473 (($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't 1474 ($x->{_es} eq '+')) ? 1 : 0; # 1e-1 => no integer 1475} 1476 1477############################################################################### 1478# Comparison methods 1479############################################################################### 1480 1481sub bcmp { 1482 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort) 1483 1484 # set up parameters 1485 my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 1486 ? (ref($_[0]), @_) 1487 : objectify(2, @_); 1488 1489 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 1490 1491 # Handle all 'nan' cases. 1492 1493 return if ($x->{sign} eq $nan) || ($y->{sign} eq $nan); 1494 1495 # Handle all '+inf' and '-inf' cases. 1496 1497 return 0 if ($x->{sign} eq '+inf' && $y->{sign} eq '+inf' || 1498 $x->{sign} eq '-inf' && $y->{sign} eq '-inf'); 1499 return +1 if $x->{sign} eq '+inf'; # x = +inf and y < +inf 1500 return -1 if $x->{sign} eq '-inf'; # x = -inf and y > -inf 1501 return -1 if $y->{sign} eq '+inf'; # x < +inf and y = +inf 1502 return +1 if $y->{sign} eq '-inf'; # x > -inf and y = -inf 1503 1504 # Handle all cases with opposite signs. 1505 1506 return +1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # also does 0 <=> -y 1507 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # also does -x <=> 0 1508 1509 # Handle all remaining zero cases. 1510 1511 my $xz = $x->is_zero(); 1512 my $yz = $y->is_zero(); 1513 return 0 if $xz && $yz; # 0 <=> 0 1514 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y 1515 return +1 if $yz && $x->{sign} eq '+'; # +x <=> 0 1516 1517 # Both arguments are now finite, non-zero numbers with the same sign. 1518 1519 my $cmp; 1520 1521 # The next step is to compare the exponents, but since each mantissa is an 1522 # integer of arbitrary value, the exponents must be normalized by the length 1523 # of the mantissas before we can compare them. 1524 1525 my $mxl = $LIB->_len($x->{_m}); 1526 my $myl = $LIB->_len($y->{_m}); 1527 1528 # If the mantissas have the same length, there is no point in normalizing 1529 # the exponents by the length of the mantissas, so treat that as a special 1530 # case. 1531 1532 if ($mxl == $myl) { 1533 1534 # First handle the two cases where the exponents have different signs. 1535 1536 if ($x->{_es} eq '+' && $y->{_es} eq '-') { 1537 $cmp = +1; 1538 } elsif ($x->{_es} eq '-' && $y->{_es} eq '+') { 1539 $cmp = -1; 1540 } 1541 1542 # Then handle the case where the exponents have the same sign. 1543 1544 else { 1545 $cmp = $LIB->_acmp($x->{_e}, $y->{_e}); 1546 $cmp = -$cmp if $x->{_es} eq '-'; 1547 } 1548 1549 # Adjust for the sign, which is the same for x and y, and bail out if 1550 # we're done. 1551 1552 $cmp = -$cmp if $x->{sign} eq '-'; # 124 > 123, but -124 < -123 1553 return $cmp if $cmp; 1554 1555 } 1556 1557 # We must normalize each exponent by the length of the corresponding 1558 # mantissa. Life is a lot easier if we first make both exponents 1559 # non-negative. We do this by adding the same positive value to both 1560 # exponent. This is safe, because when comparing the exponents, only the 1561 # relative difference is important. 1562 1563 my $ex; 1564 my $ey; 1565 1566 if ($x->{_es} eq '+') { 1567 1568 # If the exponent of x is >= 0 and the exponent of y is >= 0, there is 1569 # no need to do anything special. 1570 1571 if ($y->{_es} eq '+') { 1572 $ex = $LIB->_copy($x->{_e}); 1573 $ey = $LIB->_copy($y->{_e}); 1574 } 1575 1576 # If the exponent of x is >= 0 and the exponent of y is < 0, add the 1577 # absolute value of the exponent of y to both. 1578 1579 else { 1580 $ex = $LIB->_copy($x->{_e}); 1581 $ex = $LIB->_add($ex, $y->{_e}); # ex + |ey| 1582 $ey = $LIB->_zero(); # -ex + |ey| = 0 1583 } 1584 1585 } else { 1586 1587 # If the exponent of x is < 0 and the exponent of y is >= 0, add the 1588 # absolute value of the exponent of x to both. 1589 1590 if ($y->{_es} eq '+') { 1591 $ex = $LIB->_zero(); # -ex + |ex| = 0 1592 $ey = $LIB->_copy($y->{_e}); 1593 $ey = $LIB->_add($ey, $x->{_e}); # ey + |ex| 1594 } 1595 1596 # If the exponent of x is < 0 and the exponent of y is < 0, add the 1597 # absolute values of both exponents to both exponents. 1598 1599 else { 1600 $ex = $LIB->_copy($y->{_e}); # -ex + |ey| + |ex| = |ey| 1601 $ey = $LIB->_copy($x->{_e}); # -ey + |ex| + |ey| = |ex| 1602 } 1603 1604 } 1605 1606 # Now we can normalize the exponents by adding lengths of the mantissas. 1607 1608 $ex = $LIB->_add($ex, $LIB->_new($mxl)); 1609 $ey = $LIB->_add($ey, $LIB->_new($myl)); 1610 1611 # We're done if the exponents are different. 1612 1613 $cmp = $LIB->_acmp($ex, $ey); 1614 $cmp = -$cmp if $x->{sign} eq '-'; # 124 > 123, but -124 < -123 1615 return $cmp if $cmp; 1616 1617 # Compare the mantissas, but first normalize them by padding the shorter 1618 # mantissa with zeros (shift left) until it has the same length as the 1619 # longer mantissa. 1620 1621 my $mx = $x->{_m}; 1622 my $my = $y->{_m}; 1623 1624 if ($mxl > $myl) { 1625 $my = $LIB->_lsft($LIB->_copy($my), $LIB->_new($mxl - $myl), 10); 1626 } elsif ($mxl < $myl) { 1627 $mx = $LIB->_lsft($LIB->_copy($mx), $LIB->_new($myl - $mxl), 10); 1628 } 1629 1630 $cmp = $LIB->_acmp($mx, $my); 1631 $cmp = -$cmp if $x->{sign} eq '-'; # 124 > 123, but -124 < -123 1632 return $cmp; 1633 1634} 1635 1636sub bacmp { 1637 # Compares 2 values, ignoring their signs. 1638 # Returns one of undef, <0, =0, >0. (suitable for sort) 1639 1640 # set up parameters 1641 my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 1642 ? (ref($_[0]), @_) 1643 : objectify(2, @_); 1644 1645 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 1646 1647 # handle +-inf and NaN's 1648 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/) { 1649 return if (($x->{sign} eq $nan) || ($y->{sign} eq $nan)); 1650 return 0 if ($x->is_inf() && $y->is_inf()); 1651 return 1 if ($x->is_inf() && !$y->is_inf()); 1652 return -1; 1653 } 1654 1655 # shortcut 1656 my $xz = $x->is_zero(); 1657 my $yz = $y->is_zero(); 1658 return 0 if $xz && $yz; # 0 <=> 0 1659 return -1 if $xz && !$yz; # 0 <=> +y 1660 return 1 if $yz && !$xz; # +x <=> 0 1661 1662 # adjust so that exponents are equal 1663 my $lxm = $LIB->_len($x->{_m}); 1664 my $lym = $LIB->_len($y->{_m}); 1665 my ($xes, $yes) = (1, 1); 1666 $xes = -1 if $x->{_es} ne '+'; 1667 $yes = -1 if $y->{_es} ne '+'; 1668 # the numify somewhat limits our length, but makes it much faster 1669 my $lx = $lxm + $xes * $LIB->_num($x->{_e}); 1670 my $ly = $lym + $yes * $LIB->_num($y->{_e}); 1671 my $l = $lx - $ly; 1672 return $l <=> 0 if $l != 0; 1673 1674 # lengths (corrected by exponent) are equal 1675 # so make mantissa equal-length by padding with zero (shift left) 1676 my $diff = $lxm - $lym; 1677 my $xm = $x->{_m}; # not yet copy it 1678 my $ym = $y->{_m}; 1679 if ($diff > 0) { 1680 $ym = $LIB->_copy($y->{_m}); 1681 $ym = $LIB->_lsft($ym, $LIB->_new($diff), 10); 1682 } elsif ($diff < 0) { 1683 $xm = $LIB->_copy($x->{_m}); 1684 $xm = $LIB->_lsft($xm, $LIB->_new(-$diff), 10); 1685 } 1686 $LIB->_acmp($xm, $ym); 1687} 1688 1689############################################################################### 1690# Arithmetic methods 1691############################################################################### 1692 1693sub bneg { 1694 # (BINT or num_str) return BINT 1695 # negate number or make a negated number from string 1696 my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_); 1697 1698 return $x if $x->modify('bneg'); 1699 1700 return $x -> bnan(@r) if $x -> is_nan(); 1701 1702 # For +0 do not negate (to have always normalized +0). 1703 $x->{sign} =~ tr/+-/-+/ 1704 unless $x->{sign} eq '+' && $LIB->_is_zero($x->{_m}); 1705 1706 return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade) 1707 && ($x -> is_int() || $x -> is_inf() || $x -> is_nan()); 1708 return $x -> round(@r); 1709} 1710 1711sub bnorm { 1712 # bnorm() can't support rounding, because bround() and bfround() call 1713 # bnorm(), which would recurse indefinitely. 1714 1715 # adjust m and e so that m is smallest possible 1716 my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_); 1717 1718 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 1719 1720 # inf, nan etc 1721 if ($x->{sign} !~ /^[+-]$/) { 1722 return $downgrade -> new($x) if defined $downgrade; 1723 return $x; 1724 } 1725 1726 my $zeros = $LIB->_zeros($x->{_m}); # correct for trailing zeros 1727 if ($zeros != 0) { 1728 my $z = $LIB->_new($zeros); 1729 $x->{_m} = $LIB->_rsft($x->{_m}, $z, 10); 1730 if ($x->{_es} eq '-') { 1731 if ($LIB->_acmp($x->{_e}, $z) >= 0) { 1732 $x->{_e} = $LIB->_sub($x->{_e}, $z); 1733 $x->{_es} = '+' if $LIB->_is_zero($x->{_e}); 1734 } else { 1735 $x->{_e} = $LIB->_sub($LIB->_copy($z), $x->{_e}); 1736 $x->{_es} = '+'; 1737 } 1738 } else { 1739 $x->{_e} = $LIB->_add($x->{_e}, $z); 1740 } 1741 } else { 1742 # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing 1743 # zeros). So, for something like 0Ey, set y to 0, and -0 => +0 1744 if ($LIB->_is_zero($x->{_m})) { 1745 $x->{sign} = '+'; 1746 $x->{_es} = '+'; 1747 $x->{_e} = $LIB->_zero(); 1748 } 1749 } 1750 1751 return $downgrade -> new($x) 1752 if defined($downgrade) && $x->is_int(); 1753 return $x; 1754} 1755 1756sub binc { 1757 # increment arg by one 1758 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 1759 1760 return $x if $x->modify('binc'); 1761 1762 # Inf and NaN 1763 1764 return $x -> bnan(@r) if $x -> is_nan(); 1765 return $x -> binf($x->{sign}, @r) if $x -> is_inf(); 1766 1767 # Non-integer 1768 1769 if ($x->{_es} eq '-') { 1770 return $x->badd($class->bone(), @r); 1771 } 1772 1773 # If the exponent is non-zero, convert the internal representation, so that, 1774 # e.g., 12e+3 becomes 12000e+0 and we can easily increment the mantissa. 1775 1776 if (!$LIB->_is_zero($x->{_e})) { 1777 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # 1e2 => 100 1778 $x->{_e} = $LIB->_zero(); # normalize 1779 $x->{_es} = '+'; 1780 # we know that the last digit of $x will be '1' or '9', depending on the 1781 # sign 1782 } 1783 1784 # now $x->{_e} == 0 1785 if ($x->{sign} eq '+') { 1786 $x->{_m} = $LIB->_inc($x->{_m}); 1787 return $x->bnorm()->bround(@r); 1788 } elsif ($x->{sign} eq '-') { 1789 $x->{_m} = $LIB->_dec($x->{_m}); 1790 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # -1 +1 => -0 => +0 1791 return $x->bnorm()->bround(@r); 1792 } 1793 1794 return $downgrade -> new($x -> bdstr(), @r) 1795 if defined($downgrade) && $x -> is_int(); 1796 return $x; 1797} 1798 1799sub bdec { 1800 # decrement arg by one 1801 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 1802 1803 return $x if $x->modify('bdec'); 1804 1805 # Inf and NaN 1806 1807 return $x -> bnan(@r) if $x -> is_nan(); 1808 return $x -> binf($x->{sign}, @r) if $x -> is_inf(); 1809 1810 # Non-integer 1811 1812 if ($x->{_es} eq '-') { 1813 return $x->badd($class->bone('-'), @r); 1814 } 1815 1816 # If the exponent is non-zero, convert the internal representation, so that, 1817 # e.g., 12e+3 becomes 12000e+0 and we can easily increment the mantissa. 1818 1819 if (!$LIB->_is_zero($x->{_e})) { 1820 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # 1e2 => 100 1821 $x->{_e} = $LIB->_zero(); # normalize 1822 $x->{_es} = '+'; 1823 } 1824 1825 # now $x->{_e} == 0 1826 my $zero = $x->is_zero(); 1827 if (($x->{sign} eq '-') || $zero) { # x <= 0 1828 $x->{_m} = $LIB->_inc($x->{_m}); 1829 $x->{sign} = '-' if $zero; # 0 => 1 => -1 1830 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # -1 +1 => -0 => +0 1831 return $x->bnorm()->round(@r); 1832 } 1833 elsif ($x->{sign} eq '+') { # x > 0 1834 $x->{_m} = $LIB->_dec($x->{_m}); 1835 return $x->bnorm()->round(@r); 1836 } 1837 1838 return $downgrade -> new($x -> bdstr(), @r) 1839 if defined($downgrade) && $x -> is_int(); 1840 return $x -> round(@r); 1841} 1842 1843sub badd { 1844 # set up parameters 1845 my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 1846 ? (ref($_[0]), @_) 1847 : objectify(2, @_); 1848 1849 return $x if $x->modify('badd'); 1850 1851 # inf and NaN handling 1852 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/) { 1853 1854 # $x is NaN and/or $y is NaN 1855 if ($x->{sign} eq $nan || $y->{sign} eq $nan) { 1856 $x = $x->bnan(); 1857 } 1858 1859 # $x is Inf and $y is Inf 1860 elsif ($x->{sign} =~ /^[+-]inf$/ && $y->{sign} =~ /^[+-]inf$/) { 1861 # +Inf + +Inf or -Inf + -Inf => same, rest is NaN 1862 $x = $x->bnan() if $x->{sign} ne $y->{sign}; 1863 } 1864 1865 # +-inf + something => +-inf; something +-inf => +-inf 1866 elsif ($y->{sign} =~ /^[+-]inf$/) { 1867 $x->{sign} = $y->{sign}; 1868 } 1869 1870 return $downgrade -> new($x -> bdstr(), @r) if defined $downgrade; 1871 return $x -> round(@r); 1872 } 1873 1874 return $upgrade->badd($x, $y, @r) if defined $upgrade; 1875 1876 $r[3] = $y; # no push! 1877 1878 # for speed: no add for $x + 0 1879 if ($y->is_zero()) { 1880 $x = $x->round(@r); 1881 } 1882 1883 # for speed: no add for 0 + $y 1884 elsif ($x->is_zero()) { 1885 # make copy, clobbering up x (modify in place!) 1886 $x->{_e} = $LIB->_copy($y->{_e}); 1887 $x->{_es} = $y->{_es}; 1888 $x->{_m} = $LIB->_copy($y->{_m}); 1889 $x->{sign} = $y->{sign} || $nan; 1890 $x = $x->round(@r); 1891 } 1892 1893 # both $x and $y are non-zero 1894 else { 1895 1896 # take lower of the two e's and adapt m1 to it to match m2 1897 my $e = $y->{_e}; 1898 $e = $LIB->_zero() if !defined $e; # if no BFLOAT? 1899 $e = $LIB->_copy($e); # make copy (didn't do it yet) 1900 1901 my $es; 1902 1903 ($e, $es) = $LIB -> _ssub($e, $y->{_es} || '+', $x->{_e}, $x->{_es}); 1904 1905 my $add = $LIB->_copy($y->{_m}); 1906 1907 if ($es eq '-') { # < 0 1908 $x->{_m} = $LIB->_lsft($x->{_m}, $e, 10); 1909 ($x->{_e}, $x->{_es}) = $LIB -> _sadd($x->{_e}, $x->{_es}, $e, $es); 1910 } elsif (!$LIB->_is_zero($e)) { # > 0 1911 $add = $LIB->_lsft($add, $e, 10); 1912 } 1913 1914 # else: both e are the same, so just leave them 1915 1916 if ($x->{sign} eq $y->{sign}) { 1917 $x->{_m} = $LIB->_add($x->{_m}, $add); 1918 } else { 1919 ($x->{_m}, $x->{sign}) = 1920 $LIB -> _sadd($x->{_m}, $x->{sign}, $add, $y->{sign}); 1921 } 1922 1923 # delete trailing zeros, then round 1924 $x = $x->bnorm()->round(@r); 1925 } 1926 1927 return $downgrade -> new($x -> bdstr(), @r) 1928 if defined($downgrade) && $x -> is_int(); 1929 return $x; # rounding already done above 1930} 1931 1932sub bsub { 1933 # set up parameters 1934 my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 1935 ? (ref($_[0]), @_) 1936 : objectify(2, @_); 1937 1938 return $x if $x -> modify('bsub'); 1939 1940 if ($y -> is_zero()) { 1941 $x = $x -> round(@r); 1942 } else { 1943 1944 # To correctly handle the special case $x -> bsub($x), we note the sign 1945 # of $x, then flip the sign of $y, and if the sign of $x changed too, 1946 # then we know that $x and $y are the same object. 1947 1948 my $xsign = $x -> {sign}; 1949 $y -> {sign} =~ tr/+-/-+/; # does nothing for NaN 1950 if ($xsign ne $x -> {sign}) { 1951 # special case of $x -> bsub($x) results in 0 1952 if ($xsign =~ /^[+-]$/) { 1953 $x = $x -> bzero(@r); 1954 } else { 1955 $x = $x -> bnan(); # NaN, -inf, +inf 1956 } 1957 return $downgrade -> new($x -> bdstr(), @r) if defined $downgrade; 1958 return $x -> round(@r); 1959 } 1960 $x = $x -> badd($y, @r); # badd does not leave internal zeros 1961 $y -> {sign} =~ tr/+-/-+/; # reset $y (does nothing for NaN) 1962 } 1963 1964 return $downgrade -> new($x -> bdstr(), @r) 1965 if defined($downgrade) && ($x->is_int() || $x->is_inf() || $x->is_nan()); 1966 $x; # already rounded by badd() or no rounding 1967} 1968 1969sub bmul { 1970 # multiply two numbers 1971 1972 # set up parameters 1973 my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 1974 ? (ref($_[0]), @_) 1975 : objectify(2, @_); 1976 1977 return $x if $x->modify('bmul'); 1978 1979 return $x->bnan(@r) if ($x->{sign} eq $nan) || ($y->{sign} eq $nan); 1980 1981 # inf handling 1982 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/)) { 1983 return $x->bnan(@r) if $x->is_zero() || $y->is_zero(); 1984 # result will always be +-inf: 1985 # +inf * +/+inf => +inf, -inf * -/-inf => +inf 1986 # +inf * -/-inf => -inf, -inf * +/+inf => -inf 1987 return $x->binf(@r) if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/); 1988 return $x->binf(@r) if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/); 1989 return $x->binf('-', @r); 1990 } 1991 1992 return $upgrade->bmul($x, $y, @r) if defined $upgrade; 1993 1994 # aEb * cEd = (a*c)E(b+d) 1995 $x->{_m} = $LIB->_mul($x->{_m}, $y->{_m}); 1996 ($x->{_e}, $x->{_es}) 1997 = $LIB -> _sadd($x->{_e}, $x->{_es}, $y->{_e}, $y->{_es}); 1998 1999 $r[3] = $y; # no push! 2000 2001 # adjust sign: 2002 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+'; 2003 $x = $x->bnorm->round(@r); 2004 2005 return $downgrade -> new($x -> bdstr(), @r) 2006 if defined($downgrade) && ($x->is_int() || $x->is_inf() || $x->is_nan()); 2007 return $x; 2008} 2009 2010sub bmuladd { 2011 # multiply two numbers and add the third to the result 2012 2013 # set up parameters 2014 my ($class, $x, $y, $z, @r) 2015 = ref($_[0]) && ref($_[0]) eq ref($_[1]) && ref($_[1]) eq ref($_[2]) 2016 ? (ref($_[0]), @_) 2017 : objectify(3, @_); 2018 2019 return $x if $x->modify('bmuladd'); 2020 2021 return $x->bnan(@r) if (($x->{sign} eq $nan) || 2022 ($y->{sign} eq $nan) || 2023 ($z->{sign} eq $nan)); 2024 2025 # inf handling 2026 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/)) { 2027 return $x->bnan(@r) if $x->is_zero() || $y->is_zero(); 2028 # result will always be +-inf: 2029 # +inf * +/+inf => +inf, -inf * -/-inf => +inf 2030 # +inf * -/-inf => -inf, -inf * +/+inf => -inf 2031 return $x->binf(@r) if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/); 2032 return $x->binf(@r) if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/); 2033 return $x->binf('-', @r); 2034 } 2035 2036 # aEb * cEd = (a*c)E(b+d) 2037 $x->{_m} = $LIB->_mul($x->{_m}, $y->{_m}); 2038 ($x->{_e}, $x->{_es}) 2039 = $LIB -> _sadd($x->{_e}, $x->{_es}, $y->{_e}, $y->{_es}); 2040 2041 $r[3] = $y; # no push! 2042 2043 # adjust sign: 2044 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+'; 2045 2046 # z=inf handling (z=NaN handled above) 2047 if ($z->{sign} =~ /^[+-]inf$/) { 2048 $x->{sign} = $z->{sign}; 2049 return $downgrade -> new($x -> bdstr(), @r) if defined $downgrade; 2050 return $x -> round(@r); 2051 } 2052 2053 # take lower of the two e's and adapt m1 to it to match m2 2054 my $e = $z->{_e}; 2055 $e = $LIB->_zero() if !defined $e; # if no BFLOAT? 2056 $e = $LIB->_copy($e); # make copy (didn't do it yet) 2057 2058 my $es; 2059 2060 ($e, $es) = $LIB -> _ssub($e, $z->{_es} || '+', $x->{_e}, $x->{_es}); 2061 2062 my $add = $LIB->_copy($z->{_m}); 2063 2064 if ($es eq '-') # < 0 2065 { 2066 $x->{_m} = $LIB->_lsft($x->{_m}, $e, 10); 2067 ($x->{_e}, $x->{_es}) = $LIB -> _sadd($x->{_e}, $x->{_es}, $e, $es); 2068 } elsif (!$LIB->_is_zero($e)) # > 0 2069 { 2070 $add = $LIB->_lsft($add, $e, 10); 2071 } 2072 # else: both e are the same, so just leave them 2073 2074 if ($x->{sign} eq $z->{sign}) { 2075 # add 2076 $x->{_m} = $LIB->_add($x->{_m}, $add); 2077 } else { 2078 ($x->{_m}, $x->{sign}) = 2079 $LIB -> _sadd($x->{_m}, $x->{sign}, $add, $z->{sign}); 2080 } 2081 2082 # delete trailing zeros, then round 2083 $x = $x->bnorm()->round(@r); 2084 2085 return $downgrade -> new($x -> bdstr(), @r) 2086 if defined($downgrade) && ($x->is_int() || $x->is_inf() || $x->is_nan()); 2087 return $x; 2088} 2089 2090sub bdiv { 2091 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return 2092 # (BFLOAT, BFLOAT) (quo, rem) or BFLOAT (only quo) 2093 2094 # set up parameters 2095 my ($class, $x, $y, @r) = (ref($_[0]), @_); 2096 # objectify is costly, so avoid it 2097 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 2098 ($class, $x, $y, @r) = objectify(2, @_); 2099 } 2100 2101 return $x if $x->modify('bdiv'); 2102 2103 my $wantarray = wantarray; # call only once 2104 2105 # At least one argument is NaN. This is handled the same way as in 2106 # Math::BigInt -> bdiv(). 2107 2108 if ($x -> is_nan() || $y -> is_nan()) { 2109 return $wantarray ? ($x -> bnan(@r), $class -> bnan(@r)) 2110 : $x -> bnan(@r); 2111 } 2112 2113 # Divide by zero and modulo zero. This is handled the same way as in 2114 # Math::BigInt -> bdiv(). See the comment in the code for Math::BigInt -> 2115 # bdiv() for further details. 2116 2117 if ($y -> is_zero()) { 2118 my ($quo, $rem); 2119 if ($wantarray) { 2120 $rem = $x -> copy() -> round(@r); 2121 $rem = $downgrade -> new($rem, @r) 2122 if defined($downgrade) && $rem -> is_int(); 2123 } 2124 if ($x -> is_zero()) { 2125 $quo = $x -> bnan(@r); 2126 } else { 2127 $quo = $x -> binf($x -> {sign}, @r); 2128 } 2129 return $wantarray ? ($quo, $rem) : $quo; 2130 } 2131 2132 # Numerator (dividend) is +/-inf. This is handled the same way as in 2133 # Math::BigInt -> bdiv(). See the comment in the code for Math::BigInt -> 2134 # bdiv() for further details. 2135 2136 if ($x -> is_inf()) { 2137 my ($quo, $rem); 2138 $rem = $class -> bnan(@r) if $wantarray; 2139 if ($y -> is_inf()) { 2140 $quo = $x -> bnan(@r); 2141 } else { 2142 my $sign = $x -> bcmp(0) == $y -> bcmp(0) ? '+' : '-'; 2143 $quo = $x -> binf($sign, @r); 2144 } 2145 return $wantarray ? ($quo, $rem) : $quo; 2146 } 2147 2148 # Denominator (divisor) is +/-inf. This is handled the same way as in 2149 # Math::BigInt -> bdiv(), with one exception: In scalar context, 2150 # Math::BigFloat does true division (although rounded), not floored division 2151 # (F-division), so a finite number divided by +/-inf is always zero. See the 2152 # comment in the code for Math::BigInt -> bdiv() for further details. 2153 2154 if ($y -> is_inf()) { 2155 my ($quo, $rem); 2156 if ($wantarray) { 2157 if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) { 2158 $rem = $x -> copy() -> round(@r); 2159 $rem = $downgrade -> new($rem, @r) 2160 if defined($downgrade) && $rem -> is_int(); 2161 $quo = $x -> bzero(@r); 2162 } else { 2163 $rem = $class -> binf($y -> {sign}, @r); 2164 $quo = $x -> bone('-', @r); 2165 } 2166 return ($quo, $rem); 2167 } else { 2168 if ($y -> is_inf()) { 2169 if ($x -> is_nan() || $x -> is_inf()) { 2170 return $x -> bnan(@r); 2171 } else { 2172 return $x -> bzero(@r); 2173 } 2174 } 2175 } 2176 } 2177 2178 # At this point, both the numerator and denominator are finite numbers, and 2179 # the denominator (divisor) is non-zero. 2180 2181 # x == 0? 2182 if ($x->is_zero()) { 2183 my ($quo, $rem); 2184 $quo = $x->round(@r); 2185 $quo = $downgrade -> new($quo, @r) 2186 if defined($downgrade) && $quo -> is_int(); 2187 if ($wantarray) { 2188 $rem = $class -> bzero(@r); 2189 return $quo, $rem; 2190 } 2191 return $quo; 2192 } 2193 2194 # Division might return a value that we can not represent exactly, so 2195 # upgrade, if upgrading is enabled. 2196 2197 return $upgrade -> bdiv($x, $y, @r) 2198 if defined($upgrade) && !wantarray && !$LIB -> _is_one($y -> {_m}); 2199 2200 # we need to limit the accuracy to protect against overflow 2201 my $fallback = 0; 2202 my (@params, $scale); 2203 ($x, @params) = $x->_find_round_parameters($r[0], $r[1], $r[2], $y); 2204 2205 return $x -> round(@r) if $x->is_nan(); # error in _find_round_parameters? 2206 2207 # no rounding at all, so must use fallback 2208 if (scalar @params == 0) { 2209 # simulate old behaviour 2210 $params[0] = $class->div_scale(); # and round to it as accuracy 2211 $scale = $params[0]+4; # at least four more for proper round 2212 $params[2] = $r[2]; # round mode by caller or undef 2213 $fallback = 1; # to clear a/p afterwards 2214 } else { 2215 # the 4 below is empirical, and there might be cases where it is not 2216 # enough... 2217 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 2218 } 2219 2220 my $rem; 2221 $rem = $class -> bzero() if wantarray; 2222 2223 $y = $class->new($y) unless $y->isa('Math::BigFloat'); 2224 2225 my $lx = $LIB -> _len($x->{_m}); 2226 my $ly = $LIB -> _len($y->{_m}); 2227 $scale = $lx if $lx > $scale; 2228 $scale = $ly if $ly > $scale; 2229 my $diff = $ly - $lx; 2230 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx! 2231 2232 # check that $y is not 1 nor -1 and cache the result: 2233 my $y_not_one = !($LIB->_is_zero($y->{_e}) && $LIB->_is_one($y->{_m})); 2234 2235 # flipping the sign of $y will also flip the sign of $x for the special 2236 # case of $x->bsub($x); so we can catch it below: 2237 my $xsign = $x->{sign}; 2238 $y->{sign} =~ tr/+-/-+/; 2239 2240 if ($xsign ne $x->{sign}) { 2241 # special case of $x /= $x results in 1 2242 $x = $x->bone(); # "fixes" also sign of $y, since $x is $y 2243 } else { 2244 # correct $y's sign again 2245 $y->{sign} =~ tr/+-/-+/; 2246 # continue with normal div code: 2247 2248 # make copy of $x in case of list context for later remainder 2249 # calculation 2250 if (wantarray && $y_not_one) { 2251 $rem = $x->copy(); 2252 } 2253 2254 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+'; 2255 2256 # check for / +-1 (+/- 1E0) 2257 if ($y_not_one) { 2258 # promote Math::BigInt and its subclasses (except when already a 2259 # Math::BigFloat) 2260 $y = $class->new($y) unless $y->isa('Math::BigFloat'); 2261 2262 # calculate the result to $scale digits and then round it 2263 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d) 2264 $x->{_m} = $LIB->_lsft($x->{_m}, $LIB->_new($scale), 10); 2265 $x->{_m} = $LIB->_div($x->{_m}, $y->{_m}); # a/c 2266 2267 # correct exponent of $x 2268 ($x->{_e}, $x->{_es}) 2269 = $LIB -> _ssub($x->{_e}, $x->{_es}, $y->{_e}, $y->{_es}); 2270 # correct for 10**scale 2271 ($x->{_e}, $x->{_es}) 2272 = $LIB -> _ssub($x->{_e}, $x->{_es}, $LIB->_new($scale), '+'); 2273 $x = $x->bnorm(); # remove trailing 0's 2274 } 2275 } # end else $x != $y 2276 2277 # shortcut to not run through _find_round_parameters again 2278 if (defined $params[0]) { 2279 $x->{accuracy} = undef; # clear before round 2280 $x = $x->bround($params[0], $params[2]); # then round accordingly 2281 } else { 2282 $x->{precision} = undef; # clear before round 2283 $x = $x->bfround($params[1], $params[2]); # then round accordingly 2284 } 2285 if ($fallback) { 2286 # clear a/p after round, since user did not request it 2287 $x->{accuracy} = undef; 2288 $x->{precision} = undef; 2289 } 2290 2291 if (wantarray) { 2292 if ($y_not_one) { 2293 $x = $x -> bfloor(); 2294 $rem = $rem->bmod($y, @params); # copy already done 2295 } 2296 if ($fallback) { 2297 # clear a/p after round, since user did not request it 2298 $rem->{accuracy} = undef; 2299 $rem->{precision} = undef; 2300 } 2301 $x = $downgrade -> new($x -> bdstr(), @r) 2302 if defined($downgrade) && $x -> is_int(); 2303 $rem = $downgrade -> new($rem -> bdstr(), @r) 2304 if defined($downgrade) && $rem -> is_int(); 2305 return ($x, $rem); 2306 } 2307 2308 $x = $downgrade -> new($x, @r) 2309 if defined($downgrade) && $x -> is_int(); 2310 $x; # rounding already done above 2311} 2312 2313sub bmod { 2314 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return remainder 2315 2316 # set up parameters 2317 my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 2318 ? (ref($_[0]), @_) 2319 : objectify(2, @_); 2320 2321 return $x if $x->modify('bmod'); 2322 2323 # At least one argument is NaN. This is handled the same way as in 2324 # Math::BigInt -> bmod(). 2325 2326 return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan(); 2327 2328 # Modulo zero. This is handled the same way as in Math::BigInt -> bmod(). 2329 2330 if ($y -> is_zero()) { 2331 return $x -> round(@r); 2332 } 2333 2334 # Numerator (dividend) is +/-inf. This is handled the same way as in 2335 # Math::BigInt -> bmod(). 2336 2337 if ($x -> is_inf()) { 2338 return $x -> bnan(@r); 2339 } 2340 2341 # Denominator (divisor) is +/-inf. This is handled the same way as in 2342 # Math::BigInt -> bmod(). 2343 2344 if ($y -> is_inf()) { 2345 if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) { 2346 return $x -> round(@r); 2347 } else { 2348 return $x -> binf($y -> sign(), @r); 2349 } 2350 } 2351 2352 return $x->bzero(@r) if $x->is_zero() 2353 || ($x->is_int() && 2354 # check that $y == +1 or $y == -1: 2355 ($LIB->_is_zero($y->{_e}) && $LIB->_is_one($y->{_m}))); 2356 2357 my $cmp = $x->bacmp($y); # equal or $x < $y? 2358 if ($cmp == 0) { # $x == $y => result 0 2359 return $x -> bzero(@r); 2360 } 2361 2362 # only $y of the operands negative? 2363 my $neg = $x->{sign} ne $y->{sign} ? 1 : 0; 2364 2365 $x->{sign} = $y->{sign}; # calc sign first 2366 if ($cmp < 0 && $neg == 0) { # $x < $y => result $x 2367 return $x -> round(@r); 2368 } 2369 2370 my $ym = $LIB->_copy($y->{_m}); 2371 2372 # 2e1 => 20 2373 $ym = $LIB->_lsft($ym, $y->{_e}, 10) 2374 if $y->{_es} eq '+' && !$LIB->_is_zero($y->{_e}); 2375 2376 # if $y has digits after dot 2377 my $shifty = 0; # correct _e of $x by this 2378 if ($y->{_es} eq '-') # has digits after dot 2379 { 2380 # 123 % 2.5 => 1230 % 25 => 5 => 0.5 2381 $shifty = $LIB->_num($y->{_e}); # no more digits after dot 2382 # 123 => 1230, $y->{_m} is already 25 2383 $x->{_m} = $LIB->_lsft($x->{_m}, $y->{_e}, 10); 2384 } 2385 # $ym is now mantissa of $y based on exponent 0 2386 2387 my $shiftx = 0; # correct _e of $x by this 2388 if ($x->{_es} eq '-') # has digits after dot 2389 { 2390 # 123.4 % 20 => 1234 % 200 2391 $shiftx = $LIB->_num($x->{_e}); # no more digits after dot 2392 $ym = $LIB->_lsft($ym, $x->{_e}, 10); # 123 => 1230 2393 } 2394 # 123e1 % 20 => 1230 % 20 2395 if ($x->{_es} eq '+' && !$LIB->_is_zero($x->{_e})) { 2396 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # es => '+' here 2397 } 2398 2399 $x->{_e} = $LIB->_new($shiftx); 2400 $x->{_es} = '+'; 2401 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0; 2402 $x->{_e} = $LIB->_add($x->{_e}, $LIB->_new($shifty)) if $shifty != 0; 2403 2404 # now mantissas are equalized, exponent of $x is adjusted, so calc result 2405 2406 $x->{_m} = $LIB->_mod($x->{_m}, $ym); 2407 2408 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # fix sign for -0 2409 $x = $x->bnorm(); 2410 2411 # if one of them negative => correct in place 2412 if ($neg != 0 && ! $x -> is_zero()) { 2413 my $r = $y - $x; 2414 $x->{_m} = $r->{_m}; 2415 $x->{_e} = $r->{_e}; 2416 $x->{_es} = $r->{_es}; 2417 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # fix sign for -0 2418 $x = $x->bnorm(); 2419 } 2420 2421 $x = $x->round($r[0], $r[1], $r[2], $y); 2422 return $downgrade -> new($x -> bdstr(), @r) 2423 if defined($downgrade) && ($x->is_int() || $x->is_inf() || $x->is_nan()); 2424 return $x; 2425} 2426 2427sub bmodpow { 2428 # takes a very large number to a very large exponent in a given very 2429 # large modulus, quickly, thanks to binary exponentiation. Supports 2430 # negative exponents. 2431 my ($class, $num, $exp, $mod, @r) 2432 = ref($_[0]) && ref($_[0]) eq ref($_[1]) && ref($_[1]) eq ref($_[2]) 2433 ? (ref($_[0]), @_) 2434 : objectify(3, @_); 2435 2436 return $num if $num->modify('bmodpow'); 2437 2438 return $num -> bnan(@r) 2439 if $mod->is_nan() || $exp->is_nan() || $mod->is_nan(); 2440 2441 # check modulus for valid values 2442 return $num->bnan(@r) if $mod->{sign} ne '+' || $mod->is_zero(); 2443 2444 # check exponent for valid values 2445 if ($exp->{sign} =~ /\w/) { 2446 # i.e., if it's NaN, +inf, or -inf... 2447 return $num->bnan(@r); 2448 } 2449 2450 $num = $num->bmodinv($mod, @r) if $exp->{sign} eq '-'; 2451 2452 # check num for valid values (also NaN if there was no inverse but $exp < 0) 2453 return $num->bnan(@r) if $num->{sign} !~ /^[+-]$/; 2454 2455 # $mod is positive, sign on $exp is ignored, result also positive 2456 2457 # XXX TODO: speed it up when all three numbers are integers 2458 $num = $num->bpow($exp)->bmod($mod); 2459 2460 return $downgrade -> new($num -> bdstr(), @r) if defined($downgrade) 2461 && ($num->is_int() || $num->is_inf() || $num->is_nan()); 2462 return $num -> round(@r); 2463} 2464 2465sub bpow { 2466 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT 2467 # compute power of two numbers, second arg is used as integer 2468 # modifies first argument 2469 2470 # set up parameters 2471 my ($class, $x, $y, $a, $p, $r) = (ref($_[0]), @_); 2472 # objectify is costly, so avoid it 2473 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) { 2474 ($class, $x, $y, $a, $p, $r) = objectify(2, @_); 2475 } 2476 2477 return $x if $x -> modify('bpow'); 2478 2479 # $x and/or $y is a NaN 2480 return $x -> bnan() if $x -> is_nan() || $y -> is_nan(); 2481 2482 # $x and/or $y is a +/-Inf 2483 if ($x -> is_inf("-")) { 2484 return $x -> bzero() if $y -> is_negative(); 2485 return $x -> bnan() if $y -> is_zero(); 2486 return $x if $y -> is_odd(); 2487 return $x -> bneg(); 2488 } elsif ($x -> is_inf("+")) { 2489 return $x -> bzero() if $y -> is_negative(); 2490 return $x -> bnan() if $y -> is_zero(); 2491 return $x; 2492 } elsif ($y -> is_inf("-")) { 2493 return $x -> bnan() if $x -> is_one("-"); 2494 return $x -> binf("+") if $x > -1 && $x < 1; 2495 return $x -> bone() if $x -> is_one("+"); 2496 return $x -> bzero(); 2497 } elsif ($y -> is_inf("+")) { 2498 return $x -> bnan() if $x -> is_one("-"); 2499 return $x -> bzero() if $x > -1 && $x < 1; 2500 return $x -> bone() if $x -> is_one("+"); 2501 return $x -> binf("+"); 2502 } 2503 2504 if ($x -> is_zero()) { 2505 return $x -> bone() if $y -> is_zero(); 2506 return $x -> binf() if $y -> is_negative(); 2507 return $x; 2508 } 2509 2510 # We don't support complex numbers, so upgrade or return NaN. 2511 2512 if ($x -> is_negative() && !$y -> is_int()) { 2513 return $upgrade -> bpow($x, $y, $a, $p, $r) if defined $upgrade; 2514 return $x -> bnan(); 2515 } 2516 2517 if ($x -> is_one("+") || $y -> is_one()) { 2518 return $x; 2519 } 2520 2521 if ($x -> is_one("-")) { 2522 return $x if $y -> is_odd(); 2523 return $x -> bneg(); 2524 } 2525 2526 return $x -> _pow($y, $a, $p, $r) if !$y -> is_int(); 2527 2528 my $y1 = $y -> as_int()->{value}; # make MBI part 2529 2530 my $new_sign = '+'; 2531 $new_sign = $LIB -> _is_odd($y1) ? '-' : '+' if $x->{sign} ne '+'; 2532 2533 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster) 2534 $x->{_m} = $LIB -> _pow($x->{_m}, $y1); 2535 $x->{_e} = $LIB -> _mul($x->{_e}, $y1); 2536 2537 $x->{sign} = $new_sign; 2538 $x = $x -> bnorm(); 2539 2540 # x ** (-y) = 1 / (x ** y) 2541 2542 if ($y->{sign} eq '-') { 2543 # modify $x in place! 2544 my $z = $x -> copy(); 2545 $x = $x -> bone(); 2546 # round in one go (might ignore y's A!) 2547 return scalar $x -> bdiv($z, $a, $p, $r); 2548 } 2549 2550 $x = $x -> round($a, $p, $r, $y); 2551 2552 return $downgrade -> new($x) 2553 if defined($downgrade) && ($x->is_int() || $x->is_inf() || $x->is_nan()); 2554 return $x; 2555} 2556 2557sub binv { 2558 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 2559 2560 return $x if $x->modify('binv'); 2561 2562 my $inv = $class -> bdiv($class -> bone(), $x, @r); 2563 2564 return $downgrade -> new($inv, @r) if defined($downgrade) 2565 && ($inv -> is_int() || $inv -> is_inf() || $inv -> is_nan()); 2566 2567 for my $key (qw/ sign _m _es _e /) { 2568 $x -> {$key} = $inv -> {$key}; 2569 } 2570 2571 $x; 2572} 2573 2574sub blog { 2575 # Return the logarithm of the operand. If a second operand is defined, that 2576 # value is used as the base, otherwise the base is assumed to be Euler's 2577 # constant. 2578 2579 my ($class, $x, $base, @r); 2580 2581 # Only objectify the base if it is defined, since an undefined base, as in 2582 # $x->blog() or $x->blog(undef) signals that the base is Euler's number = 2583 # 2.718281828... 2584 2585 if (!ref($_[0]) && $_[0] =~ /^[A-Za-z]|::/) { 2586 # E.g., Math::BigFloat->blog(256, 2) 2587 ($class, $x, $base, @r) = 2588 defined $_[2] ? objectify(2, @_) : objectify(1, @_); 2589 } else { 2590 # E.g., $x->blog(2) or the deprecated Math::BigFloat::blog(256, 2) 2591 ($class, $x, $base, @r) = 2592 defined $_[1] ? objectify(2, @_) : objectify(1, @_); 2593 } 2594 2595 return $x if $x->modify('blog'); 2596 2597 # Handle all exception cases and all trivial cases. I have used Wolfram 2598 # Alpha (http://www.wolframalpha.com) as the reference for these cases. 2599 2600 return $x -> bnan(@r) if $x -> is_nan(); 2601 2602 if (defined $base) { 2603 $base = $class -> new($base) 2604 unless defined(blessed($base)) && $base -> isa(__PACKAGE__); 2605 if ($base -> is_nan() || $base -> is_one()) { 2606 return $x -> bnan(@r); 2607 } elsif ($base -> is_inf() || $base -> is_zero()) { 2608 return $x -> bnan(@r) if $x -> is_inf() || $x -> is_zero(); 2609 return $x -> bzero(@r); 2610 } elsif ($base -> is_negative()) { # -inf < base < 0 2611 return $x -> bzero(@r) if $x -> is_one(); # x = 1 2612 return $x -> bone('+', @r) if $x == $base; # x = base 2613 # we can't handle these cases, so upgrade, if we can 2614 return $upgrade -> blog($x, $base, @r) if defined $upgrade; 2615 return $x -> bnan(@r); 2616 } 2617 return $x -> bone(@r) if $x == $base; # 0 < base && 0 < x < inf 2618 } 2619 2620 if ($x -> is_inf()) { # x = +/-inf 2621 my $sign = defined($base) && $base < 1 ? '-' : '+'; 2622 return $x -> binf($sign, @r); 2623 } elsif ($x -> is_neg()) { # -inf < x < 0 2624 return $upgrade -> blog($x, $base, @r) if defined $upgrade; 2625 return $x -> bnan(@r); 2626 } elsif ($x -> is_one()) { # x = 1 2627 return $x -> bzero(@r); 2628 } elsif ($x -> is_zero()) { # x = 0 2629 my $sign = defined($base) && $base < 1 ? '+' : '-'; 2630 return $x -> binf($sign, @r); 2631 } 2632 2633 # we need to limit the accuracy to protect against overflow 2634 my $fallback = 0; 2635 my ($scale, @params); 2636 ($x, @params) = $x->_find_round_parameters(@r); 2637 2638 # no rounding at all, so must use fallback 2639 if (scalar @params == 0) { 2640 # simulate old behaviour 2641 $params[0] = $class->div_scale(); # and round to it as accuracy 2642 $params[1] = undef; # P = undef 2643 $scale = $params[0]+4; # at least four more for proper round 2644 $params[2] = $r[2]; # round mode by caller or undef 2645 $fallback = 1; # to clear a/p afterwards 2646 } else { 2647 # the 4 below is empirical, and there might be cases where it is not 2648 # enough... 2649 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 2650 } 2651 2652 # When user set globals, they would interfere with our calculation, so 2653 # disable them and later re-enable them. 2654 2655 my $ab = $class -> accuracy(); 2656 my $pb = $class -> precision(); 2657 $class -> accuracy(undef); 2658 $class -> precision(undef); 2659 2660 # Disabling upgrading and downgrading is no longer necessary to avoid an 2661 # infinite recursion, but it avoids unnecessary upgrading and downgrading in 2662 # the intermediate computations. 2663 2664 my $upg = $class -> upgrade(); 2665 my $dng = $class -> downgrade(); 2666 $class -> upgrade(undef); 2667 $class -> downgrade(undef); 2668 2669 # We also need to disable any set A or P on $x (_find_round_parameters took 2670 # them already into account), since these would interfere, too. 2671 2672 $x->{accuracy} = undef; 2673 $x->{precision} = undef; 2674 2675 my $done = 0; 2676 2677 # If both $x and $base are integers, try to calculate an integer result 2678 # first. This is very fast, and if the exact result was found, we are done. 2679 2680 if (defined($base) && $base -> is_int() && $x -> is_int()) { 2681 my $x_lib = $LIB -> _new($x -> bdstr()); 2682 my $b_lib = $LIB -> _new($base -> bdstr()); 2683 ($x_lib, my $exact) = $LIB -> _log_int($x_lib, $b_lib); 2684 if ($exact) { 2685 $x->{_m} = $x_lib; 2686 $x->{_e} = $LIB -> _zero(); 2687 $x = $x -> bnorm(); 2688 $done = 1; 2689 } 2690 } 2691 2692 # If the integer result was not accurate, compute the natural logarithm 2693 # log($x) (using reduction by 10 and possibly also by 2), and if a 2694 # different base was requested, convert the result with log($x)/log($base). 2695 2696 unless ($done) { 2697 $x = $x -> _log_10($scale); 2698 if (defined $base) { 2699 # log_b(x) = ln(x) / ln(b), so compute ln(b) 2700 my $base_log_e = $base -> copy() -> _log_10($scale); 2701 $x = $x -> bdiv($base_log_e, $scale); 2702 } 2703 } 2704 2705 # shortcut to not run through _find_round_parameters again 2706 2707 if (defined $params[0]) { 2708 $x = $x -> bround($params[0], $params[2]); # then round accordingly 2709 } else { 2710 $x = $x -> bfround($params[1], $params[2]); # then round accordingly 2711 } 2712 if ($fallback) { 2713 # clear a/p after round, since user did not request it 2714 $x->{accuracy} = undef; 2715 $x->{precision} = undef; 2716 } 2717 2718 # Restore globals. We need to do it like this, because setting one 2719 # undefines the other. 2720 2721 if (defined $ab) { 2722 $class -> accuracy($ab); 2723 } else { 2724 $class -> precision($pb); 2725 } 2726 2727 $class -> upgrade($upg); 2728 $class -> downgrade($dng); 2729 2730 return $downgrade -> new($x -> bdstr(), @r) 2731 if defined($downgrade) && $x -> is_int(); 2732 return $x; 2733} 2734 2735sub bexp { 2736 # Calculate e ** X (Euler's number to the power of X) 2737 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 2738 2739 return $x if $x -> modify('bexp'); 2740 2741 return $x -> bnan(@r) if $x -> is_nan(); 2742 return $x -> binf(@r) if $x->{sign} eq '+inf'; 2743 return $x -> bzero(@r) if $x->{sign} eq '-inf'; 2744 2745 # Get the rounding parameters, if any. 2746 2747 my $fallback = 0; 2748 my ($scale, @params); 2749 ($x, @params) = $x -> _find_round_parameters(@r); 2750 2751 # Error in _find_round_parameters? 2752 2753 return $x -> bnan(@r) if $x->{sign} eq 'NaN'; 2754 2755 return $x -> bone(@r) if $x -> is_zero(); 2756 2757 # If no rounding parameters are give, use fallback. 2758 2759 if (!@params) { 2760 $params[0] = $class -> div_scale(); # fallback accuracy 2761 $params[1] = undef; # no precision 2762 $params[2] = $r[2]; # rounding mode 2763 $scale = $params[0]; 2764 $fallback = 1; # to clear a/p afterwards 2765 } else { 2766 if (defined($params[0])) { 2767 $scale = $params[0]; 2768 } else { 2769 # We perform the computations below using accuracy only, not 2770 # precision, so when precision is given, we need to "convert" this 2771 # to accuracy. To do that, we need to know, at least approximately, 2772 # how many digits there will be in the final result. 2773 # 2774 # log10(exp($x)) = log(exp($x)) / log(10) = $x / log(10) 2775 2776 #$scale = 1 + int(log($ms) / log(10) + $es) - $params[1]; 2777 my $ndig = $x -> numify() / log(10); 2778 $scale = 1 + int($ndig) - $params[1]; 2779 } 2780 } 2781 2782 # Add extra digits to reduce the consequence of round-off errors in the 2783 # intermediate computations. 2784 2785 $scale += 4; 2786 2787 if (!$x -> isa('Math::BigFloat')) { 2788 $x = Math::BigFloat -> new($x); 2789 $class = ref($x); 2790 } 2791 2792 # When user set globals, they would interfere with our calculation, so 2793 # disable them and later re-enable them. 2794 2795 my $ab = $class -> accuracy(); 2796 my $pb = $class -> precision(); 2797 $class -> accuracy(undef); 2798 $class -> precision(undef); 2799 2800 # Disabling upgrading and downgrading is no longer necessary to avoid an 2801 # infinite recursion, but it avoids unnecessary upgrading and downgrading in 2802 # the intermediate computations. 2803 2804 my $upg = $class -> upgrade(); 2805 my $dng = $class -> downgrade(); 2806 $class -> upgrade(undef); 2807 $class -> downgrade(undef); 2808 2809 # We also need to disable any set A or P on $x (_find_round_parameters took 2810 # them already into account), since these would interfere, too. 2811 2812 $x->{accuracy} = undef; 2813 $x->{precision} = undef; 2814 2815 my $x_orig = $x -> copy(); 2816 2817 # We use the following Taylor series: 2818 2819 # x x^2 x^3 x^4 2820 # e = 1 + --- + --- + --- + --- ... 2821 # 1! 2! 3! 4! 2822 2823 # The difference for each term is X and N, which would result in: 2824 # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term 2825 2826 # But it is faster to compute exp(1) and then raising it to the 2827 # given power, esp. if $x is really big and an integer because: 2828 2829 # * The numerator is always 1, making the computation faster 2830 # * the series converges faster in the case of x == 1 2831 # * We can also easily check when we have reached our limit: when the 2832 # term to be added is smaller than "1E$scale", we can stop - f.i. 2833 # scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5. 2834 # * we can compute the *exact* result by simulating bigrat math: 2835 2836 # 1 1 gcd(3, 4) = 1 1*24 + 1*6 5 2837 # - + - = ---------- = -- 2838 # 6 24 6*24 24 2839 2840 # We do not compute the gcd() here, but simple do: 2841 # 1 1 1*24 + 1*6 30 2842 # - + - = --------- = -- 2843 # 6 24 6*24 144 2844 2845 # In general: 2846 # a c a*d + c*b and note that c is always 1 and d = (b*f) 2847 # - + - = --------- 2848 # b d b*d 2849 2850 # This leads to: which can be reduced by b to: 2851 # a 1 a*b*f + b a*f + 1 2852 # - + - = --------- = ------- 2853 # b b*f b*b*f b*f 2854 2855 # The first terms in the series are: 2856 2857 # 1 1 1 1 1 1 1 1 13700 2858 # -- + -- + -- + -- + -- + --- + --- + ---- = ----- 2859 # 1 1 2 6 24 120 720 5040 5040 2860 2861 # Note that we cannot simply reduce 13700/5040 to 685/252, but must keep 2862 # the numerator and the denominator! 2863 2864 if ($scale <= 75) { 2865 # set $x directly from a cached string form 2866 $x->{_m} = $LIB->_new("2718281828459045235360287471352662497757" . 2867 "2470936999595749669676277240766303535476"); 2868 $x->{sign} = '+'; 2869 $x->{_es} = '-'; 2870 $x->{_e} = $LIB->_new(79); 2871 } else { 2872 # compute A and B so that e = A / B. 2873 2874 # After some terms we end up with this, so we use it as a starting 2875 # point: 2876 my $A = $LIB->_new("9093339520860578540197197" . 2877 "0164779391644753259799242"); 2878 my $F = $LIB->_new(42); 2879 my $step = 42; 2880 2881 # Compute number of steps needed to get $A and $B sufficiently large. 2882 2883 my $steps = _len_to_steps($scale - 4); 2884 # print STDERR "# Doing $steps steps for ", $scale-4, " digits\n"; 2885 2886 while ($step++ <= $steps) { 2887 # calculate $a * $f + 1 2888 $A = $LIB -> _mul($A, $F); 2889 $A = $LIB -> _inc($A); 2890 # increment f 2891 $F = $LIB -> _inc($F); 2892 } 2893 2894 # Compute $B as factorial of $steps (this is faster than doing it 2895 # manually) 2896 my $B = $LIB->_fac($LIB->_new($steps)); 2897 2898 # print "A ", $LIB->_str($A), "\nB ", $LIB->_str($B), "\n"; 2899 2900 # compute A/B with $scale digits in the result (truncate, not round) 2901 $A = $LIB->_lsft($A, $LIB->_new($scale), 10); 2902 $A = $LIB->_div($A, $B); 2903 2904 $x->{_m} = $A; 2905 $x->{sign} = '+'; 2906 $x->{_es} = '-'; 2907 $x->{_e} = $LIB->_new($scale); 2908 } 2909 2910 # Now $x contains now an estimate of e, with some additional digits. 2911 2912 if ($x_orig -> is_one()) { 2913 2914 # else just round the already computed result 2915 2916 $x->{accuracy} = undef; 2917 $x->{precision} = undef; 2918 2919 # shortcut to not run through _find_round_parameters again 2920 2921 if (defined $params[0]) { 2922 $x = $x -> bround($params[0], $params[2]); # then round accordingly 2923 } else { 2924 $x = $x -> bfround($params[1], $params[2]); # then round accordingly 2925 } 2926 2927 } else { 2928 2929 # Use the fact exp(x) = exp(x/n)**n. In our case, n = 2**i for some 2930 # integer i. We use this to compute exp(y) where y = x / (2**i) and 2931 # 1 <= |y| < 2. 2932 # 2933 # The code below is similar to the code found in to_ieee754(). 2934 2935 # We need to find the base 2 exponent. First make an estimate of the 2936 # base 2 exponent, before adjusting it below. We could skip this 2937 # estimation and go straight to the while-loops below, but the loops 2938 # are slow, especially when the final exponent is far from zero and 2939 # even more so if the number of digits is large. This initial 2940 # estimation speeds up the computation dramatically. 2941 # 2942 # log2($m * 10**$e) = log10($m + 10**$e) * log(10)/log(2) 2943 # = (log10($m) + $e) * log(10)/log(2) 2944 # = (log($m)/log(10) + $e) * log(10)/log(2) 2945 2946 my ($m, $e) = $x_orig -> nparts(); 2947 my $ms = $m -> numify(); 2948 my $es = $e -> numify(); 2949 2950 # We start off by initializing the exponent to zero and the mantissa to 2951 # the input value. Then we increase the mantissa and decrease the 2952 # exponent, or vice versa, until the mantissa is in the desired range 2953 # or we hit one of the limits for the exponent. 2954 2955 my $mant = $x_orig -> copy() -> babs(); 2956 my $expo; 2957 2958 my $one = $class -> bone(); 2959 my $two = $class -> new("2"); 2960 my $half = $class -> new("0.5"); 2961 2962 my $expo_est = (log(abs($ms))/log(10) + $es) * log(10)/log(2); 2963 $expo_est = int($expo_est); 2964 2965 # Don't multiply by a number raised to a negative exponent. This will 2966 # cause a division, whose result is truncated to some fixed number of 2967 # digits. Instead, multiply by the inverse number raised to a positive 2968 # exponent. 2969 2970 $expo = $class -> new($expo_est); 2971 if ($expo_est > 0) { 2972 $mant = $mant -> bmul($half -> copy() -> bpow($expo)); 2973 } elsif ($expo_est < 0) { 2974 my $expo_abs = $expo -> copy() -> bneg(); 2975 $mant = $mant -> bmul($two -> copy() -> bpow($expo_abs)); 2976 } 2977 2978 # Final adjustment of the estimate above. 2979 2980 while ($mant -> bcmp($two) >= 0) { # $mant <= $two 2981 $mant = $mant -> bmul($half); 2982 $expo = $expo -> binc(); 2983 } 2984 2985 while ($mant -> bcmp($one) < 0) { # $mant > $one 2986 $mant = $mant -> bmul($two); 2987 $expo = $expo -> bdec(); 2988 } 2989 2990 # Because of the upscaling, we need some additional digits. 2991 2992 my $rescale = int($scale + abs($expo) * log(2) / log(10) + 1); 2993 $rescale = 4 if $rescale < 4; 2994 2995 $x = $x -> bpow($mant, $rescale); 2996 my $pow2 = $two -> bpow($expo, $rescale); 2997 $pow2 -> bneg() if $x_orig -> is_negative(); 2998 2999 # The bpow() below fails with the GMP and GMPz libraries if abs($pow2) 3000 # >= 2**30 = 1073741824. With the Pari library, it fails already when 3001 # abs($pow) >= 2**13 = 8192. With the Calc library, it is rediculously 3002 # slow when abs($pow2) is large. Fixme? 3003 3004 croak "cannot compute bexp(); input value is too large" 3005 if $pow2 -> copy() -> babs() -> bcmp("1073741824") >= 0; 3006 3007 $x = $x -> bpow($pow2, $rescale); 3008 3009 # Rounding parameters given as arguments currently don't override 3010 # instance variables, so accuracy (which is set in the computations 3011 # above) must be undefined before rounding. Fixme. 3012 3013 $x->{accuracy} = undef; 3014 $x -> round(@params); 3015 } 3016 3017 if ($fallback) { 3018 # clear a/p after round, since user did not request it 3019 $x->{accuracy} = undef; 3020 $x->{precision} = undef; 3021 } 3022 3023 # Restore globals. We need to do it like this, because setting one 3024 # undefines the other. 3025 3026 if (defined $ab) { 3027 $class -> accuracy($ab); 3028 } else { 3029 $class -> precision($pb); 3030 } 3031 3032 $class -> upgrade($upg); 3033 $class -> downgrade($dng); 3034 3035 # If downgrading, remember to preserve the relevant instance parameters. 3036 # There should be a more elegant way to do this. Fixme. 3037 3038 if ($downgrade && $x -> is_int()) { 3039 @r = ($x->{accuracy}, $x->{_r}); 3040 my $tmp = $downgrade -> new($x, @r); 3041 %$x = %$tmp; 3042 return bless $x, $downgrade; 3043 } 3044 3045 $x; 3046} 3047 3048sub bilog2 { 3049 croak "Method ", (caller(0))[3], "() not implemented yet"; 3050} 3051 3052sub bilog10 { 3053 croak "Method ", (caller(0))[3], "() not implemented yet"; 3054} 3055 3056sub bclog2 { 3057 croak "Method ", (caller(0))[3], "() not implemented yet"; 3058} 3059 3060sub bclog10 { 3061 croak "Method ", (caller(0))[3], "() not implemented yet"; 3062} 3063 3064sub bnok { 3065 # Calculate n over k (binomial coefficient or "choose" function) as integer. 3066 # set up parameters 3067 my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 3068 ? (ref($_[0]), @_) 3069 : objectify(2, @_); 3070 3071 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 3072 3073 return $x if $x->modify('bnok'); 3074 3075 return $x->bnan() if $x->is_nan() || $y->is_nan(); 3076 return $x->bnan() if (($x->is_finite() && !$x->is_int()) || 3077 ($y->is_finite() && !$y->is_int())); 3078 3079 my $xint = Math::BigInt -> new($x -> bsstr()); 3080 my $yint = Math::BigInt -> new($y -> bsstr()); 3081 $xint = $xint -> bnok($yint); 3082 3083 return $xint if defined $downgrade; 3084 3085 my $xflt = Math::BigFloat -> new($xint); 3086 3087 $x->{_m} = $xflt->{_m}; 3088 $x->{_e} = $xflt->{_e}; 3089 $x->{_es} = $xflt->{_es}; 3090 $x->{sign} = $xflt->{sign}; 3091 3092 return $x; 3093} 3094 3095sub bsin { 3096 # Calculate a sinus of x. 3097 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 3098 3099 # First we apply range reduction to x. This is because if x is large, the 3100 # Taylor series converges slowly and requires higher accuracy in the 3101 # intermediate computation. The Taylor series is: 3102 # 3103 # x^3 x^5 x^7 x^9 3104 # sin(x) = x - --- + --- - --- + --- ... 3105 # 3! 5! 7! 9! 3106 3107 return $x if $x -> modify('bsin'); 3108 3109 return $x -> bzero(@r) if $x -> is_zero(); 3110 return $x -> bnan(@r) if $x -> is_nan() || $x -> is_inf(); 3111 3112 # Get the rounding parameters, if any. 3113 3114 my $fallback = 0; 3115 my ($scale, @params); 3116 ($x, @params) = $x -> _find_round_parameters(@r); 3117 3118 # Error in _find_round_parameters? 3119 3120 return $x -> bnan(@r) if $x -> is_nan(); 3121 3122 # If no rounding parameters are given, use fallback. 3123 3124 if (!@params) { 3125 $params[0] = $class -> div_scale(); # fallback accuracy 3126 $params[1] = undef; # no precision 3127 $params[2] = $r[2]; # rounding mode 3128 $scale = $params[0]; 3129 $fallback = 1; # to clear a/p afterwards 3130 } else { 3131 if (defined($params[0])) { 3132 $scale = $params[0]; 3133 } else { 3134 # We perform the computations below using accuracy only, not 3135 # precision, so when precision is given, we need to "convert" this 3136 # to accuracy. 3137 $scale = 1 - $params[1]; 3138 } 3139 } 3140 3141 # Add more digits to the scale if the magnitude of $x is large. 3142 3143 my ($m, $e) = $x -> nparts(); 3144 $scale += $e if $x >= 10; 3145 $scale = 4 if $scale < 4; 3146 3147 # When user set globals, they would interfere with our calculation, so 3148 # disable them and later re-enable them 3149 3150 my $ab = $class -> accuracy(); 3151 my $pb = $class -> precision(); 3152 $class -> accuracy(undef); 3153 $class -> precision(undef); 3154 3155 # Disabling upgrading and downgrading is no longer necessary to avoid an 3156 # infinite recursion, but it avoids unnecessary upgrading and downgrading in 3157 # the intermediate computations. 3158 3159 my $upg = $class -> upgrade(); 3160 my $dng = $class -> downgrade(); 3161 $class -> upgrade(undef); 3162 $class -> downgrade(undef); 3163 3164 # We also need to disable any set A or P on $x (_find_round_parameters took 3165 # them already into account), since these would interfere, too. 3166 3167 $x->{accuracy} = undef; 3168 $x->{precision} = undef; 3169 3170 my $sin_prev; # the previous approximation of sin(x) 3171 my $sin; # the current approximation of sin(x) 3172 3173 while (1) { 3174 3175 # Compute constants to the current scale. 3176 3177 my $pi = $class -> bpi($scale); # �� 3178 my $twopi = $pi -> copy() -> bmul("2"); # 2�� 3179 my $halfpi = $pi -> copy() -> bmul("0.5"); # ��/2 3180 3181 # Use the fact that sin(-x) = -sin(x) to reduce the range to the 3182 # interval to [0,∞). 3183 3184 my $xsgn = $x < 0 ? -1 : 1; 3185 my $x = $x -> copy() -> babs(); 3186 3187 # Use the fact that sin(2��x) = sin(x) to reduce the range to the 3188 # interval to [0, 2��). 3189 3190 $x = $x -> bmod($twopi, $scale); 3191 3192 # Use the fact that sin(x+��) = -sin(x) to reduce the range to the 3193 # interval to [0,��). 3194 3195 if ($x -> bcmp($pi) > 0) { 3196 $xsgn = -$xsgn; 3197 $x = $x -> bsub($pi); 3198 } 3199 3200 # Use the fact that sin(��-x) = sin(x) to reduce the range to the 3201 # interval [0,��/2). 3202 3203 if ($x -> bcmp($halfpi) > 0) { 3204 $x = $x -> bsub($pi) -> bneg(); # �� - x 3205 } 3206 3207 my $tol = $class -> new("1E-". ($scale-1)); 3208 3209 my $xsq = $x -> copy() -> bmul($x, $scale) -> bneg(); 3210 my $term = $x -> copy(); 3211 my $fac = $class -> bone(); 3212 my $n = $class -> bone(); 3213 3214 $sin = $x -> copy(); # initialize sin(x) to the first term 3215 3216 while (1) { 3217 $n -> binc(); 3218 $fac = $n -> copy(); 3219 $n -> binc(); 3220 $fac -> bmul($n); 3221 3222 $term -> bmul($xsq, $scale) -> bdiv($fac, $scale); 3223 3224 $sin -> badd($term, $scale); 3225 last if $term -> copy() -> babs() -> bcmp($tol) < 0; 3226 } 3227 3228 $sin -> bneg() if $xsgn < 0; 3229 3230 # Rounding parameters given as arguments currently don't override 3231 # instance variables, so accuracy (which is set in the computations 3232 # above) must be undefined before rounding. Fixme. 3233 3234 $sin->{accuracy} = undef; 3235 $sin -> round(@params); 3236 3237 # Compare the current approximation of sin(x) with the previous one, 3238 # and if they are identical, we're done. 3239 3240 if (defined $sin_prev) { 3241 last if $sin -> bcmp($sin_prev) == 0; 3242 } 3243 3244 # If the current approximation of sin(x) is different from the previous 3245 # approximation, double the scale (accuracy) and retry. 3246 3247 $sin_prev = $sin; 3248 $scale *= 2; 3249 } 3250 3251 # Assign the result to the invocand. 3252 3253 %$x = %$sin; 3254 3255 if ($fallback) { 3256 # clear a/p after round, since user did not request it 3257 $x->{accuracy} = undef; 3258 $x->{precision} = undef; 3259 } 3260 3261 # Restore globals. We need to do it like this, because setting one 3262 # undefines the other. 3263 3264 if (defined $ab) { 3265 $class -> accuracy($ab); 3266 } else { 3267 $class -> precision($pb); 3268 } 3269 3270 $class -> upgrade($upg); 3271 $class -> downgrade($dng); 3272 3273 # If downgrading, remember to preserve the relevant instance parameters. 3274 # There should be a more elegant way to do this. Fixme. 3275 3276 if ($downgrade && $x -> is_int()) { 3277 @r = ($x->{accuracy}, $x->{_r}); 3278 my $tmp = $downgrade -> new($x, @r); 3279 %$x = %$tmp; 3280 return bless $x, $downgrade; 3281 } 3282 3283 $x; 3284} 3285 3286sub bcos { 3287 # Calculate a cosinus of x. 3288 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 3289 3290 # Taylor: x^2 x^4 x^6 x^8 3291 # cos = 1 - --- + --- - --- + --- ... 3292 # 2! 4! 6! 8! 3293 3294 # we need to limit the accuracy to protect against overflow 3295 my $fallback = 0; 3296 my ($scale, @params); 3297 ($x, @params) = $x->_find_round_parameters(@r); 3298 3299 # constant object or error in _find_round_parameters? 3300 return $x if $x->modify('bcos') || $x->is_nan(); 3301 return $x->bnan() if $x->is_inf(); 3302 return $x->bone(@r) if $x->is_zero(); 3303 3304 # no rounding at all, so must use fallback 3305 if (scalar @params == 0) { 3306 # simulate old behaviour 3307 $params[0] = $class->div_scale(); # and round to it as accuracy 3308 $params[1] = undef; # disable P 3309 $scale = $params[0]+4; # at least four more for proper round 3310 $params[2] = $r[2]; # round mode by caller or undef 3311 $fallback = 1; # to clear a/p afterwards 3312 } else { 3313 # the 4 below is empirical, and there might be cases where it is not 3314 # enough... 3315 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 3316 } 3317 3318 # When user set globals, they would interfere with our calculation, so 3319 # disable them and later re-enable them. 3320 3321 my $ab = $class -> accuracy(); 3322 my $pb = $class -> precision(); 3323 $class -> accuracy(undef); 3324 $class -> precision(undef); 3325 3326 # Disabling upgrading and downgrading is no longer necessary to avoid an 3327 # infinite recursion, but it avoids unnecessary upgrading and downgrading in 3328 # the intermediate computations. 3329 3330 my $upg = $class -> upgrade(); 3331 my $dng = $class -> downgrade(); 3332 $class -> upgrade(undef); 3333 $class -> downgrade(undef); 3334 3335 # We also need to disable any set A or P on $x (_find_round_parameters took 3336 # them already into account), since these would interfere, too. 3337 3338 $x->{accuracy} = undef; 3339 $x->{precision} = undef; 3340 3341 my $over = $x * $x; # X ^ 2 3342 my $x2 = $over->copy(); # X ^ 2; difference between terms 3343 my $sign = 1; # start with -= 3344 my $below = $class->new(2); 3345 my $factorial = $class->new(3); 3346 $x = $x->bone(); 3347 $x->{accuracy} = undef; 3348 $x->{precision} = undef; 3349 3350 my $limit = $class->new("1E-". ($scale-1)); 3351 #my $steps = 0; 3352 while (3 < 5) { 3353 # we calculate the next term, and add it to the last 3354 # when the next term is below our limit, it won't affect the outcome 3355 # anymore, so we stop: 3356 my $next = $over->copy()->bdiv($below, $scale); 3357 last if $next->bacmp($limit) <= 0; 3358 3359 if ($sign == 0) { 3360 $x = $x->badd($next); 3361 } else { 3362 $x = $x->bsub($next); 3363 } 3364 $sign = 1-$sign; # alternate 3365 # calculate things for the next term 3366 $over = $over->bmul($x2); # $x*$x 3367 $below = $below->bmul($factorial); # n*(n+1) 3368 $factorial = $factorial -> binc(); 3369 $below = $below->bmul($factorial); # n*(n+1) 3370 $factorial = $factorial -> binc(); 3371 } 3372 3373 # shortcut to not run through _find_round_parameters again 3374 if (defined $params[0]) { 3375 $x = $x->bround($params[0], $params[2]); # then round accordingly 3376 } else { 3377 $x = $x->bfround($params[1], $params[2]); # then round accordingly 3378 } 3379 if ($fallback) { 3380 # clear a/p after round, since user did not request it 3381 $x->{accuracy} = undef; 3382 $x->{precision} = undef; 3383 } 3384 3385 # Restore globals. We need to do it like this, because setting one 3386 # undefines the other. 3387 3388 if (defined $ab) { 3389 $class -> accuracy($ab); 3390 } else { 3391 $class -> precision($pb); 3392 } 3393 3394 $class -> upgrade($upg); 3395 $class -> downgrade($dng); 3396 3397 return $downgrade -> new($x -> bdstr(), @r) 3398 if defined($downgrade) && $x -> is_int(); 3399 $x; 3400} 3401 3402sub batan { 3403 # Calculate a arcus tangens of x. 3404 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 3405 3406 # taylor: x^3 x^5 x^7 x^9 3407 # atan = x - --- + --- - --- + --- ... 3408 # 3 5 7 9 3409 3410 return $x if $x->modify('batan'); 3411 3412 return $x -> bnan(@r) if $x->is_nan(); 3413 3414 # We need to limit the accuracy to protect against overflow. 3415 3416 my $fallback = 0; 3417 my ($scale, @params); 3418 ($x, @params) = $x->_find_round_parameters(@r); 3419 3420 # Error in _find_round_parameters? 3421 3422 return $x -> bnan(@r) if $x->is_nan(); 3423 3424 if ($x->{sign} =~ /^[+-]inf\z/) { 3425 # +inf result is PI/2 3426 # -inf result is -PI/2 3427 # calculate PI/2 3428 my $pi = $class->bpi(@r); 3429 # modify $x in place 3430 $x->{_m} = $pi->{_m}; 3431 $x->{_e} = $pi->{_e}; 3432 $x->{_es} = $pi->{_es}; 3433 # -y => -PI/2, +y => PI/2 3434 $x->{sign} = substr($x->{sign}, 0, 1); # "+inf" => "+" 3435 $x -> {_m} = $LIB->_div($x->{_m}, $LIB->_new(2)); 3436 return $x; 3437 } 3438 3439 return $x->bzero(@r) if $x->is_zero(); 3440 3441 # no rounding at all, so must use fallback 3442 if (scalar @params == 0) { 3443 # simulate old behaviour 3444 $params[0] = $class->div_scale(); # and round to it as accuracy 3445 $params[1] = undef; # disable P 3446 $scale = $params[0]+4; # at least four more for proper round 3447 $params[2] = $r[2]; # round mode by caller or undef 3448 $fallback = 1; # to clear a/p afterwards 3449 } else { 3450 # the 4 below is empirical, and there might be cases where it is not 3451 # enough... 3452 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 3453 } 3454 3455 # 1 or -1 => PI/4 3456 # inlined is_one() && is_one('-') 3457 if ($LIB->_is_one($x->{_m}) && $LIB->_is_zero($x->{_e})) { 3458 my $pi = $class->bpi($scale - 3); 3459 # modify $x in place 3460 $x->{_m} = $pi->{_m}; 3461 $x->{_e} = $pi->{_e}; 3462 $x->{_es} = $pi->{_es}; 3463 # leave the sign of $x alone (+1 => +PI/4, -1 => -PI/4) 3464 $x->{_m} = $LIB->_div($x->{_m}, $LIB->_new(4)); 3465 return $x; 3466 } 3467 3468 # This series is only valid if -1 < x < 1, so for other x we need to 3469 # calculate PI/2 - atan(1/x): 3470 my $pi = undef; 3471 if ($x->bacmp($x->copy()->bone) >= 0) { 3472 # calculate PI/2 3473 $pi = $class->bpi($scale - 3); 3474 $pi->{_m} = $LIB->_div($pi->{_m}, $LIB->_new(2)); 3475 # calculate 1/$x: 3476 my $x_copy = $x->copy(); 3477 # modify $x in place 3478 $x = $x->bone(); 3479 $x = $x->bdiv($x_copy, $scale); 3480 } 3481 3482 my $fmul = 1; 3483 foreach (0 .. int($scale / 20)) { 3484 $fmul *= 2; 3485 $x = $x->bdiv($x->copy()->bmul($x)->binc()->bsqrt($scale + 4)->binc(), 3486 $scale + 4); 3487 } 3488 3489 # When user set globals, they would interfere with our calculation, so 3490 # disable them and later re-enable them. 3491 3492 my $ab = $class -> accuracy(); 3493 my $pb = $class -> precision(); 3494 $class -> accuracy(undef); 3495 $class -> precision(undef); 3496 3497 # Disabling upgrading and downgrading is no longer necessary to avoid an 3498 # infinite recursion, but it avoids unnecessary upgrading and downgrading in 3499 # the intermediate computations. 3500 3501 my $upg = $class -> upgrade(); 3502 my $dng = $class -> downgrade(); 3503 $class -> upgrade(undef); 3504 $class -> downgrade(undef); 3505 3506 # We also need to disable any set A or P on $x (_find_round_parameters took 3507 # them already into account), since these would interfere, too. 3508 3509 $x->{accuracy} = undef; 3510 $x->{precision} = undef; 3511 3512 my $over = $x * $x; # X ^ 2 3513 my $x2 = $over->copy(); # X ^ 2; difference between terms 3514 $over = $over->bmul($x); # X ^ 3 as starting value 3515 my $sign = 1; # start with -= 3516 my $below = $class->new(3); 3517 my $two = $class->new(2); 3518 $x->{accuracy} = undef; 3519 $x->{precision} = undef; 3520 3521 my $limit = $class->new("1E-". ($scale-1)); 3522 #my $steps = 0; 3523 while (1) { 3524 # We calculate the next term, and add it to the last. When the next 3525 # term is below our limit, it won't affect the outcome anymore, so we 3526 # stop: 3527 my $next = $over->copy()->bdiv($below, $scale); 3528 last if $next->bacmp($limit) <= 0; 3529 3530 if ($sign == 0) { 3531 $x = $x->badd($next); 3532 } else { 3533 $x = $x->bsub($next); 3534 } 3535 $sign = 1-$sign; # alternatex 3536 # calculate things for the next term 3537 $over = $over->bmul($x2); # $x*$x 3538 $below = $below->badd($two); # n += 2 3539 } 3540 $x = $x->bmul($fmul); 3541 3542 if (defined $pi) { 3543 my $x_copy = $x->copy(); 3544 # modify $x in place 3545 $x->{_m} = $pi->{_m}; 3546 $x->{_e} = $pi->{_e}; 3547 $x->{_es} = $pi->{_es}; 3548 # PI/2 - $x 3549 $x = $x->bsub($x_copy); 3550 } 3551 3552 # Shortcut to not run through _find_round_parameters again. 3553 if (defined $params[0]) { 3554 $x = $x->bround($params[0], $params[2]); # then round accordingly 3555 } else { 3556 $x = $x->bfround($params[1], $params[2]); # then round accordingly 3557 } 3558 if ($fallback) { 3559 # Clear a/p after round, since user did not request it. 3560 $x->{accuracy} = undef; 3561 $x->{precision} = undef; 3562 } 3563 3564 # Restore globals. We need to do it like this, because setting one 3565 # undefines the other. 3566 3567 if (defined $ab) { 3568 $class -> accuracy($ab); 3569 } else { 3570 $class -> precision($pb); 3571 } 3572 3573 $class -> upgrade($upg); 3574 $class -> downgrade($dng); 3575 3576 return $downgrade -> new($x -> bdstr(), @r) 3577 if defined($downgrade) && ($x -> is_int() || $x -> is_inf()); 3578 $x; 3579} 3580 3581sub batan2 { 3582 # $y -> batan2($x) returns the arcus tangens of $y / $x. 3583 3584 # Set up parameters. 3585 my ($class, $y, $x, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 3586 ? (ref($_[0]), @_) 3587 : objectify(2, @_); 3588 3589 # Quick exit if $y is read-only. 3590 return $y if $y -> modify('batan2'); 3591 3592 # Handle all NaN cases. 3593 return $y -> bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan; 3594 3595 # We need to limit the accuracy to protect against overflow. 3596 my $fallback = 0; 3597 my ($scale, @params); 3598 ($y, @params) = $y -> _find_round_parameters(@r); 3599 3600 # Error in _find_round_parameters? 3601 return $y if $y->is_nan(); 3602 3603 # No rounding at all, so must use fallback. 3604 if (scalar @params == 0) { 3605 # Simulate old behaviour 3606 $params[0] = $class -> div_scale(); # and round to it as accuracy 3607 $params[1] = undef; # disable P 3608 $scale = $params[0] + 4; # at least four more for proper round 3609 $params[2] = $r[2]; # round mode by caller or undef 3610 $fallback = 1; # to clear a/p afterwards 3611 } else { 3612 # The 4 below is empirical, and there might be cases where it is not 3613 # enough ... 3614 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 3615 } 3616 3617 if ($x -> is_inf("+")) { # x = inf 3618 if ($y -> is_inf("+")) { # y = inf 3619 $y = $y -> bpi($scale) -> bmul("0.25"); # pi/4 3620 } elsif ($y -> is_inf("-")) { # y = -inf 3621 $y = $y -> bpi($scale) -> bmul("-0.25"); # -pi/4 3622 } else { # -inf < y < inf 3623 return $y -> bzero(@r); # 0 3624 } 3625 } elsif ($x -> is_inf("-")) { # x = -inf 3626 if ($y -> is_inf("+")) { # y = inf 3627 $y = $y -> bpi($scale) -> bmul("0.75"); # 3/4 pi 3628 } elsif ($y -> is_inf("-")) { # y = -inf 3629 $y = $y -> bpi($scale) -> bmul("-0.75"); # -3/4 pi 3630 } elsif ($y >= 0) { # y >= 0 3631 $y = $y -> bpi($scale); # pi 3632 } else { # y < 0 3633 $y = $y -> bpi($scale) -> bneg(); # -pi 3634 } 3635 } elsif ($x > 0) { # 0 < x < inf 3636 if ($y -> is_inf("+")) { # y = inf 3637 $y = $y -> bpi($scale) -> bmul("0.5"); # pi/2 3638 } elsif ($y -> is_inf("-")) { # y = -inf 3639 $y = $y -> bpi($scale) -> bmul("-0.5"); # -pi/2 3640 } else { # -inf < y < inf 3641 $y = $y -> bdiv($x, $scale) -> batan($scale); # atan(y/x) 3642 } 3643 } elsif ($x < 0) { # -inf < x < 0 3644 my $pi = $class -> bpi($scale); 3645 if ($y >= 0) { # y >= 0 3646 $y = $y -> bdiv($x, $scale) -> batan() # atan(y/x) + pi 3647 -> badd($pi); 3648 } else { # y < 0 3649 $y = $y -> bdiv($x, $scale) -> batan() # atan(y/x) - pi 3650 -> bsub($pi); 3651 } 3652 } else { # x = 0 3653 if ($y > 0) { # y > 0 3654 $y = $y -> bpi($scale) -> bmul("0.5"); # pi/2 3655 } elsif ($y < 0) { # y < 0 3656 $y = $y -> bpi($scale) -> bmul("-0.5"); # -pi/2 3657 } else { # y = 0 3658 return $y -> bzero(@r); # 0 3659 } 3660 } 3661 3662 $y = $y -> round(@r); 3663 3664 if ($fallback) { 3665 $y->{accuracy} = undef; 3666 $y->{precision} = undef; 3667 } 3668 3669 return $y; 3670} 3671 3672sub bsqrt { 3673 # calculate square root 3674 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 3675 3676 return $x if $x->modify('bsqrt'); 3677 3678 # Handle trivial cases. 3679 3680 return $x -> bnan(@r) if $x->is_nan(); 3681 return $x -> binf("+", @r) if $x->{sign} eq '+inf'; 3682 return $x -> round(@r) if $x->is_zero() || $x->is_one(); 3683 3684 # We don't support complex numbers. 3685 3686 if ($x -> is_neg()) { 3687 return $upgrade -> bsqrt($x, @r) if defined($upgrade); 3688 return $x -> bnan(@r); 3689 } 3690 3691 # we need to limit the accuracy to protect against overflow 3692 my $fallback = 0; 3693 my (@params, $scale); 3694 ($x, @params) = $x->_find_round_parameters(@r); 3695 3696 # error in _find_round_parameters? 3697 return $x -> bnan(@r) if $x->is_nan(); 3698 3699 # no rounding at all, so must use fallback 3700 if (scalar @params == 0) { 3701 # simulate old behaviour 3702 $params[0] = $class->div_scale(); # and round to it as accuracy 3703 $scale = $params[0]+4; # at least four more for proper round 3704 $params[2] = $r[2]; # round mode by caller or undef 3705 $fallback = 1; # to clear a/p afterwards 3706 } else { 3707 # the 4 below is empirical, and there might be cases where it is not 3708 # enough... 3709 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 3710 } 3711 3712 # Shift the significand left or right to get the desired number of digits, 3713 # which is 2*$scale with possibly one extra digit to ensure that the 3714 # exponent is an even number. 3715 3716 my $l = $LIB -> _len($x->{_m}); 3717 my $n = 2 * $scale - $l; # how much should we shift? 3718 $n++ if ($l % 2 xor $LIB -> _is_odd($x->{_e})); 3719 my ($na, $ns) = $n < 0 ? (abs($n), "-") : ($n, "+"); 3720 $na = $LIB -> _new($na); 3721 3722 $x->{_m} = $ns eq "+" ? $LIB -> _lsft($x->{_m}, $na, 10) 3723 : $LIB -> _rsft($x->{_m}, $na, 10); 3724 3725 $x->{_m} = $LIB -> _sqrt($x->{_m}); 3726 3727 # Adjust the exponent by the amount that we shifted the significand. The 3728 # square root of the exponent is simply half of it: sqrt(10^(2*a)) = 10^a. 3729 3730 ($x->{_e}, $x->{_es}) = $LIB -> _ssub($x->{_e}, $x->{_es}, $na, $ns); 3731 $x->{_e} = $LIB -> _div($x->{_e}, $LIB -> _new("2")); 3732 3733 # Normalize to get rid of any trailing zeros in the significand. 3734 3735 $x -> bnorm(); 3736 3737 # shortcut to not run through _find_round_parameters again 3738 if (defined $params[0]) { 3739 $x = $x->bround($params[0], $params[2]); # then round accordingly 3740 } else { 3741 $x = $x->bfround($params[1], $params[2]); # then round accordingly 3742 } 3743 3744 if ($fallback) { 3745 # clear a/p after round, since user did not request it 3746 $x->{accuracy} = undef; 3747 $x->{precision} = undef; 3748 } 3749 3750 return $downgrade -> new($x, @r) 3751 if defined($downgrade) && $x -> is_int(); 3752 $x; 3753} 3754 3755sub broot { 3756 # calculate $y'th root of $x 3757 3758 # set up parameters 3759 my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 3760 ? (ref($_[0]), @_) 3761 : objectify(2, @_); 3762 3763 return $x if $x->modify('broot'); 3764 3765 # Handle trivial cases. 3766 3767 return $x -> bnan(@r) if $x->is_nan() || $y->is_nan(); 3768 3769 if ($x -> is_neg()) { 3770 # -27 ** (1/3) = -3 3771 return $x -> broot($y -> copy() -> bneg(), @r) -> bneg() 3772 if $x -> is_int() && $y -> is_int() && $y -> is_neg(); 3773 return $upgrade -> broot($x, $y, @r) if defined $upgrade; 3774 return $x -> bnan(@r); 3775 } 3776 3777 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0 3778 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() || 3779 $y->{sign} !~ /^\+$/; 3780 3781 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one(); 3782 3783 # we need to limit the accuracy to protect against overflow 3784 my $fallback = 0; 3785 my (@params, $scale); 3786 ($x, @params) = $x->_find_round_parameters(@r); 3787 3788 return $x if $x->is_nan(); # error in _find_round_parameters? 3789 3790 # no rounding at all, so must use fallback 3791 if (scalar @params == 0) { 3792 # simulate old behaviour 3793 $params[0] = $class->div_scale(); # and round to it as accuracy 3794 $scale = $params[0]+4; # at least four more for proper round 3795 $params[2] = $r[2]; # round mode by caller or undef 3796 $fallback = 1; # to clear a/p afterwards 3797 } else { 3798 # the 4 below is empirical, and there might be cases where it is not 3799 # enough... 3800 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 3801 } 3802 3803 # When user set globals, they would interfere with our calculation, so 3804 # disable them and later re-enable them. 3805 3806 my $ab = $class -> accuracy(); 3807 my $pb = $class -> precision(); 3808 $class -> accuracy(undef); 3809 $class -> precision(undef); 3810 3811 # Disabling upgrading and downgrading is no longer necessary to avoid an 3812 # infinite recursion, but it avoids unnecessary upgrading and downgrading in 3813 # the intermediate computations. 3814 3815 my $upg = $class -> upgrade(); 3816 my $dng = $class -> downgrade(); 3817 $class -> upgrade(undef); 3818 $class -> downgrade(undef); 3819 3820 # We also need to disable any set A or P on $x (_find_round_parameters took 3821 # them already into account), since these would interfere, too. 3822 3823 $x->{accuracy} = undef; 3824 $x->{precision} = undef; 3825 3826 # remember sign and make $x positive, since -4 ** (1/2) => -2 3827 my $sign = 0; 3828 $sign = 1 if $x->{sign} eq '-'; 3829 $x->{sign} = '+'; 3830 3831 my $is_two = 0; 3832 if ($y->isa('Math::BigFloat')) { 3833 $is_two = $y->{sign} eq '+' && $LIB->_is_two($y->{_m}) 3834 && $LIB->_is_zero($y->{_e}); 3835 } else { 3836 $is_two = $y == 2; 3837 } 3838 3839 # normal square root if $y == 2: 3840 if ($is_two) { 3841 $x = $x->bsqrt($scale+4); 3842 } elsif ($y->is_one('-')) { 3843 # $x ** -1 => 1/$x 3844 my $u = $class->bone()->bdiv($x, $scale); 3845 # copy private parts over 3846 $x->{_m} = $u->{_m}; 3847 $x->{_e} = $u->{_e}; 3848 $x->{_es} = $u->{_es}; 3849 } else { 3850 # calculate the broot() as integer result first, and if it fits, return 3851 # it rightaway (but only if $x and $y are integer): 3852 3853 my $done = 0; # not yet 3854 if ($y->is_int() && $x->is_int()) { 3855 my $i = $LIB->_copy($x->{_m}); 3856 $i = $LIB->_lsft($i, $x->{_e}, 10) unless $LIB->_is_zero($x->{_e}); 3857 my $int = Math::BigInt->bzero(); 3858 $int->{value} = $i; 3859 $int = $int->broot($y->as_number()); 3860 # if ($exact) 3861 if ($int->copy()->bpow($y) == $x) { 3862 # found result, return it 3863 $x->{_m} = $int->{value}; 3864 $x->{_e} = $LIB->_zero(); 3865 $x->{_es} = '+'; 3866 $x = $x->bnorm(); 3867 $done = 1; 3868 } 3869 } 3870 if ($done == 0) { 3871 my $u = $class->bone()->bdiv($y, $scale+4); 3872 $u->{accuracy} = undef; 3873 $u->{precision} = undef; 3874 $x = $x->bpow($u, $scale+4); # el cheapo 3875 } 3876 } 3877 $x = $x->bneg() if $sign == 1; 3878 3879 # shortcut to not run through _find_round_parameters again 3880 if (defined $params[0]) { 3881 $x = $x->bround($params[0], $params[2]); # then round accordingly 3882 } else { 3883 $x = $x->bfround($params[1], $params[2]); # then round accordingly 3884 } 3885 if ($fallback) { 3886 # clear a/p after round, since user did not request it 3887 $x->{accuracy} = undef; 3888 $x->{precision} = undef; 3889 } 3890 3891 # Restore globals. We need to do it like this, because setting one 3892 # undefines the other. 3893 3894 if (defined $ab) { 3895 $class -> accuracy($ab); 3896 } else { 3897 $class -> precision($pb); 3898 } 3899 3900 $class -> upgrade($upg); 3901 $class -> downgrade($dng); 3902 3903 return $downgrade -> new($x -> bdstr(), @r) 3904 if defined($downgrade) && ($x -> is_int() || $x -> is_inf()); 3905 $x; 3906} 3907 3908sub bfac { 3909 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT 3910 # compute factorial number, modifies first argument 3911 3912 # set up parameters 3913 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 3914 3915 # inf => inf 3916 return $x if $x->modify('bfac'); 3917 3918 return $x -> bnan(@r) if $x->is_nan() || $x->is_inf("-"); 3919 return $x -> binf("+", @r) if $x->is_inf("+"); 3920 return $x -> bone(@r) if $x->is_zero() || $x->is_one(); 3921 3922 if ($x -> is_neg() || !$x -> is_int()) { 3923 return $upgrade -> bfac($x, @r) if defined($upgrade); 3924 return $x -> bnan(@r); 3925 } 3926 3927 if (! $LIB->_is_zero($x->{_e})) { 3928 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # change 12e1 to 120e0 3929 $x->{_e} = $LIB->_zero(); # normalize 3930 $x->{_es} = '+'; 3931 } 3932 $x->{_m} = $LIB->_fac($x->{_m}); # calculate factorial 3933 3934 $x = $x->bnorm()->round(@r); # norm again and round result 3935 3936 return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade) 3937 && ($x -> is_int() || $x -> is_inf()); 3938 $x; 3939} 3940 3941sub bdfac { 3942 # compute double factorial 3943 3944 # set up parameters 3945 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 3946 3947 return $x if $x->modify('bdfac'); 3948 3949 return $x -> bnan(@r) if $x->is_nan() || $x->is_inf("-"); 3950 return $x -> binf("+", @r) if $x->is_inf("+"); 3951 3952 if ($x <= -2 || !$x -> is_int()) { 3953 return $upgrade -> bdfac($x, @r) if defined($upgrade); 3954 return $x -> bnan(@r); 3955 } 3956 3957 return $x->bone() if $x <= 1; 3958 3959 croak("bdfac() requires a newer version of the $LIB library.") 3960 unless $LIB->can('_dfac'); 3961 3962 if (! $LIB->_is_zero($x->{_e})) { 3963 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # change 12e1 to 120e0 3964 $x->{_e} = $LIB->_zero(); # normalize 3965 $x->{_es} = '+'; 3966 } 3967 $x->{_m} = $LIB->_dfac($x->{_m}); # calculate factorial 3968 3969 $x = $x->bnorm()->round(@r); # norm again and round result 3970 3971 return $downgrade -> new($x -> bdstr(), @r) 3972 if defined($downgrade) && $x -> is_int(); 3973 return $x; 3974} 3975 3976sub btfac { 3977 # compute triple factorial 3978 3979 # set up parameters 3980 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 3981 3982 return $x if $x->modify('btfac'); 3983 3984 return $x -> bnan(@r) if $x->is_nan() || $x->is_inf("-"); 3985 return $x -> binf("+", @r) if $x->is_inf("+"); 3986 3987 if ($x <= -3 || !$x -> is_int()) { 3988 return $upgrade -> btfac($x, @r) if defined($upgrade); 3989 return $x -> bnan(@r); 3990 } 3991 3992 my $k = $class -> new("3"); 3993 return $x->bnan(@r) if $x <= -$k; 3994 3995 my $one = $class -> bone(); 3996 return $x->bone(@r) if $x <= $one; 3997 3998 my $f = $x -> copy(); 3999 while ($f -> bsub($k) > $one) { 4000 $x = $x -> bmul($f); 4001 } 4002 4003 $x = $x->round(@r); 4004 4005 return $downgrade -> new($x -> bdstr(), @r) 4006 if defined($downgrade) && $x -> is_int(); 4007 return $x; 4008} 4009 4010sub bmfac { 4011 my ($class, $x, $k, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 4012 ? (ref($_[0]), @_) 4013 : objectify(2, @_); 4014 4015 return $x if $x->modify('bmfac'); 4016 4017 return $x -> bnan(@r) if $x->is_nan() || $x->is_inf("-") || !$k->is_pos(); 4018 return $x -> binf("+", @r) if $x->is_inf("+"); 4019 4020 if ($x <= -$k || !$x -> is_int() || 4021 ($k -> is_finite() && !$k -> is_int())) 4022 { 4023 return $upgrade -> bmfac($x, $k, @r) if defined($upgrade); 4024 return $x -> bnan(@r); 4025 } 4026 4027 my $one = $class -> bone(); 4028 return $x->bone(@r) if $x <= $one; 4029 4030 my $f = $x -> copy(); 4031 while ($f -> bsub($k) > $one) { 4032 $x = $x -> bmul($f); 4033 } 4034 4035 $x = $x->round(@r); 4036 4037 return $downgrade -> new($x -> bdstr(), @r) 4038 if defined($downgrade) && $x -> is_int(); 4039 return $x; 4040} 4041 4042sub blsft { 4043 # shift left by $y in base $b, i.e., multiply by $b ** $y 4044 4045 # set up parameters 4046 my ($class, $x, $y, $b, @r) 4047 = ref($_[0]) && ref($_[0]) eq ref($_[1]) && ref($_[1]) eq ref($_[2]) 4048 ? (ref($_[0]), @_) 4049 : objectify(2, @_); 4050 4051 return $x if $x -> modify('blsft'); 4052 4053 return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan(); 4054 4055 $b = 2 if !defined $b; 4056 $b = $class -> new($b) 4057 unless defined(blessed($b)) && $b -> isa(__PACKAGE__); 4058 return $x -> bnan(@r) if $b -> is_nan(); 4059 4060 # There needs to be more checking for special cases here. Fixme! 4061 4062 # shift by a negative amount? 4063 return $x -> brsft($y -> copy() -> babs(), $b) if $y -> {sign} =~ /^-/; 4064 4065 $x = $x -> bmul($b -> bpow($y), $r[0], $r[1], $r[2], $y); 4066 4067 return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade) 4068 && ($x -> is_int() || $x -> is_inf() || $x -> is_nan()); 4069 return $x; 4070} 4071 4072sub brsft { 4073 # shift right by $y in base $b, i.e., divide by $b ** $y 4074 4075 # set up parameters 4076 my ($class, $x, $y, $b, @r) 4077 = ref($_[0]) && ref($_[0]) eq ref($_[1]) && ref($_[1]) eq ref($_[2]) 4078 ? (ref($_[0]), @_) 4079 : objectify(2, @_); 4080 4081 return $x if $x -> modify('brsft'); 4082 4083 return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan(); 4084 4085 # There needs to be more checking for special cases here. Fixme! 4086 4087 $b = 2 if !defined $b; 4088 $b = $class -> new($b) 4089 unless defined(blessed($b)) && $b -> isa(__PACKAGE__); 4090 return $x -> bnan(@r) if $b -> is_nan(); 4091 4092 # shift by a negative amount? 4093 return $x -> blsft($y -> copy() -> babs(), $b) if $y -> {sign} =~ /^-/; 4094 4095 # call bdiv() 4096 $x = $x -> bdiv($b -> bpow($y), $r[0], $r[1], $r[2], $y); 4097 4098 return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade) 4099 && ($x -> is_int() || $x -> is_inf() || $x -> is_nan()); 4100 return $x; 4101} 4102 4103############################################################################### 4104# Bitwise methods 4105############################################################################### 4106 4107# Bitwise left shift. 4108 4109sub bblsft { 4110 # We don't call objectify(), because the bitwise methods should not 4111 # upgrade/downgrade, even when upgrading/downgrading is enabled. 4112 4113 my ($class, $x, $y, @r) = ref($_[0]) ? (ref($_[0]), @_) : @_; 4114 4115 my $xint = Math::BigInt -> bblsft($x, $y, @r); 4116 4117 # Temporarily disable downgrading. 4118 4119 my $dng = $class -> downgrade(); 4120 $class -> downgrade(undef); 4121 4122 # convert to our class without downgrading. 4123 4124 my $xflt = $class -> new($xint); 4125 4126 # Reset downgrading. 4127 4128 $class -> downgrade($dng); 4129 4130 # If we are called as a class method, the first operand might not be an 4131 # object of this class, so check. 4132 4133 if (defined(blessed($x)) && $x -> isa(__PACKAGE__)) { 4134 $x -> {sign} = $xflt -> {sign}; 4135 $x -> {_m} = $xflt -> {_m}; 4136 $x -> {_es} = $xflt -> {_es}; 4137 $x -> {_e} = $xflt -> {_e}; 4138 } else { 4139 $x = $xflt; 4140 } 4141 4142 # Now we might downgrade. 4143 4144 return $downgrade -> new($x) if defined($downgrade); 4145 $x -> round(@r); 4146} 4147 4148# Bitwise right shift. 4149 4150sub bbrsft { 4151 # We don't call objectify(), because the bitwise methods should not 4152 # upgrade/downgrade, even when upgrading/downgrading is enabled. 4153 4154 my ($class, $x, $y, @r) = ref($_[0]) ? (ref($_[0]), @_) : @_; 4155 4156 my $xint = Math::BigInt -> bbrsft($x, $y, @r); 4157 4158 # Temporarily disable downgrading. 4159 4160 my $dng = $class -> downgrade(); 4161 $class -> downgrade(undef); 4162 4163 # Convert to our class without downgrading. 4164 4165 my $xflt = $class -> new($xint); 4166 4167 # Reset downgrading. 4168 4169 $class -> downgrade($dng); 4170 4171 # If we are called as a class method, the first operand might not be an 4172 # object of this class, so check. 4173 4174 if (defined(blessed($x)) && $x -> isa(__PACKAGE__)) { 4175 $x -> {sign} = $xflt -> {sign}; 4176 $x -> {_m} = $xflt -> {_m}; 4177 $x -> {_es} = $xflt -> {_es}; 4178 $x -> {_e} = $xflt -> {_e}; 4179 } else { 4180 $x = $xflt; 4181 } 4182 4183 # Now we might downgrade. 4184 4185 return $downgrade -> new($x) if defined($downgrade); 4186 $x -> round(@r); 4187} 4188 4189sub band { 4190 my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 4191 ? (ref($_[0]), @_) 4192 : objectify(2, @_); 4193 4194 return if $x -> modify('band'); 4195 4196 return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan(); 4197 4198 my $xint = $x -> as_int(); # to Math::BigInt 4199 my $yint = $y -> as_int(); # to Math::BigInt 4200 4201 $xint = $xint -> band($yint); 4202 4203 return $xint -> round(@r) if defined $downgrade; 4204 4205 my $xflt = $class -> new($xint); # back to Math::BigFloat 4206 $x -> {sign} = $xflt -> {sign}; 4207 $x -> {_m} = $xflt -> {_m}; 4208 $x -> {_es} = $xflt -> {_es}; 4209 $x -> {_e} = $xflt -> {_e}; 4210 4211 return $x -> round(@r); 4212} 4213 4214sub bior { 4215 my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 4216 ? (ref($_[0]), @_) 4217 : objectify(2, @_); 4218 4219 return if $x -> modify('bior'); 4220 4221 return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan(); 4222 4223 my $xint = $x -> as_int(); # to Math::BigInt 4224 my $yint = $y -> as_int(); # to Math::BigInt 4225 4226 $xint = $xint -> bior($yint); 4227 4228 return $xint -> round(@r) if defined $downgrade; 4229 4230 my $xflt = $class -> new($xint); # back to Math::BigFloat 4231 $x -> {sign} = $xflt -> {sign}; 4232 $x -> {_m} = $xflt -> {_m}; 4233 $x -> {_es} = $xflt -> {_es}; 4234 $x -> {_e} = $xflt -> {_e}; 4235 4236 return $x -> round(@r); 4237} 4238 4239sub bxor { 4240 my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1]) 4241 ? (ref($_[0]), @_) 4242 : objectify(2, @_); 4243 4244 return if $x -> modify('bxor'); 4245 4246 return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan(); 4247 4248 my $xint = $x -> as_int(); # to Math::BigInt 4249 my $yint = $y -> as_int(); # to Math::BigInt 4250 4251 $xint = $xint -> bxor($yint); 4252 4253 return $xint -> round(@r) if defined $downgrade; 4254 4255 my $xflt = $class -> new($xint); # back to Math::BigFloat 4256 $x -> {sign} = $xflt -> {sign}; 4257 $x -> {_m} = $xflt -> {_m}; 4258 $x -> {_es} = $xflt -> {_es}; 4259 $x -> {_e} = $xflt -> {_e}; 4260 4261 return $x -> round(@r); 4262} 4263 4264sub bnot { 4265 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4266 4267 return if $x -> modify('bnot'); 4268 4269 return $x -> bnan(@r) if $x -> is_nan(); 4270 4271 my $xint = $x -> as_int(); # to Math::BigInt 4272 $xint = $xint -> bnot(); 4273 4274 return $xint -> round(@r) if defined $downgrade; 4275 4276 my $xflt = $class -> new($xint); # back to Math::BigFloat 4277 $x -> {sign} = $xflt -> {sign}; 4278 $x -> {_m} = $xflt -> {_m}; 4279 $x -> {_es} = $xflt -> {_es}; 4280 $x -> {_e} = $xflt -> {_e}; 4281 4282 return $x -> round(@r); 4283} 4284 4285############################################################################### 4286# Rounding methods 4287############################################################################### 4288 4289sub bround { 4290 # accuracy: preserve $N digits, and overwrite the rest with 0's 4291 4292 my ($class, $x, @a) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4293 4294 if (($a[0] || 0) < 0) { 4295 croak('bround() needs positive accuracy'); 4296 } 4297 4298 return $x if $x->modify('bround'); 4299 4300 my ($scale, $mode) = $x->_scale_a(@a); 4301 if (!defined $scale) { # no-op 4302 return $downgrade -> new($x) if defined($downgrade) 4303 && ($x->is_int() || $x->is_inf() || $x->is_nan()); 4304 return $x; 4305 } 4306 4307 # Scale is now either $x->{accuracy}, $accuracy, or the input argument. Test 4308 # whether $x already has lower accuracy, do nothing in this case but do 4309 # round if the accuracy is the same, since a math operation might want to 4310 # round a number with A=5 to 5 digits afterwards again 4311 4312 if (defined $x->{accuracy} && $x->{accuracy} < $scale) { 4313 return $downgrade -> new($x) if defined($downgrade) 4314 && ($x->is_int() || $x->is_inf() || $x->is_nan()); 4315 return $x; 4316 } 4317 4318 # scale < 0 makes no sense 4319 # scale == 0 => keep all digits 4320 # never round a +-inf, NaN 4321 4322 if ($scale <= 0 || $x->{sign} !~ /^[+-]$/) { 4323 return $downgrade -> new($x) if defined($downgrade) 4324 && ($x->is_int() || $x->is_inf() || $x->is_nan()); 4325 return $x; 4326 } 4327 4328 # 1: never round a 0 4329 # 2: if we should keep more digits than the mantissa has, do nothing 4330 if ($x->is_zero() || $LIB->_len($x->{_m}) <= $scale) { 4331 $x->{accuracy} = $scale if !defined $x->{accuracy} || $x->{accuracy} > $scale; 4332 return $downgrade -> new($x) if defined($downgrade) 4333 && ($x->is_int() || $x->is_inf() || $x->is_nan()); 4334 return $x; 4335 } 4336 4337 # pass sign to bround for '+inf' and '-inf' rounding modes 4338 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt'; 4339 4340 $m = $m->bround($scale, $mode); # round mantissa 4341 $x->{_m} = $m->{value}; # get our mantissa back 4342 $x->{accuracy} = $scale; # remember rounding 4343 $x->{precision} = undef; # and clear P 4344 4345 # bnorm() downgrades if necessary, so no need to check whether to downgrade. 4346 $x->bnorm(); # del trailing zeros gen. by bround() 4347} 4348 4349sub bfround { 4350 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.' 4351 # $n == 0 means round to integer 4352 # expects and returns normalized numbers! 4353 4354 my ($class, $x, @p) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4355 4356 return $x if $x->modify('bfround'); # no-op 4357 4358 my ($scale, $mode) = $x->_scale_p(@p); 4359 if (!defined $scale) { 4360 return $downgrade -> new($x) if defined($downgrade) 4361 && ($x->is_int() || $x->is_inf() || $x->is_nan()); 4362 return $x; 4363 } 4364 4365 # never round a 0, +-inf, NaN 4366 4367 if ($x->is_zero()) { 4368 $x->{precision} = $scale if !defined $x->{precision} || $x->{precision} < $scale; # -3 < -2 4369 return $downgrade -> new($x) if defined($downgrade) 4370 && ($x->is_int() || $x->is_inf() || $x->is_nan()); 4371 return $x; 4372 } 4373 4374 if ($x->{sign} !~ /^[+-]$/) { 4375 return $downgrade -> new($x) if defined($downgrade) 4376 && ($x->is_int() || $x->is_inf() || $x->is_nan()); 4377 return $x; 4378 } 4379 4380 # don't round if x already has lower precision 4381 if (defined $x->{precision} && $x->{precision} < 0 && $scale < $x->{precision}) { 4382 return $downgrade -> new($x) if defined($downgrade) 4383 && ($x->is_int() || $x->is_inf() || $x->is_nan()); 4384 return $x; 4385 } 4386 4387 $x->{precision} = $scale; # remember round in any case 4388 $x->{accuracy} = undef; # and clear A 4389 if ($scale < 0) { 4390 # round right from the '.' 4391 4392 if ($x->{_es} eq '+') { # e >= 0 => nothing to round 4393 return $downgrade -> new($x) if defined($downgrade) 4394 && ($x->is_int() || $x->is_inf() || $x->is_nan()); 4395 return $x; 4396 } 4397 4398 $scale = -$scale; # positive for simplicity 4399 my $len = $LIB->_len($x->{_m}); # length of mantissa 4400 4401 # the following poses a restriction on _e, but if _e is bigger than a 4402 # scalar, you got other problems (memory etc) anyway 4403 my $dad = -(0+ ($x->{_es}.$LIB->_num($x->{_e}))); # digits after dot 4404 my $zad = 0; # zeros after dot 4405 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style 4406 4407 # print "scale $scale dad $dad zad $zad len $len\n"; 4408 # number bsstr len zad dad 4409 # 0.123 123e-3 3 0 3 4410 # 0.0123 123e-4 3 1 4 4411 # 0.001 1e-3 1 2 3 4412 # 1.23 123e-2 3 0 2 4413 # 1.2345 12345e-4 5 0 4 4414 4415 # do not round after/right of the $dad 4416 4417 if ($scale > $dad) { # 0.123, scale >= 3 => exit 4418 return $downgrade -> new($x) if defined($downgrade) 4419 && ($x->is_int() || $x->is_inf() || $x->is_nan()); 4420 return $x; 4421 } 4422 4423 # round to zero if rounding inside the $zad, but not for last zero like: 4424 # 0.0065, scale -2, round last '0' with following '65' (scale == zad 4425 # case) 4426 if ($scale < $zad) { 4427 return $downgrade -> new($x) if defined($downgrade) 4428 && ($x->is_int() || $x->is_inf() || $x->is_nan()); 4429 return $x->bzero(); 4430 } 4431 4432 if ($scale == $zad) { # for 0.006, scale -3 and trunc 4433 $scale = -$len; 4434 } else { 4435 # adjust round-point to be inside mantissa 4436 if ($zad != 0) { 4437 $scale = $scale-$zad; 4438 } else { 4439 my $dbd = $len - $dad; 4440 $dbd = 0 if $dbd < 0; # digits before dot 4441 $scale = $dbd+$scale; 4442 } 4443 } 4444 } else { 4445 # round left from the '.' 4446 4447 # 123 => 100 means length(123) = 3 - $scale (2) => 1 4448 4449 my $dbt = $LIB->_len($x->{_m}); 4450 # digits before dot 4451 my $dbd = $dbt + ($x->{_es} . $LIB->_num($x->{_e})); 4452 # should be the same, so treat it as this 4453 $scale = 1 if $scale == 0; 4454 # shortcut if already integer 4455 if ($scale == 1 && $dbt <= $dbd) { 4456 return $downgrade -> new($x) if defined($downgrade) 4457 && ($x->is_int() || $x->is_inf() || $x->is_nan()); 4458 return $x; 4459 } 4460 # maximum digits before dot 4461 ++$dbd; 4462 4463 if ($scale > $dbd) { 4464 # not enough digits before dot, so round to zero 4465 return $downgrade -> new($x) if defined($downgrade); 4466 return $x->bzero; 4467 } elsif ($scale == $dbd) { 4468 # maximum 4469 $scale = -$dbt; 4470 } else { 4471 $scale = $dbd - $scale; 4472 } 4473 } 4474 4475 # pass sign to bround for rounding modes '+inf' and '-inf' 4476 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt'; 4477 $m = $m->bround($scale, $mode); 4478 $x->{_m} = $m->{value}; # get our mantissa back 4479 4480 # bnorm() downgrades if necessary, so no need to check whether to downgrade. 4481 $x->bnorm(); 4482} 4483 4484sub bfloor { 4485 # round towards minus infinity 4486 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4487 4488 return $x if $x->modify('bfloor'); 4489 4490 return $x -> bnan(@r) if $x -> is_nan(); 4491 4492 if ($x->{sign} =~ /^[+-]$/) { 4493 # if $x has digits after dot, remove them 4494 if ($x->{_es} eq '-') { 4495 $x->{_m} = $LIB->_rsft($x->{_m}, $x->{_e}, 10); 4496 $x->{_e} = $LIB->_zero(); 4497 $x->{_es} = '+'; 4498 # increment if negative 4499 $x->{_m} = $LIB->_inc($x->{_m}) if $x->{sign} eq '-'; 4500 } 4501 $x = $x->round(@r); 4502 } 4503 return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade); 4504 return $x; 4505} 4506 4507sub bceil { 4508 # round towards plus infinity 4509 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4510 4511 return $x if $x->modify('bceil'); 4512 4513 return $x -> bnan(@r) if $x -> is_nan(); 4514 4515 # if $x has digits after dot, remove them 4516 if ($x->{sign} =~ /^[+-]$/) { 4517 if ($x->{_es} eq '-') { 4518 $x->{_m} = $LIB->_rsft($x->{_m}, $x->{_e}, 10); 4519 $x->{_e} = $LIB->_zero(); 4520 $x->{_es} = '+'; 4521 if ($x->{sign} eq '+') { 4522 $x->{_m} = $LIB->_inc($x->{_m}); # increment if positive 4523 } else { 4524 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # avoid -0 4525 } 4526 } 4527 $x = $x->round(@r); 4528 } 4529 4530 return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade); 4531 return $x; 4532} 4533 4534sub bint { 4535 # round towards zero 4536 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4537 4538 return $x if $x->modify('bint'); 4539 4540 return $x -> bnan(@r) if $x -> is_nan(); 4541 4542 if ($x->{sign} =~ /^[+-]$/) { 4543 # if $x has digits after the decimal point 4544 if ($x->{_es} eq '-') { 4545 $x->{_m} = $LIB->_rsft($x->{_m}, $x->{_e}, 10); # remove frac part 4546 $x->{_e} = $LIB->_zero(); # truncate/normalize 4547 $x->{_es} = '+'; # abs e 4548 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # avoid -0 4549 } 4550 $x = $x->round(@r); 4551 } 4552 4553 return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade); 4554 return $x; 4555} 4556 4557############################################################################### 4558# Other mathematical methods 4559############################################################################### 4560 4561sub bgcd { 4562 # (BINT or num_str, BINT or num_str) return BINT 4563 # does not modify arguments, but returns new object 4564 4565 # Class::method(...) -> Class->method(...) 4566 unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) || 4567 $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i)) 4568 { 4569 #carp "Using ", (caller(0))[3], "() as a function is deprecated;", 4570 # " use is as a method instead"; 4571 unshift @_, __PACKAGE__; 4572 } 4573 4574 my ($class, @args) = objectify(0, @_); 4575 4576 my $x = shift @args; 4577 $x = defined(blessed($x)) && $x -> isa(__PACKAGE__) ? $x -> copy() 4578 : $class -> new($x); 4579 return $class->bnan() unless $x -> is_int(); 4580 4581 while (@args) { 4582 my $y = shift @args; 4583 $y = $class->new($y) 4584 unless defined(blessed($y)) && $y -> isa(__PACKAGE__); 4585 return $class->bnan() unless $y -> is_int(); 4586 4587 # greatest common divisor 4588 while (! $y->is_zero()) { 4589 ($x, $y) = ($y->copy(), $x->copy()->bmod($y)); 4590 } 4591 4592 last if $x -> is_one(); 4593 } 4594 $x = $x -> babs(); 4595 4596 return $downgrade -> new($x) 4597 if defined $downgrade && $x->is_int(); 4598 return $x; 4599} 4600 4601sub blcm { 4602 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT 4603 # does not modify arguments, but returns new object 4604 # Least Common Multiple 4605 4606 # Class::method(...) -> Class->method(...) 4607 unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) || 4608 $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i)) 4609 { 4610 #carp "Using ", (caller(0))[3], "() as a function is deprecated;", 4611 # " use is as a method instead"; 4612 unshift @_, __PACKAGE__; 4613 } 4614 4615 my ($class, @args) = objectify(0, @_); 4616 4617 my $x = shift @args; 4618 $x = defined(blessed($x)) && $x -> isa(__PACKAGE__) ? $x -> copy() 4619 : $class -> new($x); 4620 return $class->bnan() if $x->{sign} !~ /^[+-]$/; # x NaN? 4621 4622 while (@args) { 4623 my $y = shift @args; 4624 $y = $class -> new($y) 4625 unless defined(blessed($y)) && $y -> isa(__PACKAGE__); 4626 return $x->bnan() unless $y -> is_int(); 4627 my $gcd = $x -> bgcd($y); 4628 $x = $x -> bdiv($gcd) -> bmul($y); 4629 } 4630 4631 $x = $x -> babs(); 4632 4633 return $downgrade -> new($x) 4634 if defined $downgrade && $x->is_int(); 4635 return $x; 4636} 4637 4638############################################################################### 4639# Object property methods 4640############################################################################### 4641 4642sub length { 4643 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4644 4645 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 4646 4647 return 1 if $LIB->_is_zero($x->{_m}); 4648 4649 my $len = $LIB->_len($x->{_m}); 4650 $len += $LIB->_num($x->{_e}) if $x->{_es} eq '+'; 4651 if (wantarray()) { 4652 my $t = 0; 4653 $t = $LIB->_num($x->{_e}) if $x->{_es} eq '-'; 4654 return ($len, $t); 4655 } 4656 $len; 4657} 4658 4659sub mantissa { 4660 # return a copy of the mantissa 4661 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4662 4663 # The following line causes a lot of noise in the test suits for 4664 # the Math-BigRat and bignum distributions. Fixme! 4665 #carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 4666 4667 return $x -> bnan(@r) if $x -> is_nan(); 4668 4669 if ($x->{sign} !~ /^[+-]$/) { 4670 my $s = $x->{sign}; 4671 $s =~ s/^\+//; 4672 return Math::BigInt->new($s, undef, undef); # -inf, +inf => +inf 4673 } 4674 my $m = Math::BigInt->new($LIB->_str($x->{_m}), undef, undef); 4675 $m = $m->bneg() if $x->{sign} eq '-'; 4676 $m; 4677} 4678 4679sub exponent { 4680 # return a copy of the exponent 4681 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4682 4683 # The following line causes a lot of noise in the test suits for 4684 # the Math-BigRat and bignum distributions. Fixme! 4685 #carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 4686 4687 return $x -> bnan(@r) if $x -> is_nan(); 4688 4689 if ($x->{sign} !~ /^[+-]$/) { 4690 my $s = $x->{sign}; 4691 $s =~ s/^[+-]//; 4692 return Math::BigInt->new($s, undef, undef); # -inf, +inf => +inf 4693 } 4694 Math::BigInt->new($x->{_es} . $LIB->_str($x->{_e}), undef, undef); 4695} 4696 4697sub parts { 4698 # return a copy of both the exponent and the mantissa 4699 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4700 4701 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 4702 4703 if ($x->{sign} !~ /^[+-]$/) { 4704 my $s = $x->{sign}; 4705 $s =~ s/^\+//; 4706 my $se = $s; 4707 $se =~ s/^-//; 4708 # +inf => inf and -inf, +inf => inf 4709 return ($class->new($s), $class->new($se)); 4710 } 4711 my $m = Math::BigInt->bzero(); 4712 $m->{value} = $LIB->_copy($x->{_m}); 4713 $m = $m->bneg() if $x->{sign} eq '-'; 4714 ($m, Math::BigInt->new($x->{_es} . $LIB->_num($x->{_e}))); 4715} 4716 4717# Parts used for scientific notation with significand/mantissa and exponent as 4718# integers. E.g., "12345.6789" is returned as "123456789" (mantissa) and "-4" 4719# (exponent). 4720 4721sub sparts { 4722 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4723 4724 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 4725 4726 # Not-a-number. 4727 4728 if ($x -> is_nan()) { 4729 my $mant = $class -> bnan(); # mantissa 4730 return $mant unless wantarray; # scalar context 4731 my $expo = $class -> bnan(); # exponent 4732 return ($mant, $expo); # list context 4733 } 4734 4735 # Infinity. 4736 4737 if ($x -> is_inf()) { 4738 my $mant = $class -> binf($x->{sign}); # mantissa 4739 return $mant unless wantarray; # scalar context 4740 my $expo = $class -> binf('+'); # exponent 4741 return ($mant, $expo); # list context 4742 } 4743 4744 # Finite number. 4745 4746 my $mant = $class -> new($x); 4747 $mant->{_es} = '+'; 4748 $mant->{_e} = $LIB->_zero(); 4749 $mant = $downgrade -> new($mant) if defined $downgrade; 4750 return $mant unless wantarray; 4751 4752 my $expo = $class -> new($x -> {_es} . $LIB->_str($x -> {_e})); 4753 $expo = $downgrade -> new($expo) if defined $downgrade; 4754 return ($mant, $expo); 4755} 4756 4757# Parts used for normalized notation with significand/mantissa as either 0 or a 4758# number in the semi-open interval [1,10). E.g., "12345.6789" is returned as 4759# "1.23456789" and "4". 4760 4761sub nparts { 4762 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4763 4764 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 4765 4766 # Not-a-number and Infinity. 4767 4768 return $x -> sparts() if $x -> is_nan() || $x -> is_inf(); 4769 4770 # Finite number. 4771 4772 my ($mant, $expo) = $x -> sparts(); 4773 4774 if ($mant -> bcmp(0)) { 4775 my ($ndigtot, $ndigfrac) = $mant -> length(); 4776 my $expo10adj = $ndigtot - $ndigfrac - 1; 4777 4778 if ($expo10adj > 0) { # if mantissa is not an integer 4779 $mant = $mant -> brsft($expo10adj, 10); 4780 return $mant unless wantarray; 4781 $expo = $expo -> badd($expo10adj); 4782 return ($mant, $expo); 4783 } 4784 } 4785 4786 return $mant unless wantarray; 4787 return ($mant, $expo); 4788} 4789 4790# Parts used for engineering notation with significand/mantissa as either 0 or a 4791# number in the semi-open interval [1,1000) and the exponent is a multiple of 3. 4792# E.g., "12345.6789" is returned as "12.3456789" and "3". 4793 4794sub eparts { 4795 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4796 4797 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 4798 4799 # Not-a-number and Infinity. 4800 4801 return $x -> sparts() if $x -> is_nan() || $x -> is_inf(); 4802 4803 # Finite number. 4804 4805 my ($mant, $expo) = $x -> nparts(); 4806 4807 my $c = $expo -> copy() -> bmod(3); 4808 $mant = $mant -> blsft($c, 10); 4809 return $mant unless wantarray; 4810 4811 $expo = $expo -> bsub($c); 4812 return ($mant, $expo); 4813} 4814 4815# Parts used for decimal notation, e.g., "12345.6789" is returned as "12345" 4816# (integer part) and "0.6789" (fraction part). 4817 4818sub dparts { 4819 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4820 4821 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 4822 4823 # Not-a-number. 4824 4825 if ($x -> is_nan()) { 4826 my $int = $class -> bnan(); 4827 return $int unless wantarray; 4828 my $frc = $class -> bzero(); # or NaN? 4829 return ($int, $frc); 4830 } 4831 4832 # Infinity. 4833 4834 if ($x -> is_inf()) { 4835 my $int = $class -> binf($x->{sign}); 4836 return $int unless wantarray; 4837 my $frc = $class -> bzero(); 4838 return ($int, $frc); 4839 } 4840 4841 # Finite number. 4842 4843 my $int = $x -> copy(); 4844 my $frc; 4845 4846 # If the input is an integer. 4847 4848 if ($int->{_es} eq '+') { 4849 $frc = $class -> bzero(); 4850 } 4851 4852 # If the input has a fraction part 4853 4854 else { 4855 $int->{_m} = $LIB -> _rsft($int->{_m}, $int->{_e}, 10); 4856 $int->{_e} = $LIB -> _zero(); 4857 $int->{_es} = '+'; 4858 $int->{sign} = '+' if $LIB->_is_zero($int->{_m}); # avoid -0 4859 return $int unless wantarray; 4860 $frc = $x -> copy() -> bsub($int); 4861 return ($int, $frc); 4862 } 4863 4864 $int = $downgrade -> new($int) if defined $downgrade; 4865 return $int unless wantarray; 4866 return $int, $frc; 4867} 4868 4869# Fractional parts with the numerator and denominator as integers. E.g., 4870# "123.4375" is returned as "1975" and "16". 4871 4872sub fparts { 4873 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4874 4875 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 4876 4877 # NaN => NaN/NaN 4878 4879 if ($x -> is_nan()) { 4880 return $class -> bnan() unless wantarray; 4881 return $class -> bnan(), $class -> bnan(); 4882 } 4883 4884 # ±Inf => ±Inf/1 4885 4886 if ($x -> is_inf()) { 4887 my $numer = $class -> binf($x->{sign}); 4888 return $numer unless wantarray; 4889 my $denom = $class -> bone(); 4890 return $numer, $denom; 4891 } 4892 4893 # Finite number. 4894 4895 # If we get here, we know that the output is an integer. 4896 4897 $class = $downgrade if defined $downgrade; 4898 4899 my @flt_parts = ($x->{sign}, $x->{_m}, $x->{_es}, $x->{_e}); 4900 my @rat_parts = $class -> _flt_lib_parts_to_rat_lib_parts(@flt_parts); 4901 my $num = $class -> new($LIB -> _str($rat_parts[1])); 4902 my $den = $class -> new($LIB -> _str($rat_parts[2])); 4903 $num = $num -> bneg() if $rat_parts[0] eq "-"; 4904 return $num unless wantarray; 4905 return $num, $den; 4906} 4907 4908# Given "123.4375", returns "1975", since "123.4375" is "1975/16". 4909 4910sub numerator { 4911 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4912 4913 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 4914 4915 return $class -> bnan() if $x -> is_nan(); 4916 return $class -> binf($x -> sign()) if $x -> is_inf(); 4917 return $class -> bzero() if $x -> is_zero(); 4918 4919 # If we get here, we know that the output is an integer. 4920 4921 $class = $downgrade if defined $downgrade; 4922 4923 if ($x -> {_es} eq '-') { # exponent < 0 4924 my $numer_lib = $LIB -> _copy($x -> {_m}); 4925 my $denom_lib = $LIB -> _1ex($x -> {_e}); 4926 my $gcd_lib = $LIB -> _gcd($LIB -> _copy($numer_lib), $denom_lib); 4927 $numer_lib = $LIB -> _div($numer_lib, $gcd_lib); 4928 return $class -> new($x -> {sign} . $LIB -> _str($numer_lib)); 4929 } 4930 4931 elsif (! $LIB -> _is_zero($x -> {_e})) { # exponent > 0 4932 my $numer_lib = $LIB -> _copy($x -> {_m}); 4933 $numer_lib = $LIB -> _lsft($numer_lib, $x -> {_e}, 10); 4934 return $class -> new($x -> {sign} . $LIB -> _str($numer_lib)); 4935 } 4936 4937 else { # exponent = 0 4938 return $class -> new($x -> {sign} . $LIB -> _str($x -> {_m})); 4939 } 4940} 4941 4942# Given "123.4375", returns "16", since "123.4375" is "1975/16". 4943 4944sub denominator { 4945 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4946 4947 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 4948 4949 return $class -> bnan() if $x -> is_nan(); 4950 4951 # If we get here, we know that the output is an integer. 4952 4953 $class = $downgrade if defined $downgrade; 4954 4955 if ($x -> {_es} eq '-') { # exponent < 0 4956 my $numer_lib = $LIB -> _copy($x -> {_m}); 4957 my $denom_lib = $LIB -> _1ex($x -> {_e}); 4958 my $gcd_lib = $LIB -> _gcd($LIB -> _copy($numer_lib), $denom_lib); 4959 $denom_lib = $LIB -> _div($denom_lib, $gcd_lib); 4960 return $class -> new($LIB -> _str($denom_lib)); 4961 } 4962 4963 else { # exponent >= 0 4964 return $class -> bone(); 4965 } 4966} 4967 4968############################################################################### 4969# String conversion methods 4970############################################################################### 4971 4972sub bstr { 4973 # (ref to BFLOAT or num_str) return num_str 4974 # Convert number from internal format to (non-scientific) string format. 4975 # internal format is always normalized (no leading zeros, "-0" => "+0") 4976 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 4977 4978 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 4979 4980 # Inf and NaN 4981 4982 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 4983 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 4984 return 'inf'; # +inf 4985 } 4986 4987 # Finite number 4988 4989 my $es = '0'; 4990 my $len = 1; 4991 my $cad = 0; 4992 my $dot = '.'; 4993 4994 # $x is zero? 4995 my $not_zero = !($x->{sign} eq '+' && $LIB->_is_zero($x->{_m})); 4996 if ($not_zero) { 4997 $es = $LIB->_str($x->{_m}); 4998 $len = CORE::length($es); 4999 my $e = $LIB->_num($x->{_e}); 5000 $e = -$e if $x->{_es} eq '-'; 5001 if ($e < 0) { 5002 $dot = ''; 5003 # if _e is bigger than a scalar, the following will blow your memory 5004 if ($e <= -$len) { 5005 my $r = abs($e) - $len; 5006 $es = '0.'. ('0' x $r) . $es; 5007 $cad = -($len+$r); 5008 } else { 5009 substr($es, $e, 0) = '.'; 5010 $cad = $LIB->_num($x->{_e}); 5011 $cad = -$cad if $x->{_es} eq '-'; 5012 } 5013 } elsif ($e > 0) { 5014 # expand with zeros 5015 $es .= '0' x $e; 5016 $len += $e; 5017 $cad = 0; 5018 } 5019 } # if not zero 5020 5021 $es = '-'.$es if $x->{sign} eq '-'; 5022 # if set accuracy or precision, pad with zeros on the right side 5023 if ((defined $x->{accuracy}) && ($not_zero)) { 5024 # 123400 => 6, 0.1234 => 4, 0.001234 => 4 5025 my $zeros = $x->{accuracy} - $cad; # cad == 0 => 12340 5026 $zeros = $x->{accuracy} - $len if $cad != $len; 5027 $es .= $dot.'0' x $zeros if $zeros > 0; 5028 } elsif ((($x->{precision} || 0) < 0)) { 5029 # 123400 => 6, 0.1234 => 4, 0.001234 => 6 5030 my $zeros = -$x->{precision} + $cad; 5031 $es .= $dot.'0' x $zeros if $zeros > 0; 5032 } 5033 $es; 5034} 5035 5036# Decimal notation, e.g., "12345.6789" (no exponent). 5037 5038sub bdstr { 5039 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 5040 5041 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5042 5043 # Inf and NaN 5044 5045 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 5046 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 5047 return 'inf'; # +inf 5048 } 5049 5050 # Upgrade? 5051 5052 return $upgrade -> bdstr($x, @r) 5053 if defined($upgrade) && !$x -> isa(__PACKAGE__); 5054 5055 # Finite number 5056 5057 my $mant = $LIB->_str($x->{_m}); 5058 my $esgn = $x->{_es}; 5059 my $eabs = $LIB -> _num($x->{_e}); 5060 5061 my $uintmax = ~0; 5062 5063 my $str = $mant; 5064 if ($esgn eq '+') { 5065 5066 croak("The absolute value of the exponent is too large") 5067 if $eabs > $uintmax; 5068 5069 $str .= "0" x $eabs; 5070 5071 } else { 5072 my $mlen = CORE::length($mant); 5073 my $c = $mlen - $eabs; 5074 5075 my $intmax = ($uintmax - 1) / 2; 5076 croak("The absolute value of the exponent is too large") 5077 if (1 - $c) > $intmax; 5078 5079 $str = "0" x (1 - $c) . $str if $c <= 0; 5080 substr($str, -$eabs, 0) = '.'; 5081 } 5082 5083 return $x->{sign} eq '-' ? '-' . $str : $str; 5084} 5085 5086# Scientific notation with significand/mantissa and exponent as integers, e.g., 5087# "12345.6789" is written as "123456789e-4". 5088 5089sub bsstr { 5090 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 5091 5092 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5093 5094 # Inf and NaN 5095 5096 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 5097 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 5098 return 'inf'; # +inf 5099 } 5100 5101 # Upgrade? 5102 5103 return $upgrade -> bsstr($x, @r) 5104 if defined($upgrade) && !$x -> isa(__PACKAGE__); 5105 5106 # Finite number 5107 5108 ($x->{sign} eq '-' ? '-' : '') . $LIB->_str($x->{_m}) 5109 . 'e' . $x->{_es} . $LIB->_str($x->{_e}); 5110} 5111 5112# Normalized notation, e.g., "12345.6789" is written as "1.23456789e+4". 5113 5114sub bnstr { 5115 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 5116 5117 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5118 5119 # Inf and NaN 5120 5121 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 5122 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 5123 return 'inf'; # +inf 5124 } 5125 5126 # Upgrade? 5127 5128 return $upgrade -> bnstr($x, @r) 5129 if defined($upgrade) && !$x -> isa(__PACKAGE__); 5130 5131 # Finite number 5132 5133 my $str = $x->{sign} eq '-' ? '-' : ''; 5134 5135 # Get the mantissa and the length of the mantissa. 5136 5137 my $mant = $LIB->_str($x->{_m}); 5138 my $mantlen = CORE::length($mant); 5139 5140 if ($mantlen == 1) { 5141 5142 # Not decimal point when the mantissa has length one, i.e., return the 5143 # number 2 as the string "2", not "2.". 5144 5145 $str .= $mant . 'e' . $x->{_es} . $LIB->_str($x->{_e}); 5146 5147 } else { 5148 5149 # Compute new exponent where the original exponent is adjusted by the 5150 # length of the mantissa minus one (because the decimal point is after 5151 # one digit). 5152 5153 my ($eabs, $esgn) = $LIB -> _sadd($LIB -> _copy($x->{_e}), $x->{_es}, 5154 $LIB -> _new($mantlen - 1), "+"); 5155 substr $mant, 1, 0, "."; 5156 $str .= $mant . 'e' . $esgn . $LIB->_str($eabs); 5157 5158 } 5159 5160 return $str; 5161} 5162 5163# Engineering notation, e.g., "12345.6789" is written as "12.3456789e+3". 5164 5165sub bestr { 5166 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 5167 5168 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5169 5170 # Inf and NaN 5171 5172 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 5173 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 5174 return 'inf'; # +inf 5175 } 5176 5177 # Upgrade? 5178 5179 return $upgrade -> bestr($x, @r) 5180 if defined($upgrade) && !$x -> isa(__PACKAGE__); 5181 5182 # Finite number 5183 5184 my $str = $x->{sign} eq '-' ? '-' : ''; 5185 5186 # Get the mantissa, the length of the mantissa, and adjust the exponent by 5187 # the length of the mantissa minus 1 (because the dot is after one digit). 5188 5189 my $mant = $LIB->_str($x->{_m}); 5190 my $mantlen = CORE::length($mant); 5191 my ($eabs, $esgn) = $LIB -> _sadd($LIB -> _copy($x->{_e}), $x->{_es}, 5192 $LIB -> _new($mantlen - 1), "+"); 5193 5194 my $dotpos = 1; 5195 my $mod = $LIB -> _mod($LIB -> _copy($eabs), $LIB -> _new("3")); 5196 unless ($LIB -> _is_zero($mod)) { 5197 if ($esgn eq '+') { 5198 $eabs = $LIB -> _sub($eabs, $mod); 5199 $dotpos += $LIB -> _num($mod); 5200 } else { 5201 my $delta = $LIB -> _sub($LIB -> _new("3"), $mod); 5202 $eabs = $LIB -> _add($eabs, $delta); 5203 $dotpos += $LIB -> _num($delta); 5204 } 5205 } 5206 5207 if ($dotpos < $mantlen) { 5208 substr $mant, $dotpos, 0, "."; 5209 } elsif ($dotpos > $mantlen) { 5210 $mant .= "0" x ($dotpos - $mantlen); 5211 } 5212 5213 $str .= $mant . 'e' . $esgn . $LIB->_str($eabs); 5214 5215 return $str; 5216} 5217 5218# Fractional notation, e.g., "123.4375" is written as "1975/16". 5219 5220sub bfstr { 5221 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 5222 5223 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5224 5225 # Inf and NaN 5226 5227 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 5228 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 5229 return 'inf'; # +inf 5230 } 5231 5232 # Upgrade? 5233 5234 return $upgrade -> bfstr($x, @r) 5235 if defined($upgrade) && !$x -> isa(__PACKAGE__); 5236 5237 # Finite number 5238 5239 my $str = $x->{sign} eq '-' ? '-' : ''; 5240 5241 if ($x->{_es} eq '+') { 5242 $str .= $LIB -> _str($x->{_m}) . ("0" x $LIB -> _num($x->{_e})); 5243 } else { 5244 my @flt_parts = ($x->{sign}, $x->{_m}, $x->{_es}, $x->{_e}); 5245 my @rat_parts = $class -> _flt_lib_parts_to_rat_lib_parts(@flt_parts); 5246 $str = $LIB -> _str($rat_parts[1]) . "/" . $LIB -> _str($rat_parts[2]); 5247 $str = "-" . $str if $rat_parts[0] eq "-"; 5248 } 5249 5250 return $str; 5251} 5252 5253sub to_hex { 5254 # return number as hexadecimal string (only for integers defined) 5255 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 5256 5257 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5258 5259 # Inf and NaN 5260 5261 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 5262 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 5263 return 'inf'; # +inf 5264 } 5265 5266 # Upgrade? 5267 5268 return $upgrade -> to_hex($x, @r) 5269 if defined($upgrade) && !$x -> isa(__PACKAGE__); 5270 5271 # Finite number 5272 5273 return '0' if $x->is_zero(); 5274 5275 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex? 5276 5277 my $z = $LIB->_copy($x->{_m}); 5278 if (! $LIB->_is_zero($x->{_e})) { # > 0 5279 $z = $LIB->_lsft($z, $x->{_e}, 10); 5280 } 5281 my $str = $LIB->_to_hex($z); 5282 return $x->{sign} eq '-' ? "-$str" : $str; 5283} 5284 5285sub to_oct { 5286 # return number as octal digit string (only for integers defined) 5287 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 5288 5289 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5290 5291 # Inf and NaN 5292 5293 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 5294 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 5295 return 'inf'; # +inf 5296 } 5297 5298 # Upgrade? 5299 5300 return $upgrade -> to_oct($x, @r) 5301 if defined($upgrade) && !$x -> isa(__PACKAGE__); 5302 5303 # Finite number 5304 5305 return '0' if $x->is_zero(); 5306 5307 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in octal? 5308 5309 my $z = $LIB->_copy($x->{_m}); 5310 if (! $LIB->_is_zero($x->{_e})) { # > 0 5311 $z = $LIB->_lsft($z, $x->{_e}, 10); 5312 } 5313 my $str = $LIB->_to_oct($z); 5314 return $x->{sign} eq '-' ? "-$str" : $str; 5315} 5316 5317sub to_bin { 5318 # return number as binary digit string (only for integers defined) 5319 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_); 5320 5321 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5322 5323 # Inf and NaN 5324 5325 if ($x->{sign} ne '+' && $x->{sign} ne '-') { 5326 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN 5327 return 'inf'; # +inf 5328 } 5329 5330 # Upgrade? 5331 5332 return $upgrade -> to_bin($x, @r) 5333 if defined($upgrade) && !$x -> isa(__PACKAGE__); 5334 5335 # Finite number 5336 5337 return '0' if $x->is_zero(); 5338 5339 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in binary? 5340 5341 my $z = $LIB->_copy($x->{_m}); 5342 if (! $LIB->_is_zero($x->{_e})) { # > 0 5343 $z = $LIB->_lsft($z, $x->{_e}, 10); 5344 } 5345 my $str = $LIB->_to_bin($z); 5346 return $x->{sign} eq '-' ? "-$str" : $str; 5347} 5348 5349sub to_ieee754 { 5350 my ($class, $x, $format, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_); 5351 5352 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5353 5354 my $enc; # significand encoding (applies only to decimal) 5355 my $k; # storage width in bits 5356 my $b; # base 5357 5358 if ($format =~ /^binary(\d+)\z/) { 5359 $k = $1; 5360 $b = 2; 5361 } elsif ($format =~ /^decimal(\d+)(dpd|bcd)?\z/) { 5362 $k = $1; 5363 $b = 10; 5364 $enc = $2 || 'dpd'; # default is dencely-packed decimals (DPD) 5365 } elsif ($format eq 'half') { 5366 $k = 16; 5367 $b = 2; 5368 } elsif ($format eq 'single') { 5369 $k = 32; 5370 $b = 2; 5371 } elsif ($format eq 'double') { 5372 $k = 64; 5373 $b = 2; 5374 } elsif ($format eq 'quadruple') { 5375 $k = 128; 5376 $b = 2; 5377 } elsif ($format eq 'octuple') { 5378 $k = 256; 5379 $b = 2; 5380 } elsif ($format eq 'sexdecuple') { 5381 $k = 512; 5382 $b = 2; 5383 } 5384 5385 if ($b == 2) { 5386 5387 # Get the parameters for this format. 5388 5389 my $p; # precision (in bits) 5390 my $t; # number of bits in significand 5391 my $w; # number of bits in exponent 5392 5393 if ($k == 16) { # binary16 (half-precision) 5394 $p = 11; 5395 $t = 10; 5396 $w = 5; 5397 } elsif ($k == 32) { # binary32 (single-precision) 5398 $p = 24; 5399 $t = 23; 5400 $w = 8; 5401 } elsif ($k == 64) { # binary64 (double-precision) 5402 $p = 53; 5403 $t = 52; 5404 $w = 11; 5405 } else { # binaryN (quadruple-precition and above) 5406 if ($k < 128 || $k != 32 * sprintf('%.0f', $k / 32)) { 5407 croak "Number of bits must be 16, 32, 64, or >= 128 and", 5408 " a multiple of 32"; 5409 } 5410 $p = $k - sprintf('%.0f', 4 * log($k) / log(2)) + 13; 5411 $t = $p - 1; 5412 $w = $k - $t - 1; 5413 } 5414 5415 # The maximum exponent, minimum exponent, and exponent bias. 5416 5417 my $emax = $class -> new(2) -> bpow($w - 1) -> bdec(); 5418 my $emin = 1 - $emax; 5419 my $bias = $emax; 5420 5421 # Get numerical sign, exponent, and mantissa/significand for bit 5422 # string. 5423 5424 my $sign = 0; 5425 my $expo; 5426 my $mant; 5427 5428 if ($x -> is_nan()) { # nan 5429 $sign = 1; 5430 $expo = $emax -> copy() -> binc(); 5431 $mant = $class -> new(2) -> bpow($t - 1); 5432 } elsif ($x -> is_inf()) { # inf 5433 $sign = 1 if $x -> is_neg(); 5434 $expo = $emax -> copy() -> binc(); 5435 $mant = $class -> bzero(); 5436 } elsif ($x -> is_zero()) { # zero 5437 $expo = $emin -> copy() -> bdec(); 5438 $mant = $class -> bzero(); 5439 } else { # normal and subnormal 5440 5441 $sign = 1 if $x -> is_neg(); 5442 5443 # Now we need to compute the mantissa and exponent in base $b. 5444 5445 my $binv = $class -> new("0.5"); 5446 my $b = $class -> new(2); 5447 my $one = $class -> bone(); 5448 5449 # We start off by initializing the exponent to zero and the 5450 # mantissa to the input value. Then we increase the mantissa and 5451 # decrease the exponent, or vice versa, until the mantissa is in 5452 # the desired range or we hit one of the limits for the exponent. 5453 5454 $mant = $x -> copy() -> babs(); 5455 5456 # We need to find the base 2 exponent. First make an estimate of 5457 # the base 2 exponent, before adjusting it below. We could skip 5458 # this estimation and go straight to the while-loops below, but the 5459 # loops are slow, especially when the final exponent is far from 5460 # zero and even more so if the number of digits is large. This 5461 # initial estimation speeds up the computation dramatically. 5462 # 5463 # log2($m * 10**$e) = log10($m + 10**$e) * log(10)/log(2) 5464 # = (log10($m) + $e) * log(10)/log(2) 5465 # = (log($m)/log(10) + $e) * log(10)/log(2) 5466 5467 my ($m, $e) = $x -> nparts(); 5468 my $ms = $m -> numify(); 5469 my $es = $e -> numify(); 5470 5471 my $expo_est = (log(abs($ms))/log(10) + $es) * log(10)/log(2); 5472 $expo_est = int($expo_est); 5473 5474 # Limit the exponent. 5475 5476 if ($expo_est > $emax) { 5477 $expo_est = $emax; 5478 } elsif ($expo_est < $emin) { 5479 $expo_est = $emin; 5480 } 5481 5482 # Don't multiply by a number raised to a negative exponent. This 5483 # will cause a division, whose result is truncated to some fixed 5484 # number of digits. Instead, multiply by the inverse number raised 5485 # to a positive exponent. 5486 5487 $expo = $class -> new($expo_est); 5488 if ($expo_est > 0) { 5489 $mant = $mant -> bmul($binv -> copy() -> bpow($expo)); 5490 } elsif ($expo_est < 0) { 5491 my $expo_abs = $expo -> copy() -> bneg(); 5492 $mant = $mant -> bmul($b -> copy() -> bpow($expo_abs)); 5493 } 5494 5495 # Final adjustment of the estimate above. 5496 5497 while ($mant >= $b && $expo <= $emax) { 5498 $mant = $mant -> bmul($binv); 5499 $expo = $expo -> binc(); 5500 } 5501 5502 while ($mant < $one && $expo >= $emin) { 5503 $mant = $mant -> bmul($b); 5504 $expo = $expo -> bdec(); 5505 } 5506 5507 # This is when the magnitude is larger than what can be represented 5508 # in this format. Encode as infinity. 5509 5510 if ($expo > $emax) { 5511 $mant = $class -> bzero(); 5512 $expo = $emax -> copy() -> binc(); 5513 } 5514 5515 # This is when the magnitude is so small that the number is encoded 5516 # as a subnormal number. 5517 # 5518 # If the magnitude is smaller than that of the smallest subnormal 5519 # number, and rounded downwards, it is encoded as zero. This works 5520 # transparently and does not need to be treated as a special case. 5521 # 5522 # If the number is between the largest subnormal number and the 5523 # smallest normal number, and the value is rounded upwards, the 5524 # value must be encoded as a normal number. This must be treated as 5525 # a special case. 5526 5527 elsif ($expo < $emin) { 5528 5529 # Scale up the mantissa (significand), and round to integer. 5530 5531 my $const = $class -> new($b) -> bpow($t - 1); 5532 $mant = $mant -> bmul($const); 5533 $mant = $mant -> bfround(0); 5534 5535 # If the mantissa overflowed, encode as the smallest normal 5536 # number. 5537 5538 if ($mant == $const -> bmul($b)) { 5539 $mant = $mant -> bzero(); 5540 $expo = $expo -> binc(); 5541 } 5542 } 5543 5544 # This is when the magnitude is within the range of what can be 5545 # encoded as a normal number. 5546 5547 else { 5548 5549 # Remove implicit leading bit, scale up the mantissa 5550 # (significand) to an integer, and round. 5551 5552 $mant = $mant -> bdec(); 5553 my $const = $class -> new($b) -> bpow($t); 5554 $mant = $mant -> bmul($const) -> bfround(0); 5555 5556 # If the mantissa overflowed, encode as the next larger value. 5557 # This works correctly also when the next larger value is 5558 # infinity. 5559 5560 if ($mant == $const) { 5561 $mant = $mant -> bzero(); 5562 $expo = $expo -> binc(); 5563 } 5564 } 5565 } 5566 5567 $expo = $expo -> badd($bias); # add bias 5568 5569 my $signbit = "$sign"; 5570 5571 my $mantbits = $mant -> to_bin(); 5572 $mantbits = ("0" x ($t - CORE::length($mantbits))) . $mantbits; 5573 5574 my $expobits = $expo -> to_bin(); 5575 $expobits = ("0" x ($w - CORE::length($expobits))) . $expobits; 5576 5577 my $bin = $signbit . $expobits . $mantbits; 5578 return pack "B*", $bin; 5579 } 5580 5581 croak("The format '$format' is not yet supported."); 5582} 5583 5584sub as_hex { 5585 # return number as hexadecimal string (only for integers defined) 5586 5587 my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_); 5588 5589 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5590 5591 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc 5592 return '0x0' if $x->is_zero(); 5593 5594 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex? 5595 5596 my $z = $LIB->_copy($x->{_m}); 5597 if (! $LIB->_is_zero($x->{_e})) { # > 0 5598 $z = $LIB->_lsft($z, $x->{_e}, 10); 5599 } 5600 my $str = $LIB->_as_hex($z); 5601 return $x->{sign} eq '-' ? "-$str" : $str; 5602} 5603 5604sub as_oct { 5605 # return number as octal digit string (only for integers defined) 5606 5607 my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_); 5608 5609 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5610 5611 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc 5612 return '00' if $x->is_zero(); 5613 5614 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in octal? 5615 5616 my $z = $LIB->_copy($x->{_m}); 5617 if (! $LIB->_is_zero($x->{_e})) { # > 0 5618 $z = $LIB->_lsft($z, $x->{_e}, 10); 5619 } 5620 my $str = $LIB->_as_oct($z); 5621 return $x->{sign} eq '-' ? "-$str" : $str; 5622} 5623 5624sub as_bin { 5625 # return number as binary digit string (only for integers defined) 5626 5627 my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_); 5628 5629 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5630 5631 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc 5632 return '0b0' if $x->is_zero(); 5633 5634 return $nan if $x->{_es} ne '+'; # how to do 1e-1 in binary? 5635 5636 my $z = $LIB->_copy($x->{_m}); 5637 if (! $LIB->_is_zero($x->{_e})) { # > 0 5638 $z = $LIB->_lsft($z, $x->{_e}, 10); 5639 } 5640 my $str = $LIB->_as_bin($z); 5641 return $x->{sign} eq '-' ? "-$str" : $str; 5642} 5643 5644sub numify { 5645 # Make a Perl scalar number from a Math::BigFloat object. 5646 5647 my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_); 5648 5649 carp "Rounding is not supported for ", (caller(0))[3], "()" if @r; 5650 5651 if ($x -> is_nan()) { 5652 require Math::Complex; 5653 my $inf = $Math::Complex::Inf; 5654 return $inf - $inf; 5655 } 5656 5657 if ($x -> is_inf()) { 5658 require Math::Complex; 5659 my $inf = $Math::Complex::Inf; 5660 return $x -> is_negative() ? -$inf : $inf; 5661 } 5662 5663 # Create a string and let Perl's atoi()/atof() handle the rest. 5664 5665 return 0 + $x -> bnstr(); 5666} 5667 5668############################################################################### 5669# Private methods and functions. 5670############################################################################### 5671 5672sub import { 5673 my $class = shift; 5674 $IMPORT++; # remember we did import() 5675 my @a; # unrecognized arguments 5676 5677 my @import = (); 5678 5679 while (@_) { 5680 my $param = shift; 5681 5682 # Enable overloading of constants. 5683 5684 if ($param eq ':constant') { 5685 overload::constant 5686 5687 integer => sub { 5688 $class -> new(shift); 5689 }, 5690 5691 float => sub { 5692 $class -> new(shift); 5693 }, 5694 5695 binary => sub { 5696 # E.g., a literal 0377 shall result in an object whose value 5697 # is decimal 255, but new("0377") returns decimal 377. 5698 return $class -> from_oct($_[0]) if $_[0] =~ /^0_*[0-7]/; 5699 $class -> new(shift); 5700 }; 5701 next; 5702 } 5703 5704 # Upgrading. 5705 5706 if ($param eq 'upgrade') { 5707 $class -> upgrade(shift); 5708 next; 5709 } 5710 5711 # Downgrading. 5712 5713 if ($param eq 'downgrade') { 5714 $class -> downgrade(shift); 5715 next; 5716 } 5717 5718 # Accuracy. 5719 5720 if ($param eq 'accuracy') { 5721 $class -> accuracy(shift); 5722 next; 5723 } 5724 5725 # Precision. 5726 5727 if ($param eq 'precision') { 5728 $class -> precision(shift); 5729 next; 5730 } 5731 5732 # Rounding mode. 5733 5734 if ($param eq 'round_mode') { 5735 $class -> round_mode(shift); 5736 next; 5737 } 5738 5739 # Fall-back accuracy. 5740 5741 if ($param eq 'div_scale') { 5742 $class -> div_scale(shift); 5743 next; 5744 } 5745 5746 # Backend library. 5747 5748 if ($param =~ /^(lib|try|only)\z/) { 5749 push @import, $param; 5750 push @import, shift() if @_; 5751 next; 5752 } 5753 5754 if ($param eq 'with') { 5755 # alternative class for our private parts() 5756 # XXX: no longer supported 5757 # $LIB = shift() || 'Calc'; 5758 # carp "'with' is no longer supported, use 'lib', 'try', or 'only'"; 5759 shift; 5760 next; 5761 } 5762 5763 # Unrecognized parameter. 5764 5765 push @a, $param; 5766 } 5767 5768 Math::BigInt -> import(@import); 5769 5770 # find out which library was actually loaded 5771 $LIB = Math::BigInt -> config('lib'); 5772 5773 $class -> SUPER::import(@a); # for subclasses 5774 $class -> export_to_level(1, $class, @a) if @a; # need this, too 5775} 5776 5777sub _len_to_steps { 5778 # Given D (digits in decimal), compute N so that N! (N factorial) is 5779 # at least D digits long. D should be at least 50. 5780 my $d = shift; 5781 5782 # two constants for the Ramanujan estimate of ln(N!) 5783 my $lg2 = log(2 * 3.14159265) / 2; 5784 my $lg10 = log(10); 5785 5786 # D = 50 => N => 42, so L = 40 and R = 50 5787 my $l = 40; 5788 my $r = $d; 5789 5790 # Otherwise this does not work under -Mbignum and we do not yet have "no 5791 # bignum;" :( 5792 $l = $l->numify if ref($l); 5793 $r = $r->numify if ref($r); 5794 $lg2 = $lg2->numify if ref($lg2); 5795 $lg10 = $lg10->numify if ref($lg10); 5796 5797 # binary search for the right value (could this be written as the reverse of 5798 # lg(n!)?) 5799 while ($r - $l > 1) { 5800 my $n = int(($r - $l) / 2) + $l; 5801 my $ramanujan 5802 = int(($n * log($n) - $n + log($n * (1 + 4*$n*(1+2*$n))) / 6 + $lg2) 5803 / $lg10); 5804 $ramanujan > $d ? $r = $n : $l = $n; 5805 } 5806 $l; 5807} 5808 5809sub _log { 5810 # internal log function to calculate ln() based on Taylor series. 5811 # Modifies $x in place. 5812 my ($x, $scale) = @_; 5813 my $class = ref $x; 5814 5815 # in case of $x == 1, result is 0 5816 return $x -> bzero() if $x -> is_one(); 5817 5818 # XXX TODO: rewrite this in a similar manner to bexp() 5819 5820 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log 5821 5822 # u = x-1, v = x+1 5823 # _ _ 5824 # Taylor: | u 1 u^3 1 u^5 | 5825 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0 5826 # |_ v 3 v^3 5 v^5 _| 5827 5828 # This takes much more steps to calculate the result and is thus not used 5829 # u = x-1 5830 # _ _ 5831 # Taylor: | u 1 u^2 1 u^3 | 5832 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2 5833 # |_ x 2 x^2 3 x^3 _| 5834 5835 # scale used in intermediate computations 5836 my $scaleup = $scale + 4; 5837 5838 my ($v, $u, $numer, $denom, $factor, $f); 5839 5840 $v = $x -> copy(); 5841 $v = $v -> binc(); # v = x+1 5842 $x = $x -> bdec(); 5843 $u = $x -> copy(); # u = x-1; x = x-1 5844 5845 $x = $x -> bdiv($v, $scaleup); # first term: u/v 5846 5847 $numer = $u -> copy(); # numerator 5848 $denom = $v -> copy(); # denominator 5849 5850 $u = $u -> bmul($u); # u^2 5851 $v = $v -> bmul($v); # v^2 5852 5853 $numer = $numer -> bmul($u); # u^3 5854 $denom = $denom -> bmul($v); # v^3 5855 5856 $factor = $class -> new(3); 5857 $f = $class -> new(2); 5858 5859 while (1) { 5860 my $next = $numer -> copy() -> bround($scaleup) 5861 -> bdiv($denom -> copy() -> bmul($factor) -> bround($scaleup), $scaleup); 5862 5863 $next->{accuracy} = undef; 5864 $next->{precision} = undef; 5865 my $x_prev = $x -> copy(); 5866 $x = $x -> badd($next); 5867 5868 last if $x -> bacmp($x_prev) == 0; 5869 5870 # calculate things for the next term 5871 $numer = $numer -> bmul($u); 5872 $denom = $denom -> bmul($v); 5873 $factor = $factor -> badd($f); 5874 } 5875 5876 $x = $x -> bmul($f); # $x *= 2 5877 $x = $x -> bround($scale); 5878} 5879 5880sub _log_10 { 5881 # Internal log function based on reducing input to the range of 0.1 .. 9.99 5882 # and then "correcting" the result to the proper one. Modifies $x in place. 5883 my ($x, $scale) = @_; 5884 my $class = ref $x; 5885 5886 # Taking blog() from numbers greater than 10 takes a *very long* time, so we 5887 # break the computation down into parts based on the observation that: 5888 # blog(X*Y) = blog(X) + blog(Y) 5889 # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller 5890 # $x is the faster it gets. Since 2*$x takes about 10 times as 5891 # long, we make it faster by about a factor of 100 by dividing $x by 10. 5892 5893 # The same observation is valid for numbers smaller than 0.1, e.g. computing 5894 # log(1) is fastest, and the further away we get from 1, the longer it 5895 # takes. So we also 'break' this down by multiplying $x with 10 and subtract 5896 # the log(10) afterwards to get the correct result. 5897 5898 # To get $x even closer to 1, we also divide by 2 and then use log(2) to 5899 # correct for this. For instance if $x is 2.4, we use the formula: 5900 # blog(2.4 * 2) == blog(1.2) + blog(2) 5901 # and thus calculate only blog(1.2) and blog(2), which is faster in total 5902 # than calculating blog(2.4). 5903 5904 # In addition, the values for blog(2) and blog(10) are cached. 5905 5906 # Calculate the number of digits before the dot, i.e., 1 + floor(log10(x)): 5907 # x = 123 => dbd = 3 5908 # x = 1.23 => dbd = 1 5909 # x = 0.0123 => dbd = -1 5910 # x = 0.000123 => dbd = -3 5911 # etc. 5912 5913 my $dbd = $LIB->_num($x->{_e}); 5914 $dbd = -$dbd if $x->{_es} eq '-'; 5915 $dbd += $LIB->_len($x->{_m}); 5916 5917 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid 5918 # infinite recursion 5919 5920 my $calc = 1; # do some calculation? 5921 5922 # No upgrading or downgrading in the intermediate computations. 5923 5924 my $upg = $class -> upgrade(); 5925 my $dng = $class -> downgrade(); 5926 $class -> upgrade(undef); 5927 $class -> downgrade(undef); 5928 5929 # disable the shortcut for 10, since we need log(10) and this would recurse 5930 # infinitely deep 5931 if ($x->{_es} eq '+' && # $x == 10 5932 ($LIB->_is_one($x->{_e}) && 5933 $LIB->_is_one($x->{_m}))) 5934 { 5935 $dbd = 0; # disable shortcut 5936 # we can use the cached value in these cases 5937 if ($scale <= $LOG_10_A) { 5938 $x = $x->bzero(); 5939 $x = $x->badd($LOG_10); # modify $x in place 5940 $calc = 0; # no need to calc, but round 5941 } 5942 # if we can't use the shortcut, we continue normally 5943 } else { 5944 # disable the shortcut for 2, since we maybe have it cached 5945 if (($LIB->_is_zero($x->{_e}) && # $x == 2 5946 $LIB->_is_two($x->{_m}))) 5947 { 5948 $dbd = 0; # disable shortcut 5949 # we can use the cached value in these cases 5950 if ($scale <= $LOG_2_A) { 5951 $x = $x->bzero(); 5952 $x = $x->badd($LOG_2); # modify $x in place 5953 $calc = 0; # no need to calc, but round 5954 } 5955 # if we can't use the shortcut, we continue normally 5956 } 5957 } 5958 5959 # if $x = 0.1, we know the result must be 0-log(10) 5960 if ($calc != 0 && 5961 ($x->{_es} eq '-' && # $x == 0.1 5962 ($LIB->_is_one($x->{_e}) && 5963 $LIB->_is_one($x->{_m})))) 5964 { 5965 $dbd = 0; # disable shortcut 5966 # we can use the cached value in these cases 5967 if ($scale <= $LOG_10_A) { 5968 $x = $x->bzero(); 5969 $x = $x->bsub($LOG_10); 5970 $calc = 0; # no need to calc, but round 5971 } 5972 } 5973 5974 return $x if $calc == 0; # already have the result 5975 5976 # default: these correction factors are undef and thus not used 5977 my $l_10; # value of ln(10) to A of $scale 5978 my $l_2; # value of ln(2) to A of $scale 5979 5980 my $two = $class->new(2); 5981 5982 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1 5983 # so don't do this shortcut for 1 or 0 5984 if (($dbd > 1) || ($dbd < 0)) { 5985 # convert our cached value to an object if not already (avoid doing this 5986 # at import() time, since not everybody needs this) 5987 $LOG_10 = $class->new($LOG_10, undef, undef) unless ref $LOG_10; 5988 5989 # got more than one digit before the dot, or more than one zero after 5990 # the dot, so do: 5991 # log(123) == log(1.23) + log(10) * 2 5992 # log(0.0123) == log(1.23) - log(10) * 2 5993 5994 if ($scale <= $LOG_10_A) { 5995 # use cached value 5996 $l_10 = $LOG_10->copy(); # copy for mul 5997 } else { 5998 # else: slower, compute and cache result 5999 6000 # shorten the time to calculate log(10) based on the following: 6001 # log(1.25 * 8) = log(1.25) + log(8) 6002 # = log(1.25) + log(2) + log(2) + log(2) 6003 6004 # first get $l_2 (and possible compute and cache log(2)) 6005 $LOG_2 = $class->new($LOG_2, undef, undef) unless ref $LOG_2; 6006 if ($scale <= $LOG_2_A) { 6007 # use cached value 6008 $l_2 = $LOG_2->copy(); # copy() for the mul below 6009 } else { 6010 # else: slower, compute and cache result 6011 $l_2 = $two->copy(); 6012 $l_2 = $l_2->_log($scale); # scale+4, actually 6013 $LOG_2 = $l_2->copy(); # cache the result for later 6014 # the copy() is for mul below 6015 $LOG_2_A = $scale; 6016 } 6017 6018 # now calculate log(1.25): 6019 $l_10 = $class->new('1.25'); 6020 $l_10 = $l_10->_log($scale); # scale+4, actually 6021 6022 # log(1.25) + log(2) + log(2) + log(2): 6023 $l_10 = $l_10->badd($l_2); 6024 $l_10 = $l_10->badd($l_2); 6025 $l_10 = $l_10->badd($l_2); 6026 $LOG_10 = $l_10->copy(); # cache the result for later 6027 # the copy() is for mul below 6028 $LOG_10_A = $scale; 6029 } 6030 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1 6031 $l_10 = $l_10->bmul($class->new($dbd)); # log(10) * (digits_before_dot-1) 6032 my $dbd_sign = '+'; 6033 if ($dbd < 0) { 6034 $dbd = -$dbd; 6035 $dbd_sign = '-'; 6036 } 6037 ($x->{_e}, $x->{_es}) = 6038 $LIB -> _ssub($x->{_e}, $x->{_es}, $LIB->_new($dbd), $dbd_sign); 6039 } 6040 6041 # Now: 0.1 <= $x < 10 (and possible correction in l_10) 6042 6043 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div 6044 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1) 6045 6046 $HALF = $class->new($HALF) unless ref($HALF); 6047 6048 my $twos = 0; # default: none (0 times) 6049 while ($x->bacmp($HALF) <= 0) { # X <= 0.5 6050 $twos--; 6051 $x = $x->bmul($two); 6052 } 6053 while ($x->bacmp($two) >= 0) { # X >= 2 6054 $twos++; 6055 $x = $x->bdiv($two, $scale+4); # keep all digits 6056 } 6057 $x = $x->bround($scale+4); 6058 # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both) 6059 # So calculate correction factor based on ln(2): 6060 if ($twos != 0) { 6061 $LOG_2 = $class->new($LOG_2, undef, undef) unless ref $LOG_2; 6062 if ($scale <= $LOG_2_A) { 6063 # use cached value 6064 $l_2 = $LOG_2->copy(); # copy() for the mul below 6065 } else { 6066 # else: slower, compute and cache result 6067 $l_2 = $two->copy(); 6068 $l_2 = $l_2->_log($scale); # scale+4, actually 6069 $LOG_2 = $l_2->copy(); # cache the result for later 6070 # the copy() is for mul below 6071 $LOG_2_A = $scale; 6072 } 6073 $l_2 = $l_2->bmul($twos); # * -2 => subtract, * 2 => add 6074 } else { 6075 undef $l_2; 6076 } 6077 6078 $x = $x->_log($scale); # need to do the "normal" way 6079 $x = $x->badd($l_10) if defined $l_10; # correct it by ln(10) 6080 $x = $x->badd($l_2) if defined $l_2; # and maybe by ln(2) 6081 6082 # Restore globals 6083 6084 $class -> upgrade($upg); 6085 $class -> downgrade($dng); 6086 6087 # all done, $x contains now the result 6088 $x; 6089} 6090 6091sub _pow { 6092 # Calculate a power where $y is a non-integer, like 2 ** 0.3 6093 my ($x, $y, @r) = @_; 6094 my $class = ref($x); 6095 6096 # if $y == 0.5, it is sqrt($x) 6097 $HALF = $class->new($HALF) unless ref($HALF); 6098 return $x->bsqrt(@r, $y) if $y->bcmp($HALF) == 0; 6099 6100 # Using: 6101 # a ** x == e ** (x * ln a) 6102 6103 # u = y * ln x 6104 # _ _ 6105 # Taylor: | u u^2 u^3 | 6106 # x ** y = 1 + | --- + --- + ----- + ... | 6107 # |_ 1 1*2 1*2*3 _| 6108 6109 # we need to limit the accuracy to protect against overflow 6110 my $fallback = 0; 6111 my ($scale, @params); 6112 ($x, @params) = $x->_find_round_parameters(@r); 6113 6114 return $x if $x->is_nan(); # error in _find_round_parameters? 6115 6116 # no rounding at all, so must use fallback 6117 if (scalar @params == 0) { 6118 # simulate old behaviour 6119 $params[0] = $class->div_scale(); # and round to it as accuracy 6120 $params[1] = undef; # disable P 6121 $scale = $params[0]+4; # at least four more for proper round 6122 $params[2] = $r[2]; # round mode by caller or undef 6123 $fallback = 1; # to clear a/p afterwards 6124 } else { 6125 # the 4 below is empirical, and there might be cases where it is not 6126 # enough... 6127 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined 6128 } 6129 6130 # When user set globals, they would interfere with our calculation, so 6131 # disable them and later re-enable them. 6132 6133 my $ab = $class -> accuracy(); 6134 my $pb = $class -> precision(); 6135 $class -> accuracy(undef); 6136 $class -> precision(undef); 6137 6138 # Disabling upgrading and downgrading is no longer necessary to avoid an 6139 # infinite recursion, but it avoids unnecessary upgrading and downgrading in 6140 # the intermediate computations. 6141 6142 my $upg = $class -> upgrade(); 6143 my $dng = $class -> downgrade(); 6144 $class -> upgrade(undef); 6145 $class -> downgrade(undef); 6146 6147 # We also need to disable any set A or P on $x (_find_round_parameters took 6148 # them already into account), since these would interfere, too. 6149 6150 $x->{accuracy} = undef; 6151 $x->{precision} = undef; 6152 6153 my ($limit, $v, $u, $below, $factor, $next, $over); 6154 6155 $u = $x->copy()->blog(undef, $scale)->bmul($y); 6156 my $do_invert = ($u->{sign} eq '-'); 6157 $u = $u->bneg() if $do_invert; 6158 $v = $class->bone(); # 1 6159 $factor = $class->new(2); # 2 6160 $x = $x->bone(); # first term: 1 6161 6162 $below = $v->copy(); 6163 $over = $u->copy(); 6164 6165 $limit = $class->new("1E-". ($scale-1)); 6166 while (3 < 5) { 6167 # we calculate the next term, and add it to the last 6168 # when the next term is below our limit, it won't affect the outcome 6169 # anymore, so we stop: 6170 $next = $over->copy()->bdiv($below, $scale); 6171 last if $next->bacmp($limit) <= 0; 6172 $x = $x->badd($next); 6173 # calculate things for the next term 6174 $over *= $u; 6175 $below *= $factor; 6176 $factor = $factor->binc(); 6177 6178 last if $x->{sign} !~ /^[-+]$/; 6179 } 6180 6181 if ($do_invert) { 6182 my $x_copy = $x->copy(); 6183 $x = $x->bone->bdiv($x_copy, $scale); 6184 } 6185 6186 # shortcut to not run through _find_round_parameters again 6187 if (defined $params[0]) { 6188 $x = $x->bround($params[0], $params[2]); # then round accordingly 6189 } else { 6190 $x = $x->bfround($params[1], $params[2]); # then round accordingly 6191 } 6192 if ($fallback) { 6193 # clear a/p after round, since user did not request it 6194 $x->{accuracy} = undef; 6195 $x->{precision} = undef; 6196 } 6197 6198 # Restore globals. We need to do it like this, because setting one 6199 # undefines the other. 6200 6201 if (defined $ab) { 6202 $class -> accuracy($ab); 6203 } else { 6204 $class -> precision($pb); 6205 } 6206 6207 $class -> upgrade($upg); 6208 $class -> downgrade($dng); 6209 6210 $x; 6211} 6212 6213# These functions are only provided for backwards compabibility so that old 6214# version of Math::BigRat etc. don't complain about missing them. 6215 6216sub _e_add { 6217 my ($x, $y, $xs, $ys) = @_; 6218 return $LIB -> _sadd($x, $xs, $y, $ys); 6219} 6220 6221sub _e_sub { 6222 my ($x, $y, $xs, $ys) = @_; 6223 return $LIB -> _ssub($x, $xs, $y, $ys); 6224} 6225 62261; 6227 6228__END__ 6229 6230=pod 6231 6232=head1 NAME 6233 6234Math::BigFloat - arbitrary size floating point math package 6235 6236=head1 SYNOPSIS 6237 6238 use Math::BigFloat; 6239 6240 # Configuration methods (may be used as class methods and instance methods) 6241 6242 Math::BigFloat->accuracy(); # get class accuracy 6243 Math::BigFloat->accuracy($n); # set class accuracy 6244 Math::BigFloat->precision(); # get class precision 6245 Math::BigFloat->precision($n); # set class precision 6246 Math::BigFloat->round_mode(); # get class rounding mode 6247 Math::BigFloat->round_mode($m); # set global round mode, must be one of 6248 # 'even', 'odd', '+inf', '-inf', 'zero', 6249 # 'trunc', or 'common' 6250 Math::BigFloat->config("lib"); # name of backend math library 6251 6252 # Constructor methods (when the class methods below are used as instance 6253 # methods, the value is assigned the invocand) 6254 6255 $x = Math::BigFloat->new($str); # defaults to 0 6256 $x = Math::BigFloat->new('0x123'); # from hexadecimal 6257 $x = Math::BigFloat->new('0o377'); # from octal 6258 $x = Math::BigFloat->new('0b101'); # from binary 6259 $x = Math::BigFloat->from_hex('0xc.afep+3'); # from hex 6260 $x = Math::BigFloat->from_hex('cafe'); # ditto 6261 $x = Math::BigFloat->from_oct('1.3267p-4'); # from octal 6262 $x = Math::BigFloat->from_oct('01.3267p-4'); # ditto 6263 $x = Math::BigFloat->from_oct('0o1.3267p-4'); # ditto 6264 $x = Math::BigFloat->from_oct('0377'); # ditto 6265 $x = Math::BigFloat->from_bin('0b1.1001p-4'); # from binary 6266 $x = Math::BigFloat->from_bin('0101'); # ditto 6267 $x = Math::BigFloat->from_ieee754($b, "binary64"); # from IEEE-754 bytes 6268 $x = Math::BigFloat->bzero(); # create a +0 6269 $x = Math::BigFloat->bone(); # create a +1 6270 $x = Math::BigFloat->bone('-'); # create a -1 6271 $x = Math::BigFloat->binf(); # create a +inf 6272 $x = Math::BigFloat->binf('-'); # create a -inf 6273 $x = Math::BigFloat->bnan(); # create a Not-A-Number 6274 $x = Math::BigFloat->bpi(); # returns pi 6275 6276 $y = $x->copy(); # make a copy (unlike $y = $x) 6277 $y = $x->as_int(); # return as BigInt 6278 $y = $x->as_float(); # return as a Math::BigFloat 6279 $y = $x->as_rat(); # return as a Math::BigRat 6280 6281 # Boolean methods (these don't modify the invocand) 6282 6283 $x->is_zero(); # if $x is 0 6284 $x->is_one(); # if $x is +1 6285 $x->is_one("+"); # ditto 6286 $x->is_one("-"); # if $x is -1 6287 $x->is_inf(); # if $x is +inf or -inf 6288 $x->is_inf("+"); # if $x is +inf 6289 $x->is_inf("-"); # if $x is -inf 6290 $x->is_nan(); # if $x is NaN 6291 6292 $x->is_positive(); # if $x > 0 6293 $x->is_pos(); # ditto 6294 $x->is_negative(); # if $x < 0 6295 $x->is_neg(); # ditto 6296 6297 $x->is_odd(); # if $x is odd 6298 $x->is_even(); # if $x is even 6299 $x->is_int(); # if $x is an integer 6300 6301 # Comparison methods 6302 6303 $x->bcmp($y); # compare numbers (undef, < 0, == 0, > 0) 6304 $x->bacmp($y); # compare absolutely (undef, < 0, == 0, > 0) 6305 $x->beq($y); # true if and only if $x == $y 6306 $x->bne($y); # true if and only if $x != $y 6307 $x->blt($y); # true if and only if $x < $y 6308 $x->ble($y); # true if and only if $x <= $y 6309 $x->bgt($y); # true if and only if $x > $y 6310 $x->bge($y); # true if and only if $x >= $y 6311 6312 # Arithmetic methods 6313 6314 $x->bneg(); # negation 6315 $x->babs(); # absolute value 6316 $x->bsgn(); # sign function (-1, 0, 1, or NaN) 6317 $x->bnorm(); # normalize (no-op) 6318 $x->binc(); # increment $x by 1 6319 $x->bdec(); # decrement $x by 1 6320 $x->badd($y); # addition (add $y to $x) 6321 $x->bsub($y); # subtraction (subtract $y from $x) 6322 $x->bmul($y); # multiplication (multiply $x by $y) 6323 $x->bmuladd($y,$z); # $x = $x * $y + $z 6324 $x->bdiv($y); # division (floored), set $x to quotient 6325 # return (quo,rem) or quo if scalar 6326 $x->btdiv($y); # division (truncated), set $x to quotient 6327 # return (quo,rem) or quo if scalar 6328 $x->bmod($y); # modulus (x % y) 6329 $x->btmod($y); # modulus (truncated) 6330 $x->bmodinv($mod); # modular multiplicative inverse 6331 $x->bmodpow($y,$mod); # modular exponentiation (($x ** $y) % $mod) 6332 $x->bpow($y); # power of arguments (x ** y) 6333 $x->blog(); # logarithm of $x to base e (Euler's number) 6334 $x->blog($base); # logarithm of $x to base $base (e.g., base 2) 6335 $x->bexp(); # calculate e ** $x where e is Euler's number 6336 $x->bnok($y); # x over y (binomial coefficient n over k) 6337 $x->bsin(); # sine 6338 $x->bcos(); # cosine 6339 $x->batan(); # inverse tangent 6340 $x->batan2($y); # two-argument inverse tangent 6341 $x->bsqrt(); # calculate square root 6342 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root) 6343 $x->bfac(); # factorial of $x (1*2*3*4*..$x) 6344 6345 $x->blsft($n); # left shift $n places in base 2 6346 $x->blsft($n,$b); # left shift $n places in base $b 6347 # returns (quo,rem) or quo (scalar context) 6348 $x->brsft($n); # right shift $n places in base 2 6349 $x->brsft($n,$b); # right shift $n places in base $b 6350 # returns (quo,rem) or quo (scalar context) 6351 6352 # Bitwise methods 6353 6354 $x->bblsft($y); # bitwise left shift 6355 $x->bbrsft($y); # bitwise right shift 6356 $x->band($y); # bitwise and 6357 $x->bior($y); # bitwise inclusive or 6358 $x->bxor($y); # bitwise exclusive or 6359 $x->bnot(); # bitwise not (two's complement) 6360 6361 # Rounding methods 6362 $x->round($A,$P,$mode); # round to accuracy or precision using 6363 # rounding mode $mode 6364 $x->bround($n); # accuracy: preserve $n digits 6365 $x->bfround($n); # $n > 0: round to $nth digit left of dec. point 6366 # $n < 0: round to $nth digit right of dec. point 6367 $x->bfloor(); # round towards minus infinity 6368 $x->bceil(); # round towards plus infinity 6369 $x->bint(); # round towards zero 6370 6371 # Other mathematical methods 6372 6373 $x->bgcd($y); # greatest common divisor 6374 $x->blcm($y); # least common multiple 6375 6376 # Object property methods (do not modify the invocand) 6377 6378 $x->sign(); # the sign, either +, - or NaN 6379 $x->digit($n); # the nth digit, counting from the right 6380 $x->digit(-$n); # the nth digit, counting from the left 6381 $x->length(); # return number of digits in number 6382 ($xl,$f) = $x->length(); # length of number and length of fraction 6383 # part, latter is always 0 digits long 6384 # for Math::BigInt objects 6385 $x->mantissa(); # return (signed) mantissa as BigInt 6386 $x->exponent(); # return exponent as BigInt 6387 $x->parts(); # return (mantissa,exponent) as BigInt 6388 $x->sparts(); # mantissa and exponent (as integers) 6389 $x->nparts(); # mantissa and exponent (normalised) 6390 $x->eparts(); # mantissa and exponent (engineering notation) 6391 $x->dparts(); # integer and fraction part 6392 $x->fparts(); # numerator and denominator 6393 $x->numerator(); # numerator 6394 $x->denominator(); # denominator 6395 6396 # Conversion methods (do not modify the invocand) 6397 6398 $x->bstr(); # decimal notation, possibly zero padded 6399 $x->bsstr(); # string in scientific notation with integers 6400 $x->bnstr(); # string in normalized notation 6401 $x->bestr(); # string in engineering notation 6402 $x->bdstr(); # string in decimal notation 6403 $x->bfstr(); # string in fractional notation 6404 6405 $x->as_hex(); # as signed hexadecimal string with prefixed 0x 6406 $x->as_bin(); # as signed binary string with prefixed 0b 6407 $x->as_oct(); # as signed octal string with prefixed 0 6408 $x->to_ieee754($format); # to bytes encoded according to IEEE 754-2008 6409 6410 # Other conversion methods 6411 6412 $x->numify(); # return as scalar (might overflow or underflow) 6413 6414=head1 DESCRIPTION 6415 6416Math::BigFloat provides support for arbitrary precision floating point. 6417Overloading is also provided for Perl operators. 6418 6419All operators (including basic math operations) are overloaded if you 6420declare your big floating point numbers as 6421 6422 $x = Math::BigFloat -> new('12_3.456_789_123_456_789E-2'); 6423 6424Operations with overloaded operators preserve the arguments, which is 6425exactly what you expect. 6426 6427=head2 Input 6428 6429Input values to these routines may be any scalar number or string that looks 6430like a number. Anything that is accepted by Perl as a literal numeric constant 6431should be accepted by this module. 6432 6433=over 6434 6435=item * 6436 6437Leading and trailing whitespace is ignored. 6438 6439=item * 6440 6441Leading zeros are ignored, except for floating point numbers with a binary 6442exponent, in which case the number is interpreted as an octal floating point 6443number. For example, "01.4p+0" gives 1.5, "00.4p+0" gives 0.5, but "0.4p+0" 6444gives a NaN. And while "0377" gives 255, "0377p0" gives 255. 6445 6446=item * 6447 6448If the string has a "0x" or "0X" prefix, it is interpreted as a hexadecimal 6449number. 6450 6451=item * 6452 6453If the string has a "0o" or "0O" prefix, it is interpreted as an octal number. A 6454floating point literal with a "0" prefix is also interpreted as an octal number. 6455 6456=item * 6457 6458If the string has a "0b" or "0B" prefix, it is interpreted as a binary number. 6459 6460=item * 6461 6462Underline characters are allowed in the same way as they are allowed in literal 6463numerical constants. 6464 6465=item * 6466 6467If the string can not be interpreted, NaN is returned. 6468 6469=item * 6470 6471For hexadecimal, octal, and binary floating point numbers, the exponent must be 6472separated from the significand (mantissa) by the letter "p" or "P", not "e" or 6473"E" as with decimal numbers. 6474 6475=back 6476 6477Some examples of valid string input 6478 6479 Input string Resulting value 6480 6481 123 123 6482 1.23e2 123 6483 12300e-2 123 6484 6485 67_538_754 67538754 6486 -4_5_6.7_8_9e+0_1_0 -4567890000000 6487 6488 0x13a 314 6489 0x13ap0 314 6490 0x1.3ap+8 314 6491 0x0.00013ap+24 314 6492 0x13a000p-12 314 6493 6494 0o472 314 6495 0o1.164p+8 314 6496 0o0.0001164p+20 314 6497 0o1164000p-10 314 6498 6499 0472 472 Note! 6500 01.164p+8 314 6501 00.0001164p+20 314 6502 01164000p-10 314 6503 6504 0b100111010 314 6505 0b1.0011101p+8 314 6506 0b0.00010011101p+12 314 6507 0b100111010000p-3 314 6508 6509 0x1.921fb5p+1 3.14159262180328369140625e+0 6510 0o1.2677025p1 2.71828174591064453125 6511 01.2677025p1 2.71828174591064453125 6512 0b1.1001p-4 9.765625e-2 6513 6514=head2 Output 6515 6516Output values are usually Math::BigFloat objects. 6517 6518Boolean operators C<is_zero()>, C<is_one()>, C<is_inf()>, etc. return true or 6519false. 6520 6521Comparison operators C<bcmp()> and C<bacmp()>) return -1, 0, 1, or 6522undef. 6523 6524=head1 METHODS 6525 6526Math::BigFloat supports all methods that Math::BigInt supports, except it 6527calculates non-integer results when possible. Please see L<Math::BigInt> for a 6528full description of each method. Below are just the most important differences: 6529 6530=head2 Configuration methods 6531 6532=over 6533 6534=item accuracy() 6535 6536 $x->accuracy(5); # local for $x 6537 CLASS->accuracy(5); # global for all members of CLASS 6538 # Note: This also applies to new()! 6539 6540 $A = $x->accuracy(); # read out accuracy that affects $x 6541 $A = CLASS->accuracy(); # read out global accuracy 6542 6543Set or get the global or local accuracy, aka how many significant digits the 6544results have. If you set a global accuracy, then this also applies to new()! 6545 6546Warning! The accuracy I<sticks>, e.g. once you created a number under the 6547influence of C<< CLASS->accuracy($A) >>, all results from math operations with 6548that number will also be rounded. 6549 6550In most cases, you should probably round the results explicitly using one of 6551L<Math::BigInt/round()>, L<Math::BigInt/bround()> or L<Math::BigInt/bfround()> 6552or by passing the desired accuracy to the math operation as additional 6553parameter: 6554 6555 my $x = Math::BigInt->new(30000); 6556 my $y = Math::BigInt->new(7); 6557 print scalar $x->copy()->bdiv($y, 2); # print 4300 6558 print scalar $x->copy()->bdiv($y)->bround(2); # print 4300 6559 6560=item precision() 6561 6562 $x->precision(-2); # local for $x, round at the second 6563 # digit right of the dot 6564 $x->precision(2); # ditto, round at the second digit 6565 # left of the dot 6566 6567 CLASS->precision(5); # Global for all members of CLASS 6568 # This also applies to new()! 6569 CLASS->precision(-5); # ditto 6570 6571 $P = CLASS->precision(); # read out global precision 6572 $P = $x->precision(); # read out precision that affects $x 6573 6574Note: You probably want to use L</accuracy()> instead. With L</accuracy()> you 6575set the number of digits each result should have, with L</precision()> you 6576set the place where to round! 6577 6578=back 6579 6580=head2 Constructor methods 6581 6582=over 6583 6584=item from_hex() 6585 6586 $x -> from_hex("0x1.921fb54442d18p+1"); 6587 $x = Math::BigFloat -> from_hex("0x1.921fb54442d18p+1"); 6588 6589Interpret input as a hexadecimal string.A prefix ("0x", "x", ignoring case) is 6590optional. A single underscore character ("_") may be placed between any two 6591digits. If the input is invalid, a NaN is returned. The exponent is in base 2 6592using decimal digits. 6593 6594If called as an instance method, the value is assigned to the invocand. 6595 6596=item from_oct() 6597 6598 $x -> from_oct("1.3267p-4"); 6599 $x = Math::BigFloat -> from_oct("1.3267p-4"); 6600 6601Interpret input as an octal string. A single underscore character ("_") may be 6602placed between any two digits. If the input is invalid, a NaN is returned. The 6603exponent is in base 2 using decimal digits. 6604 6605If called as an instance method, the value is assigned to the invocand. 6606 6607=item from_bin() 6608 6609 $x -> from_bin("0b1.1001p-4"); 6610 $x = Math::BigFloat -> from_bin("0b1.1001p-4"); 6611 6612Interpret input as a hexadecimal string. A prefix ("0b" or "b", ignoring case) 6613is optional. A single underscore character ("_") may be placed between any two 6614digits. If the input is invalid, a NaN is returned. The exponent is in base 2 6615using decimal digits. 6616 6617If called as an instance method, the value is assigned to the invocand. 6618 6619=item from_ieee754() 6620 6621Interpret the input as a value encoded as described in IEEE754-2008. The input 6622can be given as a byte string, hex string or binary string. The input is 6623assumed to be in big-endian byte-order. 6624 6625 # both $dbl and $mbf are 3.141592... 6626 $bytes = "\x40\x09\x21\xfb\x54\x44\x2d\x18"; 6627 $dbl = unpack "d>", $bytes; 6628 $mbf = Math::BigFloat -> from_ieee754($bytes, "binary64"); 6629 6630=item bpi() 6631 6632 print Math::BigFloat->bpi(100), "\n"; 6633 6634Calculate PI to N digits (including the 3 before the dot). The result is 6635rounded according to the current rounding mode, which defaults to "even". 6636 6637This method was added in v1.87 of Math::BigInt (June 2007). 6638 6639=back 6640 6641=head2 Arithmetic methods 6642 6643=over 6644 6645=item bmuladd() 6646 6647 $x->bmuladd($y,$z); 6648 6649Multiply $x by $y, and then add $z to the result. 6650 6651This method was added in v1.87 of Math::BigInt (June 2007). 6652 6653=item binv() 6654 6655 $x->binv(); 6656 6657Invert the value of $x, i.e., compute 1/$x. 6658 6659=item bdiv() 6660 6661 $q = $x->bdiv($y); 6662 ($q, $r) = $x->bdiv($y); 6663 6664In scalar context, divides $x by $y and returns the result to the given or 6665default accuracy/precision. In list context, does floored division 6666(F-division), returning an integer $q and a remainder $r so that $x = $q * $y + 6667$r. The remainer (modulo) is equal to what is returned by C<< $x->bmod($y) >>. 6668 6669=item bmod() 6670 6671 $x->bmod($y); 6672 6673Returns $x modulo $y. When $x is finite, and $y is finite and non-zero, the 6674result is identical to the remainder after floored division (F-division). If, 6675in addition, both $x and $y are integers, the result is identical to the result 6676from Perl's % operator. 6677 6678=item bexp() 6679 6680 $x->bexp($accuracy); # calculate e ** X 6681 6682Calculates the expression C<e ** $x> where C<e> is Euler's number. 6683 6684This method was added in v1.82 of Math::BigInt (April 2007). 6685 6686=item bnok() 6687 6688 $x->bnok($y); # x over y (binomial coefficient n over k) 6689 6690Calculates the binomial coefficient n over k, also called the "choose" 6691function. The result is equivalent to: 6692 6693 ( n ) n! 6694 | - | = ------- 6695 ( k ) k!(n-k)! 6696 6697This method was added in v1.84 of Math::BigInt (April 2007). 6698 6699=item bsin() 6700 6701 my $x = Math::BigFloat->new(1); 6702 print $x->bsin(100), "\n"; 6703 6704Calculate the sinus of $x, modifying $x in place. 6705 6706This method was added in v1.87 of Math::BigInt (June 2007). 6707 6708=item bcos() 6709 6710 my $x = Math::BigFloat->new(1); 6711 print $x->bcos(100), "\n"; 6712 6713Calculate the cosinus of $x, modifying $x in place. 6714 6715This method was added in v1.87 of Math::BigInt (June 2007). 6716 6717=item batan() 6718 6719 my $x = Math::BigFloat->new(1); 6720 print $x->batan(100), "\n"; 6721 6722Calculate the arcus tanges of $x, modifying $x in place. See also L</batan2()>. 6723 6724This method was added in v1.87 of Math::BigInt (June 2007). 6725 6726=item batan2() 6727 6728 my $y = Math::BigFloat->new(2); 6729 my $x = Math::BigFloat->new(3); 6730 print $y->batan2($x), "\n"; 6731 6732Calculate the arcus tanges of C<$y> divided by C<$x>, modifying $y in place. 6733See also L</batan()>. 6734 6735This method was added in v1.87 of Math::BigInt (June 2007). 6736 6737=item as_float() 6738 6739This method is called when Math::BigFloat encounters an object it doesn't know 6740how to handle. For instance, assume $x is a Math::BigFloat, or subclass 6741thereof, and $y is defined, but not a Math::BigFloat, or subclass thereof. If 6742you do 6743 6744 $x -> badd($y); 6745 6746$y needs to be converted into an object that $x can deal with. This is done by 6747first checking if $y is something that $x might be upgraded to. If that is the 6748case, no further attempts are made. The next is to see if $y supports the 6749method C<as_float()>. The method C<as_float()> is expected to return either an 6750object that has the same class as $x, a subclass thereof, or a string that 6751C<ref($x)-E<gt>new()> can parse to create an object. 6752 6753In Math::BigFloat, C<as_float()> has the same effect as C<copy()>. 6754 6755=item to_ieee754() 6756 6757Encodes the invocand as a byte string in the given format as specified in IEEE 6758754-2008. Note that the encoded value is the nearest possible representation of 6759the value. This value might not be exactly the same as the value in the 6760invocand. 6761 6762 # $x = 3.1415926535897932385 6763 $x = Math::BigFloat -> bpi(30); 6764 6765 $b = $x -> to_ieee754("binary64"); # encode as 8 bytes 6766 $h = unpack "H*", $b; # "400921fb54442d18" 6767 6768 # 3.141592653589793115997963... 6769 $y = Math::BigFloat -> from_ieee754($h, "binary64"); 6770 6771All binary formats in IEEE 754-2008 are accepted. For convenience, som aliases 6772are recognized: "half" for "binary16", "single" for "binary32", "double" for 6773"binary64", "quadruple" for "binary128", "octuple" for "binary256", and 6774"sexdecuple" for "binary512". 6775 6776See also L<https://en.wikipedia.org/wiki/IEEE_754>. 6777 6778=back 6779 6780=head2 ACCURACY AND PRECISION 6781 6782See also: L<Rounding|/Rounding>. 6783 6784Math::BigFloat supports both precision (rounding to a certain place before or 6785after the dot) and accuracy (rounding to a certain number of digits). For a 6786full documentation, examples and tips on these topics please see the large 6787section about rounding in L<Math::BigInt>. 6788 6789Since things like C<sqrt(2)> or C<1 / 3> must presented with a limited 6790accuracy lest a operation consumes all resources, each operation produces 6791no more than the requested number of digits. 6792 6793If there is no global precision or accuracy set, B<and> the operation in 6794question was not called with a requested precision or accuracy, B<and> the 6795input $x has no accuracy or precision set, then a fallback parameter will 6796be used. For historical reasons, it is called C<div_scale> and can be accessed 6797via: 6798 6799 $d = Math::BigFloat->div_scale(); # query 6800 Math::BigFloat->div_scale($n); # set to $n digits 6801 6802The default value for C<div_scale> is 40. 6803 6804In case the result of one operation has more digits than specified, 6805it is rounded. The rounding mode taken is either the default mode, or the one 6806supplied to the operation after the I<scale>: 6807 6808 $x = Math::BigFloat->new(2); 6809 Math::BigFloat->accuracy(5); # 5 digits max 6810 $y = $x->copy()->bdiv(3); # gives 0.66667 6811 $y = $x->copy()->bdiv(3,6); # gives 0.666667 6812 $y = $x->copy()->bdiv(3,6,undef,'odd'); # gives 0.666667 6813 Math::BigFloat->round_mode('zero'); 6814 $y = $x->copy()->bdiv(3,6); # will also give 0.666667 6815 6816Note that C<< Math::BigFloat->accuracy() >> and 6817C<< Math::BigFloat->precision() >> set the global variables, and thus B<any> 6818newly created number will be subject to the global rounding B<immediately>. This 6819means that in the examples above, the C<3> as argument to C<bdiv()> will also 6820get an accuracy of B<5>. 6821 6822It is less confusing to either calculate the result fully, and afterwards 6823round it explicitly, or use the additional parameters to the math 6824functions like so: 6825 6826 use Math::BigFloat; 6827 $x = Math::BigFloat->new(2); 6828 $y = $x->copy()->bdiv(3); 6829 print $y->bround(5),"\n"; # gives 0.66667 6830 6831 or 6832 6833 use Math::BigFloat; 6834 $x = Math::BigFloat->new(2); 6835 $y = $x->copy()->bdiv(3,5); # gives 0.66667 6836 print "$y\n"; 6837 6838=head2 Rounding 6839 6840=over 6841 6842=item bfround ( +$scale ) 6843 6844Rounds to the $scale'th place left from the '.', counting from the dot. 6845The first digit is numbered 1. 6846 6847=item bfround ( -$scale ) 6848 6849Rounds to the $scale'th place right from the '.', counting from the dot. 6850 6851=item bfround ( 0 ) 6852 6853Rounds to an integer. 6854 6855=item bround ( +$scale ) 6856 6857Preserves accuracy to $scale digits from the left (aka significant digits) and 6858pads the rest with zeros. If the number is between 1 and -1, the significant 6859digits count from the first non-zero after the '.' 6860 6861=item bround ( -$scale ) and bround ( 0 ) 6862 6863These are effectively no-ops. 6864 6865=back 6866 6867All rounding functions take as a second parameter a rounding mode from one of 6868the following: 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'. 6869 6870The default rounding mode is 'even'. By using 6871C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default 6872mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is 6873no longer supported. 6874The second parameter to the round functions then overrides the default 6875temporarily. 6876 6877The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses 6878'trunc' as rounding mode to make it equivalent to: 6879 6880 $x = 2.5; 6881 $y = int($x) + 2; 6882 6883You can override this by passing the desired rounding mode as parameter to 6884C<as_number()>: 6885 6886 $x = Math::BigFloat->new(2.5); 6887 $y = $x->as_number('odd'); # $y = 3 6888 6889=head1 NUMERIC LITERALS 6890 6891After C<use Math::BigFloat ':constant'> all numeric literals in the given scope 6892are converted to C<Math::BigFloat> objects. This conversion happens at compile 6893time. 6894 6895For example, 6896 6897 perl -MMath::BigFloat=:constant -le 'print 2e-150' 6898 6899prints the exact value of C<2e-150>. Note that without conversion of constants 6900the expression C<2e-150> is calculated using Perl scalars, which leads to an 6901inaccuracte result. 6902 6903Note that strings are not affected, so that 6904 6905 use Math::BigFloat qw/:constant/; 6906 6907 $y = "1234567890123456789012345678901234567890" 6908 + "123456789123456789"; 6909 6910does not give you what you expect. You need an explicit Math::BigFloat->new() 6911around at least one of the operands. You should also quote large constants to 6912prevent loss of precision: 6913 6914 use Math::BigFloat; 6915 6916 $x = Math::BigFloat->new("1234567889123456789123456789123456789"); 6917 6918Without the quotes Perl converts the large number to a floating point constant 6919at compile time, and then converts the result to a Math::BigFloat object at 6920runtime, which results in an inaccurate result. 6921 6922=head2 Hexadecimal, octal, and binary floating point literals 6923 6924Perl (and this module) accepts hexadecimal, octal, and binary floating point 6925literals, but use them with care with Perl versions before v5.32.0, because some 6926versions of Perl silently give the wrong result. Below are some examples of 6927different ways to write the number decimal 314. 6928 6929Hexadecimal floating point literals: 6930 6931 0x1.3ap+8 0X1.3AP+8 6932 0x1.3ap8 0X1.3AP8 6933 0x13a0p-4 0X13A0P-4 6934 6935Octal floating point literals (with "0" prefix): 6936 6937 01.164p+8 01.164P+8 6938 01.164p8 01.164P8 6939 011640p-4 011640P-4 6940 6941Octal floating point literals (with "0o" prefix) (requires v5.34.0): 6942 6943 0o1.164p+8 0O1.164P+8 6944 0o1.164p8 0O1.164P8 6945 0o11640p-4 0O11640P-4 6946 6947Binary floating point literals: 6948 6949 0b1.0011101p+8 0B1.0011101P+8 6950 0b1.0011101p8 0B1.0011101P8 6951 0b10011101000p-2 0B10011101000P-2 6952 6953=head2 Math library 6954 6955Math with the numbers is done (by default) by a module called 6956Math::BigInt::Calc. This is equivalent to saying: 6957 6958 use Math::BigFloat lib => "Calc"; 6959 6960You can change this by using: 6961 6962 use Math::BigFloat lib => "GMP"; 6963 6964B<Note>: General purpose packages should not be explicit about the library to 6965use; let the script author decide which is best. 6966 6967Note: The keyword 'lib' will warn when the requested library could not be 6968loaded. To suppress the warning use 'try' instead: 6969 6970 use Math::BigFloat try => "GMP"; 6971 6972If your script works with huge numbers and Calc is too slow for them, you can 6973also for the loading of one of these libraries and if none of them can be used, 6974the code will die: 6975 6976 use Math::BigFloat only => "GMP,Pari"; 6977 6978The following would first try to find Math::BigInt::Foo, then Math::BigInt::Bar, 6979and when this also fails, revert to Math::BigInt::Calc: 6980 6981 use Math::BigFloat lib => "Foo,Math::BigInt::Bar"; 6982 6983See the respective low-level library documentation for further details. 6984 6985See L<Math::BigInt> for more details about using a different low-level library. 6986 6987=head1 EXPORTS 6988 6989C<Math::BigFloat> exports nothing by default, but can export the C<bpi()> 6990method: 6991 6992 use Math::BigFloat qw/bpi/; 6993 6994 print bpi(10), "\n"; 6995 6996=over 6997 6998=item stringify, bstr() 6999 7000Both stringify and bstr() now drop the leading '+'. The old code would return 7001'+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for 7002reasoning and details. 7003 7004=item brsft() 7005 7006The following will probably not print what you expect: 7007 7008 my $c = Math::BigFloat->new('3.14159'); 7009 print $c->brsft(3,10),"\n"; # prints 0.00314153.1415 7010 7011It prints both quotient and remainder, since print calls C<brsft()> in list 7012context. Also, C<< $c->brsft() >> will modify $c, so be careful. 7013You probably want to use 7014 7015 print scalar $c->copy()->brsft(3,10),"\n"; 7016 # or if you really want to modify $c 7017 print scalar $c->brsft(3,10),"\n"; 7018 7019instead. 7020 7021=item Modifying and = 7022 7023Beware of: 7024 7025 $x = Math::BigFloat->new(5); 7026 $y = $x; 7027 7028It will not do what you think, e.g. making a copy of $x. Instead it just makes 7029a second reference to the B<same> object and stores it in $y. Thus anything 7030that modifies $x will modify $y (except overloaded math operators), and vice 7031versa. See L<Math::BigInt> for details and how to avoid that. 7032 7033=item precision() vs. accuracy() 7034 7035A common pitfall is to use L</precision()> when you want to round a result to 7036a certain number of digits: 7037 7038 use Math::BigFloat; 7039 7040 Math::BigFloat->precision(4); # does not do what you 7041 # think it does 7042 my $x = Math::BigFloat->new(12345); # rounds $x to "12000"! 7043 print "$x\n"; # print "12000" 7044 my $y = Math::BigFloat->new(3); # rounds $y to "0"! 7045 print "$y\n"; # print "0" 7046 $z = $x / $y; # 12000 / 0 => NaN! 7047 print "$z\n"; 7048 print $z->precision(),"\n"; # 4 7049 7050Replacing L</precision()> with L</accuracy()> is probably not what you want, 7051either: 7052 7053 use Math::BigFloat; 7054 7055 Math::BigFloat->accuracy(4); # enables global rounding: 7056 my $x = Math::BigFloat->new(123456); # rounded immediately 7057 # to "12350" 7058 print "$x\n"; # print "123500" 7059 my $y = Math::BigFloat->new(3); # rounded to "3 7060 print "$y\n"; # print "3" 7061 print $z = $x->copy()->bdiv($y),"\n"; # 41170 7062 print $z->accuracy(),"\n"; # 4 7063 7064What you want to use instead is: 7065 7066 use Math::BigFloat; 7067 7068 my $x = Math::BigFloat->new(123456); # no rounding 7069 print "$x\n"; # print "123456" 7070 my $y = Math::BigFloat->new(3); # no rounding 7071 print "$y\n"; # print "3" 7072 print $z = $x->copy()->bdiv($y,4),"\n"; # 41150 7073 print $z->accuracy(),"\n"; # undef 7074 7075In addition to computing what you expected, the last example also does B<not> 7076"taint" the result with an accuracy or precision setting, which would 7077influence any further operation. 7078 7079=back 7080 7081=head1 BUGS 7082 7083Please report any bugs or feature requests to 7084C<bug-math-bigint at rt.cpan.org>, or through the web interface at 7085L<https://rt.cpan.org/Ticket/Create.html?Queue=Math-BigInt> (requires login). 7086We will be notified, and then you'll automatically be notified of progress on 7087your bug as I make changes. 7088 7089=head1 SUPPORT 7090 7091You can find documentation for this module with the perldoc command. 7092 7093 perldoc Math::BigFloat 7094 7095You can also look for information at: 7096 7097=over 4 7098 7099=item * GitHub 7100 7101L<https://github.com/pjacklam/p5-Math-BigInt> 7102 7103=item * RT: CPAN's request tracker 7104 7105L<https://rt.cpan.org/Dist/Display.html?Name=Math-BigInt> 7106 7107=item * MetaCPAN 7108 7109L<https://metacpan.org/release/Math-BigInt> 7110 7111=item * CPAN Testers Matrix 7112 7113L<http://matrix.cpantesters.org/?dist=Math-BigInt> 7114 7115=back 7116 7117=head1 LICENSE 7118 7119This program is free software; you may redistribute it and/or modify it under 7120the same terms as Perl itself. 7121 7122=head1 SEE ALSO 7123 7124L<Math::BigInt> and L<Math::BigRat> as well as the backend libraries 7125L<Math::BigInt::FastCalc>, L<Math::BigInt::GMP>, and L<Math::BigInt::Pari>, 7126L<Math::BigInt::GMPz>, and L<Math::BigInt::BitVect>. 7127 7128The pragmas L<bigint>, L<bigfloat>, and L<bigrat> might also be of interest. In 7129addition there is the L<bignum> pragma which does upgrading and downgrading. 7130 7131=head1 AUTHORS 7132 7133=over 4 7134 7135=item * 7136 7137Mark Biggar, overloaded interface by Ilya Zakharevich, 1996-2001. 7138 7139=item * 7140 7141Completely rewritten by Tels L<http://bloodgate.com> in 2001-2008. 7142 7143=item * 7144 7145Florian Ragwitz E<lt>flora@cpan.orgE<gt>, 2010. 7146 7147=item * 7148 7149Peter John Acklam E<lt>pjacklam@gmail.comE<gt>, 2011-. 7150 7151=back 7152 7153=cut 7154