#line 1 "Math/BigFloat.pm"
package Math::BigFloat;

# 
# Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
#

# The following hash values are internally used:
#   _e	: exponent (ref to $CALC object)
#   _m	: mantissa (ref to $CALC object)
#   _es	: sign of _e
# sign	: +,-,+inf,-inf, or "NaN" if not a number
#   _a	: accuracy
#   _p	: precision

use 5.006001;
use strict;
use warnings;

use Carp ();

our $VERSION = '1.999722';
$VERSION = eval $VERSION;

require Exporter;
our @ISA	= qw/Math::BigInt/;
our @EXPORT_OK	= qw/bpi/;

# $_trap_inf/$_trap_nan are internal and should never be accessed from outside
our ($AUTOLOAD, $accuracy, $precision, $div_scale, $round_mode, $rnd_mode,
     $upgrade, $downgrade, $_trap_nan, $_trap_inf);

my $class = "Math::BigFloat";

use overload

  # overload key: with_assign

  '+'     =>      sub { $_[0] -> copy() -> badd($_[1]); },

  '-'     =>      sub { my $c = $_[0] -> copy;
                        $_[2] ? $c -> bneg() -> badd( $_[1])
                              : $c -> bsub($_[1]); },

  '*'     =>      sub { $_[0] -> copy() -> bmul($_[1]); },

  '/'     =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bdiv($_[0])
                              : $_[0] -> copy -> bdiv($_[1]); },


  '%'     =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bmod($_[0])
                              : $_[0] -> copy -> bmod($_[1]); },

  '**'    =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bpow($_[0])
                              : $_[0] -> copy -> bpow($_[1]); },

  '<<'    =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> blsft($_[0])
                              : $_[0] -> copy -> blsft($_[1]); },

  '>>'    =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> brsft($_[0])
                              : $_[0] -> copy -> brsft($_[1]); },

  # overload key: assign

  '+='    =>      sub { $_[0]->badd($_[1]); },

  '-='    =>      sub { $_[0]->bsub($_[1]); },

  '*='    =>      sub { $_[0]->bmul($_[1]); },

  '/='    =>      sub { scalar $_[0]->bdiv($_[1]); },

  '%='    =>      sub { $_[0]->bmod($_[1]); },

  '**='   =>      sub { $_[0]->bpow($_[1]); },


  '<<='   =>      sub { $_[0]->blsft($_[1]); },

  '>>='   =>      sub { $_[0]->brsft($_[1]); },

#  'x='    =>      sub { },

#  '.='    =>      sub { },

  # overload key: num_comparison

  '<'     =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> blt($_[0])
                              : $_[0] -> blt($_[1]); },

  '<='    =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> ble($_[0])
                              : $_[0] -> ble($_[1]); },

  '>'     =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bgt($_[0])
                              : $_[0] -> bgt($_[1]); },

  '>='    =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bge($_[0])
                              : $_[0] -> bge($_[1]); },

  '=='    =>      sub { $_[0] -> beq($_[1]); },

  '!='    =>      sub { $_[0] -> bne($_[1]); },

  # overload key: 3way_comparison

  '<=>'   =>      sub { my $cmp = $_[0] -> bcmp($_[1]);
                        defined($cmp) && $_[2] ? -$cmp : $cmp; },

  'cmp'   =>      sub { $_[2] ? "$_[1]" cmp $_[0] -> bstr()
                              : $_[0] -> bstr() cmp "$_[1]"; },

  # overload key: str_comparison

#  'lt'    =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrlt($_[0])
#                              : $_[0] -> bstrlt($_[1]); },
#
#  'le'    =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrle($_[0])
#                              : $_[0] -> bstrle($_[1]); },
#
#  'gt'    =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrgt($_[0])
#                              : $_[0] -> bstrgt($_[1]); },
#
#  'ge'    =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrge($_[0])
#                              : $_[0] -> bstrge($_[1]); },
#
#  'eq'    =>      sub { $_[0] -> bstreq($_[1]); },
#
#  'ne'    =>      sub { $_[0] -> bstrne($_[1]); },

  # overload key: binary

  '&'     =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> band($_[0])
                              : $_[0] -> copy -> band($_[1]); },

  '&='    =>      sub { $_[0] -> band($_[1]); },

  '|'     =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bior($_[0])
                              : $_[0] -> copy -> bior($_[1]); },

  '|='    =>      sub { $_[0] -> bior($_[1]); },

  '^'     =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bxor($_[0])
                              : $_[0] -> copy -> bxor($_[1]); },

  '^='    =>      sub { $_[0] -> bxor($_[1]); },

#  '&.'    =>      sub { },

#  '&.='   =>      sub { },

#  '|.'    =>      sub { },

#  '|.='   =>      sub { },

#  '^.'    =>      sub { },

#  '^.='   =>      sub { },

  # overload key: unary

  'neg'   =>      sub { $_[0] -> copy() -> bneg(); },

#  '!'     =>      sub { },

  '~'     =>      sub { $_[0] -> copy() -> bnot(); },

#  '~.'    =>      sub { },

  # overload key: mutators

  '++'    =>      sub { $_[0] -> binc() },

  '--'    =>      sub { $_[0] -> bdec() },

  # overload key: func

  'atan2' =>      sub { $_[2] ? ref($_[0]) -> new($_[1]) -> batan2($_[0])
                              : $_[0] -> copy() -> batan2($_[1]); },

  'cos'   =>      sub { $_[0] -> copy -> bcos(); },

  'sin'   =>      sub { $_[0] -> copy -> bsin(); },

  'exp'   =>      sub { $_[0] -> copy() -> bexp($_[1]); },

  'abs'   =>      sub { $_[0] -> copy() -> babs(); },

  'log'   =>      sub { $_[0] -> copy() -> blog(); },

  'sqrt'  =>      sub { $_[0] -> copy() -> bsqrt(); },

  'int'   =>      sub { $_[0] -> copy() -> bint(); },

  # overload key: conversion

  'bool'  =>      sub { $_[0] -> is_zero() ? '' : 1; },

  '""'    =>      sub { $_[0] -> bstr(); },

  '0+'    =>      sub { $_[0] -> numify(); },

  '='     =>      sub { $_[0]->copy(); },

  ;

##############################################################################
# global constants, flags and assorted stuff

# the following are public, but their usage is not recommended. Use the
# accessor methods instead.

# class constants, use Class->constant_name() to access
# one of 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'
$round_mode = 'even';
$accuracy   = undef;
$precision  = undef;
$div_scale  = 40;

$upgrade = undef;
$downgrade = undef;
# the package we are using for our private parts, defaults to:
# Math::BigInt->config()->{lib}
my $MBI = 'Math::BigInt::Calc';

# are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
$_trap_nan = 0;
# the same for infinity
$_trap_inf = 0;

# constant for easier life
my $nan = 'NaN'; 

my $IMPORT = 0;	# was import() called yet? used to make require work

# some digits of accuracy for blog(undef,10); which we use in blog() for speed
my $LOG_10 = 
 '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
my $LOG_10_A = length($LOG_10)-1;
# ditto for log(2)
my $LOG_2 = 
 '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
my $LOG_2_A = length($LOG_2)-1;
my $HALF = '0.5';			# made into an object if nec.

##############################################################################
# the old code had $rnd_mode, so we need to support it, too

sub TIESCALAR   { my ($class) = @_; bless \$round_mode, $class; }
sub FETCH       { return $round_mode; }
sub STORE       { $rnd_mode = $_[0]->round_mode($_[1]); }

BEGIN
  {
  # when someone sets $rnd_mode, we catch this and check the value to see
  # whether it is valid or not. 
  $rnd_mode   = 'even'; tie $rnd_mode, 'Math::BigFloat';

  # we need both of them in this package:
  *as_int = \&as_number;
  }
 
##############################################################################

{
  # valid method aliases for AUTOLOAD
  my %methods = map { $_ => 1 }  
   qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
        fint facmp fcmp fzero fnan finf finc fdec ffac fneg
	fceil ffloor frsft flsft fone flog froot fexp
      /;
  # valid methods that can be handed up (for AUTOLOAD)
  my %hand_ups = map { $_ => 1 }  
   qw / is_nan is_inf is_negative is_positive is_pos is_neg
        accuracy precision div_scale round_mode fabs fnot
        objectify upgrade downgrade
	bone binf bnan bzero
	bsub
      /;

  sub _method_alias { exists $methods{$_[0]||''}; } 
  sub _method_hand_up { exists $hand_ups{$_[0]||''}; } 
}

##############################################################################
# constructors

sub new {
    # Create a new BigFloat object from a string or another bigfloat object.
    # _e: exponent
    # _m: mantissa
    # sign  => sign ("+", "-", "+inf", "-inf", or "NaN"

    my $self    = shift;
    my $selfref = ref $self;
    my $class   = $selfref || $self;

    my ($wanted, @r) = @_;

    # avoid numify-calls by not using || on $wanted!

    unless (defined $wanted) {
        #Carp::carp("Use of uninitialized value in new");
        return $self->bzero(@r);
    }

    # Using $wanted->isa("Math::BigFloat") here causes a 'Deep recursion on
    # subroutine "Math::BigFloat::as_number"' in some tests. Fixme!

    if (UNIVERSAL::isa($wanted, 'Math::BigFloat')) {
        my $copy = $wanted -> copy();
        if ($selfref) {                 # if new() called as instance method
            %$self = %$copy;
        } else {                        # if new() called as class method
            $self = $copy;
        }
        return $copy;
    }

    $class->import() if $IMPORT == 0;             # make require work

    # If called as a class method, initialize a new object.

    $self = bless {}, $class unless $selfref;

    # shortcut for bigints and its subclasses
    if ((ref($wanted)) && $wanted -> can("as_number")) {
        $self->{_m} = $wanted->as_number()->{value};  # get us a bigint copy
        $self->{_e} = $MBI->_zero();
        $self->{_es} = '+';
        $self->{sign} = $wanted->sign();
        return $self->bnorm();
    }

    # else: got a string or something masquerading as number (with overload)

    # Handle Infs.

    if ($wanted =~ /^\s*([+-]?)inf(inity)?\s*\z/i) {
        return $downgrade->new($wanted) if $downgrade;
        my $sgn = $1 || '+';
        $self->{sign} = $sgn . 'inf';   # set a default sign for bstr()
        return $self->binf($sgn);
    }

    # Handle explicit NaNs (not the ones returned due to invalid input).

    if ($wanted =~ /^\s*([+-]?)nan\s*\z/i) {
        return $downgrade->new($wanted) if $downgrade;
        $self = $class -> bnan();
        $self->round(@r) unless @r >= 2 && !defined $r[0] && !defined $r[1];
        return $self;
    }

    # Handle hexadecimal numbers.

    if ($wanted =~ /^\s*[+-]?0[Xx]/) {
        $self = $class -> from_hex($wanted);
        $self->round(@r) unless @r >= 2 && !defined $r[0] && !defined $r[1];
        return $self;
    }

    # Handle binary numbers.

    if ($wanted =~ /^\s*[+-]?0[Bb]/) {
        $self = $class -> from_bin($wanted);
        $self->round(@r) unless @r >= 2 && !defined $r[0] && !defined $r[1];
        return $self;
    }

    # Shortcut for simple forms like '12' that have no trailing zeros.
    if ($wanted =~ /^([+-]?)0*([1-9][0-9]*[1-9])$/) {
        $self->{_e}   = $MBI -> _zero();
        $self->{_es}  = '+';
        $self->{sign} = $1 || '+';
        $self->{_m}   = $MBI -> _new($2);
        if (!$downgrade) {
            $self->round(@r) unless @r >= 2 && !defined $r[0] && !defined $r[1];
            return $self;
        }
    }

  my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
  if (!ref $mis)
    {
    if ($_trap_nan)
      {
      Carp::croak("$wanted is not a number initialized to $class");
      }
    
    return $downgrade->bnan() if $downgrade;
    
    $self->{_e} = $MBI->_zero();
    $self->{_es} = '+';
    $self->{_m} = $MBI->_zero();
    $self->{sign} = $nan;
    }
  else
    {
    # make integer from mantissa by adjusting exp, then convert to int
    $self->{_e} = $MBI->_new($$ev);		# exponent
    $self->{_es} = $$es || '+';
    my $mantissa = "$$miv$$mfv"; 		# create mant.
    $mantissa =~ s/^0+(\d)/$1/;			# strip leading zeros
    $self->{_m} = $MBI->_new($mantissa); 	# create mant.

    # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
    if (CORE::length($$mfv) != 0)
      {
      my $len = $MBI->_new( CORE::length($$mfv));
      ($self->{_e}, $self->{_es}) =
	_e_sub ($self->{_e}, $len, $self->{_es}, '+');
      }
    # we can only have trailing zeros on the mantissa if $$mfv eq ''
    else
      {
      # Use a regexp to count the trailing zeros in $$miv instead of _zeros()
      # because that is faster, especially when _m is not stored in base 10.
      my $zeros = 0; $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/; 
      if ($zeros != 0)
        {
        my $z = $MBI->_new($zeros);
        # turn '120e2' into '12e3'
        $MBI->_rsft ( $self->{_m}, $z, 10);
        ($self->{_e}, $self->{_es}) =
	  _e_add ( $self->{_e}, $z, $self->{_es}, '+');
        }
      }
    $self->{sign} = $$mis;

    # for something like 0Ey, set y to 0, and -0 => +0
    # Check $$miv for being '0' and $$mfv eq '', because otherwise _m could not
    # have become 0. That's faster than to call $MBI->_is_zero().
    $self->{sign} = '+', $self->{_e} = $MBI->_zero()
     if $$miv eq '0' and $$mfv eq '';

    if (!$downgrade) {
        $self->round(@r) unless @r >= 2 && !defined $r[0] && !defined $r[1];
        return $self;
    }
    }
  # if downgrade, inf, NaN or integers go down

  if ($downgrade && $self->{_es} eq '+')
    {
    if ($MBI->_is_zero( $self->{_e} ))
      {
      return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
      }
    return $downgrade->new($self->bsstr()); 
    }
    $self->bnorm();
    $self->round(@r) unless @r >= 2 && !defined $r[0] && !defined $r[1];
    return $self;
}

sub bnan {
    # create/assign a 'NaN'

    if (@_ == 0) {
        Carp::carp("Using bnan() as a function is deprecated;",
                   " use bnan() as a method instead");
        unshift @_, __PACKAGE__;
    }

    my $self    = shift;
    my $selfref = ref $self;
    my $class   = $selfref || $self;

    {
        no strict 'refs';
        if (${"${class}::_trap_nan"}) {
            Carp::croak("Tried to create NaN in $class->bnan()");
        }
    }

    $self->import() if $IMPORT == 0;            # make require work
    return if $self->modify('bnan');

    $self = bless {}, $class unless $selfref;

    $self -> {sign} = $nan;
    $self -> {_m}   = $MBI -> _zero();
    $self -> {_es}  = '+';
    $self -> {_e}   = $MBI -> _zero();

    return $self;
}

sub binf {
    # create/assign a '+inf' or '-inf'

    if (@_ == 0 || (defined($_[0]) && !ref($_[0]) &&
                    $_[0] =~ /^\s*[+-](inf(inity)?)?\s*$/))
    {
        Carp::carp("Using binf() as a function is deprecated;",
                   " use binf() as a method instead");
        unshift @_, __PACKAGE__;
    }

    my $self    = shift;
    my $selfref = ref $self;
    my $class   = $selfref || $self;

    {
        no strict 'refs';
        if (${"${class}::_trap_inf"}) {
            Carp::croak("Tried to create +-inf in $class->binf()");
        }
    }

    $self->import() if $IMPORT == 0;            # make require work
    return if $self->modify('binf');

    my $sign = shift;
    $sign = defined $sign && $sign =~ /^\s*-/ ? "-" : "+";

    $self = bless {}, $class unless $selfref;

    $self -> {sign} = $sign . 'inf';
    $self -> {_m}   = $MBI -> _zero();
    $self -> {_es}  = '+';
    $self -> {_e}   = $MBI -> _zero();

    return $self;
}

sub bzero {
    # create/assign '+0'

    if (@_ == 0) {
        Carp::carp("Using bone() as a function is deprecated;",
                   " use bone() as a method instead");
        unshift @_, __PACKAGE__;
    }

    my $self    = shift;
    my $selfref = ref $self;
    my $class   = $selfref || $self;

    $self->import() if $IMPORT == 0;            # make require work
    return if $self->modify('bzero');

    $self = bless {}, $class unless $selfref;

    $self -> {sign} = '+';
    $self -> {_m}   = $MBI -> _zero();
    $self -> {_es}  = '+';
    $self -> {_e}   = $MBI -> _zero();

    if (@_ > 0) {
        if (@_ > 3) {
            # call like: $x->bzero($a, $p, $r, $y);
            ($self, $self->{_a}, $self->{_p}) = $self->_find_round_parameters(@_);
        } else {
            # call like: $x->bzero($a, $p, $r);
            $self->{_a} = $_[0]
              if !defined $self->{_a} || (defined $_[0] && $_[0] > $self->{_a});
            $self->{_p} = $_[1]
              if !defined $self->{_p} || (defined $_[1] && $_[1] > $self->{_p});
        }
    }

    return $self;
}

