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