cpan/Math-BigInt/lib/Math/BigFloat.pm

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