sub bone {
    # Create or assign '+1' (or -1 if given sign '-').

    if (@_ == 0 || (defined($_[0]) && ($_[0] eq '+' || $_[0] eq '-'))) {
        Carp::carp("Using bone() as a function is deprecated;",
                   " use bone() as a method instead");
        unshift @_, __PACKAGE__;
    }

    my $self    = shift;
    my $selfref = ref $self;
    my $class   = $selfref || $self;

    $self->import() if $IMPORT == 0;            # make require work
    return if $self->modify('bone');

    my $sign = shift;
    $sign = defined $sign && $sign =~ /^\s*-/ ? "-" : "+";

    $self = bless {}, $class unless $selfref;

    $self -> {sign} = $sign;
    $self -> {_m}   = $MBI -> _one();
    $self -> {_es}  = '+';
    $self -> {_e}   = $MBI -> _zero();

    if (@_ > 0) {
        if (@_ > 3) {
            # call like: $x->bone($sign, $a, $p, $r, $y, ...);
            ($self, $self->{_a}, $self->{_p}) = $self->_find_round_parameters(@_);
        } else {
            # call like: $x->bone($sign, $a, $p, $r);
            $self->{_a} = $_[0]
              if ( (!defined $self->{_a}) || (defined $_[0] && $_[0] > $self->{_a}));
            $self->{_p} = $_[1]
              if ( (!defined $self->{_p}) || (defined $_[1] && $_[1] > $self->{_p}));
        }
    }

    return $self;
}

sub copy {
    my $self    = shift;
    my $selfref = ref $self;
    my $class   = $selfref || $self;

    # If called as a class method, the object to copy is the next argument.

    $self = shift() unless $selfref;

    my $copy = bless {}, $class;

    $copy->{sign} = $self->{sign};
    $copy->{_es}  = $self->{_es};
    $copy->{_m}   = $MBI->_copy($self->{_m});
    $copy->{_e}   = $MBI->_copy($self->{_e});
    $copy->{_a}   = $self->{_a} if exists $self->{_a};
    $copy->{_p}   = $self->{_p} if exists $self->{_p};

    return $copy;
}

sub isa
  {
  my ($self,$class) = @_;
  return if $class =~ /^Math::BigInt/;		# we aren't one of these
  UNIVERSAL::isa($self,$class);
  }

sub config
  {
  # return (later set?) configuration data as hash ref
  my $class = shift || 'Math::BigFloat';

  if (@_ == 1 && ref($_[0]) ne 'HASH')
    {
    my $cfg = $class->SUPER::config();
    return $cfg->{$_[0]};
    }

  my $cfg = $class->SUPER::config(@_);

  # now we need only to override the ones that are different from our parent
  $cfg->{class} = $class;
  $cfg->{with} = $MBI;
  $cfg;
  }

##############################################################################
# string conversion

sub bstr 
  {
  # (ref to BFLOAT or num_str ) return num_str
  # Convert number from internal format to (non-scientific) string format.
  # internal format is always normalized (no leading zeros, "-0" => "+0")
  my ($class,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);

  if ($x->{sign} !~ /^[+-]$/)
    {
    return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
    return 'inf';                                       # +inf
    }

  my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';

  # $x is zero?
  my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
  if ($not_zero)
    {
    $es = $MBI->_str($x->{_m});
    $len = CORE::length($es);
    my $e = $MBI->_num($x->{_e});	
    $e = -$e if $x->{_es} eq '-';
    if ($e < 0)
      {
      $dot = '';
      # if _e is bigger than a scalar, the following will blow your memory
      if ($e <= -$len)
        {
        my $r = abs($e) - $len;
        $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
        }
      else
        {
        substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
        $cad = -$cad if $x->{_es} eq '-';
        }
      }
    elsif ($e > 0)
      {
      # expand with zeros
      $es .= '0' x $e; $len += $e; $cad = 0;
      }
    } # if not zero

  $es = '-'.$es if $x->{sign} eq '-';
  # if set accuracy or precision, pad with zeros on the right side
  if ((defined $x->{_a}) && ($not_zero))
    {
    # 123400 => 6, 0.1234 => 4, 0.001234 => 4
    my $zeros = $x->{_a} - $cad;		# cad == 0 => 12340
    $zeros = $x->{_a} - $len if $cad != $len;
    $es .= $dot.'0' x $zeros if $zeros > 0;
    }
  elsif ((($x->{_p} || 0) < 0))
    {
    # 123400 => 6, 0.1234 => 4, 0.001234 => 6
    my $zeros = -$x->{_p} + $cad;
    $es .= $dot.'0' x $zeros if $zeros > 0;
    }
  $es;
  }

sub bsstr
  {
  # (ref to BFLOAT or num_str ) return num_str
  # Convert number from internal format to scientific string format.
  # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
  my ($class,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);

  if ($x->{sign} !~ /^[+-]$/)
    {
    return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
    return 'inf';                                       # +inf
    }
  my $sep = 'e'.$x->{_es};
  my $sign = $x->{sign}; $sign = '' if $sign eq '+';
  $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
  }
    
sub numify 
  {
  # Make a Perl scalar number from a Math::BigFloat object.
  my ($class,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);

  if ($x -> is_nan()) {
      require Math::Complex;
      my $inf = Math::Complex::Inf();
      return $inf - $inf;
  }

  if ($x -> is_inf()) {
      require Math::Complex;
      my $inf = Math::Complex::Inf();
      return $x -> is_negative() ? -$inf : $inf;
  }

  # Create a string and let Perl's atoi()/atof() handle the rest.
  return 0 + $x -> bsstr();
  }

##############################################################################
# public stuff (usually prefixed with "b")

sub bneg
  {
  # (BINT or num_str) return BINT
  # negate number or make a negated number from string
  my ($class,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);

  return $x if $x->modify('bneg');

  # for +0 do not negate (to have always normalized +0). Does nothing for 'NaN'
  $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
  $x;
  }

sub bcmp 
  {
  # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)

  # set up parameters
  my ($class,$x,$y) = (ref($_[0]),@_);

  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($class,$x,$y) = objectify(2,@_);
    }

  return $upgrade->bcmp($x,$y) if defined $upgrade &&
    ((!$x->isa($class)) || (!$y->isa($class)));

  # Handle all 'nan' cases.

  return undef if ($x->{sign} eq $nan) || ($y->{sign} eq $nan);

  # Handle all '+inf' and '-inf' cases.

  return  0 if ($x->{sign} eq '+inf' && $y->{sign} eq '+inf' ||
                $x->{sign} eq '-inf' && $y->{sign} eq '-inf');
  return +1 if $x->{sign} eq '+inf';    # x = +inf and y < +inf
  return -1 if $x->{sign} eq '-inf';    # x = -inf and y > -inf
  return -1 if $y->{sign} eq '+inf';    # x < +inf and y = +inf
  return +1 if $y->{sign} eq '-inf';    # x > -inf and y = -inf

  # Handle all cases with opposite signs.

  return +1 if $x->{sign} eq '+' && $y->{sign} eq '-';  # also does 0 <=> -y
  return -1 if $x->{sign} eq '-' && $y->{sign} eq '+';  # also does -x <=> 0

  # Handle all remaining zero cases.

  my $xz = $x->is_zero();
  my $yz = $y->is_zero();
  return  0 if $xz && $yz;                              # 0 <=> 0
  return -1 if $xz && $y->{sign} eq '+';                # 0 <=> +y
  return +1 if $yz && $x->{sign} eq '+';                # +x <=> 0

  # Both arguments are now finite, non-zero numbers with the same sign.

  my $cmp;

  # The next step is to compare the exponents, but since each mantissa is an
  # integer of arbitrary value, the exponents must be normalized by the length
  # of the mantissas before we can compare them.

  my $mxl = $MBI->_len($x->{_m});
  my $myl = $MBI->_len($y->{_m});

  # If the mantissas have the same length, there is no point in normalizing the
  # exponents by the length of the mantissas, so treat that as a special case.

  if ($mxl == $myl) {

      # First handle the two cases where the exponents have different signs.

      if ($x->{_es} eq '+' && $y->{_es} eq '-') {
          $cmp = +1;
      }

      elsif ($x->{_es} eq '-' && $y->{_es} eq '+') {
          $cmp = -1;
      }

      # Then handle the case where the exponents have the same sign.

      else {
          $cmp = $MBI->_acmp($x->{_e}, $y->{_e});
          $cmp = -$cmp if $x->{_es} eq '-';
      }

      # Adjust for the sign, which is the same for x and y, and bail out if
      # we're done.

      $cmp = -$cmp if $x->{sign} eq '-';        # 124 > 123, but -124 < -123
      return $cmp if $cmp;

  }

  # We must normalize each exponent by the length of the corresponding
  # mantissa. Life is a lot easier if we first make both exponents
  # non-negative. We do this by adding the same positive value to both
  # exponent. This is safe, because when comparing the exponents, only the
  # relative difference is important.

  my $ex;
  my $ey;

  if ($x->{_es} eq '+') {

      # If the exponent of x is >= 0 and the exponent of y is >= 0, there is no
      # need to do anything special.

      if ($y->{_es} eq '+') {
          $ex = $MBI->_copy($x->{_e});
          $ey = $MBI->_copy($y->{_e});
      }

      # If the exponent of x is >= 0 and the exponent of y is < 0, add the
      # absolute value of the exponent of y to both.

      else {
          $ex = $MBI->_copy($x->{_e});
          $ex = $MBI->_add($ex, $y->{_e});      # ex + |ey|
          $ey = $MBI->_zero();                  # -ex + |ey| = 0
      }

  } else {

      # If the exponent of x is < 0 and the exponent of y is >= 0, add the
      # absolute value of the exponent of x to both.

      if ($y->{_es} eq '+') {
          $ex = $MBI->_zero();                  # -ex + |ex| = 0
          $ey = $MBI->_copy($y->{_e});
          $ey = $MBI->_add($ey, $x->{_e});      # ey + |ex|
      }

      # If the exponent of x is < 0 and the exponent of y is < 0, add the
      # absolute values of both exponents to both exponents.

      else {
          $ex = $MBI->_copy($y->{_e});          # -ex + |ey| + |ex| = |ey|
          $ey = $MBI->_copy($x->{_e});          # -ey + |ex| + |ey| = |ex|
      }

  }

  # Now we can normalize the exponents by adding lengths of the mantissas.

  $MBI->_add($ex, $MBI->_new($mxl));
  $MBI->_add($ey, $MBI->_new($myl));

  # We're done if the exponents are different.

  $cmp = $MBI->_acmp($ex, $ey);
  $cmp = -$cmp if $x->{sign} eq '-';            # 124 > 123, but -124 < -123
  return $cmp if $cmp;

  # Compare the mantissas, but first normalize them by padding the shorter
  # mantissa with zeros (shift left) until it has the same length as the longer
  # mantissa.

  my $mx = $x->{_m};
  my $my = $y->{_m};

  if ($mxl > $myl) {
      $my = $MBI->_lsft($MBI->_copy($my), $MBI->_new($mxl - $myl), 10);
  } elsif ($mxl < $myl) {
      $mx = $MBI->_lsft($MBI->_copy($mx), $MBI->_new($myl - $mxl), 10);
  }

  $cmp = $MBI->_acmp($mx, $my);
  $cmp = -$cmp if $x->{sign} eq '-';            # 124 > 123, but -124 < -123
  return $cmp;

  }

sub bacmp 
  {
  # Compares 2 values, ignoring their signs. 
  # Returns one of undef, <0, =0, >0. (suitable for sort)
  
  # set up parameters
  my ($class,$x,$y) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($class,$x,$y) = objectify(2,@_);
    }

  return $upgrade->bacmp($x,$y) if defined $upgrade &&
    ((!$x->isa($class)) || (!$y->isa($class)));

  # handle +-inf and NaN's
  if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
    {
    return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
    return 0 if ($x->is_inf() && $y->is_inf());
    return 1 if ($x->is_inf() && !$y->is_inf());
    return -1;
    }

  # shortcut 
  my $xz = $x->is_zero();
  my $yz = $y->is_zero();
  return 0 if $xz && $yz;				# 0 <=> 0
  return -1 if $xz && !$yz;				# 0 <=> +y
  return 1 if $yz && !$xz;				# +x <=> 0

  # adjust so that exponents are equal
  my $lxm = $MBI->_len($x->{_m});
  my $lym = $MBI->_len($y->{_m});
  my ($xes,$yes) = (1,1);
  $xes = -1 if $x->{_es} ne '+';
  $yes = -1 if $y->{_es} ne '+';
  # the numify somewhat limits our length, but makes it much faster
  my $lx = $lxm + $xes * $MBI->_num($x->{_e});
  my $ly = $lym + $yes * $MBI->_num($y->{_e});
  my $l = $lx - $ly;
  return $l <=> 0 if $l != 0;
  
  # lengths (corrected by exponent) are equal
  # so make mantissa equal-length by padding with zero (shift left)
  my $diff = $lxm - $lym;
  my $xm = $x->{_m};		# not yet copy it
  my $ym = $y->{_m};
  if ($diff > 0)
    {
    $ym = $MBI->_copy($y->{_m});
    $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
    }
  elsif ($diff < 0)
    {
    $xm = $MBI->_copy($x->{_m});
    $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
    }
  $MBI->_acmp($xm,$ym);
  }

sub badd 
  {
  # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
  # return result as BFLOAT

  # set up parameters
  my ($class,$x,$y,@r) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($class,$x,$y,@r) = objectify(2,@_);
    }
 
  return $x if $x->modify('badd');

  # inf and NaN handling
  if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
    {
    # NaN first
    return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
    # inf handling
    if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
      {
      # +inf++inf or -inf+-inf => same, rest is NaN
      return $x if $x->{sign} eq $y->{sign};
      return $x->bnan();
      }
    # +-inf + something => +inf; something +-inf => +-inf
    $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
    return $x;
    }

  return $upgrade->badd($x,$y,@r) if defined $upgrade &&
   ((!$x->isa($class)) || (!$y->isa($class)));

  $r[3] = $y;						# no push!

  # speed: no add for 0+y or x+0
  return $x->bround(@r) if $y->is_zero();		# x+0
  if ($x->is_zero())					# 0+y
    {
    # make copy, clobbering up x (modify in place!)
    $x->{_e} = $MBI->_copy($y->{_e});
    $x->{_es} = $y->{_es};
    $x->{_m} = $MBI->_copy($y->{_m});
    $x->{sign} = $y->{sign} || $nan;
    return $x->round(@r);
    }
 
  # take lower of the two e's and adapt m1 to it to match m2
  my $e = $y->{_e};
  $e = $MBI->_zero() if !defined $e;		# if no BFLOAT?
  $e = $MBI->_copy($e);				# make copy (didn't do it yet)

  my $es;

  ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});

  my $add = $MBI->_copy($y->{_m});

  if ($es eq '-')				# < 0
    {
    $MBI->_lsft( $x->{_m}, $e, 10);
    ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
    }
  elsif (!$MBI->_is_zero($e))			# > 0
    {
    $MBI->_lsft($add, $e, 10);
    }
  # else: both e are the same, so just leave them

  if ($x->{sign} eq $y->{sign})
    {
    # add
    $x->{_m} = $MBI->_add($x->{_m}, $add);
    }
  else
    {
    ($x->{_m}, $x->{sign}) = 
     _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
    }

  # delete trailing zeros, then round
  $x->bnorm()->round(@r);
  }

# sub bsub is inherited from Math::BigInt!

sub bsub {
    # (BINT or num_str, BINT or num_str) return BINT
    # subtract second arg from first, modify first

    # set up parameters
    my ($class, $x, $y, @r) = (ref($_[0]), @_);

    # objectify is costly, so avoid it
    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) {
        ($class, $x, $y, @r) = objectify(2, @_);
    }

    return $x if $x -> modify('bsub');

    return $upgrade -> new($x) -> bsub($upgrade -> new($y), @r)
      if defined $upgrade && ((!$x -> isa($class)) || (!$y -> isa($class)));

    return $x -> round(@r) if $y -> is_zero();

    # To correctly handle the lone special case $x -> bsub($x), we note the
    # sign of $x, then flip the sign from $y, and if the sign of $x did change,
    # too, then we caught the special case:

    my $xsign = $x -> {sign};
    $y -> {sign} =~ tr/+-/-+/;  # does nothing for NaN
    if ($xsign ne $x -> {sign}) {
        # special case of $x -> bsub($x) results in 0
        return $x -> bzero(@r) if $xsign =~ /^[+-]$/;
        return $x -> bnan();    # NaN, -inf, +inf
    }
    $x -> badd($y, @r);         # badd does not leave internal zeros
    $y -> {sign} =~ tr/+-/-+/;  # refix $y (does nothing for NaN)
    $x;                         # already rounded by badd() or no rounding
}

sub binc
  {
  # increment arg by one
  my ($class,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);

  return $x if $x->modify('binc');

  if ($x->{_es} eq '-')
    {
    return $x->badd($class->bone(),@r);	#  digits after dot
    }

  if (!$MBI->_is_zero($x->{_e}))		# _e == 0 for NaN, inf, -inf
    {
    # 1e2 => 100, so after the shift below _m has a '0' as last digit
    $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10);	# 1e2 => 100
    $x->{_e} = $MBI->_zero();				# normalize
    $x->{_es} = '+';
    # we know that the last digit of $x will be '1' or '9', depending on the
    # sign
    }
  # now $x->{_e} == 0
  if ($x->{sign} eq '+')
    {
    $MBI->_inc($x->{_m});
    return $x->bnorm()->bround(@r);
    }
  elsif ($x->{sign} eq '-')
    {
    $MBI->_dec($x->{_m});
    $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
    return $x->bnorm()->bround(@r);
    }
  # inf, nan handling etc
  $x->badd($class->bone(),@r);			# badd() does round 
  }

sub bdec
  {
  # decrement arg by one
  my ($class,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);

  return $x if $x->modify('bdec');

  if ($x->{_es} eq '-')
    {
    return $x->badd($class->bone('-'),@r);	#  digits after dot
    }

  if (!$MBI->_is_zero($x->{_e}))
    {
    $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10);	# 1e2 => 100
    $x->{_e} = $MBI->_zero();				# normalize
    $x->{_es} = '+';
    }
  # now $x->{_e} == 0
  my $zero = $x->is_zero();
  # <= 0
  if (($x->{sign} eq '-') || $zero)
    {
    $MBI->_inc($x->{_m});
    $x->{sign} = '-' if $zero;				# 0 => 1 => -1
    $x->{sign} = '+' if $MBI->_is_zero($x->{_m});	# -1 +1 => -0 => +0
    return $x->bnorm()->round(@r);
    }
  # > 0
  elsif ($x->{sign} eq '+')
    {
    $MBI->_dec($x->{_m});
    return $x->bnorm()->round(@r);
    }
  # inf, nan handling etc
  $x->badd($class->bone('-'),@r);		# does round
  } 

