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.9994';
  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 do not 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   else
1502     {
1503       undef $l_2;
1504     }
1505
1506   $self->_log($x,$scale);                       # need to do the "normal" way
1507   $x->badd($l_10) if defined $l_10;             # correct it by ln(10)
1508   $x->badd($l_2) if defined $l_2;               # and maybe by ln(2)
1509
1510   # all done, $x contains now the result
1511   $x;
1512   }
1513
1514 sub blcm
1515   {
1516   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1517   # does not modify arguments, but returns new object
1518   # Lowest Common Multiplicator
1519
1520   my ($self,@arg) = objectify(0,@_);
1521   my $x = $self->new(shift @arg);
1522   while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); }
1523   $x;
1524   }
1525
1526 sub bgcd
1527   {
1528   # (BINT or num_str, BINT or num_str) return BINT
1529   # does not modify arguments, but returns new object
1530
1531   my $y = shift;
1532   $y = __PACKAGE__->new($y) if !ref($y);
1533   my $self = ref($y);
1534   my $x = $y->copy()->babs();                   # keep arguments
1535
1536   return $x->bnan() if $x->{sign} !~ /^[+-]$/   # x NaN?
1537         || !$x->is_int();                       # only for integers now
1538
1539   while (@_)
1540     {
1541     my $t = shift; $t = $self->new($t) if !ref($t);
1542     $y = $t->copy()->babs();
1543
1544     return $x->bnan() if $y->{sign} !~ /^[+-]$/ # y NaN?
1545         || !$y->is_int();                       # only for integers now
1546
1547     # greatest common divisor
1548     while (! $y->is_zero())
1549       {
1550       ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
1551       }
1552
1553     last if $x->is_one();
1554     }
1555   $x;
1556   }
1557
1558 ##############################################################################
1559
1560 sub _e_add
1561   {
1562   # Internal helper sub to take two positive integers and their signs and
1563   # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1564   # output ($CALC,('+'|'-'))
1565   my ($x,$y,$xs,$ys) = @_;
1566
1567   # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1568   if ($xs eq $ys)
1569     {
1570     $x = $MBI->_add ($x, $y );          # a+b
1571     # the sign follows $xs
1572     return ($x, $xs);
1573     }
1574
1575   my $a = $MBI->_acmp($x,$y);
1576   if ($a > 0)
1577     {
1578     $x = $MBI->_sub ($x , $y);                          # abs sub
1579     }
1580   elsif ($a == 0)
1581     {
1582     $x = $MBI->_zero();                                 # result is 0
1583     $xs = '+';
1584     }
1585   else # a < 0
1586     {
1587     $x = $MBI->_sub ( $y, $x, 1 );                      # abs sub
1588     $xs = $ys;
1589     }
1590   ($x,$xs);
1591   }
1592
1593 sub _e_sub
1594   {
1595   # Internal helper sub to take two positive integers and their signs and
1596   # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1597   # output ($CALC,('+'|'-'))
1598   my ($x,$y,$xs,$ys) = @_;
1599
1600   # flip sign
1601   $ys =~ tr/+-/-+/;
1602   _e_add($x,$y,$xs,$ys);                # call add (does subtract now)
1603   }
1604
1605 ###############################################################################
1606 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1607
1608 sub is_int
1609   {
1610   # return true if arg (BFLOAT or num_str) is an integer
1611   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1612
1613   (($x->{sign} =~ /^[+-]$/) &&                  # NaN and +-inf aren't
1614    ($x->{_es} eq '+')) ? 1 : 0;                 # 1e-1 => no integer
1615   }
1616
1617 sub is_zero
1618   {
1619   # return true if arg (BFLOAT or num_str) is zero
1620   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1621
1622   ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m})) ? 1 : 0;
1623   }
1624
1625 sub is_one
1626   {
1627   # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1628   my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1629
1630   $sign = '+' if !defined $sign || $sign ne '-';
1631
1632   ($x->{sign} eq $sign &&
1633    $MBI->_is_zero($x->{_e}) &&
1634    $MBI->_is_one($x->{_m}) ) ? 1 : 0;
1635   }
1636
1637 sub is_odd
1638   {
1639   # return true if arg (BFLOAT or num_str) is odd or false if even
1640   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1641
1642   (($x->{sign} =~ /^[+-]$/) &&          # NaN & +-inf aren't
1643    ($MBI->_is_zero($x->{_e})) &&
1644    ($MBI->_is_odd($x->{_m}))) ? 1 : 0;
1645   }
1646
1647 sub is_even
1648   {
1649   # return true if arg (BINT or num_str) is even or false if odd
1650   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1651
1652   (($x->{sign} =~ /^[+-]$/) &&                  # NaN & +-inf aren't
1653    ($x->{_es} eq '+') &&                        # 123.45 isn't
1654    ($MBI->_is_even($x->{_m}))) ? 1 : 0;         # but 1200 is
1655   }
1656
1657 sub bmul
1658   {
1659   # multiply two numbers
1660
1661   # set up parameters
1662   my ($self,$x,$y,@r) = (ref($_[0]),@_);
1663   # objectify is costly, so avoid it
1664   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1665     {
1666     ($self,$x,$y,@r) = objectify(2,@_);
1667     }
1668
1669   return $x if $x->modify('bmul');
1670
1671   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1672
1673   # inf handling
1674   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1675     {
1676     return $x->bnan() if $x->is_zero() || $y->is_zero();
1677     # result will always be +-inf:
1678     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1679     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1680     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1681     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1682     return $x->binf('-');
1683     }
1684
1685   return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1686    ((!$x->isa($self)) || (!$y->isa($self)));
1687
1688   # aEb * cEd = (a*c)E(b+d)
1689   $MBI->_mul($x->{_m},$y->{_m});
1690   ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1691
1692   $r[3] = $y;                           # no push!
1693
1694   # adjust sign:
1695   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1696   $x->bnorm->round(@r);
1697   }
1698
1699 sub bmuladd
1700   {
1701   # multiply two numbers and add the third to the result
1702
1703   # set up parameters
1704   my ($self,$x,$y,$z,@r) = objectify(3,@_);
1705
1706   return $x if $x->modify('bmuladd');
1707
1708   return $x->bnan() if (($x->{sign} eq $nan) ||
1709                         ($y->{sign} eq $nan) ||
1710                         ($z->{sign} eq $nan));
1711
1712   # inf handling
1713   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1714     {
1715     return $x->bnan() if $x->is_zero() || $y->is_zero();
1716     # result will always be +-inf:
1717     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1718     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1719     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1720     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1721     return $x->binf('-');
1722     }
1723
1724   return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1725    ((!$x->isa($self)) || (!$y->isa($self)));
1726
1727   # aEb * cEd = (a*c)E(b+d)
1728   $MBI->_mul($x->{_m},$y->{_m});
1729   ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1730
1731   $r[3] = $y;                           # no push!
1732
1733   # adjust sign:
1734   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1735
1736   # z=inf handling (z=NaN handled above)
1737   $x->{sign} = $z->{sign}, return $x if $z->{sign} =~ /^[+-]inf$/;
1738
1739   # take lower of the two e's and adapt m1 to it to match m2
1740   my $e = $z->{_e};
1741   $e = $MBI->_zero() if !defined $e;            # if no BFLOAT?
1742   $e = $MBI->_copy($e);                         # make copy (didn't do it yet)
1743
1744   my $es;
1745
1746   ($e,$es) = _e_sub($e, $x->{_e}, $z->{_es} || '+', $x->{_es});
1747
1748   my $add = $MBI->_copy($z->{_m});
1749
1750   if ($es eq '-')                               # < 0
1751     {
1752     $MBI->_lsft( $x->{_m}, $e, 10);
1753     ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
1754     }
1755   elsif (!$MBI->_is_zero($e))                   # > 0
1756     {
1757     $MBI->_lsft($add, $e, 10);
1758     }
1759   # else: both e are the same, so just leave them
1760
1761   if ($x->{sign} eq $z->{sign})
1762     {
1763     # add
1764     $x->{_m} = $MBI->_add($x->{_m}, $add);
1765     }
1766   else
1767     {
1768     ($x->{_m}, $x->{sign}) =
1769      _e_add($x->{_m}, $add, $x->{sign}, $z->{sign});
1770     }
1771
1772   # delete trailing zeros, then round
1773   $x->bnorm()->round(@r);
1774   }
1775
1776 sub bdiv
1777   {
1778   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1779   # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1780
1781   # set up parameters
1782   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1783   # objectify is costly, so avoid it
1784   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1785     {
1786     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1787     }
1788
1789   return $x if $x->modify('bdiv');
1790
1791   return $self->_div_inf($x,$y)
1792    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1793
1794   # x== 0 # also: or y == 1 or y == -1
1795   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1796
1797   # upgrade ?
1798   return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1799
1800   # we need to limit the accuracy to protect against overflow
1801   my $fallback = 0;
1802   my (@params,$scale);
1803   ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1804
1805   return $x if $x->is_nan();            # error in _find_round_parameters?
1806
1807   # no rounding at all, so must use fallback
1808   if (scalar @params == 0)
1809     {
1810     # simulate old behaviour
1811     $params[0] = $self->div_scale();    # and round to it as accuracy
1812     $scale = $params[0]+4;              # at least four more for proper round
1813     $params[2] = $r;                    # round mode by caller or undef
1814     $fallback = 1;                      # to clear a/p afterwards
1815     }
1816   else
1817     {
1818     # the 4 below is empirical, and there might be cases where it is not
1819     # enough...
1820     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1821     }
1822
1823   my $rem; $rem = $self->bzero() if wantarray;
1824
1825   $y = $self->new($y) unless $y->isa('Math::BigFloat');
1826
1827   my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
1828   $scale = $lx if $lx > $scale;
1829   $scale = $ly if $ly > $scale;
1830   my $diff = $ly - $lx;
1831   $scale += $diff if $diff > 0;         # if lx << ly, but not if ly << lx!
1832
1833   # already handled inf/NaN/-inf above:
1834
1835   # check that $y is not 1 nor -1 and cache the result:
1836   my $y_not_one = !($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m}));
1837
1838   # flipping the sign of $y will also flip the sign of $x for the special
1839   # case of $x->bsub($x); so we can catch it below:
1840   my $xsign = $x->{sign};
1841   $y->{sign} =~ tr/+-/-+/;
1842
1843   if ($xsign ne $x->{sign})
1844     {
1845     # special case of $x /= $x results in 1
1846     $x->bone();                 # "fixes" also sign of $y, since $x is $y
1847     }
1848   else
1849     {
1850     # correct $y's sign again
1851     $y->{sign} =~ tr/+-/-+/;
1852     # continue with normal div code:
1853
1854     # make copy of $x in case of list context for later remainder calculation
1855     if (wantarray && $y_not_one)
1856       {
1857       $rem = $x->copy();
1858       }
1859
1860     $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1861
1862     # check for / +-1 ( +/- 1E0)
1863     if ($y_not_one)
1864       {
1865       # promote BigInts and it's subclasses (except when already a BigFloat)
1866       $y = $self->new($y) unless $y->isa('Math::BigFloat');
1867
1868       # calculate the result to $scale digits and then round it
1869       # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1870       $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1871       $MBI->_div ($x->{_m},$y->{_m});   # a/c
1872
1873       # correct exponent of $x
1874       ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1875       # correct for 10**scale
1876       ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1877       $x->bnorm();              # remove trailing 0's
1878       }
1879     } # end else $x != $y
1880
1881   # shortcut to not run through _find_round_parameters again
1882   if (defined $params[0])
1883     {
1884     delete $x->{_a};                            # clear before round
1885     $x->bround($params[0],$params[2]);          # then round accordingly
1886     }
1887   else
1888     {
1889     delete $x->{_p};                            # clear before round
1890     $x->bfround($params[1],$params[2]);         # then round accordingly
1891     }
1892   if ($fallback)
1893     {
1894     # clear a/p after round, since user did not request it
1895     delete $x->{_a}; delete $x->{_p};
1896     }
1897
1898   if (wantarray)
1899     {
1900     if ($y_not_one)
1901       {
1902       $x -> bint();
1903       $rem->bmod($y,@params);                   # copy already done
1904       }
1905     if ($fallback)
1906       {
1907       # clear a/p after round, since user did not request it
1908       delete $rem->{_a}; delete $rem->{_p};
1909       }
1910     return ($x,$rem);
1911     }
1912   $x;
1913   }
1914
1915 sub bmod
1916   {
1917   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return remainder
1918
1919   # set up parameters
1920   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1921   # objectify is costly, so avoid it
1922   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1923     {
1924     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1925     }
1926
1927   return $x if $x->modify('bmod');
1928
1929   # handle NaN, inf, -inf
1930   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1931     {
1932     my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1933     $x->{sign} = $re->{sign};
1934     $x->{_e} = $re->{_e};
1935     $x->{_m} = $re->{_m};
1936     return $x->round($a,$p,$r,$y);
1937     }
1938   if ($y->is_zero())
1939     {
1940     return $x->bnan() if $x->is_zero();
1941     return $x;
1942     }
1943
1944   return $x->bzero() if $x->is_zero()
1945  || ($x->is_int() &&
1946   # check that $y == +1 or $y == -1:
1947     ($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m})));
1948
1949   my $cmp = $x->bacmp($y);                      # equal or $x < $y?
1950   return $x->bzero($a,$p) if $cmp == 0;         # $x == $y => result 0
1951
1952   # only $y of the operands negative?
1953   my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1954
1955   $x->{sign} = $y->{sign};                              # calc sign first
1956   return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0;  # $x < $y => result $x
1957
1958   my $ym = $MBI->_copy($y->{_m});
1959
1960   # 2e1 => 20
1961   $MBI->_lsft( $ym, $y->{_e}, 10)
1962    if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
1963
1964   # if $y has digits after dot
1965   my $shifty = 0;                       # correct _e of $x by this
1966   if ($y->{_es} eq '-')                 # has digits after dot
1967     {
1968     # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1969     $shifty = $MBI->_num($y->{_e});     # no more digits after dot
1970     $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
1971     }
1972   # $ym is now mantissa of $y based on exponent 0
1973
1974   my $shiftx = 0;                       # correct _e of $x by this
1975   if ($x->{_es} eq '-')                 # has digits after dot
1976     {
1977     # 123.4 % 20 => 1234 % 200
1978     $shiftx = $MBI->_num($x->{_e});     # no more digits after dot
1979     $MBI->_lsft($ym, $x->{_e}, 10);     # 123 => 1230
1980     }
1981   # 123e1 % 20 => 1230 % 20
1982   if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
1983     {
1984     $MBI->_lsft( $x->{_m}, $x->{_e},10);        # es => '+' here
1985     }
1986
1987   $x->{_e} = $MBI->_new($shiftx);
1988   $x->{_es} = '+';
1989   $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1990   $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
1991
1992   # now mantissas are equalized, exponent of $x is adjusted, so calc result
1993
1994   $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
1995
1996   $x->{sign} = '+' if $MBI->_is_zero($x->{_m});         # fix sign for -0
1997   $x->bnorm();
1998
1999   if ($neg != 0)        # one of them negative => correct in place
2000     {
2001     my $r = $y - $x;
2002     $x->{_m} = $r->{_m};
2003     $x->{_e} = $r->{_e};
2004     $x->{_es} = $r->{_es};
2005     $x->{sign} = '+' if $MBI->_is_zero($x->{_m});       # fix sign for -0
2006     $x->bnorm();
2007     }
2008
2009   $x->round($a,$p,$r,$y);       # round and return
2010   }
2011
2012 sub broot
2013   {
2014   # calculate $y'th root of $x
2015
2016   # set up parameters
2017   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2018   # objectify is costly, so avoid it
2019   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2020     {
2021     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2022     }
2023
2024   return $x if $x->modify('broot');
2025
2026   # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
2027   return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
2028          $y->{sign} !~ /^\+$/;
2029
2030   return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
2031
2032   # we need to limit the accuracy to protect against overflow
2033   my $fallback = 0;
2034   my (@params,$scale);
2035   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2036
2037   return $x if $x->is_nan();            # error in _find_round_parameters?
2038
2039   # no rounding at all, so must use fallback
2040   if (scalar @params == 0)
2041     {
2042     # simulate old behaviour
2043     $params[0] = $self->div_scale();    # and round to it as accuracy
2044     $scale = $params[0]+4;              # at least four more for proper round
2045     $params[2] = $r;                    # round mode by caller or undef
2046     $fallback = 1;                      # to clear a/p afterwards
2047     }
2048   else
2049     {
2050     # the 4 below is empirical, and there might be cases where it is not
2051     # enough...
2052     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2053     }
2054
2055   # when user set globals, they would interfere with our calculation, so
2056   # disable them and later re-enable them
2057   no strict 'refs';
2058   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2059   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2060   # we also need to disable any set A or P on $x (_find_round_parameters took
2061   # them already into account), since these would interfere, too
2062   delete $x->{_a}; delete $x->{_p};
2063   # need to disable $upgrade in BigInt, to avoid deep recursion
2064   local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
2065
2066   # remember sign and make $x positive, since -4 ** (1/2) => -2
2067   my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
2068
2069   my $is_two = 0;
2070   if ($y->isa('Math::BigFloat'))
2071     {
2072     $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
2073     }
2074   else
2075     {
2076     $is_two = ($y == 2);
2077     }
2078
2079   # normal square root if $y == 2:
2080   if ($is_two)
2081     {
2082     $x->bsqrt($scale+4);
2083     }
2084   elsif ($y->is_one('-'))
2085     {
2086     # $x ** -1 => 1/$x
2087     my $u = $self->bone()->bdiv($x,$scale);
2088     # copy private parts over
2089     $x->{_m} = $u->{_m};
2090     $x->{_e} = $u->{_e};
2091     $x->{_es} = $u->{_es};
2092     }
2093   else
2094     {
2095     # calculate the broot() as integer result first, and if it fits, return
2096     # it rightaway (but only if $x and $y are integer):
2097
2098     my $done = 0;                               # not yet
2099     if ($y->is_int() && $x->is_int())
2100       {
2101       my $i = $MBI->_copy( $x->{_m} );
2102       $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2103       my $int = Math::BigInt->bzero();
2104       $int->{value} = $i;
2105       $int->broot($y->as_number());
2106       # if ($exact)
2107       if ($int->copy()->bpow($y) == $x)
2108         {
2109         # found result, return it
2110         $x->{_m} = $int->{value};
2111         $x->{_e} = $MBI->_zero();
2112         $x->{_es} = '+';
2113         $x->bnorm();
2114         $done = 1;
2115         }
2116       }
2117     if ($done == 0)
2118       {
2119       my $u = $self->bone()->bdiv($y,$scale+4);
2120       delete $u->{_a}; delete $u->{_p};         # otherwise it conflicts
2121       $x->bpow($u,$scale+4);                    # el cheapo
2122       }
2123     }
2124   $x->bneg() if $sign == 1;
2125
2126   # shortcut to not run through _find_round_parameters again
2127   if (defined $params[0])
2128     {
2129     $x->bround($params[0],$params[2]);          # then round accordingly
2130     }
2131   else
2132     {
2133     $x->bfround($params[1],$params[2]);         # then round accordingly
2134     }
2135   if ($fallback)
2136     {
2137     # clear a/p after round, since user did not request it
2138     delete $x->{_a}; delete $x->{_p};
2139     }
2140   # restore globals
2141   $$abr = $ab; $$pbr = $pb;
2142   $x;
2143   }
2144
2145 sub bsqrt
2146   {
2147   # calculate square root
2148   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2149
2150   return $x if $x->modify('bsqrt');
2151
2152   return $x->bnan() if $x->{sign} !~ /^[+]/;    # NaN, -inf or < 0
2153   return $x if $x->{sign} eq '+inf';            # sqrt(inf) == inf
2154   return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
2155
2156   # we need to limit the accuracy to protect against overflow
2157   my $fallback = 0;
2158   my (@params,$scale);
2159   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2160
2161   return $x if $x->is_nan();            # error in _find_round_parameters?
2162
2163   # no rounding at all, so must use fallback
2164   if (scalar @params == 0)
2165     {
2166     # simulate old behaviour
2167     $params[0] = $self->div_scale();    # and round to it as accuracy
2168     $scale = $params[0]+4;              # at least four more for proper round
2169     $params[2] = $r;                    # round mode by caller or undef
2170     $fallback = 1;                      # to clear a/p afterwards
2171     }
2172   else
2173     {
2174     # the 4 below is empirical, and there might be cases where it is not
2175     # enough...
2176     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2177     }
2178
2179   # when user set globals, they would interfere with our calculation, so
2180   # disable them and later re-enable them
2181   no strict 'refs';
2182   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2183   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2184   # we also need to disable any set A or P on $x (_find_round_parameters took
2185   # them already into account), since these would interfere, too
2186   delete $x->{_a}; delete $x->{_p};
2187   # need to disable $upgrade in BigInt, to avoid deep recursion
2188   local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
2189
2190   my $i = $MBI->_copy( $x->{_m} );
2191   $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2192   my $xas = Math::BigInt->bzero();
2193   $xas->{value} = $i;
2194
2195   my $gs = $xas->copy()->bsqrt();       # some guess
2196
2197   if (($x->{_es} ne '-')                # guess can't be accurate if there are
2198                                         # digits after the dot
2199    && ($xas->bacmp($gs * $gs) == 0))    # guess hit the nail on the head?
2200     {
2201     # exact result, copy result over to keep $x
2202     $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
2203     $x->bnorm();
2204     # shortcut to not run through _find_round_parameters again
2205     if (defined $params[0])
2206       {
2207       $x->bround($params[0],$params[2]);        # then round accordingly
2208       }
2209     else
2210       {
2211       $x->bfround($params[1],$params[2]);       # then round accordingly
2212       }
2213     if ($fallback)
2214       {
2215       # clear a/p after round, since user did not request it
2216       delete $x->{_a}; delete $x->{_p};
2217       }
2218     # re-enable A and P, upgrade is taken care of by "local"
2219     ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
2220     return $x;
2221     }
2222
2223   # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
2224   # of the result by multiplying the input by 100 and then divide the integer
2225   # result of sqrt(input) by 10. Rounding afterwards returns the real result.
2226
2227   # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
2228   my $y1 = $MBI->_copy($x->{_m});
2229
2230   my $length = $MBI->_len($y1);
2231
2232   # Now calculate how many digits the result of sqrt(y1) would have
2233   my $digits = int($length / 2);
2234
2235   # But we need at least $scale digits, so calculate how many are missing
2236   my $shift = $scale - $digits;
2237
2238   # This happens if the input had enough digits
2239   # (we take care of integer guesses above)
2240   $shift = 0 if $shift < 0;
2241
2242   # Multiply in steps of 100, by shifting left two times the "missing" digits
2243   my $s2 = $shift * 2;
2244
2245   # We now make sure that $y1 has the same odd or even number of digits than
2246   # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
2247   # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
2248   # steps of 10. The length of $x does not count, since an even or odd number
2249   # of digits before the dot is not changed by adding an even number of digits
2250   # after the dot (the result is still odd or even digits long).
2251   $s2++ if $MBI->_is_odd($x->{_e});
2252
2253   $MBI->_lsft( $y1, $MBI->_new($s2), 10);
2254
2255   # now take the square root and truncate to integer
2256   $y1 = $MBI->_sqrt($y1);
2257
2258   # By "shifting" $y1 right (by creating a negative _e) we calculate the final
2259   # result, which is than later rounded to the desired scale.
2260
2261   # calculate how many zeros $x had after the '.' (or before it, depending
2262   # on sign of $dat, the result should have half as many:
2263   my $dat = $MBI->_num($x->{_e});
2264   $dat = -$dat if $x->{_es} eq '-';
2265   $dat += $length;
2266
2267   if ($dat > 0)
2268     {
2269     # no zeros after the dot (e.g. 1.23, 0.49 etc)
2270     # preserve half as many digits before the dot than the input had
2271     # (but round this "up")
2272     $dat = int(($dat+1)/2);
2273     }
2274   else
2275     {
2276     $dat = int(($dat)/2);
2277     }
2278   $dat -= $MBI->_len($y1);
2279   if ($dat < 0)
2280     {
2281     $dat = abs($dat);
2282     $x->{_e} = $MBI->_new( $dat );
2283     $x->{_es} = '-';
2284     }
2285   else
2286     {
2287     $x->{_e} = $MBI->_new( $dat );
2288     $x->{_es} = '+';
2289     }
2290   $x->{_m} = $y1;
2291   $x->bnorm();
2292
2293   # shortcut to not run through _find_round_parameters again
2294   if (defined $params[0])
2295     {
2296     $x->bround($params[0],$params[2]);          # then round accordingly
2297     }
2298   else
2299     {
2300     $x->bfround($params[1],$params[2]);         # then round accordingly
2301     }
2302   if ($fallback)
2303     {
2304     # clear a/p after round, since user did not request it
2305     delete $x->{_a}; delete $x->{_p};
2306     }
2307   # restore globals
2308   $$abr = $ab; $$pbr = $pb;
2309   $x;
2310   }
2311
2312 sub bfac
2313   {
2314   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2315   # compute factorial number, modifies first argument
2316
2317   # set up parameters
2318   my ($self,$x,@r) = (ref($_[0]),@_);
2319   # objectify is costly, so avoid it
2320   ($self,$x,@r) = objectify(1,@_) if !ref($x);
2321
2322   # inf => inf
2323   return $x if $x->modify('bfac') || $x->{sign} eq '+inf';
2324
2325   return $x->bnan()
2326     if (($x->{sign} ne '+') ||          # inf, NaN, <0 etc => NaN
2327      ($x->{_es} ne '+'));               # digits after dot?
2328
2329   # use BigInt's bfac() for faster calc
2330   if (! $MBI->_is_zero($x->{_e}))
2331     {
2332     $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
2333     $x->{_e} = $MBI->_zero();           # normalize
2334     $x->{_es} = '+';
2335     }
2336   $MBI->_fac($x->{_m});                 # calculate factorial
2337   $x->bnorm()->round(@r);               # norm again and round result
2338   }
2339
2340 sub _pow
2341   {
2342   # Calculate a power where $y is a non-integer, like 2 ** 0.3
2343   my ($x,$y,@r) = @_;
2344   my $self = ref($x);
2345
2346   # if $y == 0.5, it is sqrt($x)
2347   $HALF = $self->new($HALF) unless ref($HALF);
2348   return $x->bsqrt(@r,$y) if $y->bcmp($HALF) == 0;
2349
2350   # Using:
2351   # a ** x == e ** (x * ln a)
2352
2353   # u = y * ln x
2354   #                _                         _
2355   # Taylor:       |   u    u^2    u^3         |
2356   # x ** y  = 1 + |  --- + --- + ----- + ...  |
2357   #               |_  1    1*2   1*2*3       _|
2358
2359   # we need to limit the accuracy to protect against overflow
2360   my $fallback = 0;
2361   my ($scale,@params);
2362   ($x,@params) = $x->_find_round_parameters(@r);
2363
2364   return $x if $x->is_nan();            # error in _find_round_parameters?
2365
2366   # no rounding at all, so must use fallback
2367   if (scalar @params == 0)
2368     {
2369     # simulate old behaviour
2370     $params[0] = $self->div_scale();    # and round to it as accuracy
2371     $params[1] = undef;                 # disable P
2372     $scale = $params[0]+4;              # at least four more for proper round
2373     $params[2] = $r[2];                 # round mode by caller or undef
2374     $fallback = 1;                      # to clear a/p afterwards
2375     }
2376   else
2377     {
2378     # the 4 below is empirical, and there might be cases where it is not
2379     # enough...
2380     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2381     }
2382
2383   # when user set globals, they would interfere with our calculation, so
2384   # disable them and later re-enable them
2385   no strict 'refs';
2386   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2387   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2388   # we also need to disable any set A or P on $x (_find_round_parameters took
2389   # them already into account), since these would interfere, too
2390   delete $x->{_a}; delete $x->{_p};
2391   # need to disable $upgrade in BigInt, to avoid deep recursion
2392   local $Math::BigInt::upgrade = undef;
2393
2394   my ($limit,$v,$u,$below,$factor,$next,$over);
2395
2396   $u = $x->copy()->blog(undef,$scale)->bmul($y);
2397   $v = $self->bone();                           # 1
2398   $factor = $self->new(2);                      # 2
2399   $x->bone();                                   # first term: 1
2400
2401   $below = $v->copy();
2402   $over = $u->copy();
2403
2404   $limit = $self->new("1E-". ($scale-1));
2405   #my $steps = 0;
2406   while (3 < 5)
2407     {
2408     # we calculate the next term, and add it to the last
2409     # when the next term is below our limit, it won't affect the outcome
2410     # anymore, so we stop:
2411     $next = $over->copy()->bdiv($below,$scale);
2412     last if $next->bacmp($limit) <= 0;
2413     $x->badd($next);
2414     # calculate things for the next term
2415     $over *= $u; $below *= $factor; $factor->binc();
2416
2417     last if $x->{sign} !~ /^[-+]$/;
2418
2419     #$steps++;
2420     }
2421
2422   # shortcut to not run through _find_round_parameters again
2423   if (defined $params[0])
2424     {
2425     $x->bround($params[0],$params[2]);          # then round accordingly
2426     }
2427   else
2428     {
2429     $x->bfround($params[1],$params[2]);         # then round accordingly
2430     }
2431   if ($fallback)
2432     {
2433     # clear a/p after round, since user did not request it
2434     delete $x->{_a}; delete $x->{_p};
2435     }
2436   # restore globals
2437   $$abr = $ab; $$pbr = $pb;
2438   $x;
2439   }
2440
2441 sub bpow
2442   {
2443   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2444   # compute power of two numbers, second arg is used as integer
2445   # modifies first argument
2446
2447   # set up parameters
2448   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2449   # objectify is costly, so avoid it
2450   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2451     {
2452     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2453     }
2454
2455   return $x if $x->modify('bpow');
2456
2457   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
2458   return $x if $x->{sign} =~ /^[+-]inf$/;
2459
2460   # cache the result of is_zero
2461   my $y_is_zero = $y->is_zero();
2462   return $x->bone() if $y_is_zero;
2463   return $x         if $x->is_one() || $y->is_one();
2464
2465   my $x_is_zero = $x->is_zero();
2466   return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int();         # non-integer power
2467
2468   my $y1 = $y->as_number()->{value};                    # make MBI part
2469
2470   # if ($x == -1)
2471   if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
2472     {
2473     # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
2474     return $MBI->_is_odd($y1) ? $x : $x->babs(1);
2475     }
2476   if ($x_is_zero)
2477     {
2478     return $x if $y->{sign} eq '+';     # 0**y => 0 (if not y <= 0)
2479     # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
2480     return $x->binf();
2481     }
2482
2483   my $new_sign = '+';
2484   $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
2485
2486   # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
2487   $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
2488   $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
2489
2490   $x->{sign} = $new_sign;
2491   $x->bnorm();
2492   if ($y->{sign} eq '-')
2493     {
2494     # modify $x in place!
2495     my $z = $x->copy(); $x->bone();
2496     return scalar $x->bdiv($z,$a,$p,$r);        # round in one go (might ignore y's A!)
2497     }
2498   $x->round($a,$p,$r,$y);
2499   }
2500
2501 sub bmodpow
2502   {
2503   # takes a very large number to a very large exponent in a given very
2504   # large modulus, quickly, thanks to binary exponentiation. Supports
2505   # negative exponents.
2506   my ($self,$num,$exp,$mod,@r) = objectify(3,@_);
2507
2508   return $num if $num->modify('bmodpow');
2509
2510   # check modulus for valid values
2511   return $num->bnan() if ($mod->{sign} ne '+'           # NaN, - , -inf, +inf
2512                        || $mod->is_zero());
2513
2514   # check exponent for valid values
2515   if ($exp->{sign} =~ /\w/)
2516     {
2517     # i.e., if it's NaN, +inf, or -inf...
2518     return $num->bnan();
2519     }
2520
2521   $num->bmodinv ($mod) if ($exp->{sign} eq '-');
2522
2523   # check num for valid values (also NaN if there was no inverse but $exp < 0)
2524   return $num->bnan() if $num->{sign} !~ /^[+-]$/;
2525
2526   # $mod is positive, sign on $exp is ignored, result also positive
2527
2528   # XXX TODO: speed it up when all three numbers are integers
2529   $num->bpow($exp)->bmod($mod);
2530   }
2531
2532 ###############################################################################
2533 # trigonometric functions
2534
2535 # helper function for bpi() and batan2(), calculates arcus tanges (1/x)
2536
2537 sub _atan_inv
2538   {
2539   # return a/b so that a/b approximates atan(1/x) to at least limit digits
2540   my ($self, $x, $limit) = @_;
2541
2542   # Taylor:       x^3   x^5   x^7   x^9
2543   #    atan = x - --- + --- - --- + --- - ...
2544   #                3     5     7     9
2545
2546   #               1      1         1        1
2547   #    atan 1/x = - - ------- + ------- - ------- + ...
2548   #               x   x^3 * 3   x^5 * 5   x^7 * 7
2549
2550   #               1      1         1            1
2551   #    atan 1/x = - - --------- + ---------- - ----------- + ...
2552   #               5    3 * 125     5 * 3125     7 * 78125
2553
2554   # Subtraction/addition of a rational:
2555
2556   #  5    7    5*3 +- 7*4
2557   #  - +- -  = ----------
2558   #  4    3       4*3
2559
2560   # Term:  N        N+1
2561   #
2562   #        a             1                  a * d * c +- b
2563   #        ----- +- ------------------  =  ----------------
2564   #        b           d * c                b * d * c
2565
2566   #  since b1 = b0 * (d-2) * c
2567
2568   #        a             1                  a * d +- b / c
2569   #        ----- +- ------------------  =  ----------------
2570   #        b           d * c                b * d
2571
2572   # and  d = d + 2
2573   # and  c = c * x * x
2574
2575   #        u = d * c
2576   #        stop if length($u) > limit
2577   #        a = a * u +- b
2578   #        b = b * u
2579   #        d = d + 2
2580   #        c = c * x * x
2581   #        sign = 1 - sign
2582
2583   my $a = $MBI->_one();
2584   my $b = $MBI->_copy($x);
2585
2586   my $x2  = $MBI->_mul( $MBI->_copy($x), $b);           # x2 = x * x
2587   my $d   = $MBI->_new( 3 );                            # d = 3
2588   my $c   = $MBI->_mul( $MBI->_copy($x), $x2);          # c = x ^ 3
2589   my $two = $MBI->_new( 2 );
2590
2591   # run the first step unconditionally
2592   my $u = $MBI->_mul( $MBI->_copy($d), $c);
2593   $a = $MBI->_mul($a, $u);
2594   $a = $MBI->_sub($a, $b);
2595   $b = $MBI->_mul($b, $u);
2596   $d = $MBI->_add($d, $two);
2597   $c = $MBI->_mul($c, $x2);
2598
2599   # a is now a * (d-3) * c
2600   # b is now b * (d-2) * c
2601
2602   # run the second step unconditionally
2603   $u = $MBI->_mul( $MBI->_copy($d), $c);
2604   $a = $MBI->_mul($a, $u);
2605   $a = $MBI->_add($a, $b);
2606   $b = $MBI->_mul($b, $u);
2607   $d = $MBI->_add($d, $two);
2608   $c = $MBI->_mul($c, $x2);
2609
2610   # a is now a * (d-3) * (d-5) * c * c
2611   # b is now b * (d-2) * (d-4) * c * c
2612
2613   # so we can remove c * c from both a and b to shorten the numbers involved:
2614   $a = $MBI->_div($a, $x2);
2615   $b = $MBI->_div($b, $x2);
2616   $a = $MBI->_div($a, $x2);
2617   $b = $MBI->_div($b, $x2);
2618
2619 #  my $step = 0;
2620   my $sign = 0;                                         # 0 => -, 1 => +
2621   while (3 < 5)
2622     {
2623 #    $step++;
2624 #    if (($i++ % 100) == 0)
2625 #      {
2626 #    print "a=",$MBI->_str($a),"\n";
2627 #    print "b=",$MBI->_str($b),"\n";
2628 #      }
2629 #    print "d=",$MBI->_str($d),"\n";
2630 #    print "x2=",$MBI->_str($x2),"\n";
2631 #    print "c=",$MBI->_str($c),"\n";
2632
2633     my $u = $MBI->_mul( $MBI->_copy($d), $c);
2634     # use _alen() for libs like GMP where _len() would be O(N^2)
2635     last if $MBI->_alen($u) > $limit;
2636     my ($bc,$r) = $MBI->_div( $MBI->_copy($b), $c);
2637     if ($MBI->_is_zero($r))
2638       {
2639       # b / c is an integer, so we can remove c from all terms
2640       # this happens almost every time:
2641       $a = $MBI->_mul($a, $d);
2642       $a = $MBI->_sub($a, $bc) if $sign == 0;
2643       $a = $MBI->_add($a, $bc) if $sign == 1;
2644       $b = $MBI->_mul($b, $d);
2645       }
2646     else
2647       {
2648       # b / c is not an integer, so we keep c in the terms
2649       # this happens very rarely, for instance for x = 5, this happens only
2650       # at the following steps:
2651       # 1, 5, 14, 32, 72, 157, 340, ...
2652       $a = $MBI->_mul($a, $u);
2653       $a = $MBI->_sub($a, $b) if $sign == 0;
2654       $a = $MBI->_add($a, $b) if $sign == 1;
2655       $b = $MBI->_mul($b, $u);
2656       }
2657     $d = $MBI->_add($d, $two);
2658     $c = $MBI->_mul($c, $x2);
2659     $sign = 1 - $sign;
2660
2661     }
2662
2663 #  print "Took $step steps for ", $MBI->_str($x),"\n";
2664 #  print "a=",$MBI->_str($a),"\n"; print "b=",$MBI->_str($b),"\n";
2665   # return a/b so that a/b approximates atan(1/x)
2666   ($a,$b);
2667   }
2668
2669 sub bpi
2670   {
2671   my ($self,$n) = @_;
2672   if (@_ == 0)
2673     {
2674     $self = $class;
2675     }
2676   if (@_ == 1)
2677     {
2678     # called like Math::BigFloat::bpi(10);
2679     $n = $self; $self = $class;
2680     # called like Math::BigFloat->bpi();
2681     $n = undef if $n eq 'Math::BigFloat';
2682     }
2683   $self = ref($self) if ref($self);
2684   my $fallback = defined $n ? 0 : 1;
2685   $n = 40 if !defined $n || $n < 1;
2686
2687   # after 黃見利 (Hwang Chien-Lih) (1997)
2688   # pi/4 = 183 * atan(1/239) + 32 * atan(1/1023) – 68 * atan(1/5832)
2689   #      + 12 * atan(1/110443) - 12 * atan(1/4841182) - 100 * atan(1/6826318)
2690
2691   # a few more to prevent rounding errors
2692   $n += 4;
2693
2694   my ($a,$b) = $self->_atan_inv( $MBI->_new(239),$n);
2695   my ($c,$d) = $self->_atan_inv( $MBI->_new(1023),$n);
2696   my ($e,$f) = $self->_atan_inv( $MBI->_new(5832),$n);
2697   my ($g,$h) = $self->_atan_inv( $MBI->_new(110443),$n);
2698   my ($i,$j) = $self->_atan_inv( $MBI->_new(4841182),$n);
2699   my ($k,$l) = $self->_atan_inv( $MBI->_new(6826318),$n);
2700
2701   $MBI->_mul($a, $MBI->_new(732));
2702   $MBI->_mul($c, $MBI->_new(128));
2703   $MBI->_mul($e, $MBI->_new(272));
2704   $MBI->_mul($g, $MBI->_new(48));
2705   $MBI->_mul($i, $MBI->_new(48));
2706   $MBI->_mul($k, $MBI->_new(400));
2707
2708   my $x = $self->bone(); $x->{_m} = $a; my $x_d = $self->bone(); $x_d->{_m} = $b;
2709   my $y = $self->bone(); $y->{_m} = $c; my $y_d = $self->bone(); $y_d->{_m} = $d;
2710   my $z = $self->bone(); $z->{_m} = $e; my $z_d = $self->bone(); $z_d->{_m} = $f;
2711   my $u = $self->bone(); $u->{_m} = $g; my $u_d = $self->bone(); $u_d->{_m} = $h;
2712   my $v = $self->bone(); $v->{_m} = $i; my $v_d = $self->bone(); $v_d->{_m} = $j;
2713   my $w = $self->bone(); $w->{_m} = $k; my $w_d = $self->bone(); $w_d->{_m} = $l;
2714   $x->bdiv($x_d, $n);
2715   $y->bdiv($y_d, $n);
2716   $z->bdiv($z_d, $n);
2717   $u->bdiv($u_d, $n);
2718   $v->bdiv($v_d, $n);
2719   $w->bdiv($w_d, $n);
2720
2721   delete $x->{_a}; delete $y->{_a}; delete $z->{_a};
2722   delete $u->{_a}; delete $v->{_a}; delete $w->{_a};
2723   $x->badd($y)->bsub($z)->badd($u)->bsub($v)->bsub($w);
2724
2725   $x->bround($n-4);
2726   delete $x->{_a} if $fallback == 1;
2727   $x;
2728   }
2729
2730 sub bcos
2731   {
2732   # Calculate a cosinus of x.
2733   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2734
2735   # Taylor:      x^2   x^4   x^6   x^8
2736   #    cos = 1 - --- + --- - --- + --- ...
2737   #               2!    4!    6!    8!
2738
2739   # we need to limit the accuracy to protect against overflow
2740   my $fallback = 0;
2741   my ($scale,@params);
2742   ($x,@params) = $x->_find_round_parameters(@r);
2743
2744   #         constant object       or error in _find_round_parameters?
2745   return $x if $x->modify('bcos') || $x->is_nan();
2746
2747   return $x->bone(@r) if $x->is_zero();
2748
2749   # no rounding at all, so must use fallback
2750   if (scalar @params == 0)
2751     {
2752     # simulate old behaviour
2753     $params[0] = $self->div_scale();    # and round to it as accuracy
2754     $params[1] = undef;                 # disable P
2755     $scale = $params[0]+4;              # at least four more for proper round
2756     $params[2] = $r[2];                 # round mode by caller or undef
2757     $fallback = 1;                      # to clear a/p afterwards
2758     }
2759   else
2760     {
2761     # the 4 below is empirical, and there might be cases where it is not
2762     # enough...
2763     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2764     }
2765
2766   # when user set globals, they would interfere with our calculation, so
2767   # disable them and later re-enable them
2768   no strict 'refs';
2769   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2770   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2771   # we also need to disable any set A or P on $x (_find_round_parameters took
2772   # them already into account), since these would interfere, too
2773   delete $x->{_a}; delete $x->{_p};
2774   # need to disable $upgrade in BigInt, to avoid deep recursion
2775   local $Math::BigInt::upgrade = undef;
2776
2777   my $last = 0;
2778   my $over = $x * $x;                   # X ^ 2
2779   my $x2 = $over->copy();               # X ^ 2; difference between terms
2780   my $sign = 1;                         # start with -=
2781   my $below = $self->new(2); my $factorial = $self->new(3);
2782   $x->bone(); delete $x->{_a}; delete $x->{_p};
2783
2784   my $limit = $self->new("1E-". ($scale-1));
2785   #my $steps = 0;
2786   while (3 < 5)
2787     {
2788     # we calculate the next term, and add it to the last
2789     # when the next term is below our limit, it won't affect the outcome
2790     # anymore, so we stop:
2791     my $next = $over->copy()->bdiv($below,$scale);
2792     last if $next->bacmp($limit) <= 0;
2793
2794     if ($sign == 0)
2795       {
2796       $x->badd($next);
2797       }
2798     else
2799       {
2800       $x->bsub($next);
2801       }
2802     $sign = 1-$sign;                                    # alternate
2803     # calculate things for the next term
2804     $over->bmul($x2);                                   # $x*$x
2805     $below->bmul($factorial); $factorial->binc();       # n*(n+1)
2806     $below->bmul($factorial); $factorial->binc();       # n*(n+1)
2807     }
2808
2809   # shortcut to not run through _find_round_parameters again
2810   if (defined $params[0])
2811     {
2812     $x->bround($params[0],$params[2]);          # then round accordingly
2813     }
2814   else
2815     {
2816     $x->bfround($params[1],$params[2]);         # then round accordingly
2817     }
2818   if ($fallback)
2819     {
2820     # clear a/p after round, since user did not request it
2821     delete $x->{_a}; delete $x->{_p};
2822     }
2823   # restore globals
2824   $$abr = $ab; $$pbr = $pb;
2825   $x;
2826   }
2827
2828 sub bsin
2829   {
2830   # Calculate a sinus of x.
2831   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2832
2833   # taylor:      x^3   x^5   x^7   x^9
2834   #    sin = x - --- + --- - --- + --- ...
2835   #               3!    5!    7!    9!
2836
2837   # we need to limit the accuracy to protect against overflow
2838   my $fallback = 0;
2839   my ($scale,@params);
2840   ($x,@params) = $x->_find_round_parameters(@r);
2841
2842   #         constant object       or error in _find_round_parameters?
2843   return $x if $x->modify('bsin') || $x->is_nan();
2844
2845   return $x->bzero(@r) if $x->is_zero();
2846
2847   # no rounding at all, so must use fallback
2848   if (scalar @params == 0)
2849     {
2850     # simulate old behaviour
2851     $params[0] = $self->div_scale();    # and round to it as accuracy
2852     $params[1] = undef;                 # disable P
2853     $scale = $params[0]+4;              # at least four more for proper round
2854     $params[2] = $r[2];                 # round mode by caller or undef
2855     $fallback = 1;                      # to clear a/p afterwards
2856     }
2857   else
2858     {
2859     # the 4 below is empirical, and there might be cases where it is not
2860     # enough...
2861     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2862     }
2863
2864   # when user set globals, they would interfere with our calculation, so
2865   # disable them and later re-enable them
2866   no strict 'refs';
2867   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2868   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2869   # we also need to disable any set A or P on $x (_find_round_parameters took
2870   # them already into account), since these would interfere, too
2871   delete $x->{_a}; delete $x->{_p};
2872   # need to disable $upgrade in BigInt, to avoid deep recursion
2873   local $Math::BigInt::upgrade = undef;
2874
2875   my $last = 0;
2876   my $over = $x * $x;                   # X ^ 2
2877   my $x2 = $over->copy();               # X ^ 2; difference between terms
2878   $over->bmul($x);                      # X ^ 3 as starting value
2879   my $sign = 1;                         # start with -=
2880   my $below = $self->new(6); my $factorial = $self->new(4);
2881   delete $x->{_a}; delete $x->{_p};
2882
2883   my $limit = $self->new("1E-". ($scale-1));
2884   #my $steps = 0;
2885   while (3 < 5)
2886     {
2887     # we calculate the next term, and add it to the last
2888     # when the next term is below our limit, it won't affect the outcome
2889     # anymore, so we stop:
2890     my $next = $over->copy()->bdiv($below,$scale);
2891     last if $next->bacmp($limit) <= 0;
2892
2893     if ($sign == 0)
2894       {
2895       $x->badd($next);
2896       }
2897     else
2898       {
2899       $x->bsub($next);
2900       }
2901     $sign = 1-$sign;                                    # alternate
2902     # calculate things for the next term
2903     $over->bmul($x2);                                   # $x*$x
2904     $below->bmul($factorial); $factorial->binc();       # n*(n+1)
2905     $below->bmul($factorial); $factorial->binc();       # n*(n+1)
2906     }
2907
2908   # shortcut to not run through _find_round_parameters again
2909   if (defined $params[0])
2910     {
2911     $x->bround($params[0],$params[2]);          # then round accordingly
2912     }
2913   else
2914     {
2915     $x->bfround($params[1],$params[2]);         # then round accordingly
2916     }
2917   if ($fallback)
2918     {
2919     # clear a/p after round, since user did not request it
2920     delete $x->{_a}; delete $x->{_p};
2921     }
2922   # restore globals
2923   $$abr = $ab; $$pbr = $pb;
2924   $x;
2925   }
2926
2927 sub batan2
2928   {
2929   # calculate arcus tangens of ($y/$x)
2930
2931   # set up parameters
2932   my ($self,$y,$x,@r) = (ref($_[0]),@_);
2933   # objectify is costly, so avoid it
2934   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2935     {
2936     ($self,$y,$x,@r) = objectify(2,@_);
2937     }
2938
2939   return $y if $y->modify('batan2');
2940
2941   return $y->bnan() if ($y->{sign} eq $nan) || ($x->{sign} eq $nan);
2942
2943   # Y X
2944   # 0 0 result is 0
2945   # 0 +x result is 0
2946   # ? inf result is 0
2947   return $y->bzero(@r) if ($x->is_inf('+') && !$y->is_inf()) || ($y->is_zero() && $x->{sign} eq '+');
2948
2949   # Y    X
2950   # != 0 -inf result is +- pi
2951   if ($x->is_inf() || $y->is_inf())
2952     {
2953     # calculate PI
2954     my $pi = $self->bpi(@r);
2955     if ($y->is_inf())
2956       {
2957       # upgrade to BigRat etc.
2958       return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
2959       if ($x->{sign} eq '-inf')
2960         {
2961         # calculate 3 pi/4
2962         $MBI->_mul($pi->{_m}, $MBI->_new(3));
2963         $MBI->_div($pi->{_m}, $MBI->_new(4));
2964         }
2965       elsif ($x->{sign} eq '+inf')
2966         {
2967         # calculate pi/4
2968         $MBI->_div($pi->{_m}, $MBI->_new(4));
2969         }
2970       else
2971         {
2972         # calculate pi/2
2973         $MBI->_div($pi->{_m}, $MBI->_new(2));
2974         }
2975       $y->{sign} = substr($y->{sign},0,1); # keep +/-
2976       }
2977     # modify $y in place
2978     $y->{_m} = $pi->{_m};
2979     $y->{_e} = $pi->{_e};
2980     $y->{_es} = $pi->{_es};
2981     # keep the sign of $y
2982     return $y;
2983     }
2984
2985   return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
2986
2987   # Y X
2988   # 0 -x result is PI
2989   if ($y->is_zero())
2990     {
2991     # calculate PI
2992     my $pi = $self->bpi(@r);
2993     # modify $y in place
2994     $y->{_m} = $pi->{_m};
2995     $y->{_e} = $pi->{_e};
2996     $y->{_es} = $pi->{_es};
2997     $y->{sign} = '+';
2998     return $y;
2999     }
3000
3001   # Y X
3002   # +y 0 result is PI/2
3003   # -y 0 result is -PI/2
3004   if ($x->is_zero())
3005     {
3006     # calculate PI/2
3007     my $pi = $self->bpi(@r);
3008     # modify $y in place
3009     $y->{_m} = $pi->{_m};
3010     $y->{_e} = $pi->{_e};
3011     $y->{_es} = $pi->{_es};
3012     # -y => -PI/2, +y => PI/2
3013     $MBI->_div($y->{_m}, $MBI->_new(2));
3014     return $y;
3015     }
3016
3017   # we need to limit the accuracy to protect against overflow
3018   my $fallback = 0;
3019   my ($scale,@params);
3020   ($y,@params) = $y->_find_round_parameters(@r);
3021
3022   # error in _find_round_parameters?
3023   return $y if $y->is_nan();
3024
3025   # no rounding at all, so must use fallback
3026   if (scalar @params == 0)
3027     {
3028     # simulate old behaviour
3029     $params[0] = $self->div_scale();    # and round to it as accuracy
3030     $params[1] = undef;                 # disable P
3031     $scale = $params[0]+4;              # at least four more for proper round
3032     $params[2] = $r[2];                 # round mode by caller or undef
3033     $fallback = 1;                      # to clear a/p afterwards
3034     }
3035   else
3036     {
3037     # the 4 below is empirical, and there might be cases where it is not
3038     # enough...
3039     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3040     }
3041
3042   # inlined is_one() && is_one('-')
3043   if ($MBI->_is_one($y->{_m}) && $MBI->_is_zero($y->{_e}))
3044     {
3045     # shortcut: 1 1 result is PI/4
3046     # inlined is_one() && is_one('-')
3047     if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
3048       {
3049       # 1,1 => PI/4
3050       my $pi_4 = $self->bpi( $scale - 3);
3051       # modify $y in place
3052       $y->{_m} = $pi_4->{_m};
3053       $y->{_e} = $pi_4->{_e};
3054       $y->{_es} = $pi_4->{_es};
3055       # 1 1 => +
3056       # -1 1 => -
3057       # 1 -1 => -
3058       # -1 -1 => +
3059       $y->{sign} = $x->{sign} eq $y->{sign} ? '+' : '-';
3060       $MBI->_div($y->{_m}, $MBI->_new(4));
3061       return $y;
3062       }
3063     # shortcut: 1 int(X) result is _atan_inv(X)
3064
3065     # is integer
3066     if ($x->{_es} eq '+')
3067       {
3068       my $x1 = $MBI->_copy($x->{_m});
3069       $MBI->_lsft($x1, $x->{_e},10) unless $MBI->_is_zero($x->{_e});
3070
3071       my ($a,$b) = $self->_atan_inv($x1, $scale);
3072       my $y_sign = $y->{sign};
3073       # calculate A/B
3074       $y->bone(); $y->{_m} = $a; my $y_d = $self->bone(); $y_d->{_m} = $b;
3075       $y->bdiv($y_d, @r);
3076       $y->{sign} = $y_sign;
3077       return $y;
3078       }
3079     }
3080
3081   # handle all other cases
3082   #  X  Y
3083   # +x +y 0 to PI/2
3084   # -x +y PI/2 to PI
3085   # +x -y 0 to -PI/2
3086   # -x -y -PI/2 to -PI
3087
3088   my $y_sign = $y->{sign};
3089
3090   # divide $x by $y
3091   $y->bdiv($x, $scale) unless $x->is_one();
3092   $y->batan(@r);
3093
3094   # restore sign
3095   $y->{sign} = $y_sign;
3096
3097   $y;
3098   }
3099
3100 sub batan
3101   {
3102   # Calculate a arcus tangens of x.
3103   my ($x,@r) = @_;
3104   my $self = ref($x);
3105
3106   # taylor:       x^3   x^5   x^7   x^9
3107   #    atan = x - --- + --- - --- + --- ...
3108   #                3     5     7     9
3109
3110   # we need to limit the accuracy to protect against overflow
3111   my $fallback = 0;
3112   my ($scale,@params);
3113   ($x,@params) = $x->_find_round_parameters(@r);
3114
3115   #         constant object       or error in _find_round_parameters?
3116   return $x if $x->modify('batan') || $x->is_nan();
3117
3118   if ($x->{sign} =~ /^[+-]inf\z/)
3119     {
3120     # +inf result is PI/2
3121     # -inf result is -PI/2
3122     # calculate PI/2
3123     my $pi = $self->bpi(@r);
3124     # modify $x in place
3125     $x->{_m} = $pi->{_m};
3126     $x->{_e} = $pi->{_e};
3127     $x->{_es} = $pi->{_es};
3128     # -y => -PI/2, +y => PI/2
3129     $x->{sign} = substr($x->{sign},0,1);                # +inf => +
3130     $MBI->_div($x->{_m}, $MBI->_new(2));
3131     return $x;
3132     }
3133
3134   return $x->bzero(@r) if $x->is_zero();
3135
3136   # no rounding at all, so must use fallback
3137   if (scalar @params == 0)
3138     {
3139     # simulate old behaviour
3140     $params[0] = $self->div_scale();    # and round to it as accuracy
3141     $params[1] = undef;                 # disable P
3142     $scale = $params[0]+4;              # at least four more for proper round
3143     $params[2] = $r[2];                 # round mode by caller or undef
3144     $fallback = 1;                      # to clear a/p afterwards
3145     }
3146   else
3147     {
3148     # the 4 below is empirical, and there might be cases where it is not
3149     # enough...
3150     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3151     }
3152
3153   # 1 or -1 => PI/4
3154   # inlined is_one() && is_one('-')
3155   if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
3156     {
3157     my $pi = $self->bpi($scale - 3);
3158     # modify $x in place
3159     $x->{_m} = $pi->{_m};
3160     $x->{_e} = $pi->{_e};
3161     $x->{_es} = $pi->{_es};
3162     # leave the sign of $x alone (+1 => +PI/4, -1 => -PI/4)
3163     $MBI->_div($x->{_m}, $MBI->_new(4));
3164     return $x;
3165     }
3166
3167   # This series is only valid if -1 < x < 1, so for other x we need to
3168   # to calculate PI/2 - atan(1/x):
3169   my $one = $MBI->_new(1);
3170   my $pi = undef;
3171   if ($x->{_es} eq '+' && ($MBI->_acmp($x->{_m},$one) >= 0))
3172     {
3173     # calculate PI/2
3174     $pi = $self->bpi($scale - 3);
3175     $MBI->_div($pi->{_m}, $MBI->_new(2));
3176     # calculate 1/$x:
3177     my $x_copy = $x->copy();
3178     # modify $x in place
3179     $x->bone(); $x->bdiv($x_copy,$scale);
3180     }
3181
3182   # when user set globals, they would interfere with our calculation, so
3183   # disable them and later re-enable them
3184   no strict 'refs';
3185   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
3186   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
3187   # we also need to disable any set A or P on $x (_find_round_parameters took
3188   # them already into account), since these would interfere, too
3189   delete $x->{_a}; delete $x->{_p};
3190   # need to disable $upgrade in BigInt, to avoid deep recursion
3191   local $Math::BigInt::upgrade = undef;
3192
3193   my $last = 0;
3194   my $over = $x * $x;                   # X ^ 2
3195   my $x2 = $over->copy();               # X ^ 2; difference between terms
3196   $over->bmul($x);                      # X ^ 3 as starting value
3197   my $sign = 1;                         # start with -=
3198   my $below = $self->new(3);
3199   my $two = $self->new(2);
3200   delete $x->{_a}; delete $x->{_p};
3201
3202   my $limit = $self->new("1E-". ($scale-1));
3203   #my $steps = 0;
3204   while (3 < 5)
3205     {
3206     # we calculate the next term, and add it to the last
3207     # when the next term is below our limit, it won't affect the outcome
3208     # anymore, so we stop:
3209     my $next = $over->copy()->bdiv($below,$scale);
3210     last if $next->bacmp($limit) <= 0;
3211
3212     if ($sign == 0)
3213       {
3214       $x->badd($next);
3215       }
3216     else
3217       {
3218       $x->bsub($next);
3219       }
3220     $sign = 1-$sign;                                    # alternate
3221     # calculate things for the next term
3222     $over->bmul($x2);                                   # $x*$x
3223     $below->badd($two);                                 # n += 2
3224     }
3225
3226   if (defined $pi)
3227     {
3228     my $x_copy = $x->copy();
3229     # modify $x in place
3230     $x->{_m} = $pi->{_m};
3231     $x->{_e} = $pi->{_e};
3232     $x->{_es} = $pi->{_es};
3233     # PI/2 - $x
3234     $x->bsub($x_copy);
3235     }
3236
3237   # shortcut to not run through _find_round_parameters again
3238   if (defined $params[0])
3239     {
3240     $x->bround($params[0],$params[2]);          # then round accordingly
3241     }
3242   else
3243     {
3244     $x->bfround($params[1],$params[2]);         # then round accordingly
3245     }
3246   if ($fallback)
3247     {
3248     # clear a/p after round, since user did not request it
3249     delete $x->{_a}; delete $x->{_p};
3250     }
3251   # restore globals
3252   $$abr = $ab; $$pbr = $pb;
3253   $x;
3254   }
3255
3256 ###############################################################################
3257 # rounding functions
3258
3259 sub bfround
3260   {
3261   # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
3262   # $n == 0 means round to integer
3263   # expects and returns normalized numbers!
3264   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
3265
3266   my ($scale,$mode) = $x->_scale_p(@_);
3267   return $x if !defined $scale || $x->modify('bfround'); # no-op
3268
3269   # never round a 0, +-inf, NaN
3270   if ($x->is_zero())
3271     {
3272     $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
3273     return $x;
3274     }
3275   return $x if $x->{sign} !~ /^[+-]$/;
3276
3277   # don't round if x already has lower precision
3278   return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
3279
3280   $x->{_p} = $scale;                    # remember round in any case
3281   delete $x->{_a};                      # and clear A
3282   if ($scale < 0)
3283     {
3284     # round right from the '.'
3285
3286     return $x if $x->{_es} eq '+';              # e >= 0 => nothing to round
3287
3288     $scale = -$scale;                           # positive for simplicity
3289     my $len = $MBI->_len($x->{_m});             # length of mantissa
3290
3291     # the following poses a restriction on _e, but if _e is bigger than a
3292     # scalar, you got other problems (memory etc) anyway
3293     my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e})));   # digits after dot
3294     my $zad = 0;                                # zeros after dot
3295     $zad = $dad - $len if (-$dad < -$len);      # for 0.00..00xxx style
3296
3297     # print "scale $scale dad $dad zad $zad len $len\n";
3298     # number  bsstr   len zad dad
3299     # 0.123   123e-3    3   0 3
3300     # 0.0123  123e-4    3   1 4
3301     # 0.001   1e-3      1   2 3
3302     # 1.23    123e-2    3   0 2
3303     # 1.2345  12345e-4  5   0 4
3304
3305     # do not round after/right of the $dad
3306     return $x if $scale > $dad;                 # 0.123, scale >= 3 => exit
3307
3308     # round to zero if rounding inside the $zad, but not for last zero like:
3309     # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
3310     return $x->bzero() if $scale < $zad;
3311     if ($scale == $zad)                 # for 0.006, scale -3 and trunc
3312       {
3313       $scale = -$len;
3314       }
3315     else
3316       {
3317       # adjust round-point to be inside mantissa
3318       if ($zad != 0)
3319         {
3320         $scale = $scale-$zad;
3321         }
3322       else
3323         {
3324         my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;    # digits before dot
3325         $scale = $dbd+$scale;
3326         }
3327       }
3328     }
3329   else
3330     {
3331     # round left from the '.'
3332
3333     # 123 => 100 means length(123) = 3 - $scale (2) => 1
3334
3335     my $dbt = $MBI->_len($x->{_m});
3336     # digits before dot
3337     my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
3338     # should be the same, so treat it as this
3339     $scale = 1 if $scale == 0;
3340     # shortcut if already integer
3341     return $x if $scale == 1 && $dbt <= $dbd;
3342     # maximum digits before dot
3343     ++$dbd;
3344
3345     if ($scale > $dbd)
3346        {
3347        # not enough digits before dot, so round to zero
3348        return $x->bzero;
3349        }
3350     elsif ( $scale == $dbd )
3351        {
3352        # maximum
3353        $scale = -$dbt;
3354        }
3355     else
3356        {
3357        $scale = $dbd - $scale;
3358        }
3359     }
3360   # pass sign to bround for rounding modes '+inf' and '-inf'
3361   my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
3362   $m->bround($scale,$mode);
3363   $x->{_m} = $m->{value};                       # get our mantissa back
3364   $x->bnorm();
3365   }
3366
3367 sub bround
3368   {
3369   # accuracy: preserve $N digits, and overwrite the rest with 0's
3370   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
3371
3372   if (($_[0] || 0) < 0)
3373     {
3374     require Carp; Carp::croak ('bround() needs positive accuracy');
3375     }
3376
3377   my ($scale,$mode) = $x->_scale_a(@_);
3378   return $x if !defined $scale || $x->modify('bround'); # no-op
3379
3380   # scale is now either $x->{_a}, $accuracy, or the user parameter
3381   # test whether $x already has lower accuracy, do nothing in this case
3382   # but do round if the accuracy is the same, since a math operation might
3383   # want to round a number with A=5 to 5 digits afterwards again
3384   return $x if defined $x->{_a} && $x->{_a} < $scale;
3385
3386   # scale < 0 makes no sense
3387   # scale == 0 => keep all digits
3388   # never round a +-inf, NaN
3389   return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
3390
3391   # 1: never round a 0
3392   # 2: if we should keep more digits than the mantissa has, do nothing
3393   if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
3394     {
3395     $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
3396     return $x;
3397     }
3398
3399   # pass sign to bround for '+inf' and '-inf' rounding modes
3400   my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
3401
3402   $m->bround($scale,$mode);             # round mantissa
3403   $x->{_m} = $m->{value};               # get our mantissa back
3404   $x->{_a} = $scale;                    # remember rounding
3405   delete $x->{_p};                      # and clear P
3406   $x->bnorm();                          # del trailing zeros gen. by bround()
3407   }
3408
3409 sub bfloor
3410   {
3411   # round towards minus infinity
3412   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3413
3414   return $x if $x->modify('bfloor');
3415
3416   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3417
3418   # if $x has digits after dot
3419   if ($x->{_es} eq '-')
3420     {
3421     $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3422     $x->{_e} = $MBI->_zero();                   # trunc/norm
3423     $x->{_es} = '+';                            # abs e
3424     $MBI->_inc($x->{_m}) if $x->{sign} eq '-';  # increment if negative
3425     }
3426   $x->round($a,$p,$r);
3427   }
3428
3429 sub bceil
3430   {
3431   # round towards plus infinity
3432   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3433
3434   return $x if $x->modify('bceil');
3435   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3436
3437   # if $x has digits after dot
3438   if ($x->{_es} eq '-')
3439     {
3440     $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3441     $x->{_e} = $MBI->_zero();                   # trunc/norm
3442     $x->{_es} = '+';                            # abs e
3443     $MBI->_inc($x->{_m}) if $x->{sign} eq '+';  # increment if positive
3444     }
3445   $x->round($a,$p,$r);
3446   }
3447
3448 sub bint
3449   {
3450   # round towards zero
3451   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3452
3453   return $x if $x->modify('bint');
3454   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3455
3456   # if $x has digits after the decimal point
3457   if ($x->{_es} eq '-')
3458     {
3459     $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3460     $x->{_e} = $MBI->_zero();                     # truncate/normalize
3461     $x->{_es} = '+';                              # abs e
3462     }
3463   $x->round($a,$p,$r);
3464   }
3465
3466 sub brsft
3467   {
3468   # shift right by $y (divide by power of $n)
3469
3470   # set up parameters
3471   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3472   # objectify is costly, so avoid it
3473   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3474     {
3475     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3476     }
3477
3478   return $x if $x->modify('brsft');
3479   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3480
3481   $n = 2 if !defined $n; $n = $self->new($n);
3482
3483   # negative amount?
3484   return $x->blsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3485
3486   # the following call to bdiv() will return either quo or (quo,remainder):
3487   $x->bdiv($n->bpow($y),$a,$p,$r,$y);
3488   }
3489
3490 sub blsft
3491   {
3492   # shift left by $y (multiply by power of $n)
3493
3494   # set up parameters
3495   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3496   # objectify is costly, so avoid it
3497   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3498     {
3499     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3500     }
3501
3502   return $x if $x->modify('blsft');
3503   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3504
3505   $n = 2 if !defined $n; $n = $self->new($n);
3506
3507   # negative amount?
3508   return $x->brsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3509
3510   $x->bmul($n->bpow($y),$a,$p,$r,$y);
3511   }
3512
3513 ###############################################################################
3514
3515 sub DESTROY
3516   {
3517   # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
3518   }
3519
3520 sub AUTOLOAD
3521   {
3522   # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
3523   # or falling back to MBI::bxxx()
3524   my $name = $AUTOLOAD;
3525
3526   $name =~ s/(.*):://;  # split package
3527   my $c = $1 || $class;
3528   no strict 'refs';
3529   $c->import() if $IMPORT == 0;
3530   if (!_method_alias($name))
3531     {
3532     if (!defined $name)
3533       {
3534       # delayed load of Carp and avoid recursion
3535       require Carp;
3536       Carp::croak ("$c: Can't call a method without name");
3537       }
3538     if (!_method_hand_up($name))
3539       {
3540       # delayed load of Carp and avoid recursion
3541       require Carp;
3542       Carp::croak ("Can't call $c\-\>$name, not a valid method");
3543       }
3544     # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
3545     $name =~ s/^f/b/;
3546     return &{"Math::BigInt"."::$name"}(@_);
3547     }
3548   my $bname = $name; $bname =~ s/^f/b/;
3549   $c .= "::$name";
3550   *{$c} = \&{$bname};
3551   &{$c};        # uses @_
3552   }
3553
3554 sub exponent
3555   {
3556   # return a copy of the exponent
3557   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3558
3559   if ($x->{sign} !~ /^[+-]$/)
3560     {
3561     my $s = $x->{sign}; $s =~ s/^[+-]//;
3562     return Math::BigInt->new($s);               # -inf, +inf => +inf
3563     }
3564   Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
3565   }
3566
3567 sub mantissa
3568   {
3569   # return a copy of the mantissa
3570   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3571
3572   if ($x->{sign} !~ /^[+-]$/)
3573     {
3574     my $s = $x->{sign}; $s =~ s/^[+]//;
3575     return Math::BigInt->new($s);               # -inf, +inf => +inf
3576     }
3577   my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
3578   $m->bneg() if $x->{sign} eq '-';
3579
3580   $m;
3581   }
3582
3583 sub parts
3584   {
3585   # return a copy of both the exponent and the mantissa
3586   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3587
3588   if ($x->{sign} !~ /^[+-]$/)
3589     {
3590     my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
3591     return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
3592     }
3593   my $m = Math::BigInt->bzero();
3594   $m->{value} = $MBI->_copy($x->{_m});
3595   $m->bneg() if $x->{sign} eq '-';
3596   ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
3597   }
3598
3599 ##############################################################################
3600 # private stuff (internal use only)
3601
3602 sub import
3603   {
3604   my $self = shift;
3605   my $l = scalar @_;
3606   my $lib = ''; my @a;
3607   my $lib_kind = 'try';
3608   $IMPORT=1;
3609   for ( my $i = 0; $i < $l ; $i++)
3610     {
3611     if ( $_[$i] eq ':constant' )
3612       {
3613       # This causes overlord er load to step in. 'binary' and 'integer'
3614       # are handled by BigInt.
3615       overload::constant float => sub { $self->new(shift); };
3616       }
3617     elsif ($_[$i] eq 'upgrade')
3618       {
3619       # this causes upgrading
3620       $upgrade = $_[$i+1];              # or undef to disable
3621       $i++;
3622       }
3623     elsif ($_[$i] eq 'downgrade')
3624       {
3625       # this causes downgrading
3626       $downgrade = $_[$i+1];            # or undef to disable
3627       $i++;
3628       }
3629     elsif ($_[$i] =~ /^(lib|try|only)\z/)
3630       {
3631       # alternative library
3632       $lib = $_[$i+1] || '';            # default Calc
3633       $lib_kind = $1;                   # lib, try or only
3634       $i++;
3635       }
3636     elsif ($_[$i] eq 'with')
3637       {
3638       # alternative class for our private parts()
3639       # XXX: no longer supported
3640       # $MBI = $_[$i+1] || 'Math::BigInt';
3641       $i++;
3642       }
3643     else
3644       {
3645       push @a, $_[$i];
3646       }
3647     }
3648
3649   $lib =~ tr/a-zA-Z0-9,://cd;           # restrict to sane characters
3650   # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
3651   my $mbilib = eval { Math::BigInt->config()->{lib} };
3652   if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
3653     {
3654     # MBI already loaded
3655     Math::BigInt->import( $lib_kind, "$lib,$mbilib", 'objectify');
3656     }
3657   else
3658     {
3659     # MBI not loaded, or with ne "Math::BigInt::Calc"
3660     $lib .= ",$mbilib" if defined $mbilib;
3661     $lib =~ s/^,//;                             # don't leave empty
3662
3663     # replacement library can handle lib statement, but also could ignore it
3664
3665     # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
3666     # used in the same script, or eval inside import(). So we require MBI:
3667     require Math::BigInt;
3668     Math::BigInt->import( $lib_kind => $lib, 'objectify' );
3669     }
3670   if ($@)
3671     {
3672     require Carp; Carp::croak ("Couldn't load $lib: $! $@");
3673     }
3674   # find out which one was actually loaded
3675   $MBI = Math::BigInt->config()->{lib};
3676
3677   # register us with MBI to get notified of future lib changes
3678   Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
3679
3680   $self->export_to_level(1,$self,@a);           # export wanted functions
3681   }
3682
3683 sub bnorm
3684   {
3685   # adjust m and e so that m is smallest possible
3686   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
3687
3688   return $x if $x->{sign} !~ /^[+-]$/;          # inf, nan etc
3689
3690   my $zeros = $MBI->_zeros($x->{_m});           # correct for trailing zeros
3691   if ($zeros != 0)
3692     {
3693     my $z = $MBI->_new($zeros);
3694     $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
3695     if ($x->{_es} eq '-')
3696       {
3697       if ($MBI->_acmp($x->{_e},$z) >= 0)
3698         {
3699         $x->{_e} = $MBI->_sub ($x->{_e}, $z);
3700         $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
3701         }
3702       else
3703         {
3704         $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
3705         $x->{_es} = '+';
3706         }
3707       }
3708     else
3709       {
3710       $x->{_e} = $MBI->_add ($x->{_e}, $z);
3711       }
3712     }
3713   else
3714     {
3715     # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
3716     # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
3717     $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
3718      if $MBI->_is_zero($x->{_m});
3719     }
3720
3721   $x;                                   # MBI bnorm is no-op, so do not call it
3722   }
3723
3724 ##############################################################################
3725
3726 sub as_hex
3727   {
3728   # return number as hexadecimal string (only for integers defined)
3729   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3730
3731   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
3732   return '0x0' if $x->is_zero();
3733
3734   return $nan if $x->{_es} ne '+';              # how to do 1e-1 in hex!?
3735
3736   my $z = $MBI->_copy($x->{_m});
3737   if (! $MBI->_is_zero($x->{_e}))               # > 0
3738     {
3739     $MBI->_lsft( $z, $x->{_e},10);
3740     }
3741   $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3742   $z->as_hex();
3743   }
3744
3745 sub as_bin
3746   {
3747   # return number as binary digit string (only for integers defined)
3748   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3749
3750   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
3751   return '0b0' if $x->is_zero();
3752
3753   return $nan if $x->{_es} ne '+';              # how to do 1e-1 in hex!?
3754
3755   my $z = $MBI->_copy($x->{_m});
3756   if (! $MBI->_is_zero($x->{_e}))               # > 0
3757     {
3758     $MBI->_lsft( $z, $x->{_e},10);
3759     }
3760   $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3761   $z->as_bin();
3762   }
3763
3764 sub as_oct
3765   {
3766   # return number as octal digit string (only for integers defined)
3767   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3768
3769   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
3770   return '0' if $x->is_zero();
3771
3772   return $nan if $x->{_es} ne '+';              # how to do 1e-1 in hex!?
3773
3774   my $z = $MBI->_copy($x->{_m});
3775   if (! $MBI->_is_zero($x->{_e}))               # > 0
3776     {
3777     $MBI->_lsft( $z, $x->{_e},10);
3778     }
3779   $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3780   $z->as_oct();
3781   }
3782
3783 sub as_number
3784   {
3785   # return copy as a bigint representation of this BigFloat number
3786   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3787
3788   return $x if $x->modify('as_number');
3789
3790   if (!$x->isa('Math::BigFloat'))
3791     {
3792     # if the object can as_number(), use it
3793     return $x->as_number() if $x->can('as_number');
3794     # otherwise, get us a float and then a number
3795     $x = $x->can('as_float') ? $x->as_float() : $self->new(0+"$x");
3796     }
3797
3798   return Math::BigInt->binf($x->sign()) if $x->is_inf();
3799   return Math::BigInt->bnan()           if $x->is_nan();
3800
3801   my $z = $MBI->_copy($x->{_m});
3802   if ($x->{_es} eq '-')                 # < 0
3803     {
3804     $MBI->_rsft( $z, $x->{_e},10);
3805     }
3806   elsif (! $MBI->_is_zero($x->{_e}))    # > 0
3807     {
3808     $MBI->_lsft( $z, $x->{_e},10);
3809     }
3810   $z = Math::BigInt->new( $x->{sign} . $MBI->_str($z));
3811   $z;
3812   }
3813
3814 sub length
3815   {
3816   my $x = shift;
3817   my $class = ref($x) || $x;
3818   $x = $class->new(shift) unless ref($x);
3819
3820   return 1 if $MBI->_is_zero($x->{_m});
3821
3822   my $len = $MBI->_len($x->{_m});
3823   $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
3824   if (wantarray())
3825     {
3826     my $t = 0;
3827     $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
3828     return ($len, $t);
3829     }
3830   $len;
3831   }
3832
3833 1;
3834
3835 __END__
3836
3837 =pod
3838
3839 =head1 NAME
3840
3841 Math::BigFloat - Arbitrary size floating point math package
3842
3843 =head1 SYNOPSIS
3844
3845  use Math::BigFloat;
3846
3847  # Number creation
3848  my $x = Math::BigFloat->new($str);     # defaults to 0
3849  my $y = $x->copy();                    # make a true copy
3850  my $nan  = Math::BigFloat->bnan();     # create a NotANumber
3851  my $zero = Math::BigFloat->bzero();    # create a +0
3852  my $inf = Math::BigFloat->binf();      # create a +inf
3853  my $inf = Math::BigFloat->binf('-');   # create a -inf
3854  my $one = Math::BigFloat->bone();      # create a +1
3855  my $mone = Math::BigFloat->bone('-');  # create a -1
3856
3857  my $pi = Math::BigFloat->bpi(100);     # PI to 100 digits
3858
3859  # the following examples compute their result to 100 digits accuracy:
3860  my $cos  = Math::BigFloat->new(1)->bcos(100);        # cosinus(1)
3861  my $sin  = Math::BigFloat->new(1)->bsin(100);        # sinus(1)
3862  my $atan = Math::BigFloat->new(1)->batan(100);       # arcus tangens(1)
3863
3864  my $atan2 = Math::BigFloat->new(  1 )->batan2( 1 ,100); # batan(1)
3865  my $atan2 = Math::BigFloat->new(  1 )->batan2( 8 ,100); # batan(1/8)
3866  my $atan2 = Math::BigFloat->new( -2 )->batan2( 1 ,100); # batan(-2)
3867
3868  # Testing
3869  $x->is_zero();          # true if arg is +0
3870  $x->is_nan();           # true if arg is NaN
3871  $x->is_one();           # true if arg is +1
3872  $x->is_one('-');        # true if arg is -1
3873  $x->is_odd();           # true if odd, false for even
3874  $x->is_even();          # true if even, false for odd
3875  $x->is_pos();           # true if >= 0
3876  $x->is_neg();           # true if <  0
3877  $x->is_inf(sign);       # true if +inf, or -inf (default is '+')
3878
3879  $x->bcmp($y);           # compare numbers (undef,<0,=0,>0)
3880  $x->bacmp($y);          # compare absolutely (undef,<0,=0,>0)
3881  $x->sign();             # return the sign, either +,- or NaN
3882  $x->digit($n);          # return the nth digit, counting from right
3883  $x->digit(-$n);         # return the nth digit, counting from left
3884
3885  # The following all modify their first argument. If you want to pre-
3886  # serve $x, use $z = $x->copy()->bXXX($y); See under L</CAVEATS> for
3887  # necessary when mixing $a = $b assignments with non-overloaded math.
3888
3889  # set
3890  $x->bzero();            # set $i to 0
3891  $x->bnan();             # set $i to NaN
3892  $x->bone();             # set $x to +1
3893  $x->bone('-');          # set $x to -1
3894  $x->binf();             # set $x to inf
3895  $x->binf('-');          # set $x to -inf
3896
3897  $x->bneg();             # negation
3898  $x->babs();             # absolute value
3899  $x->bnorm();            # normalize (no-op)
3900  $x->bnot();             # two's complement (bit wise not)
3901  $x->binc();             # increment x by 1
3902  $x->bdec();             # decrement x by 1
3903
3904  $x->badd($y);           # addition (add $y to $x)
3905  $x->bsub($y);           # subtraction (subtract $y from $x)
3906  $x->bmul($y);           # multiplication (multiply $x by $y)
3907  $x->bdiv($y);           # divide, set $x to quotient
3908                          # return (quo,rem) or quo if scalar
3909
3910  $x->bmod($y);           # modulus ($x % $y)
3911  $x->bpow($y);           # power of arguments ($x ** $y)
3912  $x->bmodpow($exp,$mod); # modular exponentiation (($num**$exp) % $mod))
3913  $x->blsft($y, $n);      # left shift by $y places in base $n
3914  $x->brsft($y, $n);      # right shift by $y places in base $n
3915                          # returns (quo,rem) or quo if in scalar context
3916
3917  $x->blog();             # logarithm of $x to base e (Euler's number)
3918  $x->blog($base);        # logarithm of $x to base $base (f.i. 2)
3919  $x->bexp();             # calculate e ** $x where e is Euler's number
3920
3921  $x->band($y);           # bit-wise and
3922  $x->bior($y);           # bit-wise inclusive or
3923  $x->bxor($y);           # bit-wise exclusive or
3924  $x->bnot();             # bit-wise not (two's complement)
3925
3926  $x->bsqrt();            # calculate square-root
3927  $x->broot($y);          # $y'th root of $x (e.g. $y == 3 => cubic root)
3928  $x->bfac();             # factorial of $x (1*2*3*4*..$x)
3929
3930  $x->bround($N);         # accuracy: preserve $N digits
3931  $x->bfround($N);        # precision: round to the $Nth digit
3932
3933  $x->bfloor();           # return integer less or equal than $x
3934  $x->bceil();            # return integer greater or equal than $x
3935  $x->bint();             # round towards zero
3936
3937   # The following do not modify their arguments:
3938
3939  bgcd(@values);          # greatest common divisor
3940  blcm(@values);          # lowest common multiplicator
3941
3942  $x->bstr();             # return string
3943  $x->bsstr();            # return string in scientific notation
3944
3945  $x->as_int();           # return $x as BigInt
3946  $x->exponent();         # return exponent as BigInt
3947  $x->mantissa();         # return mantissa as BigInt
3948  $x->parts();            # return (mantissa,exponent) as BigInt
3949
3950  $x->length();           # number of digits (w/o sign and '.')
3951  ($l,$f) = $x->length(); # number of digits, and length of fraction
3952
3953  $x->precision();        # return P of $x (or global, if P of $x undef)
3954  $x->precision($n);      # set P of $x to $n
3955  $x->accuracy();         # return A of $x (or global, if A of $x undef)
3956  $x->accuracy($n);       # set A $x to $n
3957
3958  # these get/set the appropriate global value for all BigFloat objects
3959  Math::BigFloat->precision();   # Precision
3960  Math::BigFloat->accuracy();    # Accuracy
3961  Math::BigFloat->round_mode();  # rounding mode
3962
3963 =head1 DESCRIPTION
3964
3965 All operators (including basic math operations) are overloaded if you
3966 declare your big floating point numbers as
3967
3968   $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
3969
3970 Operations with overloaded operators preserve the arguments, which is
3971 exactly what you expect.
3972
3973 =head2 Input
3974
3975 Input to these routines are either BigFloat objects, or strings of the
3976 following four forms:
3977
3978 =over
3979
3980 =item *
3981
3982 C</^[+-]\d+$/>
3983
3984 =item *
3985
3986 C</^[+-]\d+\.\d*$/>
3987
3988 =item *
3989
3990 C</^[+-]\d+E[+-]?\d+$/>
3991
3992 =item *
3993
3994 C</^[+-]\d*\.\d+E[+-]?\d+$/>
3995
3996 =back
3997
3998 all with optional leading and trailing zeros and/or spaces. Additionally,
3999 numbers are allowed to have an underscore between any two digits.
4000
4001 Empty strings as well as other illegal numbers results in 'NaN'.
4002
4003 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
4004 are always stored in normalized form. On a string, it creates a BigFloat
4005 object.
4006
4007 =head2 Output
4008
4009 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
4010
4011 The string output will always have leading and trailing zeros stripped and drop
4012 a plus sign. C<bstr()> will give you always the form with a decimal point,
4013 while C<bsstr()> (s for scientific) gives you the scientific notation.
4014
4015         Input                   bstr()          bsstr()
4016         '-0'                    '0'             '0E1'
4017         '  -123 123 123'        '-123123123'    '-123123123E0'
4018         '00.0123'               '0.0123'        '123E-4'
4019         '123.45E-2'             '1.2345'        '12345E-4'
4020         '10E+3'                 '10000'         '1E4'
4021
4022 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
4023 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
4024 return either undef, <0, 0 or >0 and are suited for sort.
4025
4026 Actual math is done by using the class defined with C<< with => Class; >>
4027 (which defaults to BigInts) to represent the mantissa and exponent.
4028
4029 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
4030 represent the result when input arguments are not numbers, as well as
4031 the result of dividing by zero.
4032
4033 =head2 mantissa(), exponent() and parts()
4034
4035 mantissa() and exponent() return the said parts of the BigFloat
4036 as BigInts such that:
4037
4038         $m = $x->mantissa();
4039         $e = $x->exponent();
4040         $y = $m * ( 10 ** $e );
4041         print "ok\n" if $x == $y;
4042
4043 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
4044
4045 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
4046
4047 Currently the mantissa is reduced as much as possible, favouring higher
4048 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
4049 This might change in the future, so do not depend on it.
4050
4051 =head2 Accuracy vs. Precision
4052
4053 See also: L<Rounding|/Rounding>.
4054
4055 Math::BigFloat supports both precision (rounding to a certain place before or
4056 after the dot) and accuracy (rounding to a certain number of digits). For a
4057 full documentation, examples and tips on these topics please see the large
4058 section about rounding in L<Math::BigInt>.
4059
4060 Since things like C<sqrt(2)> or C<1 / 3> must presented with a limited
4061 accuracy lest a operation consumes all resources, each operation produces
4062 no more than the requested number of digits.
4063
4064 If there is no global precision or accuracy set, B<and> the operation in
4065 question was not called with a requested precision or accuracy, B<and> the
4066 input $x has no accuracy or precision set, then a fallback parameter will
4067 be used. For historical reasons, it is called C<div_scale> and can be accessed
4068 via:
4069
4070         $d = Math::BigFloat->div_scale();       # query
4071         Math::BigFloat->div_scale($n);          # set to $n digits
4072
4073 The default value for C<div_scale> is 40.
4074
4075 In case the result of one operation has more digits than specified,
4076 it is rounded. The rounding mode taken is either the default mode, or the one
4077 supplied to the operation after the I<scale>:
4078
4079     $x = Math::BigFloat->new(2);
4080     Math::BigFloat->accuracy(5);              # 5 digits max
4081     $y = $x->copy()->bdiv(3);                 # will give 0.66667
4082     $y = $x->copy()->bdiv(3,6);               # will give 0.666667
4083     $y = $x->copy()->bdiv(3,6,undef,'odd');   # will give 0.666667
4084     Math::BigFloat->round_mode('zero');
4085     $y = $x->copy()->bdiv(3,6);               # will also give 0.666667
4086
4087 Note that C<< Math::BigFloat->accuracy() >> and C<< Math::BigFloat->precision() >>
4088 set the global variables, and thus B<any> newly created number will be subject
4089 to the global rounding B<immediately>. This means that in the examples above, the
4090 C<3> as argument to C<bdiv()> will also get an accuracy of B<5>.
4091
4092 It is less confusing to either calculate the result fully, and afterwards
4093 round it explicitly, or use the additional parameters to the math
4094 functions like so:
4095
4096         use Math::BigFloat;
4097         $x = Math::BigFloat->new(2);
4098         $y = $x->copy()->bdiv(3);
4099         print $y->bround(5),"\n";               # will give 0.66667
4100
4101         or
4102
4103         use Math::BigFloat;
4104         $x = Math::BigFloat->new(2);
4105         $y = $x->copy()->bdiv(3,5);             # will give 0.66667
4106         print "$y\n";
4107
4108 =head2 Rounding
4109
4110 =over
4111
4112 =item ffround ( +$scale )
4113
4114 Rounds to the $scale'th place left from the '.', counting from the dot.
4115 The first digit is numbered 1.
4116
4117 =item ffround ( -$scale )
4118
4119 Rounds to the $scale'th place right from the '.', counting from the dot.
4120
4121 =item ffround ( 0 )
4122
4123 Rounds to an integer.
4124
4125 =item fround  ( +$scale )
4126
4127 Preserves accuracy to $scale digits from the left (aka significant digits)
4128 and pads the rest with zeros. If the number is between 1 and -1, the
4129 significant digits count from the first non-zero after the '.'
4130
4131 =item fround  ( -$scale ) and fround ( 0 )
4132
4133 These are effectively no-ops.
4134
4135 =back
4136
4137 All rounding functions take as a second parameter a rounding mode from one of
4138 the following: 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'.
4139
4140 The default rounding mode is 'even'. By using
4141 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
4142 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
4143 no longer supported.
4144 The second parameter to the round functions then overrides the default
4145 temporarily.
4146
4147 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
4148 'trunc' as rounding mode to make it equivalent to:
4149
4150         $x = 2.5;
4151         $y = int($x) + 2;
4152
4153 You can override this by passing the desired rounding mode as parameter to
4154 C<as_number()>:
4155
4156         $x = Math::BigFloat->new(2.5);
4157         $y = $x->as_number('odd');      # $y = 3
4158
4159 =head1 METHODS
4160
4161 Math::BigFloat supports all methods that Math::BigInt supports, except it
4162 calculates non-integer results when possible. Please see L<Math::BigInt>
4163 for a full description of each method. Below are just the most important
4164 differences:
4165
4166 =over
4167
4168 =item accuracy()
4169
4170       $x->accuracy(5);           # local for $x
4171       CLASS->accuracy(5);        # global for all members of CLASS
4172                                  # Note: This also applies to new()!
4173
4174       $A = $x->accuracy();       # read out accuracy that affects $x
4175       $A = CLASS->accuracy();    # read out global accuracy
4176
4177 Set or get the global or local accuracy, aka how many significant digits the
4178 results have. If you set a global accuracy, then this also applies to new()!
4179
4180 Warning! The accuracy I<sticks>, e.g. once you created a number under the
4181 influence of C<< CLASS->accuracy($A) >>, all results from math operations with
4182 that number will also be rounded.
4183
4184 In most cases, you should probably round the results explicitly using one of
4185 L<Math::BigInt/round()>, L<Math::BigInt/bround()> or L<Math::BigInt/bfround()> or by passing the desired accuracy
4186 to the math operation as additional parameter:
4187
4188         my $x = Math::BigInt->new(30000);
4189         my $y = Math::BigInt->new(7);
4190         print scalar $x->copy()->bdiv($y, 2);           # print 4300
4191         print scalar $x->copy()->bdiv($y)->bround(2);   # print 4300
4192
4193 =item precision()
4194
4195       $x->precision(-2);      # local for $x, round at the second
4196                               # digit right of the dot
4197       $x->precision(2);       # ditto, round at the second digit
4198                               # left of the dot
4199
4200       CLASS->precision(5);    # Global for all members of CLASS
4201                               # This also applies to new()!
4202       CLASS->precision(-5);   # ditto
4203
4204       $P = CLASS->precision();  # read out global precision
4205       $P = $x->precision();     # read out precision that affects $x
4206
4207 Note: You probably want to use L</accuracy()> instead. With L</accuracy()> you
4208 set the number of digits each result should have, with L</precision()> you
4209 set the place where to round!
4210
4211 =item bexp()
4212
4213         $x->bexp($accuracy);            # calculate e ** X
4214
4215 Calculates the expression C<e ** $x> where C<e> is Euler's number.
4216
4217 This method was added in v1.82 of Math::BigInt (April 2007).
4218
4219 =item bnok()
4220
4221         $x->bnok($y);   # x over y (binomial coefficient n over k)
4222
4223 Calculates the binomial coefficient n over k, also called the "choose"
4224 function. The result is equivalent to:
4225
4226         ( n )      n!
4227         | - |  = -------
4228         ( k )    k!(n-k)!
4229
4230 This method was added in v1.84 of Math::BigInt (April 2007).
4231
4232 =item bpi()
4233
4234         print Math::BigFloat->bpi(100), "\n";
4235
4236 Calculate PI to N digits (including the 3 before the dot). The result is
4237 rounded according to the current rounding mode, which defaults to "even".
4238
4239 This method was added in v1.87 of Math::BigInt (June 2007).
4240
4241 =item bcos()
4242
4243         my $x = Math::BigFloat->new(1);
4244         print $x->bcos(100), "\n";
4245
4246 Calculate the cosinus of $x, modifying $x in place.
4247
4248 This method was added in v1.87 of Math::BigInt (June 2007).
4249
4250 =item bsin()
4251
4252         my $x = Math::BigFloat->new(1);
4253         print $x->bsin(100), "\n";
4254
4255 Calculate the sinus of $x, modifying $x in place.
4256
4257 This method was added in v1.87 of Math::BigInt (June 2007).
4258
4259 =item batan2()
4260
4261         my $y = Math::BigFloat->new(2);
4262         my $x = Math::BigFloat->new(3);
4263         print $y->batan2($x), "\n";
4264
4265 Calculate the arcus tanges of C<$y> divided by C<$x>, modifying $y in place.
4266 See also L</batan()>.
4267
4268 This method was added in v1.87 of Math::BigInt (June 2007).
4269
4270 =item batan()
4271
4272         my $x = Math::BigFloat->new(1);
4273         print $x->batan(100), "\n";
4274
4275 Calculate the arcus tanges of $x, modifying $x in place. See also L</batan2()>.
4276
4277 This method was added in v1.87 of Math::BigInt (June 2007).
4278
4279 =item bmuladd()
4280
4281         $x->bmuladd($y,$z);
4282
4283 Multiply $x by $y, and then add $z to the result.
4284
4285 This method was added in v1.87 of Math::BigInt (June 2007).
4286
4287 =back
4288
4289 =head1 Autocreating constants
4290
4291 After C<use Math::BigFloat ':constant'> all the floating point constants
4292 in the given scope are converted to C<Math::BigFloat>. This conversion
4293 happens at compile time.
4294
4295 In particular
4296
4297   perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
4298
4299 prints the value of C<2E-100>. Note that without conversion of
4300 constants the expression 2E-100 will be calculated as normal floating point
4301 number.
4302
4303 Please note that ':constant' does not affect integer constants, nor binary
4304 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
4305 work.
4306
4307 =head2 Math libr