lib/Math/BigFloat.pm

   1 package Math::BigFloat;
   2
   3 #
   4 # Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
   5 #
   6
   7 # The following hash values are internally used:
   8 #   _e: exponent (BigInt)
   9 #   _m: mantissa (absolute BigInt)
  10 # sign: +,-,"NaN" if not a number
  11 #   _a: accuracy
  12 #   _p: precision
  13 #   _f: flags, used to signal MBI not to touch our private parts
  14
  15 $VERSION = '1.38';
  16 require 5.005;
  17 use Exporter;
  18 @ISA =       qw(Exporter Math::BigInt);
  19
  20 use strict;
  21 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode/;
  22 use vars qw/$upgrade $downgrade/;
  23 my $class = "Math::BigFloat";
  24
  25 use overload
  26 '<=>'   =>      sub { $_[2] ?
  27                       ref($_[0])->bcmp($_[1],$_[0]) :
  28                       ref($_[0])->bcmp($_[0],$_[1])},
  29 'int'   =>      sub { $_[0]->as_number() },             # 'trunc' to bigint
  30 ;
  31
  32 ##############################################################################
  33 # global constants, flags and accessory
  34
  35 use constant MB_NEVER_ROUND => 0x0001;
  36
  37 # are NaNs ok?
  38 my $NaNOK=1;
  39 # constant for easier life
  40 my $nan = 'NaN';
  41
  42 # class constants, use Class->constant_name() to access
  43 $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
  44 $accuracy   = undef;
  45 $precision  = undef;
  46 $div_scale  = 40;
  47
  48 $upgrade = undef;
  49 $downgrade = undef;
  50 my $MBI = 'Math::BigInt'; # the package we are using for our private parts
  51                           # changable by use Math::BigFloat with => 'package'
  52
  53 ##############################################################################
  54 # the old code had $rnd_mode, so we need to support it, too
  55
  56 sub TIESCALAR   { my ($class) = @_; bless \$round_mode, $class; }
  57 sub FETCH       { return $round_mode; }
  58 sub STORE       { $rnd_mode = $_[0]->round_mode($_[1]); }
  59
  60 BEGIN
  61   {
  62   $rnd_mode   = 'even';
  63   tie $rnd_mode, 'Math::BigFloat';
  64   }
  65
  66 ##############################################################################
  67
  68 # in case we call SUPER::->foo() and this wants to call modify()
  69 # sub modify () { 0; }
  70
  71 {
  72   # valid method aliases for AUTOLOAD
  73   my %methods = map { $_ => 1 }
  74    qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
  75         fint facmp fcmp fzero fnan finf finc fdec flog ffac
  76         fceil ffloor frsft flsft fone flog
  77       /;
  78   # valid method's that can be hand-ed up (for AUTOLOAD)
  79   my %hand_ups = map { $_ => 1 }
  80    qw / is_nan is_inf is_negative is_positive
  81         accuracy precision div_scale round_mode fneg fabs babs fnot
  82         objectify upgrade downgrade
  83         bone binf bnan bzero
  84       /;
  85
  86   sub method_alias { return exists $methods{$_[0]||''}; }
  87   sub method_hand_up { return exists $hand_ups{$_[0]||''}; }
  88 }
  89
  90 ##############################################################################
  91 # constructors
  92
  93 sub new
  94   {
  95   # create a new BigFloat object from a string or another bigfloat object.
  96   # _e: exponent
  97   # _m: mantissa
  98   # sign  => sign (+/-), or "NaN"
  99
 100   my ($class,$wanted,@r) = @_;
 101
 102   # avoid numify-calls by not using || on $wanted!
 103   return $class->bzero() if !defined $wanted;   # default to 0
 104   return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
 105
 106   my $self = {}; bless $self, $class;
 107   # shortcut for bigints and its subclasses
 108   if ((ref($wanted)) && (ref($wanted) ne $class))
 109     {
 110     $self->{_m} = $wanted->as_number();         # get us a bigint copy
 111     $self->{_e} = $MBI->bzero();
 112     $self->{_m}->babs();
 113     $self->{sign} = $wanted->sign();
 114     return $self->bnorm();
 115     }
 116   # got string
 117   # handle '+inf', '-inf' first
 118   if ($wanted =~ /^[+-]?inf$/)
 119     {
 120     return $downgrade->new($wanted) if $downgrade;
 121
 122     $self->{_e} = $MBI->bzero();
 123     $self->{_m} = $MBI->bzero();
 124     $self->{sign} = $wanted;
 125     $self->{sign} = '+inf' if $self->{sign} eq 'inf';
 126     return $self->bnorm();
 127     }
 128   #print "new string '$wanted'\n";
 129   my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
 130   if (!ref $mis)
 131     {
 132     die "$wanted is not a number initialized to $class" if !$NaNOK;
 133
 134     return $downgrade->bnan() if $downgrade;
 135
 136     $self->{_e} = $MBI->bzero();
 137     $self->{_m} = $MBI->bzero();
 138     $self->{sign} = $nan;
 139     }
 140   else
 141     {
 142     # make integer from mantissa by adjusting exp, then convert to bigint
 143     # undef,undef to signal MBI that we don't need no bloody rounding
 144     $self->{_e} = $MBI->new("$$es$$ev",undef,undef);    # exponent
 145     $self->{_m} = $MBI->new("$$miv$$mfv",undef,undef);  # create mant.
 146     # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
 147     $self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0;
 148     $self->{sign} = $$mis;
 149     }
 150   # if downgrade, inf, NaN or integers go down
 151
 152   if ($downgrade && $self->{_e}->{sign} eq '+')
 153     {
 154 #   print "downgrading $$miv$$mfv"."E$$es$$ev";
 155     if ($self->{_e}->is_zero())
 156       {
 157       $self->{_m}->{sign} = $$mis;              # negative if wanted
 158       return $downgrade->new($self->{_m});
 159       }
 160     return $downgrade->new("$$mis$$miv$$mfv"."E$$es$$ev");
 161     }
 162   # print "mbf new $self->{sign} $self->{_m} e $self->{_e} ",ref($self),"\n";
 163   $self->bnorm()->round(@r);            # first normalize, then round
 164   }
 165
 166 sub _bnan
 167   {
 168   # used by parent class bone() to initialize number to 1
 169   my $self = shift;
 170   $self->{_m} = $MBI->bzero();
 171   $self->{_e} = $MBI->bzero();
 172   }
 173
 174 sub _binf
 175   {
 176   # used by parent class bone() to initialize number to 1
 177   my $self = shift;
 178   $self->{_m} = $MBI->bzero();
 179   $self->{_e} = $MBI->bzero();
 180   }
 181
 182 sub _bone
 183   {
 184   # used by parent class bone() to initialize number to 1
 185   my $self = shift;
 186   $self->{_m} = $MBI->bone();
 187   $self->{_e} = $MBI->bzero();
 188   }
 189
 190 sub _bzero
 191   {
 192   # used by parent class bone() to initialize number to 1
 193   my $self = shift;
 194   $self->{_m} = $MBI->bzero();
 195   $self->{_e} = $MBI->bone();
 196   }
 197
 198 sub isa
 199   {
 200   my ($self,$class) = @_;
 201   return if $class =~ /^Math::BigInt/;          # we aren't one of these
 202   UNIVERSAL::isa($self,$class);
 203   }
 204
 205 sub config
 206   {
 207   # return (later set?) configuration data as hash ref
 208   my $class = shift || 'Math::BigFloat';
 209
 210   my $cfg = $MBI->config();
 211
 212   no strict 'refs';
 213   $cfg->{class} = $class;
 214   $cfg->{with} = $MBI;
 215   foreach (
 216    qw/upgrade downgrade precision accuracy round_mode VERSION div_scale/)
 217     {
 218     $cfg->{lc($_)} = ${"${class}::$_"};
 219     };
 220   $cfg;
 221   }
 222
 223 ##############################################################################
 224 # string conversation
 225
 226 sub bstr
 227   {
 228   # (ref to BFLOAT or num_str ) return num_str
 229   # Convert number from internal format to (non-scientific) string format.
 230   # internal format is always normalized (no leading zeros, "-0" => "+0")
 231   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 232   #my $x = shift; my $class = ref($x) || $x;
 233   #$x = $class->new(shift) unless ref($x);
 234
 235   #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
 236   #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
 237   if ($x->{sign} !~ /^[+-]$/)
 238     {
 239     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
 240     return 'inf';                                       # +inf
 241     }
 242
 243   my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
 244
 245   my $not_zero = ! $x->is_zero();
 246   if ($not_zero)
 247     {
 248     $es = $x->{_m}->bstr();
 249     $len = CORE::length($es);
 250     if (!$x->{_e}->is_zero())
 251       {
 252       if ($x->{_e}->sign() eq '-')
 253         {
 254         $dot = '';
 255         if ($x->{_e} <= -$len)
 256           {
 257           # print "style: 0.xxxx\n";
 258           my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
 259           $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
 260           }
 261         else
 262           {
 263           # print "insert '.' at $x->{_e} in '$es'\n";
 264           substr($es,$x->{_e},0) = '.'; $cad = $x->{_e};
 265           }
 266         }
 267       else
 268         {
 269         # expand with zeros
 270         $es .= '0' x $x->{_e}; $len += $x->{_e}; $cad = 0;
 271         }
 272       }
 273     } # if not zero
 274   $es = $x->{sign}.$es if $x->{sign} eq '-';
 275   # if set accuracy or precision, pad with zeros
 276   if ((defined $x->{_a}) && ($not_zero))
 277     {
 278     # 123400 => 6, 0.1234 => 4, 0.001234 => 4
 279     my $zeros = $x->{_a} - $cad;                # cad == 0 => 12340
 280     $zeros = $x->{_a} - $len if $cad != $len;
 281     $es .= $dot.'0' x $zeros if $zeros > 0;
 282     }
 283   elsif ($x->{_p} || 0 < 0)
 284     {
 285     # 123400 => 6, 0.1234 => 4, 0.001234 => 6
 286     my $zeros = -$x->{_p} + $cad;
 287     $es .= $dot.'0' x $zeros if $zeros > 0;
 288     }
 289   $es;
 290   }
 291
 292 sub bsstr
 293   {
 294   # (ref to BFLOAT or num_str ) return num_str
 295   # Convert number from internal format to scientific string format.
 296   # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
 297   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 298   #my $x = shift; my $class = ref($x) || $x;
 299   #$x = $class->new(shift) unless ref($x);
 300
 301   #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
 302   #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
 303   if ($x->{sign} !~ /^[+-]$/)
 304     {
 305     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
 306     return 'inf';                                       # +inf
 307     }
 308   my $esign = $x->{_e}->{sign}; $esign = '' if $esign eq '-';
 309   my $sep = 'e'.$esign;
 310   my $sign = $x->{sign}; $sign = '' if $sign eq '+';
 311   $sign . $x->{_m}->bstr() . $sep . $x->{_e}->bstr();
 312   }
 313
 314 sub numify
 315   {
 316   # Make a number from a BigFloat object
 317   # simple return string and let Perl's atoi()/atof() handle the rest
 318   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 319   $x->bsstr();
 320   }
 321
 322 ##############################################################################
 323 # public stuff (usually prefixed with "b")
 324
 325 # tels 2001-08-04
 326 # todo: this must be overwritten and return NaN for non-integer values
 327 # band(), bior(), bxor(), too
 328 #sub bnot
 329 #  {
 330 #  $class->SUPER::bnot($class,@_);
 331 #  }
 332
 333 sub bcmp
 334   {
 335   # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
 336   # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
 337
 338   # set up parameters
 339   my ($self,$x,$y) = (ref($_[0]),@_);
 340   # objectify is costly, so avoid it
 341   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 342     {
 343     ($self,$x,$y) = objectify(2,@_);
 344     }
 345
 346   return $upgrade->bcmp($x,$y) if defined $upgrade &&
 347     ((!$x->isa($self)) || (!$y->isa($self)));
 348
 349   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
 350     {
 351     # handle +-inf and NaN
 352     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 353     return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
 354     return +1 if $x->{sign} eq '+inf';
 355     return -1 if $x->{sign} eq '-inf';
 356     return -1 if $y->{sign} eq '+inf';
 357     return +1;
 358     }
 359
 360   # check sign for speed first
 361   return 1 if $x->{sign} eq '+' && $y->{sign} eq '-';   # does also 0 <=> -y
 362   return -1 if $x->{sign} eq '-' && $y->{sign} eq '+';  # does also -x <=> 0
 363
 364   # shortcut
 365   my $xz = $x->is_zero();
 366   my $yz = $y->is_zero();
 367   return 0 if $xz && $yz;                               # 0 <=> 0
 368   return -1 if $xz && $y->{sign} eq '+';                # 0 <=> +y
 369   return 1 if $yz && $x->{sign} eq '+';                 # +x <=> 0
 370
 371   # adjust so that exponents are equal
 372   my $lxm = $x->{_m}->length();
 373   my $lym = $y->{_m}->length();
 374   # the numify somewhat limits our length, but makes it much faster
 375   my $lx = $lxm + $x->{_e}->numify();
 376   my $ly = $lym + $y->{_e}->numify();
 377   my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
 378   return $l <=> 0 if $l != 0;
 379
 380   # lengths (corrected by exponent) are equal
 381   # so make mantissa equal length by padding with zero (shift left)
 382   my $diff = $lxm - $lym;
 383   my $xm = $x->{_m};            # not yet copy it
 384   my $ym = $y->{_m};
 385   if ($diff > 0)
 386     {
 387     $ym = $y->{_m}->copy()->blsft($diff,10);
 388     }
 389   elsif ($diff < 0)
 390     {
 391     $xm = $x->{_m}->copy()->blsft(-$diff,10);
 392     }
 393   my $rc = $xm->bacmp($ym);
 394   $rc = -$rc if $x->{sign} eq '-';              # -124 < -123
 395   $rc <=> 0;
 396   }
 397
 398 sub bacmp
 399   {
 400   # Compares 2 values, ignoring their signs.
 401   # Returns one of undef, <0, =0, >0. (suitable for sort)
 402   # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
 403
 404   # set up parameters
 405   my ($self,$x,$y) = (ref($_[0]),@_);
 406   # objectify is costly, so avoid it
 407   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 408     {
 409     ($self,$x,$y) = objectify(2,@_);
 410     }
 411
 412   return $upgrade->bacmp($x,$y) if defined $upgrade &&
 413     ((!$x->isa($self)) || (!$y->isa($self)));
 414
 415   # handle +-inf and NaN's
 416   if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
 417     {
 418     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 419     return 0 if ($x->is_inf() && $y->is_inf());
 420     return 1 if ($x->is_inf() && !$y->is_inf());
 421     return -1;
 422     }
 423
 424   # shortcut
 425   my $xz = $x->is_zero();
 426   my $yz = $y->is_zero();
 427   return 0 if $xz && $yz;                               # 0 <=> 0
 428   return -1 if $xz && !$yz;                             # 0 <=> +y
 429   return 1 if $yz && !$xz;                              # +x <=> 0
 430
 431   # adjust so that exponents are equal
 432   my $lxm = $x->{_m}->length();
 433   my $lym = $y->{_m}->length();
 434   # the numify somewhat limits our length, but makes it much faster
 435   my $lx = $lxm + $x->{_e}->numify();
 436   my $ly = $lym + $y->{_e}->numify();
 437   my $l = $lx - $ly;
 438   return $l <=> 0 if $l != 0;
 439
 440   # lengths (corrected by exponent) are equal
 441   # so make mantissa equal-length by padding with zero (shift left)
 442   my $diff = $lxm - $lym;
 443   my $xm = $x->{_m};            # not yet copy it
 444   my $ym = $y->{_m};
 445   if ($diff > 0)
 446     {
 447     $ym = $y->{_m}->copy()->blsft($diff,10);
 448     }
 449   elsif ($diff < 0)
 450     {
 451     $xm = $x->{_m}->copy()->blsft(-$diff,10);
 452     }
 453   $xm->bacmp($ym) <=> 0;
 454   }
 455
 456 sub badd
 457   {
 458   # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
 459   # return result as BFLOAT
 460
 461   # set up parameters
 462   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
 463   # objectify is costly, so avoid it
 464   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 465     {
 466     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 467     }
 468
 469   # inf and NaN handling
 470   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
 471     {
 472     # NaN first
 473     return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 474     # inf handling
 475     if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
 476       {
 477       # +inf++inf or -inf+-inf => same, rest is NaN
 478       return $x if $x->{sign} eq $y->{sign};
 479       return $x->bnan();
 480       }
 481     # +-inf + something => +inf; something +-inf => +-inf
 482     $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
 483     return $x;
 484     }
 485
 486   return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
 487    ((!$x->isa($self)) || (!$y->isa($self)));
 488
 489   # speed: no add for 0+y or x+0
 490   return $x->bround($a,$p,$r) if $y->is_zero();         # x+0
 491   if ($x->is_zero())                                    # 0+y
 492     {
 493     # make copy, clobbering up x (modify in place!)
 494     $x->{_e} = $y->{_e}->copy();
 495     $x->{_m} = $y->{_m}->copy();
 496     $x->{sign} = $y->{sign} || $nan;
 497     return $x->round($a,$p,$r,$y);
 498     }
 499
 500   # take lower of the two e's and adapt m1 to it to match m2
 501   my $e = $y->{_e};
 502   $e = $MBI->bzero() if !defined $e;    # if no BFLOAT ?
 503   $e = $e->copy();                      # make copy (didn't do it yet)
 504   $e->bsub($x->{_e});
 505   my $add = $y->{_m}->copy();
 506   if ($e->{sign} eq '-')                # < 0
 507     {
 508     my $e1 = $e->copy()->babs();
 509     #$x->{_m} *= (10 ** $e1);
 510     $x->{_m}->blsft($e1,10);
 511     $x->{_e} += $e;                     # need the sign of e
 512     }
 513   elsif (!$e->is_zero())                # > 0
 514     {
 515     #$add *= (10 ** $e);
 516     $add->blsft($e,10);
 517     }
 518   # else: both e are the same, so just leave them
 519   $x->{_m}->{sign} = $x->{sign};                # fiddle with signs
 520   $add->{sign} = $y->{sign};
 521   $x->{_m} += $add;                             # finally do add/sub
 522   $x->{sign} = $x->{_m}->{sign};                # re-adjust signs
 523   $x->{_m}->{sign} = '+';                       # mantissa always positiv
 524   # delete trailing zeros, then round
 525   return $x->bnorm()->round($a,$p,$r,$y);
 526   }
 527
 528 sub bsub
 529   {
 530   # (BigFloat or num_str, BigFloat or num_str) return BigFloat
 531   # subtract second arg from first, modify first
 532
 533   # set up parameters
 534   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
 535   # objectify is costly, so avoid it
 536   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 537     {
 538     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 539     }
 540
 541   if ($y->is_zero())            # still round for not adding zero
 542     {
 543     return $x->round($a,$p,$r);
 544     }
 545
 546   $y->{sign} =~ tr/+\-/-+/;     # does nothing for NaN
 547   $x->badd($y,$a,$p,$r);        # badd does not leave internal zeros
 548   $y->{sign} =~ tr/+\-/-+/;     # refix $y (does nothing for NaN)
 549   $x;                           # already rounded by badd()
 550   }
 551
 552 sub binc
 553   {
 554   # increment arg by one
 555   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 556
 557   if ($x->{_e}->sign() eq '-')
 558     {
 559     return $x->badd($self->bone(),$a,$p,$r);    #  digits after dot
 560     }
 561
 562   if (!$x->{_e}->is_zero())
 563     {
 564     $x->{_m}->blsft($x->{_e},10);               # 1e2 => 100
 565     $x->{_e}->bzero();
 566     }
 567   # now $x->{_e} == 0
 568   if ($x->{sign} eq '+')
 569     {
 570     $x->{_m}->binc();
 571     return $x->bnorm()->bround($a,$p,$r);
 572     }
 573   elsif ($x->{sign} eq '-')
 574     {
 575     $x->{_m}->bdec();
 576     $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
 577     return $x->bnorm()->bround($a,$p,$r);
 578     }
 579   # inf, nan handling etc
 580   $x->badd($self->__one(),$a,$p,$r);            # does round
 581   }
 582
 583 sub bdec
 584   {
 585   # decrement arg by one
 586   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 587
 588   if ($x->{_e}->sign() eq '-')
 589     {
 590     return $x->badd($self->bone('-'),$a,$p,$r); #  digits after dot
 591     }
 592
 593   if (!$x->{_e}->is_zero())
 594     {
 595     $x->{_m}->blsft($x->{_e},10);               # 1e2 => 100
 596     $x->{_e}->bzero();
 597     }
 598   # now $x->{_e} == 0
 599   my $zero = $x->is_zero();
 600   # <= 0
 601   if (($x->{sign} eq '-') || $zero)
 602     {
 603     $x->{_m}->binc();
 604     $x->{sign} = '-' if $zero;                  # 0 => 1 => -1
 605     $x->{sign} = '+' if $x->{_m}->is_zero();    # -1 +1 => -0 => +0
 606     return $x->bnorm()->round($a,$p,$r);
 607     }
 608   # > 0
 609   elsif ($x->{sign} eq '+')
 610     {
 611     $x->{_m}->bdec();
 612     return $x->bnorm()->round($a,$p,$r);
 613     }
 614   # inf, nan handling etc
 615   $x->badd($self->bone('-'),$a,$p,$r);          # does round
 616   }
 617
 618 sub blog
 619   {
 620   my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(2,@_);
 621
 622   # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
 623
 624   # u = x-1, v = x+1
 625   #              _                               _
 626   # Taylor:     |    u    1   u^3   1   u^5       |
 627   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 0
 628   #             |_   v    3   v^3   5   v^5      _|
 629
 630   # This takes much more steps to calculate the result:
 631   # u = x-1
 632   #              _                               _
 633   # Taylor:     |    u    1   u^2   1   u^3       |
 634   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 1/2
 635   #             |_   x    2   x^2   3   x^3      _|
 636
 637   # we need to limit the accuracy to protect against overflow
 638   my $fallback = 0;
 639   my $scale = 0;
 640   my @params = $x->_find_round_parameters($a,$p,$r);
 641
 642   # no rounding at all, so must use fallback
 643   if (scalar @params == 1)
 644     {
 645     # simulate old behaviour
 646     $params[1] = $self->div_scale();    # and round to it as accuracy
 647     $params[0] = undef;
 648     $scale = $params[1]+4;              # at least four more for proper round
 649     $params[3] = $r;                    # round mode by caller or undef
 650     $fallback = 1;                      # to clear a/p afterwards
 651     }
 652   else
 653     {
 654     # the 4 below is empirical, and there might be cases where it is not
 655     # enough...
 656     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
 657     }
 658
 659   return $x->bzero(@params) if $x->is_one();
 660   return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
 661   return $x->bone('+',@params) if $x->bcmp($base) == 0;
 662
 663   # when user set globals, they would interfere with our calculation, so
 664   # disable them and later re-enable them
 665   no strict 'refs';
 666   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
 667   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
 668   # we also need to disable any set A or P on $x (_find_round_parameters took
 669   # them already into account), since these would interfere, too
 670   delete $x->{_a}; delete $x->{_p};
 671   # need to disable $upgrade in BigInt, to avoid deep recursion
 672   local $Math::BigInt::upgrade = undef;
 673
 674   my ($case,$limit,$v,$u,$below,$factor,$two,$next,$over,$f);
 675
 676   if (3 < 5)
 677   #if ($x <= Math::BigFloat->new("0.5"))
 678     {
 679     $case = 0;
 680   #  print "case $case $x < 0.5\n";
 681     $v = $x->copy(); $v->binc();                # v = x+1
 682     $x->bdec(); $u = $x->copy();                # u = x-1; x = x-1
 683     $x->bdiv($v,$scale);                        # first term: u/v
 684     $below = $v->copy();
 685     $over = $u->copy();
 686     $u *= $u; $v *= $v;                         # u^2, v^2
 687     $below->bmul($v);                           # u^3, v^3
 688     $over->bmul($u);
 689     $factor = $self->new(3); $f = $self->new(2);
 690     }
 691   #else
 692   #  {
 693   #  $case = 1;
 694   #  print "case 1 $x > 0.5\n";
 695   #  $v = $x->copy();                           # v = x
 696   #  $u = $x->copy(); $u->bdec();               # u = x-1;
 697   #  $x->bdec(); $x->bdiv($v,$scale);           # first term: x-1/x
 698   #  $below = $v->copy();
 699   #  $over = $u->copy();
 700   #  $below->bmul($v);                          # u^2, v^2
 701   #  $over->bmul($u);
 702   #  $factor = $self->new(2); $f = $self->bone();
 703   #  }
 704   $limit = $self->new("1E-". ($scale-1));
 705   #my $steps = 0;
 706   while (3 < 5)
 707     {
 708     # we calculate the next term, and add it to the last
 709     # when the next term is below our limit, it won't affect the outcome
 710     # anymore, so we stop
 711     $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
 712     last if $next->bcmp($limit) <= 0;
 713     $x->badd($next);
 714     # print "step  $x\n";
 715     # calculate things for the next term
 716     $over *= $u; $below *= $v; $factor->badd($f);
 717     #$steps++;
 718     }
 719   $x->bmul(2) if $case == 0;
 720   #print "took $steps steps\n";
 721
 722   # shortcut to not run trough _find_round_parameters again
 723   if (defined $params[1])
 724     {
 725     $x->bround($params[1],$params[3]);          # then round accordingly
 726     }
 727   else
 728     {
 729     $x->bfround($params[2],$params[3]);         # then round accordingly
 730     }
 731   if ($fallback)
 732     {
 733     # clear a/p after round, since user did not request it
 734     $x->{_a} = undef; $x->{_p} = undef;
 735     }
 736   # restore globals
 737   $$abr = $ab; $$pbr = $pb;
 738
 739   $x;
 740   }
 741
 742 sub blcm
 743   {
 744   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
 745   # does not modify arguments, but returns new object
 746   # Lowest Common Multiplicator
 747
 748   my ($self,@arg) = objectify(0,@_);
 749   my $x = $self->new(shift @arg);
 750   while (@arg) { $x = _lcm($x,shift @arg); }
 751   $x;
 752   }
 753
 754 sub bgcd
 755   {
 756   # (BFLOAT or num_str, BFLOAT or num_str) return BINT
 757   # does not modify arguments, but returns new object
 758   # GCD -- Euclids algorithm Knuth Vol 2 pg 296
 759
 760   my ($self,@arg) = objectify(0,@_);
 761   my $x = $self->new(shift @arg);
 762   while (@arg) { $x = _gcd($x,shift @arg); }
 763   $x;
 764   }
 765
 766 ###############################################################################
 767 # is_foo methods (is_negative, is_positive are inherited from BigInt)
 768
 769 sub is_int
 770   {
 771   # return true if arg (BFLOAT or num_str) is an integer
 772   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 773
 774   return 1 if ($x->{sign} =~ /^[+-]$/) &&       # NaN and +-inf aren't
 775     $x->{_e}->{sign} eq '+';                    # 1e-1 => no integer
 776   0;
 777   }
 778
 779 sub is_zero
 780   {
 781   # return true if arg (BFLOAT or num_str) is zero
 782   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 783
 784   return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
 785   0;
 786   }
 787
 788 sub is_one
 789   {
 790   # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
 791   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 792
 793   my $sign = shift || ''; $sign = '+' if $sign ne '-';
 794   return 1
 795    if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
 796   0;
 797   }
 798
 799 sub is_odd
 800   {
 801   # return true if arg (BFLOAT or num_str) is odd or false if even
 802   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 803
 804   return 1 if ($x->{sign} =~ /^[+-]$/) &&               # NaN & +-inf aren't
 805     ($x->{_e}->is_zero() && $x->{_m}->is_odd());
 806   0;
 807   }
 808
 809 sub is_even
 810   {
 811   # return true if arg (BINT or num_str) is even or false if odd
 812   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 813
 814   return 0 if $x->{sign} !~ /^[+-]$/;                   # NaN & +-inf aren't
 815   return 1 if ($x->{_e}->{sign} eq '+'                  # 123.45 is never
 816      && $x->{_m}->is_even());                           # but 1200 is
 817   0;
 818   }
 819
 820 sub bmul
 821   {
 822   # multiply two numbers -- stolen from Knuth Vol 2 pg 233
 823   # (BINT or num_str, BINT or num_str) return BINT
 824
 825   # set up parameters
 826   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
 827   # objectify is costly, so avoid it
 828   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 829     {
 830     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 831     }
 832
 833   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 834
 835   # inf handling
 836   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
 837     {
 838     return $x->bnan() if $x->is_zero() || $y->is_zero();
 839     # result will always be +-inf:
 840     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
 841     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
 842     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
 843     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
 844     return $x->binf('-');
 845     }
 846   # handle result = 0
 847   return $x->bzero() if $x->is_zero() || $y->is_zero();
 848
 849   return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
 850    ((!$x->isa($self)) || (!$y->isa($self)));
 851
 852   # aEb * cEd = (a*c)E(b+d)
 853   $x->{_m}->bmul($y->{_m});
 854   $x->{_e}->badd($y->{_e});
 855   # adjust sign:
 856   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
 857   return $x->bnorm()->round($a,$p,$r,$y);
 858   }
 859
 860 sub bdiv
 861   {
 862   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
 863   # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
 864
 865   # set up parameters
 866   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
 867   # objectify is costly, so avoid it
 868   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 869     {
 870     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 871     }
 872
 873   return $self->_div_inf($x,$y)
 874    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
 875
 876   # x== 0 # also: or y == 1 or y == -1
 877   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
 878
 879   # upgrade ?
 880   return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
 881
 882   # we need to limit the accuracy to protect against overflow
 883   my $fallback = 0;
 884   my $scale = 0;
 885   my @params = $x->_find_round_parameters($a,$p,$r,$y);
 886
 887   # no rounding at all, so must use fallback
 888   if (scalar @params == 1)
 889     {
 890     # simulate old behaviour
 891     $params[1] = $self->div_scale();    # and round to it as accuracy
 892     $scale = $params[1]+4;              # at least four more for proper round
 893     $params[3] = $r;                    # round mode by caller or undef
 894     $fallback = 1;                      # to clear a/p afterwards
 895     }
 896   else
 897     {
 898     # the 4 below is empirical, and there might be cases where it is not
 899     # enough...
 900     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
 901     }
 902   my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
 903   $scale = $lx if $lx > $scale;
 904   $scale = $ly if $ly > $scale;
 905   my $diff = $ly - $lx;
 906   $scale += $diff if $diff > 0;         # if lx << ly, but not if ly << lx!
 907
 908   # make copy of $x in case of list context for later reminder calculation
 909   my $rem;
 910   if (wantarray && !$y->is_one())
 911     {
 912     $rem = $x->copy();
 913     }
 914
 915   $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
 916
 917   # check for / +-1 ( +/- 1E0)
 918   if (!$y->is_one())
 919     {
 920     # promote BigInts and it's subclasses (except when already a BigFloat)
 921     $y = $self->new($y) unless $y->isa('Math::BigFloat');
 922
 923     # need to disable $upgrade in BigInt, to avoid deep recursion
 924     local $Math::BigInt::upgrade = undef;       # should be parent class vs MBI
 925
 926     # calculate the result to $scale digits and then round it
 927     # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
 928     $x->{_m}->blsft($scale,10);
 929     $x->{_m}->bdiv( $y->{_m} ); # a/c
 930     $x->{_e}->bsub( $y->{_e} ); # b-d
 931     $x->{_e}->bsub($scale);     # correct for 10**scale
 932     $x->bnorm();                # remove trailing 0's
 933     }
 934
 935   # shortcut to not run trough _find_round_parameters again
 936   if (defined $params[1])
 937     {
 938     $x->{_a} = undef;                           # clear before round
 939     $x->bround($params[1],$params[3]);          # then round accordingly
 940     }
 941   else
 942     {
 943     $x->{_p} = undef;                           # clear before round
 944     $x->bfround($params[2],$params[3]);         # then round accordingly
 945     }
 946   if ($fallback)
 947     {
 948     # clear a/p after round, since user did not request it
 949     $x->{_a} = undef; $x->{_p} = undef;
 950     }
 951
 952   if (wantarray)
 953     {
 954     if (!$y->is_one())
 955       {
 956       $rem->bmod($y,$params[1],$params[2],$params[3]);  # copy already done
 957       }
 958     else
 959       {
 960       $rem = $self->bzero();
 961       }
 962     if ($fallback)
 963       {
 964       # clear a/p after round, since user did not request it
 965       $rem->{_a} = undef; $rem->{_p} = undef;
 966       }
 967     return ($x,$rem);
 968     }
 969   $x;
 970   }
 971
 972 sub bmod
 973   {
 974   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
 975
 976   # set up parameters
 977   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
 978   # objectify is costly, so avoid it
 979   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
 980     {
 981     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 982     }
 983
 984   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
 985     {
 986     my ($d,$re) = $self->SUPER::_div_inf($x,$y);
 987     $x->{sign} = $re->{sign};
 988     $x->{_e} = $re->{_e};
 989     $x->{_m} = $re->{_m};
 990     return $x->round($a,$p,$r,$y);
 991     }
 992   return $x->bnan() if $x->is_zero() && $y->is_zero();
 993   return $x if $y->is_zero();
 994   return $x->bnan() if $x->is_nan() || $y->is_nan();
 995   return $x->bzero() if $y->is_one() || $x->is_zero();
 996
 997   # inf handling is missing here
 998
 999   my $cmp = $x->bacmp($y);                      # equal or $x < $y?
1000   return $x->bzero($a,$p) if $cmp == 0;         # $x == $y => result 0
1001
1002   # only $y of the operands negative?
1003   my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1004
1005   $x->{sign} = $y->{sign};                              # calc sign first
1006   return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0;  # $x < $y => result $x
1007
1008   my $ym = $y->{_m}->copy();
1009
1010   # 2e1 => 20
1011   $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
1012
1013   # if $y has digits after dot
1014   my $shifty = 0;                       # correct _e of $x by this
1015   if ($y->{_e}->{sign} eq '-')          # has digits after dot
1016     {
1017     # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1018     $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
1019     $x->blsft($shifty,10);              # 123 => 1230, $y->{_m} is already 25
1020     }
1021   # $ym is now mantissa of $y based on exponent 0
1022
1023   my $shiftx = 0;                       # correct _e of $x by this
1024   if ($x->{_e}->{sign} eq '-')          # has digits after dot
1025     {
1026     # 123.4 % 20 => 1234 % 200
1027     $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
1028     $ym->blsft($shiftx,10);
1029     }
1030   # 123e1 % 20 => 1230 % 20
1031   if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
1032     {
1033     $x->{_m}->blsft($x->{_e},10);
1034     }
1035   $x->{_e} = $MBI->bzero() unless $x->{_e}->is_zero();
1036
1037   $x->{_e}->bsub($shiftx) if $shiftx != 0;
1038   $x->{_e}->bsub($shifty) if $shifty != 0;
1039
1040   # now mantissas are equalized, exponent of $x is adjusted, so calc result
1041
1042   $x->{_m}->bmod($ym);
1043
1044   $x->{sign} = '+' if $x->{_m}->is_zero();              # fix sign for -0
1045   $x->bnorm();
1046
1047   if ($neg != 0)        # one of them negative => correct in place
1048     {
1049     my $r = $y - $x;
1050     $x->{_m} = $r->{_m};
1051     $x->{_e} = $r->{_e};
1052     $x->{sign} = '+' if $x->{_m}->is_zero();            # fix sign for -0
1053     $x->bnorm();
1054     }
1055
1056   $x->round($a,$p,$r,$y);       # round and return
1057   }
1058
1059 sub bsqrt
1060   {
1061   # calculate square root; this should probably
1062   # use a different test to see whether the accuracy we want is...
1063   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1064
1065   return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
1066   return $x if $x->{sign} eq '+inf';                              # +inf
1067   return $x if $x->is_zero() || $x->is_one();
1068
1069   # we need to limit the accuracy to protect against overflow
1070   my $fallback = 0;
1071   my $scale = 0;
1072   my @params = $x->_find_round_parameters($a,$p,$r);
1073
1074   # no rounding at all, so must use fallback
1075   if ((scalar @params == 1) ||
1076       (!defined($params[1] || $params[2])))
1077     {
1078     # simulate old behaviour
1079     $params[1] = $self->div_scale();    # and round to it as accuracy
1080     $scale = $params[1]+4;              # at least four more for proper round
1081     $params[3] = $r;                    # round mode by caller or undef
1082     $fallback = 1;                      # to clear a/p afterwards
1083     }
1084   else
1085     {
1086     # the 4 below is empirical, and there might be cases where it is not
1087     # enough...
1088     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1089     }
1090
1091   # when user set globals, they would interfere with our calculation, so
1092   # disable them and later re-enable them
1093   no strict 'refs';
1094   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1095   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1096   # we also need to disable any set A or P on $x (_find_round_parameters took
1097   # them already into account), since these would interfere, too
1098   delete $x->{_a}; delete $x->{_p};
1099   # need to disable $upgrade in BigInt, to avoid deep recursion
1100   local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1101
1102   my $xas = $x->as_number();
1103   my $gs = $xas->copy()->bsqrt();       # some guess
1104
1105   if (($x->{_e}->{sign} ne '-')         # guess can't be accurate if there are
1106                                         # digits after the dot
1107    && ($xas->bacmp($gs * $gs) == 0))    # guess hit the nail on the head?
1108     {
1109     # exact result
1110     $x->{_m} = $gs; $x->{_e} = $MBI->bzero(); $x->bnorm();
1111     # shortcut to not run trough _find_round_parameters again
1112     if (defined $params[1])
1113       {
1114       $x->bround($params[1],$params[3]);        # then round accordingly
1115       }
1116     else
1117       {
1118       $x->bfround($params[2],$params[3]);       # then round accordingly
1119       }
1120     if ($fallback)
1121       {
1122       # clear a/p after round, since user did not request it
1123       $x->{_a} = undef; $x->{_p} = undef;
1124       }
1125     # re-enable A and P, upgrade is taken care of by "local"
1126     ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1127     return $x;
1128     }
1129
1130   # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
1131   # of the result by multipyling the input by 100 and then divide the integer
1132   # result of sqrt(input) by 10. Rounding afterwards returns the real result.
1133   # this will transform 123.456 (in $x) into 123456 (in $y1)
1134   my $y1 = $x->{_m}->copy();
1135   # We now make sure that $y1 has the same odd or even number of digits than
1136   # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
1137   # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
1138   # steps of 10. The length of $x does not count, since an even or odd number
1139   # of digits before the dot is not changed by adding an even number of digits
1140   # after the dot (the result is still odd or even digits long).
1141   $y1->bmul(10) if $x->{_e}->is_odd();
1142   # now calculate how many digits the result of sqrt(y1) would have
1143   my $digits = int($y1->length() / 2);
1144   # but we need at least $scale digits, so calculate how many are missing
1145   my $shift = $scale - $digits;
1146   # that should never happen (we take care of integer guesses above)
1147   # $shift = 0 if $shift < 0;
1148   # multiply in steps of 100, by shifting left two times the "missing" digits
1149   $y1->blsft($shift*2,10);
1150   # now take the square root and truncate to integer
1151   $y1->bsqrt();
1152   # By "shifting" $y1 right (by creating a negative _e) we calculate the final
1153   # result, which is than later rounded to the desired scale.
1154   $x->{_m} = $y1;
1155   # gs->length() is the number of digits before the dot. Since gs is always
1156   # truncated (9.99 => 9), it is always right (if gs was rounded, it would be
1157   # '10' and thus gs->length() == 2, which would be wrong).
1158   $x->{_e} = $MBI->new(- $y1->length() + $gs->length());
1159
1160   # shortcut to not run trough _find_round_parameters again
1161   if (defined $params[1])
1162     {
1163     $x->bround($params[1],$params[3]);          # then round accordingly
1164     }
1165   else
1166     {
1167     $x->bfround($params[2],$params[3]);         # then round accordingly
1168     }
1169   if ($fallback)
1170     {
1171     # clear a/p after round, since user did not request it
1172     $x->{_a} = undef; $x->{_p} = undef;
1173     }
1174   # restore globals
1175   $$abr = $ab; $$pbr = $pb;
1176   $x;
1177   }
1178
1179 sub bfac
1180   {
1181   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1182   # compute factorial numbers
1183   # modifies first argument
1184   my ($self,$x,@r) = objectify(1,@_);
1185
1186   return $x->bnan()
1187     if (($x->{sign} ne '+') ||          # inf, NaN, <0 etc => NaN
1188      ($x->{_e}->{sign} ne '+'));        # digits after dot?
1189
1190   return $x->bone('+',@r) if $x->is_zero() || $x->is_one();     # 0 or 1 => 1
1191
1192   # use BigInt's bfac() for faster calc
1193   $x->{_m}->blsft($x->{_e},10);         # un-norm m
1194   $x->{_e}->bzero();                    # norm $x again
1195   $x->{_m}->bfac();                     # factorial
1196   $x->bnorm()->round(@r);
1197   }
1198
1199 sub _pow2
1200   {
1201   # Calculate a power where $y is a non-integer, like 2 ** 0.5
1202   my ($x,$y,$a,$p,$r) = @_;
1203   my $self = ref($x);
1204
1205   # we need to limit the accuracy to protect against overflow
1206   my $fallback = 0;
1207   my $scale = 0;
1208   my @params = $x->_find_round_parameters($a,$p,$r);
1209
1210   # no rounding at all, so must use fallback
1211   if (scalar @params == 1)
1212     {
1213     # simulate old behaviour
1214     $params[1] = $self->div_scale();    # and round to it as accuracy
1215     $scale = $params[1]+4;              # at least four more for proper round
1216     $params[3] = $r;                    # round mode by caller or undef
1217     $fallback = 1;                      # to clear a/p afterwards
1218     }
1219   else
1220     {
1221     # the 4 below is empirical, and there might be cases where it is not
1222     # enough...
1223     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1224     }
1225
1226   # when user set globals, they would interfere with our calculation, so
1227   # disable them and later re-enable them
1228   no strict 'refs';
1229   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1230   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1231   # we also need to disable any set A or P on $x (_find_round_parameters took
1232   # them already into account), since these would interfere, too
1233   delete $x->{_a}; delete $x->{_p};
1234   # need to disable $upgrade in BigInt, to avoid deep recursion
1235   local $Math::BigInt::upgrade = undef;
1236
1237   # split the second argument into its integer and fraction part
1238   # we calculate the result then from these two parts, like in
1239   # 2 ** 2.4 == (2 ** 2) * (2 ** 0.4)
1240   my $c = $self->new($y->as_number());  # integer part
1241   my $d = $y-$c;                        # fractional part
1242   my $xc = $x->copy();                  # a temp. copy
1243
1244   # now calculate binary fraction from the decimal fraction on the fly
1245   # f.i. 0.654:
1246   # 0.654 * 2 = 1.308 > 1 => 0.1        ( 1.308 - 1 = 0.308)
1247   # 0.308 * 2 = 0.616 < 1 => 0.10
1248   # 0.616 * 2 = 1.232 > 1 => 0.101      ( 1.232 - 1 = 0.232)
1249   # and so on...
1250   # The process stops when the result is exactly one, or when we have
1251   # enough accuracy
1252
1253   # From the binary fraction we calculate the result as follows:
1254   # we assume the fraction ends in 1, and we remove this one first.
1255   # For each digit after the dot, assume 1 eq R and 0 eq XR, where R means
1256   # take square root and X multiply with the original X.
1257
1258   my $i = 0;
1259   while ($i++ < 50)
1260     {
1261     $d->badd($d);                                               # * 2
1262     last if $d->is_one();                                       # == 1
1263     $x->bsqrt();                                                # 0
1264     if ($d > 1)
1265       {
1266       $x->bsqrt(); $x->bmul($xc); $d->bdec();                   # 1
1267       }
1268     }
1269   # assume fraction ends in 1
1270   $x->bsqrt();                                                  # 1
1271   if (!$c->is_one())
1272     {
1273     $x->bmul( $xc->bpow($c) );
1274     }
1275   elsif (!$c->is_zero())
1276     {
1277     $x->bmul( $xc );
1278     }
1279   # done
1280
1281   # shortcut to not run trough _find_round_parameters again
1282   if (defined $params[1])
1283     {
1284     $x->bround($params[1],$params[3]);          # then round accordingly
1285     }
1286   else
1287     {
1288     $x->bfround($params[2],$params[3]);         # then round accordingly
1289     }
1290   if ($fallback)
1291     {
1292     # clear a/p after round, since user did not request it
1293     $x->{_a} = undef; $x->{_p} = undef;
1294     }
1295   # restore globals
1296   $$abr = $ab; $$pbr = $pb;
1297   $x;
1298   }
1299
1300 sub _pow
1301   {
1302   # Calculate a power where $y is a non-integer, like 2 ** 0.5
1303   my ($x,$y,$a,$p,$r) = @_;
1304   my $self = ref($x);
1305
1306   # if $y == 0.5, it is sqrt($x)
1307   return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
1308
1309   # u = y * ln x
1310   #                _                             _
1311   # Taylor:       |    u     u^2      u^3         |
1312   # x ** y  = 1 + |   --- +  --- + * ----- + ...  |
1313   #               |_   1     1*2     1*2*3       _|
1314
1315   # we need to limit the accuracy to protect against overflow
1316   my $fallback = 0;
1317   my $scale = 0;
1318   my @params = $x->_find_round_parameters($a,$p,$r);
1319
1320   # no rounding at all, so must use fallback
1321   if (scalar @params == 1)
1322     {
1323     # simulate old behaviour
1324     $params[1] = $self->div_scale();    # and round to it as accuracy
1325     $scale = $params[1]+4;              # at least four more for proper round
1326     $params[3] = $r;                    # round mode by caller or undef
1327     $fallback = 1;                      # to clear a/p afterwards
1328     }
1329   else
1330     {
1331     # the 4 below is empirical, and there might be cases where it is not
1332     # enough...
1333     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1334     }
1335
1336   # when user set globals, they would interfere with our calculation, so
1337   # disable them and later re-enable them
1338   no strict 'refs';
1339   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1340   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1341   # we also need to disable any set A or P on $x (_find_round_parameters took
1342   # them already into account), since these would interfere, too
1343   delete $x->{_a}; delete $x->{_p};
1344   # need to disable $upgrade in BigInt, to avoid deep recursion
1345   local $Math::BigInt::upgrade = undef;
1346
1347   my ($limit,$v,$u,$below,$factor,$next,$over);
1348
1349   $u = $x->copy()->blog($scale)->bmul($y);
1350   $v = $self->bone();                           # 1
1351   $factor = $self->new(2);                      # 2
1352   $x->bone();                                   # first term: 1
1353
1354   $below = $v->copy();
1355   $over = $u->copy();
1356
1357   $limit = $self->new("1E-". ($scale-1));
1358   #my $steps = 0;
1359   while (3 < 5)
1360     {
1361     # we calculate the next term, and add it to the last
1362     # when the next term is below our limit, it won't affect the outcome
1363     # anymore, so we stop
1364     $next = $over->copy()->bdiv($below,$scale);
1365     last if $next->bcmp($limit) <= 0;
1366     $x->badd($next);
1367 #    print "at $x\n";
1368     # calculate things for the next term
1369     $over *= $u; $below *= $factor; $factor->binc();
1370     #$steps++;
1371     }
1372
1373   # shortcut to not run trough _find_round_parameters again
1374   if (defined $params[1])
1375     {
1376     $x->bround($params[1],$params[3]);          # then round accordingly
1377     }
1378   else
1379     {
1380     $x->bfround($params[2],$params[3]);         # then round accordingly
1381     }
1382   if ($fallback)
1383     {
1384     # clear a/p after round, since user did not request it
1385     $x->{_a} = undef; $x->{_p} = undef;
1386     }
1387   # restore globals
1388   $$abr = $ab; $$pbr = $pb;
1389   $x;
1390   }
1391
1392 sub bpow
1393   {
1394   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1395   # compute power of two numbers, second arg is used as integer
1396   # modifies first argument
1397
1398   # set up parameters
1399   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1400   # objectify is costly, so avoid it
1401   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1402     {
1403     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1404     }
1405
1406   return $x if $x->{sign} =~ /^[+-]inf$/;
1407   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1408   return $x->bone() if $y->is_zero();
1409   return $x         if $x->is_one() || $y->is_one();
1410
1411   return $x->_pow($y,$a,$p,$r) if !$y->is_int();        # non-integer power
1412
1413   my $y1 = $y->as_number();             # make bigint
1414   # if ($x == -1)
1415   if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
1416     {
1417     # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
1418     return $y1->is_odd() ? $x : $x->babs(1);
1419     }
1420   if ($x->is_zero())
1421     {
1422     return $x if $y->{sign} eq '+';     # 0**y => 0 (if not y <= 0)
1423     # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1424     $x->binf();
1425     }
1426
1427   # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1428   $y1->babs();
1429   $x->{_m}->bpow($y1);
1430   $x->{_e}->bmul($y1);
1431   $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1432   $x->bnorm();
1433   if ($y->{sign} eq '-')
1434     {
1435     # modify $x in place!
1436     my $z = $x->copy(); $x->bzero()->binc();
1437     return $x->bdiv($z,$a,$p,$r);       # round in one go (might ignore y's A!)
1438     }
1439   $x->round($a,$p,$r,$y);
1440   }
1441
1442 ###############################################################################
1443 # rounding functions
1444
1445 sub bfround
1446   {
1447   # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1448   # $n == 0 means round to integer
1449   # expects and returns normalized numbers!
1450   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1451
1452   return $x if $x->modify('bfround');
1453
1454   my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1455   return $x if !defined $scale;                 # no-op
1456
1457   # never round a 0, +-inf, NaN
1458   if ($x->is_zero())
1459     {
1460     $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1461     return $x;
1462     }
1463   return $x if $x->{sign} !~ /^[+-]$/;
1464
1465   # don't round if x already has lower precision
1466   return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1467
1468   $x->{_p} = $scale;                    # remember round in any case
1469   $x->{_a} = undef;                     # and clear A
1470   if ($scale < 0)
1471     {
1472     # round right from the '.'
1473
1474     return $x if $x->{_e}->{sign} eq '+';       # e >= 0 => nothing to round
1475
1476     $scale = -$scale;                           # positive for simplicity
1477     my $len = $x->{_m}->length();               # length of mantissa
1478
1479     # the following poses a restriction on _e, but if _e is bigger than a
1480     # scalar, you got other problems (memory etc) anyway
1481     my $dad = -($x->{_e}->numify());            # digits after dot
1482     my $zad = 0;                                # zeros after dot
1483     $zad = $dad - $len if (-$dad < -$len);      # for 0.00..00xxx style
1484
1485     #print "scale $scale dad $dad zad $zad len $len\n";
1486     # number  bsstr   len zad dad
1487     # 0.123   123e-3    3   0 3
1488     # 0.0123  123e-4    3   1 4
1489     # 0.001   1e-3      1   2 3
1490     # 1.23    123e-2    3   0 2
1491     # 1.2345  12345e-4  5   0 4
1492
1493     # do not round after/right of the $dad
1494     return $x if $scale > $dad;                 # 0.123, scale >= 3 => exit
1495
1496     # round to zero if rounding inside the $zad, but not for last zero like:
1497     # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1498     return $x->bzero() if $scale < $zad;
1499     if ($scale == $zad)                 # for 0.006, scale -3 and trunc
1500       {
1501       $scale = -$len;
1502       }
1503     else
1504       {
1505       # adjust round-point to be inside mantissa
1506       if ($zad != 0)
1507         {
1508         $scale = $scale-$zad;
1509         }
1510       else
1511         {
1512         my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;    # digits before dot
1513         $scale = $dbd+$scale;
1514         }
1515       }
1516     }
1517   else
1518     {
1519     # round left from the '.'
1520
1521     # 123 => 100 means length(123) = 3 - $scale (2) => 1
1522
1523     my $dbt = $x->{_m}->length();
1524     # digits before dot
1525     my $dbd = $dbt + $x->{_e}->numify();
1526     # should be the same, so treat it as this
1527     $scale = 1 if $scale == 0;
1528     # shortcut if already integer
1529     return $x if $scale == 1 && $dbt <= $dbd;
1530     # maximum digits before dot
1531     ++$dbd;
1532
1533     if ($scale > $dbd)
1534        {
1535        # not enough digits before dot, so round to zero
1536        return $x->bzero;
1537        }
1538     elsif ( $scale == $dbd )
1539        {
1540        # maximum
1541        $scale = -$dbt;
1542        }
1543     else
1544        {
1545        $scale = $dbd - $scale;
1546        }
1547     }
1548   # pass sign to bround for rounding modes '+inf' and '-inf'
1549   $x->{_m}->{sign} = $x->{sign};
1550   $x->{_m}->bround($scale,$mode);
1551   $x->{_m}->{sign} = '+';               # fix sign back
1552   $x->bnorm();
1553   }
1554
1555 sub bround
1556   {
1557   # accuracy: preserve $N digits, and overwrite the rest with 0's
1558   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1559
1560   die ('bround() needs positive accuracy') if ($_[0] || 0) < 0;
1561
1562   my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1563   return $x if !defined $scale;                         # no-op
1564
1565   return $x if $x->modify('bround');
1566
1567   # scale is now either $x->{_a}, $accuracy, or the user parameter
1568   # test whether $x already has lower accuracy, do nothing in this case
1569   # but do round if the accuracy is the same, since a math operation might
1570   # want to round a number with A=5 to 5 digits afterwards again
1571   return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
1572
1573   # scale < 0 makes no sense
1574   # never round a +-inf, NaN
1575   return $x if ($scale < 0) ||  $x->{sign} !~ /^[+-]$/;
1576
1577   # 1: $scale == 0 => keep all digits
1578   # 2: never round a 0
1579   # 3: if we should keep more digits than the mantissa has, do nothing
1580   if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1581     {
1582     $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1583     return $x;
1584     }
1585
1586   # pass sign to bround for '+inf' and '-inf' rounding modes
1587   $x->{_m}->{sign} = $x->{sign};
1588   $x->{_m}->bround($scale,$mode);       # round mantissa
1589   $x->{_m}->{sign} = '+';               # fix sign back
1590   # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
1591   $x->{_a} = $scale;                    # remember rounding
1592   $x->{_p} = undef;                     # and clear P
1593   $x->bnorm();                          # del trailing zeros gen. by bround()
1594   }
1595
1596 sub bfloor
1597   {
1598   # return integer less or equal then $x
1599   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1600
1601   return $x if $x->modify('bfloor');
1602
1603   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1604
1605   # if $x has digits after dot
1606   if ($x->{_e}->{sign} eq '-')
1607     {
1608     $x->{_e}->{sign} = '+';                     # negate e
1609     $x->{_m}->brsft($x->{_e},10);               # cut off digits after dot
1610     $x->{_e}->bzero();                          # trunc/norm
1611     $x->{_m}->binc() if $x->{sign} eq '-';      # decrement if negative
1612     }
1613   $x->round($a,$p,$r);
1614   }
1615
1616 sub bceil
1617   {
1618   # return integer greater or equal then $x
1619   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1620
1621   return $x if $x->modify('bceil');
1622   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1623
1624   # if $x has digits after dot
1625   if ($x->{_e}->{sign} eq '-')
1626     {
1627     #$x->{_m}->brsft(-$x->{_e},10);
1628     #$x->{_e}->bzero();
1629     #$x++ if $x->{sign} eq '+';
1630
1631     $x->{_e}->{sign} = '+';                     # negate e
1632     $x->{_m}->brsft($x->{_e},10);               # cut off digits after dot
1633     $x->{_e}->bzero();                          # trunc/norm
1634     $x->{_m}->binc() if $x->{sign} eq '+';      # decrement if negative
1635     }
1636   $x->round($a,$p,$r);
1637   }
1638
1639 sub brsft
1640   {
1641   # shift right by $y (divide by power of $n)
1642
1643   # set up parameters
1644   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1645   # objectify is costly, so avoid it
1646   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1647     {
1648     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1649     }
1650
1651   return $x if $x->modify('brsft');
1652   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1653
1654   $n = 2 if !defined $n; $n = $self->new($n);
1655   $x->bdiv($n->bpow($y),$a,$p,$r,$y);
1656   }
1657
1658 sub blsft
1659   {
1660   # shift left by $y (multiply by power of $n)
1661
1662   # set up parameters
1663   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1664   # objectify is costly, so avoid it
1665   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1666     {
1667     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1668     }
1669
1670   return $x if $x->modify('blsft');
1671   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1672
1673   $n = 2 if !defined $n; $n = $self->new($n);
1674   $x->bmul($n->bpow($y),$a,$p,$r,$y);
1675   }
1676
1677 ###############################################################################
1678
1679 sub DESTROY
1680   {
1681   # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
1682   }
1683
1684 sub AUTOLOAD
1685   {
1686   # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
1687   # or falling back to MBI::bxxx()
1688   my $name = $AUTOLOAD;
1689
1690   $name =~ s/.*:://;    # split package
1691   no strict 'refs';
1692   if (!method_alias($name))
1693     {
1694     if (!defined $name)
1695       {
1696       # delayed load of Carp and avoid recursion
1697       require Carp;
1698       Carp::croak ("Can't call a method without name");
1699       }
1700     if (!method_hand_up($name))
1701       {
1702       # delayed load of Carp and avoid recursion
1703       require Carp;
1704       Carp::croak ("Can't call $class\-\>$name, not a valid method");
1705       }
1706     # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1707     $name =~ s/^f/b/;
1708     return &{"$MBI"."::$name"}(@_);
1709     }
1710   my $bname = $name; $bname =~ s/^f/b/;
1711   *{$class."::$name"} = \&$bname;
1712   &$bname;      # uses @_
1713   }
1714
1715 sub exponent
1716   {
1717   # return a copy of the exponent
1718   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1719
1720   if ($x->{sign} !~ /^[+-]$/)
1721     {
1722     my $s = $x->{sign}; $s =~ s/^[+-]//;
1723     return $self->new($s);                      # -inf, +inf => +inf
1724     }
1725   return $x->{_e}->copy();
1726   }
1727
1728 sub mantissa
1729   {
1730   # return a copy of the mantissa
1731   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1732
1733   if ($x->{sign} !~ /^[+-]$/)
1734     {
1735     my $s = $x->{sign}; $s =~ s/^[+]//;
1736     return $self->new($s);                      # -inf, +inf => +inf
1737     }
1738   my $m = $x->{_m}->copy();             # faster than going via bstr()
1739   $m->bneg() if $x->{sign} eq '-';
1740
1741   $m;
1742   }
1743
1744 sub parts
1745   {
1746   # return a copy of both the exponent and the mantissa
1747   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1748
1749   if ($x->{sign} !~ /^[+-]$/)
1750     {
1751     my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
1752     return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
1753     }
1754   my $m = $x->{_m}->copy();     # faster than going via bstr()
1755   $m->bneg() if $x->{sign} eq '-';
1756   return ($m,$x->{_e}->copy());
1757   }
1758
1759 ##############################################################################
1760 # private stuff (internal use only)
1761
1762 sub import
1763   {
1764   my $self = shift;
1765   my $l = scalar @_;
1766   my $lib = ''; my @a;
1767   for ( my $i = 0; $i < $l ; $i++)
1768     {
1769     if ( $_[$i] eq ':constant' )
1770       {
1771       # this rest causes overlord er load to step in
1772       # print "overload @_\n";
1773       overload::constant float => sub { $self->new(shift); };
1774       }
1775     elsif ($_[$i] eq 'upgrade')
1776       {
1777       # this causes upgrading
1778       $upgrade = $_[$i+1];              # or undef to disable
1779       $i++;
1780       }
1781     elsif ($_[$i] eq 'downgrade')
1782       {
1783       # this causes downgrading
1784       $downgrade = $_[$i+1];            # or undef to disable
1785       $i++;
1786       }
1787     elsif ($_[$i] eq 'lib')
1788       {
1789       $lib = $_[$i+1] || '';            # default Calc
1790       $i++;
1791       }
1792     elsif ($_[$i] eq 'with')
1793       {
1794       $MBI = $_[$i+1] || 'Math::BigInt';        # default Math::BigInt
1795       $i++;
1796       }
1797     else
1798       {
1799       push @a, $_[$i];
1800       }
1801     }
1802
1803   # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
1804   my $mbilib = eval { Math::BigInt->config()->{lib} };
1805   if ((defined $mbilib) && ($MBI eq 'Math::BigInt'))
1806     {
1807     # MBI already loaded
1808     $MBI->import('lib',"$lib,$mbilib", 'objectify');
1809     }
1810   else
1811     {
1812     # MBI not loaded, or with ne "Math::BigInt"
1813     $lib .= ",$mbilib" if defined $mbilib;
1814     $lib =~ s/^,//;                             # don't leave empty
1815     if ($] < 5.006)
1816       {
1817       # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
1818       # used in the same script, or eval inside import().
1819       my @parts = split /::/, $MBI;             # Math::BigInt => Math BigInt
1820       my $file = pop @parts; $file .= '.pm';    # BigInt => BigInt.pm
1821       require File::Spec;
1822       $file = File::Spec->catfile (@parts, $file);
1823       eval { require "$file"; };
1824       $MBI->import( lib => $lib, 'objectify' );
1825       }
1826     else
1827       {
1828       my $rc = "use $MBI lib => '$lib', 'objectify';";
1829       eval $rc;
1830       }
1831     }
1832   die ("Couldn't load $MBI: $! $@") if $@;
1833
1834   # any non :constant stuff is handled by our parent, Exporter
1835   # even if @_ is empty, to give it a chance
1836   $self->SUPER::import(@a);             # for subclasses
1837   $self->export_to_level(1,$self,@a);   # need this, too
1838   }
1839
1840 sub bnorm
1841   {
1842   # adjust m and e so that m is smallest possible
1843   # round number according to accuracy and precision settings
1844   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1845
1846   return $x if $x->{sign} !~ /^[+-]$/;          # inf, nan etc
1847
1848 #  if (!$x->{_m}->is_odd())
1849 #    {
1850     my $zeros = $x->{_m}->_trailing_zeros();    # correct for trailing zeros
1851     if ($zeros != 0)
1852       {
1853       $x->{_m}->brsft($zeros,10); $x->{_e}->badd($zeros);
1854       }
1855     # for something like 0Ey, set y to 1, and -0 => +0
1856     $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
1857 #    }
1858   # this is to prevent automatically rounding when MBI's globals are set
1859   $x->{_m}->{_f} = MB_NEVER_ROUND;
1860   $x->{_e}->{_f} = MB_NEVER_ROUND;
1861   # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
1862   $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
1863   $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
1864   $x;                                   # MBI bnorm is no-op, so dont call it
1865   }
1866
1867 ##############################################################################
1868
1869 sub as_hex
1870   {
1871   # return number as hexadecimal string (only for integers defined)
1872   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1873
1874   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
1875   return '0x0' if $x->is_zero();
1876
1877   return 'NaN' if $x->{_e}->{sign} ne '+';      # how to do 1e-1 in hex!?
1878
1879   my $z = $x->{_m}->copy();
1880   if (!$x->{_e}->is_zero())             # > 0
1881     {
1882     $z->blsft($x->{_e},10);
1883     }
1884   $z->{sign} = $x->{sign};
1885   $z->as_hex();
1886   }
1887
1888 sub as_bin
1889   {
1890   # return number as binary digit string (only for integers defined)
1891   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1892
1893   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
1894   return '0b0' if $x->is_zero();
1895
1896   return 'NaN' if $x->{_e}->{sign} ne '+';      # how to do 1e-1 in hex!?
1897
1898   my $z = $x->{_m}->copy();
1899   if (!$x->{_e}->is_zero())             # > 0
1900     {
1901     $z->blsft($x->{_e},10);
1902     }
1903   $z->{sign} = $x->{sign};
1904   $z->as_bin();
1905   }
1906
1907 sub as_number
1908   {
1909   # return copy as a bigint representation of this BigFloat number
1910   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1911
1912   my $z = $x->{_m}->copy();
1913   if ($x->{_e}->{sign} eq '-')          # < 0
1914     {
1915     $x->{_e}->{sign} = '+';             # flip
1916     $z->brsft($x->{_e},10);
1917     $x->{_e}->{sign} = '-';             # flip back
1918     }
1919   elsif (!$x->{_e}->is_zero())          # > 0
1920     {
1921     $z->blsft($x->{_e},10);
1922     }
1923   $z->{sign} = $x->{sign};
1924   $z;
1925   }
1926
1927 sub length
1928   {
1929   my $x = shift;
1930   my $class = ref($x) || $x;
1931   $x = $class->new(shift) unless ref($x);
1932
1933   return 1 if $x->{_m}->is_zero();
1934   my $len = $x->{_m}->length();
1935   $len += $x->{_e} if $x->{_e}->sign() eq '+';
1936   if (wantarray())
1937     {
1938     my $t = $MBI->bzero();
1939     $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1940     return ($len,$t);
1941     }
1942   $len;
1943   }
1944
1945 1;
1946 __END__
1947
1948 =head1 NAME
1949
1950 Math::BigFloat - Arbitrary size floating point math package
1951
1952 =head1 SYNOPSIS
1953
1954   use Math::BigFloat;
1955
1956   # Number creation
1957   $x = Math::BigFloat->new($str);       # defaults to 0
1958   $nan  = Math::BigFloat->bnan();       # create a NotANumber
1959   $zero = Math::BigFloat->bzero();      # create a +0
1960   $inf = Math::BigFloat->binf();        # create a +inf
1961   $inf = Math::BigFloat->binf('-');     # create a -inf
1962   $one = Math::BigFloat->bone();        # create a +1
1963   $one = Math::BigFloat->bone('-');     # create a -1
1964
1965   # Testing
1966   $x->is_zero();                # true if arg is +0
1967   $x->is_nan();                 # true if arg is NaN
1968   $x->is_one();                 # true if arg is +1
1969   $x->is_one('-');              # true if arg is -1
1970   $x->is_odd();                 # true if odd, false for even
1971   $x->is_even();                # true if even, false for odd
1972   $x->is_positive();            # true if >= 0
1973   $x->is_negative();            # true if <  0
1974   $x->is_inf(sign);             # true if +inf, or -inf (default is '+')
1975
1976   $x->bcmp($y);                 # compare numbers (undef,<0,=0,>0)
1977   $x->bacmp($y);                # compare absolutely (undef,<0,=0,>0)
1978   $x->sign();                   # return the sign, either +,- or NaN
1979   $x->digit($n);                # return the nth digit, counting from right
1980   $x->digit(-$n);               # return the nth digit, counting from left
1981
1982   # The following all modify their first argument:
1983
1984   # set
1985   $x->bzero();                  # set $i to 0
1986   $x->bnan();                   # set $i to NaN
1987   $x->bone();                   # set $x to +1
1988   $x->bone('-');                # set $x to -1
1989   $x->binf();                   # set $x to inf
1990   $x->binf('-');                # set $x to -inf
1991
1992   $x->bneg();                   # negation
1993   $x->babs();                   # absolute value
1994   $x->bnorm();                  # normalize (no-op)
1995   $x->bnot();                   # two's complement (bit wise not)
1996   $x->binc();                   # increment x by 1
1997   $x->bdec();                   # decrement x by 1
1998
1999   $x->badd($y);                 # addition (add $y to $x)
2000   $x->bsub($y);                 # subtraction (subtract $y from $x)
2001   $x->bmul($y);                 # multiplication (multiply $x by $y)
2002   $x->bdiv($y);                 # divide, set $i to quotient
2003                                 # return (quo,rem) or quo if scalar
2004
2005   $x->bmod($y);                 # modulus
2006   $x->bpow($y);                 # power of arguments (a**b)
2007   $x->blsft($y);                # left shift
2008   $x->brsft($y);                # right shift
2009                                 # return (quo,rem) or quo if scalar
2010
2011   $x->blog($base);              # logarithm of $x, base defaults to e
2012                                 # (other bases than e not supported yet)
2013
2014   $x->band($y);                 # bit-wise and
2015   $x->bior($y);                 # bit-wise inclusive or
2016   $x->bxor($y);                 # bit-wise exclusive or
2017   $x->bnot();                   # bit-wise not (two's complement)
2018
2019   $x->bsqrt();                  # calculate square-root
2020   $x->bfac();                   # factorial of $x (1*2*3*4*..$x)
2021
2022   $x->bround($N);               # accuracy: preserver $N digits
2023   $x->bfround($N);              # precision: round to the $Nth digit
2024
2025   # The following do not modify their arguments:
2026   bgcd(@values);                # greatest common divisor
2027   blcm(@values);                # lowest common multiplicator
2028
2029   $x->bstr();                   # return string
2030   $x->bsstr();                  # return string in scientific notation
2031
2032   $x->bfloor();                 # return integer less or equal than $x
2033   $x->bceil();                  # return integer greater or equal than $x
2034
2035   $x->exponent();               # return exponent as BigInt
2036   $x->mantissa();               # return mantissa as BigInt
2037   $x->parts();                  # return (mantissa,exponent) as BigInt
2038
2039   $x->length();                 # number of digits (w/o sign and '.')
2040   ($l,$f) = $x->length();       # number of digits, and length of fraction
2041
2042   $x->precision();              # return P of $x (or global, if P of $x undef)
2043   $x->precision($n);            # set P of $x to $n
2044   $x->accuracy();               # return A of $x (or global, if A of $x undef)
2045   $x->accuracy($n);             # set A $x to $n
2046
2047   Math::BigFloat->precision();  # get/set global P for all BigFloat objects
2048   Math::BigFloat->accuracy();   # get/set global A for all BigFloat objects
2049
2050 =head1 DESCRIPTION
2051
2052 All operators (inlcuding basic math operations) are overloaded if you
2053 declare your big floating point numbers as
2054
2055   $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2056
2057 Operations with overloaded operators preserve the arguments, which is
2058 exactly what you expect.
2059
2060 =head2 Canonical notation
2061
2062 Input to these routines are either BigFloat objects, or strings of the
2063 following four forms:
2064
2065 =over 2
2066
2067 =item *
2068
2069 C</^[+-]\d+$/>
2070
2071 =item *
2072
2073 C</^[+-]\d+\.\d*$/>
2074
2075 =item *
2076
2077 C</^[+-]\d+E[+-]?\d+$/>
2078
2079 =item *
2080
2081 C</^[+-]\d*\.\d+E[+-]?\d+$/>
2082
2083 =back
2084
2085 all with optional leading and trailing zeros and/or spaces. Additonally,
2086 numbers are allowed to have an underscore between any two digits.
2087
2088 Empty strings as well as other illegal numbers results in 'NaN'.
2089
2090 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2091 are always stored in normalized form. On a string, it creates a BigFloat
2092 object.
2093
2094 =head2 Output
2095
2096 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2097
2098 The string output will always have leading and trailing zeros stripped and drop
2099 a plus sign. C<bstr()> will give you always the form with a decimal point,
2100 while C<bsstr()> (for scientific) gives you the scientific notation.
2101
2102         Input                   bstr()          bsstr()
2103         '-0'                    '0'             '0E1'
2104         '  -123 123 123'        '-123123123'    '-123123123E0'
2105         '00.0123'               '0.0123'        '123E-4'
2106         '123.45E-2'             '1.2345'        '12345E-4'
2107         '10E+3'                 '10000'         '1E4'
2108
2109 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2110 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2111 return either undef, <0, 0 or >0 and are suited for sort.
2112
2113 Actual math is done by using BigInts to represent the mantissa and exponent.
2114 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2115 represent the result when input arguments are not numbers, as well as
2116 the result of dividing by zero.
2117
2118 =head2 C<mantissa()>, C<exponent()> and C<parts()>
2119
2120 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2121 as BigInts such that:
2122
2123         $m = $x->mantissa();
2124         $e = $x->exponent();
2125         $y = $m * ( 10 ** $e );
2126         print "ok\n" if $x == $y;
2127
2128 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2129
2130 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2131
2132 Currently the mantissa is reduced as much as possible, favouring higher
2133 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2134 This might change in the future, so do not depend on it.
2135
2136 =head2 Accuracy vs. Precision
2137
2138 See also: L<Rounding|Rounding>.
2139
2140 Math::BigFloat supports both precision and accuracy. For a full documentation,
2141 examples and tips on these topics please see the large section in
2142 L<Math::BigInt>.
2143
2144 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2145 a operation consumes all resources, each operation produces no more than
2146 C<Math::BigFloat::precision()> digits.
2147
2148 In case the result of one operation has more precision than specified,
2149 it is rounded. The rounding mode taken is either the default mode, or the one
2150 supplied to the operation after the I<scale>:
2151
2152         $x = Math::BigFloat->new(2);
2153         Math::BigFloat::precision(5);           # 5 digits max
2154         $y = $x->copy()->bdiv(3);               # will give 0.66666
2155         $y = $x->copy()->bdiv(3,6);             # will give 0.666666
2156         $y = $x->copy()->bdiv(3,6,'odd');       # will give 0.666667
2157         Math::BigFloat::round_mode('zero');
2158         $y = $x->copy()->bdiv(3,6);             # will give 0.666666
2159
2160 =head2 Rounding
2161
2162 =over 2
2163
2164 =item ffround ( +$scale )
2165
2166 Rounds to the $scale'th place left from the '.', counting from the dot.
2167 The first digit is numbered 1.
2168
2169 =item ffround ( -$scale )
2170
2171 Rounds to the $scale'th place right from the '.', counting from the dot.
2172
2173 =item ffround ( 0 )
2174
2175 Rounds to an integer.
2176
2177 =item fround  ( +$scale )
2178
2179 Preserves accuracy to $scale digits from the left (aka significant digits)
2180 and pads the rest with zeros. If the number is between 1 and -1, the
2181 significant digits count from the first non-zero after the '.'
2182
2183 =item fround  ( -$scale ) and fround ( 0 )
2184
2185 These are effetively no-ops.
2186
2187 =back
2188
2189 All rounding functions take as a second parameter a rounding mode from one of
2190 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2191
2192 The default rounding mode is 'even'. By using
2193 C<< Math::BigFloat::round_mode($round_mode); >> you can get and set the default
2194 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2195 no longer supported.
2196 The second parameter to the round functions then overrides the default
2197 temporarily.
2198
2199 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
2200 'trunc' as rounding mode to make it equivalent to:
2201
2202         $x = 2.5;
2203         $y = int($x) + 2;
2204
2205 You can override this by passing the desired rounding mode as parameter to
2206 C<as_number()>:
2207
2208         $x = Math::BigFloat->new(2.5);
2209         $y = $x->as_number('odd');      # $y = 3
2210
2211 =head1 EXAMPLES
2212
2213   # not ready yet
2214
2215 =head1 Autocreating constants
2216
2217 After C<use Math::BigFloat ':constant'> all the floating point constants
2218 in the given scope are converted to C<Math::BigFloat>. This conversion
2219 happens at compile time.
2220
2221 In particular
2222
2223   perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2224
2225 prints the value of C<2E-100>. Note that without conversion of
2226 constants the expression 2E-100 will be calculated as normal floating point
2227 number.
2228
2229 Please note that ':constant' does not affect integer constants, nor binary
2230 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2231 work.
2232
2233 =head2 Math library
2234
2235 Math with the numbers is done (by default) by a module called
2236 Math::BigInt::Calc. This is equivalent to saying:
2237
2238         use Math::BigFloat lib => 'Calc';
2239
2240 You can change this by using:
2241
2242         use Math::BigFloat lib => 'BitVect';
2243
2244 The following would first try to find Math::BigInt::Foo, then
2245 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2246
2247         use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2248
2249 Calc.pm uses as internal format an array of elements of some decimal base
2250 (usually 1e7, but this might be differen for some systems) with the least
2251 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2252 significant bit first. Other modules might use even different means of
2253 representing the numbers. See the respective module documentation for further
2254 details.
2255
2256 Please note that Math::BigFloat does B<not> use the denoted library itself,
2257 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2258 to do:
2259
2260         use Math::BigInt lib => 'GMP';
2261         use Math::BigFloat;
2262
2263 you can roll it all into one line:
2264
2265         use Math::BigFloat lib => 'GMP';
2266
2267 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details.
2268
2269 =head2 Using Math::BigInt::Lite
2270
2271 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2272
2273         # 1
2274         use Math::BigFloat with => 'Math::BigInt::Lite';
2275
2276 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2277 can combine these if you want. For instance, you may want to use
2278 Math::BigInt objects in your main script, too.
2279
2280         # 2
2281         use Math::BigInt;
2282         use Math::BigFloat with => 'Math::BigInt::Lite';
2283
2284 Of course, you can combine this with the C<lib> parameter.
2285
2286         # 3
2287         use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2288
2289 If you want to use Math::BigInt's, too, simple add a Math::BigInt B<before>:
2290
2291         # 4
2292         use Math::BigInt;
2293         use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2294
2295 Notice that the module with the last C<lib> will "win" and thus
2296 it's lib will be used if the lib is available:
2297
2298         # 5
2299         use Math::BigInt lib => 'Bar,Baz';
2300         use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2301
2302 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2303 words, Math::BigFloat will try to retain previously loaded libs when you
2304 don't specify it one.
2305
2306 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2307 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2308 same as trying the latter load alone, except for the fact that Bar or Baz
2309 might be loaded needlessly in an intermidiate step
2310
2311 The old way still works though:
2312
2313         # 6
2314         use Math::BigInt lib => 'Bar,Baz';
2315         use Math::BigFloat;
2316
2317 But B<examples #3 and #4 are recommended> for usage.
2318
2319 =head1 BUGS
2320
2321 =over 2
2322
2323 =item *
2324
2325 The following does not work yet:
2326
2327         $m = $x->mantissa();
2328         $e = $x->exponent();
2329         $y = $m * ( 10 ** $e );
2330         print "ok\n" if $x == $y;
2331
2332 =item *
2333
2334 There is no fmod() function yet.
2335
2336 =back
2337
2338 =head1 CAVEAT
2339
2340 =over 1
2341
2342 =item stringify, bstr()
2343
2344 Both stringify and bstr() now drop the leading '+'. The old code would return
2345 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2346 reasoning and details.
2347
2348 =item bdiv
2349
2350 The following will probably not do what you expect:
2351
2352         print $c->bdiv(123.456),"\n";
2353
2354 It prints both quotient and reminder since print works in list context. Also,
2355 bdiv() will modify $c, so be carefull. You probably want to use
2356
2357         print $c / 123.456,"\n";
2358         print scalar $c->bdiv(123.456),"\n";  # or if you want to modify $c
2359
2360 instead.
2361
2362 =item Modifying and =
2363
2364 Beware of:
2365
2366         $x = Math::BigFloat->new(5);
2367         $y = $x;
2368
2369 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2370 a second reference to the B<same> object and stores it in $y. Thus anything
2371 that modifies $x will modify $y, and vice versa.
2372
2373         $x->bmul(2);
2374         print "$x, $y\n";       # prints '10, 10'
2375
2376 If you want a true copy of $x, use:
2377
2378         $y = $x->copy();
2379
2380 See also the documentation in L<overload> regarding C<=>.
2381
2382 =item bpow
2383
2384 C<bpow()> now modifies the first argument, unlike the old code which left
2385 it alone and only returned the result. This is to be consistent with
2386 C<badd()> etc. The first will modify $x, the second one won't:
2387
2388         print bpow($x,$i),"\n";         # modify $x
2389         print $x->bpow($i),"\n";        # ditto
2390         print $x ** $i,"\n";            # leave $x alone
2391
2392 =back
2393
2394 =head1 LICENSE
2395
2396 This program is free software; you may redistribute it and/or modify it under
2397 the same terms as Perl itself.
2398
2399 =head1 AUTHORS
2400
2401 Mark Biggar, overloaded interface by Ilya Zakharevich.
2402 Completely rewritten by Tels http://bloodgate.com in 2001.
2403
2404 =cut