sub DEBUG () { 0; }

sub blog
  {
  my ($class,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);

  # If called as $x -> blog() or $x -> blog(undef), don't objectify the
  # undefined base, since undef signals that the base is Euler's number.
  #unless (ref($x) && !defined($base)) {
  #    # objectify is costly, so avoid it
  #    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) {
  #        ($class,$x,$base,$a,$p,$r) = objectify(2,@_);
  #    }
  #}

  return $x if $x->modify('blog');

  return $x -> bnan() if $x -> is_nan();

  # we need to limit the accuracy to protect against overflow
  my $fallback = 0;
  my ($scale,@params);
  ($x,@params) = $x->_find_round_parameters($a,$p,$r);

  # no rounding at all, so must use fallback
  if (scalar @params == 0)
    {
    # simulate old behaviour
    $params[0] = $class->div_scale();	# and round to it as accuracy
    $params[1] = undef;			# P = undef
    $scale = $params[0]+4; 		# at least four more for proper round
    $params[2] = $r;			# round mode by caller or undef
    $fallback = 1;			# to clear a/p afterwards
    }
  else
    {
    # the 4 below is empirical, and there might be cases where it is not
    # enough...
    $scale = abs($params[0] || $params[1]) + 4;	# take whatever is defined
    }

  my $done = 0;
  if (defined $base) {
      $base = $class -> new($base) unless ref $base;
      if ($base -> is_nan() || $base -> is_one()) {
          $x -> bnan();
          $done = 1;
      } elsif ($base -> is_inf() || $base -> is_zero()) {
          if ($x -> is_inf() || $x -> is_zero()) {
              $x -> bnan();
          } else {
              $x -> bzero(@params);
          }
          $done = 1;
      } elsif ($base -> is_negative()) {        # -inf < base < 0
          if ($x -> is_one()) {                 #     x = 1
              $x -> bzero(@params);
          } elsif ($x == $base) {
              $x -> bone('+', @params);         #     x = base
          } else {
              $x -> bnan();                     #     otherwise
          }
          $done = 1;
      } elsif ($x == $base) {
          $x -> bone('+', @params);             # 0 < base && 0 < x < inf
          $done = 1;
      }
  }

  # We now know that the base is either undefined or positive and finite.

  unless ($done) {
      if ($x -> is_inf()) {             #   x = +/-inf
          my $sign = defined $base && $base < 1 ? '-' : '+';
          $x -> binf($sign);
          $done = 1;
      } elsif ($x -> is_neg()) {        #   -inf < x < 0
          $x -> bnan();
          $done = 1;
      } elsif ($x -> is_one()) {        #   x = 1
          $x -> bzero(@params);
          $done = 1;
      } elsif ($x -> is_zero()) {       #   x = 0
          my $sign = defined $base && $base < 1 ? '+' : '-';
          $x -> binf($sign);
          $done = 1;
      }
  }

  if ($done) {
      if ($fallback) {
        # clear a/p after round, since user did not request it
          delete $x->{_a};
          delete $x->{_p};
        }
      return $x;
      }

  # when user set globals, they would interfere with our calculation, so
  # disable them and later re-enable them
  no strict 'refs';
  my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef;
  my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef;
  # we also need to disable any set A or P on $x (_find_round_parameters took
  # them already into account), since these would interfere, too
  delete $x->{_a}; delete $x->{_p};
  # need to disable $upgrade in BigInt, to avoid deep recursion
  local $Math::BigInt::upgrade = undef;
  local $Math::BigFloat::downgrade = undef;

  # upgrade $x if $x is not a BigFloat (handle BigInt input)
  # XXX TODO: rebless!
  if (!$x->isa('Math::BigFloat'))
    {
    $x = Math::BigFloat->new($x);
    $class = ref($x);
    }
  
  $done = 0;

  # If the base is defined and an integer, try to calculate integer result
  # first. This is very fast, and in case the real result was found, we can
  # stop right here.
  if (defined $base && $base->is_int() && $x->is_int())
    {
    my $i = $MBI->_copy( $x->{_m} );
    $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
    my $int = Math::BigInt->bzero();
    $int->{value} = $i;
    $int->blog($base->as_number());
    # if ($exact)
    if ($base->as_number()->bpow($int) == $x)
      {
      # found result, return it
      $x->{_m} = $int->{value};
      $x->{_e} = $MBI->_zero();
      $x->{_es} = '+';
      $x->bnorm();
      $done = 1;
      }
    }

  if ($done == 0)
    {
    # base is undef, so base should be e (Euler's number), so first calculate the
    # log to base e (using reduction by 10 (and probably 2)):
    $class->_log_10($x,$scale);

    # and if a different base was requested, convert it
    if (defined $base)
      {
      $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
      # not ln, but some other base (don't modify $base)
      $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
      }
    }
 
  # shortcut to not run through _find_round_parameters again
  if (defined $params[0])
    {
    $x->bround($params[0],$params[2]);		# then round accordingly
    }
  else
    {
    $x->bfround($params[1],$params[2]);		# then round accordingly
    }
  if ($fallback)
    {
    # clear a/p after round, since user did not request it
    delete $x->{_a}; delete $x->{_p};
    }
  # restore globals
  $$abr = $ab; $$pbr = $pb;

  $x;
  }

sub _len_to_steps
  {
  # Given D (digits in decimal), compute N so that N! (N factorial) is
  # at least D digits long. D should be at least 50.
  my $d = shift;

  # two constants for the Ramanujan estimate of ln(N!)
  my $lg2 = log(2 * 3.14159265) / 2;
  my $lg10 = log(10);

  # D = 50 => N => 42, so L = 40 and R = 50
  my $l = 40; my $r = $d;

  # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
  $l = $l->numify if ref($l);
  $r = $r->numify if ref($r);
  $lg2 = $lg2->numify if ref($lg2);
  $lg10 = $lg10->numify if ref($lg10);

  # binary search for the right value (could this be written as the reverse of lg(n!)?)
  while ($r - $l > 1)
    {
    my $n = int(($r - $l) / 2) + $l;
    my $ramanujan = 
      int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
    $ramanujan > $d ? $r = $n : $l = $n;
    }
  $l;
  }

sub bnok
  {
  # Calculate n over k (binomial coefficient or "choose" function) as integer.
  # set up parameters
  my ($class,$x,$y,@r) = (ref($_[0]),@_);

  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($class,$x,$y,@r) = objectify(2,@_);
    }

  return $x if $x->modify('bnok');

  return $x->bnan() if $x->is_nan() || $y->is_nan();
  return $x->binf() if $x->is_inf();

  my $u = $x->as_int();
  $u->bnok($y->as_int());

  $x->{_m} = $u->{value};
  $x->{_e} = $MBI->_zero();
  $x->{_es} = '+';
  $x->{sign} = '+';
  $x->bnorm(@r);
  }

sub bexp
  {
  # Calculate e ** X (Euler's number to the power of X)
  my ($class,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);

  return $x if $x->modify('bexp');

  return $x->binf() if $x->{sign} eq '+inf';
  return $x->bzero() if $x->{sign} eq '-inf';

  # we need to limit the accuracy to protect against overflow
  my $fallback = 0;
  my ($scale,@params);
  ($x,@params) = $x->_find_round_parameters($a,$p,$r);

  # also takes care of the "error in _find_round_parameters?" case
  return $x if $x->{sign} eq 'NaN';

  # no rounding at all, so must use fallback
  if (scalar @params == 0)
    {
    # simulate old behaviour
    $params[0] = $class->div_scale();	# and round to it as accuracy
    $params[1] = undef;			# P = undef
    $scale = $params[0]+4; 		# at least four more for proper round
    $params[2] = $r;			# round mode by caller or undef
    $fallback = 1;			# to clear a/p afterwards
    }
  else
    {
    # the 4 below is empirical, and there might be cases where it's not enough...
    $scale = abs($params[0] || $params[1]) + 4;	# take whatever is defined
    }

  return $x->bone(@params) if $x->is_zero();

  if (!$x->isa('Math::BigFloat'))
    {
    $x = Math::BigFloat->new($x);
    $class = ref($x);
    }
  
  # when user set globals, they would interfere with our calculation, so
  # disable them and later re-enable them
  no strict 'refs';
  my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef;
  my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef;
  # we also need to disable any set A or P on $x (_find_round_parameters took
  # them already into account), since these would interfere, too
  delete $x->{_a}; delete $x->{_p};
  # need to disable $upgrade in BigInt, to avoid deep recursion
  local $Math::BigInt::upgrade = undef;
  local $Math::BigFloat::downgrade = undef;

  my $x_org = $x->copy();

  # We use the following Taylor series:

  #           x    x^2   x^3   x^4
  #  e = 1 + --- + --- + --- + --- ...
  #           1!    2!    3!    4!

  # The difference for each term is X and N, which would result in:
  # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term

  # But it is faster to compute exp(1) and then raising it to the
  # given power, esp. if $x is really big and an integer because:

  #  * The numerator is always 1, making the computation faster
  #  * the series converges faster in the case of x == 1
  #  * We can also easily check when we have reached our limit: when the
  #    term to be added is smaller than "1E$scale", we can stop - f.i.
  #    scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
  #  * we can compute the *exact* result by simulating bigrat math:

  #  1   1    gcd(3,4) = 1    1*24 + 1*6    5
  #  - + -                  = ---------- =  --                 
  #  6   24                      6*24       24

  # We do not compute the gcd() here, but simple do:
  #  1   1    1*24 + 1*6   30
  #  - + -  = --------- =  --                 
  #  6   24       6*24     144

  # In general:
  #  a   c    a*d + c*b 	and note that c is always 1 and d = (b*f)
  #  - + -  = ---------
  #  b   d       b*d

  # This leads to:         which can be reduced by b to:
  #  a   1     a*b*f + b    a*f + 1
  #  - + -   = --------- =  -------
  #  b   b*f     b*b*f        b*f

  # The first terms in the series are:

  # 1     1    1    1    1    1     1     1     13700
  # -- + -- + -- + -- + -- + --- + --- + ---- = -----
  # 1     1    2    6   24   120   720   5040   5040

  # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!

  if ($scale <= 75)
    {
    # set $x directly from a cached string form
    $x->{_m} = $MBI->_new(
    "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
    $x->{sign} = '+';
    $x->{_es} = '-';
    $x->{_e} = $MBI->_new(79);
    }
  else
    {
    # compute A and B so that e = A / B.
 
    # After some terms we end up with this, so we use it as a starting point:
    my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
    my $F = $MBI->_new(42); my $step = 42;

    # Compute how many steps we need to take to get $A and $B sufficiently big
    my $steps = _len_to_steps($scale - 4);
#    print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
    while ($step++ <= $steps)
      {
      # calculate $a * $f + 1
      $A = $MBI->_mul($A, $F);
      $A = $MBI->_inc($A);
      # increment f
      $F = $MBI->_inc($F);
      }
    # compute $B as factorial of $steps (this is faster than doing it manually)
    my $B = $MBI->_fac($MBI->_new($steps));
    
#  print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";

    # compute A/B with $scale digits in the result (truncate, not round)
    $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
    $A = $MBI->_div( $A, $B );

    $x->{_m} = $A;
    $x->{sign} = '+';
    $x->{_es} = '-';
    $x->{_e} = $MBI->_new($scale);
    }

  # $x contains now an estimate of e, with some surplus digits, so we can round
  if (!$x_org->is_one())
    {
    # Reduce size of fractional part, followup with integer power of two.
    my $lshift = 0;
    while ($lshift < 30 && $x_org->bacmp(2 << $lshift) > 0)
      {
        $lshift++;
      }
    # Raise $x to the wanted power and round it.
    if ($lshift == 0)
      {
        $x->bpow($x_org, @params);
      }
    else
      {
        my($mul, $rescale) = (1 << $lshift, $scale+1+$lshift);
        $x->bpow(scalar $x_org->bdiv($mul,$rescale),$rescale)->bpow($mul, @params);
      }
    }
  else
    {
    # else just round the already computed result
    delete $x->{_a}; delete $x->{_p};
    # shortcut to not run through _find_round_parameters again
    if (defined $params[0])
      {
      $x->bround($params[0],$params[2]);		# then round accordingly
      }
    else
      {
      $x->bfround($params[1],$params[2]);		# then round accordingly
      }
    }
  if ($fallback)
    {
    # clear a/p after round, since user did not request it
    delete $x->{_a}; delete $x->{_p};
    }
  # restore globals
  $$abr = $ab; $$pbr = $pb;

  $x;						# return modified $x
  }

sub _log
  {
  # internal log function to calculate ln() based on Taylor series.
  # Modifies $x in place.
  my ($class,$x,$scale) = @_;

  # in case of $x == 1, result is 0
  return $x->bzero() if $x->is_one();

  # XXX TODO: rewrite this in a similar manner to bexp()

  # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log

  # u = x-1, v = x+1
  #              _                               _
  # Taylor:     |    u    1   u^3   1   u^5       |
  # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 0
  #             |_   v    3   v^3   5   v^5      _|

  # This takes much more steps to calculate the result and is thus not used
  # u = x-1
  #              _                               _
  # Taylor:     |    u    1   u^2   1   u^3       |
  # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 1/2
  #             |_   x    2   x^2   3   x^3      _|

  my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);

  $v = $x->copy(); $v->binc();		# v = x+1
  $x->bdec(); $u = $x->copy();		# u = x-1; x = x-1
  $x->bdiv($v,$scale);			# first term: u/v
  $below = $v->copy();
  $over = $u->copy();
  $u *= $u; $v *= $v;				# u^2, v^2
  $below->bmul($v);				# u^3, v^3
  $over->bmul($u);
  $factor = $class->new(3); $f = $class->new(2);

  my $steps = 0;
  $limit = $class->new("1E-". ($scale-1));
  while (3 < 5)
    {
    # we calculate the next term, and add it to the last
    # when the next term is below our limit, it won't affect the outcome
    # anymore, so we stop

    # calculating the next term simple from over/below will result in quite
    # a time hog if the input has many digits, since over and below will
    # accumulate more and more digits, and the result will also have many
    # digits, but in the end it is rounded to $scale digits anyway. So if we
    # round $over and $below first, we save a lot of time for the division
    # (not with log(1.2345), but try log (123**123) to see what I mean. This
    # can introduce a rounding error if the division result would be f.i.
    # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
    # if we truncated $over and $below we might get 0.12345. Does this matter
    # for the end result? So we give $over and $below 4 more digits to be
    # on the safe side (unscientific error handling as usual... :+D

    $next = $over->copy->bround($scale+4)->bdiv(
      $below->copy->bmul($factor)->bround($scale+4), 
      $scale);

## old version:    
##    $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);

    last if $next->bacmp($limit) <= 0;

    delete $next->{_a}; delete $next->{_p};
    $x->badd($next);
    # calculate things for the next term
    $over *= $u; $below *= $v; $factor->badd($f);
    if (DEBUG)
      {
      $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
      }
    }
  print "took $steps steps\n" if DEBUG;
  $x->bmul($f);					# $x *= 2
  }

