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