[Summary view] [Print] [Text view]

   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 internally used:
   8  #   _e    : exponent (ref to $CALC object)
   9  #   _m    : mantissa (ref to $CALC object)
  10  #   _es    : sign of _e
  11  # sign    : +,-,+inf,-inf, or "NaN" if not a number
  12  #   _a    : accuracy
  13  #   _p    : precision
  14  
  15  $VERSION = '1.59';
  16  require 5.006;
  17  
  18  require Exporter;
  19  @ISA        = qw/Math::BigInt/;
  20  @EXPORT_OK    = qw/bpi/;
  21  
  22  use strict;
  23  # $_trap_inf/$_trap_nan are internal and should never be accessed from outside
  24  use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode
  25          $upgrade $downgrade $_trap_nan $_trap_inf/;
  26  my $class = "Math::BigFloat";
  27  
  28  use overload
  29  '<=>'    =>    sub { my $rc = $_[2] ?
  30                        ref($_[0])->bcmp($_[1],$_[0]) : 
  31                        ref($_[0])->bcmp($_[0],$_[1]); 
  32                $rc = 1 unless defined $rc;
  33                $rc <=> 0;
  34          },
  35  # we need '>=' to get things like "1 >= NaN" right:
  36  '>='    =>    sub { my $rc = $_[2] ?
  37                        ref($_[0])->bcmp($_[1],$_[0]) : 
  38                        ref($_[0])->bcmp($_[0],$_[1]);
  39                # if there was a NaN involved, return false
  40                return '' unless defined $rc;
  41                $rc >= 0;
  42          },
  43  'int'    =>    sub { $_[0]->as_number() },        # 'trunc' to bigint
  44  ;
  45  
  46  ##############################################################################
  47  # global constants, flags and assorted stuff
  48  
  49  # the following are public, but their usage is not recommended. Use the
  50  # accessor methods instead.
  51  
  52  # class constants, use Class->constant_name() to access
  53  # one of 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'
  54  $round_mode = 'even';
  55  $accuracy   = undef;
  56  $precision  = undef;
  57  $div_scale  = 40;
  58  
  59  $upgrade = undef;
  60  $downgrade = undef;
  61  # the package we are using for our private parts, defaults to:
  62  # Math::BigInt->config()->{lib}
  63  my $MBI = 'Math::BigInt::FastCalc';
  64  
  65  # are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
  66  $_trap_nan = 0;
  67  # the same for infinity
  68  $_trap_inf = 0;
  69  
  70  # constant for easier life
  71  my $nan = 'NaN'; 
  72  
  73  my $IMPORT = 0;    # was import() called yet? used to make require work
  74  
  75  # some digits of accuracy for blog(undef,10); which we use in blog() for speed
  76  my $LOG_10 = 
  77   '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
  78  my $LOG_10_A = length($LOG_10)-1;
  79  # ditto for log(2)
  80  my $LOG_2 = 
  81   '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
  82  my $LOG_2_A = length($LOG_2)-1;
  83  my $HALF = '0.5';            # made into an object if nec.
  84  
  85  ##############################################################################
  86  # the old code had $rnd_mode, so we need to support it, too
  87  
  88  sub TIESCALAR   { my ($class) = @_; bless \$round_mode, $class; }
  89  sub FETCH       { return $round_mode; }
  90  sub STORE       { $rnd_mode = $_[0]->round_mode($_[1]); }
  91  
  92  BEGIN
  93    {
  94    # when someone sets $rnd_mode, we catch this and check the value to see
  95    # whether it is valid or not. 
  96    $rnd_mode   = 'even'; tie $rnd_mode, 'Math::BigFloat';
  97  
  98    # we need both of them in this package:
  99    *as_int = \&as_number;
 100    }
 101   
 102  ##############################################################################
 103  
 104  {
 105    # valid method aliases for AUTOLOAD
 106    my %methods = map { $_ => 1 }  
 107     qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
 108          fint facmp fcmp fzero fnan finf finc fdec ffac fneg
 109      fceil ffloor frsft flsft fone flog froot fexp
 110        /;
 111    # valid methods that can be handed up (for AUTOLOAD)
 112    my %hand_ups = map { $_ => 1 }  
 113     qw / is_nan is_inf is_negative is_positive is_pos is_neg
 114          accuracy precision div_scale round_mode fabs fnot
 115          objectify upgrade downgrade
 116      bone binf bnan bzero
 117      bsub
 118        /;
 119  
 120    sub _method_alias { exists $methods{$_[0]||''}; } 
 121    sub _method_hand_up { exists $hand_ups{$_[0]||''}; } 
 122  }
 123  
 124  ##############################################################################
 125  # constructors
 126  
 127  sub new 
 128    {
 129    # create a new BigFloat object from a string or another bigfloat object. 
 130    # _e: exponent
 131    # _m: mantissa
 132    # sign  => sign (+/-), or "NaN"
 133  
 134    my ($class,$wanted,@r) = @_;
 135  
 136    # avoid numify-calls by not using || on $wanted!
 137    return $class->bzero() if !defined $wanted;    # default to 0
 138    return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
 139  
 140    $class->import() if $IMPORT == 0;             # make require work
 141  
 142    my $self = {}; bless $self, $class;
 143    # shortcut for bigints and its subclasses
 144    if ((ref($wanted)) && UNIVERSAL::can( $wanted, "as_number"))
 145      {
 146      $self->{_m} = $wanted->as_number()->{value}; # get us a bigint copy
 147      $self->{_e} = $MBI->_zero();
 148      $self->{_es} = '+';
 149      $self->{sign} = $wanted->sign();
 150      return $self->bnorm();
 151      }
 152    # else: got a string or something maskerading as number (with overload)
 153  
 154    # handle '+inf', '-inf' first
 155    if ($wanted =~ /^[+-]?inf\z/)
 156      {
 157      return $downgrade->new($wanted) if $downgrade;
 158  
 159      $self->{sign} = $wanted;        # set a default sign for bstr()
 160      return $self->binf($wanted);
 161      }
 162  
 163    # shortcut for simple forms like '12' that neither have trailing nor leading
 164    # zeros
 165    if ($wanted =~ /^([+-]?)([1-9][0-9]*[1-9])$/)
 166      {
 167      $self->{_e} = $MBI->_zero();
 168      $self->{_es} = '+';
 169      $self->{sign} = $1 || '+';
 170      $self->{_m} = $MBI->_new($2);
 171      return $self->round(@r) if !$downgrade;
 172      }
 173  
 174    my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
 175    if (!ref $mis)
 176      {
 177      if ($_trap_nan)
 178        {
 179        require Carp;
 180        Carp::croak ("$wanted is not a number initialized to $class");
 181        }
 182      
 183      return $downgrade->bnan() if $downgrade;
 184      
 185      $self->{_e} = $MBI->_zero();
 186      $self->{_es} = '+';
 187      $self->{_m} = $MBI->_zero();
 188      $self->{sign} = $nan;
 189      }
 190    else
 191      {
 192      # make integer from mantissa by adjusting exp, then convert to int
 193      $self->{_e} = $MBI->_new($$ev);        # exponent
 194      $self->{_es} = $$es || '+';
 195      my $mantissa = "$$miv$$mfv";         # create mant.
 196      $mantissa =~ s/^0+(\d)/$1/;            # strip leading zeros
 197      $self->{_m} = $MBI->_new($mantissa);     # create mant.
 198  
 199      # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
 200      if (CORE::length($$mfv) != 0)
 201        {
 202        my $len = $MBI->_new( CORE::length($$mfv));
 203        ($self->{_e}, $self->{_es}) =
 204      _e_sub ($self->{_e}, $len, $self->{_es}, '+');
 205        }
 206      # we can only have trailing zeros on the mantissa if $$mfv eq ''
 207      else
 208        {
 209        # Use a regexp to count the trailing zeros in $$miv instead of _zeros()
 210        # because that is faster, especially when _m is not stored in base 10.
 211        my $zeros = 0; $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/; 
 212        if ($zeros != 0)
 213          {
 214          my $z = $MBI->_new($zeros);
 215          # turn '120e2' into '12e3'
 216          $MBI->_rsft ( $self->{_m}, $z, 10);
 217          ($self->{_e}, $self->{_es}) =
 218        _e_add ( $self->{_e}, $z, $self->{_es}, '+');
 219          }
 220        }
 221      $self->{sign} = $$mis;
 222  
 223      # for something like 0Ey, set y to 1, and -0 => +0
 224      # Check $$miv for being '0' and $$mfv eq '', because otherwise _m could not
 225      # have become 0. That's faster than to call $MBI->_is_zero().
 226      $self->{sign} = '+', $self->{_e} = $MBI->_one()
 227       if $$miv eq '0' and $$mfv eq '';
 228  
 229      return $self->round(@r) if !$downgrade;
 230      }
 231    # if downgrade, inf, NaN or integers go down
 232  
 233    if ($downgrade && $self->{_es} eq '+')
 234      {
 235      if ($MBI->_is_zero( $self->{_e} ))
 236        {
 237        return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
 238        }
 239      return $downgrade->new($self->bsstr()); 
 240      }
 241    $self->bnorm()->round(@r);            # first normalize, then round
 242    }
 243  
 244  sub copy
 245    {
 246    # if two arguments, the first one is the class to "swallow" subclasses
 247    if (@_ > 1)
 248      {
 249      my  $self = bless {
 250      sign => $_[1]->{sign}, 
 251      _es => $_[1]->{_es}, 
 252      _m => $MBI->_copy($_[1]->{_m}),
 253      _e => $MBI->_copy($_[1]->{_e}),
 254      }, $_[0] if @_ > 1;
 255  
 256      $self->{_a} = $_[1]->{_a} if defined $_[1]->{_a};
 257      $self->{_p} = $_[1]->{_p} if defined $_[1]->{_p};
 258      return $self;
 259      }
 260  
 261    my $self = bless {
 262      sign => $_[0]->{sign}, 
 263      _es => $_[0]->{_es}, 
 264      _m => $MBI->_copy($_[0]->{_m}),
 265      _e => $MBI->_copy($_[0]->{_e}),
 266      }, ref($_[0]);
 267  
 268    $self->{_a} = $_[0]->{_a} if defined $_[0]->{_a};
 269    $self->{_p} = $_[0]->{_p} if defined $_[0]->{_p};
 270    $self;
 271    }
 272  
 273  sub _bnan
 274    {
 275    # used by parent class bone() to initialize number to NaN
 276    my $self = shift;
 277    
 278    if ($_trap_nan)
 279      {
 280      require Carp;
 281      my $class = ref($self);
 282      Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
 283      }
 284  
 285    $IMPORT=1;                    # call our import only once
 286    $self->{_m} = $MBI->_zero();
 287    $self->{_e} = $MBI->_zero();
 288    $self->{_es} = '+';
 289    }
 290  
 291  sub _binf
 292    {
 293    # used by parent class bone() to initialize number to +-inf
 294    my $self = shift;
 295    
 296    if ($_trap_inf)
 297      {
 298      require Carp;
 299      my $class = ref($self);
 300      Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
 301      }
 302  
 303    $IMPORT=1;                    # call our import only once
 304    $self->{_m} = $MBI->_zero();
 305    $self->{_e} = $MBI->_zero();
 306    $self->{_es} = '+';
 307    }
 308  
 309  sub _bone
 310    {
 311    # used by parent class bone() to initialize number to 1
 312    my $self = shift;
 313    $IMPORT=1;                    # call our import only once
 314    $self->{_m} = $MBI->_one();
 315    $self->{_e} = $MBI->_zero();
 316    $self->{_es} = '+';
 317    }
 318  
 319  sub _bzero
 320    {
 321    # used by parent class bone() to initialize number to 0
 322    my $self = shift;
 323    $IMPORT=1;                    # call our import only once
 324    $self->{_m} = $MBI->_zero();
 325    $self->{_e} = $MBI->_one();
 326    $self->{_es} = '+';
 327    }
 328  
 329  sub isa
 330    {
 331    my ($self,$class) = @_;
 332    return if $class =~ /^Math::BigInt/;        # we aren't one of these
 333    UNIVERSAL::isa($self,$class);
 334    }
 335  
 336  sub config
 337    {
 338    # return (later set?) configuration data as hash ref
 339    my $class = shift || 'Math::BigFloat';
 340  
 341    if (@_ == 1 && ref($_[0]) ne 'HASH')
 342      {
 343      my $cfg = $class->SUPER::config();
 344      return $cfg->{$_[0]};
 345      }
 346  
 347    my $cfg = $class->SUPER::config(@_);
 348  
 349    # now we need only to override the ones that are different from our parent
 350    $cfg->{class} = $class;
 351    $cfg->{with} = $MBI;
 352    $cfg;
 353    }
 354  
 355  ##############################################################################
 356  # string conversation
 357  
 358  sub bstr 
 359    {
 360    # (ref to BFLOAT or num_str ) return num_str
 361    # Convert number from internal format to (non-scientific) string format.
 362    # internal format is always normalized (no leading zeros, "-0" => "+0")
 363    my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
 364  
 365    if ($x->{sign} !~ /^[+-]$/)
 366      {
 367      return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
 368      return 'inf';                                       # +inf
 369      }
 370  
 371    my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
 372  
 373    # $x is zero?
 374    my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
 375    if ($not_zero)
 376      {
 377      $es = $MBI->_str($x->{_m});
 378      $len = CORE::length($es);
 379      my $e = $MBI->_num($x->{_e});    
 380      $e = -$e if $x->{_es} eq '-';
 381      if ($e < 0)
 382        {
 383        $dot = '';
 384        # if _e is bigger than a scalar, the following will blow your memory
 385        if ($e <= -$len)
 386          {
 387          my $r = abs($e) - $len;
 388          $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
 389          }
 390        else
 391          {
 392          substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
 393          $cad = -$cad if $x->{_es} eq '-';
 394          }
 395        }
 396      elsif ($e > 0)
 397        {
 398        # expand with zeros
 399        $es .= '0' x $e; $len += $e; $cad = 0;
 400        }
 401      } # if not zero
 402  
 403    $es = '-'.$es if $x->{sign} eq '-';
 404    # if set accuracy or precision, pad with zeros on the right side
 405    if ((defined $x->{_a}) && ($not_zero))
 406      {
 407      # 123400 => 6, 0.1234 => 4, 0.001234 => 4
 408      my $zeros = $x->{_a} - $cad;        # cad == 0 => 12340
 409      $zeros = $x->{_a} - $len if $cad != $len;
 410      $es .= $dot.'0' x $zeros if $zeros > 0;
 411      }
 412    elsif ((($x->{_p} || 0) < 0))
 413      {
 414      # 123400 => 6, 0.1234 => 4, 0.001234 => 6
 415      my $zeros = -$x->{_p} + $cad;
 416      $es .= $dot.'0' x $zeros if $zeros > 0;
 417      }
 418    $es;
 419    }
 420  
 421  sub bsstr
 422    {
 423    # (ref to BFLOAT or num_str ) return num_str
 424    # Convert number from internal format to scientific string format.
 425    # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
 426    my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
 427  
 428    if ($x->{sign} !~ /^[+-]$/)
 429      {
 430      return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
 431      return 'inf';                                       # +inf
 432      }
 433    my $sep = 'e'.$x->{_es};
 434    my $sign = $x->{sign}; $sign = '' if $sign eq '+';
 435    $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
 436    }
 437      
 438  sub numify 
 439    {
 440    # Make a number from a BigFloat object
 441    # simple return a string and let Perl's atoi()/atof() handle the rest
 442    my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
 443    $x->bsstr(); 
 444    }
 445  
 446  ##############################################################################
 447  # public stuff (usually prefixed with "b")
 448  
 449  sub bneg
 450    {
 451    # (BINT or num_str) return BINT
 452    # negate number or make a negated number from string
 453    my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
 454  
 455    return $x if $x->modify('bneg');
 456  
 457    # for +0 dont negate (to have always normalized +0). Does nothing for 'NaN'
 458    $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
 459    $x;
 460    }
 461  
 462  # tels 2001-08-04 
 463  # XXX TODO this must be overwritten and return NaN for non-integer values
 464  # band(), bior(), bxor(), too
 465  #sub bnot
 466  #  {
 467  #  $class->SUPER::bnot($class,@_);
 468  #  }
 469  
 470  sub bcmp 
 471    {
 472    # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
 473  
 474    # set up parameters
 475    my ($self,$x,$y) = (ref($_[0]),@_);
 476    # objectify is costly, so avoid it
 477    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 478      {
 479      ($self,$x,$y) = objectify(2,@_);
 480      }
 481  
 482    return $upgrade->bcmp($x,$y) if defined $upgrade &&
 483      ((!$x->isa($self)) || (!$y->isa($self)));
 484  
 485    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
 486      {
 487      # handle +-inf and NaN
 488      return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 489      return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
 490      return +1 if $x->{sign} eq '+inf';
 491      return -1 if $x->{sign} eq '-inf';
 492      return -1 if $y->{sign} eq '+inf';
 493      return +1;
 494      }
 495  
 496    # check sign for speed first
 497    return 1 if $x->{sign} eq '+' && $y->{sign} eq '-';    # does also 0 <=> -y
 498    return -1 if $x->{sign} eq '-' && $y->{sign} eq '+';    # does also -x <=> 0
 499  
 500    # shortcut 
 501    my $xz = $x->is_zero();
 502    my $yz = $y->is_zero();
 503    return 0 if $xz && $yz;                # 0 <=> 0
 504    return -1 if $xz && $y->{sign} eq '+';        # 0 <=> +y
 505    return 1 if $yz && $x->{sign} eq '+';            # +x <=> 0
 506  
 507    # adjust so that exponents are equal
 508    my $lxm = $MBI->_len($x->{_m});
 509    my $lym = $MBI->_len($y->{_m});
 510    # the numify somewhat limits our length, but makes it much faster
 511    my ($xes,$yes) = (1,1);
 512    $xes = -1 if $x->{_es} ne '+';
 513    $yes = -1 if $y->{_es} ne '+';
 514    my $lx = $lxm + $xes * $MBI->_num($x->{_e});
 515    my $ly = $lym + $yes * $MBI->_num($y->{_e});
 516    my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
 517    return $l <=> 0 if $l != 0;
 518    
 519    # lengths (corrected by exponent) are equal
 520    # so make mantissa equal length by padding with zero (shift left)
 521    my $diff = $lxm - $lym;
 522    my $xm = $x->{_m};        # not yet copy it
 523    my $ym = $y->{_m};
 524    if ($diff > 0)
 525      {
 526      $ym = $MBI->_copy($y->{_m});
 527      $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
 528      }
 529    elsif ($diff < 0)
 530      {
 531      $xm = $MBI->_copy($x->{_m});
 532      $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
 533      }
 534    my $rc = $MBI->_acmp($xm,$ym);
 535    $rc = -$rc if $x->{sign} eq '-';        # -124 < -123
 536    $rc <=> 0;
 537    }
 538  
 539  sub bacmp 
 540    {
 541    # Compares 2 values, ignoring their signs. 
 542    # Returns one of undef, <0, =0, >0. (suitable for sort)
 543    
 544    # set up parameters
 545    my ($self,$x,$y) = (ref($_[0]),@_);
 546    # objectify is costly, so avoid it
 547    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 548      {
 549      ($self,$x,$y) = objectify(2,@_);
 550      }
 551  
 552    return $upgrade->bacmp($x,$y) if defined $upgrade &&
 553      ((!$x->isa($self)) || (!$y->isa($self)));
 554  
 555    # handle +-inf and NaN's
 556    if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
 557      {
 558      return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 559      return 0 if ($x->is_inf() && $y->is_inf());
 560      return 1 if ($x->is_inf() && !$y->is_inf());
 561      return -1;
 562      }
 563  
 564    # shortcut 
 565    my $xz = $x->is_zero();
 566    my $yz = $y->is_zero();
 567    return 0 if $xz && $yz;                # 0 <=> 0
 568    return -1 if $xz && !$yz;                # 0 <=> +y
 569    return 1 if $yz && !$xz;                # +x <=> 0
 570  
 571    # adjust so that exponents are equal
 572    my $lxm = $MBI->_len($x->{_m});
 573    my $lym = $MBI->_len($y->{_m});
 574    my ($xes,$yes) = (1,1);
 575    $xes = -1 if $x->{_es} ne '+';
 576    $yes = -1 if $y->{_es} ne '+';
 577    # the numify somewhat limits our length, but makes it much faster
 578    my $lx = $lxm + $xes * $MBI->_num($x->{_e});
 579    my $ly = $lym + $yes * $MBI->_num($y->{_e});
 580    my $l = $lx - $ly;
 581    return $l <=> 0 if $l != 0;
 582    
 583    # lengths (corrected by exponent) are equal
 584    # so make mantissa equal-length by padding with zero (shift left)
 585    my $diff = $lxm - $lym;
 586    my $xm = $x->{_m};        # not yet copy it
 587    my $ym = $y->{_m};
 588    if ($diff > 0)
 589      {
 590      $ym = $MBI->_copy($y->{_m});
 591      $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
 592      }
 593    elsif ($diff < 0)
 594      {
 595      $xm = $MBI->_copy($x->{_m});
 596      $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
 597      }
 598    $MBI->_acmp($xm,$ym);
 599    }
 600  
 601  sub badd 
 602    {
 603    # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
 604    # return result as BFLOAT
 605  
 606    # set up parameters
 607    my ($self,$x,$y,@r) = (ref($_[0]),@_);
 608    # objectify is costly, so avoid it
 609    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 610      {
 611      ($self,$x,$y,@r) = objectify(2,@_);
 612      }
 613   
 614    return $x if $x->modify('badd');
 615  
 616    # inf and NaN handling
 617    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
 618      {
 619      # NaN first
 620      return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 621      # inf handling
 622      if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
 623        {
 624        # +inf++inf or -inf+-inf => same, rest is NaN
 625        return $x if $x->{sign} eq $y->{sign};
 626        return $x->bnan();
 627        }
 628      # +-inf + something => +inf; something +-inf => +-inf
 629      $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
 630      return $x;
 631      }
 632  
 633    return $upgrade->badd($x,$y,@r) if defined $upgrade &&
 634     ((!$x->isa($self)) || (!$y->isa($self)));
 635  
 636    $r[3] = $y;                        # no push!
 637  
 638    # speed: no add for 0+y or x+0
 639    return $x->bround(@r) if $y->is_zero();        # x+0
 640    if ($x->is_zero())                    # 0+y
 641      {
 642      # make copy, clobbering up x (modify in place!)
 643      $x->{_e} = $MBI->_copy($y->{_e});
 644      $x->{_es} = $y->{_es};
 645      $x->{_m} = $MBI->_copy($y->{_m});
 646      $x->{sign} = $y->{sign} || $nan;
 647      return $x->round(@r);
 648      }
 649   
 650    # take lower of the two e's and adapt m1 to it to match m2
 651    my $e = $y->{_e};
 652    $e = $MBI->_zero() if !defined $e;        # if no BFLOAT?
 653    $e = $MBI->_copy($e);                # make copy (didn't do it yet)
 654  
 655    my $es;
 656  
 657    ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});
 658  
 659    my $add = $MBI->_copy($y->{_m});
 660  
 661    if ($es eq '-')                # < 0
 662      {
 663      $MBI->_lsft( $x->{_m}, $e, 10);
 664      ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
 665      }
 666    elsif (!$MBI->_is_zero($e))            # > 0
 667      {
 668      $MBI->_lsft($add, $e, 10);
 669      }
 670    # else: both e are the same, so just leave them
 671  
 672    if ($x->{sign} eq $y->{sign})
 673      {
 674      # add
 675      $x->{_m} = $MBI->_add($x->{_m}, $add);
 676      }
 677    else
 678      {
 679      ($x->{_m}, $x->{sign}) = 
 680       _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
 681      }
 682  
 683    # delete trailing zeros, then round
 684    $x->bnorm()->round(@r);
 685    }
 686  
 687  # sub bsub is inherited from Math::BigInt!
 688  
 689  sub binc
 690    {
 691    # increment arg by one
 692    my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 693  
 694    return $x if $x->modify('binc');
 695  
 696    if ($x->{_es} eq '-')
 697      {
 698      return $x->badd($self->bone(),@r);    #  digits after dot
 699      }
 700  
 701    if (!$MBI->_is_zero($x->{_e}))        # _e == 0 for NaN, inf, -inf
 702      {
 703      # 1e2 => 100, so after the shift below _m has a '0' as last digit
 704      $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10);    # 1e2 => 100
 705      $x->{_e} = $MBI->_zero();                # normalize
 706      $x->{_es} = '+';
 707      # we know that the last digit of $x will be '1' or '9', depending on the
 708      # sign
 709      }
 710    # now $x->{_e} == 0
 711    if ($x->{sign} eq '+')
 712      {
 713      $MBI->_inc($x->{_m});
 714      return $x->bnorm()->bround(@r);
 715      }
 716    elsif ($x->{sign} eq '-')
 717      {
 718      $MBI->_dec($x->{_m});
 719      $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
 720      return $x->bnorm()->bround(@r);
 721      }
 722    # inf, nan handling etc
 723    $x->badd($self->bone(),@r);            # badd() does round 
 724    }
 725  
 726  sub bdec
 727    {
 728    # decrement arg by one
 729    my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 730  
 731    return $x if $x->modify('bdec');
 732  
 733    if ($x->{_es} eq '-')
 734      {
 735      return $x->badd($self->bone('-'),@r);    #  digits after dot
 736      }
 737  
 738    if (!$MBI->_is_zero($x->{_e}))
 739      {
 740      $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10);    # 1e2 => 100
 741      $x->{_e} = $MBI->_zero();                # normalize
 742      $x->{_es} = '+';
 743      }
 744    # now $x->{_e} == 0
 745    my $zero = $x->is_zero();
 746    # <= 0
 747    if (($x->{sign} eq '-') || $zero)
 748      {
 749      $MBI->_inc($x->{_m});
 750      $x->{sign} = '-' if $zero;                # 0 => 1 => -1
 751      $x->{sign} = '+' if $MBI->_is_zero($x->{_m});    # -1 +1 => -0 => +0
 752      return $x->bnorm()->round(@r);
 753      }
 754    # > 0
 755    elsif ($x->{sign} eq '+')
 756      {
 757      $MBI->_dec($x->{_m});
 758      return $x->bnorm()->round(@r);
 759      }
 760    # inf, nan handling etc
 761    $x->badd($self->bone('-'),@r);        # does round
 762    } 
 763  
 764  sub DEBUG () { 0; }
 765  
 766  sub blog
 767    {
 768    my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 769  
 770    return $x if $x->modify('blog');
 771  
 772    # $base > 0, $base != 1; if $base == undef default to $base == e
 773    # $x >= 0
 774  
 775    # we need to limit the accuracy to protect against overflow
 776    my $fallback = 0;
 777    my ($scale,@params);
 778    ($x,@params) = $x->_find_round_parameters($a,$p,$r);
 779  
 780    # also takes care of the "error in _find_round_parameters?" case
 781    return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
 782  
 783    # no rounding at all, so must use fallback
 784    if (scalar @params == 0)
 785      {
 786      # simulate old behaviour
 787      $params[0] = $self->div_scale();    # and round to it as accuracy
 788      $params[1] = undef;            # P = undef
 789      $scale = $params[0]+4;         # at least four more for proper round
 790      $params[2] = $r;            # round mode by caller or undef
 791      $fallback = 1;            # to clear a/p afterwards
 792      }
 793    else
 794      {
 795      # the 4 below is empirical, and there might be cases where it is not
 796      # enough...
 797      $scale = abs($params[0] || $params[1]) + 4;    # take whatever is defined
 798      }
 799  
 800    return $x->bzero(@params) if $x->is_one();
 801    # base not defined => base == Euler's number e
 802    if (defined $base)
 803      {
 804      # make object, since we don't feed it through objectify() to still get the
 805      # case of $base == undef
 806      $base = $self->new($base) unless ref($base);
 807      # $base > 0; $base != 1
 808      return $x->bnan() if $base->is_zero() || $base->is_one() ||
 809        $base->{sign} ne '+';
 810      # if $x == $base, we know the result must be 1.0
 811      if ($x->bcmp($base) == 0)
 812        {
 813        $x->bone('+',@params);
 814        if ($fallback)
 815          {
 816          # clear a/p after round, since user did not request it
 817          delete $x->{_a}; delete $x->{_p};
 818          }
 819        return $x;
 820        }
 821      }
 822  
 823    # when user set globals, they would interfere with our calculation, so
 824    # disable them and later re-enable them
 825    no strict 'refs';
 826    my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
 827    my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
 828    # we also need to disable any set A or P on $x (_find_round_parameters took
 829    # them already into account), since these would interfere, too
 830    delete $x->{_a}; delete $x->{_p};
 831    # need to disable $upgrade in BigInt, to avoid deep recursion
 832    local $Math::BigInt::upgrade = undef;
 833    local $Math::BigFloat::downgrade = undef;
 834  
 835    # upgrade $x if $x is not a BigFloat (handle BigInt input)
 836    # XXX TODO: rebless!
 837    if (!$x->isa('Math::BigFloat'))
 838      {
 839      $x = Math::BigFloat->new($x);
 840      $self = ref($x);
 841      }
 842    
 843    my $done = 0;
 844  
 845    # If the base is defined and an integer, try to calculate integer result
 846    # first. This is very fast, and in case the real result was found, we can
 847    # stop right here.
 848    if (defined $base && $base->is_int() && $x->is_int())
 849      {
 850      my $i = $MBI->_copy( $x->{_m} );
 851      $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
 852      my $int = Math::BigInt->bzero();
 853      $int->{value} = $i;
 854      $int->blog($base->as_number());
 855      # if ($exact)
 856      if ($base->as_number()->bpow($int) == $x)
 857        {
 858        # found result, return it
 859        $x->{_m} = $int->{value};
 860        $x->{_e} = $MBI->_zero();
 861        $x->{_es} = '+';
 862        $x->bnorm();
 863        $done = 1;
 864        }
 865      }
 866  
 867    if ($done == 0)
 868      {
 869      # base is undef, so base should be e (Euler's number), so first calculate the
 870      # log to base e (using reduction by 10 (and probably 2)):
 871      $self->_log_10($x,$scale);
 872  
 873      # and if a different base was requested, convert it
 874      if (defined $base)
 875        {
 876        $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
 877        # not ln, but some other base (don't modify $base)
 878        $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
 879        }
 880      }
 881   
 882    # shortcut to not run through _find_round_parameters again
 883    if (defined $params[0])
 884      {
 885      $x->bround($params[0],$params[2]);        # then round accordingly
 886      }
 887    else
 888      {
 889      $x->bfround($params[1],$params[2]);        # then round accordingly
 890      }
 891    if ($fallback)
 892      {
 893      # clear a/p after round, since user did not request it
 894      delete $x->{_a}; delete $x->{_p};
 895      }
 896    # restore globals
 897    $$abr = $ab; $$pbr = $pb;
 898  
 899    $x;
 900    }
 901  
 902  sub _len_to_steps
 903    {
 904    # Given D (digits in decimal), compute N so that N! (N factorial) is
 905    # at least D digits long. D should be at least 50.
 906    my $d = shift;
 907  
 908    # two constants for the Ramanujan estimate of ln(N!)
 909    my $lg2 = log(2 * 3.14159265) / 2;
 910    my $lg10 = log(10);
 911  
 912    # D = 50 => N => 42, so L = 40 and R = 50
 913    my $l = 40; my $r = $d;
 914  
 915    # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
 916    $l = $l->numify if ref($l);
 917    $r = $r->numify if ref($r);
 918    $lg2 = $lg2->numify if ref($lg2);
 919    $lg10 = $lg10->numify if ref($lg10);
 920  
 921    # binary search for the right value (could this be written as the reverse of lg(n!)?)
 922    while ($r - $l > 1)
 923      {
 924      my $n = int(($r - $l) / 2) + $l;
 925      my $ramanujan = 
 926        int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
 927      $ramanujan > $d ? $r = $n : $l = $n;
 928      }
 929    $l;
 930    }
 931  
 932  sub bnok
 933    {
 934    # Calculate n over k (binomial coefficient or "choose" function) as integer.
 935    # set up parameters
 936    my ($self,$x,$y,@r) = (ref($_[0]),@_);
 937  
 938    # objectify is costly, so avoid it
 939    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 940      {
 941      ($self,$x,$y,@r) = objectify(2,@_);
 942      }
 943  
 944    return $x if $x->modify('bnok');
 945  
 946    return $x->bnan() if $x->is_nan() || $y->is_nan();
 947    return $x->binf() if $x->is_inf();
 948  
 949    my $u = $x->as_int();
 950    $u->bnok($y->as_int());
 951  
 952    $x->{_m} = $u->{value};
 953    $x->{_e} = $MBI->_zero();
 954    $x->{_es} = '+';
 955    $x->{sign} = '+';
 956    $x->bnorm(@r);
 957    }
 958  
 959  sub bexp
 960    {
 961    # Calculate e ** X (Euler's number to the power of X)
 962    my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 963  
 964    return $x if $x->modify('bexp');
 965  
 966    return $x->binf() if $x->{sign} eq '+inf';
 967    return $x->bzero() if $x->{sign} eq '-inf';
 968  
 969    # we need to limit the accuracy to protect against overflow
 970    my $fallback = 0;
 971    my ($scale,@params);
 972    ($x,@params) = $x->_find_round_parameters($a,$p,$r);
 973  
 974    # also takes care of the "error in _find_round_parameters?" case
 975    return $x if $x->{sign} eq 'NaN';
 976  
 977    # no rounding at all, so must use fallback
 978    if (scalar @params == 0)
 979      {
 980      # simulate old behaviour
 981      $params[0] = $self->div_scale();    # and round to it as accuracy
 982      $params[1] = undef;            # P = undef
 983      $scale = $params[0]+4;         # at least four more for proper round
 984      $params[2] = $r;            # round mode by caller or undef
 985      $fallback = 1;            # to clear a/p afterwards
 986      }
 987    else
 988      {
 989      # the 4 below is empirical, and there might be cases where it's not enough...
 990      $scale = abs($params[0] || $params[1]) + 4;    # take whatever is defined
 991      }
 992  
 993    return $x->bone(@params) if $x->is_zero();
 994  
 995    if (!$x->isa('Math::BigFloat'))
 996      {
 997      $x = Math::BigFloat->new($x);
 998      $self = ref($x);
 999      }
1000    
1001    # when user set globals, they would interfere with our calculation, so
1002    # disable them and later re-enable them
1003    no strict 'refs';
1004    my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1005    my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1006    # we also need to disable any set A or P on $x (_find_round_parameters took
1007    # them already into account), since these would interfere, too
1008    delete $x->{_a}; delete $x->{_p};
1009    # need to disable $upgrade in BigInt, to avoid deep recursion
1010    local $Math::BigInt::upgrade = undef;
1011    local $Math::BigFloat::downgrade = undef;
1012  
1013    my $x_org = $x->copy();
1014  
1015    # We use the following Taylor series:
1016  
1017    #           x    x^2   x^3   x^4
1018    #  e = 1 + --- + --- + --- + --- ...
1019    #           1!    2!    3!    4!
1020  
1021    # The difference for each term is X and N, which would result in:
1022    # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term
1023  
1024    # But it is faster to compute exp(1) and then raising it to the
1025    # given power, esp. if $x is really big and an integer because:
1026  
1027    #  * The numerator is always 1, making the computation faster
1028    #  * the series converges faster in the case of x == 1
1029    #  * We can also easily check when we have reached our limit: when the
1030    #    term to be added is smaller than "1E$scale", we can stop - f.i.
1031    #    scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
1032    #  * we can compute the *exact* result by simulating bigrat math:
1033  
1034    #  1   1    gcd(3,4) = 1    1*24 + 1*6    5
1035    #  - + -                  = ---------- =  --                 
1036    #  6   24                      6*24       24
1037  
1038    # We do not compute the gcd() here, but simple do:
1039    #  1   1    1*24 + 1*6   30
1040    #  - + -  = --------- =  --                 
1041    #  6   24       6*24     144
1042  
1043    # In general:
1044    #  a   c    a*d + c*b     and note that c is always 1 and d = (b*f)
1045    #  - + -  = ---------
1046    #  b   d       b*d
1047  
1048    # This leads to:         which can be reduced by b to:
1049    #  a   1     a*b*f + b    a*f + 1
1050    #  - + -   = --------- =  -------
1051    #  b   b*f     b*b*f        b*f
1052  
1053    # The first terms in the series are:
1054  
1055    # 1     1    1    1    1    1     1     1     13700
1056    # -- + -- + -- + -- + -- + --- + --- + ---- = -----
1057    # 1     1    2    6   24   120   720   5040   5040
1058  
1059    # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!
1060  
1061    if ($scale <= 75)
1062      {
1063      # set $x directly from a cached string form
1064      $x->{_m} = $MBI->_new(
1065      "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
1066      $x->{sign} = '+';
1067      $x->{_es} = '-';
1068      $x->{_e} = $MBI->_new(79);
1069      }
1070    else
1071      {
1072      # compute A and B so that e = A / B.
1073   
1074      # After some terms we end up with this, so we use it as a starting point:
1075      my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
1076      my $F = $MBI->_new(42); my $step = 42;
1077  
1078      # Compute how many steps we need to take to get $A and $B sufficiently big
1079      my $steps = _len_to_steps($scale - 4);
1080  #    print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
1081      while ($step++ <= $steps)
1082        {
1083        # calculate $a * $f + 1
1084        $A = $MBI->_mul($A, $F);
1085        $A = $MBI->_inc($A);
1086        # increment f
1087        $F = $MBI->_inc($F);
1088        }
1089      # compute $B as factorial of $steps (this is faster than doing it manually)
1090      my $B = $MBI->_fac($MBI->_new($steps));
1091      
1092  #  print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";
1093  
1094      # compute A/B with $scale digits in the result (truncate, not round)
1095      $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
1096      $A = $MBI->_div( $A, $B );
1097  
1098      $x->{_m} = $A;
1099      $x->{sign} = '+';
1100      $x->{_es} = '-';
1101      $x->{_e} = $MBI->_new($scale);
1102      }
1103  
1104    # $x contains now an estimate of e, with some surplus digits, so we can round
1105    if (!$x_org->is_one())
1106      {
1107      # raise $x to the wanted power and round it in one step:
1108      $x->bpow($x_org, @params);
1109      }
1110    else
1111      {
1112      # else just round the already computed result
1113      delete $x->{_a}; delete $x->{_p};
1114      # shortcut to not run through _find_round_parameters again
1115      if (defined $params[0])
1116        {
1117        $x->bround($params[0],$params[2]);        # then round accordingly
1118        }
1119      else
1120        {
1121        $x->bfround($params[1],$params[2]);        # then round accordingly
1122        }
1123      }
1124    if ($fallback)
1125      {
1126      # clear a/p after round, since user did not request it
1127      delete $x->{_a}; delete $x->{_p};
1128      }
1129    # restore globals
1130    $$abr = $ab; $$pbr = $pb;
1131  
1132    $x;                        # return modified $x
1133    }
1134  
1135  sub _log
1136    {
1137    # internal log function to calculate ln() based on Taylor series.
1138    # Modifies $x in place.
1139    my ($self,$x,$scale) = @_;
1140  
1141    # in case of $x == 1, result is 0
1142    return $x->bzero() if $x->is_one();
1143  
1144    # XXX TODO: rewrite this in a similiar manner to bexp()
1145  
1146    # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
1147  
1148    # u = x-1, v = x+1
1149    #              _                               _
1150    # Taylor:     |    u    1   u^3   1   u^5       |
1151    # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 0
1152    #             |_   v    3   v^3   5   v^5      _|
1153  
1154    # This takes much more steps to calculate the result and is thus not used
1155    # u = x-1
1156    #              _                               _
1157    # Taylor:     |    u    1   u^2   1   u^3       |
1158    # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 1/2
1159    #             |_   x    2   x^2   3   x^3      _|
1160  
1161    my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
1162  
1163    $v = $x->copy(); $v->binc();        # v = x+1
1164    $x->bdec(); $u = $x->copy();        # u = x-1; x = x-1
1165    $x->bdiv($v,$scale);            # first term: u/v
1166    $below = $v->copy();
1167    $over = $u->copy();
1168    $u *= $u; $v *= $v;                # u^2, v^2
1169    $below->bmul($v);                # u^3, v^3
1170    $over->bmul($u);
1171    $factor = $self->new(3); $f = $self->new(2);
1172  
1173    my $steps = 0 if DEBUG;  
1174    $limit = $self->new("1E-". ($scale-1));
1175    while (3 < 5)
1176      {
1177      # we calculate the next term, and add it to the last
1178      # when the next term is below our limit, it won't affect the outcome
1179      # anymore, so we stop
1180  
1181      # calculating the next term simple from over/below will result in quite
1182      # a time hog if the input has many digits, since over and below will
1183      # accumulate more and more digits, and the result will also have many
1184      # digits, but in the end it is rounded to $scale digits anyway. So if we
1185      # round $over and $below first, we save a lot of time for the division
1186      # (not with log(1.2345), but try log (123**123) to see what I mean. This
1187      # can introduce a rounding error if the division result would be f.i.
1188      # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
1189      # if we truncated $over and $below we might get 0.12345. Does this matter
1190      # for the end result? So we give $over and $below 4 more digits to be
1191      # on the safe side (unscientific error handling as usual... :+D
1192  
1193      $next = $over->copy->bround($scale+4)->bdiv(
1194        $below->copy->bmul($factor)->bround($scale+4), 
1195        $scale);
1196  
1197  ## old version:    
1198  ##    $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
1199  
1200      last if $next->bacmp($limit) <= 0;
1201  
1202      delete $next->{_a}; delete $next->{_p};
1203      $x->badd($next);
1204      # calculate things for the next term
1205      $over *= $u; $below *= $v; $factor->badd($f);
1206      if (DEBUG)
1207        {
1208        $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
1209        }
1210      }
1211    print "took $steps steps\n" if DEBUG;
1212    $x->bmul($f);                    # $x *= 2
1213    }
1214  
1215  sub _log_10
1216    {
1217    # Internal log function based on reducing input to the range of 0.1 .. 9.99
1218    # and then "correcting" the result to the proper one. Modifies $x in place.
1219    my ($self,$x,$scale) = @_;
1220  
1221    # Taking blog() from numbers greater than 10 takes a *very long* time, so we
1222    # break the computation down into parts based on the observation that:
1223    #  blog(X*Y) = blog(X) + blog(Y)
1224    # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
1225    # $x is the faster it gets. Since 2*$x takes about 10 times as
1226    # long, we make it faster by about a factor of 100 by dividing $x by 10.
1227  
1228    # The same observation is valid for numbers smaller than 0.1, e.g. computing
1229    # log(1) is fastest, and the further away we get from 1, the longer it takes.
1230    # So we also 'break' this down by multiplying $x with 10 and subtract the
1231    # log(10) afterwards to get the correct result.
1232  
1233    # To get $x even closer to 1, we also divide by 2 and then use log(2) to
1234    # correct for this. For instance if $x is 2.4, we use the formula:
1235    #  blog(2.4 * 2) == blog (1.2) + blog(2)
1236    # and thus calculate only blog(1.2) and blog(2), which is faster in total
1237    # than calculating blog(2.4).
1238  
1239    # In addition, the values for blog(2) and blog(10) are cached.
1240  
1241    # Calculate nr of digits before dot:
1242    my $dbd = $MBI->_num($x->{_e});
1243    $dbd = -$dbd if $x->{_es} eq '-';
1244    $dbd += $MBI->_len($x->{_m});
1245  
1246    # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
1247    # infinite recursion
1248  
1249    my $calc = 1;                    # do some calculation?
1250  
1251    # disable the shortcut for 10, since we need log(10) and this would recurse
1252    # infinitely deep
1253    if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
1254      {
1255      $dbd = 0;                    # disable shortcut
1256      # we can use the cached value in these cases
1257      if ($scale <= $LOG_10_A)
1258        {
1259        $x->bzero(); $x->badd($LOG_10);        # modify $x in place
1260        $calc = 0;                 # no need to calc, but round
1261        }
1262      # if we can't use the shortcut, we continue normally
1263      }
1264    else
1265      {
1266      # disable the shortcut for 2, since we maybe have it cached
1267      if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
1268        {
1269        $dbd = 0;                    # disable shortcut
1270        # we can use the cached value in these cases
1271        if ($scale <= $LOG_2_A)
1272          {
1273          $x->bzero(); $x->badd($LOG_2);        # modify $x in place
1274          $calc = 0;                 # no need to calc, but round
1275          }
1276        # if we can't use the shortcut, we continue normally
1277        }
1278      }
1279  
1280    # if $x = 0.1, we know the result must be 0-log(10)
1281    if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
1282        $MBI->_is_one($x->{_m}))
1283      {
1284      $dbd = 0;                    # disable shortcut
1285      # we can use the cached value in these cases
1286      if ($scale <= $LOG_10_A)
1287        {
1288        $x->bzero(); $x->bsub($LOG_10);
1289        $calc = 0;                 # no need to calc, but round
1290        }
1291      }
1292  
1293    return if $calc == 0;                # already have the result
1294  
1295    # default: these correction factors are undef and thus not used
1296    my $l_10;                # value of ln(10) to A of $scale
1297    my $l_2;                # value of ln(2) to A of $scale
1298  
1299    my $two = $self->new(2);
1300  
1301    # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
1302    # so don't do this shortcut for 1 or 0
1303    if (($dbd > 1) || ($dbd < 0))
1304      {
1305      # convert our cached value to an object if not already (avoid doing this
1306      # at import() time, since not everybody needs this)
1307      $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
1308  
1309      #print "x = $x, dbd = $dbd, calc = $calc\n";
1310      # got more than one digit before the dot, or more than one zero after the
1311      # dot, so do:
1312      #  log(123)    == log(1.23) + log(10) * 2
1313      #  log(0.0123) == log(1.23) - log(10) * 2
1314    
1315      if ($scale <= $LOG_10_A)
1316        {
1317        # use cached value
1318        $l_10 = $LOG_10->copy();        # copy for mul
1319        }
1320      else
1321        {
1322        # else: slower, compute and cache result
1323        # also disable downgrade for this code path
1324        local $Math::BigFloat::downgrade = undef;
1325  
1326        # shorten the time to calculate log(10) based on the following:
1327        # log(1.25 * 8) = log(1.25) + log(8)
1328        #               = log(1.25) + log(2) + log(2) + log(2)
1329  
1330        # first get $l_2 (and possible compute and cache log(2))
1331        $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1332        if ($scale <= $LOG_2_A)
1333          {
1334          # use cached value
1335          $l_2 = $LOG_2->copy();            # copy() for the mul below
1336          }
1337        else
1338          {
1339          # else: slower, compute and cache result
1340          $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1341          $LOG_2 = $l_2->copy();            # cache the result for later
1342                          # the copy() is for mul below
1343          $LOG_2_A = $scale;
1344          }
1345  
1346        # now calculate log(1.25):
1347        $l_10 = $self->new('1.25'); $self->_log($l_10, $scale); # scale+4, actually
1348  
1349        # log(1.25) + log(2) + log(2) + log(2):
1350        $l_10->badd($l_2);
1351        $l_10->badd($l_2);
1352        $l_10->badd($l_2);
1353        $LOG_10 = $l_10->copy();        # cache the result for later
1354                      # the copy() is for mul below
1355        $LOG_10_A = $scale;
1356        }
1357      $dbd-- if ($dbd > 1);         # 20 => dbd=2, so make it dbd=1    
1358      $l_10->bmul( $self->new($dbd));    # log(10) * (digits_before_dot-1)
1359      my $dbd_sign = '+';
1360      if ($dbd < 0)
1361        {
1362        $dbd = -$dbd;
1363        $dbd_sign = '-';
1364        }
1365      ($x->{_e}, $x->{_es}) = 
1366      _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
1367   
1368      }
1369  
1370    # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1371  
1372    ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1373    ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1374  
1375    $HALF = $self->new($HALF) unless ref($HALF);
1376  
1377    my $twos = 0;                # default: none (0 times)    
1378    while ($x->bacmp($HALF) <= 0)        # X <= 0.5
1379      {
1380      $twos--; $x->bmul($two);
1381      }
1382    while ($x->bacmp($two) >= 0)        # X >= 2
1383      {
1384      $twos++; $x->bdiv($two,$scale+4);        # keep all digits
1385      }
1386    # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
1387    # So calculate correction factor based on ln(2):
1388    if ($twos != 0)
1389      {
1390      $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1391      if ($scale <= $LOG_2_A)
1392        {
1393        # use cached value
1394        $l_2 = $LOG_2->copy();            # copy() for the mul below
1395        }
1396      else
1397        {
1398        # else: slower, compute and cache result
1399        # also disable downgrade for this code path
1400        local $Math::BigFloat::downgrade = undef;
1401        $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1402        $LOG_2 = $l_2->copy();            # cache the result for later
1403                          # the copy() is for mul below
1404        $LOG_2_A = $scale;
1405        }
1406      $l_2->bmul($twos);        # * -2 => subtract, * 2 => add
1407      }
1408    
1409    $self->_log($x,$scale);            # need to do the "normal" way
1410    $x->badd($l_10) if defined $l_10;         # correct it by ln(10)
1411    $x->badd($l_2) if defined $l_2;        # and maybe by ln(2)
1412  
1413    # all done, $x contains now the result
1414    $x;
1415    }
1416  
1417  sub blcm 
1418    { 
1419    # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1420    # does not modify arguments, but returns new object
1421    # Lowest Common Multiplicator
1422  
1423    my ($self,@arg) = objectify(0,@_);
1424    my $x = $self->new(shift @arg);
1425    while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); } 
1426    $x;
1427    }
1428  
1429  sub bgcd
1430    {
1431    # (BINT or num_str, BINT or num_str) return BINT
1432    # does not modify arguments, but returns new object
1433  
1434    my $y = shift;
1435    $y = __PACKAGE__->new($y) if !ref($y);
1436    my $self = ref($y);
1437    my $x = $y->copy()->babs();            # keep arguments
1438  
1439    return $x->bnan() if $x->{sign} !~ /^[+-]$/    # x NaN?
1440      || !$x->is_int();            # only for integers now
1441  
1442    while (@_)
1443      {
1444      my $t = shift; $t = $self->new($t) if !ref($t);
1445      $y = $t->copy()->babs();
1446      
1447      return $x->bnan() if $y->{sign} !~ /^[+-]$/    # y NaN?
1448           || !$y->is_int();            # only for integers now
1449  
1450      # greatest common divisor
1451      while (! $y->is_zero())
1452        {
1453        ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
1454        }
1455  
1456      last if $x->is_one();
1457      }
1458    $x;
1459    }
1460  
1461  ##############################################################################
1462  
1463  sub _e_add
1464    {
1465    # Internal helper sub to take two positive integers and their signs and
1466    # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')), 
1467    # output ($CALC,('+'|'-'))
1468    my ($x,$y,$xs,$ys) = @_;
1469  
1470    # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1471    if ($xs eq $ys)
1472      {
1473      $x = $MBI->_add ($x, $y );        # a+b
1474      # the sign follows $xs
1475      return ($x, $xs);
1476      }
1477  
1478    my $a = $MBI->_acmp($x,$y);
1479    if ($a > 0)
1480      {
1481      $x = $MBI->_sub ($x , $y);                # abs sub
1482      }
1483    elsif ($a == 0)
1484      {
1485      $x = $MBI->_zero();                    # result is 0
1486      $xs = '+';
1487      }
1488    else # a < 0
1489      {
1490      $x = $MBI->_sub ( $y, $x, 1 );            # abs sub
1491      $xs = $ys;
1492      }
1493    ($x,$xs);
1494    }
1495  
1496  sub _e_sub
1497    {
1498    # Internal helper sub to take two positive integers and their signs and
1499    # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')), 
1500    # output ($CALC,('+'|'-'))
1501    my ($x,$y,$xs,$ys) = @_;
1502  
1503    # flip sign
1504    $ys =~ tr/+-/-+/;
1505    _e_add($x,$y,$xs,$ys);        # call add (does subtract now)
1506    }
1507  
1508  ###############################################################################
1509  # is_foo methods (is_negative, is_positive are inherited from BigInt)
1510  
1511  sub is_int
1512    {
1513    # return true if arg (BFLOAT or num_str) is an integer
1514    my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1515  
1516    (($x->{sign} =~ /^[+-]$/) &&            # NaN and +-inf aren't
1517     ($x->{_es} eq '+')) ? 1 : 0;            # 1e-1 => no integer
1518    }
1519  
1520  sub is_zero
1521    {
1522    # return true if arg (BFLOAT or num_str) is zero
1523    my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1524  
1525    ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m})) ? 1 : 0;
1526    }
1527  
1528  sub is_one
1529    {
1530    # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1531    my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1532  
1533    $sign = '+' if !defined $sign || $sign ne '-';
1534  
1535    ($x->{sign} eq $sign && 
1536     $MBI->_is_zero($x->{_e}) &&
1537     $MBI->_is_one($x->{_m}) ) ? 1 : 0; 
1538    }
1539  
1540  sub is_odd
1541    {
1542    # return true if arg (BFLOAT or num_str) is odd or false if even
1543    my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1544    
1545    (($x->{sign} =~ /^[+-]$/) &&        # NaN & +-inf aren't
1546     ($MBI->_is_zero($x->{_e})) &&
1547     ($MBI->_is_odd($x->{_m}))) ? 1 : 0; 
1548    }
1549  
1550  sub is_even
1551    {
1552    # return true if arg (BINT or num_str) is even or false if odd
1553    my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1554  
1555    (($x->{sign} =~ /^[+-]$/) &&            # NaN & +-inf aren't
1556     ($x->{_es} eq '+') &&             # 123.45 isn't
1557     ($MBI->_is_even($x->{_m}))) ? 1 : 0;        # but 1200 is
1558    }
1559  
1560  sub bmul
1561    { 
1562    # multiply two numbers
1563    
1564    # set up parameters
1565    my ($self,$x,$y,@r) = (ref($_[0]),@_);
1566    # objectify is costly, so avoid it
1567    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1568      {
1569      ($self,$x,$y,@r) = objectify(2,@_);
1570      }
1571  
1572    return $x if $x->modify('bmul');
1573  
1574    return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1575  
1576    # inf handling
1577    if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1578      {
1579      return $x->bnan() if $x->is_zero() || $y->is_zero(); 
1580      # result will always be +-inf:
1581      # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1582      # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1583      return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1584      return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1585      return $x->binf('-');
1586      }
1587    
1588    return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1589     ((!$x->isa($self)) || (!$y->isa($self)));
1590  
1591    # aEb * cEd = (a*c)E(b+d)
1592    $MBI->_mul($x->{_m},$y->{_m});
1593    ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1594  
1595    $r[3] = $y;                # no push!
1596  
1597    # adjust sign:
1598    $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1599    $x->bnorm->round(@r);
1600    }
1601  
1602  sub bmuladd
1603    { 
1604    # multiply two numbers and add the third to the result
1605    
1606    # set up parameters
1607    my ($self,$x,$y,$z,@r) = (ref($_[0]),@_);
1608    # objectify is costly, so avoid it
1609    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1610      {
1611      ($self,$x,$y,$z,@r) = objectify(3,@_);
1612      }
1613  
1614    return $x if $x->modify('bmuladd');
1615  
1616    return $x->bnan() if (($x->{sign} eq $nan) ||
1617              ($y->{sign} eq $nan) ||
1618              ($z->{sign} eq $nan));
1619  
1620    # inf handling
1621    if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1622      {
1623      return $x->bnan() if $x->is_zero() || $y->is_zero(); 
1624      # result will always be +-inf:
1625      # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1626      # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1627      return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1628      return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1629      return $x->binf('-');
1630      }
1631  
1632    return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1633     ((!$x->isa($self)) || (!$y->isa($self)));
1634  
1635    # aEb * cEd = (a*c)E(b+d)
1636    $MBI->_mul($x->{_m},$y->{_m});
1637    ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1638  
1639    $r[3] = $y;                # no push!
1640  
1641    # adjust sign:
1642    $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1643  
1644    # z=inf handling (z=NaN handled above)
1645    $x->{sign} = $z->{sign}, return $x if $z->{sign} =~ /^[+-]inf$/;
1646  
1647    # take lower of the two e's and adapt m1 to it to match m2
1648    my $e = $z->{_e};
1649    $e = $MBI->_zero() if !defined $e;        # if no BFLOAT?
1650    $e = $MBI->_copy($e);                # make copy (didn't do it yet)
1651  
1652    my $es;
1653  
1654    ($e,$es) = _e_sub($e, $x->{_e}, $z->{_es} || '+', $x->{_es});
1655  
1656    my $add = $MBI->_copy($z->{_m});
1657  
1658    if ($es eq '-')                # < 0
1659      {
1660      $MBI->_lsft( $x->{_m}, $e, 10);
1661      ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
1662      }
1663    elsif (!$MBI->_is_zero($e))            # > 0
1664      {
1665      $MBI->_lsft($add, $e, 10);
1666      }
1667    # else: both e are the same, so just leave them
1668  
1669    if ($x->{sign} eq $z->{sign})
1670      {
1671      # add
1672      $x->{_m} = $MBI->_add($x->{_m}, $add);
1673      }
1674    else
1675      {
1676      ($x->{_m}, $x->{sign}) = 
1677       _e_add($x->{_m}, $add, $x->{sign}, $z->{sign});
1678      }
1679  
1680    # delete trailing zeros, then round
1681    $x->bnorm()->round(@r);
1682    }
1683  
1684  sub bdiv 
1685    {
1686    # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return 
1687    # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1688  
1689    # set up parameters
1690    my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1691    # objectify is costly, so avoid it
1692    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1693      {
1694      ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1695      }
1696  
1697    return $x if $x->modify('bdiv');
1698  
1699    return $self->_div_inf($x,$y)
1700     if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1701  
1702    # x== 0 # also: or y == 1 or y == -1
1703    return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1704  
1705    # upgrade ?
1706    return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1707  
1708    # we need to limit the accuracy to protect against overflow
1709    my $fallback = 0;
1710    my (@params,$scale);
1711    ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1712  
1713    return $x if $x->is_nan();        # error in _find_round_parameters?
1714  
1715    # no rounding at all, so must use fallback
1716    if (scalar @params == 0)
1717      {
1718      # simulate old behaviour
1719      $params[0] = $self->div_scale();    # and round to it as accuracy
1720      $scale = $params[0]+4;         # at least four more for proper round
1721      $params[2] = $r;            # round mode by caller or undef
1722      $fallback = 1;            # to clear a/p afterwards
1723      }
1724    else
1725      {
1726      # the 4 below is empirical, and there might be cases where it is not
1727      # enough...
1728      $scale = abs($params[0] || $params[1]) + 4;    # take whatever is defined
1729      }
1730  
1731    my $rem; $rem = $self->bzero() if wantarray;
1732  
1733    $y = $self->new($y) unless $y->isa('Math::BigFloat');
1734  
1735    my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
1736    $scale = $lx if $lx > $scale;
1737    $scale = $ly if $ly > $scale;
1738    my $diff = $ly - $lx;
1739    $scale += $diff if $diff > 0;        # if lx << ly, but not if ly << lx!
1740  
1741    # already handled inf/NaN/-inf above:
1742  
1743    # check that $y is not 1 nor -1 and cache the result:
1744    my $y_not_one = !($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m}));
1745  
1746    # flipping the sign of $y will also flip the sign of $x for the special
1747    # case of $x->bsub($x); so we can catch it below:
1748    my $xsign = $x->{sign};
1749    $y->{sign} =~ tr/+-/-+/;
1750  
1751    if ($xsign ne $x->{sign})
1752      {
1753      # special case of $x /= $x results in 1
1754      $x->bone();            # "fixes" also sign of $y, since $x is $y
1755      }
1756    else
1757      {
1758      # correct $y's sign again
1759      $y->{sign} =~ tr/+-/-+/;
1760      # continue with normal div code:
1761  
1762      # make copy of $x in case of list context for later reminder calculation
1763      if (wantarray && $y_not_one)
1764        {
1765        $rem = $x->copy();
1766        }
1767  
1768      $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+'; 
1769  
1770      # check for / +-1 ( +/- 1E0)
1771      if ($y_not_one)
1772        {
1773        # promote BigInts and it's subclasses (except when already a BigFloat)
1774        $y = $self->new($y) unless $y->isa('Math::BigFloat'); 
1775  
1776        # calculate the result to $scale digits and then round it
1777        # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1778        $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1779        $MBI->_div ($x->{_m},$y->{_m});    # a/c
1780  
1781        # correct exponent of $x
1782        ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1783        # correct for 10**scale
1784        ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1785        $x->bnorm();        # remove trailing 0's
1786        }
1787      } # ende else $x != $y
1788  
1789    # shortcut to not run through _find_round_parameters again
1790    if (defined $params[0])
1791      {
1792      delete $x->{_a};                 # clear before round
1793      $x->bround($params[0],$params[2]);        # then round accordingly
1794      }
1795    else
1796      {
1797      delete $x->{_p};                 # clear before round
1798      $x->bfround($params[1],$params[2]);        # then round accordingly
1799      }
1800    if ($fallback)
1801      {
1802      # clear a/p after round, since user did not request it
1803      delete $x->{_a}; delete $x->{_p};
1804      }
1805  
1806    if (wantarray)
1807      {
1808      if ($y_not_one)
1809        {
1810        $rem->bmod($y,@params);            # copy already done
1811        }
1812      if ($fallback)
1813        {
1814        # clear a/p after round, since user did not request it
1815        delete $rem->{_a}; delete $rem->{_p};
1816        }
1817      return ($x,$rem);
1818      }
1819    $x;
1820    }
1821  
1822  sub bmod 
1823    {
1824    # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder 
1825  
1826    # set up parameters
1827    my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1828    # objectify is costly, so avoid it
1829    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1830      {
1831      ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1832      }
1833  
1834    return $x if $x->modify('bmod');
1835  
1836    # handle NaN, inf, -inf
1837    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1838      {
1839      my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1840      $x->{sign} = $re->{sign};
1841      $x->{_e} = $re->{_e};
1842      $x->{_m} = $re->{_m};
1843      return $x->round($a,$p,$r,$y);
1844      } 
1845    if ($y->is_zero())
1846      {
1847      return $x->bnan() if $x->is_zero();
1848      return $x;
1849      }
1850  
1851    return $x->bzero() if $x->is_zero()
1852   || ($x->is_int() &&
1853    # check that $y == +1 or $y == -1:
1854      ($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m})));
1855  
1856    my $cmp = $x->bacmp($y);            # equal or $x < $y?
1857    return $x->bzero($a,$p) if $cmp == 0;        # $x == $y => result 0
1858  
1859    # only $y of the operands negative? 
1860    my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1861  
1862    $x->{sign} = $y->{sign};                # calc sign first
1863    return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0;    # $x < $y => result $x
1864    
1865    my $ym = $MBI->_copy($y->{_m});
1866    
1867    # 2e1 => 20
1868    $MBI->_lsft( $ym, $y->{_e}, 10) 
1869     if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
1870   
1871    # if $y has digits after dot
1872    my $shifty = 0;            # correct _e of $x by this
1873    if ($y->{_es} eq '-')            # has digits after dot
1874      {
1875      # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1876      $shifty = $MBI->_num($y->{_e});     # no more digits after dot
1877      $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
1878      }
1879    # $ym is now mantissa of $y based on exponent 0
1880  
1881    my $shiftx = 0;            # correct _e of $x by this
1882    if ($x->{_es} eq '-')            # has digits after dot
1883      {
1884      # 123.4 % 20 => 1234 % 200
1885      $shiftx = $MBI->_num($x->{_e});    # no more digits after dot
1886      $MBI->_lsft($ym, $x->{_e}, 10);    # 123 => 1230
1887      }
1888    # 123e1 % 20 => 1230 % 20
1889    if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
1890      {
1891      $MBI->_lsft( $x->{_m}, $x->{_e},10);    # es => '+' here
1892      }
1893  
1894    $x->{_e} = $MBI->_new($shiftx);
1895    $x->{_es} = '+'; 
1896    $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1897    $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
1898    
1899    # now mantissas are equalized, exponent of $x is adjusted, so calc result
1900  
1901    $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
1902  
1903    $x->{sign} = '+' if $MBI->_is_zero($x->{_m});        # fix sign for -0
1904    $x->bnorm();
1905  
1906    if ($neg != 0)    # one of them negative => correct in place
1907      {
1908      my $r = $y - $x;
1909      $x->{_m} = $r->{_m};
1910      $x->{_e} = $r->{_e};
1911      $x->{_es} = $r->{_es};
1912      $x->{sign} = '+' if $MBI->_is_zero($x->{_m});    # fix sign for -0
1913      $x->bnorm();
1914      }
1915  
1916    $x->round($a,$p,$r,$y);    # round and return
1917    }
1918  
1919  sub broot
1920    {
1921    # calculate $y'th root of $x
1922    
1923    # set up parameters
1924    my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1925    # objectify is costly, so avoid it
1926    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1927      {
1928      ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1929      }
1930  
1931    return $x if $x->modify('broot');
1932  
1933    # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1934    return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1935           $y->{sign} !~ /^\+$/;
1936  
1937    return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1938    
1939    # we need to limit the accuracy to protect against overflow
1940    my $fallback = 0;
1941    my (@params,$scale);
1942    ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1943  
1944    return $x if $x->is_nan();        # error in _find_round_parameters?
1945  
1946    # no rounding at all, so must use fallback
1947    if (scalar @params == 0) 
1948      {
1949      # simulate old behaviour
1950      $params[0] = $self->div_scale();    # and round to it as accuracy
1951      $scale = $params[0]+4;         # at least four more for proper round
1952      $params[2] = $r;            # iound mode by caller or undef
1953      $fallback = 1;            # to clear a/p afterwards
1954      }
1955    else
1956      {
1957      # the 4 below is empirical, and there might be cases where it is not
1958      # enough...
1959      $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1960      }
1961  
1962    # when user set globals, they would interfere with our calculation, so
1963    # disable them and later re-enable them
1964    no strict 'refs';
1965    my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1966    my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1967    # we also need to disable any set A or P on $x (_find_round_parameters took
1968    # them already into account), since these would interfere, too
1969    delete $x->{_a}; delete $x->{_p};
1970    # need to disable $upgrade in BigInt, to avoid deep recursion
1971    local $Math::BigInt::upgrade = undef;    # should be really parent class vs MBI
1972  
1973    # remember sign and make $x positive, since -4 ** (1/2) => -2
1974    my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
1975  
1976    my $is_two = 0;
1977    if ($y->isa('Math::BigFloat'))
1978      {
1979      $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
1980      }
1981    else
1982      {
1983      $is_two = ($y == 2);
1984      }
1985  
1986    # normal square root if $y == 2:
1987    if ($is_two)
1988      {
1989      $x->bsqrt($scale+4);
1990      }
1991    elsif ($y->is_one('-'))
1992      {
1993      # $x ** -1 => 1/$x
1994      my $u = $self->bone()->bdiv($x,$scale);
1995      # copy private parts over
1996      $x->{_m} = $u->{_m};
1997      $x->{_e} = $u->{_e};
1998      $x->{_es} = $u->{_es};
1999      }
2000    else
2001      {
2002      # calculate the broot() as integer result first, and if it fits, return
2003      # it rightaway (but only if $x and $y are integer):
2004  
2005      my $done = 0;                # not yet
2006      if ($y->is_int() && $x->is_int())
2007        {
2008        my $i = $MBI->_copy( $x->{_m} );
2009        $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2010        my $int = Math::BigInt->bzero();
2011        $int->{value} = $i;
2012        $int->broot($y->as_number());
2013        # if ($exact)
2014        if ($int->copy()->bpow($y) == $x)
2015          {
2016          # found result, return it
2017          $x->{_m} = $int->{value};
2018          $x->{_e} = $MBI->_zero();
2019          $x->{_es} = '+';
2020          $x->bnorm();
2021          $done = 1;
2022          }
2023        }
2024      if ($done == 0)
2025        {
2026        my $u = $self->bone()->bdiv($y,$scale+4);
2027        delete $u->{_a}; delete $u->{_p};         # otherwise it conflicts
2028        $x->bpow($u,$scale+4);                    # el cheapo
2029        }
2030      }
2031    $x->bneg() if $sign == 1;
2032    
2033    # shortcut to not run through _find_round_parameters again
2034    if (defined $params[0])
2035      {
2036      $x->bround($params[0],$params[2]);        # then round accordingly
2037      }
2038    else
2039      {
2040      $x->bfround($params[1],$params[2]);        # then round accordingly
2041      }
2042    if ($fallback)
2043      {
2044      # clear a/p after round, since user did not request it
2045      delete $x->{_a}; delete $x->{_p};
2046      }
2047    # restore globals
2048    $$abr = $ab; $$pbr = $pb;
2049    $x;
2050    }
2051  
2052  sub bsqrt
2053    { 
2054    # calculate square root
2055    my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2056  
2057    return $x if $x->modify('bsqrt');
2058  
2059    return $x->bnan() if $x->{sign} !~ /^[+]/;    # NaN, -inf or < 0
2060    return $x if $x->{sign} eq '+inf';        # sqrt(inf) == inf
2061    return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
2062  
2063    # we need to limit the accuracy to protect against overflow
2064    my $fallback = 0;
2065    my (@params,$scale);
2066    ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2067  
2068    return $x if $x->is_nan();        # error in _find_round_parameters?
2069  
2070    # no rounding at all, so must use fallback
2071    if (scalar @params == 0) 
2072      {
2073      # simulate old behaviour
2074      $params[0] = $self->div_scale();    # and round to it as accuracy
2075      $scale = $params[0]+4;         # at least four more for proper round
2076      $params[2] = $r;            # round mode by caller or undef
2077      $fallback = 1;            # to clear a/p afterwards
2078      }
2079    else
2080      {
2081      # the 4 below is empirical, and there might be cases where it is not
2082      # enough...
2083      $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2084      }
2085  
2086    # when user set globals, they would interfere with our calculation, so
2087    # disable them and later re-enable them
2088    no strict 'refs';
2089    my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2090    my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2091    # we also need to disable any set A or P on $x (_find_round_parameters took
2092    # them already into account), since these would interfere, too
2093    delete $x->{_a}; delete $x->{_p};
2094    # need to disable $upgrade in BigInt, to avoid deep recursion
2095    local $Math::BigInt::upgrade = undef;    # should be really parent class vs MBI
2096  
2097    my $i = $MBI->_copy( $x->{_m} );
2098    $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2099    my $xas = Math::BigInt->bzero();
2100    $xas->{value} = $i;
2101  
2102    my $gs = $xas->copy()->bsqrt();    # some guess
2103  
2104    if (($x->{_es} ne '-')        # guess can't be accurate if there are
2105                      # digits after the dot
2106     && ($xas->bacmp($gs * $gs) == 0))    # guess hit the nail on the head?
2107      {
2108      # exact result, copy result over to keep $x
2109      $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
2110      $x->bnorm();
2111      # shortcut to not run through _find_round_parameters again
2112      if (defined $params[0])
2113        {
2114        $x->bround($params[0],$params[2]);    # then round accordingly
2115        }
2116      else
2117        {
2118        $x->bfround($params[1],$params[2]);    # then round accordingly
2119        }
2120      if ($fallback)
2121        {
2122        # clear a/p after round, since user did not request it
2123        delete $x->{_a}; delete $x->{_p};
2124        }
2125      # re-enable A and P, upgrade is taken care of by "local"
2126      ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
2127      return $x;
2128      }
2129   
2130    # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
2131    # of the result by multipyling the input by 100 and then divide the integer
2132    # result of sqrt(input) by 10. Rounding afterwards returns the real result.
2133  
2134    # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
2135    my $y1 = $MBI->_copy($x->{_m});
2136  
2137    my $length = $MBI->_len($y1);
2138    
2139    # Now calculate how many digits the result of sqrt(y1) would have
2140    my $digits = int($length / 2);
2141  
2142    # But we need at least $scale digits, so calculate how many are missing
2143    my $shift = $scale - $digits;
2144  
2145    # That should never happen (we take care of integer guesses above)
2146    # $shift = 0 if $shift < 0; 
2147  
2148    # Multiply in steps of 100, by shifting left two times the "missing" digits
2149    my $s2 = $shift * 2;
2150  
2151    # We now make sure that $y1 has the same odd or even number of digits than
2152    # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
2153    # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
2154    # steps of 10. The length of $x does not count, since an even or odd number
2155    # of digits before the dot is not changed by adding an even number of digits
2156    # after the dot (the result is still odd or even digits long).
2157    $s2++ if $MBI->_is_odd($x->{_e});
2158  
2159    $MBI->_lsft( $y1, $MBI->_new($s2), 10);
2160  
2161    # now take the square root and truncate to integer
2162    $y1 = $MBI->_sqrt($y1);
2163  
2164    # By "shifting" $y1 right (by creating a negative _e) we calculate the final
2165    # result, which is than later rounded to the desired scale.
2166  
2167    # calculate how many zeros $x had after the '.' (or before it, depending
2168    # on sign of $dat, the result should have half as many:
2169    my $dat = $MBI->_num($x->{_e});
2170    $dat = -$dat if $x->{_es} eq '-';
2171    $dat += $length;
2172  
2173    if ($dat > 0)
2174      {
2175      # no zeros after the dot (e.g. 1.23, 0.49 etc)
2176      # preserve half as many digits before the dot than the input had 
2177      # (but round this "up")
2178      $dat = int(($dat+1)/2);
2179      }
2180    else
2181      {
2182      $dat = int(($dat)/2);
2183      }
2184    $dat -= $MBI->_len($y1);
2185    if ($dat < 0)
2186      {
2187      $dat = abs($dat);
2188      $x->{_e} = $MBI->_new( $dat );
2189      $x->{_es} = '-';
2190      }
2191    else
2192      {    
2193      $x->{_e} = $MBI->_new( $dat );
2194      $x->{_es} = '+';
2195      }
2196    $x->{_m} = $y1;
2197    $x->bnorm();
2198  
2199    # shortcut to not run through _find_round_parameters again
2200    if (defined $params[0])
2201      {
2202      $x->bround($params[0],$params[2]);        # then round accordingly
2203      }
2204    else
2205      {
2206      $x->bfround($params[1],$params[2]);        # then round accordingly
2207      }
2208    if ($fallback)
2209      {
2210      # clear a/p after round, since user did not request it
2211      delete $x->{_a}; delete $x->{_p};
2212      }
2213    # restore globals
2214    $$abr = $ab; $$pbr = $pb;
2215    $x;
2216    }
2217  
2218  sub bfac
2219    {
2220    # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2221    # compute factorial number, modifies first argument
2222  
2223    # set up parameters
2224    my ($self,$x,@r) = (ref($_[0]),@_);
2225    # objectify is costly, so avoid it
2226    ($self,$x,@r) = objectify(1,@_) if !ref($x);
2227  
2228    # inf => inf
2229    return $x if $x->modify('bfac') || $x->{sign} eq '+inf';    
2230  
2231    return $x->bnan() 
2232      if (($x->{sign} ne '+') ||        # inf, NaN, <0 etc => NaN
2233       ($x->{_es} ne '+'));        # digits after dot?
2234  
2235    # use BigInt's bfac() for faster calc
2236    if (! $MBI->_is_zero($x->{_e}))
2237      {
2238      $MBI->_lsft($x->{_m}, $x->{_e},10);    # change 12e1 to 120e0
2239      $x->{_e} = $MBI->_zero();        # normalize
2240      $x->{_es} = '+';
2241      }
2242    $MBI->_fac($x->{_m});            # calculate factorial
2243    $x->bnorm()->round(@r);         # norm again and round result
2244    }
2245  
2246  sub _pow
2247    {
2248    # Calculate a power where $y is a non-integer, like 2 ** 0.3
2249    my ($x,$y,@r) = @_;
2250    my $self = ref($x);
2251  
2252    # if $y == 0.5, it is sqrt($x)
2253    $HALF = $self->new($HALF) unless ref($HALF);
2254    return $x->bsqrt(@r,$y) if $y->bcmp($HALF) == 0;
2255  
2256    # Using:
2257    # a ** x == e ** (x * ln a)
2258  
2259    # u = y * ln x
2260    #                _                         _
2261    # Taylor:       |   u    u^2    u^3         |
2262    # x ** y  = 1 + |  --- + --- + ----- + ...  |
2263    #               |_  1    1*2   1*2*3       _|
2264  
2265    # we need to limit the accuracy to protect against overflow
2266    my $fallback = 0;
2267    my ($scale,@params);
2268    ($x,@params) = $x->_find_round_parameters(@r);
2269      
2270    return $x if $x->is_nan();        # error in _find_round_parameters?
2271  
2272    # no rounding at all, so must use fallback
2273    if (scalar @params == 0)
2274      {
2275      # simulate old behaviour
2276      $params[0] = $self->div_scale();    # and round to it as accuracy
2277      $params[1] = undef;            # disable P
2278      $scale = $params[0]+4;         # at least four more for proper round
2279      $params[2] = $r[2];            # round mode by caller or undef
2280      $fallback = 1;            # to clear a/p afterwards
2281      }
2282    else
2283      {
2284      # the 4 below is empirical, and there might be cases where it is not
2285      # enough...
2286      $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2287      }
2288  
2289    # when user set globals, they would interfere with our calculation, so
2290    # disable them and later re-enable them
2291    no strict 'refs';
2292    my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2293    my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2294    # we also need to disable any set A or P on $x (_find_round_parameters took
2295    # them already into account), since these would interfere, too
2296    delete $x->{_a}; delete $x->{_p};
2297    # need to disable $upgrade in BigInt, to avoid deep recursion
2298    local $Math::BigInt::upgrade = undef;
2299   
2300    my ($limit,$v,$u,$below,$factor,$next,$over);
2301  
2302    $u = $x->copy()->blog(undef,$scale)->bmul($y);
2303    $v = $self->bone();                # 1
2304    $factor = $self->new(2);            # 2
2305    $x->bone();                    # first term: 1
2306  
2307    $below = $v->copy();
2308    $over = $u->copy();
2309  
2310    $limit = $self->new("1E-". ($scale-1));
2311    #my $steps = 0;
2312    while (3 < 5)
2313      {
2314      # we calculate the next term, and add it to the last
2315      # when the next term is below our limit, it won't affect the outcome
2316      # anymore, so we stop:
2317      $next = $over->copy()->bdiv($below,$scale);
2318      last if $next->bacmp($limit) <= 0;
2319      $x->badd($next);
2320      # calculate things for the next term
2321      $over *= $u; $below *= $factor; $factor->binc();
2322  
2323      last if $x->{sign} !~ /^[-+]$/;
2324  
2325      #$steps++;
2326      }
2327    
2328    # shortcut to not run through _find_round_parameters again
2329    if (defined $params[0])
2330      {
2331      $x->bround($params[0],$params[2]);        # then round accordingly
2332      }
2333    else
2334      {
2335      $x->bfround($params[1],$params[2]);        # then round accordingly
2336      }
2337    if ($fallback)
2338      {
2339      # clear a/p after round, since user did not request it
2340      delete $x->{_a}; delete $x->{_p};
2341      }
2342    # restore globals
2343    $$abr = $ab; $$pbr = $pb;
2344    $x;
2345    }
2346  
2347  sub bpow 
2348    {
2349    # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2350    # compute power of two numbers, second arg is used as integer
2351    # modifies first argument
2352  
2353    # set up parameters
2354    my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2355    # objectify is costly, so avoid it
2356    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2357      {
2358      ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2359      }
2360  
2361    return $x if $x->modify('bpow');
2362  
2363    return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
2364    return $x if $x->{sign} =~ /^[+-]inf$/;
2365    
2366    # cache the result of is_zero
2367    my $y_is_zero = $y->is_zero();
2368    return $x->bone() if $y_is_zero;
2369    return $x         if $x->is_one() || $y->is_one();
2370  
2371    my $x_is_zero = $x->is_zero();
2372    return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int();        # non-integer power
2373  
2374    my $y1 = $y->as_number()->{value};            # make MBI part
2375  
2376    # if ($x == -1)
2377    if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
2378      {
2379      # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
2380      return $MBI->_is_odd($y1) ? $x : $x->babs(1);
2381      }
2382    if ($x_is_zero)
2383      {
2384      return $x if $y->{sign} eq '+';     # 0**y => 0 (if not y <= 0)
2385      # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
2386      return $x->binf();
2387      }
2388  
2389    my $new_sign = '+';
2390    $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
2391  
2392    # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
2393    $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
2394    $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
2395  
2396    $x->{sign} = $new_sign;
2397    $x->bnorm();
2398    if ($y->{sign} eq '-')
2399      {
2400      # modify $x in place!
2401      my $z = $x->copy(); $x->bone();
2402      return scalar $x->bdiv($z,$a,$p,$r);    # round in one go (might ignore y's A!)
2403      }
2404    $x->round($a,$p,$r,$y);
2405    }
2406  
2407  sub bmodpow
2408    {
2409    # takes a very large number to a very large exponent in a given very
2410    # large modulus, quickly, thanks to binary exponentation. Supports
2411    # negative exponents.
2412    my ($self,$num,$exp,$mod,@r) = objectify(3,@_);
2413  
2414    return $num if $num->modify('bmodpow');
2415  
2416    # check modulus for valid values
2417    return $num->bnan() if ($mod->{sign} ne '+'           # NaN, - , -inf, +inf
2418                         || $mod->is_zero());
2419  
2420    # check exponent for valid values
2421    if ($exp->{sign} =~ /\w/)
2422      {
2423      # i.e., if it's NaN, +inf, or -inf...
2424      return $num->bnan();
2425      }
2426  
2427    $num->bmodinv ($mod) if ($exp->{sign} eq '-');
2428  
2429    # check num for valid values (also NaN if there was no inverse but $exp < 0)
2430    return $num->bnan() if $num->{sign} !~ /^[+-]$/;
2431  
2432    # $mod is positive, sign on $exp is ignored, result also positive
2433  
2434    # XXX TODO: speed it up when all three numbers are integers
2435    $num->bpow($exp)->bmod($mod);
2436    }
2437  
2438  ###############################################################################
2439  # trigonometric functions
2440  
2441  # helper function for bpi() and batan2(), calculates arcus tanges (1/x)
2442  
2443  sub _atan_inv
2444    {
2445    # return a/b so that a/b approximates atan(1/x) to at least limit digits
2446    my ($self, $x, $limit) = @_;
2447  
2448    # Taylor:       x^3   x^5   x^7   x^9
2449    #    atan = x - --- + --- - --- + --- - ...
2450    #                3     5     7     9 
2451  
2452    #               1      1         1        1
2453    #    atan 1/x = - - ------- + ------- - ------- + ...
2454    #               x   x^3 * 3   x^5 * 5   x^7 * 7 
2455  
2456    #               1      1         1            1
2457    #    atan 1/x = - - --------- + ---------- - ----------- + ... 
2458    #               5    3 * 125     5 * 3125     7 * 78125
2459  
2460    # Subtraction/addition of a rational:
2461  
2462    #  5    7    5*3 +- 7*4
2463    #  - +- -  = ----------
2464    #  4    3       4*3
2465  
2466    # Term:  N        N+1
2467    #
2468    #        a             1                  a * d * c +- b
2469    #        ----- +- ------------------  =  ----------------
2470    #        b           d * c                b * d * c
2471  
2472    #  since b1 = b0 * (d-2) * c
2473  
2474    #        a             1                  a * d +- b / c
2475    #        ----- +- ------------------  =  ----------------
2476    #        b           d * c                b * d 
2477  
2478    # and  d = d + 2
2479    # and  c = c * x * x
2480  
2481    #        u = d * c
2482    #        stop if length($u) > limit 
2483    #        a = a * u +- b
2484    #        b = b * u
2485    #        d = d + 2
2486    #        c = c * x * x
2487    #        sign = 1 - sign
2488  
2489    my $a = $MBI->_one();
2490    my $b = $MBI->_copy($x);
2491   
2492    my $x2  = $MBI->_mul( $MBI->_copy($x), $b);        # x2 = x * x
2493    my $d   = $MBI->_new( 3 );                # d = 3
2494    my $c   = $MBI->_mul( $MBI->_copy($x), $x2);        # c = x ^ 3
2495    my $two = $MBI->_new( 2 );
2496  
2497    # run the first step unconditionally
2498    my $u = $MBI->_mul( $MBI->_copy($d), $c);
2499    $a = $MBI->_mul($a, $u);
2500    $a = $MBI->_sub($a, $b);
2501    $b = $MBI->_mul($b, $u);
2502    $d = $MBI->_add($d, $two);
2503    $c = $MBI->_mul($c, $x2);
2504  
2505    # a is now a * (d-3) * c
2506    # b is now b * (d-2) * c
2507  
2508    # run the second step unconditionally
2509    $u = $MBI->_mul( $MBI->_copy($d), $c);
2510    $a = $MBI->_mul($a, $u);
2511    $a = $MBI->_add($a, $b);
2512    $b = $MBI->_mul($b, $u);
2513    $d = $MBI->_add($d, $two);
2514    $c = $MBI->_mul($c, $x2);
2515  
2516    # a is now a * (d-3) * (d-5) * c * c  
2517    # b is now b * (d-2) * (d-4) * c * c
2518  
2519    # so we can remove c * c from both a and b to shorten the numbers involved:
2520    $a = $MBI->_div($a, $x2);
2521    $b = $MBI->_div($b, $x2);
2522    $a = $MBI->_div($a, $x2);
2523    $b = $MBI->_div($b, $x2);
2524  
2525  #  my $step = 0; 
2526    my $sign = 0;                        # 0 => -, 1 => +
2527    while (3 < 5)
2528      {
2529  #    $step++;
2530  #    if (($i++ % 100) == 0)
2531  #      {
2532  #    print "a=",$MBI->_str($a),"\n";
2533  #    print "b=",$MBI->_str($b),"\n";
2534  #      }
2535  #    print "d=",$MBI->_str($d),"\n";
2536  #    print "x2=",$MBI->_str($x2),"\n";
2537  #    print "c=",$MBI->_str($c),"\n";
2538  
2539      my $u = $MBI->_mul( $MBI->_copy($d), $c);
2540      # use _alen() for libs like GMP where _len() would be O(N^2)
2541      last if $MBI->_alen($u) > $limit;
2542      my ($bc,$r) = $MBI->_div( $MBI->_copy($b), $c);
2543      if ($MBI->_is_zero($r))
2544        {
2545        # b / c is an integer, so we can remove c from all terms
2546        # this happens almost every time:
2547        $a = $MBI->_mul($a, $d);
2548        $a = $MBI->_sub($a, $bc) if $sign == 0;
2549        $a = $MBI->_add($a, $bc) if $sign == 1;
2550        $b = $MBI->_mul($b, $d);
2551        }
2552      else
2553        {
2554        # b / c is not an integer, so we keep c in the terms
2555        # this happens very rarely, for instance for x = 5, this happens only
2556        # at the following steps:
2557        # 1, 5, 14, 32, 72, 157, 340, ...
2558        $a = $MBI->_mul($a, $u);
2559        $a = $MBI->_sub($a, $b) if $sign == 0;
2560        $a = $MBI->_add($a, $b) if $sign == 1;
2561        $b = $MBI->_mul($b, $u);
2562        }
2563      $d = $MBI->_add($d, $two);
2564      $c = $MBI->_mul($c, $x2);
2565      $sign = 1 - $sign;
2566  
2567      }
2568  
2569  #  print "Took $step steps for ", $MBI->_str($x),"\n";
2570  #  print "a=",$MBI->_str($a),"\n"; print "b=",$MBI->_str($b),"\n";
2571    # return a/b so that a/b approximates atan(1/x)
2572    ($a,$b);
2573    }
2574  
2575  sub bpi
2576    {
2577    my ($self,$n) = @_;
2578    if (@_ == 0)
2579      {
2580      $self = $class;
2581      }
2582    if (@_ == 1)
2583      {
2584      # called like Math::BigFloat::bpi(10);
2585      $n = $self; $self = $class;
2586      # called like Math::BigFloat->bpi();
2587      $n = undef if $n eq 'Math::BigFloat';
2588      }
2589    $self = ref($self) if ref($self);
2590    my $fallback = defined $n ? 0 : 1;
2591    $n = 40 if !defined $n || $n < 1;
2592  
2593    # after 黃見利 (Hwang Chien-Lih) (1997)
2594    # pi/4 = 183 * atan(1/239) + 32 * atan(1/1023) – 68 * atan(1/5832)
2595    #     + 12 * atan(1/110443) - 12 * atan(1/4841182) - 100 * atan(1/6826318)
2596  
2597    # a few more to prevent rounding errors
2598    $n += 4;
2599  
2600    my ($a,$b) = $self->_atan_inv( $MBI->_new(239),$n);
2601    my ($c,$d) = $self->_atan_inv( $MBI->_new(1023),$n);
2602    my ($e,$f) = $self->_atan_inv( $MBI->_new(5832),$n);
2603    my ($g,$h) = $self->_atan_inv( $MBI->_new(110443),$n);
2604    my ($i,$j) = $self->_atan_inv( $MBI->_new(4841182),$n);
2605    my ($k,$l) = $self->_atan_inv( $MBI->_new(6826318),$n);
2606  
2607    $MBI->_mul($a, $MBI->_new(732));
2608    $MBI->_mul($c, $MBI->_new(128));
2609    $MBI->_mul($e, $MBI->_new(272));
2610    $MBI->_mul($g, $MBI->_new(48));
2611    $MBI->_mul($i, $MBI->_new(48));
2612    $MBI->_mul($k, $MBI->_new(400));
2613  
2614    my $x = $self->bone(); $x->{_m} = $a; my $x_d = $self->bone(); $x_d->{_m} = $b;
2615    my $y = $self->bone(); $y->{_m} = $c; my $y_d = $self->bone(); $y_d->{_m} = $d;
2616    my $z = $self->bone(); $z->{_m} = $e; my $z_d = $self->bone(); $z_d->{_m} = $f;
2617    my $u = $self->bone(); $u->{_m} = $g; my $u_d = $self->bone(); $u_d->{_m} = $h;
2618    my $v = $self->bone(); $v->{_m} = $i; my $v_d = $self->bone(); $v_d->{_m} = $j;
2619    my $w = $self->bone(); $w->{_m} = $k; my $w_d = $self->bone(); $w_d->{_m} = $l;
2620    $x->bdiv($x_d, $n);
2621    $y->bdiv($y_d, $n);
2622    $z->bdiv($z_d, $n);
2623    $u->bdiv($u_d, $n);
2624    $v->bdiv($v_d, $n);
2625    $w->bdiv($w_d, $n);
2626  
2627    delete $x->{_a}; delete $y->{_a}; delete $z->{_a};
2628    delete $u->{_a}; delete $v->{_a}; delete $w->{_a};
2629    $x->badd($y)->bsub($z)->badd($u)->bsub($v)->bsub($w);
2630  
2631    $x->bround($n-4);
2632    delete $x->{_a} if $fallback == 1;
2633    $x;
2634    }
2635  
2636  sub bcos
2637    {
2638    # Calculate a cosinus of x.
2639    my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2640  
2641    # Taylor:      x^2   x^4   x^6   x^8
2642    #    cos = 1 - --- + --- - --- + --- ...
2643    #               2!    4!    6!    8!
2644  
2645    # we need to limit the accuracy to protect against overflow
2646    my $fallback = 0;
2647    my ($scale,@params);
2648    ($x,@params) = $x->_find_round_parameters(@r);
2649      
2650    #         constant object       or error in _find_round_parameters?
2651    return $x if $x->modify('bcos') || $x->is_nan();
2652  
2653    return $x->bone(@r) if $x->is_zero();
2654  
2655    # no rounding at all, so must use fallback
2656    if (scalar @params == 0)
2657      {
2658      # simulate old behaviour
2659      $params[0] = $self->div_scale();    # and round to it as accuracy
2660      $params[1] = undef;            # disable P
2661      $scale = $params[0]+4;         # at least four more for proper round
2662      $params[2] = $r[2];            # round mode by caller or undef
2663      $fallback = 1;            # to clear a/p afterwards
2664      }
2665    else
2666      {
2667      # the 4 below is empirical, and there might be cases where it is not
2668      # enough...
2669      $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2670      }
2671  
2672    # when user set globals, they would interfere with our calculation, so
2673    # disable them and later re-enable them
2674    no strict 'refs';
2675    my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2676    my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2677    # we also need to disable any set A or P on $x (_find_round_parameters took
2678    # them already into account), since these would interfere, too
2679    delete $x->{_a}; delete $x->{_p};
2680    # need to disable $upgrade in BigInt, to avoid deep recursion
2681    local $Math::BigInt::upgrade = undef;
2682   
2683    my $last = 0;
2684    my $over = $x * $x;                   # X ^ 2
2685    my $x2 = $over->copy();               # X ^ 2; difference between terms
2686    my $sign = 1;                         # start with -=
2687    my $below = $self->new(2); my $factorial = $self->new(3);
2688    $x->bone(); delete $x->{_a}; delete $x->{_p};
2689  
2690    my $limit = $self->new("1E-". ($scale-1));
2691    #my $steps = 0;
2692    while (3 < 5)
2693      {
2694      # we calculate the next term, and add it to the last
2695      # when the next term is below our limit, it won't affect the outcome
2696      # anymore, so we stop:
2697      my $next = $over->copy()->bdiv($below,$scale);
2698      last if $next->bacmp($limit) <= 0;
2699  
2700      if ($sign == 0)
2701        {
2702        $x->badd($next);
2703        }
2704      else
2705        {
2706        $x->bsub($next);
2707        }
2708      $sign = 1-$sign;                    # alternate
2709      # calculate things for the next term
2710      $over->bmul($x2);                    # $x*$x
2711      $below->bmul($factorial); $factorial->binc();    # n*(n+1)
2712      $below->bmul($factorial); $factorial->binc();    # n*(n+1)
2713      }
2714  
2715    # shortcut to not run through _find_round_parameters again
2716    if (defined $params[0])
2717      {
2718      $x->bround($params[0],$params[2]);        # then round accordingly
2719      }
2720    else
2721      {
2722      $x->bfround($params[1],$params[2]);        # then round accordingly
2723      }
2724    if ($fallback)
2725      {
2726      # clear a/p after round, since user did not request it
2727      delete $x->{_a}; delete $x->{_p};
2728      }
2729    # restore globals
2730    $$abr = $ab; $$pbr = $pb;
2731    $x;
2732    }
2733  
2734  sub bsin
2735    {
2736    # Calculate a sinus of x.
2737    my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2738  
2739    # taylor:      x^3   x^5   x^7   x^9
2740    #    sin = x - --- + --- - --- + --- ...
2741    #               3!    5!    7!    9!
2742  
2743    # we need to limit the accuracy to protect against overflow
2744    my $fallback = 0;
2745    my ($scale,@params);
2746    ($x,@params) = $x->_find_round_parameters(@r);
2747      
2748    #         constant object       or error in _find_round_parameters?
2749    return $x if $x->modify('bsin') || $x->is_nan();
2750  
2751    return $x->bzero(@r) if $x->is_zero();
2752  
2753    # no rounding at all, so must use fallback
2754    if (scalar @params == 0)
2755      {
2756      # simulate old behaviour
2757      $params[0] = $self->div_scale();    # and round to it as accuracy
2758      $params[1] = undef;            # disable P
2759      $scale = $params[0]+4;         # at least four more for proper round
2760      $params[2] = $r[2];            # round mode by caller or undef
2761      $fallback = 1;            # to clear a/p afterwards
2762      }
2763    else
2764      {
2765      # the 4 below is empirical, and there might be cases where it is not
2766      # enough...
2767      $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2768      }
2769  
2770    # when user set globals, they would interfere with our calculation, so
2771    # disable them and later re-enable them
2772    no strict 'refs';
2773    my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2774    my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2775    # we also need to disable any set A or P on $x (_find_round_parameters took
2776    # them already into account), since these would interfere, too
2777    delete $x->{_a}; delete $x->{_p};
2778    # need to disable $upgrade in BigInt, to avoid deep recursion
2779    local $Math::BigInt::upgrade = undef;
2780   
2781    my $last = 0;
2782    my $over = $x * $x;            # X ^ 2
2783    my $x2 = $over->copy();        # X ^ 2; difference between terms
2784    $over->bmul($x);            # X ^ 3 as starting value
2785    my $sign = 1;                # start with -=
2786    my $below = $self->new(6); my $factorial = $self->new(4);
2787    delete $x->{_a}; delete $x->{_p};
2788  
2789    my $limit = $self->new("1E-". ($scale-1));
2790    #my $steps = 0;
2791    while (3 < 5)
2792      {
2793      # we calculate the next term, and add it to the last
2794      # when the next term is below our limit, it won't affect the outcome
2795      # anymore, so we stop:
2796      my $next = $over->copy()->bdiv($below,$scale);
2797      last if $next->bacmp($limit) <= 0;
2798  
2799      if ($sign == 0)
2800        {
2801        $x->badd($next);
2802        }
2803      else
2804        {
2805        $x->bsub($next);
2806        }
2807      $sign = 1-$sign;                    # alternate
2808      # calculate things for the next term
2809      $over->bmul($x2);                    # $x*$x
2810      $below->bmul($factorial); $factorial->binc();    # n*(n+1)
2811      $below->bmul($factorial); $factorial->binc();    # n*(n+1)
2812      }
2813  
2814    # shortcut to not run through _find_round_parameters again
2815    if (defined $params[0])
2816      {
2817      $x->bround($params[0],$params[2]);        # then round accordingly
2818      }
2819    else
2820      {
2821      $x->bfround($params[1],$params[2]);        # then round accordingly
2822      }
2823    if ($fallback)
2824      {
2825      # clear a/p after round, since user did not request it
2826      delete $x->{_a}; delete $x->{_p};
2827      }
2828    # restore globals
2829    $$abr = $ab; $$pbr = $pb;
2830    $x;
2831    }
2832  
2833  sub batan2
2834    { 
2835    # calculate arcus tangens of ($y/$x)
2836    
2837    # set up parameters
2838    my ($self,$y,$x,@r) = (ref($_[0]),@_);
2839    # objectify is costly, so avoid it
2840    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2841      {
2842      ($self,$y,$x,@r) = objectify(2,@_);
2843      }
2844  
2845    return $y if $y->modify('batan2');
2846  
2847    return $y->bnan() if ($y->{sign} eq $nan) || ($x->{sign} eq $nan);
2848  
2849    return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
2850  
2851    # Y X
2852    # 0 0 result is 0
2853    # 0 +x result is 0
2854    return $y->bzero(@r) if $y->is_zero() && $x->{sign} eq '+';
2855  
2856    # Y X
2857    # 0 -x result is PI
2858    if ($y->is_zero())
2859      {
2860      # calculate PI
2861      my $pi = $self->bpi(@r);
2862      # modify $x in place
2863      $y->{_m} = $pi->{_m};
2864      $y->{_e} = $pi->{_e};
2865      $y->{_es} = $pi->{_es};
2866      $y->{sign} = '+';
2867      return $y;
2868      }
2869  
2870    # Y X
2871    # +y 0 result is PI/2
2872    # -y 0 result is -PI/2
2873    if ($y->is_inf() || $x->is_zero())
2874      {
2875      # calculate PI/2
2876      my $pi = $self->bpi(@r);
2877      # modify $x in place
2878      $y->{_m} = $pi->{_m};
2879      $y->{_e} = $pi->{_e};
2880      $y->{_es} = $pi->{_es};
2881      # -y => -PI/2, +y => PI/2
2882      $y->{sign} = substr($y->{sign},0,1);        # +inf => +
2883      $MBI->_div($y->{_m}, $MBI->_new(2));
2884      return $y;
2885      }
2886  
2887    # we need to limit the accuracy to protect against overflow
2888    my $fallback = 0;
2889    my ($scale,@params);
2890    ($y,@params) = $y->_find_round_parameters(@r);
2891      
2892    # error in _find_round_parameters?
2893    return $y if $y->is_nan();
2894  
2895    # no rounding at all, so must use fallback
2896    if (scalar @params == 0)
2897      {
2898      # simulate old behaviour
2899      $params[0] = $self->div_scale();    # and round to it as accuracy
2900      $params[1] = undef;            # disable P
2901      $scale = $params[0]+4;         # at least four more for proper round
2902      $params[2] = $r[2];            # round mode by caller or undef
2903      $fallback = 1;            # to clear a/p afterwards
2904      }
2905    else
2906      {
2907      # the 4 below is empirical, and there might be cases where it is not
2908      # enough...
2909      $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2910      }
2911  
2912    # inlined is_one() && is_one('-')
2913    if ($MBI->_is_one($y->{_m}) && $MBI->_is_zero($y->{_e}))
2914      {
2915      # shortcut: 1 1 result is PI/4
2916      # inlined is_one() && is_one('-')
2917      if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
2918        {
2919        # 1,1 => PI/4
2920        my $pi_4 = $self->bpi( $scale - 3);
2921        # modify $x in place
2922        $y->{_m} = $pi_4->{_m};
2923        $y->{_e} = $pi_4->{_e};
2924        $y->{_es} = $pi_4->{_es};
2925        # 1 1 => +
2926        # -1 1 => -
2927        # 1 -1 => -
2928        # -1 -1 => +
2929        $y->{sign} = $x->{sign} eq $y->{sign} ? '+' : '-';
2930        $MBI->_div($y->{_m}, $MBI->_new(4));
2931        return $y;
2932        }
2933      # shortcut: 1 int(X) result is _atan_inv(X)
2934  
2935      # is integer
2936      if ($x->{_es} eq '+')
2937        {
2938        my $x1 = $MBI->_copy($x->{_m});
2939        $MBI->_lsft($x1, $x->{_e},10) unless $MBI->_is_zero($x->{_e});
2940  
2941        my ($a,$b) = $self->_atan_inv($x1, $scale);
2942        my $y_sign = $y->{sign};
2943        # calculate A/B
2944        $y->bone(); $y->{_m} = $a; my $y_d = $self->bone(); $y_d->{_m} = $b;
2945        $y->bdiv($y_d, @r);
2946        $y->{sign} = $y_sign;
2947        return $y;
2948        }
2949      }
2950  
2951    # handle all other cases
2952    #  X  Y
2953    # +x +y 0 to PI/2
2954    # -x +y PI/2 to PI
2955    # +x -y 0 to -PI/2
2956    # -x -y -PI/2 to -PI 
2957  
2958    my $y_sign = $y->{sign};
2959  
2960    # divide $x by $y
2961    $y->bdiv($x, $scale) unless $x->is_one();
2962    $y->batan(@r);
2963  
2964    # restore sign
2965    $y->{sign} = $y_sign;
2966  
2967    $y;
2968    }
2969  
2970  sub batan
2971    {
2972    # Calculate a arcus tangens of x.
2973    my ($x,@r) = @_;
2974    my $self = ref($x);
2975  
2976    # taylor:       x^3   x^5   x^7   x^9
2977    #    atan = x - --- + --- - --- + --- ...
2978    #                3     5     7     9 
2979  
2980    # we need to limit the accuracy to protect against overflow
2981    my $fallback = 0;
2982    my ($scale,@params);
2983    ($x,@params) = $x->_find_round_parameters(@r);
2984      
2985    #         constant object       or error in _find_round_parameters?
2986    return $x if $x->modify('batan') || $x->is_nan();
2987  
2988    if ($x->{sign} =~ /^[+-]inf\z/)
2989      {
2990      # +inf result is PI/2
2991      # -inf result is -PI/2
2992      # calculate PI/2
2993      my $pi = $self->bpi(@r);
2994      # modify $x in place
2995      $x->{_m} = $pi->{_m};
2996      $x->{_e} = $pi->{_e};
2997      $x->{_es} = $pi->{_es};
2998      # -y => -PI/2, +y => PI/2
2999      $x->{sign} = substr($x->{sign},0,1);        # +inf => +
3000      $MBI->_div($x->{_m}, $MBI->_new(2));
3001      return $x;
3002      }
3003  
3004    return $x->bzero(@r) if $x->is_zero();
3005  
3006    # no rounding at all, so must use fallback
3007    if (scalar @params == 0)
3008      {
3009      # simulate old behaviour
3010      $params[0] = $self->div_scale();    # and round to it as accuracy
3011      $params[1] = undef;            # disable P
3012      $scale = $params[0]+4;         # at least four more for proper round
3013      $params[2] = $r[2];            # round mode by caller or undef
3014      $fallback = 1;            # to clear a/p afterwards
3015      }
3016    else
3017      {
3018      # the 4 below is empirical, and there might be cases where it is not
3019      # enough...
3020      $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3021      }
3022  
3023    # 1 or -1 => PI/4
3024    # inlined is_one() && is_one('-')
3025    if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
3026      {
3027      my $pi = $self->bpi($scale - 3);
3028      # modify $x in place
3029      $x->{_m} = $pi->{_m};
3030      $x->{_e} = $pi->{_e};
3031      $x->{_es} = $pi->{_es};
3032      # leave the sign of $x alone (+1 => +PI/4, -1 => -PI/4)
3033      $MBI->_div($x->{_m}, $MBI->_new(4));
3034      return $x;
3035      }
3036    
3037    # This series is only valid if -1 < x < 1, so for other x we need to
3038    # to calculate PI/2 - atan(1/x):
3039    my $one = $MBI->_new(1);
3040    my $pi = undef;
3041    if ($x->{_es} eq '+' && ($MBI->_acmp($x->{_m},$one) >= 0))
3042      {
3043      # calculate PI/2
3044      $pi = $self->bpi($scale - 3);
3045      $MBI->_div($pi->{_m}, $MBI->_new(2));
3046      # calculate 1/$x:
3047      my $x_copy = $x->copy();
3048      # modify $x in place
3049      $x->bone(); $x->bdiv($x_copy,$scale);
3050      }
3051  
3052    # when user set globals, they would interfere with our calculation, so
3053    # disable them and later re-enable them
3054    no strict 'refs';
3055    my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
3056    my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
3057    # we also need to disable any set A or P on $x (_find_round_parameters took
3058    # them already into account), since these would interfere, too
3059    delete $x->{_a}; delete $x->{_p};
3060    # need to disable $upgrade in BigInt, to avoid deep recursion
3061    local $Math::BigInt::upgrade = undef;
3062   
3063    my $last = 0;
3064    my $over = $x * $x;            # X ^ 2
3065    my $x2 = $over->copy();        # X ^ 2; difference between terms
3066    $over->bmul($x);            # X ^ 3 as starting value
3067    my $sign = 1;                # start with -=
3068    my $below = $self->new(3);
3069    my $two = $self->new(2);
3070    delete $x->{_a}; delete $x->{_p};
3071  
3072    my $limit = $self->new("1E-". ($scale-1));
3073    #my $steps = 0;
3074    while (3 < 5)
3075      {
3076      # we calculate the next term, and add it to the last
3077      # when the next term is below our limit, it won't affect the outcome
3078      # anymore, so we stop:
3079      my $next = $over->copy()->bdiv($below,$scale);
3080      last if $next->bacmp($limit) <= 0;
3081  
3082      if ($sign == 0)
3083        {
3084        $x->badd($next);
3085        }
3086      else
3087        {
3088        $x->bsub($next);
3089        }
3090      $sign = 1-$sign;                    # alternate
3091      # calculate things for the next term
3092      $over->bmul($x2);                    # $x*$x
3093      $below->badd($two);                    # n += 2
3094      }
3095  
3096    if (defined $pi)
3097      {
3098      my $x_copy = $x->copy();
3099      # modify $x in place
3100      $x->{_m} = $pi->{_m};
3101      $x->{_e} = $pi->{_e};
3102      $x->{_es} = $pi->{_es};
3103      # PI/2 - $x
3104      $x->bsub($x_copy);
3105      }
3106  
3107    # shortcut to not run through _find_round_parameters again
3108    if (defined $params[0])
3109      {
3110      $x->bround($params[0],$params[2]);        # then round accordingly
3111      }
3112    else
3113      {
3114      $x->bfround($params[1],$params[2]);        # then round accordingly
3115      }
3116    if ($fallback)
3117      {
3118      # clear a/p after round, since user did not request it
3119      delete $x->{_a}; delete $x->{_p};
3120      }
3121    # restore globals
3122    $$abr = $ab; $$pbr = $pb;
3123    $x;
3124    }
3125  
3126  ###############################################################################
3127  # rounding functions
3128  
3129  sub bfround
3130    {
3131    # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
3132    # $n == 0 means round to integer
3133    # expects and returns normalized numbers!
3134    my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
3135  
3136    my ($scale,$mode) = $x->_scale_p(@_);
3137    return $x if !defined $scale || $x->modify('bfround'); # no-op
3138  
3139    # never round a 0, +-inf, NaN
3140    if ($x->is_zero())
3141      {
3142      $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
3143      return $x; 
3144      }
3145    return $x if $x->{sign} !~ /^[+-]$/;
3146  
3147    # don't round if x already has lower precision
3148    return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
3149  
3150    $x->{_p} = $scale;            # remember round in any case
3151    delete $x->{_a};            # and clear A
3152    if ($scale < 0)
3153      {
3154      # round right from the '.'
3155  
3156      return $x if $x->{_es} eq '+';        # e >= 0 => nothing to round
3157  
3158      $scale = -$scale;                # positive for simplicity
3159      my $len = $MBI->_len($x->{_m});        # length of mantissa
3160  
3161      # the following poses a restriction on _e, but if _e is bigger than a
3162      # scalar, you got other problems (memory etc) anyway
3163      my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e})));    # digits after dot
3164      my $zad = 0;                # zeros after dot
3165      $zad = $dad - $len if (-$dad < -$len);    # for 0.00..00xxx style
3166     
3167      # p rint "scale $scale dad $dad zad $zad len $len\n";
3168      # number  bsstr   len zad dad    
3169      # 0.123   123e-3    3   0 3
3170      # 0.0123  123e-4    3   1 4
3171      # 0.001   1e-3      1   2 3
3172      # 1.23    123e-2    3   0 2
3173      # 1.2345  12345e-4    5   0 4
3174  
3175      # do not round after/right of the $dad
3176      return $x if $scale > $dad;            # 0.123, scale >= 3 => exit
3177  
3178      # round to zero if rounding inside the $zad, but not for last zero like:
3179      # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
3180      return $x->bzero() if $scale < $zad;
3181      if ($scale == $zad)            # for 0.006, scale -3 and trunc
3182        {
3183        $scale = -$len;
3184        }
3185      else
3186        {
3187        # adjust round-point to be inside mantissa
3188        if ($zad != 0)
3189          {
3190      $scale = $scale-$zad;
3191          }
3192        else
3193          {
3194          my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;    # digits before dot
3195      $scale = $dbd+$scale;
3196          }
3197        }
3198      }
3199    else
3200      {
3201      # round left from the '.'
3202  
3203      # 123 => 100 means length(123) = 3 - $scale (2) => 1
3204  
3205      my $dbt = $MBI->_len($x->{_m}); 
3206      # digits before dot 
3207      my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
3208      # should be the same, so treat it as this 
3209      $scale = 1 if $scale == 0; 
3210      # shortcut if already integer 
3211      return $x if $scale == 1 && $dbt <= $dbd; 
3212      # maximum digits before dot 
3213      ++$dbd;
3214  
3215      if ($scale > $dbd) 
3216         { 
3217         # not enough digits before dot, so round to zero 
3218         return $x->bzero; 
3219         }
3220      elsif ( $scale == $dbd )
3221         { 
3222         # maximum 
3223         $scale = -$dbt; 
3224         } 
3225      else
3226         { 
3227         $scale = $dbd - $scale; 
3228         }
3229      }
3230    # pass sign to bround for rounding modes '+inf' and '-inf'
3231    my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
3232    $m->bround($scale,$mode);
3233    $x->{_m} = $m->{value};            # get our mantissa back
3234    $x->bnorm();
3235    }
3236  
3237  sub bround
3238    {
3239    # accuracy: preserve $N digits, and overwrite the rest with 0's
3240    my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
3241  
3242    if (($_[0] || 0) < 0)
3243      {
3244      require Carp; Carp::croak ('bround() needs positive accuracy');
3245      }
3246  
3247    my ($scale,$mode) = $x->_scale_a(@_);
3248    return $x if !defined $scale || $x->modify('bround');    # no-op
3249  
3250    # scale is now either $x->{_a}, $accuracy, or the user parameter
3251    # test whether $x already has lower accuracy, do nothing in this case 
3252    # but do round if the accuracy is the same, since a math operation might
3253    # want to round a number with A=5 to 5 digits afterwards again
3254    return $x if defined $x->{_a} && $x->{_a} < $scale;
3255  
3256    # scale < 0 makes no sense
3257    # scale == 0 => keep all digits
3258    # never round a +-inf, NaN
3259    return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
3260  
3261    # 1: never round a 0
3262    # 2: if we should keep more digits than the mantissa has, do nothing
3263    if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
3264      {
3265      $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
3266      return $x; 
3267      }
3268  
3269    # pass sign to bround for '+inf' and '-inf' rounding modes
3270    my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
3271  
3272    $m->bround($scale,$mode);        # round mantissa
3273    $x->{_m} = $m->{value};        # get our mantissa back
3274    $x->{_a} = $scale;            # remember rounding
3275    delete $x->{_p};            # and clear P
3276    $x->bnorm();                # del trailing zeros gen. by bround()
3277    }
3278  
3279  sub bfloor
3280    {
3281    # return integer less or equal then $x
3282    my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3283  
3284    return $x if $x->modify('bfloor');
3285     
3286    return $x if $x->{sign} !~ /^[+-]$/;    # nan, +inf, -inf
3287  
3288    # if $x has digits after dot
3289    if ($x->{_es} eq '-')
3290      {
3291      $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3292      $x->{_e} = $MBI->_zero();            # trunc/norm    
3293      $x->{_es} = '+';                # abs e
3294      $MBI->_inc($x->{_m}) if $x->{sign} eq '-';    # increment if negative
3295      }
3296    $x->round($a,$p,$r);
3297    }
3298  
3299  sub bceil
3300    {
3301    # return integer greater or equal then $x
3302    my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3303  
3304    return $x if $x->modify('bceil');
3305    return $x if $x->{sign} !~ /^[+-]$/;    # nan, +inf, -inf
3306  
3307    # if $x has digits after dot
3308    if ($x->{_es} eq '-')
3309      {
3310      $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3311      $x->{_e} = $MBI->_zero();            # trunc/norm    
3312      $x->{_es} = '+';                # abs e
3313      $MBI->_inc($x->{_m}) if $x->{sign} eq '+';    # increment if positive
3314      }
3315    $x->round($a,$p,$r);
3316    }
3317  
3318  sub brsft
3319    {
3320    # shift right by $y (divide by power of $n)
3321    
3322    # set up parameters
3323    my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3324    # objectify is costly, so avoid it
3325    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3326      {
3327      ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3328      }
3329  
3330    return $x if $x->modify('brsft');
3331    return $x if $x->{sign} !~ /^[+-]$/;    # nan, +inf, -inf
3332  
3333    $n = 2 if !defined $n; $n = $self->new($n);
3334  
3335    # negative amount?
3336    return $x->blsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3337  
3338    # the following call to bdiv() will return either quo or (quo,reminder):
3339    $x->bdiv($n->bpow($y),$a,$p,$r,$y);
3340    }
3341  
3342  sub blsft
3343    {
3344    # shift left by $y (multiply by power of $n)
3345    
3346    # set up parameters
3347    my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3348    # objectify is costly, so avoid it
3349    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3350      {
3351      ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3352      }
3353  
3354    return $x if $x->modify('blsft');
3355    return $x if $x->{sign} !~ /^[+-]$/;    # nan, +inf, -inf
3356  
3357    $n = 2 if !defined $n; $n = $self->new($n);
3358  
3359    # negative amount?
3360    return $x->brsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3361  
3362    $x->bmul($n->bpow($y),$a,$p,$r,$y);
3363    }
3364  
3365  ###############################################################################
3366  
3367  sub DESTROY
3368    {
3369    # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
3370    }
3371  
3372  sub AUTOLOAD
3373    {
3374    # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
3375    # or falling back to MBI::bxxx()
3376    my $name = $AUTOLOAD;
3377  
3378    $name =~ s/(.*):://;    # split package
3379    my $c = $1 || $class;
3380    no strict 'refs';
3381    $c->import() if $IMPORT == 0;
3382    if (!_method_alias($name))
3383      {
3384      if (!defined $name)
3385        {
3386        # delayed load of Carp and avoid recursion    
3387        require Carp;
3388        Carp::croak ("$c: Can't call a method without name");
3389        }
3390      if (!_method_hand_up($name))
3391        {
3392        # delayed load of Carp and avoid recursion    
3393        require Carp;
3394        Carp::croak ("Can't call $c\-\>$name, not a valid method");
3395        }
3396      # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
3397      $name =~ s/^f/b/;
3398      return &{"Math::BigInt"."::$name"}(@_);
3399      }
3400    my $bname = $name; $bname =~ s/^f/b/;
3401    $c .= "::$name";
3402    *{$c} = \&{$bname};
3403    &{$c};    # uses @_
3404    }
3405  
3406  sub exponent
3407    {
3408    # return a copy of the exponent
3409    my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3410  
3411    if ($x->{sign} !~ /^[+-]$/)
3412      {
3413      my $s = $x->{sign}; $s =~ s/^[+-]//;
3414      return Math::BigInt->new($s);         # -inf, +inf => +inf
3415      }
3416    Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
3417    }
3418  
3419  sub mantissa
3420    {
3421    # return a copy of the mantissa
3422    my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3423   
3424    if ($x->{sign} !~ /^[+-]$/)
3425      {
3426      my $s = $x->{sign}; $s =~ s/^[+]//;
3427      return Math::BigInt->new($s);        # -inf, +inf => +inf
3428      }
3429    my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
3430    $m->bneg() if $x->{sign} eq '-';
3431  
3432    $m;
3433    }
3434  
3435  sub parts
3436    {
3437    # return a copy of both the exponent and the mantissa
3438    my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3439  
3440    if ($x->{sign} !~ /^[+-]$/)
3441      {
3442      my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
3443      return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
3444      }
3445    my $m = Math::BigInt->bzero();
3446    $m->{value} = $MBI->_copy($x->{_m});
3447    $m->bneg() if $x->{sign} eq '-';
3448    ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
3449    }
3450  
3451  ##############################################################################
3452  # private stuff (internal use only)
3453  
3454  sub import
3455    {
3456    my $self = shift;
3457    my $l = scalar @_;
3458    my $lib = ''; my @a;
3459    my $lib_kind = 'try';
3460    $IMPORT=1;
3461    for ( my $i = 0; $i < $l ; $i++)
3462      {
3463      if ( $_[$i] eq ':constant' )
3464        {
3465        # This causes overlord er load to step in. 'binary' and 'integer'
3466        # are handled by BigInt.
3467        overload::constant float => sub { $self->new(shift); }; 
3468        }
3469      elsif ($_[$i] eq 'upgrade')
3470        {
3471        # this causes upgrading
3472        $upgrade = $_[$i+1];        # or undef to disable
3473        $i++;
3474        }
3475      elsif ($_[$i] eq 'downgrade')
3476        {
3477        # this causes downgrading
3478        $downgrade = $_[$i+1];        # or undef to disable
3479        $i++;
3480        }
3481      elsif ($_[$i] =~ /^(lib|try|only)\z/)
3482        {
3483        # alternative library
3484        $lib = $_[$i+1] || '';        # default Calc
3485        $lib_kind = $1;            # lib, try or only
3486        $i++;
3487        }
3488      elsif ($_[$i] eq 'with')
3489        {
3490        # alternative class for our private parts()
3491        # XXX: no longer supported
3492        # $MBI = $_[$i+1] || 'Math::BigInt';
3493        $i++;
3494        }
3495      else
3496        {
3497        push @a, $_[$i];
3498        }
3499      }
3500  
3501    $lib =~ tr/a-zA-Z0-9,://cd;        # restrict to sane characters
3502    # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
3503    my $mbilib = eval { Math::BigInt->config()->{lib} };
3504    if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
3505      {
3506      # MBI already loaded
3507      Math::BigInt->import( $lib_kind, "$lib,$mbilib", 'objectify');
3508      }
3509    else
3510      {
3511      # MBI not loaded, or with ne "Math::BigInt::Calc"
3512      $lib .= ",$mbilib" if defined $mbilib;
3513      $lib =~ s/^,//;                # don't leave empty 
3514      
3515      # replacement library can handle lib statement, but also could ignore it
3516      
3517      # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
3518      # used in the same script, or eval inside import(). So we require MBI:
3519      require Math::BigInt;
3520      Math::BigInt->import( $lib_kind => $lib, 'objectify' );
3521      }
3522    if ($@)
3523      {
3524      require Carp; Carp::croak ("Couldn't load $lib: $! $@");
3525      }
3526    # find out which one was actually loaded
3527    $MBI = Math::BigInt->config()->{lib};
3528  
3529    # register us with MBI to get notified of future lib changes
3530    Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
3531  
3532    $self->export_to_level(1,$self,@a);        # export wanted functions
3533    }
3534  
3535  sub bnorm
3536    {
3537    # adjust m and e so that m is smallest possible
3538    my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
3539  
3540    return $x if $x->{sign} !~ /^[+-]$/;        # inf, nan etc
3541  
3542    my $zeros = $MBI->_zeros($x->{_m});        # correct for trailing zeros
3543    if ($zeros != 0)
3544      {
3545      my $z = $MBI->_new($zeros);
3546      $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
3547      if ($x->{_es} eq '-')
3548        {
3549        if ($MBI->_acmp($x->{_e},$z) >= 0)
3550          {
3551          $x->{_e} = $MBI->_sub ($x->{_e}, $z);
3552          $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
3553          }
3554        else
3555          {
3556          $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
3557          $x->{_es} = '+';
3558          }
3559        }
3560      else
3561        {
3562        $x->{_e} = $MBI->_add ($x->{_e}, $z);
3563        }
3564      }
3565    else
3566      {
3567      # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
3568      # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
3569      $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
3570       if $MBI->_is_zero($x->{_m});
3571      }
3572  
3573    $x;                    # MBI bnorm is no-op, so dont call it
3574    } 
3575   
3576  ##############################################################################
3577  
3578  sub as_hex
3579    {
3580    # return number as hexadecimal string (only for integers defined)
3581    my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3582  
3583    return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
3584    return '0x0' if $x->is_zero();
3585  
3586    return $nan if $x->{_es} ne '+';        # how to do 1e-1 in hex!?
3587  
3588    my $z = $MBI->_copy($x->{_m});
3589    if (! $MBI->_is_zero($x->{_e}))        # > 0 
3590      {
3591      $MBI->_lsft( $z, $x->{_e},10);
3592      }
3593    $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3594    $z->as_hex();
3595    }
3596  
3597  sub as_bin
3598    {
3599    # return number as binary digit string (only for integers defined)
3600    my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3601  
3602    return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
3603    return '0b0' if $x->is_zero();
3604  
3605    return $nan if $x->{_es} ne '+';        # how to do 1e-1 in hex!?
3606  
3607    my $z = $MBI->_copy($x->{_m});
3608    if (! $MBI->_is_zero($x->{_e}))        # > 0 
3609      {
3610      $MBI->_lsft( $z, $x->{_e},10);
3611      }
3612    $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3613    $z->as_bin();
3614    }
3615  
3616  sub as_oct
3617    {
3618    # return number as octal digit string (only for integers defined)
3619    my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3620  
3621    return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
3622    return '0' if $x->is_zero();
3623  
3624    return $nan if $x->{_es} ne '+';        # how to do 1e-1 in hex!?
3625  
3626    my $z = $MBI->_copy($x->{_m});
3627    if (! $MBI->_is_zero($x->{_e}))        # > 0 
3628      {
3629      $MBI->_lsft( $z, $x->{_e},10);
3630      }
3631    $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3632    $z->as_oct();
3633    }
3634  
3635  sub as_number
3636    {
3637    # return copy as a bigint representation of this BigFloat number
3638    my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3639  
3640    return $x if $x->modify('as_number');
3641  
3642    my $z = $MBI->_copy($x->{_m});
3643    if ($x->{_es} eq '-')            # < 0
3644      {
3645      $MBI->_rsft( $z, $x->{_e},10);
3646      } 
3647    elsif (! $MBI->_is_zero($x->{_e}))    # > 0 
3648      {
3649      $MBI->_lsft( $z, $x->{_e},10);
3650      }
3651    $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3652    $z;
3653    }
3654  
3655  sub length
3656    {
3657    my $x = shift;
3658    my $class = ref($x) || $x;
3659    $x = $class->new(shift) unless ref($x);
3660  
3661    return 1 if $MBI->_is_zero($x->{_m});
3662  
3663    my $len = $MBI->_len($x->{_m});
3664    $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
3665    if (wantarray())
3666      {
3667      my $t = 0;
3668      $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
3669      return ($len, $t);
3670      }
3671    $len;
3672    }
3673  
3674  1;
3675  __END__
3676  
3677  =head1 NAME
3678  
3679  Math::BigFloat - Arbitrary size floating point math package
3680  
3681  =head1 SYNOPSIS
3682  
3683    use Math::BigFloat;
3684  
3685    # Number creation
3686    my $x = Math::BigFloat->new($str);    # defaults to 0
3687    my $y = $x->copy();            # make a true copy
3688    my $nan  = Math::BigFloat->bnan();    # create a NotANumber
3689    my $zero = Math::BigFloat->bzero();    # create a +0
3690    my $inf = Math::BigFloat->binf();    # create a +inf
3691    my $inf = Math::BigFloat->binf('-');    # create a -inf
3692    my $one = Math::BigFloat->bone();    # create a +1
3693    my $mone = Math::BigFloat->bone('-');    # create a -1
3694  
3695    my $pi = Math::BigFloat->bpi(100);    # PI to 100 digits
3696  
3697    # the following examples compute their result to 100 digits accuracy:
3698    my $cos  = Math::BigFloat->new(1)->bcos(100);        # cosinus(1)
3699    my $sin  = Math::BigFloat->new(1)->bsin(100);        # sinus(1)
3700    my $atan = Math::BigFloat->new(1)->batan(100);    # arcus tangens(1)
3701  
3702    my $atan2 = Math::BigFloat->new(  1 )->batan2( 1 ,100); # batan(1)
3703    my $atan2 = Math::BigFloat->new(  1 )->batan2( 8 ,100); # batan(1/8)
3704    my $atan2 = Math::BigFloat->new( -2 )->batan2( 1 ,100); # batan(-2)
3705  
3706    # Testing
3707    $x->is_zero();        # true if arg is +0
3708    $x->is_nan();            # true if arg is NaN
3709    $x->is_one();            # true if arg is +1
3710    $x->is_one('-');        # true if arg is -1
3711    $x->is_odd();            # true if odd, false for even
3712    $x->is_even();        # true if even, false for odd
3713    $x->is_pos();            # true if >= 0
3714    $x->is_neg();            # true if <  0
3715    $x->is_inf(sign);        # true if +inf, or -inf (default is '+')
3716  
3717    $x->bcmp($y);            # compare numbers (undef,<0,=0,>0)
3718    $x->bacmp($y);        # compare absolutely (undef,<0,=0,>0)
3719    $x->sign();            # return the sign, either +,- or NaN
3720    $x->digit($n);        # return the nth digit, counting from right
3721    $x->digit(-$n);        # return the nth digit, counting from left 
3722  
3723    # The following all modify their first argument. If you want to preserve
3724    # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
3725    # necessary when mixing $a = $b assignments with non-overloaded math.
3726   
3727    # set 
3728    $x->bzero();            # set $i to 0
3729    $x->bnan();            # set $i to NaN
3730    $x->bone();                   # set $x to +1
3731    $x->bone('-');                # set $x to -1
3732    $x->binf();                   # set $x to inf
3733    $x->binf('-');                # set $x to -inf
3734  
3735    $x->bneg();            # negation
3736    $x->babs();            # absolute value
3737    $x->bnorm();            # normalize (no-op)
3738    $x->bnot();            # two's complement (bit wise not)
3739    $x->binc();            # increment x by 1
3740    $x->bdec();            # decrement x by 1
3741    
3742    $x->badd($y);            # addition (add $y to $x)
3743    $x->bsub($y);            # subtraction (subtract $y from $x)
3744    $x->bmul($y);            # multiplication (multiply $x by $y)
3745    $x->bdiv($y);            # divide, set $x to quotient
3746                  # return (quo,rem) or quo if scalar
3747  
3748    $x->bmod($y);            # modulus ($x % $y)
3749    $x->bpow($y);            # power of arguments ($x ** $y)
3750    $x->bmodpow($exp,$mod);    # modular exponentation (($num**$exp) % $mod))
3751    $x->blsft($y, $n);        # left shift by $y places in base $n
3752    $x->brsft($y, $n);        # right shift by $y places in base $n
3753                  # returns (quo,rem) or quo if in scalar context
3754    
3755    $x->blog();            # logarithm of $x to base e (Euler's number)
3756    $x->blog($base);        # logarithm of $x to base $base (f.i. 2)
3757    $x->bexp();            # calculate e ** $x where e is Euler's number
3758    
3759    $x->band($y);            # bit-wise and
3760    $x->bior($y);            # bit-wise inclusive or
3761    $x->bxor($y);            # bit-wise exclusive or
3762    $x->bnot();            # bit-wise not (two's complement)
3763   
3764    $x->bsqrt();            # calculate square-root
3765    $x->broot($y);        # $y'th root of $x (e.g. $y == 3 => cubic root)
3766    $x->bfac();            # factorial of $x (1*2*3*4*..$x)
3767   
3768    $x->bround($N);         # accuracy: preserve $N digits
3769    $x->bfround($N);        # precision: round to the $Nth digit
3770  
3771    $x->bfloor();            # return integer less or equal than $x
3772    $x->bceil();            # return integer greater or equal than $x
3773  
3774    # The following do not modify their arguments:
3775  
3776    bgcd(@values);        # greatest common divisor
3777    blcm(@values);        # lowest common multiplicator
3778    
3779    $x->bstr();            # return string
3780    $x->bsstr();            # return string in scientific notation
3781  
3782    $x->as_int();            # return $x as BigInt 
3783    $x->exponent();        # return exponent as BigInt
3784    $x->mantissa();        # return mantissa as BigInt
3785    $x->parts();            # return (mantissa,exponent) as BigInt
3786  
3787    $x->length();            # number of digits (w/o sign and '.')
3788    ($l,$f) = $x->length();    # number of digits, and length of fraction    
3789  
3790    $x->precision();        # return P of $x (or global, if P of $x undef)
3791    $x->precision($n);        # set P of $x to $n
3792    $x->accuracy();        # return A of $x (or global, if A of $x undef)
3793    $x->accuracy($n);        # set A $x to $n
3794  
3795    # these get/set the appropriate global value for all BigFloat objects
3796    Math::BigFloat->precision();    # Precision
3797    Math::BigFloat->accuracy();    # Accuracy
3798    Math::BigFloat->round_mode();    # rounding mode
3799  
3800  =head1 DESCRIPTION
3801  
3802  All operators (including basic math operations) are overloaded if you
3803  declare your big floating point numbers as
3804  
3805    $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
3806  
3807  Operations with overloaded operators preserve the arguments, which is
3808  exactly what you expect.
3809  
3810  =head2 Canonical notation
3811  
3812  Input to these routines are either BigFloat objects, or strings of the
3813  following four forms:
3814  
3815  =over 2
3816  
3817  =item *
3818  
3819  C</^[+-]\d+$/>
3820  
3821  =item *
3822  
3823  C</^[+-]\d+\.\d*$/>
3824  
3825  =item *
3826  
3827  C</^[+-]\d+E[+-]?\d+$/>
3828  
3829  =item *
3830  
3831  C</^[+-]\d*\.\d+E[+-]?\d+$/>
3832  
3833  =back
3834  
3835  all with optional leading and trailing zeros and/or spaces. Additionally,
3836  numbers are allowed to have an underscore between any two digits.
3837  
3838  Empty strings as well as other illegal numbers results in 'NaN'.
3839  
3840  bnorm() on a BigFloat object is now effectively a no-op, since the numbers 
3841  are always stored in normalized form. On a string, it creates a BigFloat 
3842  object.
3843  
3844  =head2 Output
3845  
3846  Output values are BigFloat objects (normalized), except for bstr() and bsstr().
3847  
3848  The string output will always have leading and trailing zeros stripped and drop
3849  a plus sign. C<bstr()> will give you always the form with a decimal point,
3850  while C<bsstr()> (s for scientific) gives you the scientific notation.
3851  
3852      Input            bstr()        bsstr()
3853      '-0'            '0'        '0E1'
3854         '  -123 123 123'    '-123123123'    '-123123123E0'
3855      '00.0123'        '0.0123'    '123E-4'
3856      '123.45E-2'        '1.2345'    '12345E-4'
3857      '10E+3'            '10000'        '1E4'
3858  
3859  Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
3860  C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
3861  return either undef, <0, 0 or >0 and are suited for sort.
3862  
3863  Actual math is done by using the class defined with C<with => Class;> (which
3864  defaults to BigInts) to represent the mantissa and exponent.
3865  
3866  The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to 
3867  represent the result when input arguments are not numbers, as well as 
3868  the result of dividing by zero.
3869  
3870  =head2 C<mantissa()>, C<exponent()> and C<parts()>
3871  
3872  C<mantissa()> and C<exponent()> return the said parts of the BigFloat 
3873  as BigInts such that:
3874  
3875      $m = $x->mantissa();
3876      $e = $x->exponent();
3877      $y = $m * ( 10 ** $e );
3878      print "ok\n" if $x == $y;
3879  
3880  C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
3881  
3882  A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
3883  
3884  Currently the mantissa is reduced as much as possible, favouring higher
3885  exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
3886  This might change in the future, so do not depend on it.
3887  
3888  =head2 Accuracy vs. Precision
3889  
3890  See also: L<Rounding|Rounding>.
3891  
3892  Math::BigFloat supports both precision (rounding to a certain place before or
3893  after the dot) and accuracy (rounding to a certain number of digits). For a
3894  full documentation, examples and tips on these topics please see the large
3895  section about rounding in L<Math::BigInt>.
3896  
3897  Since things like C<sqrt(2)> or C<1 / 3> must presented with a limited
3898  accuracy lest a operation consumes all resources, each operation produces
3899  no more than the requested number of digits.
3900  
3901  If there is no gloabl precision or accuracy set, B<and> the operation in
3902  question was not called with a requested precision or accuracy, B<and> the
3903  input $x has no accuracy or precision set, then a fallback parameter will
3904  be used. For historical reasons, it is called C<div_scale> and can be accessed
3905  via:
3906  
3907      $d = Math::BigFloat->div_scale();        # query
3908      Math::BigFloat->div_scale($n);            # set to $n digits
3909  
3910  The default value for C<div_scale> is 40.
3911  
3912  In case the result of one operation has more digits than specified,
3913  it is rounded. The rounding mode taken is either the default mode, or the one
3914  supplied to the operation after the I<scale>:
3915  
3916      $x = Math::BigFloat->new(2);
3917      Math::BigFloat->accuracy(5);        # 5 digits max
3918      $y = $x->copy()->bdiv(3);        # will give 0.66667
3919      $y = $x->copy()->bdiv(3,6);        # will give 0.666667
3920      $y = $x->copy()->bdiv(3,6,undef,'odd');    # will give 0.666667
3921      Math::BigFloat->round_mode('zero');
3922      $y = $x->copy()->bdiv(3,6);        # will also give 0.666667
3923  
3924  Note that C<< Math::BigFloat->accuracy() >> and C<< Math::BigFloat->precision() >>
3925  set the global variables, and thus B<any> newly created number will be subject
3926  to the global rounding B<immediately>. This means that in the examples above, the
3927  C<3> as argument to C<bdiv()> will also get an accuracy of B<5>.
3928  
3929  It is less confusing to either calculate the result fully, and afterwards
3930  round it explicitly, or use the additional parameters to the math
3931  functions like so:
3932  
3933      use Math::BigFloat;    
3934      $x = Math::BigFloat->new(2);
3935      $y = $x->copy()->bdiv(3);
3936      print $y->bround(5),"\n";        # will give 0.66667
3937  
3938      or
3939  
3940      use Math::BigFloat;    
3941      $x = Math::BigFloat->new(2);
3942      $y = $x->copy()->bdiv(3,5);        # will give 0.66667
3943      print "$y\n";
3944  
3945  =head2 Rounding
3946  
3947  =over 2
3948  
3949  =item ffround ( +$scale )
3950  
3951  Rounds to the $scale'th place left from the '.', counting from the dot.
3952  The first digit is numbered 1. 
3953  
3954  =item ffround ( -$scale )
3955  
3956  Rounds to the $scale'th place right from the '.', counting from the dot.
3957  
3958  =item ffround ( 0 )
3959  
3960  Rounds to an integer.
3961  
3962  =item fround  ( +$scale )
3963  
3964  Preserves accuracy to $scale digits from the left (aka significant digits)
3965  and pads the rest with zeros. If the number is between 1 and -1, the
3966  significant digits count from the first non-zero after the '.'
3967  
3968  =item fround  ( -$scale ) and fround ( 0 )
3969  
3970  These are effectively no-ops.
3971  
3972  =back
3973  
3974  All rounding functions take as a second parameter a rounding mode from one of
3975  the following: 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'.
3976  
3977  The default rounding mode is 'even'. By using
3978  C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
3979  mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
3980  no longer supported.
3981  The second parameter to the round functions then overrides the default
3982  temporarily. 
3983  
3984  The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
3985  'trunc' as rounding mode to make it equivalent to:
3986  
3987      $x = 2.5;
3988      $y = int($x) + 2;
3989  
3990  You can override this by passing the desired rounding mode as parameter to
3991  C<as_number()>:
3992  
3993      $x = Math::BigFloat->new(2.5);
3994      $y = $x->as_number('odd');    # $y = 3
3995  
3996  =head1 METHODS
3997  
3998  Math::BigFloat supports all methods that Math::BigInt supports, except it
3999  calculates non-integer results when possible. Please see L<Math::BigInt>
4000  for a full description of each method. Below are just the most important
4001  differences:
4002  
4003  =head2 accuracy
4004  
4005          $x->accuracy(5);                # local for $x
4006          CLASS->accuracy(5);             # global for all members of CLASS
4007                                          # Note: This also applies to new()!
4008  
4009          $A = $x->accuracy();            # read out accuracy that affects $x
4010          $A = CLASS->accuracy();         # read out global accuracy
4011  
4012  Set or get the global or local accuracy, aka how many significant digits the
4013  results have. If you set a global accuracy, then this also applies to new()!
4014  
4015  Warning! The accuracy I<sticks>, e.g. once you created a number under the
4016  influence of C<< CLASS->accuracy($A) >>, all results from math operations with
4017  that number will also be rounded.
4018  
4019  In most cases, you should probably round the results explicitly using one of
4020  L<round()>, L<bround()> or L<bfround()> or by passing the desired accuracy
4021  to the math operation as additional parameter:
4022  
4023          my $x = Math::BigInt->new(30000);
4024          my $y = Math::BigInt->new(7);
4025          print scalar $x->copy()->bdiv($y, 2);           # print 4300
4026          print scalar $x->copy()->bdiv($y)->bround(2);   # print 4300
4027  
4028  =head2 precision()
4029  
4030          $x->precision(-2);      # local for $x, round at the second digit right of the dot
4031          $x->precision(2);       # ditto, round at the second digit left of the dot
4032  
4033          CLASS->precision(5);    # Global for all members of CLASS
4034                                  # This also applies to new()!
4035          CLASS->precision(-5);   # ditto
4036  
4037          $P = CLASS->precision();        # read out global precision
4038          $P = $x->precision();           # read out precision that affects $x
4039  
4040  Note: You probably want to use L<accuracy()> instead. With L<accuracy> you
4041  set the number of digits each result should have, with L<precision> you
4042  set the place where to round!
4043  
4044  =head2 bexp()
4045  
4046      $x->bexp($accuracy);        # calculate e ** X
4047  
4048  Calculates the expression C<e ** $x> where C<e> is Euler's number.
4049  
4050  This method was added in v1.82 of Math::BigInt (April 2007).
4051  
4052  =head2 bnok()
4053  
4054      $x->bnok($y);           # x over y (binomial coefficient n over k)
4055  
4056  Calculates the binomial coefficient n over k, also called the "choose"
4057  function. The result is equivalent to:
4058  
4059      ( n )      n!
4060      | - |  = -------
4061      ( k )    k!(n-k)!
4062  
4063  This method was added in v1.84 of Math::BigInt (April 2007).
4064  
4065  =head2 bpi()
4066  
4067      print Math::BigFloat->bpi(100), "\n";
4068  
4069  Calculate PI to N digits (including the 3 before the dot). The result is
4070  rounded according to the current rounding mode, which defaults to "even".
4071  
4072  This method was added in v1.87 of Math::BigInt (June 2007).
4073  
4074  =head2 bcos()
4075  
4076      my $x = Math::BigFloat->new(1);
4077      print $x->bcos(100), "\n";
4078  
4079  Calculate the cosinus of $x, modifying $x in place.
4080  
4081  This method was added in v1.87 of Math::BigInt (June 2007).
4082  
4083  =head2 bsin()
4084  
4085      my $x = Math::BigFloat->new(1);
4086      print $x->bsin(100), "\n";
4087  
4088  Calculate the sinus of $x, modifying $x in place.
4089  
4090  This method was added in v1.87 of Math::BigInt (June 2007).
4091  
4092  =head2 batan2()
4093  
4094      my $y = Math::BigFloat->new(2);
4095      my $x = Math::BigFloat->new(3);
4096      print $y->batan2($x), "\n";
4097  
4098  Calculate the arcus tanges of C<$y> divided by C<$x>, modifying $y in place.
4099  See also L<batan()>.
4100  
4101  This method was added in v1.87 of Math::BigInt (June 2007).
4102  
4103  =head2 batan()
4104  
4105      my $x = Math::BigFloat->new(1);
4106      print $x->batan(100), "\n";
4107  
4108  Calculate the arcus tanges of $x, modifying $x in place. See also L<batan2()>.
4109  
4110  This method was added in v1.87 of Math::BigInt (June 2007).
4111  
4112  =head2 bmuladd()
4113  
4114      $x->bmuladd($y,$z);        
4115  
4116  Multiply $x by $y, and then add $z to the result.
4117  
4118  This method was added in v1.87 of Math::BigInt (June 2007).
4119  
4120  =head1 Autocreating constants
4121  
4122  After C<use Math::BigFloat ':constant'> all the floating point constants
4123  in the given scope are converted to C<Math::BigFloat>. This conversion
4124  happens at compile time.
4125  
4126  In particular
4127  
4128    perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
4129  
4130  prints the value of C<2E-100>. Note that without conversion of 
4131  constants the expression 2E-100 will be calculated as normal floating point 
4132  number.
4133  
4134  Please note that ':constant' does not affect integer constants, nor binary 
4135  nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
4136  work.
4137  
4138  =head2 Math library
4139  
4140  Math with the numbers is done (by default) by a module called
4141  Math::BigInt::Calc. This is equivalent to saying:
4142  
4143      use Math::BigFloat lib => 'Calc';
4144  
4145  You can change this by using:
4146  
4147      use Math::BigFloat lib => 'GMP';
4148  
4149  Note: The keyword 'lib' will warn when the requested library could not be
4150  loaded. To suppress the warning use 'try' instead:
4151  
4152      use Math::BigFloat try => 'GMP';
4153  
4154  To turn the warning into a die(), use 'only' instead:
4155  
4156      use Math::BigFloat only => 'GMP';
4157  
4158  The following would first try to find Math::BigInt::Foo, then
4159  Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
4160  
4161      use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
4162  
4163  See the respective low-level library documentation for further details.
4164  
4165  Please note that Math::BigFloat does B<not> use the denoted library itself,
4166  but it merely passes the lib argument to Math::BigInt. So, instead of the need
4167  to do:
4168  
4169      use Math::BigInt lib => 'GMP';
4170      use Math::BigFloat;
4171  
4172  you can roll it all into one line:
4173  
4174      use Math::BigFloat lib => 'GMP';
4175  
4176  It is also possible to just require Math::BigFloat:
4177  
4178      require Math::BigFloat;
4179  
4180  This will load the necessary things (like BigInt) when they are needed, and
4181  automatically.
4182  
4183  See L<Math::BigInt> for more details than you ever wanted to know about using
4184  a different low-level library.
4185  
4186  =head2 Using Math::BigInt::Lite
4187  
4188  For backwards compatibility reasons it is still possible to
4189  request a different storage class for use with Math::BigFloat:
4190  
4191          use Math::BigFloat with => 'Math::BigInt::Lite';
4192  
4193  However, this request is ignored, as the current code now uses the low-level
4194  math libary for directly storing the number parts.
4195  
4196  =head1 EXPORTS
4197  
4198  C<Math::BigFloat> exports nothing by default, but can export the C<bpi()> method:
4199  
4200      use Math::BigFloat qw/bpi/;
4201  
4202      print bpi(10), "\n";
4203  
4204  =head1 BUGS
4205  
4206  Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
4207  
4208  =head1 CAVEATS
4209  
4210  Do not try to be clever to insert some operations in between switching
4211  libraries:
4212  
4213      require Math::BigFloat;
4214      my $matter = Math::BigFloat->bone() + 4;    # load BigInt and Calc
4215      Math::BigFloat->import( lib => 'Pari' );    # load Pari, too
4216      my $anti_matter = Math::BigFloat->bone()+4;    # now use Pari
4217  
4218  This will create objects with numbers stored in two different backend libraries,
4219  and B<VERY BAD THINGS> will happen when you use these together:
4220  
4221      my $flash_and_bang = $matter + $anti_matter;    # Don't do this!
4222  
4223  =over 1
4224  
4225  =item stringify, bstr()
4226  
4227  Both stringify and bstr() now drop the leading '+'. The old code would return
4228  '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
4229  reasoning and details.
4230  
4231  =item bdiv
4232  
4233  The following will probably not print what you expect:
4234  
4235      print $c->bdiv(123.456),"\n";
4236  
4237  It prints both quotient and reminder since print works in list context. Also,
4238  bdiv() will modify $c, so be careful. You probably want to use
4239      
4240      print $c / 123.456,"\n";
4241      print scalar $c->bdiv(123.456),"\n";  # or if you want to modify $c
4242  
4243  instead.
4244  
4245  =item brsft
4246  
4247  The following will probably not print what you expect:
4248  
4249      my $c = Math::BigFloat->new('3.14159');
4250          print $c->brsft(3,10),"\n";    # prints 0.00314153.1415
4251  
4252  It prints both quotient and remainder, since print calls C<brsft()> in list
4253  context. Also, C<< $c->brsft() >> will modify $c, so be careful.
4254  You probably want to use
4255  
4256      print scalar $c->copy()->brsft(3,10),"\n";
4257      # or if you really want to modify $c
4258          print scalar $c->brsft(3,10),"\n";
4259  
4260  instead.
4261  
4262  =item Modifying and =
4263  
4264  Beware of:
4265  
4266      $x = Math::BigFloat->new(5);
4267      $y = $x;
4268  
4269  It will not do what you think, e.g. making a copy of $x. Instead it just makes
4270  a second reference to the B<same> object and stores it in $y. Thus anything
4271  that modifies $x will modify $y (except overloaded math operators), and vice
4272  versa. See L<Math::BigInt> for details and how to avoid that.
4273  
4274  =item bpow
4275  
4276  C<bpow()> now modifies the first argument, unlike the old code which left
4277  it alone and only returned the result. This is to be consistent with
4278  C<badd()> etc. The first will modify $x, the second one won't:
4279  
4280      print bpow($x,$i),"\n";     # modify $x
4281      print $x->bpow($i),"\n";     # ditto
4282      print $x ** $i,"\n";        # leave $x alone 
4283  
4284  =item precision() vs. accuracy()
4285  
4286  A common pitfall is to use L<precision()> when you want to round a result to
4287  a certain number of digits:
4288  
4289      use Math::BigFloat;
4290  
4291      Math::BigFloat->precision(4);        # does not do what you think it does
4292      my $x = Math::BigFloat->new(12345);    # rounds $x to "12000"!
4293      print "$x\n";                # print "12000"
4294      my $y = Math::BigFloat->new(3);        # rounds $y to "0"!
4295      print "$y\n";                # print "0"
4296      $z = $x / $y;                # 12000 / 0 => NaN!
4297      print "$z\n";
4298      print $z->precision(),"\n";        # 4
4299  
4300  Replacing L<precision> with L<accuracy> is probably not what you want, either:
4301  
4302      use Math::BigFloat;
4303  
4304      Math::BigFloat->accuracy(4);        # enables global rounding:
4305      my $x = Math::BigFloat->new(123456);    # rounded immediately to "12350"
4306      print "$x\n";                # print "123500"
4307      my $y = Math::BigFloat->new(3);        # rounded to "3
4308      print "$y\n";                # print "3"
4309      print $z = $x->copy()->bdiv($y),"\n";    # 41170
4310      print $z->accuracy(),"\n";        # 4
4311  
4312  What you want to use instead is:
4313  
4314      use Math::BigFloat;
4315  
4316      my $x = Math::BigFloat->new(123456);    # no rounding
4317      print "$x\n";                # print "123456"
4318      my $y = Math::BigFloat->new(3);        # no rounding
4319      print "$y\n";                # print "3"
4320      print $z = $x->copy()->bdiv($y,4),"\n";    # 41150
4321      print $z->accuracy(),"\n";        # undef
4322  
4323  In addition to computing what you expected, the last example also does B<not>
4324  "taint" the result with an accuracy or precision setting, which would
4325  influence any further operation.
4326  
4327  =back
4328  
4329  =head1 SEE ALSO
4330  
4331  L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
4332  L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and  L<Math::BigInt::GMP>.
4333  
4334  The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
4335  because they solve the autoupgrading/downgrading issue, at least partly.
4336  
4337  The package at L<http://search.cpan.org/~tels/Math-BigInt> contains
4338  more documentation including a full version history, testcases, empty
4339  subclass files and benchmarks.
4340  
4341  =head1 LICENSE
4342  
4343  This program is free software; you may redistribute it and/or modify it under
4344  the same terms as Perl itself.
4345  
4346  =head1 AUTHORS
4347  
4348  Mark Biggar, overloaded interface by Ilya Zakharevich.
4349  Completely rewritten by Tels L<http://bloodgate.com> in 2001 - 2006, and still
4350  at it in 2007.
4351  
4352  =cut