sub _log_10
  {
  # Internal log function based on reducing input to the range of 0.1 .. 9.99
  # and then "correcting" the result to the proper one. Modifies $x in place.
  my ($class,$x,$scale) = @_;

  # Taking blog() from numbers greater than 10 takes a *very long* time, so we
  # break the computation down into parts based on the observation that:
  #  blog(X*Y) = blog(X) + blog(Y)
  # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
  # $x is the faster it gets. Since 2*$x takes about 10 times as
  # long, we make it faster by about a factor of 100 by dividing $x by 10.

  # The same observation is valid for numbers smaller than 0.1, e.g. computing
  # log(1) is fastest, and the further away we get from 1, the longer it takes.
  # So we also 'break' this down by multiplying $x with 10 and subtract the
  # log(10) afterwards to get the correct result.

  # To get $x even closer to 1, we also divide by 2 and then use log(2) to
  # correct for this. For instance if $x is 2.4, we use the formula:
  #  blog(2.4 * 2) == blog (1.2) + blog(2)
  # and thus calculate only blog(1.2) and blog(2), which is faster in total
  # than calculating blog(2.4).

  # In addition, the values for blog(2) and blog(10) are cached.

  # Calculate nr of digits before dot:
  my $dbd = $MBI->_num($x->{_e});
  $dbd = -$dbd if $x->{_es} eq '-';
  $dbd += $MBI->_len($x->{_m});

  # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
  # infinite recursion

  my $calc = 1;					# do some calculation?

  # disable the shortcut for 10, since we need log(10) and this would recurse
  # infinitely deep
  if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
    {
    $dbd = 0;					# disable shortcut
    # we can use the cached value in these cases
    if ($scale <= $LOG_10_A)
      {
      $x->bzero(); $x->badd($LOG_10);		# modify $x in place
      $calc = 0; 				# no need to calc, but round
      }
    # if we can't use the shortcut, we continue normally
    }
  else
    {
    # disable the shortcut for 2, since we maybe have it cached
    if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
      {
      $dbd = 0;					# disable shortcut
      # we can use the cached value in these cases
      if ($scale <= $LOG_2_A)
        {
        $x->bzero(); $x->badd($LOG_2);		# modify $x in place
        $calc = 0; 				# no need to calc, but round
        }
      # if we can't use the shortcut, we continue normally
      }
    }

  # if $x = 0.1, we know the result must be 0-log(10)
  if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
      $MBI->_is_one($x->{_m}))
    {
    $dbd = 0;					# disable shortcut
    # we can use the cached value in these cases
    if ($scale <= $LOG_10_A)
      {
      $x->bzero(); $x->bsub($LOG_10);
      $calc = 0; 				# no need to calc, but round
      }
    }

  return if $calc == 0;				# already have the result

  # default: these correction factors are undef and thus not used
  my $l_10;				# value of ln(10) to A of $scale
  my $l_2;				# value of ln(2) to A of $scale

  my $two = $class->new(2);

  # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
  # so don't do this shortcut for 1 or 0
  if (($dbd > 1) || ($dbd < 0))
    {
    # convert our cached value to an object if not already (avoid doing this
    # at import() time, since not everybody needs this)
    $LOG_10 = $class->new($LOG_10,undef,undef) unless ref $LOG_10;

    #print "x = $x, dbd = $dbd, calc = $calc\n";
    # got more than one digit before the dot, or more than one zero after the
    # dot, so do:
    #  log(123)    == log(1.23) + log(10) * 2
    #  log(0.0123) == log(1.23) - log(10) * 2
  
    if ($scale <= $LOG_10_A)
      {
      # use cached value
      $l_10 = $LOG_10->copy();		# copy for mul
      }
    else
      {
      # else: slower, compute and cache result
      # also disable downgrade for this code path
      local $Math::BigFloat::downgrade = undef;

      # shorten the time to calculate log(10) based on the following:
      # log(1.25 * 8) = log(1.25) + log(8)
      #               = log(1.25) + log(2) + log(2) + log(2)

      # first get $l_2 (and possible compute and cache log(2))
      $LOG_2 = $class->new($LOG_2,undef,undef) unless ref $LOG_2;
      if ($scale <= $LOG_2_A)
        {
        # use cached value
        $l_2 = $LOG_2->copy();			# copy() for the mul below
        }
      else
        {
        # else: slower, compute and cache result
        $l_2 = $two->copy(); $class->_log($l_2, $scale); # scale+4, actually
        $LOG_2 = $l_2->copy();			# cache the result for later
						# the copy() is for mul below
        $LOG_2_A = $scale;
        }

      # now calculate log(1.25):
      $l_10 = $class->new('1.25'); $class->_log($l_10, $scale); # scale+4, actually

      # log(1.25) + log(2) + log(2) + log(2):
      $l_10->badd($l_2);
      $l_10->badd($l_2);
      $l_10->badd($l_2);
      $LOG_10 = $l_10->copy();		# cache the result for later
					# the copy() is for mul below
      $LOG_10_A = $scale;
      }
    $dbd-- if ($dbd > 1); 		# 20 => dbd=2, so make it dbd=1	
    $l_10->bmul( $class->new($dbd));	# log(10) * (digits_before_dot-1)
    my $dbd_sign = '+';
    if ($dbd < 0)
      {
      $dbd = -$dbd;
      $dbd_sign = '-';
      }
    ($x->{_e}, $x->{_es}) = 
	_e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
 
    }

  # Now: 0.1 <= $x < 10 (and possible correction in l_10)

  ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
  ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)

  $HALF = $class->new($HALF) unless ref($HALF);

  my $twos = 0;				# default: none (0 times)	
  while ($x->bacmp($HALF) <= 0)		# X <= 0.5
    {
    $twos--; $x->bmul($two);
    }
  while ($x->bacmp($two) >= 0)		# X >= 2
    {
    $twos++; $x->bdiv($two,$scale+4);		# keep all digits
    }
  $x->bround($scale+4);
  # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
  # So calculate correction factor based on ln(2):
  if ($twos != 0)
    {
    $LOG_2 = $class->new($LOG_2,undef,undef) unless ref $LOG_2;
    if ($scale <= $LOG_2_A)
      {
      # use cached value
      $l_2 = $LOG_2->copy();			# copy() for the mul below
      }
    else
      {
      # else: slower, compute and cache result
      # also disable downgrade for this code path
      local $Math::BigFloat::downgrade = undef;
      $l_2 = $two->copy(); $class->_log($l_2, $scale); # scale+4, actually
      $LOG_2 = $l_2->copy();			# cache the result for later
						# the copy() is for mul below
      $LOG_2_A = $scale;
      }
    $l_2->bmul($twos);		# * -2 => subtract, * 2 => add
    }
  else
    {
      undef $l_2;
    }
  
  $class->_log($x,$scale);			# need to do the "normal" way
  $x->badd($l_10) if defined $l_10; 		# correct it by ln(10)
  $x->badd($l_2) if defined $l_2;		# and maybe by ln(2)

  # all done, $x contains now the result
  $x;
  }

sub blcm 
  { 
  # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
  # does not modify arguments, but returns new object
  # Lowest Common Multiplicator

  my ($class,@arg) = objectify(0,@_);
  my $x = $class->new(shift @arg);
  while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); } 
  $x;
  }

sub bgcd
  {
  # (BINT or num_str, BINT or num_str) return BINT
  # does not modify arguments, but returns new object

  my $y = shift;
  $y = __PACKAGE__->new($y) if !ref($y);
  my $class = ref($y);
  my $x = $y->copy()->babs();			# keep arguments

  return $x->bnan() if $x->{sign} !~ /^[+-]$/	# x NaN?
	|| !$x->is_int();			# only for integers now

  while (@_)
    {
    my $t = shift; $t = $class->new($t) if !ref($t);
    $y = $t->copy()->babs();
    
    return $x->bnan() if $y->{sign} !~ /^[+-]$/	# y NaN?
     	|| !$y->is_int();			# only for integers now

    # greatest common divisor
    while (! $y->is_zero())
      {
      ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
      }

    last if $x->is_one();
    }
  $x;
  }

##############################################################################

sub _e_add {
    # Internal helper sub to take two positive integers and their signs and
    # then add them. Input ($CALC, $CALC, ('+'|'-'), ('+'|'-')), output
    # ($CALC, ('+'|'-')).

    my ($x, $y, $xs, $ys) = @_;

    # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
    if ($xs eq $ys) {
        $x = $MBI->_add($x, $y);                # +a + +b or -a + -b
    } else {
        my $a = $MBI->_acmp($x, $y);
        if ($a == 0) {
            # This does NOT modify $x in-place. TODO: Fix this?
            $x = $MBI->_zero();                 # result is 0
            $xs = '+';
            return ($x, $xs);
        }
        if ($a > 0) {
            $x = $MBI->_sub($x, $y);            # abs sub
        } else {                                # a < 0
            $x = $MBI->_sub ( $y, $x, 1 );      # abs sub
            $xs = $ys;
        }
    }

    $xs = '+' if $xs eq '-' && $MBI->_is_zero($x);      # no "-0"

    return ($x, $xs);
}

sub _e_sub {
  # Internal helper sub to take two positive integers and their signs and
  # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')), 
  # output ($CALC,('+'|'-'))
  my ($x,$y,$xs,$ys) = @_;

  # flip sign
    $ys = $ys eq '+' ? '-' : '+';       # swap sign of second operand ...
    _e_add($x, $y, $xs, $ys);           # ... and let _e_add() do the job
  }

###############################################################################
# is_foo methods (is_negative, is_positive are inherited from BigInt)

sub is_int
  {
  # return true if arg (BFLOAT or num_str) is an integer
  my ($class,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);

  (($x->{sign} =~ /^[+-]$/) &&			# NaN and +-inf aren't
   ($x->{_es} eq '+')) ? 1 : 0;			# 1e-1 => no integer
  }

sub is_zero
  {
  # return true if arg (BFLOAT or num_str) is zero
  my ($class,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);

  ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m})) ? 1 : 0;
  }

sub is_one
  {
  # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
  my ($class,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);

  $sign = '+' if !defined $sign || $sign ne '-';

  ($x->{sign} eq $sign && 
   $MBI->_is_zero($x->{_e}) &&
   $MBI->_is_one($x->{_m}) ) ? 1 : 0; 
  }

sub is_odd
  {
  # return true if arg (BFLOAT or num_str) is odd or false if even
  my ($class,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
  
  (($x->{sign} =~ /^[+-]$/) &&		# NaN & +-inf aren't
   ($MBI->_is_zero($x->{_e})) &&
   ($MBI->_is_odd($x->{_m}))) ? 1 : 0; 
  }

sub is_even
  {
  # return true if arg (BINT or num_str) is even or false if odd
  my ($class,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);

  (($x->{sign} =~ /^[+-]$/) &&			# NaN & +-inf aren't
   ($x->{_es} eq '+') &&	 		# 123.45 isn't
   ($MBI->_is_even($x->{_m}))) ? 1 : 0;		# but 1200 is
  }

sub bmul
  { 
  # multiply two numbers
  
  # set up parameters
  my ($class,$x,$y,@r) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($class,$x,$y,@r) = objectify(2,@_);
    }

  return $x if $x->modify('bmul');

  return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));

  # inf handling
  if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
    {
    return $x->bnan() if $x->is_zero() || $y->is_zero(); 
    # result will always be +-inf:
    # +inf * +/+inf => +inf, -inf * -/-inf => +inf
    # +inf * -/-inf => -inf, -inf * +/+inf => -inf
    return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
    return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
    return $x->binf('-');
    }
  
  return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
   ((!$x->isa($class)) || (!$y->isa($class)));

  # aEb * cEd = (a*c)E(b+d)
  $MBI->_mul($x->{_m},$y->{_m});
  ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});

  $r[3] = $y;				# no push!

  # adjust sign:
  $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
  $x->bnorm->round(@r);
  }

sub bmuladd
  { 
  # multiply two numbers and add the third to the result
  
  # set up parameters
  my ($class,$x,$y,$z,@r) = objectify(3,@_);

  return $x if $x->modify('bmuladd');

  return $x->bnan() if (($x->{sign} eq $nan) ||
			($y->{sign} eq $nan) ||
			($z->{sign} eq $nan));

  # inf handling
  if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
    {
    return $x->bnan() if $x->is_zero() || $y->is_zero(); 
    # result will always be +-inf:
    # +inf * +/+inf => +inf, -inf * -/-inf => +inf
    # +inf * -/-inf => -inf, -inf * +/+inf => -inf
    return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
    return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
    return $x->binf('-');
    }

  return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
   ((!$x->isa($class)) || (!$y->isa($class)));

  # aEb * cEd = (a*c)E(b+d)
  $MBI->_mul($x->{_m},$y->{_m});
  ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});

  $r[3] = $y;				# no push!

  # adjust sign:
  $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';

  # z=inf handling (z=NaN handled above)
  $x->{sign} = $z->{sign}, return $x if $z->{sign} =~ /^[+-]inf$/;

  # take lower of the two e's and adapt m1 to it to match m2
  my $e = $z->{_e};
  $e = $MBI->_zero() if !defined $e;		# if no BFLOAT?
  $e = $MBI->_copy($e);				# make copy (didn't do it yet)

  my $es;

  ($e,$es) = _e_sub($e, $x->{_e}, $z->{_es} || '+', $x->{_es});

  my $add = $MBI->_copy($z->{_m});

  if ($es eq '-')				# < 0
    {
    $MBI->_lsft( $x->{_m}, $e, 10);
    ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
    }
  elsif (!$MBI->_is_zero($e))			# > 0
    {
    $MBI->_lsft($add, $e, 10);
    }
  # else: both e are the same, so just leave them

  if ($x->{sign} eq $z->{sign})
    {
    # add
    $x->{_m} = $MBI->_add($x->{_m}, $add);
    }
  else
    {
    ($x->{_m}, $x->{sign}) = 
     _e_add($x->{_m}, $add, $x->{sign}, $z->{sign});
    }

  # delete trailing zeros, then round
  $x->bnorm()->round(@r);
  }

sub bdiv 
  {
  # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return 
    # (BFLOAT, BFLOAT) (quo, rem) or BFLOAT (only quo)

  # set up parameters
  my ($class,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($class,$x,$y,$a,$p,$r) = objectify(2,@_);
    }

  return $x if $x->modify('bdiv');

    my $wantarray = wantarray;          # call only once

    # At least one argument is NaN. This is handled the same way as in
    # Math::BigInt -> bdiv().

    if ($x -> is_nan() || $y -> is_nan()) {
        return $wantarray ? ($x -> bnan(), $class -> bnan()) : $x -> bnan();
    }

    # Divide by zero and modulo zero. This is handled the same way as in
    # Math::BigInt -> bdiv(). See the comment in the code for Math::BigInt ->
    # bdiv() for further details.

    if ($y -> is_zero()) {
        my ($quo, $rem);
        if ($wantarray) {
            $rem = $x -> copy();
        }
        if ($x -> is_zero()) {
            $quo = $x -> bnan();
        } else {
            $quo = $x -> binf($x -> {sign});
        }
        return $wantarray ? ($quo, $rem) : $quo;
    }

    # Numerator (dividend) is +/-inf. This is handled the same way as in
    # Math::BigInt -> bdiv(). See the comment in the code for Math::BigInt ->
    # bdiv() for further details.

    if ($x -> is_inf()) {
        my ($quo, $rem);
        $rem = $class -> bnan() if $wantarray;
        if ($y -> is_inf()) {
            $quo = $x -> bnan();
        } else {
            my $sign = $x -> bcmp(0) == $y -> bcmp(0) ? '+' : '-';
            $quo = $x -> binf($sign);
        }
        return $wantarray ? ($quo, $rem) : $quo;
    }

  # Denominator (divisor) is +/-inf. This is handled the same way as in
  # Math::BigInt -> bdiv(), with one exception: In scalar context,
  # Math::BigFloat does true division (although rounded), not floored division
  # (F-division), so a finite number divided by +/-inf is always zero. See the
  # comment in the code for Math::BigInt -> bdiv() for further details.

  if ($y -> is_inf()) {
      my ($quo, $rem);
      if ($wantarray) {
          if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) {
              $rem = $x -> copy();
              $quo = $x -> bzero();
          } else {
              $rem = $class -> binf($y -> {sign});
              $quo = $x -> bone('-');
          }
          return ($quo, $rem);
      } else {
        if ($y -> is_inf()) {
            if ($x -> is_nan() || $x -> is_inf()) {
                return $x -> bnan();
            } else {
                return $x -> bzero();
            }
        }
    }
    }

  # At this point, both the numerator and denominator are finite numbers, and
  # the denominator (divisor) is non-zero.

  # x == 0?
  return wantarray ? ($x,$class->bzero()) : $x if $x->is_zero();

  # upgrade ?
  return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;

  # we need to limit the accuracy to protect against overflow
  my $fallback = 0;
  my (@params,$scale);
  ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);

  return $x if $x->is_nan();		# error in _find_round_parameters?

  # no rounding at all, so must use fallback
    if (scalar @params == 0)
    {
    # simulate old behaviour
    $params[0] = $class->div_scale();	# and round to it as accuracy
    $scale = $params[0]+4; 		# at least four more for proper round
    $params[2] = $r;			# round mode by caller or undef
    $fallback = 1;			# to clear a/p afterwards
    } else {
    # the 4 below is empirical, and there might be cases where it is not
    # enough...
    $scale = abs($params[0] || $params[1]) + 4;	# take whatever is defined
    }

    my $rem;
    $rem = $class -> bzero() if wantarray;

  $y = $class->new($y) unless $y->isa('Math::BigFloat');

  my $lx = $MBI -> _len($x->{_m}); my $ly = $MBI -> _len($y->{_m});
  $scale = $lx if $lx > $scale;
  $scale = $ly if $ly > $scale;
  my $diff = $ly - $lx;
  $scale += $diff if $diff > 0;		# if lx << ly, but not if ly << lx!

  # check that $y is not 1 nor -1 and cache the result:
  my $y_not_one = !($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m}));

  # flipping the sign of $y will also flip the sign of $x for the special
  # case of $x->bsub($x); so we can catch it below:
  my $xsign = $x->{sign};
  $y->{sign} =~ tr/+-/-+/;

  if ($xsign ne $x->{sign})
    {
    # special case of $x /= $x results in 1
    $x->bone();			# "fixes" also sign of $y, since $x is $y
    }
  else
    {
    # correct $y's sign again
    $y->{sign} =~ tr/+-/-+/;
    # continue with normal div code:

    # make copy of $x in case of list context for later remainder calculation
    if (wantarray && $y_not_one)
      {
      $rem = $x->copy();
      }

    $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+'; 

    # check for / +-1 ( +/- 1E0)
    if ($y_not_one)
      {
      # promote BigInts and it's subclasses (except when already a BigFloat)
      $y = $class->new($y) unless $y->isa('Math::BigFloat'); 

      # calculate the result to $scale digits and then round it
      # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
      $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
      $MBI->_div ($x->{_m},$y->{_m});	# a/c

      # correct exponent of $x
      ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
      # correct for 10**scale
      ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
      $x->bnorm();		# remove trailing 0's
      }
    } # end else $x != $y

  # shortcut to not run through _find_round_parameters again
  if (defined $params[0])
    {
    delete $x->{_a}; 				# clear before round
    $x->bround($params[0],$params[2]);		# then round accordingly
    }
  else
    {
    delete $x->{_p}; 				# clear before round
    $x->bfround($params[1],$params[2]);		# then round accordingly
    }
  if ($fallback)
    {
    # clear a/p after round, since user did not request it
    delete $x->{_a}; delete $x->{_p};
    }

  if (wantarray)
    {
    if ($y_not_one)
      {
      $x -> bfloor();
      $rem->bmod($y,@params);			# copy already done
      }
    if ($fallback)
      {
      # clear a/p after round, since user did not request it
      delete $rem->{_a}; delete $rem->{_p};
      }
    return ($x,$rem);
    }
  $x;
  }

