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