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