sub bmod 
  {
  # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return remainder

  # set up parameters
  my ($class,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($class,$x,$y,$a,$p,$r) = objectify(2,@_);
    }

  return $x if $x->modify('bmod');

    # At least one argument is NaN. This is handled the same way as in
    # Math::BigInt -> bmod().

    if ($x -> is_nan() || $y -> is_nan()) {
        return $x -> bnan();
    } 

    # Modulo zero. This is handled the same way as in Math::BigInt -> bmod().

    if ($y -> is_zero()) {
    return $x;
    }

    # Numerator (dividend) is +/-inf. This is handled the same way as in
    # Math::BigInt -> bmod().

    if ($x -> is_inf()) {
        return $x -> bnan();
    }

    # Denominator (divisor) is +/-inf. This is handled the same way as in
    # Math::BigInt -> bmod().

    if ($y -> is_inf()) {
        if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) {
            return $x;
        } else {
            return $x -> binf($y -> sign());
        }
    }

  return $x->bzero() if $x->is_zero()
 || ($x->is_int() &&
  # check that $y == +1 or $y == -1:
    ($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m})));

  my $cmp = $x->bacmp($y);			# equal or $x < $y?
    if ($cmp == 0) {                            # $x == $y => result 0
        return $x -> bzero($a, $p);
    }

  # only $y of the operands negative? 
    my $neg = $x->{sign} ne $y->{sign} ? 1 : 0;

  $x->{sign} = $y->{sign};				# calc sign first
    if ($cmp < 0 && $neg == 0) {                # $x < $y => result $x
        return $x -> round($a, $p, $r);
    }
  
  my $ym = $MBI->_copy($y->{_m});
  
  # 2e1 => 20
  $MBI->_lsft( $ym, $y->{_e}, 10) 
   if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
 
  # if $y has digits after dot
  my $shifty = 0;			# correct _e of $x by this
  if ($y->{_es} eq '-')			# has digits after dot
    {
    # 123 % 2.5 => 1230 % 25 => 5 => 0.5
    $shifty = $MBI->_num($y->{_e}); 	# no more digits after dot
    $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
    }
  # $ym is now mantissa of $y based on exponent 0

  my $shiftx = 0;			# correct _e of $x by this
  if ($x->{_es} eq '-')			# has digits after dot
    {
    # 123.4 % 20 => 1234 % 200
    $shiftx = $MBI->_num($x->{_e});	# no more digits after dot
    $MBI->_lsft($ym, $x->{_e}, 10);	# 123 => 1230
    }
  # 123e1 % 20 => 1230 % 20
  if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
    {
    $MBI->_lsft( $x->{_m}, $x->{_e},10);	# es => '+' here
    }

  $x->{_e} = $MBI->_new($shiftx);
  $x->{_es} = '+'; 
  $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
  $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
  
  # now mantissas are equalized, exponent of $x is adjusted, so calc result

  $x->{_m} = $MBI->_mod( $x->{_m}, $ym);

  $x->{sign} = '+' if $MBI->_is_zero($x->{_m});		# fix sign for -0
  $x->bnorm();

    if ($neg != 0 && ! $x -> is_zero())   # one of them negative => correct in place
    {
    my $r = $y - $x;
    $x->{_m} = $r->{_m};
    $x->{_e} = $r->{_e};
    $x->{_es} = $r->{_es};
    $x->{sign} = '+' if $MBI->_is_zero($x->{_m});	# fix sign for -0
    $x->bnorm();
    }

  $x->round($a,$p,$r,$y);	# round and return
  }

sub broot
  {
  # calculate $y'th root of $x
  
  # set up parameters
  my ($class,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($class,$x,$y,$a,$p,$r) = objectify(2,@_);
    }

  return $x if $x->modify('broot');

  # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
  return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
         $y->{sign} !~ /^\+$/;

  return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
  
  # we need to limit the accuracy to protect against overflow
  my $fallback = 0;
  my (@params,$scale);
  ($x,@params) = $x->_find_round_parameters($a,$p,$r);

  return $x if $x->is_nan();		# error in _find_round_parameters?

  # no rounding at all, so must use fallback
  if (scalar @params == 0) 
    {
    # simulate old behaviour
    $params[0] = $class->div_scale();	# and round to it as accuracy
    $scale = $params[0]+4; 		# at least four more for proper round
    $params[2] = $r;			# round mode by caller or undef
    $fallback = 1;			# to clear a/p afterwards
    }
  else
    {
    # the 4 below is empirical, and there might be cases where it is not
    # enough...
    $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
    }

  # when user set globals, they would interfere with our calculation, so
  # disable them and later re-enable them
  no strict 'refs';
  my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef;
  my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef;
  # we also need to disable any set A or P on $x (_find_round_parameters took
  # them already into account), since these would interfere, too
  delete $x->{_a}; delete $x->{_p};
  # need to disable $upgrade in BigInt, to avoid deep recursion
  local $Math::BigInt::upgrade = undef;	# should be really parent class vs MBI

  # remember sign and make $x positive, since -4 ** (1/2) => -2
  my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';

  my $is_two = 0;
  if ($y->isa('Math::BigFloat'))
    {
    $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
    }
  else
    {
    $is_two = ($y == 2);
    }

  # normal square root if $y == 2:
  if ($is_two)
    {
    $x->bsqrt($scale+4);
    }
  elsif ($y->is_one('-'))
    {
    # $x ** -1 => 1/$x
    my $u = $class->bone()->bdiv($x,$scale);
    # copy private parts over
    $x->{_m} = $u->{_m};
    $x->{_e} = $u->{_e};
    $x->{_es} = $u->{_es};
    }
  else
    {
    # calculate the broot() as integer result first, and if it fits, return
    # it rightaway (but only if $x and $y are integer):

    my $done = 0;				# not yet
    if ($y->is_int() && $x->is_int())
      {
      my $i = $MBI->_copy( $x->{_m} );
      $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
      my $int = Math::BigInt->bzero();
      $int->{value} = $i;
      $int->broot($y->as_number());
      # if ($exact)
      if ($int->copy()->bpow($y) == $x)
        {
        # found result, return it
        $x->{_m} = $int->{value};
        $x->{_e} = $MBI->_zero();
        $x->{_es} = '+';
        $x->bnorm();
        $done = 1;
        }
      }
    if ($done == 0)
      {
      my $u = $class->bone()->bdiv($y,$scale+4);
      delete $u->{_a}; delete $u->{_p};         # otherwise it conflicts
      $x->bpow($u,$scale+4);                    # el cheapo
      }
    }
  $x->bneg() if $sign == 1;
  
  # shortcut to not run through _find_round_parameters again
  if (defined $params[0])
    {
    $x->bround($params[0],$params[2]);		# then round accordingly
    }
  else
    {
    $x->bfround($params[1],$params[2]);		# then round accordingly
    }
  if ($fallback)
    {
    # clear a/p after round, since user did not request it
    delete $x->{_a}; delete $x->{_p};
    }
  # restore globals
  $$abr = $ab; $$pbr = $pb;
  $x;
  }

sub bsqrt
  { 
  # calculate square root
  my ($class,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);

  return $x if $x->modify('bsqrt');

  return $x->bnan() if $x->{sign} !~ /^[+]/;	# NaN, -inf or < 0
  return $x if $x->{sign} eq '+inf';		# sqrt(inf) == inf
  return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();

  # we need to limit the accuracy to protect against overflow
  my $fallback = 0;
  my (@params,$scale);
  ($x,@params) = $x->_find_round_parameters($a,$p,$r);

  return $x if $x->is_nan();		# error in _find_round_parameters?

  # no rounding at all, so must use fallback
  if (scalar @params == 0) 
    {
    # simulate old behaviour
    $params[0] = $class->div_scale();	# and round to it as accuracy
    $scale = $params[0]+4; 		# at least four more for proper round
    $params[2] = $r;			# round mode by caller or undef
    $fallback = 1;			# to clear a/p afterwards
    }
  else
    {
    # the 4 below is empirical, and there might be cases where it is not
    # enough...
    $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
    }

  # when user set globals, they would interfere with our calculation, so
  # disable them and later re-enable them
  no strict 'refs';
  my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef;
  my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef;
  # we also need to disable any set A or P on $x (_find_round_parameters took
  # them already into account), since these would interfere, too
  delete $x->{_a}; delete $x->{_p};
  # need to disable $upgrade in BigInt, to avoid deep recursion
  local $Math::BigInt::upgrade = undef;	# should be really parent class vs MBI

  my $i = $MBI->_copy( $x->{_m} );
  $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
  my $xas = Math::BigInt->bzero();
  $xas->{value} = $i;

  my $gs = $xas->copy()->bsqrt();	# some guess

  if (($x->{_es} ne '-')		# guess can't be accurate if there are
					# digits after the dot
   && ($xas->bacmp($gs * $gs) == 0))	# guess hit the nail on the head?
    {
    # exact result, copy result over to keep $x
    $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
    $x->bnorm();
    # shortcut to not run through _find_round_parameters again
    if (defined $params[0])
      {
      $x->bround($params[0],$params[2]);	# then round accordingly
      }
    else
      {
      $x->bfround($params[1],$params[2]);	# then round accordingly
      }
    if ($fallback)
      {
      # clear a/p after round, since user did not request it
      delete $x->{_a}; delete $x->{_p};
      }
    # re-enable A and P, upgrade is taken care of by "local"
    ${"$class\::accuracy"} = $ab; ${"$class\::precision"} = $pb;
    return $x;
    }
 
  # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
  # of the result by multiplying the input by 100 and then divide the integer
  # result of sqrt(input) by 10. Rounding afterwards returns the real result.

  # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
  my $y1 = $MBI->_copy($x->{_m});

  my $length = $MBI->_len($y1);
  
  # Now calculate how many digits the result of sqrt(y1) would have
  my $digits = int($length / 2);

  # But we need at least $scale digits, so calculate how many are missing
  my $shift = $scale - $digits;

  # This happens if the input had enough digits
  # (we take care of integer guesses above)
  $shift = 0 if $shift < 0; 

  # Multiply in steps of 100, by shifting left two times the "missing" digits
  my $s2 = $shift * 2;

  # We now make sure that $y1 has the same odd or even number of digits than
  # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
  # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
  # steps of 10. The length of $x does not count, since an even or odd number
  # of digits before the dot is not changed by adding an even number of digits
  # after the dot (the result is still odd or even digits long).
  $s2++ if $MBI->_is_odd($x->{_e});

  $MBI->_lsft( $y1, $MBI->_new($s2), 10);

  # now take the square root and truncate to integer
  $y1 = $MBI->_sqrt($y1);

  # By "shifting" $y1 right (by creating a negative _e) we calculate the final
  # result, which is than later rounded to the desired scale.

  # calculate how many zeros $x had after the '.' (or before it, depending
  # on sign of $dat, the result should have half as many:
  my $dat = $MBI->_num($x->{_e});
  $dat = -$dat if $x->{_es} eq '-';
  $dat += $length;

  if ($dat > 0)
    {
    # no zeros after the dot (e.g. 1.23, 0.49 etc)
    # preserve half as many digits before the dot than the input had 
    # (but round this "up")
    $dat = int(($dat+1)/2);
    }
  else
    {
    $dat = int(($dat)/2);
    }
  $dat -= $MBI->_len($y1);
  if ($dat < 0)
    {
    $dat = abs($dat);
    $x->{_e} = $MBI->_new( $dat );
    $x->{_es} = '-';
    }
  else
    {    
    $x->{_e} = $MBI->_new( $dat );
    $x->{_es} = '+';
    }
  $x->{_m} = $y1;
  $x->bnorm();

  # shortcut to not run through _find_round_parameters again
  if (defined $params[0])
    {
    $x->bround($params[0],$params[2]);		# then round accordingly
    }
  else
    {
    $x->bfround($params[1],$params[2]);		# then round accordingly
    }
  if ($fallback)
    {
    # clear a/p after round, since user did not request it
    delete $x->{_a}; delete $x->{_p};
    }
  # restore globals
  $$abr = $ab; $$pbr = $pb;
  $x;
  }

sub bfac
  {
  # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
  # compute factorial number, modifies first argument

  # set up parameters
  my ($class,$x,@r) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  ($class,$x,@r) = objectify(1,@_) if !ref($x);

  # inf => inf
  return $x if $x->modify('bfac') || $x->{sign} eq '+inf';	

  return $x->bnan() 
    if (($x->{sign} ne '+') ||		# inf, NaN, <0 etc => NaN
     ($x->{_es} ne '+'));		# digits after dot?

  # use BigInt's bfac() for faster calc
  if (! $MBI->_is_zero($x->{_e}))
    {
    $MBI->_lsft($x->{_m}, $x->{_e},10);	# change 12e1 to 120e0
    $x->{_e} = $MBI->_zero();		# normalize
    $x->{_es} = '+';
    }
  $MBI->_fac($x->{_m});			# calculate factorial
  $x->bnorm()->round(@r); 		# norm again and round result
  }

sub _pow
  {
  # Calculate a power where $y is a non-integer, like 2 ** 0.3
  my ($x,$y,@r) = @_;
  my $class = ref($x);

  # if $y == 0.5, it is sqrt($x)
  $HALF = $class->new($HALF) unless ref($HALF);
  return $x->bsqrt(@r,$y) if $y->bcmp($HALF) == 0;

  # Using:
  # a ** x == e ** (x * ln a)

  # u = y * ln x
  #                _                         _
  # Taylor:       |   u    u^2    u^3         |
  # x ** y  = 1 + |  --- + --- + ----- + ...  |
  #               |_  1    1*2   1*2*3       _|

  # we need to limit the accuracy to protect against overflow
  my $fallback = 0;
  my ($scale,@params);
  ($x,@params) = $x->_find_round_parameters(@r);
    
  return $x if $x->is_nan();		# error in _find_round_parameters?

  # no rounding at all, so must use fallback
  if (scalar @params == 0)
    {
    # simulate old behaviour
    $params[0] = $class->div_scale();	# and round to it as accuracy
    $params[1] = undef;			# disable P
    $scale = $params[0]+4; 		# at least four more for proper round
    $params[2] = $r[2];			# round mode by caller or undef
    $fallback = 1;			# to clear a/p afterwards
    }
  else
    {
    # the 4 below is empirical, and there might be cases where it is not
    # enough...
    $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
    }

  # when user set globals, they would interfere with our calculation, so
  # disable them and later re-enable them
  no strict 'refs';
  my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef;
  my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef;
  # we also need to disable any set A or P on $x (_find_round_parameters took
  # them already into account), since these would interfere, too
  delete $x->{_a}; delete $x->{_p};
  # need to disable $upgrade in BigInt, to avoid deep recursion
  local $Math::BigInt::upgrade = undef;
 
  my ($limit,$v,$u,$below,$factor,$next,$over);

  $u = $x->copy()->blog(undef,$scale)->bmul($y);
  my $do_invert = ($u->{sign} eq '-');
  $u->bneg()  if $do_invert;
  $v = $class->bone();				# 1
  $factor = $class->new(2);			# 2
  $x->bone();					# first term: 1

  $below = $v->copy();
  $over = $u->copy();

  $limit = $class->new("1E-". ($scale-1));
  #my $steps = 0;
  while (3 < 5)
    {
    # we calculate the next term, and add it to the last
    # when the next term is below our limit, it won't affect the outcome
    # anymore, so we stop:
    $next = $over->copy()->bdiv($below,$scale);
    last if $next->bacmp($limit) <= 0;
    $x->badd($next);
    # calculate things for the next term
    $over *= $u; $below *= $factor; $factor->binc();

    last if $x->{sign} !~ /^[-+]$/;

    #$steps++;
    }

  if ($do_invert)
    {
    my $x_copy = $x->copy;
    $x->bone->bdiv($x_copy, $scale);
    }
  
  # shortcut to not run through _find_round_parameters again
  if (defined $params[0])
    {
    $x->bround($params[0],$params[2]);		# then round accordingly
    }
  else
    {
    $x->bfround($params[1],$params[2]);		# then round accordingly
    }
  if ($fallback)
    {
    # clear a/p after round, since user did not request it
    delete $x->{_a}; delete $x->{_p};
    }
  # restore globals
  $$abr = $ab; $$pbr = $pb;
  $x;
  }

