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