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