sub bpow 
  {
  # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
  # compute power of two numbers, second arg is used as integer
  # modifies first argument

  # set up parameters
  my ($class,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
  # objectify is costly, so avoid it
  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
    {
    ($class,$x,$y,$a,$p,$r) = objectify(2,@_);
    }

  return $x if $x->modify('bpow');

  return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
  return $x if $x->{sign} =~ /^[+-]inf$/;
  
  # cache the result of is_zero
  my $y_is_zero = $y->is_zero();
  return $x->bone() if $y_is_zero;
  return $x         if $x->is_one() || $y->is_one();

  my $x_is_zero = $x->is_zero();
  return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int();		# non-integer power

  my $y1 = $y->as_number()->{value};			# make MBI part

  # if ($x == -1)
  if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
    {
    # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
    return $MBI->_is_odd($y1) ? $x : $x->babs(1);
    }
  if ($x_is_zero)
    {
    return $x if $y->{sign} eq '+'; 	# 0**y => 0 (if not y <= 0)
    # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
    return $x->binf();
    }

  my $new_sign = '+';
  $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';

  # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
  $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
  $x->{_e} = $MBI->_mul ($x->{_e}, $y1);

  $x->{sign} = $new_sign;
  $x->bnorm();
  if ($y->{sign} eq '-')
    {
    # modify $x in place!
    my $z = $x->copy(); $x->bone();
    return scalar $x->bdiv($z,$a,$p,$r);	# round in one go (might ignore y's A!)
    }
  $x->round($a,$p,$r,$y);
  }

sub bmodpow
  {
  # takes a very large number to a very large exponent in a given very
  # large modulus, quickly, thanks to binary exponentiation. Supports
  # negative exponents.
  my ($class,$num,$exp,$mod,@r) = objectify(3,@_);

  return $num if $num->modify('bmodpow');

  # check modulus for valid values
  return $num->bnan() if ($mod->{sign} ne '+'           # NaN, - , -inf, +inf
                       || $mod->is_zero());

  # check exponent for valid values
  if ($exp->{sign} =~ /\w/)
    {
    # i.e., if it's NaN, +inf, or -inf...
    return $num->bnan();
    }

  $num->bmodinv ($mod) if ($exp->{sign} eq '-');

  # check num for valid values (also NaN if there was no inverse but $exp < 0)
  return $num->bnan() if $num->{sign} !~ /^[+-]$/;

  # $mod is positive, sign on $exp is ignored, result also positive

  # XXX TODO: speed it up when all three numbers are integers
  $num->bpow($exp)->bmod($mod);
  }

###############################################################################
# trigonometric functions

# helper function for bpi() and batan2(), calculates arcus tanges (1/x)

sub _atan_inv
  {
  # return a/b so that a/b approximates atan(1/x) to at least limit digits
  my ($class, $x, $limit) = @_;

  # Taylor:       x^3   x^5   x^7   x^9
  #    atan = x - --- + --- - --- + --- - ...
  #                3     5     7     9 

  #               1      1         1        1
  #    atan 1/x = - - ------- + ------- - ------- + ...
  #               x   x^3 * 3   x^5 * 5   x^7 * 7 

  #               1      1         1            1
  #    atan 1/x = - - --------- + ---------- - ----------- + ... 
  #               5    3 * 125     5 * 3125     7 * 78125

  # Subtraction/addition of a rational:

  #  5    7    5*3 +- 7*4
  #  - +- -  = ----------
  #  4    3       4*3

  # Term:  N        N+1
  #
  #        a             1                  a * d * c +- b
  #        ----- +- ------------------  =  ----------------
  #        b           d * c                b * d * c

  #  since b1 = b0 * (d-2) * c

  #        a             1                  a * d +- b / c
  #        ----- +- ------------------  =  ----------------
  #        b           d * c                b * d 

  # and  d = d + 2
  # and  c = c * x * x

  #        u = d * c
  #        stop if length($u) > limit 
  #        a = a * u +- b
  #        b = b * u
  #        d = d + 2
  #        c = c * x * x
  #        sign = 1 - sign

  my $a = $MBI->_one();
  my $b = $MBI->_copy($x);
 
  my $x2  = $MBI->_mul( $MBI->_copy($x), $b);		# x2 = x * x
  my $d   = $MBI->_new( 3 );				# d = 3
  my $c   = $MBI->_mul( $MBI->_copy($x), $x2);		# c = x ^ 3
  my $two = $MBI->_new( 2 );

  # run the first step unconditionally
  my $u = $MBI->_mul( $MBI->_copy($d), $c);
  $a = $MBI->_mul($a, $u);
  $a = $MBI->_sub($a, $b);
  $b = $MBI->_mul($b, $u);
  $d = $MBI->_add($d, $two);
  $c = $MBI->_mul($c, $x2);

  # a is now a * (d-3) * c
  # b is now b * (d-2) * c

  # run the second step unconditionally
  $u = $MBI->_mul( $MBI->_copy($d), $c);
  $a = $MBI->_mul($a, $u);
  $a = $MBI->_add($a, $b);
  $b = $MBI->_mul($b, $u);
  $d = $MBI->_add($d, $two);
  $c = $MBI->_mul($c, $x2);

  # a is now a * (d-3) * (d-5) * c * c  
  # b is now b * (d-2) * (d-4) * c * c

  # so we can remove c * c from both a and b to shorten the numbers involved:
  $a = $MBI->_div($a, $x2);
  $b = $MBI->_div($b, $x2);
  $a = $MBI->_div($a, $x2);
  $b = $MBI->_div($b, $x2);

#  my $step = 0; 
  my $sign = 0;						# 0 => -, 1 => +
  while (3 < 5)
    {
#    $step++;
#    if (($i++ % 100) == 0)
#      {
#    print "a=",$MBI->_str($a),"\n";
#    print "b=",$MBI->_str($b),"\n";
#      }
#    print "d=",$MBI->_str($d),"\n";
#    print "x2=",$MBI->_str($x2),"\n";
#    print "c=",$MBI->_str($c),"\n";

    my $u = $MBI->_mul( $MBI->_copy($d), $c);
    # use _alen() for libs like GMP where _len() would be O(N^2)
    last if $MBI->_alen($u) > $limit;
    my ($bc,$r) = $MBI->_div( $MBI->_copy($b), $c);
    if ($MBI->_is_zero($r))
      {
      # b / c is an integer, so we can remove c from all terms
      # this happens almost every time:
      $a = $MBI->_mul($a, $d);
      $a = $MBI->_sub($a, $bc) if $sign == 0;
      $a = $MBI->_add($a, $bc) if $sign == 1;
      $b = $MBI->_mul($b, $d);
      }
    else
      {
      # b / c is not an integer, so we keep c in the terms
      # this happens very rarely, for instance for x = 5, this happens only
      # at the following steps:
      # 1, 5, 14, 32, 72, 157, 340, ...
      $a = $MBI->_mul($a, $u);
      $a = $MBI->_sub($a, $b) if $sign == 0;
      $a = $MBI->_add($a, $b) if $sign == 1;
      $b = $MBI->_mul($b, $u);
      }
    $d = $MBI->_add($d, $two);
    $c = $MBI->_mul($c, $x2);
    $sign = 1 - $sign;

    }

#  print "Took $step steps for ", $MBI->_str($x),"\n";
#  print "a=",$MBI->_str($a),"\n"; print "b=",$MBI->_str($b),"\n";
  # return a/b so that a/b approximates atan(1/x)
  ($a,$b);
  }

sub bpi {

    # Called as                 Argument list
    # ---------                 -------------
    # Math::BigFloat->bpi()     ("Math::BigFloat")
    # Math::BigFloat->bpi(10)   ("Math::BigFloat", 10)
    # $x->bpi()                 ($x)
    # $x->bpi(10)               ($x, 10)
    # Math::BigFloat::bpi()     ()
    # Math::BigFloat::bpi(10)   (10)
    #
    # In ambiguous cases, we favour the OO-style, so the following case
    #
    #   $n = Math::BigFloat->new("10");
    #   $x = Math::BigFloat->bpi($n);
    #
    # which gives an argument list with the single element $n, is resolved as
    #
    #   $n->bpi();

    my $self    = shift;
    my $selfref = ref $self;
    my $class   = $selfref || $self;

    my $accu;                      # accuracy (number of digits)
    my $prec;                      # precision
    my $rndm;                      # round mode

    # If bpi() is called as a function ...
    #
    # This cludge is necessary because we still support bpi() as a function. If
    # bpi() is called with either no argument or one argument, and that one
    # argument is either undefined or a scalar that looks like a number, then
    # we assume bpi() is called as a function.

    if (@_ == 0 &&
        (defined($self) && !ref($self) && $self =~ /^\s*[+-]?\d/i)
          ||
        !defined($self))
    {
        $accu = $self;
        $class = __PACKAGE__;
        $self = $class -> bzero();      # initialize
    }

    # ... or if bpi() is called as a method ...

    else {
        if ($selfref) {                 # bpi() called as instance method
            return $self if $self -> modify('bpi');
        } else {                        # bpi() called as class method
            $self = $class -> bzero();  # initialize
        }
        $accu = shift;
        $prec = shift;
        $rndm = shift;
    }

    my @r = ($accu, $prec, $rndm);

    # We need to limit the accuracy to protect against overflow.
    my $fallback = 0;
    my ($scale, @params);
    ($self, @params) = $self -> _find_round_parameters(@r);

    # Error in _find_round_parameters?
    #
    # We can't return here, because that will fail if $self was a NaN when
    # bpi() was invoked, and we want to assign pi to $x. It is probably not a
    # good idea that _find_round_parameters() signals invalid round parameters
    # by silently returning a NaN. Fixme!
    #return $self if $self && $self->is_nan();

    # No rounding at all, so must use fallback.
    if (scalar @params == 0) {
        # Simulate old behaviour
        $params[0] = $self -> div_scale();  # and round to it as accuracy
        $params[1] = undef;                 # disable P
        $params[2] = $r[2];                 # round mode by caller or undef
        $fallback = 1;                      # to clear a/p afterwards
    }

    # The accuracy, i.e., the number of digits. Pi has one digit before the
    # dot, so a precision of 4 digits is equivalent to an accuracy of 5 digits.

    my $n = $params[0] || 1 - $params[1];

    if ($n < 1000) {

        # after 黃見利 (Hwang Chien-Lih) (1997)
        # pi/4 = 183 * atan(1/239) + 32 * atan(1/1023) – 68 * atan(1/5832)
        #        + 12 * atan(1/110443) - 12 * atan(1/4841182) - 100 * atan(1/6826318)

        # Use a few more digits in the intermediate computations.

        my $nextra = $n < 800 ? 4 : 5;
        $n += $nextra;

        my ($a, $b) = $class->_atan_inv($MBI->_new(239), $n);
        my ($c, $d) = $class->_atan_inv($MBI->_new(1023), $n);
        my ($e, $f) = $class->_atan_inv($MBI->_new(5832), $n);
        my ($g, $h) = $class->_atan_inv($MBI->_new(110443), $n);
        my ($i, $j) = $class->_atan_inv($MBI->_new(4841182), $n);
        my ($k, $l) = $class->_atan_inv($MBI->_new(6826318), $n);

        $MBI->_mul($a, $MBI->_new(732));
        $MBI->_mul($c, $MBI->_new(128));
        $MBI->_mul($e, $MBI->_new(272));
        $MBI->_mul($g, $MBI->_new(48));
        $MBI->_mul($i, $MBI->_new(48));
        $MBI->_mul($k, $MBI->_new(400));

        my $x = $class->bone(); $x->{_m} = $a; my $x_d = $class->bone(); $x_d->{_m} = $b;
        my $y = $class->bone(); $y->{_m} = $c; my $y_d = $class->bone(); $y_d->{_m} = $d;
        my $z = $class->bone(); $z->{_m} = $e; my $z_d = $class->bone(); $z_d->{_m} = $f;
        my $u = $class->bone(); $u->{_m} = $g; my $u_d = $class->bone(); $u_d->{_m} = $h;
        my $v = $class->bone(); $v->{_m} = $i; my $v_d = $class->bone(); $v_d->{_m} = $j;
        my $w = $class->bone(); $w->{_m} = $k; my $w_d = $class->bone(); $w_d->{_m} = $l;
        $x->bdiv($x_d, $n);
        $y->bdiv($y_d, $n);
        $z->bdiv($z_d, $n);
        $u->bdiv($u_d, $n);
        $v->bdiv($v_d, $n);
        $w->bdiv($w_d, $n);

        delete $x->{_a}; delete $y->{_a}; delete $z->{_a};
        delete $u->{_a}; delete $v->{_a}; delete $w->{_a};
        $x->badd($y)->bsub($z)->badd($u)->bsub($v)->bsub($w);

        for my $key (qw/ sign _m _es _e _a _p /) {
            $self -> {$key} = $x -> {$key} if exists $x -> {$key};
        }

    } else {

        # For large accuracy, the arctan formulas become very inefficient with
        # Math::BigFloat. Switch to Brent-Salamin (aka AGM or Gauss-Legendre).

        # Use a few more digits in the intermediate computations.
        my $nextra = 8;

        $HALF = $class -> new($HALF) unless ref($HALF);
        my ($an, $bn, $tn, $pn) = ($class -> bone, $HALF -> copy -> bsqrt($n),
                                   $HALF -> copy -> bmul($HALF), $class -> bone);
        while ($pn < $n) {
            my $prev_an = $an -> copy;
            $an -> badd($bn) -> bmul($HALF, $n);
            $bn -> bmul($prev_an) -> bsqrt($n);
            $prev_an -> bsub($an);
            $tn -> bsub($pn * $prev_an * $prev_an);
            $pn -> badd($pn);
        }
        $an -> badd($bn);
        $an -> bmul($an, $n) -> bdiv(4 * $tn, $n);

        for my $key (qw/ sign _m _es _e _a _p /) {
            $self -> {$key} = $an -> {$key} if exists $an -> {$key};;
        }
    }

    $self -> round(@params);

    if ($fallback) {
        delete $self->{_a};
        delete $self->{_p};
    }

    return $self;
}

sub bcos
  {
  # Calculate a cosinus of x.
  my ($class,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);

  # Taylor:      x^2   x^4   x^6   x^8
  #    cos = 1 - --- + --- - --- + --- ...
  #               2!    4!    6!    8!

  # we need to limit the accuracy to protect against overflow
  my $fallback = 0;
  my ($scale,@params);
  ($x,@params) = $x->_find_round_parameters(@r);
    
  #         constant object       or error in _find_round_parameters?
  return $x if $x->modify('bcos') || $x->is_nan();

  return $x->bone(@r) if $x->is_zero();

  # no rounding at all, so must use fallback
  if (scalar @params == 0)
    {
    # simulate old behaviour
    $params[0] = $class->div_scale();	# and round to it as accuracy
    $params[1] = undef;			# disable P
    $scale = $params[0]+4; 		# at least four more for proper round
    $params[2] = $r[2];			# round mode by caller or undef
    $fallback = 1;			# to clear a/p afterwards
    }
  else
    {
    # the 4 below is empirical, and there might be cases where it is not
    # enough...
    $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
    }

  # when user set globals, they would interfere with our calculation, so
  # disable them and later re-enable them
  no strict 'refs';
  my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef;
  my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef;
  # we also need to disable any set A or P on $x (_find_round_parameters took
  # them already into account), since these would interfere, too
  delete $x->{_a}; delete $x->{_p};
  # need to disable $upgrade in BigInt, to avoid deep recursion
  local $Math::BigInt::upgrade = undef;
 
  my $last = 0;
  my $over = $x * $x;                   # X ^ 2
  my $x2 = $over->copy();               # X ^ 2; difference between terms
  my $sign = 1;                         # start with -=
  my $below = $class->new(2); my $factorial = $class->new(3);
  $x->bone(); delete $x->{_a}; delete $x->{_p};

  my $limit = $class->new("1E-". ($scale-1));
  #my $steps = 0;
  while (3 < 5)
    {
    # we calculate the next term, and add it to the last
    # when the next term is below our limit, it won't affect the outcome
    # anymore, so we stop:
    my $next = $over->copy()->bdiv($below,$scale);
    last if $next->bacmp($limit) <= 0;

    if ($sign == 0)
      {
      $x->badd($next);
      }
    else
      {
      $x->bsub($next);
      }
    $sign = 1-$sign;					# alternate
    # calculate things for the next term
    $over->bmul($x2);					# $x*$x
    $below->bmul($factorial); $factorial->binc();	# n*(n+1)
    $below->bmul($factorial); $factorial->binc();	# n*(n+1)
    }

  # shortcut to not run through _find_round_parameters again
  if (defined $params[0])
    {
    $x->bround($params[0],$params[2]);		# then round accordingly
    }
  else
    {
    $x->bfround($params[1],$params[2]);		# then round accordingly
    }
  if ($fallback)
    {
    # clear a/p after round, since user did not request it
    delete $x->{_a}; delete $x->{_p};
    }
  # restore globals
  $$abr = $ab; $$pbr = $pb;
  $x;
  }

sub bsin
  {
  # Calculate a sinus of x.
  my ($class,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);

  # taylor:      x^3   x^5   x^7   x^9
  #    sin = x - --- + --- - --- + --- ...
  #               3!    5!    7!    9!

  # we need to limit the accuracy to protect against overflow
  my $fallback = 0;
  my ($scale,@params);
  ($x,@params) = $x->_find_round_parameters(@r);
    
  #         constant object       or error in _find_round_parameters?
  return $x if $x->modify('bsin') || $x->is_nan();

  return $x->bzero(@r) if $x->is_zero();

  # no rounding at all, so must use fallback
  if (scalar @params == 0)
    {
    # simulate old behaviour
    $params[0] = $class->div_scale();	# and round to it as accuracy
    $params[1] = undef;			# disable P
    $scale = $params[0]+4; 		# at least four more for proper round
    $params[2] = $r[2];			# round mode by caller or undef
    $fallback = 1;			# to clear a/p afterwards
    }
  else
    {
    # the 4 below is empirical, and there might be cases where it is not
    # enough...
    $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
    }

  # when user set globals, they would interfere with our calculation, so
  # disable them and later re-enable them
  no strict 'refs';
  my $abr = "$class\::accuracy"; my $ab = $$abr; $$abr = undef;
  my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef;
  # we also need to disable any set A or P on $x (_find_round_parameters took
  # them already into account), since these would interfere, too
  delete $x->{_a}; delete $x->{_p};
  # need to disable $upgrade in BigInt, to avoid deep recursion
  local $Math::BigInt::upgrade = undef;
 
  my $last = 0;
  my $over = $x * $x;			# X ^ 2
  my $x2 = $over->copy();		# X ^ 2; difference between terms
  $over->bmul($x);			# X ^ 3 as starting value
  my $sign = 1;				# start with -=
  my $below = $class->new(6); my $factorial = $class->new(4);
  delete $x->{_a}; delete $x->{_p};

  my $limit = $class->new("1E-". ($scale-1));
  #my $steps = 0;
  while (3 < 5)
    {
    # we calculate the next term, and add it to the last
    # when the next term is below our limit, it won't affect the outcome
    # anymore, so we stop:
    my $next = $over->copy()->bdiv($below,$scale);
    last if $next->bacmp($limit) <= 0;

    if ($sign == 0)
      {
      $x->badd($next);
      }
    else
      {
      $x->bsub($next);
      }
    $sign = 1-$sign;					# alternate
    # calculate things for the next term
    $over->bmul($x2);					# $x*$x
    $below->bmul($factorial); $factorial->binc();	# n*(n+1)
    $below->bmul($factorial); $factorial->binc();	# n*(n+1)
    }

  # shortcut to not run through _find_round_parameters again
  if (defined $params[0])
    {
    $x->bround($params[0],$params[2]);		# then round accordingly
    }
  else
    {
    $x->bfround($params[1],$params[2]);		# then round accordingly
    }
  if ($fallback)
    {
    # clear a/p after round, since user did not request it
    delete $x->{_a}; delete $x->{_p};
    }
  # restore globals
  $$abr = $ab; $$pbr = $pb;
  $x;
  }

sub batan2 {
    # $y -> batan2($x) returns the arcus tangens of $y / $x.

    # Set up parameters.
    my ($class, $y, $x, @r) = (ref($_[0]), @_);

    # Objectify is costly, so avoid it if we can.
    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) {
        ($class, $y, $x, @r) = objectify(2, @_);
    }

    # Quick exit if $y is read-only.
    return $y if $y -> modify('batan2');

    # Handle all NaN cases.
    return $y -> bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;

    # We need to limit the accuracy to protect against overflow.
    my $fallback = 0;
    my ($scale, @params);
    ($y, @params) = $y -> _find_round_parameters(@r);

    # Error in _find_round_parameters?
    return $y if $y->is_nan();

    # No rounding at all, so must use fallback.
    if (scalar @params == 0) {
        # Simulate old behaviour
        $params[0] = $class -> div_scale();  # and round to it as accuracy
        $params[1] = undef;                 # disable P
        $scale = $params[0] + 4;            # at least four more for proper round
        $params[2] = $r[2];                 # round mode by caller or undef
        $fallback = 1;                      # to clear a/p afterwards
    } else {
        # The 4 below is empirical, and there might be cases where it is not
        # enough ...
        $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
    }

    if ($x -> is_inf("+")) {                            # x = inf
        if ($y -> is_inf("+")) {                        #    y = inf
            $y -> bpi($scale) -> bmul("0.25");          #       pi/4
        } elsif ($y -> is_inf("-")) {                   #    y = -inf
            $y -> bpi($scale) -> bmul("-0.25");         #       -pi/4
        } else {                                        #    -inf < y < inf
            return $y -> bzero(@r);                     #       0
        }
    }

    elsif ($x -> is_inf("-")) {                         # x = -inf
        if ($y -> is_inf("+")) {                        #    y = inf
            $y -> bpi($scale) -> bmul("0.75");          #       3/4 pi
        } elsif ($y -> is_inf("-")) {                   #    y = -inf
            $y -> bpi($scale) -> bmul("-0.75");         #       -3/4 pi
        } elsif ($y >= 0) {                             #    y >= 0
            $y -> bpi($scale);                          #       pi
        } else {                                        #    y < 0
            $y -> bpi($scale) -> bneg();                #       -pi
        }
    }

    elsif ($x > 0) {                                    # 0 < x < inf
        if ($y -> is_inf("+")) {                        #    y = inf
            $y -> bpi($scale) -> bmul("0.5");           #       pi/2
        } elsif ($y -> is_inf("-")) {                   #    y = -inf
            $y -> bpi($scale) -> bmul("-0.5");          #       -pi/2
        } else {                                        #   -inf < y < inf
            $y -> bdiv($x, $scale) -> batan($scale);    #       atan(y/x)
        }
    }

    elsif ($x < 0) {                                    # -inf < x < 0
        my $pi = $class -> bpi($scale);
        if ($y >= 0) {                                  #    y >= 0
            $y -> bdiv($x, $scale) -> batan()           #       atan(y/x) + pi
               -> badd($pi);
        } else {                                        #    y < 0
            $y -> bdiv($x, $scale) -> batan()           #       atan(y/x) - pi
               -> bsub($pi);
        }
    }

    else {                                              # x = 0
        if ($y > 0) {                                   #    y > 0
            $y -> bpi($scale) -> bmul("0.5");           #       pi/2
        } elsif ($y < 0) {                              #    y < 0
            $y -> bpi($scale) -> bmul("-0.5");          #       -pi/2
        } else {                                        #    y = 0
            return $y -> bzero(@r);                     #       0
        }
    }

    $y -> round(@r);

    if ($fallback) {
        delete $y->{_a};
        delete $y->{_p};
    }

    return $y;
}

