lib/Math/BigFloat.pm

   1 #!/usr/bin/perl -w
   2
   3 # The following hash values are internally used:
   4 #   _e: exponent (BigInt)
   5 #   _m: mantissa (absolute BigInt)
   6 # sign: +,-,"NaN" if not a number
   7 #   _a: accuracy
   8 #   _p: precision
   9 #   _f: flags, used to signal MBI not to touch our private parts
  10 # _cow: Copy-On-Write (NRY)
  11
  12 package Math::BigFloat;
  13
  14 $VERSION = '1.23';
  15 require 5.005;
  16 use Exporter;
  17 use Math::BigInt qw/objectify/;
  18 @ISA =       qw( Exporter Math::BigInt);
  19 # can not export bneg/babs since the are only in MBI
  20 @EXPORT_OK = qw(
  21                 bcmp
  22                 badd bmul bdiv bmod bnorm bsub
  23                 bgcd blcm bround bfround
  24                 bpow bnan bzero bfloor bceil
  25                 bacmp bstr binc bdec binf
  26                 is_odd is_even is_nan is_inf is_positive is_negative
  27                 is_zero is_one sign
  28                );
  29
  30 #@EXPORT = qw( );
  31 use strict;
  32 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode/;
  33 my $class = "Math::BigFloat";
  34
  35 use overload
  36 '<=>'   =>      sub { $_[2] ?
  37                       ref($_[0])->bcmp($_[1],$_[0]) :
  38                       ref($_[0])->bcmp($_[0],$_[1])},
  39 'int'   =>      sub { $_[0]->as_number() },             # 'trunc' to bigint
  40 ;
  41
  42 ##############################################################################
  43 # global constants, flags and accessory
  44
  45 use constant MB_NEVER_ROUND => 0x0001;
  46
  47 # are NaNs ok?
  48 my $NaNOK=1;
  49 # constant for easier life
  50 my $nan = 'NaN';
  51
  52 # class constants, use Class->constant_name() to access
  53 $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
  54 $accuracy   = undef;
  55 $precision  = undef;
  56 $div_scale  = 40;
  57
  58 # in case we call SUPER::->foo() and this wants to call modify()
  59 # sub modify () { 0; }
  60
  61 {
  62   # valid method aliases for AUTOLOAD
  63   my %methods = map { $_ => 1 }
  64    qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
  65         fneg fint facmp fcmp fzero fnan finf finc fdec
  66         fceil ffloor
  67       /;
  68   # valid method's that need to be hand-ed up (for AUTOLOAD)
  69   my %hand_ups = map { $_ => 1 }
  70    qw / is_nan is_inf is_negative is_positive
  71         accuracy precision div_scale round_mode fabs babs
  72       /;
  73
  74   sub method_alias { return exists $methods{$_[0]||''}; }
  75   sub method_hand_up { return exists $hand_ups{$_[0]||''}; }
  76 }
  77
  78 ##############################################################################
  79 # constructors
  80
  81 sub new
  82   {
  83   # create a new BigFloat object from a string or another bigfloat object.
  84   # _e: exponent
  85   # _m: mantissa
  86   # sign  => sign (+/-), or "NaN"
  87
  88   my $class = shift;
  89
  90   my $wanted = shift; # avoid numify call by not using || here
  91   return $class->bzero() if !defined $wanted;      # default to 0
  92   return $wanted->copy() if ref($wanted) eq $class;
  93
  94   my $round = shift; $round = 0 if !defined $round; # no rounding as default
  95   my $self = {}; bless $self, $class;
  96   # shortcut for bigints and its subclasses
  97   if ((ref($wanted)) && (ref($wanted) ne $class))
  98     {
  99     $self->{_m} = $wanted->as_number();         # get us a bigint copy
 100     $self->{_e} = Math::BigInt->new(0);
 101     $self->{_m}->babs();
 102     $self->{sign} = $wanted->sign();
 103     return $self->bnorm();
 104     }
 105   # got string
 106   # handle '+inf', '-inf' first
 107   if ($wanted =~ /^[+-]?inf$/)
 108     {
 109     $self->{_e} = Math::BigInt->new(0);
 110     $self->{_m} = Math::BigInt->new(0);
 111     $self->{sign} = $wanted;
 112     $self->{sign} = '+inf' if $self->{sign} eq 'inf';
 113     return $self->bnorm();
 114     }
 115   #print "new string '$wanted'\n";
 116   my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
 117   if (!ref $mis)
 118     {
 119     die "$wanted is not a number initialized to $class" if !$NaNOK;
 120     $self->{_e} = Math::BigInt->new(0);
 121     $self->{_m} = Math::BigInt->new(0);
 122     $self->{sign} = $nan;
 123     }
 124   else
 125     {
 126     # make integer from mantissa by adjusting exp, then convert to bigint
 127     $self->{_e} = Math::BigInt->new("$$es$$ev");        # exponent
 128     $self->{_m} = Math::BigInt->new("$$mis$$miv$$mfv"); # create mantissa
 129     # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
 130     $self->{_e} -= CORE::length($$mfv);
 131     $self->{sign} = $self->{_m}->sign(); $self->{_m}->babs();
 132     }
 133   #print "$wanted => $self->{sign} $self->{value}\n";
 134   $self->bnorm();       # first normalize
 135   # if any of the globals is set, round to them and thus store them insid $self
 136   $self->round($accuracy,$precision,$class->round_mode)
 137    if defined $accuracy || defined $precision;
 138   return $self;
 139   }
 140
 141 sub bnan
 142   {
 143   # create a bigfloat 'NaN', if given a BigFloat, set it to 'NaN'
 144   my $self = shift;
 145   $self = $class if !defined $self;
 146   if (!ref($self))
 147     {
 148     my $c = $self; $self = {}; bless $self, $c;
 149     }
 150   $self->{_m} = Math::BigInt->bzero();
 151   $self->{_e} = Math::BigInt->bzero();
 152   $self->{sign} = $nan;
 153   return $self;
 154   }
 155
 156 sub binf
 157   {
 158   # create a bigfloat '+-inf', if given a BigFloat, set it to '+-inf'
 159   my $self = shift;
 160   my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
 161
 162   $self = $class if !defined $self;
 163   if (!ref($self))
 164     {
 165     my $c = $self; $self = {}; bless $self, $c;
 166     }
 167   $self->{_m} = Math::BigInt->bzero();
 168   $self->{_e} = Math::BigInt->bzero();
 169   $self->{sign} = $sign.'inf';
 170   return $self;
 171   }
 172
 173 sub bone
 174   {
 175   # create a bigfloat '+-1', if given a BigFloat, set it to '+-1'
 176   my $self = shift;
 177   my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
 178
 179   $self = $class if !defined $self;
 180   if (!ref($self))
 181     {
 182     my $c = $self; $self = {}; bless $self, $c;
 183     }
 184   $self->{_m} = Math::BigInt->bone();
 185   $self->{_e} = Math::BigInt->bzero();
 186   $self->{sign} = $sign;
 187   return $self;
 188   }
 189
 190 sub bzero
 191   {
 192   # create a bigfloat '+0', if given a BigFloat, set it to 0
 193   my $self = shift;
 194   $self = $class if !defined $self;
 195   if (!ref($self))
 196     {
 197     my $c = $self; $self = {}; bless $self, $c;
 198     }
 199   $self->{_m} = Math::BigInt->bzero();
 200   $self->{_e} = Math::BigInt->bone();
 201   $self->{sign} = '+';
 202   return $self;
 203   }
 204
 205 ##############################################################################
 206 # string conversation
 207
 208 sub bstr
 209   {
 210   # (ref to BFLOAT or num_str ) return num_str
 211   # Convert number from internal format to (non-scientific) string format.
 212   # internal format is always normalized (no leading zeros, "-0" => "+0")
 213   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 214   #my $x = shift; my $class = ref($x) || $x;
 215   #$x = $class->new(shift) unless ref($x);
 216
 217   #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
 218   #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
 219   if ($x->{sign} !~ /^[+-]$/)
 220     {
 221     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
 222     return 'inf';                                       # +inf
 223     }
 224
 225   my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
 226
 227   my $not_zero = !$x->is_zero();
 228   if ($not_zero)
 229     {
 230     $es = $x->{_m}->bstr();
 231     $len = CORE::length($es);
 232     if (!$x->{_e}->is_zero())
 233 #      {
 234 #      $es = $x->{sign}.$es if $x->{sign} eq '-';
 235 #      }
 236 #    else
 237       {
 238       if ($x->{_e}->sign() eq '-')
 239         {
 240         $dot = '';
 241         if ($x->{_e} <= -$len)
 242           {
 243           # print "style: 0.xxxx\n";
 244           my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
 245           $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
 246           }
 247         else
 248           {
 249           # print "insert '.' at $x->{_e} in '$es'\n";
 250           substr($es,$x->{_e},0) = '.'; $cad = $x->{_e};
 251           }
 252         }
 253       else
 254         {
 255         # expand with zeros
 256         $es .= '0' x $x->{_e}; $len += $x->{_e}; $cad = 0;
 257         }
 258       }
 259     } # if not zero
 260   $es = $x->{sign}.$es if $x->{sign} eq '-';
 261   # if set accuracy or precision, pad with zeros
 262   if ((defined $x->{_a}) && ($not_zero))
 263     {
 264     # 123400 => 6, 0.1234 => 4, 0.001234 => 4
 265     my $zeros = $x->{_a} - $cad;                # cad == 0 => 12340
 266     $zeros = $x->{_a} - $len if $cad != $len;
 267     #print "acc padd $x->{_a} $zeros (len $len cad $cad)\n";
 268     $es .= $dot.'0' x $zeros if $zeros > 0;
 269     }
 270   elsif ($x->{_p} || 0 < 0)
 271     {
 272     # 123400 => 6, 0.1234 => 4, 0.001234 => 6
 273     my $zeros = -$x->{_p} + $cad;
 274     #print "pre padd $x->{_p} $zeros (len $len cad $cad)\n";
 275     $es .= $dot.'0' x $zeros if $zeros > 0;
 276     }
 277   return $es;
 278   }
 279
 280 sub bsstr
 281   {
 282   # (ref to BFLOAT or num_str ) return num_str
 283   # Convert number from internal format to scientific string format.
 284   # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
 285   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 286   #my $x = shift; my $class = ref($x) || $x;
 287   #$x = $class->new(shift) unless ref($x);
 288
 289   #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
 290   #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
 291   if ($x->{sign} !~ /^[+-]$/)
 292     {
 293     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
 294     return 'inf';                                       # +inf
 295     }
 296   my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
 297   my $sep = 'e'.$sign;
 298   return $x->{_m}->bstr().$sep.$x->{_e}->bstr();
 299   }
 300
 301 sub numify
 302   {
 303   # Make a number from a BigFloat object
 304   # simple return string and let Perl's atoi()/atof() handle the rest
 305   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 306   return $x->bsstr();
 307   }
 308
 309 ##############################################################################
 310 # public stuff (usually prefixed with "b")
 311
 312 # really? Just for exporting them is not what I had in mind
 313 #sub babs
 314 #  {
 315 #  $class->SUPER::babs($class,@_);
 316 #  }
 317 #sub bneg
 318 #  {
 319 #  $class->SUPER::bneg($class,@_);
 320 #  }
 321
 322 # tels 2001-08-04
 323 # todo: this must be overwritten and return NaN for non-integer values
 324 # band(), bior(), bxor(), too
 325 #sub bnot
 326 #  {
 327 #  $class->SUPER::bnot($class,@_);
 328 #  }
 329
 330 sub bcmp
 331   {
 332   # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
 333   # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
 334   my ($self,$x,$y) = objectify(2,@_);
 335
 336   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
 337     {
 338     # handle +-inf and NaN
 339     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 340     return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
 341     return +1 if $x->{sign} eq '+inf';
 342     return -1 if $x->{sign} eq '-inf';
 343     return -1 if $y->{sign} eq '+inf';
 344     return +1 if $y->{sign} eq '-inf';
 345     }
 346
 347   # check sign for speed first
 348   return 1 if $x->{sign} eq '+' && $y->{sign} eq '-';   # does also 0 <=> -y
 349   return -1 if $x->{sign} eq '-' && $y->{sign} eq '+';  # does also -x <=> 0
 350
 351   # shortcut
 352   my $xz = $x->is_zero();
 353   my $yz = $y->is_zero();
 354   return 0 if $xz && $yz;                               # 0 <=> 0
 355   return -1 if $xz && $y->{sign} eq '+';                # 0 <=> +y
 356   return 1 if $yz && $x->{sign} eq '+';                 # +x <=> 0
 357
 358   # adjust so that exponents are equal
 359   my $lxm = $x->{_m}->length();
 360   my $lym = $y->{_m}->length();
 361   my $lx = $lxm + $x->{_e};
 362   my $ly = $lym + $y->{_e};
 363   # print "x $x y $y lx $lx ly $ly\n";
 364   my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
 365   # print "$l $x->{sign}\n";
 366   return $l <=> 0 if $l != 0;
 367
 368   # lengths (corrected by exponent) are equal
 369   # so make mantissa euqal length by padding with zero (shift left)
 370   my $diff = $lxm - $lym;
 371   my $xm = $x->{_m};            # not yet copy it
 372   my $ym = $y->{_m};
 373   if ($diff > 0)
 374     {
 375     $ym = $y->{_m}->copy()->blsft($diff,10);
 376     }
 377   elsif ($diff < 0)
 378     {
 379     $xm = $x->{_m}->copy()->blsft(-$diff,10);
 380     }
 381   my $rc = $xm->bcmp($ym);
 382   $rc = -$rc if $x->{sign} eq '-';              # -124 < -123
 383   return $rc <=> 0;
 384   }
 385
 386 sub bacmp
 387   {
 388   # Compares 2 values, ignoring their signs.
 389   # Returns one of undef, <0, =0, >0. (suitable for sort)
 390   # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
 391   my ($self,$x,$y) = objectify(2,@_);
 392
 393   # handle +-inf and NaN's
 394   if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]/)
 395     {
 396     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 397     return 0 if ($x->is_inf() && $y->is_inf());
 398     return 1 if ($x->is_inf() && !$y->is_inf());
 399     return -1 if (!$x->is_inf() && $y->is_inf());
 400     }
 401
 402   # shortcut
 403   my $xz = $x->is_zero();
 404   my $yz = $y->is_zero();
 405   return 0 if $xz && $yz;                               # 0 <=> 0
 406   return -1 if $xz && !$yz;                             # 0 <=> +y
 407   return 1 if $yz && !$xz;                              # +x <=> 0
 408
 409   # adjust so that exponents are equal
 410   my $lxm = $x->{_m}->length();
 411   my $lym = $y->{_m}->length();
 412   my $lx = $lxm + $x->{_e};
 413   my $ly = $lym + $y->{_e};
 414   # print "x $x y $y lx $lx ly $ly\n";
 415   my $l = $lx - $ly; # $l = -$l if $x->{sign} eq '-';
 416   # print "$l $x->{sign}\n";
 417   return $l <=> 0 if $l != 0;
 418
 419   # lengths (corrected by exponent) are equal
 420   # so make mantissa euqal length by padding with zero (shift left)
 421   my $diff = $lxm - $lym;
 422   my $xm = $x->{_m};            # not yet copy it
 423   my $ym = $y->{_m};
 424   if ($diff > 0)
 425     {
 426     $ym = $y->{_m}->copy()->blsft($diff,10);
 427     }
 428   elsif ($diff < 0)
 429     {
 430     $xm = $x->{_m}->copy()->blsft(-$diff,10);
 431     }
 432   my $rc = $xm->bcmp($ym);
 433   # $rc = -$rc if $x->{sign} eq '-';            # -124 < -123
 434   return $rc <=> 0;
 435
 436 #  # signs are ignored, so check length
 437 #  # length(x) is length(m)+e aka length of non-fraction part
 438 #  # the longer one is bigger
 439 #  my $l = $x->length() - $y->length();
 440 #  #print "$l\n";
 441 #  return $l if $l != 0;
 442 #  #print "equal lengths\n";
 443 #
 444 #  # if both are equal long, make full compare
 445 #  # first compare only the mantissa
 446 #  # if mantissa are equal, compare fractions
 447 #
 448 #  return $x->{_m} <=> $y->{_m} || $x->{_e} <=> $y->{_e};
 449   }
 450
 451 sub badd
 452   {
 453   # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
 454   # return result as BFLOAT
 455   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 456
 457   # inf and NaN handling
 458   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
 459     {
 460     # NaN first
 461     return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 462     # inf handline
 463     if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
 464       {
 465       # + and + => +, - and - => -, + and - => 0, - and + => 0
 466       return $x->bzero() if $x->{sign} ne $y->{sign};
 467       return $x;
 468       }
 469     # +-inf + something => +inf
 470     # something +-inf => +-inf
 471     $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
 472     return $x;
 473     }
 474
 475   # speed: no add for 0+y or x+0
 476   return $x if $y->is_zero();                           # x+0
 477   if ($x->is_zero())                                    # 0+y
 478     {
 479     # make copy, clobbering up x (modify in place!)
 480     $x->{_e} = $y->{_e}->copy();
 481     $x->{_m} = $y->{_m}->copy();
 482     $x->{sign} = $y->{sign} || $nan;
 483     return $x->round($a,$p,$r,$y);
 484     }
 485
 486   # take lower of the two e's and adapt m1 to it to match m2
 487   my $e = $y->{_e}; $e = Math::BigInt::bzero() if !defined $e;  # if no BFLOAT
 488   $e = $e - $x->{_e};
 489   my $add = $y->{_m}->copy();
 490   if ($e < 0)
 491     {
 492     # print "e < 0\n";
 493     #print "\$x->{_m}: $x->{_m} ";
 494     #print "\$x->{_e}: $x->{_e}\n";
 495     my $e1 = $e->copy()->babs();
 496     $x->{_m} *= (10 ** $e1);
 497     $x->{_e} += $e;                     # need the sign of e
 498     #$x->{_m} += $y->{_m};
 499     #print "\$x->{_m}: $x->{_m} ";
 500     #print "\$x->{_e}: $x->{_e}\n";
 501     }
 502   elsif ($e > 0)
 503     {
 504     # print "e > 0\n";
 505     #print "\$x->{_m}: $x->{_m} \$y->{_m}: $y->{_m} \$e: $e ",ref($e),"\n";
 506     $add *= (10 ** $e);
 507     #$x->{_m} += $y->{_m} * (10 ** $e);
 508     #print "\$x->{_m}: $x->{_m}\n";
 509     }
 510   # else: both e are same, so leave them
 511   #print "badd $x->{sign}$x->{_m} +  $y->{sign}$add\n";
 512   # fiddle with signs
 513   $x->{_m}->{sign} = $x->{sign};
 514   $add->{sign} = $y->{sign};
 515   # finally do add/sub
 516   $x->{_m} += $add;
 517   # re-adjust signs
 518   $x->{sign} = $x->{_m}->{sign};
 519   $x->{_m}->{sign} = '+';
 520   #$x->bnorm();                         # delete trailing zeros
 521   return $x->round($a,$p,$r,$y);
 522   }
 523
 524 sub bsub
 525   {
 526   # (BigFloat or num_str, BigFloat or num_str) return BigFloat
 527   # subtract second arg from first, modify first
 528   my ($self,$x,$y) = objectify(2,@_);
 529
 530   $x->badd($y->bneg()); # badd does not leave internal zeros
 531   $y->bneg();           # refix y, assumes no one reads $y in between
 532   return $x;            # badd() already normalized and rounded
 533   }
 534
 535 sub binc
 536   {
 537   # increment arg by one
 538   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 539   $x->badd($self->bone())->round($a,$p,$r);
 540   }
 541
 542 sub bdec
 543   {
 544   # decrement arg by one
 545   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 546   $x->badd($self->bone('-'))->round($a,$p,$r);
 547   }
 548
 549 sub blcm
 550   {
 551   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
 552   # does not modify arguments, but returns new object
 553   # Lowest Common Multiplicator
 554
 555   my ($self,@arg) = objectify(0,@_);
 556   my $x = $self->new(shift @arg);
 557   while (@arg) { $x = _lcm($x,shift @arg); }
 558   $x;
 559   }
 560
 561 sub bgcd
 562   {
 563   # (BFLOAT or num_str, BFLOAT or num_str) return BINT
 564   # does not modify arguments, but returns new object
 565   # GCD -- Euclids algorithm Knuth Vol 2 pg 296
 566
 567   my ($self,@arg) = objectify(0,@_);
 568   my $x = $self->new(shift @arg);
 569   while (@arg) { $x = _gcd($x,shift @arg); }
 570   $x;
 571   }
 572
 573 sub is_zero
 574   {
 575   # return true if arg (BFLOAT or num_str) is zero (array '+', '0')
 576   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 577
 578   return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
 579   return 0;
 580   }
 581
 582 sub is_one
 583   {
 584   # return true if arg (BFLOAT or num_str) is +1 (array '+', '1')
 585   # or -1 if signis given
 586   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 587
 588   my $sign = shift || ''; $sign = '+' if $sign ne '-';
 589   return 1
 590    if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one());
 591   return 0;
 592   }
 593
 594 sub is_odd
 595   {
 596   # return true if arg (BFLOAT or num_str) is odd or false if even
 597   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 598
 599   return 0 if $x->{sign} !~ /^[+-]$/;                   # NaN & +-inf aren't
 600   return 1 if ($x->{_e}->is_zero() && $x->{_m}->is_odd());
 601   return 0;
 602   }
 603
 604 sub is_even
 605   {
 606   # return true if arg (BINT or num_str) is even or false if odd
 607   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
 608
 609   return 0 if $x->{sign} !~ /^[+-]$/;                   # NaN & +-inf aren't
 610   return 1 if $x->{_m}->is_zero();                      # 0e1 is even
 611   return 1 if ($x->{_e}->is_zero() && $x->{_m}->is_even()); # 123.45 is never
 612   return 0;
 613   }
 614
 615 sub bmul
 616   {
 617   # multiply two numbers -- stolen from Knuth Vol 2 pg 233
 618   # (BINT or num_str, BINT or num_str) return BINT
 619   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 620
 621   # print "mbf bmul $x->{_m}e$x->{_e} $y->{_m}e$y->{_e}\n";
 622   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 623
 624   # handle result = 0
 625   return $x->bzero() if $x->is_zero() || $y->is_zero();
 626   # inf handling
 627   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
 628     {
 629     # result will always be +-inf:
 630     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
 631     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
 632     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
 633     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
 634     return $x->binf('-');
 635     }
 636
 637   # aEb * cEd = (a*c)E(b+d)
 638   $x->{_m} = $x->{_m} * $y->{_m};
 639   #print "m: $x->{_m}\n";
 640   $x->{_e} = $x->{_e} + $y->{_e};
 641   #print "e: $x->{_m}\n";
 642   # adjust sign:
 643   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
 644   #print "s: $x->{sign}\n";
 645   $x->bnorm();
 646   return $x->round($a,$p,$r,$y);
 647   }
 648
 649 sub bdiv
 650   {
 651   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
 652   # (BFLOAT,BFLOAT) (quo,rem) or BINT (only rem)
 653   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 654
 655
 656   # x / +-inf => 0, reminder x
 657   return wantarray ? ($x->bzero(),$x->copy()) : $x->bzero()
 658    if $y->{sign} =~ /^[+-]inf$/;
 659
 660   # NaN if x == NaN or y == NaN or x==y==0
 661   return wantarray ? ($x->bnan(),bnan()) : $x->bnan()
 662    if (($x->is_nan() || $y->is_nan()) ||
 663       ($x->is_zero() && $y->is_zero()));
 664
 665   # 5 / 0 => +inf, -6 / 0 => -inf
 666   return wantarray
 667    ? ($x->binf($x->{sign}),$self->bnan()) : $x->binf($x->{sign})
 668    if ($x->{sign} =~ /^[+-]$/ && $y->is_zero());
 669
 670   # promote BigInts and it's subclasses (except when already a BigFloat)
 671   $y = $self->new($y) unless $y->isa('Math::BigFloat');
 672
 673   # old, broken way
 674   # $y = $class->new($y) if ref($y) ne $self;           # promote bigints
 675
 676   # print "mbf bdiv $x ",ref($x)," ",$y," ",ref($y),"\n";
 677   # we need to limit the accuracy to protect against overflow
 678
 679   my $fallback = 0;
 680   my $scale = 0;
 681 #  print "s=$scale a=",$a||'undef'," p=",$p||'undef'," r=",$r||'undef',"\n";
 682   my @params = $x->_find_round_parameters($a,$p,$r,$y);
 683
 684   # no rounding at all, so must use fallback
 685   if (scalar @params == 1)
 686     {
 687     # simulate old behaviour
 688     $scale = $self->div_scale()+1;      # at least one more for proper round
 689     $params[1] = $self->div_scale();    # and round to it as accuracy
 690     $params[3] = $r;                    # round mode by caller or undef
 691     $fallback = 1;                      # to clear a/p afterwards
 692     }
 693   else
 694     {
 695     # the 4 below is empirical, and there might be cases where it is not
 696     # enough...
 697     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
 698     }
 699  # print "s=$scale a=",$params[1]||'undef'," p=",$params[2]||'undef'," f=$fallback\n";
 700   my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
 701   $scale = $lx if $lx > $scale;
 702   $scale = $ly if $ly > $scale;
 703 #  print "scale $scale $lx $ly\n";
 704   my $diff = $ly - $lx;
 705   $scale += $diff if $diff > 0;         # if lx << ly, but not if ly << lx!
 706
 707   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
 708
 709   $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
 710
 711   # check for / +-1 ( +/- 1E0)
 712   if ($y->is_one())
 713     {
 714     return wantarray ? ($x,$self->bzero()) : $x;
 715     }
 716
 717   # calculate the result to $scale digits and then round it
 718   # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
 719   #$scale = 82;
 720   #print "self: $self x: $x ref(x) ", ref($x)," m: $x->{_m}\n";
 721   $x->{_m}->blsft($scale,10);
 722   #print "m: $x->{_m} $y->{_m}\n";
 723   $x->{_m}->bdiv( $y->{_m} );   # a/c
 724   #print "m: $x->{_m}\n";
 725   #print "e: $x->{_e} $y->{_e} ",$scale,"\n";
 726   $x->{_e}->bsub($y->{_e});     # b-d
 727   #print "e: $x->{_e}\n";
 728   $x->{_e}->bsub($scale);       # correct for 10**scale
 729   #print "after div: m: $x->{_m} e: $x->{_e}\n";
 730   $x->bnorm();                  # remove trailing 0's
 731   #print "after norm: m: $x->{_m} e: $x->{_e}\n";
 732
 733   # shortcut to not run trough _find_round_parameters again
 734   if (defined $params[1])
 735     {
 736     $x->bround($params[1],undef,$params[3]);    # then round accordingly
 737     }
 738   else
 739     {
 740     $x->bfround($params[2],$params[3]);         # then round accordingly
 741     }
 742   if ($fallback)
 743     {
 744     # clear a/p after round, since user did not request it
 745     $x->{_a} = undef; $x->{_p} = undef;
 746     }
 747
 748   if (wantarray)
 749     {
 750     my $rem = $x->copy();
 751     $rem->bmod($y,$params[1],$params[2],$params[3]);
 752     if ($fallback)
 753       {
 754       # clear a/p after round, since user did not request it
 755       $rem->{_a} = undef; $rem->{_p} = undef;
 756       }
 757     return ($x,$rem);
 758     }
 759   return $x;
 760   }
 761
 762 sub bmod
 763   {
 764   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
 765   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 766
 767   return $x->bnan() if ($x->{sign} eq $nan || $y->is_nan() || $y->is_zero());
 768   return $x->bzero() if $y->is_one();
 769
 770   # XXX tels: not done yet
 771   return $x->round($a,$p,$r,$y);
 772   }
 773
 774 sub bsqrt
 775   {
 776   # calculate square root; this should probably
 777   # use a different test to see whether the accuracy we want is...
 778   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
 779
 780   return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
 781   return $x if $x->{sign} eq '+inf';                              # +inf
 782   return $x if $x->is_zero() || $x == 1;
 783
 784   # we need to limit the accuracy to protect against overflow (ignore $p)
 785   my ($scale) = $x->_scale_a($self->accuracy(),$self->round_mode,$a,$r);
 786   my $fallback = 0;
 787   if (!defined $scale)
 788     {
 789     # simulate old behaviour
 790     $scale = $self->div_scale()+1;      # one more for proper riund
 791     $a = $self->div_scale();            # and round to it
 792     $fallback = 1;                      # to clear a/p afterwards
 793     }
 794   my $lx = $x->{_m}->length();
 795   $scale = $lx if $scale < $lx;
 796   my $e = Math::BigFloat->new("1E-$scale");     # make test variable
 797   return $x->bnan() if $e->sign() eq 'NaN';
 798
 799   # start with some reasonable guess
 800   #$x *= 10 ** ($len - $org->{_e}); $x /= 2;    # !?!?
 801   $lx = $lx+$x->{_e};
 802   $lx = 1 if $lx < 1;
 803   my $gs = Math::BigFloat->new('1'. ('0' x $lx));
 804
 805 #   print "first guess: $gs (x $x) scale $scale\n";
 806
 807   my $diff = $e;
 808   my $y = $x->copy();
 809   my $two = Math::BigFloat->new(2);
 810   # promote BigInts and it's subclasses (except when already a BigFloat)
 811   $y = $self->new($y) unless $y->isa('Math::BigFloat');
 812   # old, broken way
 813   # $x = Math::BigFloat->new($x) if ref($x) ne $class;  # promote BigInts
 814   my $rem;
 815   # $scale = 2;
 816   while ($diff >= $e)
 817     {
 818     return $x->bnan() if $gs->is_zero();
 819     $rem = $y->copy(); $rem->bdiv($gs,$scale);
 820     #print "y $y gs $gs ($gs->{_a}) rem (y/gs)\n $rem\n";
 821     $x = ($rem + $gs);
 822     #print "x $x rem $rem gs $gs gsa: $gs->{_a}\n";
 823     $x->bdiv($two,$scale);
 824     #print "x $x (/2)\n";
 825     $diff = $x->copy()->bsub($gs)->babs();
 826     $gs = $x->copy();
 827     }
 828 #  print "before $x $x->{_a} ",$a||'a undef'," ",$p||'p undef',"\n";
 829   $x->round($a,$p,$r);
 830 #  print "after $x $x->{_a} ",$a||'a undef'," ",$p||'p undef',"\n";
 831   if ($fallback)
 832     {
 833     # clear a/p after round, since user did not request it
 834     $x->{_a} = undef; $x->{_p} = undef;
 835     }
 836   $x;
 837   }
 838
 839 sub bpow
 840   {
 841   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
 842   # compute power of two numbers, second arg is used as integer
 843   # modifies first argument
 844
 845   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 846
 847   return $x if $x->{sign} =~ /^[+-]inf$/;
 848   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
 849   return $x->bone() if $y->is_zero();
 850   return $x         if $x->is_one() || $y->is_one();
 851   my $y1 = $y->as_number();             # make bigint (trunc)
 852   if ($x == -1)
 853     {
 854     # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
 855     return $y1->is_odd() ? $x : $x->babs(1);
 856     }
 857   return $x if $x->is_zero() && $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
 858   # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
 859   return $x->binf() if $x->is_zero() && $y->{sign} eq '-';
 860
 861   # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
 862   $y1->babs();
 863   $x->{_m}->bpow($y1);
 864   $x->{_e}->bmul($y1);
 865   $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
 866   $x->bnorm();
 867   if ($y->{sign} eq '-')
 868     {
 869     # modify $x in place!
 870     my $z = $x->copy(); $x->bzero()->binc();
 871     return $x->bdiv($z,$a,$p,$r);       # round in one go (might ignore y's A!)
 872     }
 873   return $x->round($a,$p,$r,$y);
 874   }
 875
 876 ###############################################################################
 877 # rounding functions
 878
 879 sub bfround
 880   {
 881   # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
 882   # $n == 0 means round to integer
 883   # expects and returns normalized numbers!
 884   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
 885
 886   return $x if $x->modify('bfround');
 887
 888   my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
 889   return $x if !defined $scale;                 # no-op
 890
 891   # never round a 0, +-inf, NaN
 892   return $x if $x->{sign} !~ /^[+-]$/ || $x->is_zero();
 893   # print "MBF bfround $x to scale $scale mode $mode\n";
 894
 895   # don't round if x already has lower precision
 896   return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
 897
 898   $x->{_p} = $scale;                    # remember round in any case
 899   $x->{_a} = undef;                     # and clear A
 900   if ($scale < 0)
 901     {
 902     # print "bfround scale $scale e $x->{_e}\n";
 903     # round right from the '.'
 904     return $x if $x->{_e} >= 0;                 # nothing to round
 905     $scale = -$scale;                           # positive for simplicity
 906     my $len = $x->{_m}->length();               # length of mantissa
 907     my $dad = -$x->{_e};                        # digits after dot
 908     my $zad = 0;                                # zeros after dot
 909     $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
 910     #print "scale $scale dad $dad zad $zad len $len\n";
 911
 912     # number  bsstr   len zad dad
 913     # 0.123   123e-3    3   0 3
 914     # 0.0123  123e-4    3   1 4
 915     # 0.001   1e-3      1   2 3
 916     # 1.23    123e-2    3   0 2
 917     # 1.2345  12345e-4  5   0 4
 918
 919     # do not round after/right of the $dad
 920     return $x if $scale > $dad;                 # 0.123, scale >= 3 => exit
 921
 922     # round to zero if rounding inside the $zad, but not for last zero like:
 923     # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
 924     return $x->bzero() if $scale < $zad;
 925     if ($scale == $zad)                 # for 0.006, scale -3 and trunc
 926       {
 927       $scale = -$len-1;
 928       }
 929     else
 930       {
 931       # adjust round-point to be inside mantissa
 932       if ($zad != 0)
 933         {
 934         $scale = $scale-$zad;
 935         }
 936       else
 937         {
 938         my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;    # digits before dot
 939         $scale = $dbd+$scale;
 940         }
 941       }
 942     # print "round to $x->{_m} to $scale\n";
 943     }
 944   else
 945     {
 946     # 123 => 100 means length(123) = 3 - $scale (2) => 1
 947
 948     # calculate digits before dot
 949     my $dbt = $x->{_m}->length(); $dbt += $x->{_e} if $x->{_e}->sign() eq '-';
 950     # if not enough digits before dot, round to zero
 951     return $x->bzero() if ($scale > $dbt) && ($dbt < 0);
 952     # scale always >= 0 here
 953     if ($dbt == 0)
 954       {
 955       # 0.49->bfround(1): scale == 1, dbt == 0: => 0.0
 956       # 0.51->bfround(0): scale == 0, dbt == 0: => 1.0
 957       # 0.5->bfround(0):  scale == 0, dbt == 0: => 0
 958       # 0.05->bfround(0): scale == 0, dbt == 0: => 0
 959       # print "$scale $dbt $x->{_m}\n";
 960       $scale = -$x->{_m}->length();
 961       }
 962     elsif ($dbt > 0)
 963       {
 964       # correct by subtracting scale
 965       $scale = $dbt - $scale;
 966       }
 967     else
 968       {
 969       $scale = $x->{_m}->length() - $scale;
 970       }
 971     }
 972   # print "using $scale for $x->{_m} with '$mode'\n";
 973   # pass sign to bround for rounding modes '+inf' and '-inf'
 974   $x->{_m}->{sign} = $x->{sign};
 975   $x->{_m}->bround($scale,$mode);
 976   $x->{_m}->{sign} = '+';               # fix sign back
 977   $x->bnorm();
 978   }
 979
 980 sub bround
 981   {
 982   # accuracy: preserve $N digits, and overwrite the rest with 0's
 983   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
 984
 985   die ('bround() needs positive accuracy') if ($_[0] || 0) < 0;
 986
 987   my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
 988   return $x if !defined $scale;                         # no-op
 989
 990   return $x if $x->modify('bround');
 991
 992   # scale is now either $x->{_a}, $accuracy, or the user parameter
 993   # test whether $x already has lower accuracy, do nothing in this case
 994   # but do round if the accuracy is the same, since a math operation might
 995   # want to round a number with A=5 to 5 digits afterwards again
 996   return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
 997
 998   # print "bround $scale $mode\n";
 999   # 0 => return all digits, scale < 0 makes no sense
1000   return $x if ($scale <= 0);
1001   # never round a 0, +-inf, NaN
1002   return $x if $x->{sign} !~ /^[+-]$/ || $x->is_zero();
1003
1004   # if $e longer than $m, we have 0.0000xxxyyy style number, and must
1005   # subtract the delta from scale, to simulate keeping the zeros
1006   # -5 +5 => 1; -10 +5 => -4
1007   my $delta = $x->{_e} + $x->{_m}->length() + 1;
1008
1009   # if we should keep more digits than the mantissa has, do nothing
1010   return $x if $x->{_m}->length() <= $scale;
1011
1012   # pass sign to bround for '+inf' and '-inf' rounding modes
1013   $x->{_m}->{sign} = $x->{sign};
1014   $x->{_m}->bround($scale,$mode);       # round mantissa
1015   $x->{_m}->{sign} = '+';               # fix sign back
1016   $x->{_a} = $scale;                    # remember rounding
1017   $x->{_p} = undef;                     # and clear P
1018   $x->bnorm();                          # del trailing zeros gen. by bround()
1019   }
1020
1021 sub bfloor
1022   {
1023   # return integer less or equal then $x
1024   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1025
1026   return $x if $x->modify('bfloor');
1027
1028   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1029
1030   # if $x has digits after dot
1031   if ($x->{_e}->{sign} eq '-')
1032     {
1033     $x->{_m}->brsft(-$x->{_e},10);
1034     $x->{_e}->bzero();
1035     $x-- if $x->{sign} eq '-';
1036     }
1037   return $x->round($a,$p,$r);
1038   }
1039
1040 sub bceil
1041   {
1042   # return integer greater or equal then $x
1043   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1044
1045   return $x if $x->modify('bceil');
1046   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1047
1048   # if $x has digits after dot
1049   if ($x->{_e}->{sign} eq '-')
1050     {
1051     $x->{_m}->brsft(-$x->{_e},10);
1052     $x->{_e}->bzero();
1053     $x++ if $x->{sign} eq '+';
1054     }
1055   return $x->round($a,$p,$r);
1056   }
1057
1058 ###############################################################################
1059
1060 sub DESTROY
1061   {
1062   # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
1063   }
1064
1065 sub AUTOLOAD
1066   {
1067   # make fxxx and bxxx work
1068   # my $self = $_[0];
1069   my $name = $AUTOLOAD;
1070
1071   $name =~ s/.*:://;    # split package
1072   #print "$name\n";
1073   no strict 'refs';
1074   if (!method_alias($name))
1075     {
1076     if (!defined $name)
1077       {
1078       # delayed load of Carp and avoid recursion
1079       require Carp;
1080       Carp::croak ("Can't call a method without name");
1081       }
1082     # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1083     if (!method_hand_up($name))
1084       {
1085       # delayed load of Carp and avoid recursion
1086       require Carp;
1087       Carp::croak ("Can't call $class\-\>$name, not a valid method");
1088       }
1089     # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1090     $name =~ s/^f/b/;
1091     return &{'Math::BigInt'."::$name"}(@_);
1092     }
1093   my $bname = $name; $bname =~ s/^f/b/;
1094   *{$class."\:\:$name"} = \&$bname;
1095   &$bname;      # uses @_
1096   }
1097
1098 sub exponent
1099   {
1100   # return a copy of the exponent
1101   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1102
1103   if ($x->{sign} !~ /^[+-]$/)
1104     {
1105     my $s = $x->{sign}; $s =~ s/^[+-]//;
1106     return $self->new($s);                      # -inf, +inf => +inf
1107     }
1108   return $x->{_e}->copy();
1109   }
1110
1111 sub mantissa
1112   {
1113   # return a copy of the mantissa
1114   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1115
1116   if ($x->{sign} !~ /^[+-]$/)
1117     {
1118     my $s = $x->{sign}; $s =~ s/^[+]//;
1119     return $self->new($s);                      # -inf, +inf => +inf
1120     }
1121   my $m = $x->{_m}->copy();             # faster than going via bstr()
1122   $m->bneg() if $x->{sign} eq '-';
1123
1124   return $m;
1125   }
1126
1127 sub parts
1128   {
1129   # return a copy of both the exponent and the mantissa
1130   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1131
1132   if ($x->{sign} !~ /^[+-]$/)
1133     {
1134     my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
1135     return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
1136     }
1137   my $m = $x->{_m}->copy();     # faster than going via bstr()
1138   $m->bneg() if $x->{sign} eq '-';
1139   return ($m,$x->{_e}->copy());
1140   }
1141
1142 ##############################################################################
1143 # private stuff (internal use only)
1144
1145 sub import
1146   {
1147   my $self = shift;
1148   for ( my $i = 0; $i < @_ ; $i++ )
1149     {
1150     if ( $_[$i] eq ':constant' )
1151       {
1152       # this rest causes overlord er load to step in
1153       # print "overload @_\n";
1154       overload::constant float => sub { $self->new(shift); };
1155       splice @_, $i, 1; last;
1156       }
1157     }
1158   # any non :constant stuff is handled by our parent, Exporter
1159   # even if @_ is empty, to give it a chance
1160   #$self->SUPER::import(@_);            # does not work (would call MBI)
1161   $self->export_to_level(1,$self,@_);   # need this instead
1162   }
1163
1164 sub bnorm
1165   {
1166   # adjust m and e so that m is smallest possible
1167   # round number according to accuracy and precision settings
1168   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1169
1170   return $x if $x->{sign} !~ /^[+-]$/;          # inf, nan etc
1171
1172   my $zeros = $x->{_m}->_trailing_zeros();      # correct for trailing zeros
1173   if ($zeros != 0)
1174     {
1175     $x->{_m}->brsft($zeros,10); $x->{_e} += $zeros;
1176     }
1177   # for something like 0Ey, set y to 1, and -0 => +0
1178   $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
1179   # this is to prevent automatically rounding when MBI's globals are set
1180   $x->{_m}->{_f} = MB_NEVER_ROUND;
1181   $x->{_e}->{_f} = MB_NEVER_ROUND;
1182   # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
1183   $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
1184   $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
1185   return $x;                                    # MBI bnorm is no-op
1186   }
1187
1188 ##############################################################################
1189 # internal calculation routines
1190
1191 sub as_number
1192   {
1193   # return a bigint representation of this BigFloat number
1194   my $x = shift; my $class = ref($x) || $x; $x = $class->new(shift) unless ref($x);
1195
1196   my $z;
1197   if ($x->{_e}->is_zero())
1198     {
1199     $z = $x->{_m}->copy();
1200     $z->{sign} = $x->{sign};
1201     return $z;
1202     }
1203   $z = $x->{_m}->copy();
1204   if ($x->{_e} < 0)
1205     {
1206     $z->brsft(-$x->{_e},10);
1207     }
1208   else
1209     {
1210     $z->blsft($x->{_e},10);
1211     }
1212   $z->{sign} = $x->{sign};
1213   return $z;
1214   }
1215
1216 sub length
1217   {
1218   my $x = shift;
1219   my $class = ref($x) || $x;
1220   $x = $class->new(shift) unless ref($x);
1221
1222   return 1 if $x->{_m}->is_zero();
1223   my $len = $x->{_m}->length();
1224   $len += $x->{_e} if $x->{_e}->sign() eq '+';
1225   if (wantarray())
1226     {
1227     my $t = Math::BigInt::bzero();
1228     $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1229     return ($len,$t);
1230     }
1231   return $len;
1232   }
1233
1234 1;
1235 __END__
1236
1237 =head1 NAME
1238
1239 Math::BigFloat - Arbitrary size floating point math package
1240
1241 =head1 SYNOPSIS
1242
1243   use Math::BigFloat;
1244
1245   # Number creation
1246   $x = Math::BigInt->new($str); # defaults to 0
1247   $nan  = Math::BigInt->bnan(); # create a NotANumber
1248   $zero = Math::BigInt->bzero();# create a "+0"
1249
1250   # Testing
1251   $x->is_zero();                # return whether arg is zero or not
1252   $x->is_nan();                 # return whether arg is NaN or not
1253   $x->is_one();                 # true if arg is +1
1254   $x->is_one('-');              # true if arg is -1
1255   $x->is_odd();                 # true if odd, false for even
1256   $x->is_even();                # true if even, false for odd
1257   $x->is_positive();            # true if >= 0
1258   $x->is_negative();            # true if <  0
1259   $x->is_inf(sign)              # true if +inf or -inf (sign default '+')
1260   $x->bcmp($y);                 # compare numbers (undef,<0,=0,>0)
1261   $x->bacmp($y);                # compare absolutely (undef,<0,=0,>0)
1262   $x->sign();                   # return the sign, either +,- or NaN
1263
1264   # The following all modify their first argument:
1265
1266   # set
1267   $x->bzero();                  # set $i to 0
1268   $x->bnan();                   # set $i to NaN
1269
1270   $x->bneg();                   # negation
1271   $x->babs();                   # absolute value
1272   $x->bnorm();                  # normalize (no-op)
1273   $x->bnot();                   # two's complement (bit wise not)
1274   $x->binc();                   # increment x by 1
1275   $x->bdec();                   # decrement x by 1
1276
1277   $x->badd($y);                 # addition (add $y to $x)
1278   $x->bsub($y);                 # subtraction (subtract $y from $x)
1279   $x->bmul($y);                 # multiplication (multiply $x by $y)
1280   $x->bdiv($y);                 # divide, set $i to quotient
1281                                 # return (quo,rem) or quo if scalar
1282
1283   $x->bmod($y);                 # modulus
1284   $x->bpow($y);                 # power of arguments (a**b)
1285   $x->blsft($y);                # left shift
1286   $x->brsft($y);                # right shift
1287                                 # return (quo,rem) or quo if scalar
1288
1289   $x->band($y);                 # bit-wise and
1290   $x->bior($y);                 # bit-wise inclusive or
1291   $x->bxor($y);                 # bit-wise exclusive or
1292   $x->bnot();                   # bit-wise not (two's complement)
1293
1294   $x->bround($N);               # accuracy: preserver $N digits
1295   $x->bfround($N);              # precision: round to the $Nth digit
1296
1297   # The following do not modify their arguments:
1298
1299   bgcd(@values);                # greatest common divisor
1300   blcm(@values);                # lowest common multiplicator
1301
1302   $x->bstr();                   # return string
1303   $x->bsstr();                  # return string in scientific notation
1304
1305   $x->exponent();               # return exponent as BigInt
1306   $x->mantissa();               # return mantissa as BigInt
1307   $x->parts();                  # return (mantissa,exponent) as BigInt
1308
1309   $x->length();                 # number of digits (w/o sign and '.')
1310   ($l,$f) = $x->length();       # number of digits, and length of fraction
1311
1312 =head1 DESCRIPTION
1313
1314 All operators (inlcuding basic math operations) are overloaded if you
1315 declare your big floating point numbers as
1316
1317   $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1318
1319 Operations with overloaded operators preserve the arguments, which is
1320 exactly what you expect.
1321
1322 =head2 Canonical notation
1323
1324 Input to these routines are either BigFloat objects, or strings of the
1325 following four forms:
1326
1327 =over 2
1328
1329 =item *
1330
1331 C</^[+-]\d+$/>
1332
1333 =item *
1334
1335 C</^[+-]\d+\.\d*$/>
1336
1337 =item *
1338
1339 C</^[+-]\d+E[+-]?\d+$/>
1340
1341 =item *
1342
1343 C</^[+-]\d*\.\d+E[+-]?\d+$/>
1344
1345 =back
1346
1347 all with optional leading and trailing zeros and/or spaces. Additonally,
1348 numbers are allowed to have an underscore between any two digits.
1349
1350 Empty strings as well as other illegal numbers results in 'NaN'.
1351
1352 bnorm() on a BigFloat object is now effectively a no-op, since the numbers
1353 are always stored in normalized form. On a string, it creates a BigFloat
1354 object.
1355
1356 =head2 Output
1357
1358 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1359
1360 The string output will always have leading and trailing zeros stripped and drop
1361 a plus sign. C<bstr()> will give you always the form with a decimal point,
1362 while C<bsstr()> (for scientific) gives you the scientific notation.
1363
1364         Input                   bstr()          bsstr()
1365         '-0'                    '0'             '0E1'
1366         '  -123 123 123'        '-123123123'    '-123123123E0'
1367         '00.0123'               '0.0123'        '123E-4'
1368         '123.45E-2'             '1.2345'        '12345E-4'
1369         '10E+3'                 '10000'         '1E4'
1370
1371 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1372 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1373 return either undef, <0, 0 or >0 and are suited for sort.
1374
1375 Actual math is done by using BigInts to represent the mantissa and exponent.
1376 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
1377 represent the result when input arguments are not numbers, as well as
1378 the result of dividing by zero.
1379
1380 =head2 C<mantissa()>, C<exponent()> and C<parts()>
1381
1382 C<mantissa()> and C<exponent()> return the said parts of the BigFloat
1383 as BigInts such that:
1384
1385         $m = $x->mantissa();
1386         $e = $x->exponent();
1387         $y = $m * ( 10 ** $e );
1388         print "ok\n" if $x == $y;
1389
1390 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
1391
1392 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
1393
1394 Currently the mantissa is reduced as much as possible, favouring higher
1395 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
1396 This might change in the future, so do not depend on it.
1397
1398 =head2 Accuracy vs. Precision
1399
1400 See also: L<Rounding|Rounding>.
1401
1402 Math::BigFloat supports both precision and accuracy. (here should follow
1403 a short description of both).
1404
1405 Precision: digits after the '.', laber, schwad
1406 Accuracy: Significant digits blah blah
1407
1408 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
1409 a operation consumes all resources, each operation produces no more than
1410 C<Math::BigFloat::precision()> digits.
1411
1412 In case the result of one operation has more precision than specified,
1413 it is rounded. The rounding mode taken is either the default mode, or the one
1414 supplied to the operation after the I<scale>:
1415
1416         $x = Math::BigFloat->new(2);
1417         Math::BigFloat::precision(5);           # 5 digits max
1418         $y = $x->copy()->bdiv(3);               # will give 0.66666
1419         $y = $x->copy()->bdiv(3,6);             # will give 0.666666
1420         $y = $x->copy()->bdiv(3,6,'odd');       # will give 0.666667
1421         Math::BigFloat::round_mode('zero');
1422         $y = $x->copy()->bdiv(3,6);             # will give 0.666666
1423
1424 =head2 Rounding
1425
1426 =over 2
1427
1428 =item ffround ( +$scale )
1429
1430 Rounds to the $scale'th place left from the '.', counting from the dot.
1431 The first digit is numbered 1.
1432
1433 =item ffround ( -$scale )
1434
1435 Rounds to the $scale'th place right from the '.', counting from the dot.
1436
1437 =item ffround ( 0 )
1438
1439 Rounds to an integer.
1440
1441 =item fround  ( +$scale )
1442
1443 Preserves accuracy to $scale digits from the left (aka significant digits)
1444 and pads the rest with zeros. If the number is between 1 and -1, the
1445 significant digits count from the first non-zero after the '.'
1446
1447 =item fround  ( -$scale ) and fround ( 0 )
1448
1449 These are effetively no-ops.
1450
1451 =back
1452
1453 All rounding functions take as a second parameter a rounding mode from one of
1454 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
1455
1456 The default rounding mode is 'even'. By using
1457 C<< Math::BigFloat::round_mode($round_mode); >> you can get and set the default
1458 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
1459 no longer supported.
1460 The second parameter to the round functions then overrides the default
1461 temporarily.
1462
1463 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
1464 'trunc' as rounding mode to make it equivalent to:
1465
1466         $x = 2.5;
1467         $y = int($x) + 2;
1468
1469 You can override this by passing the desired rounding mode as parameter to
1470 C<as_number()>:
1471
1472         $x = Math::BigFloat->new(2.5);
1473         $y = $x->as_number('odd');      # $y = 3
1474
1475 =head1 EXAMPLES
1476
1477   # not ready yet
1478
1479 =head1 Autocreating constants
1480
1481 After C<use Math::BigFloat ':constant'> all the floating point constants
1482 in the given scope are converted to C<Math::BigFloat>. This conversion
1483 happens at compile time.
1484
1485 In particular
1486
1487   perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
1488
1489 prints the value of C<2E-100>.  Note that without conversion of
1490 constants the expression 2E-100 will be calculated as normal floating point
1491 number.
1492
1493 =head1 BUGS
1494
1495 =over 2
1496
1497 =item *
1498
1499 The following does not work yet:
1500
1501         $m = $x->mantissa();
1502         $e = $x->exponent();
1503         $y = $m * ( 10 ** $e );
1504         print "ok\n" if $x == $y;
1505
1506 =item *
1507
1508 There is no fmod() function yet.
1509
1510 =back
1511
1512 =head1 CAVEAT
1513
1514 =over 1
1515
1516 =item stringify, bstr()
1517
1518 Both stringify and bstr() now drop the leading '+'. The old code would return
1519 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
1520 reasoning and details.
1521
1522 =item bdiv
1523
1524 The following will probably not do what you expect:
1525
1526         print $c->bdiv(123.456),"\n";
1527
1528 It prints both quotient and reminder since print works in list context. Also,
1529 bdiv() will modify $c, so be carefull. You probably want to use
1530
1531         print $c / 123.456,"\n";
1532         print scalar $c->bdiv(123.456),"\n";  # or if you want to modify $c
1533
1534 instead.
1535
1536 =item Modifying and =
1537
1538 Beware of:
1539
1540         $x = Math::BigFloat->new(5);
1541         $y = $x;
1542
1543 It will not do what you think, e.g. making a copy of $x. Instead it just makes
1544 a second reference to the B<same> object and stores it in $y. Thus anything
1545 that modifies $x will modify $y, and vice versa.
1546
1547         $x->bmul(2);
1548         print "$x, $y\n";       # prints '10, 10'
1549
1550 If you want a true copy of $x, use:
1551
1552         $y = $x->copy();
1553
1554 See also the documentation in L<overload> regarding C<=>.
1555
1556 =item bpow
1557
1558 C<bpow()> now modifies the first argument, unlike the old code which left
1559 it alone and only returned the result. This is to be consistent with
1560 C<badd()> etc. The first will modify $x, the second one won't:
1561
1562         print bpow($x,$i),"\n";         # modify $x
1563         print $x->bpow($i),"\n";        # ditto
1564         print $x ** $i,"\n";            # leave $x alone
1565
1566 =back
1567
1568 =head1 LICENSE
1569
1570 This program is free software; you may redistribute it and/or modify it under
1571 the same terms as Perl itself.
1572
1573 =head1 AUTHORS
1574
1575 Mark Biggar, overloaded interface by Ilya Zakharevich.
1576 Completely rewritten by Tels http://bloodgate.com in 2001.
1577
1578 =cut