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