sub batan {
    # Calculate a arcus tangens of x.

    my $self    = shift;
    my $selfref = ref $self;
    my $class   = $selfref || $self;

    my (@r) = @_;

    # taylor:       x^3   x^5   x^7   x^9
    #    atan = x - --- + --- - --- + --- ...
    #                3     5     7     9

    # We need to limit the accuracy to protect against overflow.

    my $fallback = 0;
    my ($scale, @params);
    ($self, @params) = $self->_find_round_parameters(@r);

    # Constant object or error in _find_round_parameters?

    return $self if $self->modify('batan') || $self->is_nan();

    if ($self->{sign} =~ /^[+-]inf\z/) {
        # +inf result is PI/2
        # -inf result is -PI/2
        # calculate PI/2
        my $pi = $class->bpi(@r);
        # modify $self in place
        $self->{_m} = $pi->{_m};
        $self->{_e} = $pi->{_e};
        $self->{_es} = $pi->{_es};
        # -y => -PI/2, +y => PI/2
        $self->{sign} = substr($self->{sign}, 0, 1);  # "+inf" => "+"
        $MBI->_div($self->{_m}, $MBI->_new(2));
        return $self;
    }

    return $self->bzero(@r) if $self->is_zero();

    # no rounding at all, so must use fallback
    if (scalar @params == 0) {
        # simulate old behaviour
        $params[0] = $class->div_scale();  # and round to it as accuracy
        $params[1] = undef;               # disable P
        $scale = $params[0]+4;            # at least four more for proper round
        $params[2] = $r[2];               # round mode by caller or undef
        $fallback = 1;                    # to clear a/p afterwards
    } else {
        # the 4 below is empirical, and there might be cases where it is not
        # enough...
        $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
    }

    # 1 or -1 => PI/4
    # inlined is_one() && is_one('-')
    if ($MBI->_is_one($self->{_m}) && $MBI->_is_zero($self->{_e})) {
        my $pi = $class->bpi($scale - 3);
        # modify $self in place
        $self->{_m} = $pi->{_m};
        $self->{_e} = $pi->{_e};
        $self->{_es} = $pi->{_es};
        # leave the sign of $self alone (+1 => +PI/4, -1 => -PI/4)
        $MBI->_div($self->{_m}, $MBI->_new(4));
        return $self;
    }

    # This series is only valid if -1 < x < 1, so for other x we need to
    # calculate PI/2 - atan(1/x):
    my $one = $MBI->_new(1);
    my $pi = undef;
    if ($self->bacmp($self->copy->bone) >= 0) {
        # calculate PI/2
        $pi = $class->bpi($scale - 3);
        $MBI->_div($pi->{_m}, $MBI->_new(2));
        # calculate 1/$self:
        my $self_copy = $self->copy();
        # modify $self in place
        $self->bone(); $self->bdiv($self_copy, $scale);
    }

    my $fmul = 1;
    foreach my $k (0 .. int($scale / 20)) {
        $fmul *= 2;
        $self->bdiv($self->copy->bmul($self)->binc->bsqrt($scale + 4)->binc, $scale + 4);
    }

    # When user set globals, they would interfere with our calculation, so
    # disable them and later re-enable them.
    no strict 'refs';
    my $abr = "$class\::accuracy";  my $ab = $$abr; $$abr = undef;
    my $pbr = "$class\::precision"; my $pb = $$pbr; $$pbr = undef;
    # We also need to disable any set A or P on $self (_find_round_parameters
    # took them already into account), since these would interfere, too
    delete $self->{_a}; delete $self->{_p};
    # Need to disable $upgrade in BigInt, to avoid deep recursion.
    local $Math::BigInt::upgrade = undef;

    my $last = 0;
    my $over = $self * $self;		# X ^ 2
    my $self2 = $over->copy();		# X ^ 2; difference between terms
    $over->bmul($self);			# X ^ 3 as starting value
    my $sign = 1;			# start with -=
    my $below = $class->new(3);
    my $two = $class->new(2);
    delete $self->{_a}; delete $self->{_p};

    my $limit = $class->new("1E-". ($scale-1));
    #my $steps = 0;
    while (1) {
        # We calculate the next term, and add it to the last. When the next
        # term is below our limit, it won't affect the outcome anymore, so we
        # stop:
        my $next = $over->copy()->bdiv($below, $scale);
        last if $next->bacmp($limit) <= 0;

        if ($sign == 0) {
            $self->badd($next);
        } else {
            $self->bsub($next);
        }
        $sign = 1-$sign;        # alternatex
        # calculate things for the next term
        $over->bmul($self2);    # $self*$self
        $below->badd($two);     # n += 2
    }
    $self->bmul($fmul);

    if (defined $pi) {
        my $self_copy = $self->copy();
        # modify $self in place
        $self->{_m} = $pi->{_m};
        $self->{_e} = $pi->{_e};
        $self->{_es} = $pi->{_es};
        # PI/2 - $self
        $self->bsub($self_copy);
    }

    # Shortcut to not run through _find_round_parameters again.
    if (defined $params[0]) {
        $self->bround($params[0], $params[2]); # then round accordingly
    } else {
        $self->bfround($params[1], $params[2]); # then round accordingly
    }
    if ($fallback) {
        # Clear a/p after round, since user did not request it.
        delete $self->{_a}; delete $self->{_p};
    }

    # restore globals
    $$abr = $ab; $$pbr = $pb;
    $self;
}

###############################################################################
# rounding functions

sub bfround
  {
  # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
  # $n == 0 means round to integer
  # expects and returns normalized numbers!
  my $x = shift; my $class = ref($x) || $x; $x = $class->new(shift) if !ref($x);

  my ($scale,$mode) = $x->_scale_p(@_);
  return $x if !defined $scale || $x->modify('bfround'); # no-op

  # never round a 0, +-inf, NaN
  if ($x->is_zero())
    {
    $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
    return $x; 
    }
  return $x if $x->{sign} !~ /^[+-]$/;

  # don't round if x already has lower precision
  return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});

  $x->{_p} = $scale;			# remember round in any case
  delete $x->{_a};			# and clear A
  if ($scale < 0)
    {
    # round right from the '.'

    return $x if $x->{_es} eq '+';		# e >= 0 => nothing to round

    $scale = -$scale;				# positive for simplicity
    my $len = $MBI->_len($x->{_m});		# length of mantissa

    # the following poses a restriction on _e, but if _e is bigger than a
    # scalar, you got other problems (memory etc) anyway
    my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e})));	# digits after dot
    my $zad = 0;				# zeros after dot
    $zad = $dad - $len if (-$dad < -$len);	# for 0.00..00xxx style
   
    # print "scale $scale dad $dad zad $zad len $len\n";
    # number  bsstr   len zad dad	
    # 0.123   123e-3	3   0 3
    # 0.0123  123e-4	3   1 4
    # 0.001   1e-3      1   2 3
    # 1.23    123e-2	3   0 2
    # 1.2345  12345e-4	5   0 4

    # do not round after/right of the $dad
    return $x if $scale > $dad;			# 0.123, scale >= 3 => exit

    # round to zero if rounding inside the $zad, but not for last zero like:
    # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
    return $x->bzero() if $scale < $zad;
    if ($scale == $zad)			# for 0.006, scale -3 and trunc
      {
      $scale = -$len;
      }
    else
      {
      # adjust round-point to be inside mantissa
      if ($zad != 0)
        {
	$scale = $scale-$zad;
        }
      else
        {
        my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;	# digits before dot
	$scale = $dbd+$scale;
        }
      }
    }
  else
    {
    # round left from the '.'

    # 123 => 100 means length(123) = 3 - $scale (2) => 1

    my $dbt = $MBI->_len($x->{_m}); 
    # digits before dot 
    my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
    # should be the same, so treat it as this 
    $scale = 1 if $scale == 0; 
    # shortcut if already integer 
    return $x if $scale == 1 && $dbt <= $dbd; 
    # maximum digits before dot 
    ++$dbd;

    if ($scale > $dbd) 
       { 
       # not enough digits before dot, so round to zero 
       return $x->bzero; 
       }
    elsif ( $scale == $dbd )
       { 
       # maximum 
       $scale = -$dbt; 
       } 
    else
       { 
       $scale = $dbd - $scale; 
       }
    }
  # pass sign to bround for rounding modes '+inf' and '-inf'
  my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
  $m->bround($scale,$mode);
  $x->{_m} = $m->{value};			# get our mantissa back
  $x->bnorm();
  }

sub bround
  {
  # accuracy: preserve $N digits, and overwrite the rest with 0's
  my $x = shift; my $class = ref($x) || $x; $x = $class->new(shift) if !ref($x);

  if (($_[0] || 0) < 0)
    {
    Carp::croak('bround() needs positive accuracy');
    }

  my ($scale,$mode) = $x->_scale_a(@_);
  return $x if !defined $scale || $x->modify('bround');	# no-op

  # scale is now either $x->{_a}, $accuracy, or the user parameter
  # test whether $x already has lower accuracy, do nothing in this case 
  # but do round if the accuracy is the same, since a math operation might
  # want to round a number with A=5 to 5 digits afterwards again
  return $x if defined $x->{_a} && $x->{_a} < $scale;

  # scale < 0 makes no sense
  # scale == 0 => keep all digits
  # never round a +-inf, NaN
  return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;

  # 1: never round a 0
  # 2: if we should keep more digits than the mantissa has, do nothing
  if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
    {
    $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
    return $x; 
    }

  # pass sign to bround for '+inf' and '-inf' rounding modes
  my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';

  $m->bround($scale,$mode);		# round mantissa
  $x->{_m} = $m->{value};		# get our mantissa back
  $x->{_a} = $scale;			# remember rounding
  delete $x->{_p};			# and clear P
  $x->bnorm();				# del trailing zeros gen. by bround()
  }

sub bfloor
  {
  # round towards minus infinity
  my ($class,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);

  return $x if $x->modify('bfloor');
   
  return $x if $x->{sign} !~ /^[+-]$/;	# nan, +inf, -inf

  # if $x has digits after dot
  if ($x->{_es} eq '-')
    {
    $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
    $x->{_e} = $MBI->_zero();			# trunc/norm	
    $x->{_es} = '+';				# abs e
    $MBI->_inc($x->{_m}) if $x->{sign} eq '-';	# increment if negative
    }
  $x->round($a,$p,$r);
  }

sub bceil
  {
  # round towards plus infinity
  my ($class,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);

  return $x if $x->modify('bceil');
  return $x if $x->{sign} !~ /^[+-]$/;	# nan, +inf, -inf

  # if $x has digits after dot
  if ($x->{_es} eq '-')
    {
    $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
    $x->{_e} = $MBI->_zero();			# trunc/norm	
    $x->{_es} = '+';				# abs e
    if ($x->{sign} eq '+') {
        $MBI->_inc($x->{_m});                   # increment if positive
    } else {
        $x->{sign} = '+' if $MBI->_is_zero($x->{_m});   # avoid -0
    }
    }
  $x->round($a,$p,$r);
  }

sub bint
  {
  # round towards zero
  my ($class,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);

  return $x if $x->modify('bint');
  return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf

  # if $x has digits after the decimal point
  if ($x->{_es} eq '-')
    {
    $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
    $x->{_e} = $MBI->_zero();                     # truncate/normalize
    $x->{_es} = '+';                              # abs e
    $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # avoid -0
    }
  $x->round($a,$p,$r);
  }

sub brsft {
    # shift right by $y (divide $b ** $y)

    # set up parameters
    my ($class, $x, $y, $b, $a, $p, $r) = (ref($_[0]), @_);

    # objectify is costly, so avoid it
    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) {
        ($class, $x, $y, $b, $a, $p, $r) = objectify(2, @_);
    }

    return $x if $x -> modify('brsft');
    return $x if $x -> {sign} !~ /^[+-]$/; # nan, +inf, -inf

    $b = 2 if !defined $b;
    $b = $class -> new($b) unless ref($b) && $b -> isa($class);

    return $x -> bnan() if $x -> is_nan() || $y -> is_nan() || $b -> is_nan();

    # shift by a negative amount?
    return $x -> blsft($y -> copy() -> babs(), $b) if $y -> {sign} =~ /^-/;

    # the following call to bdiv() will return either quotient (scalar context)
    # or quotient and remainder (list context).
    $x -> bdiv($b -> bpow($y), $a, $p, $r, $y);
}

sub blsft {
    # shift left by $y (multiply by $b ** $y)

    # set up parameters
    my ($class, $x, $y, $b, $a, $p, $r) = (ref($_[0]), @_);

    # objectify is costly, so avoid it
    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) {
        ($class, $x, $y, $b, $a, $p, $r) = objectify(2, @_);
    }

    return $x if $x -> modify('blsft');
    return $x if $x -> {sign} !~ /^[+-]$/; # nan, +inf, -inf

    $b = 2 if !defined $b;
    $b = $class -> new($b) unless ref($b) && $b -> isa($class);

    return $x -> bnan() if $x -> is_nan() || $y -> is_nan() || $b -> is_nan();

    # shift by a negative amount?
    return $x -> brsft($y -> copy() -> babs(), $b) if $y -> {sign} =~ /^-/;

    $x -> bmul($b -> bpow($y), $a, $p, $r, $y);
}

