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

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