home *** CD-ROM | disk | FTP | other *** search
/ PC World 2005 December (Special) / PCWorld_2005-12_Special_cd.bin / Bezpecnost / lsti / lsti.exe / framework-2.5.exe / BigFloat.pm < prev    next >
Text File  |  2005-01-27  |  90KB  |  3,000 lines

  1. package Math::BigFloat;
  2.  
  3. # Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
  4. #
  5.  
  6. # The following hash values are internally used:
  7. #   _e    : exponent (ref to $CALC object)
  8. #   _m    : mantissa (ref to $CALC object)
  9. #   _es    : sign of _e
  10. # sign    : +,-,+inf,-inf, or "NaN" if not a number
  11. #   _a    : accuracy
  12. #   _p    : precision
  13.  
  14. $VERSION = '1.47';
  15. require 5.005;
  16.  
  17. require Exporter;
  18. @ISA =       qw(Exporter Math::BigInt);
  19.  
  20. use strict;
  21. # $_trap_inf/$_trap_nan are internal and should never be accessed from outside
  22. use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode
  23.         $upgrade $downgrade $_trap_nan $_trap_inf/;
  24. my $class = "Math::BigFloat";
  25.  
  26. use overload
  27. '<=>'    =>    sub { $_[2] ?
  28.                       ref($_[0])->bcmp($_[1],$_[0]) : 
  29.                       ref($_[0])->bcmp($_[0],$_[1])},
  30. 'int'    =>    sub { $_[0]->as_number() },        # 'trunc' to bigint
  31. ;
  32.  
  33. ##############################################################################
  34. # global constants, flags and assorted stuff
  35.  
  36. # the following are public, but their usage is not recommended. Use the
  37. # accessor methods instead.
  38.  
  39. # class constants, use Class->constant_name() to access
  40. $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
  41. $accuracy   = undef;
  42. $precision  = undef;
  43. $div_scale  = 40;
  44.  
  45. $upgrade = undef;
  46. $downgrade = undef;
  47. # the package we are using for our private parts, defaults to:
  48. # Math::BigInt->config()->{lib}
  49. my $MBI = 'Math::BigInt::Calc';
  50.  
  51. # are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
  52. $_trap_nan = 0;
  53. # the same for infinity
  54. $_trap_inf = 0;
  55.  
  56. # constant for easier life
  57. my $nan = 'NaN'; 
  58.  
  59. my $IMPORT = 0;    # was import() called yet? used to make require work
  60.  
  61. # some digits of accuracy for blog(undef,10); which we use in blog() for speed
  62. my $LOG_10 = 
  63.  '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
  64. my $LOG_10_A = length($LOG_10)-1;
  65. # ditto for log(2)
  66. my $LOG_2 = 
  67.  '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
  68. my $LOG_2_A = length($LOG_2)-1;
  69. my $HALF = '0.5';            # made into an object if necc.
  70.  
  71. ##############################################################################
  72. # the old code had $rnd_mode, so we need to support it, too
  73.  
  74. sub TIESCALAR   { my ($class) = @_; bless \$round_mode, $class; }
  75. sub FETCH       { return $round_mode; }
  76. sub STORE       { $rnd_mode = $_[0]->round_mode($_[1]); }
  77.  
  78. BEGIN
  79.   {
  80.   # when someone set's $rnd_mode, we catch this and check the value to see
  81.   # whether it is valid or not. 
  82.   $rnd_mode   = 'even'; tie $rnd_mode, 'Math::BigFloat'; 
  83.   }
  84.  
  85. ##############################################################################
  86.  
  87. {
  88.   # valid method aliases for AUTOLOAD
  89.   my %methods = map { $_ => 1 }  
  90.    qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
  91.         fint facmp fcmp fzero fnan finf finc fdec flog ffac
  92.     fceil ffloor frsft flsft fone flog froot
  93.       /;
  94.   # valid method's that can be hand-ed up (for AUTOLOAD)
  95.   my %hand_ups = map { $_ => 1 }  
  96.    qw / is_nan is_inf is_negative is_positive is_pos is_neg
  97.         accuracy precision div_scale round_mode fneg fabs fnot
  98.         objectify upgrade downgrade
  99.     bone binf bnan bzero
  100.       /;
  101.  
  102.   sub method_alias { exists $methods{$_[0]||''}; } 
  103.   sub method_hand_up { exists $hand_ups{$_[0]||''}; } 
  104. }
  105.  
  106. ##############################################################################
  107. # constructors
  108.  
  109. sub new 
  110.   {
  111.   # create a new BigFloat object from a string or another bigfloat object. 
  112.   # _e: exponent
  113.   # _m: mantissa
  114.   # sign  => sign (+/-), or "NaN"
  115.  
  116.   my ($class,$wanted,@r) = @_;
  117.  
  118.   # avoid numify-calls by not using || on $wanted!
  119.   return $class->bzero() if !defined $wanted;    # default to 0
  120.   return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
  121.  
  122.   $class->import() if $IMPORT == 0;             # make require work
  123.  
  124.   my $self = {}; bless $self, $class;
  125.   # shortcut for bigints and its subclasses
  126.   if ((ref($wanted)) && (ref($wanted) ne $class))
  127.     {
  128.     $self->{_m} = $wanted->as_number()->{value}; # get us a bigint copy
  129.     $self->{_e} = $MBI->_zero();
  130.     $self->{_es} = '+';
  131.     $self->{sign} = $wanted->sign();
  132.     return $self->bnorm();
  133.     }
  134.   # else: got a string
  135.  
  136.   # handle '+inf', '-inf' first
  137.   if ($wanted =~ /^[+-]?inf$/)
  138.     {
  139.     return $downgrade->new($wanted) if $downgrade;
  140.  
  141.     $self->{_e} = $MBI->_zero();
  142.     $self->{_es} = '+';
  143.     $self->{_m} = $MBI->_zero();
  144.     $self->{sign} = $wanted;
  145.     $self->{sign} = '+inf' if $self->{sign} eq 'inf';
  146.     return $self->bnorm();
  147.     }
  148.  
  149.   # shortcut for simple forms like '12' that neither have trailing nor leading
  150.   # zeros
  151.   if ($wanted =~ /^([+-]?)([1-9][0-9]*[1-9])$/)
  152.     {
  153.     $self->{_e} = $MBI->_zero();
  154.     $self->{_es} = '+';
  155.     $self->{sign} = $1 || '+';
  156.     $self->{_m} = $MBI->_new($2);
  157.     return $self->round(@r) if !$downgrade;
  158.     }
  159.  
  160.   my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
  161.   if (!ref $mis)
  162.     {
  163.     if ($_trap_nan)
  164.       {
  165.       require Carp;
  166.       Carp::croak ("$wanted is not a number initialized to $class");
  167.       }
  168.     
  169.     return $downgrade->bnan() if $downgrade;
  170.     
  171.     $self->{_e} = $MBI->_zero();
  172.     $self->{_es} = '+';
  173.     $self->{_m} = $MBI->_zero();
  174.     $self->{sign} = $nan;
  175.     }
  176.   else
  177.     {
  178.     # make integer from mantissa by adjusting exp, then convert to int
  179.     $self->{_e} = $MBI->_new($$ev);        # exponent
  180.     $self->{_es} = $$es || '+';
  181.     my $mantissa = "$$miv$$mfv";         # create mant.
  182.     $mantissa =~ s/^0+(\d)/$1/;            # strip leading zeros
  183.     $self->{_m} = $MBI->_new($mantissa);     # create mant.
  184.  
  185.     # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
  186.     if (CORE::length($$mfv) != 0)
  187.       {
  188.       my $len = $MBI->_new( CORE::length($$mfv));
  189.       ($self->{_e}, $self->{_es}) =
  190.     _e_sub ($self->{_e}, $len, $self->{_es}, '+');
  191.       }
  192.     # we can only have trailing zeros on the mantissa if $$mfv eq ''
  193.     else
  194.       {
  195.       # Use a regexp to count the trailing zeros in $$miv instead of _zeros()
  196.       # because that is faster, especially when _m is not stored in base 10.
  197.       my $zeros = 0; $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/; 
  198.       if ($zeros != 0)
  199.         {
  200.         my $z = $MBI->_new($zeros);
  201.         # turn '120e2' into '12e3'
  202.         $MBI->_rsft ( $self->{_m}, $z, 10);
  203.     _e_add ( $self->{_e}, $z, $self->{_es}, '+');
  204.         }
  205.       }
  206.     $self->{sign} = $$mis;
  207.  
  208.     # for something like 0Ey, set y to 1, and -0 => +0
  209.     # Check $$miv for beeing '0' and $$mfv eq '', because otherwise _m could not
  210.     # have become 0. That's faster than to call $MBI->_is_zero().
  211.     $self->{sign} = '+', $self->{_e} = $MBI->_one()
  212.      if $$miv eq '0' and $$mfv eq '';
  213.  
  214.     return $self->round(@r) if !$downgrade;
  215.     }
  216.   # if downgrade, inf, NaN or integers go down
  217.  
  218.   if ($downgrade && $self->{_es} eq '+')
  219.     {
  220.     if ($MBI->_is_zero( $self->{_e} ))
  221.       {
  222.       return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
  223.       }
  224.     return $downgrade->new($self->bsstr()); 
  225.     }
  226.   $self->bnorm()->round(@r);            # first normalize, then round
  227.   }
  228.  
  229. sub copy
  230.   {
  231.   my ($c,$x);
  232.   if (@_ > 1)
  233.     {
  234.     # if two arguments, the first one is the class to "swallow" subclasses
  235.     ($c,$x) = @_;
  236.     }
  237.   else
  238.     {
  239.     $x = shift;
  240.     $c = ref($x);
  241.     }
  242.   return unless ref($x); # only for objects
  243.  
  244.   my $self = {}; bless $self,$c;
  245.  
  246.   $self->{sign} = $x->{sign};
  247.   $self->{_es} = $x->{_es};
  248.   $self->{_m} = $MBI->_copy($x->{_m});
  249.   $self->{_e} = $MBI->_copy($x->{_e});
  250.   $self->{_a} = $x->{_a} if defined $x->{_a};
  251.   $self->{_p} = $x->{_p} if defined $x->{_p};
  252.   $self;
  253.   }
  254.  
  255. sub _bnan
  256.   {
  257.   # used by parent class bone() to initialize number to NaN
  258.   my $self = shift;
  259.   
  260.   if ($_trap_nan)
  261.     {
  262.     require Carp;
  263.     my $class = ref($self);
  264.     Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
  265.     }
  266.  
  267.   $IMPORT=1;                    # call our import only once
  268.   $self->{_m} = $MBI->_zero();
  269.   $self->{_e} = $MBI->_zero();
  270.   $self->{_es} = '+';
  271.   }
  272.  
  273. sub _binf
  274.   {
  275.   # used by parent class bone() to initialize number to +-inf
  276.   my $self = shift;
  277.   
  278.   if ($_trap_inf)
  279.     {
  280.     require Carp;
  281.     my $class = ref($self);
  282.     Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
  283.     }
  284.  
  285.   $IMPORT=1;                    # call our import only once
  286.   $self->{_m} = $MBI->_zero();
  287.   $self->{_e} = $MBI->_zero();
  288.   $self->{_es} = '+';
  289.   }
  290.  
  291. sub _bone
  292.   {
  293.   # used by parent class bone() to initialize number to 1
  294.   my $self = shift;
  295.   $IMPORT=1;                    # call our import only once
  296.   $self->{_m} = $MBI->_one();
  297.   $self->{_e} = $MBI->_zero();
  298.   $self->{_es} = '+';
  299.   }
  300.  
  301. sub _bzero
  302.   {
  303.   # used by parent class bone() to initialize number to 0
  304.   my $self = shift;
  305.   $IMPORT=1;                    # call our import only once
  306.   $self->{_m} = $MBI->_zero();
  307.   $self->{_e} = $MBI->_one();
  308.   $self->{_es} = '+';
  309.   }
  310.  
  311. sub isa
  312.   {
  313.   my ($self,$class) = @_;
  314.   return if $class =~ /^Math::BigInt/;        # we aren't one of these
  315.   UNIVERSAL::isa($self,$class);
  316.   }
  317.  
  318. sub config
  319.   {
  320.   # return (later set?) configuration data as hash ref
  321.   my $class = shift || 'Math::BigFloat';
  322.  
  323.   my $cfg = $class->SUPER::config(@_);
  324.  
  325.   # now we need only to override the ones that are different from our parent
  326.   $cfg->{class} = $class;
  327.   $cfg->{with} = $MBI;
  328.   $cfg;
  329.   }
  330.  
  331. ##############################################################################
  332. # string conversation
  333.  
  334. sub bstr 
  335.   {
  336.   # (ref to BFLOAT or num_str ) return num_str
  337.   # Convert number from internal format to (non-scientific) string format.
  338.   # internal format is always normalized (no leading zeros, "-0" => "+0")
  339.   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
  340.  
  341.   if ($x->{sign} !~ /^[+-]$/)
  342.     {
  343.     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
  344.     return 'inf';                                       # +inf
  345.     }
  346.  
  347.   my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
  348.  
  349.   # $x is zero?
  350.   my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
  351.   if ($not_zero)
  352.     {
  353.     $es = $MBI->_str($x->{_m});
  354.     $len = CORE::length($es);
  355.     my $e = $MBI->_num($x->{_e});    
  356.     $e = -$e if $x->{_es} eq '-';
  357.     if ($e < 0)
  358.       {
  359.       $dot = '';
  360.       # if _e is bigger than a scalar, the following will blow your memory
  361.       if ($e <= -$len)
  362.         {
  363.         my $r = abs($e) - $len;
  364.         $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
  365.         }
  366.       else
  367.         {
  368.         substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
  369.         $cad = -$cad if $x->{_es} eq '-';
  370.         }
  371.       }
  372.     elsif ($e > 0)
  373.       {
  374.       # expand with zeros
  375.       $es .= '0' x $e; $len += $e; $cad = 0;
  376.       }
  377.     } # if not zero
  378.  
  379.   $es = '-'.$es if $x->{sign} eq '-';
  380.   # if set accuracy or precision, pad with zeros on the right side
  381.   if ((defined $x->{_a}) && ($not_zero))
  382.     {
  383.     # 123400 => 6, 0.1234 => 4, 0.001234 => 4
  384.     my $zeros = $x->{_a} - $cad;        # cad == 0 => 12340
  385.     $zeros = $x->{_a} - $len if $cad != $len;
  386.     $es .= $dot.'0' x $zeros if $zeros > 0;
  387.     }
  388.   elsif ((($x->{_p} || 0) < 0))
  389.     {
  390.     # 123400 => 6, 0.1234 => 4, 0.001234 => 6
  391.     my $zeros = -$x->{_p} + $cad;
  392.     $es .= $dot.'0' x $zeros if $zeros > 0;
  393.     }
  394.   $es;
  395.   }
  396.  
  397. sub bsstr
  398.   {
  399.   # (ref to BFLOAT or num_str ) return num_str
  400.   # Convert number from internal format to scientific string format.
  401.   # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
  402.   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
  403.  
  404.   if ($x->{sign} !~ /^[+-]$/)
  405.     {
  406.     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
  407.     return 'inf';                                       # +inf
  408.     }
  409.   my $sep = 'e'.$x->{_es};
  410.   my $sign = $x->{sign}; $sign = '' if $sign eq '+';
  411.   $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
  412.   }
  413.     
  414. sub numify 
  415.   {
  416.   # Make a number from a BigFloat object
  417.   # simple return a string and let Perl's atoi()/atof() handle the rest
  418.   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
  419.   $x->bsstr(); 
  420.   }
  421.  
  422. ##############################################################################
  423. # public stuff (usually prefixed with "b")
  424.  
  425. # tels 2001-08-04 
  426. # XXX TODO this must be overwritten and return NaN for non-integer values
  427. # band(), bior(), bxor(), too
  428. #sub bnot
  429. #  {
  430. #  $class->SUPER::bnot($class,@_);
  431. #  }
  432.  
  433. sub bcmp 
  434.   {
  435.   # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
  436.  
  437.   # set up parameters
  438.   my ($self,$x,$y) = (ref($_[0]),@_);
  439.   # objectify is costly, so avoid it
  440.   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
  441.     {
  442.     ($self,$x,$y) = objectify(2,@_);
  443.     }
  444.  
  445.   return $upgrade->bcmp($x,$y) if defined $upgrade &&
  446.     ((!$x->isa($self)) || (!$y->isa($self)));
  447.  
  448.   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
  449.     {
  450.     # handle +-inf and NaN
  451.     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
  452.     return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
  453.     return +1 if $x->{sign} eq '+inf';
  454.     return -1 if $x->{sign} eq '-inf';
  455.     return -1 if $y->{sign} eq '+inf';
  456.     return +1;
  457.     }
  458.  
  459.   # check sign for speed first
  460.   return 1 if $x->{sign} eq '+' && $y->{sign} eq '-';    # does also 0 <=> -y
  461.   return -1 if $x->{sign} eq '-' && $y->{sign} eq '+';    # does also -x <=> 0
  462.  
  463.   # shortcut 
  464.   my $xz = $x->is_zero();
  465.   my $yz = $y->is_zero();
  466.   return 0 if $xz && $yz;                # 0 <=> 0
  467.   return -1 if $xz && $y->{sign} eq '+';        # 0 <=> +y
  468.   return 1 if $yz && $x->{sign} eq '+';            # +x <=> 0
  469.  
  470.   # adjust so that exponents are equal
  471.   my $lxm = $MBI->_len($x->{_m});
  472.   my $lym = $MBI->_len($y->{_m});
  473.   # the numify somewhat limits our length, but makes it much faster
  474.   my ($xes,$yes) = (1,1);
  475.   $xes = -1 if $x->{_es} ne '+';
  476.   $yes = -1 if $y->{_es} ne '+';
  477.   my $lx = $lxm + $xes * $MBI->_num($x->{_e});
  478.   my $ly = $lym + $yes * $MBI->_num($y->{_e});
  479.   my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
  480.   return $l <=> 0 if $l != 0;
  481.   
  482.   # lengths (corrected by exponent) are equal
  483.   # so make mantissa equal length by padding with zero (shift left)
  484.   my $diff = $lxm - $lym;
  485.   my $xm = $x->{_m};        # not yet copy it
  486.   my $ym = $y->{_m};
  487.   if ($diff > 0)
  488.     {
  489.     $ym = $MBI->_copy($y->{_m});
  490.     $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
  491.     }
  492.   elsif ($diff < 0)
  493.     {
  494.     $xm = $MBI->_copy($x->{_m});
  495.     $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
  496.     }
  497.   my $rc = $MBI->_acmp($xm,$ym);
  498.   $rc = -$rc if $x->{sign} eq '-';        # -124 < -123
  499.   $rc <=> 0;
  500.   }
  501.  
  502. sub bacmp 
  503.   {
  504.   # Compares 2 values, ignoring their signs. 
  505.   # Returns one of undef, <0, =0, >0. (suitable for sort)
  506.   
  507.   # set up parameters
  508.   my ($self,$x,$y) = (ref($_[0]),@_);
  509.   # objectify is costly, so avoid it
  510.   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
  511.     {
  512.     ($self,$x,$y) = objectify(2,@_);
  513.     }
  514.  
  515.   return $upgrade->bacmp($x,$y) if defined $upgrade &&
  516.     ((!$x->isa($self)) || (!$y->isa($self)));
  517.  
  518.   # handle +-inf and NaN's
  519.   if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
  520.     {
  521.     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
  522.     return 0 if ($x->is_inf() && $y->is_inf());
  523.     return 1 if ($x->is_inf() && !$y->is_inf());
  524.     return -1;
  525.     }
  526.  
  527.   # shortcut 
  528.   my $xz = $x->is_zero();
  529.   my $yz = $y->is_zero();
  530.   return 0 if $xz && $yz;                # 0 <=> 0
  531.   return -1 if $xz && !$yz;                # 0 <=> +y
  532.   return 1 if $yz && !$xz;                # +x <=> 0
  533.  
  534.   # adjust so that exponents are equal
  535.   my $lxm = $MBI->_len($x->{_m});
  536.   my $lym = $MBI->_len($y->{_m});
  537.   my ($xes,$yes) = (1,1);
  538.   $xes = -1 if $x->{_es} ne '+';
  539.   $yes = -1 if $y->{_es} ne '+';
  540.   # the numify somewhat limits our length, but makes it much faster
  541.   my $lx = $lxm + $xes * $MBI->_num($x->{_e});
  542.   my $ly = $lym + $yes * $MBI->_num($y->{_e});
  543.   my $l = $lx - $ly;
  544.   return $l <=> 0 if $l != 0;
  545.   
  546.   # lengths (corrected by exponent) are equal
  547.   # so make mantissa equal-length by padding with zero (shift left)
  548.   my $diff = $lxm - $lym;
  549.   my $xm = $x->{_m};        # not yet copy it
  550.   my $ym = $y->{_m};
  551.   if ($diff > 0)
  552.     {
  553.     $ym = $MBI->_copy($y->{_m});
  554.     $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
  555.     }
  556.   elsif ($diff < 0)
  557.     {
  558.     $xm = $MBI->_copy($x->{_m});
  559.     $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
  560.     }
  561.   $MBI->_acmp($xm,$ym);
  562.   }
  563.  
  564. sub badd 
  565.   {
  566.   # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
  567.   # return result as BFLOAT
  568.  
  569.   # set up parameters
  570.   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
  571.   # objectify is costly, so avoid it
  572.   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
  573.     {
  574.     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
  575.     }
  576.  
  577.   # inf and NaN handling
  578.   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
  579.     {
  580.     # NaN first
  581.     return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
  582.     # inf handling
  583.     if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
  584.       {
  585.       # +inf++inf or -inf+-inf => same, rest is NaN
  586.       return $x if $x->{sign} eq $y->{sign};
  587.       return $x->bnan();
  588.       }
  589.     # +-inf + something => +inf; something +-inf => +-inf
  590.     $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
  591.     return $x;
  592.     }
  593.  
  594.   return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
  595.    ((!$x->isa($self)) || (!$y->isa($self)));
  596.  
  597.   # speed: no add for 0+y or x+0
  598.   return $x->bround($a,$p,$r) if $y->is_zero();        # x+0
  599.   if ($x->is_zero())                    # 0+y
  600.     {
  601.     # make copy, clobbering up x (modify in place!)
  602.     $x->{_e} = $MBI->_copy($y->{_e});
  603.     $x->{_es} = $y->{_es};
  604.     $x->{_m} = $MBI->_copy($y->{_m});
  605.     $x->{sign} = $y->{sign} || $nan;
  606.     return $x->round($a,$p,$r,$y);
  607.     }
  608.  
  609.   # take lower of the two e's and adapt m1 to it to match m2
  610.   my $e = $y->{_e};
  611.   $e = $MBI->_zero() if !defined $e;        # if no BFLOAT?
  612.   $e = $MBI->_copy($e);                # make copy (didn't do it yet)
  613.  
  614.   my $es;
  615.  
  616.   ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});
  617.  
  618.   my $add = $MBI->_copy($y->{_m});
  619.  
  620.   if ($es eq '-')                # < 0
  621.     {
  622.     $MBI->_lsft( $x->{_m}, $e, 10);
  623.     ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
  624.     }
  625.   elsif (!$MBI->_is_zero($e))            # > 0
  626.     {
  627.     $MBI->_lsft($add, $e, 10);
  628.     }
  629.   # else: both e are the same, so just leave them
  630.  
  631.   if ($x->{sign} eq $y->{sign})
  632.     {
  633.     # add
  634.     $x->{_m} = $MBI->_add($x->{_m}, $add);
  635.     }
  636.   else
  637.     {
  638.     ($x->{_m}, $x->{sign}) = 
  639.      _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
  640.     }
  641.  
  642.   # delete trailing zeros, then round
  643.   $x->bnorm()->round($a,$p,$r,$y);
  644.   }
  645.  
  646. # sub bsub is inherited from Math::BigInt!
  647.  
  648. sub binc
  649.   {
  650.   # increment arg by one
  651.   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
  652.  
  653.   if ($x->{_es} eq '-')
  654.     {
  655.     return $x->badd($self->bone(),@r);    #  digits after dot
  656.     }
  657.  
  658.   if (!$MBI->_is_zero($x->{_e}))        # _e == 0 for NaN, inf, -inf
  659.     {
  660.     # 1e2 => 100, so after the shift below _m has a '0' as last digit
  661.     $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10);    # 1e2 => 100
  662.     $x->{_e} = $MBI->_zero();                # normalize
  663.     $x->{_es} = '+';
  664.     # we know that the last digit of $x will be '1' or '9', depending on the
  665.     # sign
  666.     }
  667.   # now $x->{_e} == 0
  668.   if ($x->{sign} eq '+')
  669.     {
  670.     $MBI->_inc($x->{_m});
  671.     return $x->bnorm()->bround(@r);
  672.     }
  673.   elsif ($x->{sign} eq '-')
  674.     {
  675.     $MBI->_dec($x->{_m});
  676.     $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
  677.     return $x->bnorm()->bround(@r);
  678.     }
  679.   # inf, nan handling etc
  680.   $x->badd($self->bone(),@r);            # badd() does round 
  681.   }
  682.  
  683. sub bdec
  684.   {
  685.   # decrement arg by one
  686.   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
  687.  
  688.   if ($x->{_es} eq '-')
  689.     {
  690.     return $x->badd($self->bone('-'),@r);    #  digits after dot
  691.     }
  692.  
  693.   if (!$MBI->_is_zero($x->{_e}))
  694.     {
  695.     $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10);    # 1e2 => 100
  696.     $x->{_e} = $MBI->_zero();                # normalize
  697.     $x->{_es} = '+';
  698.     }
  699.   # now $x->{_e} == 0
  700.   my $zero = $x->is_zero();
  701.   # <= 0
  702.   if (($x->{sign} eq '-') || $zero)
  703.     {
  704.     $MBI->_inc($x->{_m});
  705.     $x->{sign} = '-' if $zero;                # 0 => 1 => -1
  706.     $x->{sign} = '+' if $MBI->_is_zero($x->{_m});    # -1 +1 => -0 => +0
  707.     return $x->bnorm()->round(@r);
  708.     }
  709.   # > 0
  710.   elsif ($x->{sign} eq '+')
  711.     {
  712.     $MBI->_dec($x->{_m});
  713.     return $x->bnorm()->round(@r);
  714.     }
  715.   # inf, nan handling etc
  716.   $x->badd($self->bone('-'),@r);        # does round
  717.   } 
  718.  
  719. sub DEBUG () { 0; }
  720.  
  721. sub blog
  722.   {
  723.   my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
  724.  
  725.   # $base > 0, $base != 1; if $base == undef default to $base == e
  726.   # $x >= 0
  727.  
  728.   # we need to limit the accuracy to protect against overflow
  729.   my $fallback = 0;
  730.   my ($scale,@params);
  731.   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
  732.  
  733.   # also takes care of the "error in _find_round_parameters?" case
  734.   return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
  735.  
  736.  
  737.   # no rounding at all, so must use fallback
  738.   if (scalar @params == 0)
  739.     {
  740.     # simulate old behaviour
  741.     $params[0] = $self->div_scale();    # and round to it as accuracy
  742.     $params[1] = undef;            # P = undef
  743.     $scale = $params[0]+4;         # at least four more for proper round
  744.     $params[2] = $r;            # round mode by caller or undef
  745.     $fallback = 1;            # to clear a/p afterwards
  746.     }
  747.   else
  748.     {
  749.     # the 4 below is empirical, and there might be cases where it is not
  750.     # enough...
  751.     $scale = abs($params[0] || $params[1]) + 4;    # take whatever is defined
  752.     }
  753.  
  754.   return $x->bzero(@params) if $x->is_one();
  755.   # base not defined => base == Euler's constant e
  756.   if (defined $base)
  757.     {
  758.     # make object, since we don't feed it through objectify() to still get the
  759.     # case of $base == undef
  760.     $base = $self->new($base) unless ref($base);
  761.     # $base > 0; $base != 1
  762.     return $x->bnan() if $base->is_zero() || $base->is_one() ||
  763.       $base->{sign} ne '+';
  764.     # if $x == $base, we know the result must be 1.0
  765.     if ($x->bcmp($base) == 0)
  766.       {
  767.       $x->bone('+',@params);
  768.       if ($fallback)
  769.         {
  770.         # clear a/p after round, since user did not request it
  771.         delete $x->{_a}; delete $x->{_p};
  772.         }
  773.       return $x;
  774.       }
  775.     }
  776.  
  777.   # when user set globals, they would interfere with our calculation, so
  778.   # disable them and later re-enable them
  779.   no strict 'refs';
  780.   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
  781.   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
  782.   # we also need to disable any set A or P on $x (_find_round_parameters took
  783.   # them already into account), since these would interfere, too
  784.   delete $x->{_a}; delete $x->{_p};
  785.   # need to disable $upgrade in BigInt, to avoid deep recursion
  786.   local $Math::BigInt::upgrade = undef;
  787.   local $Math::BigFloat::downgrade = undef;
  788.  
  789.   # upgrade $x if $x is not a BigFloat (handle BigInt input)
  790.   if (!$x->isa('Math::BigFloat'))
  791.     {
  792.     $x = Math::BigFloat->new($x);
  793.     $self = ref($x);
  794.     }
  795.   
  796.   my $done = 0;
  797.  
  798.   # If the base is defined and an integer, try to calculate integer result
  799.   # first. This is very fast, and in case the real result was found, we can
  800.   # stop right here.
  801.   if (defined $base && $base->is_int() && $x->is_int())
  802.     {
  803.     my $i = $MBI->_copy( $x->{_m} );
  804.     $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
  805.     my $int = Math::BigInt->bzero();
  806.     $int->{value} = $i;
  807.     $int->blog($base->as_number());
  808.     # if ($exact)
  809.     if ($base->as_number()->bpow($int) == $x)
  810.       {
  811.       # found result, return it
  812.       $x->{_m} = $int->{value};
  813.       $x->{_e} = $MBI->_zero();
  814.       $x->{_es} = '+';
  815.       $x->bnorm();
  816.       $done = 1;
  817.       }
  818.     }
  819.  
  820.   if ($done == 0)
  821.     {
  822.     # first calculate the log to base e (using reduction by 10 (and probably 2))
  823.     $self->_log_10($x,$scale);
  824.  
  825.     # and if a different base was requested, convert it
  826.     if (defined $base)
  827.       {
  828.       $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
  829.       # not ln, but some other base (don't modify $base)
  830.       $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
  831.       }
  832.     }
  833.  
  834.   # shortcut to not run through _find_round_parameters again
  835.   if (defined $params[0])
  836.     {
  837.     $x->bround($params[0],$params[2]);        # then round accordingly
  838.     }
  839.   else
  840.     {
  841.     $x->bfround($params[1],$params[2]);        # then round accordingly
  842.     }
  843.   if ($fallback)
  844.     {
  845.     # clear a/p after round, since user did not request it
  846.     delete $x->{_a}; delete $x->{_p};
  847.     }
  848.   # restore globals
  849.   $$abr = $ab; $$pbr = $pb;
  850.  
  851.   $x;
  852.   }
  853.  
  854. sub _log
  855.   {
  856.   # internal log function to calculate ln() based on Taylor series.
  857.   # Modifies $x in place.
  858.   my ($self,$x,$scale) = @_;
  859.  
  860.   # in case of $x == 1, result is 0
  861.   return $x->bzero() if $x->is_one();
  862.  
  863.   # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
  864.  
  865.   # u = x-1, v = x+1
  866.   #              _                               _
  867.   # Taylor:     |    u    1   u^3   1   u^5       |
  868.   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 0
  869.   #             |_   v    3   v^3   5   v^5      _|
  870.  
  871.   # This takes much more steps to calculate the result and is thus not used
  872.   # u = x-1
  873.   #              _                               _
  874.   # Taylor:     |    u    1   u^2   1   u^3       |
  875.   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 1/2
  876.   #             |_   x    2   x^2   3   x^3      _|
  877.  
  878.   my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
  879.  
  880.   $v = $x->copy(); $v->binc();        # v = x+1
  881.   $x->bdec(); $u = $x->copy();        # u = x-1; x = x-1
  882.   $x->bdiv($v,$scale);            # first term: u/v
  883.   $below = $v->copy();
  884.   $over = $u->copy();
  885.   $u *= $u; $v *= $v;                # u^2, v^2
  886.   $below->bmul($v);                # u^3, v^3
  887.   $over->bmul($u);
  888.   $factor = $self->new(3); $f = $self->new(2);
  889.  
  890.   my $steps = 0 if DEBUG;  
  891.   $limit = $self->new("1E-". ($scale-1));
  892.   while (3 < 5)
  893.     {
  894.     # we calculate the next term, and add it to the last
  895.     # when the next term is below our limit, it won't affect the outcome
  896.     # anymore, so we stop
  897.  
  898.     # calculating the next term simple from over/below will result in quite
  899.     # a time hog if the input has many digits, since over and below will
  900.     # accumulate more and more digits, and the result will also have many
  901.     # digits, but in the end it is rounded to $scale digits anyway. So if we
  902.     # round $over and $below first, we save a lot of time for the division
  903.     # (not with log(1.2345), but try log (123**123) to see what I mean. This
  904.     # can introduce a rounding error if the division result would be f.i.
  905.     # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
  906.     # if we truncated $over and $below we might get 0.12345. Does this matter
  907.     # for the end result? So we give $over and $below 4 more digits to be
  908.     # on the safe side (unscientific error handling as usual... :+D
  909.     
  910.     $next = $over->copy->bround($scale+4)->bdiv(
  911.       $below->copy->bmul($factor)->bround($scale+4), 
  912.       $scale);
  913.  
  914. ## old version:    
  915. ##    $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
  916.  
  917.     last if $next->bacmp($limit) <= 0;
  918.  
  919.     delete $next->{_a}; delete $next->{_p};
  920.     $x->badd($next);
  921.     # calculate things for the next term
  922.     $over *= $u; $below *= $v; $factor->badd($f);
  923.     if (DEBUG)
  924.       {
  925.       $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
  926.       }
  927.     }
  928.   $x->bmul($f);                    # $x *= 2
  929.   print "took $steps steps\n" if DEBUG;
  930.   }
  931.  
  932. sub _log_10
  933.   {
  934.   # Internal log function based on reducing input to the range of 0.1 .. 9.99
  935.   # and then "correcting" the result to the proper one. Modifies $x in place.
  936.   my ($self,$x,$scale) = @_;
  937.  
  938.   # taking blog() from numbers greater than 10 takes a *very long* time, so we
  939.   # break the computation down into parts based on the observation that:
  940.   #  blog(x*y) = blog(x) + blog(y)
  941.   # We set $y here to multiples of 10 so that $x is below 1 (the smaller $x is
  942.   # the faster it get's, especially because 2*$x takes about 10 times as long,
  943.   # so by dividing $x by 10 we make it at least factor 100 faster...)
  944.  
  945.   # The same observation is valid for numbers smaller than 0.1 (e.g. computing
  946.   # log(1) is fastest, and the farther away we get from 1, the longer it takes)
  947.   # so we also 'break' this down by multiplying $x with 10 and subtract the
  948.   # log(10) afterwards to get the correct result.
  949.  
  950.   # calculate nr of digits before dot
  951.   my $dbd = $MBI->_num($x->{_e});
  952.   $dbd = -$dbd if $x->{_es} eq '-';
  953.   $dbd += $MBI->_len($x->{_m});
  954.  
  955.   # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
  956.   # infinite recursion
  957.  
  958.   my $calc = 1;                    # do some calculation?
  959.  
  960.   # disable the shortcut for 10, since we need log(10) and this would recurse
  961.   # infinitely deep
  962.   if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
  963.     {
  964.     $dbd = 0;                    # disable shortcut
  965.     # we can use the cached value in these cases
  966.     if ($scale <= $LOG_10_A)
  967.       {
  968.       $x->bzero(); $x->badd($LOG_10);
  969.       $calc = 0;                 # no need to calc, but round
  970.       }
  971.     }
  972.   else
  973.     {
  974.     # disable the shortcut for 2, since we maybe have it cached
  975.     if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
  976.       {
  977.       $dbd = 0;                    # disable shortcut
  978.       # we can use the cached value in these cases
  979.       if ($scale <= $LOG_2_A)
  980.         {
  981.         $x->bzero(); $x->badd($LOG_2);
  982.         $calc = 0;                 # no need to calc, but round
  983.         }
  984.       }
  985.     }
  986.  
  987.   # if $x = 0.1, we know the result must be 0-log(10)
  988.   if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
  989.       $MBI->_is_one($x->{_m}))
  990.     {
  991.     $dbd = 0;                    # disable shortcut
  992.     # we can use the cached value in these cases
  993.     if ($scale <= $LOG_10_A)
  994.       {
  995.       $x->bzero(); $x->bsub($LOG_10);
  996.       $calc = 0;                 # no need to calc, but round
  997.       }
  998.     }
  999.  
  1000.   return if $calc == 0;                # already have the result
  1001.  
  1002.   # default: these correction factors are undef and thus not used
  1003.   my $l_10;                # value of ln(10) to A of $scale
  1004.   my $l_2;                # value of ln(2) to A of $scale
  1005.  
  1006.   # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
  1007.   # so don't do this shortcut for 1 or 0
  1008.   if (($dbd > 1) || ($dbd < 0))
  1009.     {
  1010.     # convert our cached value to an object if not already (avoid doing this
  1011.     # at import() time, since not everybody needs this)
  1012.     $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
  1013.  
  1014.     #print "x = $x, dbd = $dbd, calc = $calc\n";
  1015.     # got more than one digit before the dot, or more than one zero after the
  1016.     # dot, so do:
  1017.     #  log(123)    == log(1.23) + log(10) * 2
  1018.     #  log(0.0123) == log(1.23) - log(10) * 2
  1019.   
  1020.     if ($scale <= $LOG_10_A)
  1021.       {
  1022.       # use cached value
  1023.       $l_10 = $LOG_10->copy();        # copy for mul
  1024.       }
  1025.     else
  1026.       {
  1027.       # else: slower, compute it (but don't cache it, because it could be big)
  1028.       # also disable downgrade for this code path
  1029.       local $Math::BigFloat::downgrade = undef;
  1030.       $l_10 = $self->new(10)->blog(undef,$scale);    # scale+4, actually
  1031.       }
  1032.     $dbd-- if ($dbd > 1);         # 20 => dbd=2, so make it dbd=1    
  1033.     $l_10->bmul( $self->new($dbd));    # log(10) * (digits_before_dot-1)
  1034.     my $dbd_sign = '+';
  1035.     if ($dbd < 0)
  1036.       {
  1037.       $dbd = -$dbd;
  1038.       $dbd_sign = '-';
  1039.       }
  1040.     ($x->{_e}, $x->{_es}) = 
  1041.     _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
  1042.  
  1043.     }
  1044.  
  1045.   # Now: 0.1 <= $x < 10 (and possible correction in l_10)
  1046.  
  1047.   ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
  1048.   ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
  1049.  
  1050.   $HALF = $self->new($HALF) unless ref($HALF);
  1051.  
  1052.   my $twos = 0;                # default: none (0 times)    
  1053.   my $two = $self->new(2);
  1054.   while ($x->bacmp($HALF) <= 0)
  1055.     {
  1056.     $twos--; $x->bmul($two);
  1057.     }
  1058.   while ($x->bacmp($two) >= 0)
  1059.     {
  1060.     $twos++; $x->bdiv($two,$scale+4);        # keep all digits
  1061.     }
  1062.   # $twos > 0 => did mul 2, < 0 => did div 2 (never both)
  1063.   # calculate correction factor based on ln(2)
  1064.   if ($twos != 0)
  1065.     {
  1066.     $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
  1067.     if ($scale <= $LOG_2_A)
  1068.       {
  1069.       # use cached value
  1070.       $l_2 = $LOG_2->copy();            # copy for mul
  1071.       }
  1072.     else
  1073.       {
  1074.       # else: slower, compute it (but don't cache it, because it could be big)
  1075.       # also disable downgrade for this code path
  1076.       local $Math::BigFloat::downgrade = undef;
  1077.       $l_2 = $two->blog(undef,$scale);    # scale+4, actually
  1078.       }
  1079.     $l_2->bmul($twos);        # * -2 => subtract, * 2 => add
  1080.     }
  1081.   
  1082.   $self->_log($x,$scale);            # need to do the "normal" way
  1083.   $x->badd($l_10) if defined $l_10;         # correct it by ln(10)
  1084.   $x->badd($l_2) if defined $l_2;        # and maybe by ln(2)
  1085.   # all done, $x contains now the result
  1086.   }
  1087.  
  1088. sub blcm 
  1089.   { 
  1090.   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
  1091.   # does not modify arguments, but returns new object
  1092.   # Lowest Common Multiplicator
  1093.  
  1094.   my ($self,@arg) = objectify(0,@_);
  1095.   my $x = $self->new(shift @arg);
  1096.   while (@arg) { $x = _lcm($x,shift @arg); } 
  1097.   $x;
  1098.   }
  1099.  
  1100. sub bgcd 
  1101.   { 
  1102.   # (BFLOAT or num_str, BFLOAT or num_str) return BINT
  1103.   # does not modify arguments, but returns new object
  1104.   # GCD -- Euclids algorithm Knuth Vol 2 pg 296
  1105.    
  1106.   my ($self,@arg) = objectify(0,@_);
  1107.   my $x = $self->new(shift @arg);
  1108.   while (@arg) { $x = _gcd($x,shift @arg); } 
  1109.   $x;
  1110.   }
  1111.  
  1112. ##############################################################################
  1113.  
  1114. sub _e_add
  1115.   {
  1116.   # Internal helper sub to take two positive integers and their signs and
  1117.   # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')), 
  1118.   # output ($CALC,('+'|'-'))
  1119.   my ($x,$y,$xs,$ys) = @_;
  1120.  
  1121.   # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
  1122.   if ($xs eq $ys)
  1123.     {
  1124.     $x = $MBI->_add ($x, $y );        # a+b
  1125.     # the sign follows $xs
  1126.     return ($x, $xs);
  1127.     }
  1128.  
  1129.   my $a = $MBI->_acmp($x,$y);
  1130.   if ($a > 0)
  1131.     {
  1132.     $x = $MBI->_sub ($x , $y);                # abs sub
  1133.     }
  1134.   elsif ($a == 0)
  1135.     {
  1136.     $x = $MBI->_zero();                    # result is 0
  1137.     $xs = '+';
  1138.     }
  1139.   else # a < 0
  1140.     {
  1141.     $x = $MBI->_sub ( $y, $x, 1 );            # abs sub
  1142.     $xs = $ys;
  1143.     }
  1144.   ($x,$xs);
  1145.   }
  1146.  
  1147. sub _e_sub
  1148.   {
  1149.   # Internal helper sub to take two positive integers and their signs and
  1150.   # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')), 
  1151.   # output ($CALC,('+'|'-'))
  1152.   my ($x,$y,$xs,$ys) = @_;
  1153.  
  1154.   # flip sign
  1155.   $ys =~ tr/+-/-+/;
  1156.   _e_add($x,$y,$xs,$ys);        # call add (does subtract now)
  1157.   }
  1158.  
  1159. ###############################################################################
  1160. # is_foo methods (is_negative, is_positive are inherited from BigInt)
  1161.  
  1162. sub is_int
  1163.   {
  1164.   # return true if arg (BFLOAT or num_str) is an integer
  1165.   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
  1166.  
  1167.   return 1 if ($x->{sign} =~ /^[+-]$/) &&    # NaN and +-inf aren't
  1168.     $x->{_es} eq '+';                # 1e-1 => no integer
  1169.   0;
  1170.   }
  1171.  
  1172. sub is_zero
  1173.   {
  1174.   # return true if arg (BFLOAT or num_str) is zero
  1175.   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
  1176.  
  1177.   return 1 if $x->{sign} eq '+' && $MBI->_is_zero($x->{_m});
  1178.   0;
  1179.   }
  1180.  
  1181. sub is_one
  1182.   {
  1183.   # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
  1184.   my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
  1185.  
  1186.   $sign = '+' if !defined $sign || $sign ne '-';
  1187.   return 1
  1188.    if ($x->{sign} eq $sign && 
  1189.     $MBI->_is_zero($x->{_e}) && $MBI->_is_one($x->{_m})); 
  1190.   0;
  1191.   }
  1192.  
  1193. sub is_odd
  1194.   {
  1195.   # return true if arg (BFLOAT or num_str) is odd or false if even
  1196.   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
  1197.   
  1198.   return 1 if ($x->{sign} =~ /^[+-]$/) &&        # NaN & +-inf aren't
  1199.     ($MBI->_is_zero($x->{_e}) && $MBI->_is_odd($x->{_m})); 
  1200.   0;
  1201.   }
  1202.  
  1203. sub is_even
  1204.   {
  1205.   # return true if arg (BINT or num_str) is even or false if odd
  1206.   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
  1207.  
  1208.   return 0 if $x->{sign} !~ /^[+-]$/;            # NaN & +-inf aren't
  1209.   return 1 if ($x->{_es} eq '+'                 # 123.45 is never
  1210.      && $MBI->_is_even($x->{_m}));            # but 1200 is
  1211.   0;
  1212.   }
  1213.  
  1214. sub bmul 
  1215.   { 
  1216.   # multiply two numbers -- stolen from Knuth Vol 2 pg 233
  1217.   # (BINT or num_str, BINT or num_str) return BINT
  1218.   
  1219.   # set up parameters
  1220.   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
  1221.   # objectify is costly, so avoid it
  1222.   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
  1223.     {
  1224.     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
  1225.     }
  1226.  
  1227.   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
  1228.  
  1229.   # inf handling
  1230.   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
  1231.     {
  1232.     return $x->bnan() if $x->is_zero() || $y->is_zero(); 
  1233.     # result will always be +-inf:
  1234.     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
  1235.     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
  1236.     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
  1237.     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
  1238.     return $x->binf('-');
  1239.     }
  1240.   # handle result = 0
  1241.   return $x->bzero() if $x->is_zero() || $y->is_zero();
  1242.   
  1243.   return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
  1244.    ((!$x->isa($self)) || (!$y->isa($self)));
  1245.  
  1246.   # aEb * cEd = (a*c)E(b+d)
  1247.   $MBI->_mul($x->{_m},$y->{_m});
  1248.   ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
  1249.  
  1250.   # adjust sign:
  1251.   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
  1252.   return $x->bnorm()->round($a,$p,$r,$y);
  1253.   }
  1254.  
  1255. sub bdiv 
  1256.   {
  1257.   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return 
  1258.   # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
  1259.  
  1260.   # set up parameters
  1261.   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
  1262.   # objectify is costly, so avoid it
  1263.   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
  1264.     {
  1265.     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
  1266.     }
  1267.  
  1268.   return $self->_div_inf($x,$y)
  1269.    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
  1270.  
  1271.   # x== 0 # also: or y == 1 or y == -1
  1272.   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
  1273.  
  1274.   # upgrade ?
  1275.   return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
  1276.  
  1277.   # we need to limit the accuracy to protect against overflow
  1278.   my $fallback = 0;
  1279.   my (@params,$scale);
  1280.   ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
  1281.  
  1282.   return $x if $x->is_nan();        # error in _find_round_parameters?
  1283.  
  1284.   # no rounding at all, so must use fallback
  1285.   if (scalar @params == 0)
  1286.     {
  1287.     # simulate old behaviour
  1288.     $params[0] = $self->div_scale();    # and round to it as accuracy
  1289.     $scale = $params[0]+4;         # at least four more for proper round
  1290.     $params[2] = $r;            # round mode by caller or undef
  1291.     $fallback = 1;            # to clear a/p afterwards
  1292.     }
  1293.   else
  1294.     {
  1295.     # the 4 below is empirical, and there might be cases where it is not
  1296.     # enough...
  1297.     $scale = abs($params[0] || $params[1]) + 4;    # take whatever is defined
  1298.     }
  1299.  
  1300.   my $rem; $rem = $self->bzero() if wantarray;
  1301.  
  1302.   $y = $self->new($y) unless $y->isa('Math::BigFloat');
  1303.  
  1304.   my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
  1305.   $scale = $lx if $lx > $scale;
  1306.   $scale = $ly if $ly > $scale;
  1307.   my $diff = $ly - $lx;
  1308.   $scale += $diff if $diff > 0;        # if lx << ly, but not if ly << lx!
  1309.   
  1310.   # cases like $x /= $x (but not $x /= $y!) were wrong due to modifying $x
  1311.   # twice below)
  1312.   require Scalar::Util;
  1313.   if (Scalar::Util::refaddr($x) == Scalar::Util::refaddr($y)) 
  1314.     {
  1315.     $x->bone();                # x/x => 1, rem 0
  1316.     }
  1317.   else
  1318.     {
  1319.  
  1320.     # make copy of $x in case of list context for later reminder calculation
  1321.     if (wantarray && !$y->is_one())
  1322.       {
  1323.       $rem = $x->copy();
  1324.       }
  1325.  
  1326.     $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+'; 
  1327.  
  1328.     # check for / +-1 ( +/- 1E0)
  1329.     if (!$y->is_one())
  1330.       {
  1331.       # promote BigInts and it's subclasses (except when already a BigFloat)
  1332.       $y = $self->new($y) unless $y->isa('Math::BigFloat'); 
  1333.  
  1334.       # calculate the result to $scale digits and then round it
  1335.       # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
  1336.       $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
  1337.       $MBI->_div ($x->{_m},$y->{_m});    # a/c
  1338.  
  1339.       # correct exponent of $x
  1340.       ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
  1341.       # correct for 10**scale
  1342.       ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
  1343.       $x->bnorm();        # remove trailing 0's
  1344.       }
  1345.     } # ende else $x != $y
  1346.  
  1347.   # shortcut to not run through _find_round_parameters again
  1348.   if (defined $params[0])
  1349.     {
  1350.     delete $x->{_a};                 # clear before round
  1351.     $x->bround($params[0],$params[2]);        # then round accordingly
  1352.     }
  1353.   else
  1354.     {
  1355.     delete $x->{_p};                 # clear before round
  1356.     $x->bfround($params[1],$params[2]);        # then round accordingly
  1357.     }
  1358.   if ($fallback)
  1359.     {
  1360.     # clear a/p after round, since user did not request it
  1361.     delete $x->{_a}; delete $x->{_p};
  1362.     }
  1363.  
  1364.   if (wantarray)
  1365.     {
  1366.     if (!$y->is_one())
  1367.       {
  1368.       $rem->bmod($y,@params);            # copy already done
  1369.       }
  1370.     if ($fallback)
  1371.       {
  1372.       # clear a/p after round, since user did not request it
  1373.       delete $rem->{_a}; delete $rem->{_p};
  1374.       }
  1375.     return ($x,$rem);
  1376.     }
  1377.   $x;
  1378.   }
  1379.  
  1380. sub bmod 
  1381.   {
  1382.   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder 
  1383.  
  1384.   # set up parameters
  1385.   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
  1386.   # objectify is costly, so avoid it
  1387.   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
  1388.     {
  1389.     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
  1390.     }
  1391.  
  1392.   # handle NaN, inf, -inf
  1393.   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
  1394.     {
  1395.     my ($d,$re) = $self->SUPER::_div_inf($x,$y);
  1396.     $x->{sign} = $re->{sign};
  1397.     $x->{_e} = $re->{_e};
  1398.     $x->{_m} = $re->{_m};
  1399.     return $x->round($a,$p,$r,$y);
  1400.     } 
  1401.   if ($y->is_zero())
  1402.     {
  1403.     return $x->bnan() if $x->is_zero();
  1404.     return $x;
  1405.     }
  1406.   return $x->bzero() if $y->is_one() || $x->is_zero();
  1407.  
  1408.   my $cmp = $x->bacmp($y);            # equal or $x < $y?
  1409.   return $x->bzero($a,$p) if $cmp == 0;        # $x == $y => result 0
  1410.  
  1411.   # only $y of the operands negative? 
  1412.   my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
  1413.  
  1414.   $x->{sign} = $y->{sign};                # calc sign first
  1415.   return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0;    # $x < $y => result $x
  1416.   
  1417.   my $ym = $MBI->_copy($y->{_m});
  1418.   
  1419.   # 2e1 => 20
  1420.   $MBI->_lsft( $ym, $y->{_e}, 10) 
  1421.    if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
  1422.  
  1423.   # if $y has digits after dot
  1424.   my $shifty = 0;            # correct _e of $x by this
  1425.   if ($y->{_es} eq '-')            # has digits after dot
  1426.     {
  1427.     # 123 % 2.5 => 1230 % 25 => 5 => 0.5
  1428.     $shifty = $MBI->_num($y->{_e});     # no more digits after dot
  1429.     $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
  1430.     }
  1431.   # $ym is now mantissa of $y based on exponent 0
  1432.  
  1433.   my $shiftx = 0;            # correct _e of $x by this
  1434.   if ($x->{_es} eq '-')            # has digits after dot
  1435.     {
  1436.     # 123.4 % 20 => 1234 % 200
  1437.     $shiftx = $MBI->_num($x->{_e});    # no more digits after dot
  1438.     $MBI->_lsft($ym, $x->{_e}, 10);    # 123 => 1230
  1439.     }
  1440.   # 123e1 % 20 => 1230 % 20
  1441.   if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
  1442.     {
  1443.     $MBI->_lsft( $x->{_m}, $x->{_e},10);    # es => '+' here
  1444.     }
  1445.  
  1446.   $x->{_e} = $MBI->_new($shiftx);
  1447.   $x->{_es} = '+'; 
  1448.   $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
  1449.   $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
  1450.   
  1451.   # now mantissas are equalized, exponent of $x is adjusted, so calc result
  1452.  
  1453.   $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
  1454.  
  1455.   $x->{sign} = '+' if $MBI->_is_zero($x->{_m});        # fix sign for -0
  1456.   $x->bnorm();
  1457.  
  1458.   if ($neg != 0)    # one of them negative => correct in place
  1459.     {
  1460.     my $r = $y - $x;
  1461.     $x->{_m} = $r->{_m};
  1462.     $x->{_e} = $r->{_e};
  1463.     $x->{_es} = $r->{_es};
  1464.     $x->{sign} = '+' if $MBI->_is_zero($x->{_m});    # fix sign for -0
  1465.     $x->bnorm();
  1466.     }
  1467.  
  1468.   $x->round($a,$p,$r,$y);    # round and return
  1469.   }
  1470.  
  1471. sub broot
  1472.   {
  1473.   # calculate $y'th root of $x
  1474.   
  1475.   # set up parameters
  1476.   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
  1477.   # objectify is costly, so avoid it
  1478.   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
  1479.     {
  1480.     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
  1481.     }
  1482.  
  1483.   # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
  1484.   return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
  1485.          $y->{sign} !~ /^\+$/;
  1486.  
  1487.   return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
  1488.   
  1489.   # we need to limit the accuracy to protect against overflow
  1490.   my $fallback = 0;
  1491.   my (@params,$scale);
  1492.   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
  1493.  
  1494.   return $x if $x->is_nan();        # error in _find_round_parameters?
  1495.  
  1496.   # no rounding at all, so must use fallback
  1497.   if (scalar @params == 0) 
  1498.     {
  1499.     # simulate old behaviour
  1500.     $params[0] = $self->div_scale();    # and round to it as accuracy
  1501.     $scale = $params[0]+4;         # at least four more for proper round
  1502.     $params[2] = $r;            # iound mode by caller or undef
  1503.     $fallback = 1;            # to clear a/p afterwards
  1504.     }
  1505.   else
  1506.     {
  1507.     # the 4 below is empirical, and there might be cases where it is not
  1508.     # enough...
  1509.     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
  1510.     }
  1511.  
  1512.   # when user set globals, they would interfere with our calculation, so
  1513.   # disable them and later re-enable them
  1514.   no strict 'refs';
  1515.   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
  1516.   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
  1517.   # we also need to disable any set A or P on $x (_find_round_parameters took
  1518.   # them already into account), since these would interfere, too
  1519.   delete $x->{_a}; delete $x->{_p};
  1520.   # need to disable $upgrade in BigInt, to avoid deep recursion
  1521.   local $Math::BigInt::upgrade = undef;    # should be really parent class vs MBI
  1522.  
  1523.   # remember sign and make $x positive, since -4 ** (1/2) => -2
  1524.   my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
  1525.  
  1526.   my $is_two = 0;
  1527.   if ($y->isa('Math::BigFloat'))
  1528.     {
  1529.     $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
  1530.     }
  1531.   else
  1532.     {
  1533.     $is_two = ($y == 2);
  1534.     }
  1535.  
  1536.   # normal square root if $y == 2:
  1537.   if ($is_two)
  1538.     {
  1539.     $x->bsqrt($scale+4);
  1540.     }
  1541.   elsif ($y->is_one('-'))
  1542.     {
  1543.     # $x ** -1 => 1/$x
  1544.     my $u = $self->bone()->bdiv($x,$scale);
  1545.     # copy private parts over
  1546.     $x->{_m} = $u->{_m};
  1547.     $x->{_e} = $u->{_e};
  1548.     $x->{_es} = $u->{_es};
  1549.     }
  1550.   else
  1551.     {
  1552.     # calculate the broot() as integer result first, and if it fits, return
  1553.     # it rightaway (but only if $x and $y are integer):
  1554.  
  1555.     my $done = 0;                # not yet
  1556.     if ($y->is_int() && $x->is_int())
  1557.       {
  1558.       my $i = $MBI->_copy( $x->{_m} );
  1559.       $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
  1560.       my $int = Math::BigInt->bzero();
  1561.       $int->{value} = $i;
  1562.       $int->broot($y->as_number());
  1563.       # if ($exact)
  1564.       if ($int->copy()->bpow($y) == $x)
  1565.         {
  1566.         # found result, return it
  1567.         $x->{_m} = $int->{value};
  1568.         $x->{_e} = $MBI->_zero();
  1569.         $x->{_es} = '+';
  1570.         $x->bnorm();
  1571.         $done = 1;
  1572.         }
  1573.       }
  1574.     if ($done == 0)
  1575.       {
  1576.       my $u = $self->bone()->bdiv($y,$scale+4);
  1577.       delete $u->{_a}; delete $u->{_p};         # otherwise it conflicts
  1578.       $x->bpow($u,$scale+4);                    # el cheapo
  1579.       }
  1580.     }
  1581.   $x->bneg() if $sign == 1;
  1582.   
  1583.   # shortcut to not run through _find_round_parameters again
  1584.   if (defined $params[0])
  1585.     {
  1586.     $x->bround($params[0],$params[2]);        # then round accordingly
  1587.     }
  1588.   else
  1589.     {
  1590.     $x->bfround($params[1],$params[2]);        # then round accordingly
  1591.     }
  1592.   if ($fallback)
  1593.     {
  1594.     # clear a/p after round, since user did not request it
  1595.     delete $x->{_a}; delete $x->{_p};
  1596.     }
  1597.   # restore globals
  1598.   $$abr = $ab; $$pbr = $pb;
  1599.   $x;
  1600.   }
  1601.  
  1602. sub bsqrt
  1603.   { 
  1604.   # calculate square root
  1605.   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
  1606.  
  1607.   return $x->bnan() if $x->{sign} !~ /^[+]/;    # NaN, -inf or < 0
  1608.   return $x if $x->{sign} eq '+inf';        # sqrt(inf) == inf
  1609.   return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
  1610.  
  1611.   # we need to limit the accuracy to protect against overflow
  1612.   my $fallback = 0;
  1613.   my (@params,$scale);
  1614.   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
  1615.  
  1616.   return $x if $x->is_nan();        # error in _find_round_parameters?
  1617.  
  1618.   # no rounding at all, so must use fallback
  1619.   if (scalar @params == 0) 
  1620.     {
  1621.     # simulate old behaviour
  1622.     $params[0] = $self->div_scale();    # and round to it as accuracy
  1623.     $scale = $params[0]+4;         # at least four more for proper round
  1624.     $params[2] = $r;            # round mode by caller or undef
  1625.     $fallback = 1;            # to clear a/p afterwards
  1626.     }
  1627.   else
  1628.     {
  1629.     # the 4 below is empirical, and there might be cases where it is not
  1630.     # enough...
  1631.     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
  1632.     }
  1633.  
  1634.   # when user set globals, they would interfere with our calculation, so
  1635.   # disable them and later re-enable them
  1636.   no strict 'refs';
  1637.   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
  1638.   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
  1639.   # we also need to disable any set A or P on $x (_find_round_parameters took
  1640.   # them already into account), since these would interfere, too
  1641.   delete $x->{_a}; delete $x->{_p};
  1642.   # need to disable $upgrade in BigInt, to avoid deep recursion
  1643.   local $Math::BigInt::upgrade = undef;    # should be really parent class vs MBI
  1644.  
  1645.   my $i = $MBI->_copy( $x->{_m} );
  1646.   $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
  1647.   my $xas = Math::BigInt->bzero();
  1648.   $xas->{value} = $i;
  1649.  
  1650.   my $gs = $xas->copy()->bsqrt();    # some guess
  1651.  
  1652.   if (($x->{_es} ne '-')        # guess can't be accurate if there are
  1653.                     # digits after the dot
  1654.    && ($xas->bacmp($gs * $gs) == 0))    # guess hit the nail on the head?
  1655.     {
  1656.     # exact result, copy result over to keep $x
  1657.     $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
  1658.     $x->bnorm();
  1659.     # shortcut to not run through _find_round_parameters again
  1660.     if (defined $params[0])
  1661.       {
  1662.       $x->bround($params[0],$params[2]);    # then round accordingly
  1663.       }
  1664.     else
  1665.       {
  1666.       $x->bfround($params[1],$params[2]);    # then round accordingly
  1667.       }
  1668.     if ($fallback)
  1669.       {
  1670.       # clear a/p after round, since user did not request it
  1671.       delete $x->{_a}; delete $x->{_p};
  1672.       }
  1673.     # re-enable A and P, upgrade is taken care of by "local"
  1674.     ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
  1675.     return $x;
  1676.     }
  1677.  
  1678.   # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
  1679.   # of the result by multipyling the input by 100 and then divide the integer
  1680.   # result of sqrt(input) by 10. Rounding afterwards returns the real result.
  1681.  
  1682.   # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
  1683.   my $y1 = $MBI->_copy($x->{_m});
  1684.  
  1685.   my $length = $MBI->_len($y1);
  1686.   
  1687.   # Now calculate how many digits the result of sqrt(y1) would have
  1688.   my $digits = int($length / 2);
  1689.  
  1690.   # But we need at least $scale digits, so calculate how many are missing
  1691.   my $shift = $scale - $digits;
  1692.  
  1693.   # That should never happen (we take care of integer guesses above)
  1694.   # $shift = 0 if $shift < 0; 
  1695.  
  1696.   # Multiply in steps of 100, by shifting left two times the "missing" digits
  1697.   my $s2 = $shift * 2;
  1698.  
  1699.   # We now make sure that $y1 has the same odd or even number of digits than
  1700.   # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
  1701.   # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
  1702.   # steps of 10. The length of $x does not count, since an even or odd number
  1703.   # of digits before the dot is not changed by adding an even number of digits
  1704.   # after the dot (the result is still odd or even digits long).
  1705.   $s2++ if $MBI->_is_odd($x->{_e});
  1706.  
  1707.   $MBI->_lsft( $y1, $MBI->_new($s2), 10);
  1708.  
  1709.   # now take the square root and truncate to integer
  1710.   $y1 = $MBI->_sqrt($y1);
  1711.  
  1712.   # By "shifting" $y1 right (by creating a negative _e) we calculate the final
  1713.   # result, which is than later rounded to the desired scale.
  1714.  
  1715.   # calculate how many zeros $x had after the '.' (or before it, depending
  1716.   # on sign of $dat, the result should have half as many:
  1717.   my $dat = $MBI->_num($x->{_e});
  1718.   $dat = -$dat if $x->{_es} eq '-';
  1719.   $dat += $length;
  1720.  
  1721.   if ($dat > 0)
  1722.     {
  1723.     # no zeros after the dot (e.g. 1.23, 0.49 etc)
  1724.     # preserve half as many digits before the dot than the input had 
  1725.     # (but round this "up")
  1726.     $dat = int(($dat+1)/2);
  1727.     }
  1728.   else
  1729.     {
  1730.     $dat = int(($dat)/2);
  1731.     }
  1732.   $dat -= $MBI->_len($y1);
  1733.   if ($dat < 0)
  1734.     {
  1735.     $dat = abs($dat);
  1736.     $x->{_e} = $MBI->_new( $dat );
  1737.     $x->{_es} = '-';
  1738.     }
  1739.   else
  1740.     {    
  1741.     $x->{_e} = $MBI->_new( $dat );
  1742.     $x->{_es} = '+';
  1743.     }
  1744.   $x->{_m} = $y1;
  1745.   $x->bnorm();
  1746.  
  1747.   # shortcut to not run through _find_round_parameters again
  1748.   if (defined $params[0])
  1749.     {
  1750.     $x->bround($params[0],$params[2]);        # then round accordingly
  1751.     }
  1752.   else
  1753.     {
  1754.     $x->bfround($params[1],$params[2]);        # then round accordingly
  1755.     }
  1756.   if ($fallback)
  1757.     {
  1758.     # clear a/p after round, since user did not request it
  1759.     delete $x->{_a}; delete $x->{_p};
  1760.     }
  1761.   # restore globals
  1762.   $$abr = $ab; $$pbr = $pb;
  1763.   $x;
  1764.   }
  1765.  
  1766. sub bfac
  1767.   {
  1768.   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
  1769.   # compute factorial number, modifies first argument
  1770.  
  1771.   # set up parameters
  1772.   my ($self,$x,@r) = (ref($_[0]),@_);
  1773.   # objectify is costly, so avoid it
  1774.   ($self,$x,@r) = objectify(1,@_) if !ref($x);
  1775.  
  1776.  return $x if $x->{sign} eq '+inf';    # inf => inf
  1777.   return $x->bnan() 
  1778.     if (($x->{sign} ne '+') ||        # inf, NaN, <0 etc => NaN
  1779.      ($x->{_es} ne '+'));        # digits after dot?
  1780.  
  1781.   # use BigInt's bfac() for faster calc
  1782.   if (! $MBI->_is_zero($x->{_e}))
  1783.     {
  1784.     $MBI->_lsft($x->{_m}, $x->{_e},10);    # change 12e1 to 120e0
  1785.     $x->{_e} = $MBI->_zero();        # normalize
  1786.     $x->{_es} = '+';
  1787.     }
  1788.   $MBI->_fac($x->{_m});            # calculate factorial
  1789.   $x->bnorm()->round(@r);         # norm again and round result
  1790.   }
  1791.  
  1792. sub _pow
  1793.   {
  1794.   # Calculate a power where $y is a non-integer, like 2 ** 0.5
  1795.   my ($x,$y,$a,$p,$r) = @_;
  1796.   my $self = ref($x);
  1797.  
  1798.   # if $y == 0.5, it is sqrt($x)
  1799.   $HALF = $self->new($HALF) unless ref($HALF);
  1800.   return $x->bsqrt($a,$p,$r,$y) if $y->bcmp($HALF) == 0;
  1801.  
  1802.   # Using:
  1803.   # a ** x == e ** (x * ln a)
  1804.  
  1805.   # u = y * ln x
  1806.   #                _                         _
  1807.   # Taylor:       |   u    u^2    u^3         |
  1808.   # x ** y  = 1 + |  --- + --- + ----- + ...  |
  1809.   #               |_  1    1*2   1*2*3       _|
  1810.  
  1811.   # we need to limit the accuracy to protect against overflow
  1812.   my $fallback = 0;
  1813.   my ($scale,@params);
  1814.   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
  1815.     
  1816.   return $x if $x->is_nan();        # error in _find_round_parameters?
  1817.  
  1818.   # no rounding at all, so must use fallback
  1819.   if (scalar @params == 0)
  1820.     {
  1821.     # simulate old behaviour
  1822.     $params[0] = $self->div_scale();    # and round to it as accuracy
  1823.     $params[1] = undef;            # disable P
  1824.     $scale = $params[0]+4;         # at least four more for proper round
  1825.     $params[2] = $r;            # round mode by caller or undef
  1826.     $fallback = 1;            # to clear a/p afterwards
  1827.     }
  1828.   else
  1829.     {
  1830.     # the 4 below is empirical, and there might be cases where it is not
  1831.     # enough...
  1832.     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
  1833.     }
  1834.  
  1835.   # when user set globals, they would interfere with our calculation, so
  1836.   # disable them and later re-enable them
  1837.   no strict 'refs';
  1838.   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
  1839.   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
  1840.   # we also need to disable any set A or P on $x (_find_round_parameters took
  1841.   # them already into account), since these would interfere, too
  1842.   delete $x->{_a}; delete $x->{_p};
  1843.   # need to disable $upgrade in BigInt, to avoid deep recursion
  1844.   local $Math::BigInt::upgrade = undef;
  1845.  
  1846.   my ($limit,$v,$u,$below,$factor,$next,$over);
  1847.  
  1848.   $u = $x->copy()->blog(undef,$scale)->bmul($y);
  1849.   $v = $self->bone();                # 1
  1850.   $factor = $self->new(2);            # 2
  1851.   $x->bone();                    # first term: 1
  1852.  
  1853.   $below = $v->copy();
  1854.   $over = $u->copy();
  1855.  
  1856.   $limit = $self->new("1E-". ($scale-1));
  1857.   #my $steps = 0;
  1858.   while (3 < 5)
  1859.     {
  1860.     # we calculate the next term, and add it to the last
  1861.     # when the next term is below our limit, it won't affect the outcome
  1862.     # anymore, so we stop
  1863.     $next = $over->copy()->bdiv($below,$scale);
  1864.     last if $next->bacmp($limit) <= 0;
  1865.     $x->badd($next);
  1866.     # calculate things for the next term
  1867.     $over *= $u; $below *= $factor; $factor->binc();
  1868.  
  1869.     last if $x->{sign} !~ /^[-+]$/;
  1870.  
  1871.     #$steps++;
  1872.     }
  1873.   
  1874.   # shortcut to not run through _find_round_parameters again
  1875.   if (defined $params[0])
  1876.     {
  1877.     $x->bround($params[0],$params[2]);        # then round accordingly
  1878.     }
  1879.   else
  1880.     {
  1881.     $x->bfround($params[1],$params[2]);        # then round accordingly
  1882.     }
  1883.   if ($fallback)
  1884.     {
  1885.     # clear a/p after round, since user did not request it
  1886.     delete $x->{_a}; delete $x->{_p};
  1887.     }
  1888.   # restore globals
  1889.   $$abr = $ab; $$pbr = $pb;
  1890.   $x;
  1891.   }
  1892.  
  1893. sub bpow 
  1894.   {
  1895.   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
  1896.   # compute power of two numbers, second arg is used as integer
  1897.   # modifies first argument
  1898.  
  1899.   # set up parameters
  1900.   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
  1901.   # objectify is costly, so avoid it
  1902.   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
  1903.     {
  1904.     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
  1905.     }
  1906.  
  1907.   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
  1908.   return $x if $x->{sign} =~ /^[+-]inf$/;
  1909.   
  1910.   # -2 ** -2 => NaN
  1911.   return $x->bnan() if $x->{sign} eq '-' && $y->{sign} eq '-';
  1912.  
  1913.   # cache the result of is_zero
  1914.   my $y_is_zero = $y->is_zero();
  1915.   return $x->bone() if $y_is_zero;
  1916.   return $x         if $x->is_one() || $y->is_one();
  1917.  
  1918.   my $x_is_zero = $x->is_zero();
  1919.   return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int();        # non-integer power
  1920.  
  1921.   my $y1 = $y->as_number()->{value};            # make MBI part
  1922.  
  1923.   # if ($x == -1)
  1924.   if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
  1925.     {
  1926.     # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
  1927.     return $MBI->_is_odd($y1) ? $x : $x->babs(1);
  1928.     }
  1929.   if ($x_is_zero)
  1930.     {
  1931.     return $x->bone() if $y_is_zero;
  1932.     return $x if $y->{sign} eq '+';     # 0**y => 0 (if not y <= 0)
  1933.     # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
  1934.     return $x->binf();
  1935.     }
  1936.  
  1937.   my $new_sign = '+';
  1938.   $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
  1939.  
  1940.   # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
  1941.   $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
  1942.   $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
  1943.  
  1944.   $x->{sign} = $new_sign;
  1945.   $x->bnorm();
  1946.   if ($y->{sign} eq '-')
  1947.     {
  1948.     # modify $x in place!
  1949.     my $z = $x->copy(); $x->bone();
  1950.     return $x->bdiv($z,$a,$p,$r);    # round in one go (might ignore y's A!)
  1951.     }
  1952.   $x->round($a,$p,$r,$y);
  1953.   }
  1954.  
  1955. ###############################################################################
  1956. # rounding functions
  1957.  
  1958. sub bfround
  1959.   {
  1960.   # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
  1961.   # $n == 0 means round to integer
  1962.   # expects and returns normalized numbers!
  1963.   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
  1964.  
  1965.   return $x if $x->modify('bfround');
  1966.  
  1967.   my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
  1968.   return $x if !defined $scale;            # no-op
  1969.  
  1970.   # never round a 0, +-inf, NaN
  1971.   if ($x->is_zero())
  1972.     {
  1973.     $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
  1974.     return $x; 
  1975.     }
  1976.   return $x if $x->{sign} !~ /^[+-]$/;
  1977.  
  1978.   # don't round if x already has lower precision
  1979.   return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
  1980.  
  1981.   $x->{_p} = $scale;            # remember round in any case
  1982.   delete $x->{_a};            # and clear A
  1983.   if ($scale < 0)
  1984.     {
  1985.     # round right from the '.'
  1986.  
  1987.     return $x if $x->{_es} eq '+';        # e >= 0 => nothing to round
  1988.  
  1989.     $scale = -$scale;                # positive for simplicity
  1990.     my $len = $MBI->_len($x->{_m});        # length of mantissa
  1991.  
  1992.     # the following poses a restriction on _e, but if _e is bigger than a
  1993.     # scalar, you got other problems (memory etc) anyway
  1994.     my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e})));    # digits after dot
  1995.     my $zad = 0;                # zeros after dot
  1996.     $zad = $dad - $len if (-$dad < -$len);    # for 0.00..00xxx style
  1997.    
  1998.     # p rint "scale $scale dad $dad zad $zad len $len\n";
  1999.     # number  bsstr   len zad dad    
  2000.     # 0.123   123e-3    3   0 3
  2001.     # 0.0123  123e-4    3   1 4
  2002.     # 0.001   1e-3      1   2 3
  2003.     # 1.23    123e-2    3   0 2
  2004.     # 1.2345  12345e-4    5   0 4
  2005.  
  2006.     # do not round after/right of the $dad
  2007.     return $x if $scale > $dad;            # 0.123, scale >= 3 => exit
  2008.  
  2009.     # round to zero if rounding inside the $zad, but not for last zero like:
  2010.     # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
  2011.     return $x->bzero() if $scale < $zad;
  2012.     if ($scale == $zad)            # for 0.006, scale -3 and trunc
  2013.       {
  2014.       $scale = -$len;
  2015.       }
  2016.     else
  2017.       {
  2018.       # adjust round-point to be inside mantissa
  2019.       if ($zad != 0)
  2020.         {
  2021.     $scale = $scale-$zad;
  2022.         }
  2023.       else
  2024.         {
  2025.         my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;    # digits before dot
  2026.     $scale = $dbd+$scale;
  2027.         }
  2028.       }
  2029.     }
  2030.   else
  2031.     {
  2032.     # round left from the '.'
  2033.  
  2034.     # 123 => 100 means length(123) = 3 - $scale (2) => 1
  2035.  
  2036.     my $dbt = $MBI->_len($x->{_m}); 
  2037.     # digits before dot 
  2038.     my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
  2039.     # should be the same, so treat it as this 
  2040.     $scale = 1 if $scale == 0; 
  2041.     # shortcut if already integer 
  2042.     return $x if $scale == 1 && $dbt <= $dbd; 
  2043.     # maximum digits before dot 
  2044.     ++$dbd;
  2045.  
  2046.     if ($scale > $dbd) 
  2047.        { 
  2048.        # not enough digits before dot, so round to zero 
  2049.        return $x->bzero; 
  2050.        }
  2051.     elsif ( $scale == $dbd )
  2052.        { 
  2053.        # maximum 
  2054.        $scale = -$dbt; 
  2055.        } 
  2056.     else
  2057.        { 
  2058.        $scale = $dbd - $scale; 
  2059.        }
  2060.     }
  2061.   # pass sign to bround for rounding modes '+inf' and '-inf'
  2062.   my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
  2063.   $m->bround($scale,$mode);
  2064.   $x->{_m} = $m->{value};            # get our mantissa back
  2065.   $x->bnorm();
  2066.   }
  2067.  
  2068. sub bround
  2069.   {
  2070.   # accuracy: preserve $N digits, and overwrite the rest with 0's
  2071.   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
  2072.  
  2073.   if (($_[0] || 0) < 0)
  2074.     {
  2075.     require Carp; Carp::croak ('bround() needs positive accuracy');
  2076.     }
  2077.  
  2078.   my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
  2079.   return $x if !defined $scale;                # no-op
  2080.  
  2081.   return $x if $x->modify('bround');
  2082.  
  2083.   # scale is now either $x->{_a}, $accuracy, or the user parameter
  2084.   # test whether $x already has lower accuracy, do nothing in this case 
  2085.   # but do round if the accuracy is the same, since a math operation might
  2086.   # want to round a number with A=5 to 5 digits afterwards again
  2087.   return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
  2088.  
  2089.   # scale < 0 makes no sense
  2090.   # never round a +-inf, NaN
  2091.   return $x if ($scale < 0) ||    $x->{sign} !~ /^[+-]$/;
  2092.  
  2093.   # 1: $scale == 0 => keep all digits
  2094.   # 2: never round a 0
  2095.   # 3: if we should keep more digits than the mantissa has, do nothing
  2096.   if ($scale == 0 || $x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
  2097.     {
  2098.     $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
  2099.     return $x; 
  2100.     }
  2101.  
  2102.   # pass sign to bround for '+inf' and '-inf' rounding modes
  2103.   my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
  2104.  
  2105.   $m->bround($scale,$mode);        # round mantissa
  2106.   $x->{_m} = $m->{value};        # get our mantissa back
  2107.   $x->{_a} = $scale;            # remember rounding
  2108.   delete $x->{_p};            # and clear P
  2109.   $x->bnorm();                # del trailing zeros gen. by bround()
  2110.   }
  2111.  
  2112. sub bfloor
  2113.   {
  2114.   # return integer less or equal then $x
  2115.   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
  2116.  
  2117.   return $x if $x->modify('bfloor');
  2118.    
  2119.   return $x if $x->{sign} !~ /^[+-]$/;    # nan, +inf, -inf
  2120.  
  2121.   # if $x has digits after dot
  2122.   if ($x->{_es} eq '-')
  2123.     {
  2124.     $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
  2125.     $x->{_e} = $MBI->_zero();            # trunc/norm    
  2126.     $x->{_es} = '+';                # abs e
  2127.     $MBI->_inc($x->{_m}) if $x->{sign} eq '-';    # increment if negative
  2128.     }
  2129.   $x->round($a,$p,$r);
  2130.   }
  2131.  
  2132. sub bceil
  2133.   {
  2134.   # return integer greater or equal then $x
  2135.   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
  2136.  
  2137.   return $x if $x->modify('bceil');
  2138.   return $x if $x->{sign} !~ /^[+-]$/;    # nan, +inf, -inf
  2139.  
  2140.   # if $x has digits after dot
  2141.   if ($x->{_es} eq '-')
  2142.     {
  2143.     $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
  2144.     $x->{_e} = $MBI->_zero();            # trunc/norm    
  2145.     $x->{_es} = '+';                # abs e
  2146.     $MBI->_inc($x->{_m}) if $x->{sign} eq '+';    # increment if positive
  2147.     }
  2148.   $x->round($a,$p,$r);
  2149.   }
  2150.  
  2151. sub brsft
  2152.   {
  2153.   # shift right by $y (divide by power of $n)
  2154.   
  2155.   # set up parameters
  2156.   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
  2157.   # objectify is costly, so avoid it
  2158.   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
  2159.     {
  2160.     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
  2161.     }
  2162.  
  2163.   return $x if $x->modify('brsft');
  2164.   return $x if $x->{sign} !~ /^[+-]$/;    # nan, +inf, -inf
  2165.  
  2166.   $n = 2 if !defined $n; $n = $self->new($n);
  2167.   $x->bdiv($n->bpow($y),$a,$p,$r,$y);
  2168.   }
  2169.  
  2170. sub blsft
  2171.   {
  2172.   # shift left by $y (multiply by power of $n)
  2173.   
  2174.   # set up parameters
  2175.   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
  2176.   # objectify is costly, so avoid it
  2177.   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
  2178.     {
  2179.     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
  2180.     }
  2181.  
  2182.   return $x if $x->modify('blsft');
  2183.   return $x if $x->{sign} !~ /^[+-]$/;    # nan, +inf, -inf
  2184.  
  2185.   $n = 2 if !defined $n; $n = $self->new($n);
  2186.   $x->bmul($n->bpow($y),$a,$p,$r,$y);
  2187.   }
  2188.  
  2189. ###############################################################################
  2190.  
  2191. sub DESTROY
  2192.   {
  2193.   # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
  2194.   }
  2195.  
  2196. sub AUTOLOAD
  2197.   {
  2198.   # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
  2199.   # or falling back to MBI::bxxx()
  2200.   my $name = $AUTOLOAD;
  2201.  
  2202.   $name =~ s/(.*):://;    # split package
  2203.   my $c = $1 || $class;
  2204.   no strict 'refs';
  2205.   $c->import() if $IMPORT == 0;
  2206.   if (!method_alias($name))
  2207.     {
  2208.     if (!defined $name)
  2209.       {
  2210.       # delayed load of Carp and avoid recursion    
  2211.       require Carp;
  2212.       Carp::croak ("$c: Can't call a method without name");
  2213.       }
  2214.     if (!method_hand_up($name))
  2215.       {
  2216.       # delayed load of Carp and avoid recursion    
  2217.       require Carp;
  2218.       Carp::croak ("Can't call $c\-\>$name, not a valid method");
  2219.       }
  2220.     # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
  2221.     $name =~ s/^f/b/;
  2222.     return &{"Math::BigInt"."::$name"}(@_);
  2223.     }
  2224.   my $bname = $name; $bname =~ s/^f/b/;
  2225.   $c .= "::$name";
  2226.   *{$c} = \&{$bname};
  2227.   &{$c};    # uses @_
  2228.   }
  2229.  
  2230. sub exponent
  2231.   {
  2232.   # return a copy of the exponent
  2233.   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
  2234.  
  2235.   if ($x->{sign} !~ /^[+-]$/)
  2236.     {
  2237.     my $s = $x->{sign}; $s =~ s/^[+-]//;
  2238.     return Math::BigInt->new($s);         # -inf, +inf => +inf
  2239.     }
  2240.   Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
  2241.   }
  2242.  
  2243. sub mantissa
  2244.   {
  2245.   # return a copy of the mantissa
  2246.   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
  2247.  
  2248.   if ($x->{sign} !~ /^[+-]$/)
  2249.     {
  2250.     my $s = $x->{sign}; $s =~ s/^[+]//;
  2251.     return Math::BigInt->new($s);        # -inf, +inf => +inf
  2252.     }
  2253.   my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
  2254.   $m->bneg() if $x->{sign} eq '-';
  2255.  
  2256.   $m;
  2257.   }
  2258.  
  2259. sub parts
  2260.   {
  2261.   # return a copy of both the exponent and the mantissa
  2262.   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
  2263.  
  2264.   if ($x->{sign} !~ /^[+-]$/)
  2265.     {
  2266.     my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
  2267.     return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
  2268.     }
  2269.   my $m = Math::BigInt->bzero();
  2270.   $m->{value} = $MBI->_copy($x->{_m});
  2271.   $m->bneg() if $x->{sign} eq '-';
  2272.   ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
  2273.   }
  2274.  
  2275. ##############################################################################
  2276. # private stuff (internal use only)
  2277.  
  2278. sub import
  2279.   {
  2280.   my $self = shift;
  2281.   my $l = scalar @_;
  2282.   my $lib = ''; my @a;
  2283.   $IMPORT=1;
  2284.   for ( my $i = 0; $i < $l ; $i++)
  2285.     {
  2286.     if ( $_[$i] eq ':constant' )
  2287.       {
  2288.       # This causes overlord er load to step in. 'binary' and 'integer'
  2289.       # are handled by BigInt.
  2290.       overload::constant float => sub { $self->new(shift); }; 
  2291.       }
  2292.     elsif ($_[$i] eq 'upgrade')
  2293.       {
  2294.       # this causes upgrading
  2295.       $upgrade = $_[$i+1];        # or undef to disable
  2296.       $i++;
  2297.       }
  2298.     elsif ($_[$i] eq 'downgrade')
  2299.       {
  2300.       # this causes downgrading
  2301.       $downgrade = $_[$i+1];        # or undef to disable
  2302.       $i++;
  2303.       }
  2304.     elsif ($_[$i] eq 'lib')
  2305.       {
  2306.       # alternative library
  2307.       $lib = $_[$i+1] || '';        # default Calc
  2308.       $i++;
  2309.       }
  2310.     elsif ($_[$i] eq 'with')
  2311.       {
  2312.       # alternative class for our private parts()
  2313.       # XXX: no longer supported
  2314.       # $MBI = $_[$i+1] || 'Math::BigInt';
  2315.       $i++;
  2316.       }
  2317.     else
  2318.       {
  2319.       push @a, $_[$i];
  2320.       }
  2321.     }
  2322.  
  2323.   # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
  2324.   my $mbilib = eval { Math::BigInt->config()->{lib} };
  2325.   if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
  2326.     {
  2327.     # MBI already loaded
  2328.     Math::BigInt->import('lib',"$lib,$mbilib", 'objectify');
  2329.     }
  2330.   else
  2331.     {
  2332.     # MBI not loaded, or with ne "Math::BigInt::Calc"
  2333.     $lib .= ",$mbilib" if defined $mbilib;
  2334.     $lib =~ s/^,//;                # don't leave empty 
  2335.     
  2336.     # replacement library can handle lib statement, but also could ignore it
  2337.     
  2338.     # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
  2339.     # used in the same script, or eval inside import(). So we require MBI:
  2340.     require Math::BigInt;
  2341.     Math::BigInt->import( lib => $lib, 'objectify' );
  2342.     }
  2343.   if ($@)
  2344.     {
  2345.     require Carp; Carp::croak ("Couldn't load $lib: $! $@");
  2346.     }
  2347.   $MBI = Math::BigInt->config()->{lib};
  2348.  
  2349.   # any non :constant stuff is handled by our parent, Exporter
  2350.   # even if @_ is empty, to give it a chance
  2351.   $self->SUPER::import(@a);          # for subclasses
  2352.   $self->export_to_level(1,$self,@a);    # need this, too
  2353.   }
  2354.  
  2355. sub bnorm
  2356.   {
  2357.   # adjust m and e so that m is smallest possible
  2358.   # round number according to accuracy and precision settings
  2359.   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
  2360.  
  2361.   return $x if $x->{sign} !~ /^[+-]$/;        # inf, nan etc
  2362.  
  2363.   my $zeros = $MBI->_zeros($x->{_m});        # correct for trailing zeros
  2364.   if ($zeros != 0)
  2365.     {
  2366.     my $z = $MBI->_new($zeros);
  2367.     $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
  2368.     if ($x->{_es} eq '-')
  2369.       {
  2370.       if ($MBI->_acmp($x->{_e},$z) >= 0)
  2371.         {
  2372.         $x->{_e} = $MBI->_sub  ($x->{_e}, $z);
  2373.         $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
  2374.         }
  2375.       else
  2376.         {
  2377.         $x->{_e} = $MBI->_sub  ( $MBI->_copy($z), $x->{_e});
  2378.         $x->{_es} = '+';
  2379.         }
  2380.       }
  2381.     else
  2382.       {
  2383.       $x->{_e} = $MBI->_add  ($x->{_e}, $z);
  2384.       }
  2385.     }
  2386.   else
  2387.     {
  2388.     # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
  2389.     # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
  2390.     $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
  2391.      if $MBI->_is_zero($x->{_m});
  2392.     }
  2393.  
  2394.   $x;                    # MBI bnorm is no-op, so dont call it
  2395.   } 
  2396.  
  2397. ##############################################################################
  2398.  
  2399. sub as_hex
  2400.   {
  2401.   # return number as hexadecimal string (only for integers defined)
  2402.   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
  2403.  
  2404.   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
  2405.   return '0x0' if $x->is_zero();
  2406.  
  2407.   return $nan if $x->{_es} ne '+';        # how to do 1e-1 in hex!?
  2408.  
  2409.   my $z = $MBI->_copy($x->{_m});
  2410.   if (! $MBI->_is_zero($x->{_e}))        # > 0 
  2411.     {
  2412.     $MBI->_lsft( $z, $x->{_e},10);
  2413.     }
  2414.   $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
  2415.   $z->as_hex();
  2416.   }
  2417.  
  2418. sub as_bin
  2419.   {
  2420.   # return number as binary digit string (only for integers defined)
  2421.   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
  2422.  
  2423.   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
  2424.   return '0b0' if $x->is_zero();
  2425.  
  2426.   return $nan if $x->{_es} ne '+';        # how to do 1e-1 in hex!?
  2427.  
  2428.   my $z = $MBI->_copy($x->{_m});
  2429.   if (! $MBI->_is_zero($x->{_e}))        # > 0 
  2430.     {
  2431.     $MBI->_lsft( $z, $x->{_e},10);
  2432.     }
  2433.   $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
  2434.   $z->as_bin();
  2435.   }
  2436.  
  2437. sub as_number
  2438.   {
  2439.   # return copy as a bigint representation of this BigFloat number
  2440.   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
  2441.  
  2442.   my $z = $MBI->_copy($x->{_m});
  2443.   if ($x->{_es} eq '-')            # < 0
  2444.     {
  2445.     $MBI->_rsft( $z, $x->{_e},10);
  2446.     } 
  2447.   elsif (! $MBI->_is_zero($x->{_e}))    # > 0 
  2448.     {
  2449.     $MBI->_lsft( $z, $x->{_e},10);
  2450.     }
  2451.   $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
  2452.   $z;
  2453.   }
  2454.  
  2455. sub length
  2456.   {
  2457.   my $x = shift;
  2458.   my $class = ref($x) || $x;
  2459.   $x = $class->new(shift) unless ref($x);
  2460.  
  2461.   return 1 if $MBI->_is_zero($x->{_m});
  2462.  
  2463.   my $len = $MBI->_len($x->{_m});
  2464.   $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
  2465.   if (wantarray())
  2466.     {
  2467.     my $t = 0;
  2468.     $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
  2469.     return ($len, $t);
  2470.     }
  2471.   $len;
  2472.   }
  2473.  
  2474. 1;
  2475. __END__
  2476.  
  2477. =head1 NAME
  2478.  
  2479. Math::BigFloat - Arbitrary size floating point math package
  2480.  
  2481. =head1 SYNOPSIS
  2482.  
  2483.   use Math::BigFloat;
  2484.  
  2485.   # Number creation
  2486.   $x = Math::BigFloat->new($str);    # defaults to 0
  2487.   $nan  = Math::BigFloat->bnan();    # create a NotANumber
  2488.   $zero = Math::BigFloat->bzero();    # create a +0
  2489.   $inf = Math::BigFloat->binf();    # create a +inf
  2490.   $inf = Math::BigFloat->binf('-');    # create a -inf
  2491.   $one = Math::BigFloat->bone();    # create a +1
  2492.   $one = Math::BigFloat->bone('-');    # create a -1
  2493.  
  2494.   # Testing
  2495.   $x->is_zero();        # true if arg is +0
  2496.   $x->is_nan();            # true if arg is NaN
  2497.   $x->is_one();            # true if arg is +1
  2498.   $x->is_one('-');        # true if arg is -1
  2499.   $x->is_odd();            # true if odd, false for even
  2500.   $x->is_even();        # true if even, false for odd
  2501.   $x->is_pos();            # true if >= 0
  2502.   $x->is_neg();            # true if <  0
  2503.   $x->is_inf(sign);        # true if +inf, or -inf (default is '+')
  2504.  
  2505.   $x->bcmp($y);            # compare numbers (undef,<0,=0,>0)
  2506.   $x->bacmp($y);        # compare absolutely (undef,<0,=0,>0)
  2507.   $x->sign();            # return the sign, either +,- or NaN
  2508.   $x->digit($n);        # return the nth digit, counting from right
  2509.   $x->digit(-$n);        # return the nth digit, counting from left 
  2510.  
  2511.   # The following all modify their first argument. If you want to preserve
  2512.   # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
  2513.   # neccessary when mixing $a = $b assigments with non-overloaded math.
  2514.  
  2515.   # set 
  2516.   $x->bzero();            # set $i to 0
  2517.   $x->bnan();            # set $i to NaN
  2518.   $x->bone();                   # set $x to +1
  2519.   $x->bone('-');                # set $x to -1
  2520.   $x->binf();                   # set $x to inf
  2521.   $x->binf('-');                # set $x to -inf
  2522.  
  2523.   $x->bneg();            # negation
  2524.   $x->babs();            # absolute value
  2525.   $x->bnorm();            # normalize (no-op)
  2526.   $x->bnot();            # two's complement (bit wise not)
  2527.   $x->binc();            # increment x by 1
  2528.   $x->bdec();            # decrement x by 1
  2529.   
  2530.   $x->badd($y);            # addition (add $y to $x)
  2531.   $x->bsub($y);            # subtraction (subtract $y from $x)
  2532.   $x->bmul($y);            # multiplication (multiply $x by $y)
  2533.   $x->bdiv($y);            # divide, set $x to quotient
  2534.                 # return (quo,rem) or quo if scalar
  2535.  
  2536.   $x->bmod($y);            # modulus ($x % $y)
  2537.   $x->bpow($y);            # power of arguments ($x ** $y)
  2538.   $x->blsft($y);        # left shift
  2539.   $x->brsft($y);        # right shift 
  2540.                 # return (quo,rem) or quo if scalar
  2541.   
  2542.   $x->blog();            # logarithm of $x to base e (Euler's number)
  2543.   $x->blog($base);        # logarithm of $x to base $base (f.i. 2)
  2544.   
  2545.   $x->band($y);            # bit-wise and
  2546.   $x->bior($y);            # bit-wise inclusive or
  2547.   $x->bxor($y);            # bit-wise exclusive or
  2548.   $x->bnot();            # bit-wise not (two's complement)
  2549.  
  2550.   $x->bsqrt();            # calculate square-root
  2551.   $x->broot($y);        # $y'th root of $x (e.g. $y == 3 => cubic root)
  2552.   $x->bfac();            # factorial of $x (1*2*3*4*..$x)
  2553.  
  2554.   $x->bround($N);         # accuracy: preserve $N digits
  2555.   $x->bfround($N);        # precision: round to the $Nth digit
  2556.  
  2557.   $x->bfloor();            # return integer less or equal than $x
  2558.   $x->bceil();            # return integer greater or equal than $x
  2559.  
  2560.   # The following do not modify their arguments:
  2561.  
  2562.   bgcd(@values);        # greatest common divisor
  2563.   blcm(@values);        # lowest common multiplicator
  2564.   
  2565.   $x->bstr();            # return string
  2566.   $x->bsstr();            # return string in scientific notation
  2567.  
  2568.   $x->as_int();            # return $x as BigInt 
  2569.   $x->exponent();        # return exponent as BigInt
  2570.   $x->mantissa();        # return mantissa as BigInt
  2571.   $x->parts();            # return (mantissa,exponent) as BigInt
  2572.  
  2573.   $x->length();            # number of digits (w/o sign and '.')
  2574.   ($l,$f) = $x->length();    # number of digits, and length of fraction    
  2575.  
  2576.   $x->precision();        # return P of $x (or global, if P of $x undef)
  2577.   $x->precision($n);        # set P of $x to $n
  2578.   $x->accuracy();        # return A of $x (or global, if A of $x undef)
  2579.   $x->accuracy($n);        # set A $x to $n
  2580.  
  2581.   # these get/set the appropriate global value for all BigFloat objects
  2582.   Math::BigFloat->precision();    # Precision
  2583.   Math::BigFloat->accuracy();    # Accuracy
  2584.   Math::BigFloat->round_mode();    # rounding mode
  2585.  
  2586. =head1 DESCRIPTION
  2587.  
  2588. All operators (inlcuding basic math operations) are overloaded if you
  2589. declare your big floating point numbers as
  2590.  
  2591.   $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
  2592.  
  2593. Operations with overloaded operators preserve the arguments, which is
  2594. exactly what you expect.
  2595.  
  2596. =head2 Canonical notation
  2597.  
  2598. Input to these routines are either BigFloat objects, or strings of the
  2599. following four forms:
  2600.  
  2601. =over 2
  2602.  
  2603. =item *
  2604.  
  2605. C</^[+-]\d+$/>
  2606.  
  2607. =item *
  2608.  
  2609. C</^[+-]\d+\.\d*$/>
  2610.  
  2611. =item *
  2612.  
  2613. C</^[+-]\d+E[+-]?\d+$/>
  2614.  
  2615. =item *
  2616.  
  2617. C</^[+-]\d*\.\d+E[+-]?\d+$/>
  2618.  
  2619. =back
  2620.  
  2621. all with optional leading and trailing zeros and/or spaces. Additonally,
  2622. numbers are allowed to have an underscore between any two digits.
  2623.  
  2624. Empty strings as well as other illegal numbers results in 'NaN'.
  2625.  
  2626. bnorm() on a BigFloat object is now effectively a no-op, since the numbers 
  2627. are always stored in normalized form. On a string, it creates a BigFloat 
  2628. object.
  2629.  
  2630. =head2 Output
  2631.  
  2632. Output values are BigFloat objects (normalized), except for bstr() and bsstr().
  2633.  
  2634. The string output will always have leading and trailing zeros stripped and drop
  2635. a plus sign. C<bstr()> will give you always the form with a decimal point,
  2636. while C<bsstr()> (s for scientific) gives you the scientific notation.
  2637.  
  2638.     Input            bstr()        bsstr()
  2639.     '-0'            '0'        '0E1'
  2640.        '  -123 123 123'    '-123123123'    '-123123123E0'
  2641.     '00.0123'        '0.0123'    '123E-4'
  2642.     '123.45E-2'        '1.2345'    '12345E-4'
  2643.     '10E+3'            '10000'        '1E4'
  2644.  
  2645. Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
  2646. C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
  2647. return either undef, <0, 0 or >0 and are suited for sort.
  2648.  
  2649. Actual math is done by using the class defined with C<with => Class;> (which
  2650. defaults to BigInts) to represent the mantissa and exponent.
  2651.  
  2652. The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to 
  2653. represent the result when input arguments are not numbers, as well as 
  2654. the result of dividing by zero.
  2655.  
  2656. =head2 C<mantissa()>, C<exponent()> and C<parts()>
  2657.  
  2658. C<mantissa()> and C<exponent()> return the said parts of the BigFloat 
  2659. as BigInts such that:
  2660.  
  2661.     $m = $x->mantissa();
  2662.     $e = $x->exponent();
  2663.     $y = $m * ( 10 ** $e );
  2664.     print "ok\n" if $x == $y;
  2665.  
  2666. C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
  2667.  
  2668. A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
  2669.  
  2670. Currently the mantissa is reduced as much as possible, favouring higher
  2671. exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
  2672. This might change in the future, so do not depend on it.
  2673.  
  2674. =head2 Accuracy vs. Precision
  2675.  
  2676. See also: L<Rounding|Rounding>.
  2677.  
  2678. Math::BigFloat supports both precision and accuracy. For a full documentation,
  2679. examples and tips on these topics please see the large section in
  2680. L<Math::BigInt>.
  2681.  
  2682. Since things like sqrt(2) or 1/3 must presented with a limited precision lest
  2683. a operation consumes all resources, each operation produces no more than
  2684. the requested number of digits.
  2685.  
  2686. Please refer to BigInt's documentation for the precedence rules of which
  2687. accuracy/precision setting will be used.
  2688.  
  2689. If there is no gloabl precision set, B<and> the operation inquestion was not
  2690. called with a requested precision or accuracy, B<and> the input $x has no
  2691. accuracy or precision set, then a fallback parameter will be used. For
  2692. historical reasons, it is called C<div_scale> and can be accessed via:
  2693.  
  2694.     $d = Math::BigFloat->div_scale();        # query
  2695.     Math::BigFloat->div_scale($n);            # set to $n digits
  2696.  
  2697. The default value is 40 digits.
  2698.  
  2699. In case the result of one operation has more precision than specified,
  2700. it is rounded. The rounding mode taken is either the default mode, or the one
  2701. supplied to the operation after the I<scale>:
  2702.  
  2703.     $x = Math::BigFloat->new(2);
  2704.     Math::BigFloat->precision(5);        # 5 digits max
  2705.     $y = $x->copy()->bdiv(3);        # will give 0.66666
  2706.     $y = $x->copy()->bdiv(3,6);        # will give 0.666666
  2707.     $y = $x->copy()->bdiv(3,6,'odd');    # will give 0.666667
  2708.     Math::BigFloat->round_mode('zero');
  2709.     $y = $x->copy()->bdiv(3,6);        # will give 0.666666
  2710.  
  2711. =head2 Rounding
  2712.  
  2713. =over 2
  2714.  
  2715. =item ffround ( +$scale )
  2716.  
  2717. Rounds to the $scale'th place left from the '.', counting from the dot.
  2718. The first digit is numbered 1. 
  2719.  
  2720. =item ffround ( -$scale )
  2721.  
  2722. Rounds to the $scale'th place right from the '.', counting from the dot.
  2723.  
  2724. =item ffround ( 0 )
  2725.  
  2726. Rounds to an integer.
  2727.  
  2728. =item fround  ( +$scale )
  2729.  
  2730. Preserves accuracy to $scale digits from the left (aka significant digits)
  2731. and pads the rest with zeros. If the number is between 1 and -1, the
  2732. significant digits count from the first non-zero after the '.'
  2733.  
  2734. =item fround  ( -$scale ) and fround ( 0 )
  2735.  
  2736. These are effectively no-ops.
  2737.  
  2738. =back
  2739.  
  2740. All rounding functions take as a second parameter a rounding mode from one of
  2741. the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
  2742.  
  2743. The default rounding mode is 'even'. By using
  2744. C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
  2745. mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
  2746. no longer supported.
  2747. The second parameter to the round functions then overrides the default
  2748. temporarily. 
  2749.  
  2750. The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
  2751. 'trunc' as rounding mode to make it equivalent to:
  2752.  
  2753.     $x = 2.5;
  2754.     $y = int($x) + 2;
  2755.  
  2756. You can override this by passing the desired rounding mode as parameter to
  2757. C<as_number()>:
  2758.  
  2759.     $x = Math::BigFloat->new(2.5);
  2760.     $y = $x->as_number('odd');    # $y = 3
  2761.  
  2762. =head1 EXAMPLES
  2763.  
  2764.   # not ready yet
  2765.  
  2766. =head1 Autocreating constants
  2767.  
  2768. After C<use Math::BigFloat ':constant'> all the floating point constants
  2769. in the given scope are converted to C<Math::BigFloat>. This conversion
  2770. happens at compile time.
  2771.  
  2772. In particular
  2773.  
  2774.   perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
  2775.  
  2776. prints the value of C<2E-100>. Note that without conversion of 
  2777. constants the expression 2E-100 will be calculated as normal floating point 
  2778. number.
  2779.  
  2780. Please note that ':constant' does not affect integer constants, nor binary 
  2781. nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
  2782. work.
  2783.  
  2784. =head2 Math library
  2785.  
  2786. Math with the numbers is done (by default) by a module called
  2787. Math::BigInt::Calc. This is equivalent to saying:
  2788.  
  2789.     use Math::BigFloat lib => 'Calc';
  2790.  
  2791. You can change this by using:
  2792.  
  2793.     use Math::BigFloat lib => 'BitVect';
  2794.  
  2795. The following would first try to find Math::BigInt::Foo, then
  2796. Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
  2797.  
  2798.     use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
  2799.  
  2800. Calc.pm uses as internal format an array of elements of some decimal base
  2801. (usually 1e7, but this might be differen for some systems) with the least
  2802. significant digit first, while BitVect.pm uses a bit vector of base 2, most
  2803. significant bit first. Other modules might use even different means of
  2804. representing the numbers. See the respective module documentation for further
  2805. details.
  2806.  
  2807. Please note that Math::BigFloat does B<not> use the denoted library itself,
  2808. but it merely passes the lib argument to Math::BigInt. So, instead of the need
  2809. to do:
  2810.  
  2811.     use Math::BigInt lib => 'GMP';
  2812.     use Math::BigFloat;
  2813.  
  2814. you can roll it all into one line:
  2815.  
  2816.     use Math::BigFloat lib => 'GMP';
  2817.  
  2818. It is also possible to just require Math::BigFloat:
  2819.  
  2820.     require Math::BigFloat;
  2821.  
  2822. This will load the neccessary things (like BigInt) when they are needed, and
  2823. automatically.
  2824.  
  2825. Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
  2826. you ever wanted to know about loading a different library.
  2827.  
  2828. =head2 Using Math::BigInt::Lite
  2829.  
  2830. It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
  2831.  
  2832.         # 1
  2833.         use Math::BigFloat with => 'Math::BigInt::Lite';
  2834.  
  2835. There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
  2836. can combine these if you want. For instance, you may want to use
  2837. Math::BigInt objects in your main script, too.
  2838.  
  2839.         # 2
  2840.         use Math::BigInt;
  2841.         use Math::BigFloat with => 'Math::BigInt::Lite';
  2842.  
  2843. Of course, you can combine this with the C<lib> parameter.
  2844.  
  2845.         # 3
  2846.         use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
  2847.  
  2848. There is no need for a "use Math::BigInt;" statement, even if you want to
  2849. use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
  2850. always loads it. But if you add it, add it B<before>:
  2851.  
  2852.         # 4
  2853.         use Math::BigInt;
  2854.         use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
  2855.  
  2856. Notice that the module with the last C<lib> will "win" and thus
  2857. it's lib will be used if the lib is available:
  2858.  
  2859.         # 5
  2860.         use Math::BigInt lib => 'Bar,Baz';
  2861.         use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
  2862.  
  2863. That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
  2864. words, Math::BigFloat will try to retain previously loaded libs when you
  2865. don't specify it onem but if you specify one, it will try to load them.
  2866.  
  2867. Actually, the lib loading order would be "Bar,Baz,Calc", and then
  2868. "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
  2869. same as trying the latter load alone, except for the fact that one of Bar or
  2870. Baz might be loaded needlessly in an intermidiate step (and thus hang around
  2871. and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
  2872. will still be tried to be loaded, but this is not as time/memory consuming as
  2873. actually loading one of them. Still, this type of usage is not recommended due
  2874. to these issues.
  2875.  
  2876. The old way (loading the lib only in BigInt) still works though:
  2877.  
  2878.         # 6
  2879.         use Math::BigInt lib => 'Bar,Baz';
  2880.         use Math::BigFloat;
  2881.  
  2882. You can even load Math::BigInt afterwards:
  2883.  
  2884.         # 7
  2885.         use Math::BigFloat;
  2886.         use Math::BigInt lib => 'Bar,Baz';
  2887.  
  2888. But this has the same problems like #5, it will first load Calc
  2889. (Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
  2890. Baz, depending on which of them works and is usable/loadable. Since this
  2891. loads Calc unnecc., it is not recommended.
  2892.  
  2893. Since it also possible to just require Math::BigFloat, this poses the question
  2894. about what libary this will use:
  2895.  
  2896.     require Math::BigFloat;
  2897.     my $x = Math::BigFloat->new(123); $x += 123;
  2898.  
  2899. It will use Calc. Please note that the call to import() is still done, but
  2900. only when you use for the first time some Math::BigFloat math (it is triggered
  2901. via any constructor, so the first time you create a Math::BigFloat, the load
  2902. will happen in the background). This means:
  2903.  
  2904.     require Math::BigFloat;
  2905.     Math::BigFloat->import ( lib => 'Foo,Bar' );
  2906.  
  2907. would be the same as:
  2908.  
  2909.     use Math::BigFloat lib => 'Foo, Bar';
  2910.  
  2911. But don't try to be clever to insert some operations in between:
  2912.  
  2913.     require Math::BigFloat;
  2914.     my $x = Math::BigFloat->bone() + 4;        # load BigInt and Calc
  2915.     Math::BigFloat->import( lib => 'Pari' );    # load Pari, too
  2916.     $x = Math::BigFloat->bone()+4;            # now use Pari
  2917.  
  2918. While this works, it loads Calc needlessly. But maybe you just wanted that?
  2919.  
  2920. B<Examples #3 is highly recommended> for daily usage.
  2921.  
  2922. =head1 BUGS
  2923.  
  2924. Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
  2925.  
  2926. =head1 CAVEATS
  2927.  
  2928. =over 1
  2929.  
  2930. =item stringify, bstr()
  2931.  
  2932. Both stringify and bstr() now drop the leading '+'. The old code would return
  2933. '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
  2934. reasoning and details.
  2935.  
  2936. =item bdiv
  2937.  
  2938. The following will probably not do what you expect:
  2939.  
  2940.     print $c->bdiv(123.456),"\n";
  2941.  
  2942. It prints both quotient and reminder since print works in list context. Also,
  2943. bdiv() will modify $c, so be carefull. You probably want to use
  2944.     
  2945.     print $c / 123.456,"\n";
  2946.     print scalar $c->bdiv(123.456),"\n";  # or if you want to modify $c
  2947.  
  2948. instead.
  2949.  
  2950. =item Modifying and =
  2951.  
  2952. Beware of:
  2953.  
  2954.     $x = Math::BigFloat->new(5);
  2955.     $y = $x;
  2956.  
  2957. It will not do what you think, e.g. making a copy of $x. Instead it just makes
  2958. a second reference to the B<same> object and stores it in $y. Thus anything
  2959. that modifies $x will modify $y (except overloaded math operators), and vice
  2960. versa. See L<Math::BigInt> for details and how to avoid that.
  2961.  
  2962. =item bpow
  2963.  
  2964. C<bpow()> now modifies the first argument, unlike the old code which left
  2965. it alone and only returned the result. This is to be consistent with
  2966. C<badd()> etc. The first will modify $x, the second one won't:
  2967.  
  2968.     print bpow($x,$i),"\n";     # modify $x
  2969.     print $x->bpow($i),"\n";     # ditto
  2970.     print $x ** $i,"\n";        # leave $x alone 
  2971.  
  2972. =back
  2973.  
  2974. =head1 SEE ALSO
  2975.  
  2976. L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
  2977. L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and  L<Math::BigInt::GMP>.
  2978.  
  2979. The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
  2980. because they solve the autoupgrading/downgrading issue, at least partly.
  2981.  
  2982. The package at
  2983. L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
  2984. more documentation including a full version history, testcases, empty
  2985. subclass files and benchmarks.
  2986.  
  2987. =head1 LICENSE
  2988.  
  2989. This program is free software; you may redistribute it and/or modify it under
  2990. the same terms as Perl itself.
  2991.  
  2992. =head1 AUTHORS
  2993.  
  2994. Mark Biggar, overloaded interface by Ilya Zakharevich.
  2995. Completely rewritten by Tels http://bloodgate.com in 2001, 2002, and still
  2996. at it in 2003.
  2997.  
  2998. =cut
  2999.