sub band {
    my $x     = shift;
    my $xref  = ref($x);
    my $class = $xref || $x;

    Carp::croak 'band() is an instance method, not a class method' unless $xref;
    Carp::croak 'Not enough arguments for band()' if @_ < 1;

    return if $x -> modify('band');

    my $y = shift;
    $y = $class -> new($y) unless ref($y);

    my @r = @_;

    my $xtmp = Math::BigInt -> new($x -> bint());   # to Math::BigInt
    $xtmp -> band($y);
    $xtmp = $class -> new($xtmp);                   # back to Math::BigFloat

    $x -> {sign} = $xtmp -> {sign};
    $x -> {_m}   = $xtmp -> {_m};
    $x -> {_es}  = $xtmp -> {_es};
    $x -> {_e}   = $xtmp -> {_e};

    return $x -> round(@r);
}

sub bior {
    my $x     = shift;
    my $xref  = ref($x);
    my $class = $xref || $x;

    Carp::croak 'bior() is an instance method, not a class method' unless $xref;
    Carp::croak 'Not enough arguments for bior()' if @_ < 1;

    return if $x -> modify('bior');

    my $y = shift;
    $y = $class -> new($y) unless ref($y);

    my @r = @_;

    my $xtmp = Math::BigInt -> new($x -> bint());   # to Math::BigInt
    $xtmp -> bior($y);
    $xtmp = $class -> new($xtmp);                   # back to Math::BigFloat

    $x -> {sign} = $xtmp -> {sign};
    $x -> {_m}   = $xtmp -> {_m};
    $x -> {_es}  = $xtmp -> {_es};
    $x -> {_e}   = $xtmp -> {_e};

    return $x -> round(@r);
}

sub bxor {
    my $x     = shift;
    my $xref  = ref($x);
    my $class = $xref || $x;

    Carp::croak 'bxor() is an instance method, not a class method' unless $xref;
    Carp::croak 'Not enough arguments for bxor()' if @_ < 1;

    return if $x -> modify('bxor');

    my $y = shift;
    $y = $class -> new($y) unless ref($y);

    my @r = @_;

    my $xtmp = Math::BigInt -> new($x -> bint());   # to Math::BigInt
    $xtmp -> bxor($y);
    $xtmp = $class -> new($xtmp);                   # back to Math::BigFloat

    $x -> {sign} = $xtmp -> {sign};
    $x -> {_m}   = $xtmp -> {_m};
    $x -> {_es}  = $xtmp -> {_es};
    $x -> {_e}   = $xtmp -> {_e};

    return $x -> round(@r);
}

sub bnot {
    my $x     = shift;
    my $xref  = ref($x);
    my $class = $xref || $x;

    Carp::croak 'bnot() is an instance method, not a class method' unless $xref;

    return if $x -> modify('bnot');

    my @r = @_;

    my $xtmp = Math::BigInt -> new($x -> bint());   # to Math::BigInt
    $xtmp -> bnot();
    $xtmp = $class -> new($xtmp);                   # back to Math::BigFloat

    $x -> {sign} = $xtmp -> {sign};
    $x -> {_m}   = $xtmp -> {_m};
    $x -> {_es}  = $xtmp -> {_es};
    $x -> {_e}   = $xtmp -> {_e};

    return $x -> round(@r);
}

###############################################################################

sub DESTROY
  {
  # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
  }

sub AUTOLOAD
  {
  # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
  # or falling back to MBI::bxxx()
  my $name = $AUTOLOAD;

  $name =~ s/(.*):://;	# split package
  my $c = $1 || $class;
  no strict 'refs';
  $c->import() if $IMPORT == 0;
  if (!_method_alias($name))
    {
    if (!defined $name)
      {
      # delayed load of Carp and avoid recursion	
      Carp::croak("$c: Can't call a method without name");
      }
    if (!_method_hand_up($name))
      {
      # delayed load of Carp and avoid recursion	
      Carp::croak("Can't call $c\-\>$name, not a valid method");
      }
    # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
    $name =~ s/^f/b/;
    return &{"Math::BigInt"."::$name"}(@_);
    }
  my $bname = $name; $bname =~ s/^f/b/;
  $c .= "::$name";
  *{$c} = \&{$bname};
  &{$c};	# uses @_
  }

sub exponent
  {
  # return a copy of the exponent
  my ($class,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);

  if ($x->{sign} !~ /^[+-]$/)
    {
    my $s = $x->{sign}; $s =~ s/^[+-]//;
    return Math::BigInt->new($s, undef, undef); 		# -inf, +inf => +inf
    }
  Math::BigInt->new($x->{_es} . $MBI->_str($x->{_e}), undef, undef);
  }

sub mantissa
  {
  # return a copy of the mantissa
  my ($class,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 
  if ($x->{sign} !~ /^[+-]$/)
    {
    my $s = $x->{sign}; $s =~ s/^[+]//;
    return Math::BigInt->new($s, undef, undef);	# -inf, +inf => +inf
    }
  my $m = Math::BigInt->new($MBI->_str($x->{_m}), undef, undef);
  $m->bneg() if $x->{sign} eq '-';

  $m;
  }

sub parts
  {
  # return a copy of both the exponent and the mantissa
  my ($class,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);

  if ($x->{sign} !~ /^[+-]$/)
    {
    my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
    return ($class->new($s),$class->new($se)); # +inf => inf and -inf,+inf => inf
    }
  my $m = Math::BigInt->bzero();
  $m->{value} = $MBI->_copy($x->{_m});
  $m->bneg() if $x->{sign} eq '-';
  ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
  }

##############################################################################
# private stuff (internal use only)

sub import
  {
  my $class = shift;
  my $l = scalar @_;
  my $lib = ''; my @a;
  my $lib_kind = 'try';
  $IMPORT=1;
  for ( my $i = 0; $i < $l ; $i++)
    {
    if ( $_[$i] eq ':constant' )
      {
      # This causes overlord er load to step in. 'binary' and 'integer'
      # are handled by BigInt.
      overload::constant float => sub { $class->new(shift); }; 
      }
    elsif ($_[$i] eq 'upgrade')
      {
      # this causes upgrading
      $upgrade = $_[$i+1];		# or undef to disable
      $i++;
      }
    elsif ($_[$i] eq 'downgrade')
      {
      # this causes downgrading
      $downgrade = $_[$i+1];		# or undef to disable
      $i++;
      }
    elsif ($_[$i] =~ /^(lib|try|only)\z/)
      {
      # alternative library
      $lib = $_[$i+1] || '';		# default Calc
      $lib_kind = $1;			# lib, try or only
      $i++;
      }
    elsif ($_[$i] eq 'with')
      {
      # alternative class for our private parts()
      # XXX: no longer supported
      # $MBI = $_[$i+1] || 'Math::BigInt';
      $i++;
      }
    else
      {
      push @a, $_[$i];
      }
    }

  $lib =~ tr/a-zA-Z0-9,://cd;		# restrict to sane characters
  # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
  my $mbilib = eval { Math::BigInt->config()->{lib} };
  if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
    {
    # MBI already loaded
    Math::BigInt->import( $lib_kind, "$lib,$mbilib", 'objectify');
    }
  else
    {
    # MBI not loaded, or with ne "Math::BigInt::Calc"
    $lib .= ",$mbilib" if defined $mbilib;
    $lib =~ s/^,//;				# don't leave empty 
    
    # replacement library can handle lib statement, but also could ignore it
    
    # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
    # used in the same script, or eval inside import(). So we require MBI:
    require Math::BigInt;
    Math::BigInt->import( $lib_kind => $lib, 'objectify' );
    }
  if ($@)
    {
    Carp::croak("Couldn't load $lib: $! $@");
    }
  # find out which one was actually loaded
  $MBI = Math::BigInt->config()->{lib};

  # register us with MBI to get notified of future lib changes
  Math::BigInt::_register_callback( $class, sub { $MBI = $_[0]; } );

  $class->export_to_level(1,$class,@a);		# export wanted functions
  }

sub bnorm
  {
  # adjust m and e so that m is smallest possible
  my ($class,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);

  return $x if $x->{sign} !~ /^[+-]$/;		# inf, nan etc

  my $zeros = $MBI->_zeros($x->{_m});		# correct for trailing zeros
  if ($zeros != 0)
    {
    my $z = $MBI->_new($zeros);
    $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
    if ($x->{_es} eq '-')
      {
      if ($MBI->_acmp($x->{_e},$z) >= 0)
        {
        $x->{_e} = $MBI->_sub ($x->{_e}, $z);
        $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
        }
      else
        {
        $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
        $x->{_es} = '+';
        }
      }
    else
      {
      $x->{_e} = $MBI->_add ($x->{_e}, $z);
      }
    }
  else
    {
    # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
    # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
    $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
     if $MBI->_is_zero($x->{_m});
    }

  $x;					# MBI bnorm is no-op, so do not call it
  } 
 
##############################################################################

sub as_hex
  {
  # return number as hexadecimal string (only for integers defined)
  my ($class,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);

  return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
  return '0x0' if $x->is_zero();

  return $nan if $x->{_es} ne '+';		# how to do 1e-1 in hex!?

  my $z = $MBI->_copy($x->{_m});
  if (! $MBI->_is_zero($x->{_e}))		# > 0 
    {
    $MBI->_lsft( $z, $x->{_e},10);
    }
  $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
  $z->as_hex();
  }

sub as_bin
  {
  # return number as binary digit string (only for integers defined)
  my ($class,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);

  return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
  return '0b0' if $x->is_zero();

  return $nan if $x->{_es} ne '+';		# how to do 1e-1 in hex!?

  my $z = $MBI->_copy($x->{_m});
  if (! $MBI->_is_zero($x->{_e}))		# > 0 
    {
    $MBI->_lsft( $z, $x->{_e},10);
    }
  $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
  $z->as_bin();
  }

sub as_oct
  {
  # return number as octal digit string (only for integers defined)
  my ($class,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);

  return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
  return '0' if $x->is_zero();

  return $nan if $x->{_es} ne '+';		# how to do 1e-1 in hex!?

  my $z = $MBI->_copy($x->{_m});
  if (! $MBI->_is_zero($x->{_e}))		# > 0 
    {
    $MBI->_lsft( $z, $x->{_e},10);
    }
  $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
  $z->as_oct();
  }

sub as_number
  {
  # return copy as a bigint representation of this BigFloat number
  my ($class,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);

  return $x if $x->modify('as_number');

  if (!$x->isa('Math::BigFloat'))
    {
    # if the object can as_number(), use it
    return $x->as_number() if $x->can('as_number');
    # otherwise, get us a float and then a number
    $x = $x->can('as_float') ? $x->as_float() : $class->new(0+"$x");
    }

  return Math::BigInt->binf($x->sign()) if $x->is_inf();
  return Math::BigInt->bnan()           if $x->is_nan();

  my $z = $MBI->_copy($x->{_m});
  if ($x->{_es} eq '-')			# < 0
    {
    $MBI->_rsft( $z, $x->{_e},10);
    } 
  elsif (! $MBI->_is_zero($x->{_e}))	# > 0 
    {
    $MBI->_lsft( $z, $x->{_e},10);
    }
  $z = Math::BigInt->new( $x->{sign} . $MBI->_str($z));
  $z;
  }

sub length
  {
  my $x = shift;
  my $class = ref($x) || $x;
  $x = $class->new(shift) unless ref($x);

  return 1 if $MBI->_is_zero($x->{_m});

  my $len = $MBI->_len($x->{_m});
  $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
  if (wantarray())
    {
    my $t = 0;
    $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
    return ($len, $t);
    }
  $len;
  }

sub from_hex {
    my $self    = shift;
    my $selfref = ref $self;
    my $class   = $selfref || $self;

    my $str = shift;

    # If called as a class method, initialize a new object.

    $self = $class -> bzero() unless $selfref;

    if ($str =~ s/
                     ^

                     # sign
                     ( [+-]? )

                     # optional "hex marker"
                     (?: 0? x )?

                     # significand using the hex digits 0..9 and a..f
                     (
                         [0-9a-fA-F]+ (?: _ [0-9a-fA-F]+ )*
                         (?:
                             \.
                             (?: [0-9a-fA-F]+ (?: _ [0-9a-fA-F]+ )* )?
                         )?
                     |
                         \.
                         [0-9a-fA-F]+ (?: _ [0-9a-fA-F]+ )*
                     )

                     # exponent (power of 2) using decimal digits
                     (?:
                         [Pp]
                         ( [+-]? )
                         ( \d+ (?: _ \d+ )* )
                     )?

                     $
                 //x)
    {
        my $s_sign  = $1 || '+';
        my $s_value = $2;
        my $e_sign  = $3 || '+';
        my $e_value = $4 || '0';
        $s_value =~ tr/_//d;
        $e_value =~ tr/_//d;

        # The significand must be multiplied by 2 raised to this exponent.

        my $two_expon = $class -> new($e_value);
        $two_expon -> bneg() if $e_sign eq '-';

        # If there is a dot in the significand, remove it and adjust the
        # exponent according to the number of digits in the fraction part of
        # the significand. Since the digits in the significand are in base 16,
        # but the exponent is only in base 2, multiply the exponent adjustment
        # value by log(16) / log(2) = 4.

        my $idx = index($s_value, '.');
        if ($idx >= 0) {
            substr($s_value, $idx, 1) = '';
            $two_expon -= $class -> new(CORE::length($s_value))
                                 -> bsub($idx)
                                 -> bmul("4");
        }

        $self -> {sign} = $s_sign;
        $self -> {_m}   = $MBI -> _from_hex('0x' . $s_value);

        if ($two_expon > 0) {
            my $factor = $class -> new("2") -> bpow($two_expon);
            $self -> bmul($factor);
        } elsif ($two_expon < 0) {
            my $factor = $class -> new("0.5") -> bpow(-$two_expon);
            $self -> bmul($factor);
        }

        return $self;
    }

    return $self->bnan();
}

sub from_oct {
    my $self    = shift;
    my $selfref = ref $self;
    my $class   = $selfref || $self;

    my $str = shift;

    # If called as a class method, initialize a new object.

    $self = $class -> bzero() unless $selfref;

    if ($str =~ s/
                     ^

                     # sign
                     ( [+-]? )

                     # significand using the octal digits 0..7
                     (
                         [0-7]+ (?: _ [0-7]+ )*
                         (?:
                             \.
                             (?: [0-7]+ (?: _ [0-7]+ )* )?
                         )?
                     |
                         \.
                         [0-7]+ (?: _ [0-7]+ )*
                     )

                     # exponent (power of 2) using decimal digits
                     (?:
                         [Pp]
                         ( [+-]? )
                         ( \d+ (?: _ \d+ )* )
                     )?

                     $
                 //x)
    {
        my $s_sign  = $1 || '+';
        my $s_value = $2;
        my $e_sign  = $3 || '+';
        my $e_value = $4 || '0';
        $s_value =~ tr/_//d;
        $e_value =~ tr/_//d;

        # The significand must be multiplied by 2 raised to this exponent.

        my $two_expon = $class -> new($e_value);
        $two_expon -> bneg() if $e_sign eq '-';

        # If there is a dot in the significand, remove it and adjust the
        # exponent according to the number of digits in the fraction part of
        # the significand. Since the digits in the significand are in base 8,
        # but the exponent is only in base 2, multiply the exponent adjustment
        # value by log(8) / log(2) = 3.

        my $idx = index($s_value, '.');
        if ($idx >= 0) {
            substr($s_value, $idx, 1) = '';
            $two_expon -= $class -> new(CORE::length($s_value))
                                 -> bsub($idx)
                                 -> bmul("3");
        }

        $self -> {sign} = $s_sign;
        $self -> {_m}   = $MBI -> _from_oct($s_value);

        if ($two_expon > 0) {
            my $factor = $class -> new("2") -> bpow($two_expon);
            $self -> bmul($factor);
        } elsif ($two_expon < 0) {
            my $factor = $class -> new("0.5") -> bpow(-$two_expon);
            $self -> bmul($factor);
        }

        return $self;
    }

    return $self->bnan();
}

sub from_bin {
    my $self    = shift;
    my $selfref = ref $self;
    my $class   = $selfref || $self;

    my $str = shift;

    # If called as a class method, initialize a new object.

    $self = $class -> bzero() unless $selfref;

    if ($str =~ s/
                     ^

                     # sign
                     ( [+-]? )

                     # optional "bin marker"
                     (?: 0? b )?

                     # significand using the binary digits 0 and 1
                     (
                         [01]+ (?: _ [01]+ )*
                         (?:
                             \.
                             (?: [01]+ (?: _ [01]+ )* )?
                         )?
                     |
                         \.
                         [01]+ (?: _ [01]+ )*
                     )

                     # exponent (power of 2) using decimal digits
                     (?:
                         [Pp]
                         ( [+-]? )
                         ( \d+ (?: _ \d+ )* )
                     )?

                     $
                 //x)
    {
        my $s_sign  = $1 || '+';
        my $s_value = $2;
        my $e_sign  = $3 || '+';
        my $e_value = $4 || '0';
        $s_value =~ tr/_//d;
        $e_value =~ tr/_//d;

        # The significand must be multiplied by 2 raised to this exponent.

        my $two_expon = $class -> new($e_value);
        $two_expon -> bneg() if $e_sign eq '-';

        # If there is a dot in the significand, remove it and adjust the
        # exponent according to the number of digits in the fraction part of
        # the significand.

        my $idx = index($s_value, '.');
        if ($idx >= 0) {
            substr($s_value, $idx, 1) = '';
            $two_expon -= $class -> new(CORE::length($s_value))
                                 -> bsub($idx);
        }

        $self -> {sign} = $s_sign;
        $self -> {_m}   = $MBI -> _from_bin('0b' . $s_value);

        if ($two_expon > 0) {
            my $factor = $class -> new("2") -> bpow($two_expon);
            $self -> bmul($factor);
        } elsif ($two_expon < 0) {
            my $factor = $class -> new("0.5") -> bpow(-$two_expon);
            $self -> bmul($factor);
        }

        return $self;
    }

    return $self->bnan();
}

1;

__END__

#line 5557
