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