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