lib/Math/BigInt.pm

   1 #!/usr/bin/perl -w
   2
   3 # mark.biggar@TrustedSysLabs.com
   4 # eay@mincom.com is dead (math::BigInteger)
   5 # see: http://www.cypherspace.org/~adam/rsa/pureperl.html (contacted c. adam
   6 # on 2000/11/13 - but email is dead
   7
   8 # todo:
   9 # - fully remove funky $# stuff (maybe)
  10 # - use integer; vs 1e7 as base
  11 # - speed issues (XS? Bit::Vector?)
  12 # - split out actual math code to Math::BigNumber
  13
  14 # Qs: what exactly happens on numify of HUGE numbers? overflow?
  15 #     $a = -$a is much slower (making copy of $a) than $a->bneg(), hm!?
  16 #     (copy_on_write will help there, but that is not yet implemented)
  17
  18 # The following hash values are used:
  19 #   value: the internal array, base 100000
  20 #   sign : +,-,NaN,+inf,-inf
  21 #   _a   : accuracy
  22 #   _p   : precision
  23 #   _cow : copy on write: number of objects that share the data (NRY)
  24 # Internally the numbers are stored in an array with at least 1 element, no
  25 # leading zero parts (except the first) and in base 100000
  26
  27 # USE_MUL: due to problems on certain os (os390, posix-bc) "* 1e-5" is used
  28 # instead of "/ 1e5" at some places, (marked with USE_MUL). But instead of
  29 # using the reverse only on problematic machines, I used it everytime to avoid
  30 # the costly comparisations. This _should_ work everywhere. Thanx Peter Prymmer
  31
  32 package Math::BigInt;
  33 my $class = "Math::BigInt";
  34
  35 $VERSION = 1.35;
  36 use Exporter;
  37 @ISA =       qw( Exporter );
  38 @EXPORT_OK = qw( bneg babs bcmp badd bmul bdiv bmod bnorm bsub
  39                  bgcd blcm
  40                  bround
  41                  blsft brsft band bior bxor bnot bpow bnan bzero
  42                  bacmp bstr bsstr binc bdec bint binf bfloor bceil
  43                  is_odd is_even is_zero is_one is_nan is_inf sign
  44                  length as_number
  45                  trace objectify _swap
  46                );
  47
  48 #@EXPORT = qw( );
  49 use vars qw/$rnd_mode $accuracy $precision $div_scale/;
  50 use strict;
  51
  52 # Inside overload, the first arg is always an object. If the original code had
  53 # it reversed (like $x = 2 * $y), then the third paramater indicates this
  54 # swapping. To make it work, we use a helper routine which not only reswaps the
  55 # params, but also makes a new object in this case. See _swap() for details,
  56 # especially the cases of operators with different classes.
  57
  58 # For overloaded ops with only one argument we simple use $_[0]->copy() to
  59 # preserve the argument.
  60
  61 # Thus inheritance of overload operators becomes possible and transparent for
  62 # our subclasses without the need to repeat the entire overload section there.
  63
  64 use overload
  65 '='     =>      sub { $_[0]->copy(); },
  66
  67 # '+' and '-' do not use _swap, since it is a triffle slower. If you want to
  68 # override _swap (if ever), then override overload of '+' and '-', too!
  69 # for sub it is a bit tricky to keep b: b-a => -a+b
  70 '-'     =>      sub { my $c = $_[0]->copy; $_[2] ?
  71                    $c->bneg()->badd($_[1]) :
  72                    $c->bsub( $_[1]) },
  73 '+'     =>      sub { $_[0]->copy()->badd($_[1]); },
  74
  75 # some shortcuts for speed (assumes that reversed order of arguments is routed
  76 # to normal '+' and we thus can always modify first arg. If this is changed,
  77 # this breaks and must be adjusted.)
  78 '+='    =>      sub { $_[0]->badd($_[1]); },
  79 '-='    =>      sub { $_[0]->bsub($_[1]); },
  80 '*='    =>      sub { $_[0]->bmul($_[1]); },
  81 '/='    =>      sub { scalar $_[0]->bdiv($_[1]); },
  82 '**='   =>      sub { $_[0]->bpow($_[1]); },
  83
  84 '<=>'   =>      sub { $_[2] ?
  85                       $class->bcmp($_[1],$_[0]) :
  86                       $class->bcmp($_[0],$_[1])},
  87 'cmp'   =>      sub {
  88          $_[2] ?
  89                $_[1] cmp $_[0]->bstr() :
  90                $_[0]->bstr() cmp $_[1] },
  91
  92 'int'   =>      sub { $_[0]->copy(); },
  93 'neg'   =>      sub { $_[0]->copy()->bneg(); },
  94 'abs'   =>      sub { $_[0]->copy()->babs(); },
  95 '~'     =>      sub { $_[0]->copy()->bnot(); },
  96
  97 '*'     =>      sub { my @a = ref($_[0])->_swap(@_); $a[0]->bmul($a[1]); },
  98 '/'     =>      sub { my @a = ref($_[0])->_swap(@_);scalar $a[0]->bdiv($a[1]);},
  99 '%'     =>      sub { my @a = ref($_[0])->_swap(@_); $a[0]->bmod($a[1]); },
 100 '**'    =>      sub { my @a = ref($_[0])->_swap(@_); $a[0]->bpow($a[1]); },
 101 '<<'    =>      sub { my @a = ref($_[0])->_swap(@_); $a[0]->blsft($a[1]); },
 102 '>>'    =>      sub { my @a = ref($_[0])->_swap(@_); $a[0]->brsft($a[1]); },
 103
 104 '&'     =>      sub { my @a = ref($_[0])->_swap(@_); $a[0]->band($a[1]); },
 105 '|'     =>      sub { my @a = ref($_[0])->_swap(@_); $a[0]->bior($a[1]); },
 106 '^'     =>      sub { my @a = ref($_[0])->_swap(@_); $a[0]->bxor($a[1]); },
 107
 108 # can modify arg of ++ and --, so avoid a new-copy for speed, but don't
 109 # use $_[0]->_one(), it modifies $_[0] to be 1!
 110 '++'    =>      sub { $_[0]->binc() },
 111 '--'    =>      sub { $_[0]->bdec() },
 112
 113 # if overloaded, O(1) instead of O(N) and twice as fast for small numbers
 114 'bool'  =>      sub {
 115   # this kludge is needed for perl prior 5.6.0 since returning 0 here fails :-/
 116   # v5.6.1 dumps on that: return !$_[0]->is_zero() || undef;                :-(
 117   my $t = !$_[0]->is_zero();
 118   undef $t if $t == 0;
 119   return $t;
 120   },
 121
 122 qw(
 123 ""      bstr
 124 0+      numify),                # Order of arguments unsignificant
 125 ;
 126
 127 ##############################################################################
 128 # global constants, flags and accessory
 129
 130 # are NaNs ok?
 131 my $NaNOK=1;
 132 # set to 1 for tracing
 133 my $trace = 0;
 134 # constants for easier life
 135 my $nan = 'NaN';
 136 my $BASE_LEN = 5;
 137 my $BASE = int("1e".$BASE_LEN);         # var for trying to change it to 1e7
 138 my $RBASE = 1e-5;                       # see USE_MUL
 139
 140 # Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
 141 $rnd_mode = 'even';
 142 $accuracy = undef;
 143 $precision = undef;
 144 $div_scale = 40;
 145
 146 sub round_mode
 147   {
 148   # make Class->round_mode() work
 149   my $self = shift || $class;
 150   # shift @_ if defined $_[0] && $_[0] eq $class;
 151   if (defined $_[0])
 152     {
 153     my $m = shift;
 154     die "Unknown round mode $m"
 155      if $m !~ /^(even|odd|\+inf|\-inf|zero|trunc)$/;
 156     $rnd_mode = $m; return;
 157     }
 158   return $rnd_mode;
 159   }
 160
 161 sub accuracy
 162   {
 163   # $x->accuracy($a);           ref($x) a
 164   # $x->accuracy();             ref($x);
 165   # Class::accuracy();          # not supported
 166   #print "MBI @_ ($class)\n";
 167   my $x = shift;
 168
 169   die ("accuracy() needs reference to object as first parameter.")
 170    if !ref $x;
 171
 172   if (@_ > 0)
 173     {
 174     $x->{_a} = shift;
 175     $x->round() if defined $x->{_a};
 176     }
 177   return $x->{_a};
 178   }
 179
 180 sub precision
 181   {
 182   my $x = shift;
 183
 184   die ("precision() needs reference to object as first parameter.")
 185    unless ref $x;
 186
 187   if (@_ > 0)
 188     {
 189     $x->{_p} = shift;
 190     $x->round() if defined $x->{_p};
 191     }
 192   return $x->{_p};
 193   }
 194
 195 sub _scale_a
 196   {
 197   # select accuracy parameter based on precedence,
 198   # used by bround() and bfround(), may return undef for scale (means no op)
 199   my ($x,$s,$m,$scale,$mode) = @_;
 200   $scale = $x->{_a} if !defined $scale;
 201   $scale = $s if (!defined $scale);
 202   $mode = $m if !defined $mode;
 203   return ($scale,$mode);
 204   }
 205
 206 sub _scale_p
 207   {
 208   # select precision parameter based on precedence,
 209   # used by bround() and bfround(), may return undef for scale (means no op)
 210   my ($x,$s,$m,$scale,$mode) = @_;
 211   $scale = $x->{_p} if !defined $scale;
 212   $scale = $s if (!defined $scale);
 213   $mode = $m if !defined $mode;
 214   return ($scale,$mode);
 215   }
 216
 217 ##############################################################################
 218 # constructors
 219
 220 sub copy
 221   {
 222   my ($c,$x);
 223   if (@_ > 1)
 224     {
 225     # if two arguments, the first one is the class to "swallow" subclasses
 226     ($c,$x) = @_;
 227     }
 228   else
 229     {
 230     $x = shift;
 231     $c = ref($x);
 232     }
 233   return unless ref($x); # only for objects
 234
 235   my $self = {}; bless $self,$c;
 236   foreach my $k (keys %$x)
 237     {
 238     if (ref($x->{$k}) eq 'ARRAY')
 239       {
 240       $self->{$k} = [ @{$x->{$k}} ];
 241       }
 242     elsif (ref($x->{$k}) eq 'HASH')
 243       {
 244       # only one level deep!
 245       foreach my $h (keys %{$x->{$k}})
 246         {
 247         $self->{$k}->{$h} = $x->{$k}->{$h};
 248         }
 249       }
 250     elsif (ref($x->{$k}))
 251       {
 252       my $c = ref($x->{$k});
 253       $self->{$k} = $c->new($x->{$k}); # no copy() due to deep rec
 254       }
 255     else
 256       {
 257       $self->{$k} = $x->{$k};
 258       }
 259     }
 260   $self;
 261   }
 262
 263 sub new
 264   {
 265   # create a new BigInts object from a string or another bigint object.
 266   # value => internal array representation
 267   # sign  => sign (+/-), or "NaN"
 268
 269   # the argument could be an object, so avoid ||, && etc on it, this would
 270   # cause costly overloaded code to be called. The only allowed op are ref()
 271   # and definend.
 272
 273   trace (@_);
 274   my $class = shift;
 275
 276   my $wanted = shift; # avoid numify call by not using || here
 277   return $class->bzero() if !defined $wanted;   # default to 0
 278   return $class->copy($wanted) if ref($wanted);
 279
 280   my $self = {}; bless $self, $class;
 281   # handle '+inf', '-inf' first
 282   if ($wanted =~ /^[+-]inf$/)
 283     {
 284     $self->{value} = [ 0 ];
 285     $self->{sign} = $wanted;
 286     return $self;
 287     }
 288   # split str in m mantissa, e exponent, i integer, f fraction, v value, s sign
 289   my ($mis,$miv,$mfv,$es,$ev) = _split(\$wanted);
 290   if (ref $mis && !ref $miv)
 291     {
 292     # _from_hex
 293     $self->{value} = $mis->{value};
 294     $self->{sign} = $mis->{sign};
 295     return $self;
 296     }
 297   if (!ref $mis)
 298     {
 299     die "$wanted is not a number initialized to $class" if !$NaNOK;
 300     #print "NaN 1\n";
 301     $self->{value} = [ 0 ];
 302     $self->{sign} = $nan;
 303     return $self;
 304     }
 305   # make integer from mantissa by adjusting exp, then convert to bigint
 306   $self->{sign} = $$mis;                        # store sign
 307   $self->{value} = [ 0 ];                       # for all the NaN cases
 308   my $e = int("$$es$$ev");                      # exponent (avoid recursion)
 309   if ($e > 0)
 310     {
 311     my $diff = $e - CORE::length($$mfv);
 312     if ($diff < 0)                              # Not integer
 313       {
 314       #print "NOI 1\n";
 315       $self->{sign} = $nan;
 316       }
 317     else                                        # diff >= 0
 318       {
 319       # adjust fraction and add it to value
 320       # print "diff > 0 $$miv\n";
 321       $$miv = $$miv . ($$mfv . '0' x $diff);
 322       }
 323     }
 324   else
 325     {
 326     if ($$mfv ne '')                            # e <= 0
 327       {
 328       # fraction and negative/zero E => NOI
 329       #print "NOI 2 \$\$mfv '$$mfv'\n";
 330       $self->{sign} = $nan;
 331       }
 332     elsif ($e < 0)
 333       {
 334       # xE-y, and empty mfv
 335       #print "xE-y\n";
 336       $e = abs($e);
 337       if ($$miv !~ s/0{$e}$//)          # can strip so many zero's?
 338         {
 339         #print "NOI 3\n";
 340         $self->{sign} = $nan;
 341         }
 342       }
 343     }
 344   $self->{sign} = '+' if $$miv eq '0';                  # normalize -0 => +0
 345   $self->_internal($miv) if $self->{sign} ne $nan;      # as internal array
 346   #print "$wanted => $self->{sign} $self->{value}->[0]\n";
 347   # if any of the globals is set, round to them and thus store them insid $self
 348   $self->round($accuracy,$precision,$rnd_mode)
 349    if defined $accuracy || defined $precision;
 350   return $self;
 351   }
 352
 353 # some shortcuts for easier life
 354 sub bint
 355   {
 356   # exportable version of new
 357   trace(@_);
 358   return $class->new(@_);
 359   }
 360
 361 sub bnan
 362   {
 363   # create a bigint 'NaN', if given a BigInt, set it to 'NaN'
 364   my $self = shift;
 365   $self = $class if !defined $self;
 366   if (!ref($self))
 367     {
 368     my $c = $self; $self = {}; bless $self, $c;
 369     }
 370   return if $self->modify('bnan');
 371   $self->{value} = [ 0 ];
 372   $self->{sign} = $nan;
 373   trace('NaN');
 374   return $self;
 375   }
 376
 377 sub binf
 378   {
 379   # create a bigint '+-inf', if given a BigInt, set it to '+-inf'
 380   # the sign is either '+', or if given, used from there
 381   my $self = shift;
 382   my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
 383   $self = $class if !defined $self;
 384   if (!ref($self))
 385     {
 386     my $c = $self; $self = {}; bless $self, $c;
 387     }
 388   return if $self->modify('binf');
 389   $self->{value} = [ 0 ];
 390   $self->{sign} = $sign.'inf';
 391   trace('inf');
 392   return $self;
 393   }
 394
 395 sub bzero
 396   {
 397   # create a bigint '+0', if given a BigInt, set it to 0
 398   my $self = shift;
 399   $self = $class if !defined $self;
 400   if (!ref($self))
 401     {
 402     my $c = $self; $self = {}; bless $self, $c;
 403     }
 404   return if $self->modify('bzero');
 405   $self->{value} = [ 0 ];
 406   $self->{sign} = '+';
 407   trace('0');
 408   return $self;
 409   }
 410
 411 ##############################################################################
 412 # string conversation
 413
 414 sub bsstr
 415   {
 416   # (ref to BFLOAT or num_str ) return num_str
 417   # Convert number from internal format to scientific string format.
 418   # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
 419   trace(@_);
 420   my ($self,$x) = objectify(1,@_);
 421
 422   return $x->{sign} if $x->{sign} !~ /^[+-]$/;
 423   my ($m,$e) = $x->parts();
 424   # can be only '+', so
 425   my $sign = 'e+';
 426   # MBF: my $s = $e->{sign}; $s = '' if $s eq '-'; my $sep = 'e'.$s;
 427   return $m->bstr().$sign.$e->bstr();
 428   }
 429
 430 sub bstr
 431   {
 432   # (ref to BINT or num_str ) return num_str
 433   # Convert number from internal base 100000 format to string format.
 434   # internal format is always normalized (no leading zeros, "-0" => "+0")
 435   trace(@_);
 436   my $x = shift; $x = $class->new($x) unless ref $x;
 437   # my ($self,$x) = objectify(1,@_);
 438
 439   return $x->{sign} if $x->{sign} !~ /^[+-]$/;
 440   my $ar = $x->{value} || return $nan;          # should not happen
 441   my $es = "";
 442   $es = $x->{sign} if $x->{sign} eq '-';        # get sign, but not '+'
 443   my $l = scalar @$ar;         # number of parts
 444   return $nan if $l < 1;       # should not happen
 445   # handle first one different to strip leading zeros from it (there are no
 446   # leading zero parts in internal representation)
 447   $l --; $es .= $ar->[$l]; $l--;
 448   # Interestingly, the pre-padd method uses more time
 449   # the old grep variant takes longer (14 to 10 sec)
 450   while ($l >= 0)
 451     {
 452     $es .= substr('0000'.$ar->[$l],-5);   # fastest way I could think of
 453     $l--;
 454     }
 455   return $es;
 456   }
 457
 458 sub numify
 459   {
 460   # Make a number from a BigInt object
 461   # old: simple return string and let Perl's atoi() handle the rest
 462   # new: calc because it is faster than bstr()+atoi()
 463   #trace (@_);
 464   #my ($self,$x) = objectify(1,@_);
 465   #return $x->bstr(); # ref($x);
 466   my $x = shift; $x = $class->new($x) unless ref $x;
 467
 468   return $nan if $x->{sign} eq $nan;
 469   my $fac = 1; $fac = -1 if $x->{sign} eq '-';
 470   return $fac*$x->{value}->[0] if @{$x->{value}} == 1;  # below $BASE
 471   my $num = 0;
 472   foreach (@{$x->{value}})
 473     {
 474     $num += $fac*$_; $fac *= $BASE;
 475     }
 476   return $num;
 477   }
 478
 479 ##############################################################################
 480 # public stuff (usually prefixed with "b")
 481
 482 sub sign
 483   {
 484   # return the sign of the number: +/-/NaN
 485   my ($self,$x) = objectify(1,@_);
 486   return $x->{sign};
 487   }
 488
 489 sub round
 490   {
 491   # After any operation or when calling round(), the result is rounded by
 492   # regarding the A & P from arguments, local parameters, or globals.
 493   # The result's A or P are set by the rounding, but not inspected beforehand
 494   # (aka only the arguments enter into it). This works because the given
 495   # 'first' argument is both the result and true first argument with unchanged
 496   # A and P settings.
 497   # This does not yet handle $x with A, and $y with P (which should be an
 498   # error).
 499   my $self = shift;
 500   my $a    = shift;     # accuracy, if given by caller
 501   my $p    = shift;     # precision, if given by caller
 502   my $r    = shift;     # round_mode, if given by caller
 503   my @args = @_;        # all 'other' arguments (0 for unary, 1 for binary ops)
 504
 505   unshift @args,$self;  # add 'first' argument
 506
 507   $self = new($self) unless ref($self); # if not object, make one
 508
 509   # find out class of argument to round
 510   my $c = ref($args[0]);
 511
 512   # now pick $a or $p, but only if we have got "arguments"
 513   if ((!defined $a) && (!defined $p) && (@args > 0))
 514     {
 515     foreach (@args)
 516       {
 517       # take the defined one, or if both defined, the one that is smaller
 518       $a = $_->{_a} if (defined $_->{_a}) && (!defined $a || $_->{_a} < $a);
 519       }
 520     if (!defined $a)            # if it still is not defined, take p
 521       {
 522       foreach (@args)
 523         {
 524         # take the defined one, or if both defined, the one that is smaller
 525         $p = $_->{_p} if (defined $_->{_p}) && (!defined $p || $_->{_p} < $p);
 526         }
 527       # if none defined, use globals (#2)
 528       if (!defined $p)
 529         {
 530         no strict 'refs';
 531         my $z = "$c\::accuracy"; $a = $$z;
 532         if (!defined $a)
 533           {
 534           $z = "$c\::precision"; $p = $$z;
 535           }
 536         }
 537       } # endif !$a
 538     } # endif !$a || !$P && args > 0
 539   # for clearity, this is not merged at place (#2)
 540   # now round, by calling fround or ffround:
 541   if (defined $a)
 542     {
 543     $self->{_a} = $a; $self->bround($a,$r);
 544     }
 545   elsif (defined $p)
 546     {
 547     $self->{_p} = $p; $self->bfround($p,$r);
 548     }
 549   return $self->bnorm();
 550   }
 551
 552 sub bnorm
 553   {
 554   # (num_str or BINT) return BINT
 555   # Normalize number -- no-op here
 556   my $self = shift;
 557
 558   return $self;
 559   }
 560
 561 sub babs
 562   {
 563   # (BINT or num_str) return BINT
 564   # make number absolute, or return absolute BINT from string
 565   #my ($self,$x) = objectify(1,@_);
 566   my $x = shift; $x = $class->new($x) unless ref $x;
 567   return $x if $x->modify('babs');
 568   # post-normalized abs for internal use (does nothing for NaN)
 569   $x->{sign} =~ s/^-/+/;
 570   $x;
 571   }
 572
 573 sub bneg
 574   {
 575   # (BINT or num_str) return BINT
 576   # negate number or make a negated number from string
 577   my ($self,$x,$a,$p,$r) = objectify(1,@_);
 578   return $x if $x->modify('bneg');
 579   # for +0 dont negate (to have always normalized)
 580   return $x if $x->is_zero();
 581   $x->{sign} =~ tr/+\-/-+/; # does nothing for NaN
 582   # $x->round($a,$p,$r);        # changing this makes $x - $y modify $y!!
 583   $x;
 584   }
 585
 586 sub bcmp
 587   {
 588   # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
 589   # (BINT or num_str, BINT or num_str) return cond_code
 590   my ($self,$x,$y) = objectify(2,@_);
 591   return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 592   &cmp($x->{value},$y->{value},$x->{sign},$y->{sign}) <=> 0;
 593   }
 594
 595 sub bacmp
 596   {
 597   # Compares 2 values, ignoring their signs.
 598   # Returns one of undef, <0, =0, >0. (suitable for sort)
 599   # (BINT, BINT) return cond_code
 600   my ($self,$x,$y) = objectify(2,@_);
 601   return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 602   acmp($x->{value},$y->{value}) <=> 0;
 603   }
 604
 605 sub badd
 606   {
 607   # add second arg (BINT or string) to first (BINT) (modifies first)
 608   # return result as BINT
 609   trace(@_);
 610   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 611
 612   return $x if $x->modify('badd');
 613   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 614
 615   # for round calls, make array
 616   my @bn = ($a,$p,$r,$y);
 617   # speed: no add for 0+y or x+0
 618   return $x->bnorm(@bn) if $y->is_zero();                       # x+0
 619   if ($x->is_zero())                                            # 0+y
 620     {
 621     # make copy, clobbering up x
 622     $x->{value} = [ @{$y->{value}} ];
 623     $x->{sign} = $y->{sign} || $nan;
 624     return $x->round(@bn);
 625     }
 626
 627   # shortcuts
 628   my $xv = $x->{value};
 629   my $yv = $y->{value};
 630   my ($sx, $sy) = ( $x->{sign}, $y->{sign} ); # get signs
 631
 632   if ($sx eq $sy)
 633     {
 634     add($xv,$yv);                       # if same sign, absolute add
 635     $x->{sign} = $sx;
 636     }
 637   else
 638     {
 639     my $a = acmp ($yv,$xv);             # absolute compare
 640     if ($a > 0)
 641       {
 642       #print "swapped sub (a=$a)\n";
 643       &sub($yv,$xv,1);                  # absolute sub w/ swapped params
 644       $x->{sign} = $sy;
 645       }
 646     elsif ($a == 0)
 647       {
 648       # speedup, if equal, set result to 0
 649       $x->{value} = [ 0 ];
 650       $x->{sign} = '+';
 651       }
 652     else # a < 0
 653       {
 654       #print "unswapped sub (a=$a)\n";
 655       &sub($xv, $yv);                   # absolute sub
 656       $x->{sign} = $sx;
 657       }
 658     }
 659   return $x->round(@bn);
 660   }
 661
 662 sub bsub
 663   {
 664   # (BINT or num_str, BINT or num_str) return num_str
 665   # subtract second arg from first, modify first
 666   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 667
 668   trace(@_);
 669   return $x if $x->modify('bsub');
 670   $x->badd($y->bneg()); # badd does not leave internal zeros
 671   $y->bneg();           # refix y, assumes no one reads $y in between
 672   return $x->round($a,$p,$r,$y);
 673   }
 674
 675 sub binc
 676   {
 677   # increment arg by one
 678   my ($self,$x,$a,$p,$r) = objectify(1,@_);
 679   # my $x = shift; $x = $class->new($x) unless ref $x; my $self = ref($x);
 680   trace(@_);
 681   return $x if $x->modify('binc');
 682   $x->badd($self->_one())->round($a,$p,$r);
 683   }
 684
 685 sub bdec
 686   {
 687   # decrement arg by one
 688   my ($self,$x,$a,$p,$r) = objectify(1,@_);
 689   trace(@_);
 690   return $x if $x->modify('bdec');
 691   $x->badd($self->_one('-'))->round($a,$p,$r);
 692   }
 693
 694 sub blcm
 695   {
 696   # (BINT or num_str, BINT or num_str) return BINT
 697   # does not modify arguments, but returns new object
 698   # Lowest Common Multiplicator
 699   trace(@_);
 700
 701   my ($self,@arg) = objectify(0,@_);
 702   my $x = $self->new(shift @arg);
 703   while (@arg) { $x = _lcm($x,shift @arg); }
 704   $x;
 705   }
 706
 707 sub bgcd
 708   {
 709   # (BINT or num_str, BINT or num_str) return BINT
 710   # does not modify arguments, but returns new object
 711   # GCD -- Euclids algorithm, variant C (Knuth Vol 3, pg 341 ff)
 712   trace(@_);
 713
 714   my ($self,@arg) = objectify(0,@_);
 715   my $x = $self->new(shift @arg);
 716   while (@arg)
 717     {
 718     #$x = _gcd($x,shift @arg); last if $x->is_one(); # new fast, but is slower
 719     $x = _gcd_old($x,shift @arg); last if $x->is_one(); # old, slow, but faster
 720     }
 721   $x;
 722   }
 723
 724 sub bmod
 725   {
 726   # modulus
 727   # (BINT or num_str, BINT or num_str) return BINT
 728   my ($self,$x,$y) = objectify(2,@_);
 729
 730   return $x if $x->modify('bmod');
 731   (&bdiv($self,$x,$y))[1];
 732   }
 733
 734 sub bnot
 735   {
 736   # (num_str or BINT) return BINT
 737   # represent ~x as twos-complement number
 738   my ($self,$x) = objectify(1,@_);
 739   return $x if $x->modify('bnot');
 740   $x->bneg(); $x->bdec(); # was: bsub(-1,$x);, time it someday
 741   $x;
 742   }
 743
 744 sub is_zero
 745   {
 746   # return true if arg (BINT or num_str) is zero (array '+', '0')
 747   #my ($self,$x) = objectify(1,@_);
 748   #trace(@_);
 749   my $x = shift; $x = $class->new($x) unless ref $x;
 750   return (@{$x->{value}} == 1) && ($x->{sign} eq '+')
 751    && ($x->{value}->[0] == 0);
 752   }
 753
 754 sub is_nan
 755   {
 756   # return true if arg (BINT or num_str) is NaN
 757   #my ($self,$x) = objectify(1,@_);
 758   #trace(@_);
 759   my $x = shift; $x = $class->new($x) unless ref $x;
 760   return ($x->{sign} eq $nan);
 761   }
 762
 763 sub is_inf
 764   {
 765   # return true if arg (BINT or num_str) is +-inf
 766   #my ($self,$x) = objectify(1,@_);
 767   #trace(@_);
 768   my $x = shift; $x = $class->new($x) unless ref $x;
 769   my $sign = shift || '';
 770
 771   return $x->{sign} =~ /^[+-]inf/ if $sign eq '';
 772   return $x->{sign} =~ /^[$sign]inf/;
 773   }
 774
 775 sub is_one
 776   {
 777   # return true if arg (BINT or num_str) is +1 (array '+', '1')
 778   # or -1 if signis given
 779   #my ($self,$x) = objectify(1,@_);
 780   my $x = shift; $x = $class->new($x) unless ref $x;
 781   my $sign = shift || '+'; #$_[2] || '+';
 782   return (@{$x->{value}} == 1) && ($x->{sign} eq $sign)
 783    && ($x->{value}->[0] == 1);
 784   }
 785
 786 sub is_odd
 787   {
 788   # return true when arg (BINT or num_str) is odd, false for even
 789   my $x = shift; $x = $class->new($x) unless ref $x;
 790   #my ($self,$x) = objectify(1,@_);
 791   return (($x->{sign} ne $nan) && ($x->{value}->[0] & 1));
 792   }
 793
 794 sub is_even
 795   {
 796   # return true when arg (BINT or num_str) is even, false for odd
 797   my $x = shift; $x = $class->new($x) unless ref $x;
 798   #my ($self,$x) = objectify(1,@_);
 799   return (($x->{sign} ne $nan) && (!($x->{value}->[0] & 1)));
 800   }
 801
 802 sub bmul
 803   {
 804   # multiply two numbers -- stolen from Knuth Vol 2 pg 233
 805   # (BINT or num_str, BINT or num_str) return BINT
 806   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 807   #print "$self bmul $x ",ref($x)," $y ",ref($y),"\n";
 808   trace(@_);
 809   return $x if $x->modify('bmul');
 810   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
 811
 812   mul($x,$y);  # do actual math
 813   return $x->round($a,$p,$r,$y);
 814   }
 815
 816 sub bdiv
 817   {
 818   # (dividend: BINT or num_str, divisor: BINT or num_str) return
 819   # (BINT,BINT) (quo,rem) or BINT (only rem)
 820   trace(@_);
 821   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 822
 823   return $x if $x->modify('bdiv');
 824
 825   # NaN?
 826   return wantarray ? ($x->bnan(),bnan()) : $x->bnan()
 827    if ($x->{sign} eq $nan || $y->{sign} eq $nan || $y->is_zero());
 828
 829   # 0 / something
 830   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
 831
 832   # Is $x in the interval [0, $y) ?
 833   my $cmp = acmp($x->{value},$y->{value});
 834   if (($cmp < 0) and ($x->{sign} eq $y->{sign}))
 835     {
 836     return $x->bzero() unless wantarray;
 837     my $t = $x->copy();      # make copy first, because $x->bzero() clobbers $x
 838     return ($x->bzero(),$t);
 839     }
 840   elsif ($cmp == 0)
 841     {
 842     # shortcut, both are the same, so set to +/- 1
 843     $x->_one( ($x->{sign} ne $y->{sign} ? '-' : '+') );
 844     return $x unless wantarray;
 845     return ($x,$self->bzero());
 846     }
 847
 848   # calc new sign and in case $y == +/- 1, return $x
 849   $x->{sign} = ($x->{sign} ne $y->{sign} ? '-' : '+');
 850   # check for / +-1 (cant use $y->is_one due to '-'
 851   if ((@{$y->{value}} == 1) && ($y->{value}->[0] == 1))
 852     {
 853     return wantarray ? ($x,$self->bzero()) : $x;
 854     }
 855
 856   # call div here
 857   my $rem = $self->bzero();
 858   $rem->{sign} = $y->{sign};
 859   ($x->{value},$rem->{value}) = div($x->{value},$y->{value});
 860   # do not leave rest "-0";
 861   $rem->{sign} = '+' if (@{$rem->{value}} == 1) && ($rem->{value}->[0] == 0);
 862   if (($x->{sign} eq '-') and (!$rem->is_zero()))
 863     {
 864     $x->bdec();
 865     }
 866   $x->round($a,$p,$r,$y);
 867   if (wantarray)
 868     {
 869     $rem->round($a,$p,$r,$x,$y);
 870     return ($x,$y-$rem) if $x->{sign} eq '-';   # was $x,$rem
 871     return ($x,$rem);
 872     }
 873   return $x;
 874   }
 875
 876 sub bpow
 877   {
 878   # (BINT or num_str, BINT or num_str) return BINT
 879   # compute power of two numbers -- stolen from Knuth Vol 2 pg 233
 880   # modifies first argument
 881   #trace(@_);
 882   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
 883
 884   return $x if $x->modify('bpow');
 885
 886   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
 887   return $x->_one() if $y->is_zero();
 888   return $x         if $x->is_one() || $y->is_one();
 889   if ($x->{sign} eq '-' && @{$x->{value}} == 1 && $x->{value}->[0] == 1)
 890     {
 891     # if $x == -1 and odd/even y => +1/-1
 892     return $y->is_odd() ? $x : $x->_set(1); # $x->babs() would work to
 893     # my Casio FX-5500L has here a bug, -1 ** 2 is -1, but -1 * -1 is 1 LOL
 894     }
 895   # shortcut for $x ** 2
 896   if ($y->{sign} eq '+' && @{$y->{value}} == 1 && $y->{value}->[0] == 2)
 897     {
 898     return $x->bmul($x)->bround($a,$p,$r);
 899     }
 900   # 1 ** -y => 1 / (1**y), so do test for negative $y after above's clause
 901   return $x->bnan() if $y->{sign} eq '-';
 902   return $x         if $x->is_zero();  # 0**y => 0 (if not y <= 0)
 903
 904   # tels: 10**x is special (actually 100**x etc is special, too) but not here
 905   #if ((@{$x->{value}} == 1) && ($x->{value}->[0] == 10))
 906   #  {
 907   #  # 10**2
 908   #  my $yi = int($y); my $yi5 = int($yi/5);
 909   #  $x->{value} = [];
 910   #  my $v = $x->{value};
 911   #  if ($yi5 > 0)
 912   #    {
 913   #    # $x->{value}->[$yi5-1] = 0;             # pre-padd array (no use)
 914   #    for (my $i = 0; $i < $yi5; $i++)
 915   #      {
 916   #      $v->[$i] = 0;
 917   #      }
 918   #    }
 919   #  push @{$v}, int( '1'.'0' x ($yi % 5));
 920   #  if ($x->{sign} eq '-')
 921   #    {
 922   #    $x->{sign} = $y->is_odd() ? '-' : '+';   # -10**2 = 100, -10**3 = -1000
 923   #    }
 924   #  return $x;
 925   #  }
 926
 927   # based on the assumption that shifting in base 10 is fast, and that bpow()
 928   # works faster if numbers are small: we count trailing zeros (this step is
 929   # O(1)..O(N), but in case of O(N) we save much more time), stripping them
 930   # out of the multiplication, and add $count * $y zeros afterwards:
 931   # 300 ** 3 == 300*300*300 == 3*3*3 . '0' x 2 * 3 == 27 . '0' x 6
 932   my $zeros = $x->_trailing_zeros();
 933   if ($zeros > 0)
 934     {
 935     $x->brsft($zeros,10);       # remove zeros
 936     $x->bpow($y);               # recursion (will not branch into here again)
 937     $zeros = $y * $zeros;       # real number of zeros to add
 938     $x->blsft($zeros,10);
 939     return $x;
 940     }
 941
 942   my $pow2 = $self->_one();
 943   my $y1 = $class->new($y);
 944   my ($res);
 945   while (!$y1->is_one())
 946     {
 947     #print "bpow: p2: $pow2 x: $x y: $y1 r: $res\n";
 948     #print "len ",$x->length(),"\n";
 949     ($y1,$res)=&bdiv($y1,2);
 950     if (!$res->is_zero()) { &bmul($pow2,$x); }
 951     if (!$y1->is_zero())  { &bmul($x,$x); }
 952     }
 953   #print "bpow: e p2: $pow2 x: $x y: $y1 r: $res\n";
 954   &bmul($x,$pow2) if (!$pow2->is_one());
 955   #print "bpow: e p2: $pow2 x: $x y: $y1 r: $res\n";
 956   return $x->round($a,$p,$r);
 957   }
 958
 959 sub blsft
 960   {
 961   # (BINT or num_str, BINT or num_str) return BINT
 962   # compute x << y, base n, y >= 0
 963   my ($self,$x,$y,$n) = objectify(2,@_);
 964
 965   return $x if $x->modify('blsft');
 966   return $x->bnan() if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/);
 967
 968   $n = 2 if !defined $n; return $x if $n == 0;
 969   return $x->bnan() if $n < 0 || $y->{sign} eq '-';
 970   if ($n != 10)
 971     {
 972     $x->bmul( $self->bpow($n, $y) );
 973     }
 974   else
 975     {
 976     # shortcut (faster) for shifting by 10) since we are in base 10eX
 977     # multiples of 5:
 978     my $src = scalar @{$x->{value}};            # source
 979     my $len = $y->numify();                     # shift-len as normal int
 980     my $rem = $len % 5;                         # reminder to shift
 981     my $dst = $src + int($len/5);               # destination
 982
 983     my $v = $x->{value};                        # speed-up
 984     my $vd;                                     # further speedup
 985     #print "src $src:",$v->[$src]||0," dst $dst:",$v->[$dst]||0," rem $rem\n";
 986     $v->[$src] = 0;                             # avoid first ||0 for speed
 987     while ($src >= 0)
 988       {
 989       $vd = $v->[$src]; $vd = '00000'.$vd;
 990       #print "s $src d $dst '$vd' ";
 991       $vd = substr($vd,-5+$rem,5-$rem);
 992       #print "'$vd' ";
 993       $vd .= $src > 0 ? substr('00000'.$v->[$src-1],-5,$rem) : '0' x $rem;
 994       #print "'$vd' ";
 995       $vd = substr($vd,-5,5) if length($vd) > 5;
 996       #print "'$vd'\n";
 997       $v->[$dst] = int($vd);
 998       $dst--; $src--;
 999       }
1000     # set lowest parts to 0
1001     while ($dst >= 0) { $v->[$dst--] = 0; }
1002     # fix spurios last zero element
1003     splice @$v,-1 if $v->[-1] == 0;
1004     #print "elems: "; my $i = 0;
1005     #foreach (reverse @$v) { print "$i $_ "; $i++; } print "\n";
1006     # old way: $x->bmul( $self->bpow($n, $y) );
1007     }
1008   return $x;
1009   }
1010
1011 sub brsft
1012   {
1013   # (BINT or num_str, BINT or num_str) return BINT
1014   # compute x >> y, base n, y >= 0
1015   my ($self,$x,$y,$n) = objectify(2,@_);
1016
1017   return $x if $x->modify('brsft');
1018   return $x->bnan() if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/);
1019
1020   $n = 2 if !defined $n; return $x->bnan() if $n <= 0 || $y->{sign} eq '-';
1021   if ($n != 10)
1022     {
1023     scalar bdiv($x, $self->bpow($n, $y));
1024     }
1025   else
1026     {
1027     # shortcut (faster) for shifting by 10)
1028     # multiples of 5:
1029     my $dst = 0;                                # destination
1030     my $src = $y->numify();                     # as normal int
1031     my $rem = $src % 5;                         # reminder to shift
1032     $src = int($src / 5);                       # source
1033     my $len = scalar @{$x->{value}} - $src;     # elems to go
1034     my $v = $x->{value};                        # speed-up
1035     if ($rem == 0)
1036       {
1037       splice (@$v,0,$src);                      # even faster, 38.4 => 39.3
1038       }
1039     else
1040       {
1041       my $vd;
1042       $v->[scalar @$v] = 0;                     # avoid || 0 test inside loop
1043       while ($dst < $len)
1044         {
1045         $vd = '00000'.$v->[$src];
1046         #print "$dst $src '$vd' ";
1047         $vd = substr($vd,-5,5-$rem);
1048         #print "'$vd' ";
1049         $src++;
1050         $vd = substr('00000'.$v->[$src],-$rem,$rem) . $vd;
1051         #print "'$vd1' ";
1052         #print "'$vd'\n";
1053         $vd = substr($vd,-5,5) if length($vd) > 5;
1054         $v->[$dst] = int($vd);
1055         $dst++;
1056         }
1057       splice (@$v,$dst) if $dst > 0;            # kill left-over array elems
1058       pop @$v if $v->[-1] == 0;                 # kill last element
1059       } # else rem == 0
1060     # old way: scalar bdiv($x, $self->bpow($n, $y));
1061     }
1062   return $x;
1063   }
1064
1065 sub band
1066   {
1067   #(BINT or num_str, BINT or num_str) return BINT
1068   # compute x & y
1069   trace(@_);
1070   my ($self,$x,$y) = objectify(2,@_);
1071
1072   return $x if $x->modify('band');
1073
1074   return $x->bnan() if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/);
1075   return $x->bzero() if $y->is_zero();
1076   my $r = $self->bzero(); my $m = new Math::BigInt 1; my ($xr,$yr);
1077   my $x10000 = new Math::BigInt (0x10000);
1078   my $y1 = copy(ref($x),$y);                            # make copy
1079   while (!$x->is_zero() && !$y1->is_zero())
1080     {
1081     ($x, $xr) = bdiv($x, $x10000);
1082     ($y1, $yr) = bdiv($y1, $x10000);
1083     $r->badd( bmul( new Math::BigInt ( $xr->numify() & $yr->numify()), $m ));
1084     $m->bmul($x10000);
1085     }
1086   $x = $r;
1087   }
1088
1089 sub bior
1090   {
1091   #(BINT or num_str, BINT or num_str) return BINT
1092   # compute x | y
1093   trace(@_);
1094   my ($self,$x,$y) = objectify(2,@_);
1095
1096   return $x if $x->modify('bior');
1097
1098   return $x->bnan() if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/);
1099   return $x if $y->is_zero();
1100   my $r = $self->bzero(); my $m = new Math::BigInt 1; my ($xr,$yr);
1101   my $x10000 = new Math::BigInt (0x10000);
1102   my $y1 = copy(ref($x),$y);                            # make copy
1103   while (!$x->is_zero() || !$y1->is_zero())
1104     {
1105     ($x, $xr) = bdiv($x,$x10000);
1106     ($y1, $yr) = bdiv($y1,$x10000);
1107     $r->badd( bmul( new Math::BigInt ( $xr->numify() | $yr->numify()), $m ));
1108     $m->bmul($x10000);
1109     }
1110   $x = $r;
1111   }
1112
1113 sub bxor
1114   {
1115   #(BINT or num_str, BINT or num_str) return BINT
1116   # compute x ^ y
1117   my ($self,$x,$y) = objectify(2,@_);
1118
1119   return $x if $x->modify('bxor');
1120
1121   return $x->bnan() if ($x->{sign} eq $nan || $y->{sign} eq $nan);
1122   return $x if $y->is_zero();
1123   return $x->bzero() if $x == $y; # shortcut
1124   my $r = $self->bzero(); my $m = new Math::BigInt 1; my ($xr,$yr);
1125   my $x10000 = new Math::BigInt (0x10000);
1126   my $y1 = copy(ref($x),$y);                    # make copy
1127   while (!$x->is_zero() || !$y1->is_zero())
1128     {
1129     ($x, $xr) = bdiv($x, $x10000);
1130     ($y1, $yr) = bdiv($y1, $x10000);
1131     $r->badd( bmul( new Math::BigInt ( $xr->numify() ^ $yr->numify()), $m ));
1132     $m->bmul($x10000);
1133     }
1134   $x = $r;
1135   }
1136
1137 sub length
1138   {
1139   my ($self,$x) = objectify(1,@_);
1140
1141   return (_digits($x->{value}), 0) if wantarray;
1142   _digits($x->{value});
1143   }
1144
1145 sub digit
1146   {
1147   # return the nth digit, negative values count backward
1148   my $x = shift;
1149   my $n = shift || 0;
1150
1151   my $len = $x->length();
1152
1153   $n = $len+$n if $n < 0;               # -1 last, -2 second-to-last
1154   $n = abs($n);                         # if negatives are to big
1155   $len--; $n = $len if $n > $len;       # n to big?
1156
1157   my $elem = int($n / 5);               # which array element
1158   my $digit = $n % 5;                   # which digit in this element
1159   $elem = '0000'.$x->{value}->[$elem];  # get element padded with 0's
1160   return substr($elem,-$digit-1,1);
1161   }
1162
1163 sub _trailing_zeros
1164   {
1165   # return the amount of trailing zeros in $x
1166   my $x = shift;
1167   $x = $class->new($x) unless ref $x;
1168
1169   return 0 if $x->is_zero() || $x->is_nan();
1170   # check each array elem in _m for having 0 at end as long as elem == 0
1171   # Upon finding a elem != 0, stop
1172   my $zeros = 0; my $elem;
1173   foreach my $e (@{$x->{value}})
1174     {
1175     if ($e != 0)
1176       {
1177       $elem = "$e";                             # preserve x
1178       $elem =~ s/.*?(0*$)/$1/;                  # strip anything not zero
1179       $zeros *= 5;                              # elems * 5
1180       $zeros += CORE::length($elem);            # count trailing zeros
1181       last;                                     # early out
1182       }
1183     $zeros ++;                                  # real else branch: 50% slower!
1184     }
1185   return $zeros;
1186   }
1187
1188 sub bsqrt
1189   {
1190   my ($self,$x) = objectify(1,@_);
1191
1192   return $x->bnan() if $x->{sign} =~ /\-|$nan/; # -x or NaN => NaN
1193   return $x->bzero() if $x->is_zero();          # 0 => 0
1194   return $x if $x == 1;                         # 1 => 1
1195
1196   my $y = $x->copy();                           # give us one more digit accur.
1197   my $l = int($x->length()/2);
1198
1199   $x->bzero();
1200   $x->binc();           # keep ref($x), but modify it
1201   $x *= 10 ** $l;
1202
1203   # print "x: $y guess $x\n";
1204
1205   my $last = $self->bzero();
1206   while ($last != $x)
1207     {
1208     $last = $x;
1209     $x += $y / $x;
1210     $x /= 2;
1211     }
1212   return $x;
1213   }
1214
1215 sub exponent
1216   {
1217   # return a copy of the exponent (here always 0, NaN or 1 for $m == 0)
1218   my ($self,$x) = objectify(1,@_);
1219
1220   return bnan() if $x->is_nan();
1221   my $e = $class->bzero();
1222   return $e->binc() if $x->is_zero();
1223   $e += $x->_trailing_zeros();
1224   return $e;
1225   }
1226
1227 sub mantissa
1228   {
1229   # return a copy of the mantissa (here always $self)
1230   my ($self,$x) = objectify(1,@_);
1231
1232   return bnan() if $x->is_nan();
1233   my $m = $x->copy();
1234   # that's inefficient
1235   my $zeros = $m->_trailing_zeros();
1236   $m /= 10 ** $zeros if $zeros != 0;
1237   return $m;
1238   }
1239
1240 sub parts
1241   {
1242   # return a copy of both the exponent and the mantissa (here 0 and self)
1243   my $self = shift;
1244   $self = $class->new($self) unless ref $self;
1245
1246   return ($self->mantissa(),$self->exponent());
1247   }
1248
1249 ##############################################################################
1250 # rounding functions
1251
1252 sub bfround
1253   {
1254   # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1255   # $n == 0 => round to integer
1256   my $x = shift; $x = $class->new($x) unless ref $x;
1257   my ($scale,$mode) = $x->_scale_p($precision,$rnd_mode,@_);
1258   return $x if !defined $scale;         # no-op
1259
1260   # no-op for BigInts if $n <= 0
1261   return $x if $scale <= 0;
1262
1263   $x->bround( $x->length()-$scale, $mode);
1264   }
1265
1266 sub _scan_for_nonzero
1267   {
1268   my $x = shift;
1269   my $pad = shift;
1270
1271   my $len = $x->length();
1272   return 0 if $len == 1;                # '5' is trailed by invisible zeros
1273   my $follow = $pad - 1;
1274   return 0 if $follow > $len || $follow < 1;
1275   #print "checking $x $r\n";
1276   # old, slow way checking string for non-zero characters
1277   my $r = substr ("$x",-$follow);
1278   return 1 if $r =~ /[^0]/; return 0;
1279
1280   # faster way checking array contents; it is actually not faster (even in a
1281   # rounding-only-shoutout, so I leave the simpler code in)
1282   #my $rem = $follow % 5; my $div = $follow / 5; my $v = $x->{value};
1283   # pad with zeros and extract
1284   #print "last part : ",'00000'.$v->[$div]," $rem = '";
1285   #print substr('00000'.$v->[$div],-$rem,5),"'\n";
1286   #my $r1 = substr ('00000'.$v->[$div],-$rem,5);
1287   #print "$r1\n";
1288   #return 1 if $r1 =~ /[^0]/;
1289   #
1290   #for (my $j = $div-1; $j >= 0; $j --)
1291   #  {
1292   #  #print "part $v->[$j]\n";
1293   #  return 1 if $v->[$j] != 0;
1294   #  }
1295   #return 0;
1296   }
1297
1298 sub fround
1299   {
1300   # to make life easier for switch between MBF and MBI (autoload fxxx()
1301   # like MBF does for bxxx()?)
1302   my $x = shift;
1303   return $x->bround(@_);
1304   }
1305
1306 sub bround
1307   {
1308   # accuracy: +$n preserve $n digits from left,
1309   #           -$n preserve $n digits from right (f.i. for 0.1234 style in MBF)
1310   # no-op for $n == 0
1311   # and overwrite the rest with 0's, return normalized number
1312   # do not return $x->bnorm(), but $x
1313   my $x = shift; $x = $class->new($x) unless ref $x;
1314   my ($scale,$mode) = $x->_scale_a($accuracy,$rnd_mode,@_);
1315   return $x if !defined $scale;         # no-op
1316
1317   # print "MBI round: $x to $scale $mode\n";
1318   # -scale means what? tom? hullo? -$scale needed by MBF round, but what for?
1319   return $x if $x->is_nan() || $x->is_zero() || $scale == 0;
1320
1321   # we have fewer digits than we want to scale to
1322   my $len = $x->length();
1323   # print "$len $scale\n";
1324   return $x if $len < abs($scale);
1325
1326   # count of 0's to pad, from left (+) or right (-): 9 - +6 => 3, or |-6| => 6
1327   my ($pad,$digit_round,$digit_after);
1328   $pad = $len - $scale;
1329   $pad = abs($scale)+1 if $scale < 0;
1330   $digit_round = '0'; $digit_round = $x->digit($pad) if $pad < $len;
1331   $digit_after = '0'; $digit_after = $x->digit($pad-1) if $pad > 0;
1332   # print "r $x: pos:$pad l:$len s:$scale r:$digit_round a:$digit_after m: $mode\n";
1333
1334   # in case of 01234 we round down, for 6789 up, and only in case 5 we look
1335   # closer at the remaining digits of the original $x, remember decision
1336   my $round_up = 1;                                     # default round up
1337   $round_up -- if
1338     ($mode eq 'trunc')                          ||      # trunc by round down
1339     ($digit_after =~ /[01234]/)                 ||      # round down anyway,
1340                                                         # 6789 => round up
1341     ($digit_after eq '5')                       &&      # not 5000...0000
1342     ($x->_scan_for_nonzero($pad) == 0)          &&
1343     (
1344      ($mode eq 'even') && ($digit_round =~ /[24680]/) ||
1345      ($mode eq 'odd')  && ($digit_round =~ /[13579]/) ||
1346      ($mode eq '+inf') && ($x->{sign} eq '-')   ||
1347      ($mode eq '-inf') && ($x->{sign} eq '+')   ||
1348      ($mode eq 'zero')          # round down if zero, sign adjusted below
1349     );
1350   # allow rounding one place left of mantissa
1351   #print "$pad $len $scale\n";
1352   # this is triggering warnings, and buggy for $scale < 0
1353   #if (-$scale != $len)
1354     {
1355     # split mantissa at $scale and then pad with zeros
1356     my $s5 = int($pad / 5);
1357     my $i = 0;
1358     while ($i < $s5)
1359       {
1360       $x->{value}->[$i++] = 0;                          # replace with 5 x 0
1361       }
1362     $x->{value}->[$s5] = '00000'.$x->{value}->[$s5];    # pad with 0
1363     my $rem = $pad % 5;                                 # so much left over
1364     if ($rem > 0)
1365       {
1366       #print "remainder $rem\n";
1367       #print "elem      $x->{value}->[$s5]\n";
1368       substr($x->{value}->[$s5],-$rem,$rem) = '0' x $rem;       # stamp w/ '0'
1369       }
1370     $x->{value}->[$s5] = int ($x->{value}->[$s5]);      # str '05' => int '5'
1371     }
1372   if ($round_up)                                        # what gave test above?
1373     {
1374     $pad = $len if $scale < 0;                          # tlr: whack 0.51=>1.0
1375     # modify $x in place, undef, undef to avoid rounding
1376     $x->badd( Math::BigInt->new($x->{sign}.'1'.'0'x$pad),
1377      undef,undef );
1378     # str creation much faster than 10 ** something
1379     }
1380   $x;
1381   }
1382
1383 sub bfloor
1384   {
1385   # return integer less or equal then number, since it is already integer,
1386   # always returns $self
1387   my ($self,$x,$a,$p,$r) = objectify(1,@_);
1388
1389   # not needed: return $x if $x->modify('bfloor');
1390
1391   return $x->round($a,$p,$r);
1392   }
1393
1394 sub bceil
1395   {
1396   # return integer greater or equal then number, since it is already integer,
1397   # always returns $self
1398   my ($self,$x,$a,$p,$r) = objectify(1,@_);
1399
1400   # not needed: return $x if $x->modify('bceil');
1401
1402   return $x->round($a,$p,$r);
1403   }
1404
1405 ##############################################################################
1406 # private stuff (internal use only)
1407
1408 sub trace
1409   {
1410   # print out a number without using bstr (avoid deep recurse) for trace/debug
1411   return unless $trace;
1412
1413   my ($package,$file,$line,$sub) = caller(1);
1414   print "'$sub' called from '$package' line $line:\n ";
1415
1416   foreach my $x (@_)
1417     {
1418     if (!defined $x)
1419       {
1420       print "undef, "; next;
1421       }
1422     if (!ref($x))
1423       {
1424       print "'$x' "; next;
1425       }
1426     next if (ref($x) ne "HASH");
1427     print "$x->{sign} ";
1428     foreach (@{$x->{value}})
1429       {
1430       print "$_ ";
1431       }
1432     print ", ";
1433     }
1434   print "\n";
1435   }
1436
1437 sub _set
1438   {
1439   # internal set routine to set X fast to an integer value < [+-]100000
1440   my $self = shift;
1441   my $wanted = shift || 0;
1442
1443   $self->{sign} = $nan, return if $wanted !~ /^[+-]?[0-9]+$/;
1444   $self->{sign} = '-'; $self->{sign} = '+' if $wanted >= 0;
1445   $self->{value} = [ abs($wanted) ];
1446   return $self;
1447   }
1448
1449 sub _one
1450   {
1451   # internal speedup, set argument to 1, or create a +/- 1
1452   my $self = shift;
1453   my $x = $self->bzero(); $x->{value} = [ 1 ]; $x->{sign} = shift || '+'; $x;
1454   }
1455
1456 sub _swap
1457   {
1458   # Overload will swap params if first one is no object ref so that the first
1459   # one is always an object ref. In this case, third param is true.
1460   # This routine is to overcome the effect of scalar,$object creating an object
1461   # of the class of this package, instead of the second param $object. This
1462   # happens inside overload, when the overload section of this package is
1463   # inherited by sub classes.
1464   # For overload cases (and this is used only there), we need to preserve the
1465   # args, hence the copy().
1466   # You can override this method in a subclass, the overload section will call
1467   # $object->_swap() to make sure it arrives at the proper subclass, with some
1468   # exceptions like '+' and '-'.
1469
1470   # object, (object|scalar) => preserve first and make copy
1471   # scalar, object          => swapped, re-swap and create new from first
1472   #                            (using class of second object, not $class!!)
1473   my $self = shift;                     # for override in subclass
1474   #print "swap $self 0:$_[0] 1:$_[1] 2:$_[2]\n";
1475   if ($_[2])
1476     {
1477     my $c = ref ($_[0]) || $class;      # fallback $class should not happen
1478     return ( $c->new($_[1]), $_[0] );
1479     }
1480   else
1481     {
1482     return ( $_[0]->copy(), $_[1] );
1483     }
1484   }
1485
1486 sub objectify
1487   {
1488   # check for strings, if yes, return objects instead
1489
1490   # the first argument is number of args objectify() should look at it will
1491   # return $count+1 elements, the first will be a classname. This is because
1492   # overloaded '""' calls bstr($object,undef,undef) and this would result in
1493   # useless objects beeing created and thrown away. So we cannot simple loop
1494   # over @_. If the given count is 0, all arguments will be used.
1495
1496   # If the second arg is a ref, use it as class.
1497   # If not, try to use it as classname, unless undef, then use $class
1498   # (aka Math::BigInt). The latter shouldn't happen,though.
1499
1500   # caller:                        gives us:
1501   # $x->badd(1);                => ref x, scalar y
1502   # Class->badd(1,2);           => classname x (scalar), scalar x, scalar y
1503   # Class->badd( Class->(1),2); => classname x (scalar), ref x, scalar y
1504   # Math::BigInt::badd(1,2);    => scalar x, scalar y
1505   # In the last case we check number of arguments to turn it silently into
1506   # $class,1,2. (We can not take '1' as class ;o)
1507   # badd($class,1) is not supported (it should, eventually, try to add undef)
1508   # currently it tries 'Math::BigInt' + 1, which will not work.
1509
1510   trace(@_);
1511   my $count = abs(shift || 0);
1512
1513   #print caller(),"\n";
1514
1515   my @a;                        # resulting array
1516   if (ref $_[0])
1517     {
1518     # okay, got object as first
1519     $a[0] = ref $_[0];
1520     }
1521   else
1522     {
1523     # nope, got 1,2 (Class->xxx(1) => Class,1 and not supported)
1524     $a[0] = $class;
1525     #print "@_\n"; sleep(1);
1526     $a[0] = shift if $_[0] =~ /^[A-Z].*::/;     # classname as first?
1527     }
1528   #print caller(),"\n";
1529   # print "Now in objectify, my class is today $a[0]\n";
1530   my $k;
1531   if ($count == 0)
1532     {
1533     while (@_)
1534       {
1535       $k = shift;
1536       if (!ref($k))
1537         {
1538         $k = $a[0]->new($k);
1539         }
1540       elsif (ref($k) ne $a[0])
1541         {
1542         # foreign object, try to convert to integer
1543         $k->can('as_number') ?  $k = $k->as_number() : $k = $a[0]->new($k);
1544         }
1545       push @a,$k;
1546       }
1547     }
1548   else
1549     {
1550     while ($count > 0)
1551       {
1552       #print "$count\n";
1553       $count--;
1554       $k = shift;
1555       if (!ref($k))
1556         {
1557         $k = $a[0]->new($k);
1558         }
1559       elsif (ref($k) ne $a[0])
1560         {
1561         # foreign object, try to convert to integer
1562         $k->can('as_number') ?  $k = $k->as_number() : $k = $a[0]->new($k);
1563         }
1564       push @a,$k;
1565       }
1566     push @a,@_;         # return other params, too
1567     }
1568   #my $i = 0;
1569   #foreach (@a)
1570   #  {
1571   #  print "o $i $a[0]\n" if $i == 0;
1572   #  print "o $i ",ref($_),"\n" if $i != 0; $i++;
1573   #  }
1574   #print "objectify done: would return ",scalar @a," values\n";
1575   #print caller(1),"\n" unless wantarray;
1576   die "$class objectify needs list context" unless wantarray;
1577   @a;
1578   }
1579
1580 sub import
1581   {
1582   my $self = shift;
1583   #print "import $self @_\n";
1584   for ( my $i = 0; $i < @_ ; $i++ )
1585     {
1586     if ( $_[$i] eq ':constant' )
1587       {
1588       # this rest causes overlord er load to step in
1589       overload::constant integer => sub { $self->new(shift) };
1590       splice @_, $i, 1; last;
1591       }
1592     }
1593   # any non :constant stuff is handled by our parent, Exporter
1594   # even if @_ is empty, to give it a chance
1595   #$self->SUPER::import(@_);                    # does not work
1596   $self->export_to_level(1,$self,@_);           # need this instead
1597   }
1598
1599 sub _internal
1600   {
1601   # (ref to self, ref to string) return ref to num_array
1602   # Convert a number from string format to internal base 100000 format.
1603   # Assumes normalized value as input.
1604   my ($s,$d) = @_;
1605   my $il = CORE::length($$d)-1;
1606   # these leaves '00000' instead of int 0 and will be corrected after any op
1607   $s->{value} = [ reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $$d)) ];
1608   $s;
1609   }
1610
1611 sub _strip_zeros
1612   {
1613   # internal normalization function that strips leading zeros from the array
1614   # args: ref to array
1615   #trace(@_);
1616   my $s = shift;
1617
1618   my $cnt = scalar @$s; # get count of parts
1619   my $i = $cnt-1;
1620   #print "strip: cnt $cnt i $i\n";
1621   # '0', '3', '4', '0', '0',
1622   #  0    1    2    3    4
1623   # cnt = 5, i = 4
1624   # i = 4
1625   # i = 3
1626   # => fcnt = cnt - i (5-2 => 3, cnt => 5-1 = 4, throw away from 4th pos)
1627   # >= 1: skip first part (this can be zero)
1628   while ($i > 0) { last if $s->[$i] != 0; $i--; }
1629   $i++; splice @$s,$i if ($i < $cnt); # $i cant be 0
1630   return $s;
1631   }
1632
1633 sub _from_hex
1634   {
1635   # convert a (ref to) big hex string to BigInt, return undef for error
1636   my $hs = shift;
1637
1638   my $x = Math::BigInt->bzero();
1639   return $x->bnan() if $$hs !~ /^[\-\+]?0x[0-9A-Fa-f]+$/;
1640
1641   my $mul = Math::BigInt->bzero(); $mul++;
1642   my $x65536 = Math::BigInt->new(65536);
1643
1644   my $len = CORE::length($$hs)-2; my $sign = '+';
1645   if ($$hs =~ /^\-/)
1646     {
1647     $sign = '-'; $len--;
1648     }
1649   $len = int($len/4);                           # 4-digit parts, w/o '0x'
1650   my $val; my $i = -4;
1651   while ($len >= 0)
1652     {
1653     $val = substr($$hs,$i,4);
1654     $val =~ s/^[\-\+]?0x// if $len == 0;        # for last part only because
1655     $val = hex($val);                           # hex does not like wrong chars
1656     # print "$val ",substr($$hs,$i,4),"\n";
1657     $i -= 4; $len --;
1658     $x += $mul * $val if $val != 0;
1659     $mul *= $x65536 if $len >= 0;               # skip last mul
1660     }
1661   $x->{sign} = $sign if !$x->is_zero();
1662   return $x;
1663   }
1664
1665 sub _from_bin
1666   {
1667   # convert a (ref to) big binary string to BigInt, return undef for error
1668   my $bs = shift;
1669
1670   my $x = Math::BigInt->bzero();
1671   return $x->bnan() if $$bs !~ /^[\-\+]?0b[01]+$/;
1672
1673   my $mul = Math::BigInt->bzero(); $mul++;
1674   my $x256 = Math::BigInt->new(256);
1675
1676   my $len = CORE::length($$bs)-2; my $sign = '+';
1677   if ($$bs =~ /^\-/)
1678     {
1679     $sign = '-'; $len--;
1680     }
1681   $len = int($len/8);                           # 8-digit parts, w/o '0b'
1682   my $val; my $i = -8;
1683   while ($len >= 0)
1684     {
1685     $val = substr($$bs,$i,8);
1686     $val =~ s/^[\-\+]?0b// if $len == 0;        # for last part only
1687     #$val = oct('0b'.$val);     # does not work on Perl prior 5.6.0
1688     $val = ('0' x (8-CORE::length($val))).$val if CORE::length($val) < 8;
1689     $val = ord(pack('B8',$val));
1690     # print "$val ",substr($$bs,$i,16),"\n";
1691     $i -= 8; $len --;
1692     $x += $mul * $val if $val != 0;
1693     $mul *= $x256 if $len >= 0;                 # skip last mul
1694     }
1695   $x->{sign} = $sign if !$x->is_zero();
1696   return $x;
1697   }
1698
1699 sub _split
1700   {
1701   # (ref to num_str) return num_str
1702   # internal, take apart a string and return the pieces
1703   my $x = shift;
1704
1705   # pre-parse input
1706   $$x =~ s/^\s+//g;                     # strip white space at front
1707   $$x =~ s/\s+$//g;                     # strip white space at end
1708   #$$x =~ s/\s+//g;                     # strip white space (no longer)
1709   return if $$x eq "";
1710
1711   return _from_hex($x) if $$x =~ /^[\-\+]?0x/;  # hex string
1712   return _from_bin($x) if $$x =~ /^[\-\+]?0b/;  # binary string
1713
1714   return if $$x !~ /^[\-\+]?\.?[0-9]/;
1715
1716   $$x =~ s/(\d)_(\d)/$1$2/g;            # strip underscores between digits
1717   $$x =~ s/(\d)_(\d)/$1$2/g;            # do twice for 1_2_3
1718
1719   # some possible inputs:
1720   # 2.1234 # 0.12        # 1          # 1E1 # 2.134E1 # 434E-10 # 1.02009E-2
1721   # .2     # 1_2_3.4_5_6 # 1.4E1_2_3  # 1e3 # +.2
1722
1723   #print "input: '$$x' ";
1724   my ($m,$e) = split /[Ee]/,$$x;
1725   $e = '0' if !defined $e || $e eq "";
1726   # print "m '$m' e '$e'\n";
1727   # sign,value for exponent,mantint,mantfrac
1728   my ($es,$ev,$mis,$miv,$mfv);
1729   # valid exponent?
1730   if ($e =~ /^([+-]?)0*(\d+)$/) # strip leading zeros
1731     {
1732     $es = $1; $ev = $2;
1733     #print "'$m' '$e' e: $es $ev ";
1734     # valid mantissa?
1735     return if $m eq '.' || $m eq '';
1736     my ($mi,$mf) = split /\./,$m;
1737     $mi = '0' if !defined $mi;
1738     $mi .= '0' if $mi =~ /^[\-\+]?$/;
1739     $mf = '0' if !defined $mf || $mf eq '';
1740     if ($mi =~ /^([+-]?)0*(\d+)$/) # strip leading zeros
1741       {
1742       $mis = $1||'+'; $miv = $2;
1743       #print "$mis $miv";
1744       # valid, existing fraction part of mantissa?
1745       return unless ($mf =~ /^(\d*?)0*$/);      # strip trailing zeros
1746       $mfv = $1;
1747       #print " split: $mis $miv . $mfv E $es $ev\n";
1748       return (\$mis,\$miv,\$mfv,\$es,\$ev);
1749       }
1750     }
1751   return; # NaN, not a number
1752   }
1753
1754 sub _digits
1755   {
1756   # computer number of digits in bigint, minus the sign
1757   # int() because add/sub leaves sometimes strings (like '00005') instead of
1758   # int ('5') in this place, causing length to fail
1759   my $cx = shift;
1760
1761   #print "len: ",(@$cx-1)*5+CORE::length(int($cx->[-1])),"\n";
1762   return (@$cx-1)*5+CORE::length(int($cx->[-1]));
1763   }
1764
1765 sub as_number
1766   {
1767   # an object might be asked to return itself as bigint on certain overloaded
1768   # operations, this does exactly this, so that sub classes can simple inherit
1769   # it or override with their own integer conversion routine
1770   my $self = shift;
1771
1772   return $self->copy();
1773   }
1774
1775 ##############################################################################
1776 # internal calculation routines
1777
1778 sub acmp
1779   {
1780   # internal absolute post-normalized compare (ignore signs)
1781   # ref to array, ref to array, return <0, 0, >0
1782   # arrays must have at least on entry, this is not checked for
1783
1784   my ($cx, $cy) = @_;
1785
1786   #print "$cx $cy\n";
1787   my ($i,$a,$x,$y,$k);
1788   # calculate length based on digits, not parts
1789   $x = _digits($cx); $y = _digits($cy);
1790   # print "length: ",($x-$y),"\n";
1791   return $x-$y if ($x - $y);              # if different in length
1792   #print "full compare\n";
1793   $i = 0; $a = 0;
1794   # first way takes 5.49 sec instead of 4.87, but has the early out advantage
1795   # so grep is slightly faster, but more unflexible. hm. $_ instead if $k
1796   # yields 5.6 instead of 5.5 sec huh?
1797   # manual way (abort if unequal, good for early ne)
1798   my $j = scalar @$cx - 1;
1799   while ($j >= 0)
1800    {
1801    # print "$cx->[$j] $cy->[$j] $a",$cx->[$j]-$cy->[$j],"\n";
1802    last if ($a = $cx->[$j] - $cy->[$j]); $j--;
1803    }
1804   return $a;
1805   # while it early aborts, it is even slower than the manual variant
1806   #grep { return $a if ($a = $_ - $cy->[$i++]); } @$cx;
1807   # grep way, go trough all (bad for early ne)
1808   #grep { $a = $_ - $cy->[$i++]; } @$cx;
1809   #return $a;
1810   }
1811
1812 sub cmp
1813   {
1814   # post-normalized compare for internal use (honors signs)
1815   # ref to array, ref to array, return < 0, 0, >0
1816   my ($cx,$cy,$sx,$sy) = @_;
1817
1818   #return 0 if (is0($cx,$sx) && is0($cy,$sy));
1819
1820   if ($sx eq '+')
1821     {
1822     return 1 if $sy eq '-'; # 0 check handled above
1823     return acmp($cx,$cy);
1824     }
1825   else
1826     {
1827     # $sx eq '-'
1828     return -1 if ($sy eq '+');
1829     return acmp($cy,$cx);
1830     }
1831   return 0; # equal
1832   }
1833
1834 sub add
1835   {
1836   # (ref to int_num_array, ref to int_num_array)
1837   # routine to add two base 1e5 numbers
1838   # stolen from Knuth Vol 2 Algorithm A pg 231
1839   # there are separate routines to add and sub as per Kunth pg 233
1840   # This routine clobbers up array x, but not y.
1841
1842   my ($x,$y) = @_;
1843
1844   # for each in Y, add Y to X and carry. If after that, something is left in
1845   # X, foreach in X add carry to X and then return X, carry
1846   # Trades one "$j++" for having to shift arrays, $j could be made integer
1847   # but this would impose a limit to number-length to 2**32.
1848   my $i; my $car = 0; my $j = 0;
1849   for $i (@$y)
1850     {
1851     $x->[$j] -= $BASE
1852       if $car = (($x->[$j] += $i + $car) >= $BASE) ? 1 : 0;
1853     $j++;
1854     }
1855   while ($car != 0)
1856     {
1857     $x->[$j] -= $BASE if $car = (($x->[$j] += $car) >= $BASE) ? 1 : 0; $j++;
1858     }
1859   }
1860
1861 sub sub
1862   {
1863   # (ref to int_num_array, ref to int_num_array)
1864   # subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
1865   # subtract Y from X (X is always greater/equal!) by modifiyng x in place
1866   my ($sx,$sy,$s) = @_;
1867
1868   my $car = 0; my $i; my $j = 0;
1869   if (!$s)
1870     {
1871     #print "case 2\n";
1872     for $i (@$sx)
1873       {
1874       last unless defined $sy->[$j] || $car;
1875       #print "x: $i y: $sy->[$j] c: $car\n";
1876       $i += $BASE if $car = (($i -= ($sy->[$j] || 0) + $car) < 0); $j++;
1877       #print "x: $i y: $sy->[$j-1] c: $car\n";
1878       }
1879     # might leave leading zeros, so fix that
1880     _strip_zeros($sx);
1881     return $sx;
1882     }
1883   else
1884     {
1885     #print "case 1 (swap)\n";
1886     for $i (@$sx)
1887       {
1888       last unless defined $sy->[$j] || $car;
1889       #print "$sy->[$j] $i $car => $sx->[$j]\n";
1890       $sy->[$j] += $BASE
1891        if $car = (($sy->[$j] = $i-($sy->[$j]||0) - $car) < 0);
1892       #print "$sy->[$j] $i $car => $sy->[$j]\n";
1893       $j++;
1894       }
1895     # might leave leading zeros, so fix that
1896     _strip_zeros($sy);
1897     return $sy;
1898     }
1899   }
1900
1901 sub mul
1902   {
1903   # (BINT, BINT) return nothing
1904   # multiply two numbers in internal representation
1905   # modifies first arg, second needs not be different from first
1906   my ($x,$y) = @_;
1907
1908   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1909   my @prod = (); my ($prod,$car,$cty,$xi,$yi);
1910   my $xv = $x->{value};
1911   my $yv = $y->{value};
1912   # since multiplying $x with $x fails, make copy in this case
1913   $yv = [@$xv] if "$xv" eq "$yv";
1914   for $xi (@$xv)
1915     {
1916     $car = 0; $cty = 0;
1917     for $yi (@$yv)
1918       {
1919       $prod = $xi * $yi + ($prod[$cty] || 0) + $car;
1920       $prod[$cty++] =
1921        $prod - ($car = int($prod * 1e-5)) * $BASE;      # see USE_MUL
1922       }
1923     $prod[$cty] += $car if $car; # need really to check for 0?
1924     $xi = shift @prod;
1925     }
1926   push @$xv, @prod;
1927   _strip_zeros($x->{value});
1928   # normalize (handled last to save check for $y->is_zero()
1929   $x->{sign} = '+' if @$xv == 1 && $xv->[0] == 0; # not is_zero due to '-'
1930   }
1931
1932 sub div
1933   {
1934   # ref to array, ref to array, modify first array and return reminder if
1935   # in list context
1936   # does no longer handle sign
1937   my ($x,$yorg) = @_;
1938   my ($car,$bar,$prd,$dd,$xi,$yi,@q,$v2,$v1);
1939
1940   my (@d,$tmp,$q,$u2,$u1,$u0);
1941
1942   $car = $bar = $prd = 0;
1943
1944   my $y = [ @$yorg ];
1945   if (($dd = int($BASE/($y->[-1]+1))) != 1)
1946     {
1947     for $xi (@$x)
1948       {
1949       $xi = $xi * $dd + $car;
1950       $xi -= ($car = int($xi * $RBASE)) * $BASE;        # see USE_MUL
1951       }
1952     push(@$x, $car); $car = 0;
1953     for $yi (@$y)
1954       {
1955       $yi = $yi * $dd + $car;
1956       $yi -= ($car = int($yi * $RBASE)) * $BASE;        # see USE_MUL
1957       }
1958     }
1959   else
1960     {
1961     push(@$x, 0);
1962     }
1963   @q = (); ($v2,$v1) = @$y[-2,-1];
1964   $v2 = 0 unless $v2;
1965   while ($#$x > $#$y)
1966     {
1967     ($u2,$u1,$u0) = @$x[-3..-1];
1968     $u2 = 0 unless $u2;
1969     print "oups v1 is 0, u0: $u0 $y->[-2] $y->[-1] l ",scalar @$y,"\n"
1970      if $v1 == 0;
1971     $q = (($u0 == $v1) ? 99999 : int(($u0*$BASE+$u1)/$v1));
1972     --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*$BASE+$u2);
1973     if ($q)
1974       {
1975       ($car, $bar) = (0,0);
1976       for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
1977         {
1978         $prd = $q * $y->[$yi] + $car;
1979         $prd -= ($car = int($prd * $RBASE)) * $BASE;    # see USE_MUL
1980         $x->[$xi] += 1e5 if ($bar = (($x->[$xi] -= $prd + $bar) < 0));
1981         }
1982       if ($x->[-1] < $car + $bar)
1983         {
1984         $car = 0; --$q;
1985         for ($yi = 0, $xi = $#$x-$#$y-1; $yi <= $#$y; ++$yi,++$xi)
1986           {
1987           $x->[$xi] -= 1e5
1988            if ($car = (($x->[$xi] += $y->[$yi] + $car) > $BASE));
1989           }
1990         }
1991       }
1992       pop(@$x); unshift(@q, $q);
1993     }
1994   if (wantarray)
1995     {
1996     @d = ();
1997     if ($dd != 1)
1998       {
1999       $car = 0;
2000       for $xi (reverse @$x)
2001         {
2002         $prd = $car * $BASE + $xi;
2003         $car = $prd - ($tmp = int($prd / $dd)) * $dd; # see USE_MUL
2004         unshift(@d, $tmp);
2005         }
2006       }
2007     else
2008       {
2009       @d = @$x;
2010       }
2011     @$x = @q;
2012     _strip_zeros($x);
2013     _strip_zeros(\@d);
2014     return ($x,\@d);
2015     }
2016   @$x = @q;
2017   _strip_zeros($x);
2018   return $x;
2019   }
2020
2021 sub _lcm
2022   {
2023   # (BINT or num_str, BINT or num_str) return BINT
2024   # does modify first argument
2025   # LCM
2026
2027   my $x = shift; my $ty = shift;
2028   return $x->bnan() if ($x->{sign} eq $nan) || ($ty->{sign} eq $nan);
2029   return $x * $ty / bgcd($x,$ty);
2030   }
2031
2032 sub _gcd_old
2033   {
2034   # (BINT or num_str, BINT or num_str) return BINT
2035   # does modify first arg
2036   # GCD -- Euclids algorithm E, Knuth Vol 2 pg 296
2037   trace(@_);
2038
2039   my $x = shift; my $ty = $class->new(shift); # preserve y
2040   return $x->bnan() if ($x->{sign} eq $nan) || ($ty->{sign} eq $nan);
2041
2042   while (!$ty->is_zero())
2043     {
2044     ($x, $ty) = ($ty,bmod($x,$ty));
2045     }
2046   $x;
2047   }
2048
2049 sub _gcd
2050   {
2051   # (BINT or num_str, BINT or num_str) return BINT
2052   # does not modify arguments
2053   # GCD -- Euclids algorithm, variant L (Lehmer), Knuth Vol 3 pg 347 ff
2054   # unfortunately, it is slower and also seems buggy (the A=0, B=1, C=1, D=0
2055   # case..)
2056   trace(@_);
2057
2058   my $u = $class->new(shift); my $v = $class->new(shift);       # preserve u,v
2059   return $u->bnan() if ($u->{sign} eq $nan) || ($v->{sign} eq $nan);
2060
2061   $u->babs(); $v->babs();               # Euclid is valid for |u| and |v|
2062
2063   my ($U,$V,$A,$B,$C,$D,$T,$Q);         # single precision variables
2064   my ($t);                              # multiprecision variables
2065
2066   while ($v > $BASE)
2067     {
2068     #sleep 1;
2069     ($u,$v) = ($v,$u) if ($u < $v);             # make sure that u >= v
2070     #print "gcd: $u $v\n";
2071     # step L1, initialize
2072     $A = 1; $B = 0; $C = 0; $D = 1;
2073     $U = $u->{value}->[-1];                     # leading digits of u
2074     $V = $v->{value}->[-1];                     # leading digits of v
2075
2076     # step L2, test quotient
2077     if (($V + $C != 0) && ($V + $D != 0))       # div by zero => go to L4
2078       {
2079       $Q = int(($U + $A)/($V + $C));            # quotient
2080       #print "L1 A=$A B=$B C=$C D=$D U=$U V=$V Q=$Q\n";
2081       # do L3? (div by zero => go to L4)
2082       while ($Q == int(($U + $B)/($V + $D)))
2083         {
2084         # step L3, emulate Euclid
2085         #print "L3a A=$A B=$B C=$C D=$D U=$U V=$V Q=$Q\n";
2086         $T = $A - $Q*$C; $A = $C; $C = $T;
2087         $T = $B - $Q*$D; $B = $D; $D = $T;
2088         $T = $U - $Q*$V; $U = $V; $V = $T;
2089         last if ($V + $D == 0) || ($V + $C == 0);       # div by zero => L4
2090         $Q = int(($U + $A)/($V + $C));  # quotient for next test
2091         #print "L3b A=$A B=$B C=$C D=$D U=$U V=$V Q=$Q\n";
2092         }
2093       }
2094     # step L4, multiprecision step
2095     # was if ($B == 0)
2096     # in case A = 0, B = 1, C = 0 and D = 1, this case would simple swap u & v
2097     # and loop endless. Not sure why this happens, Knuth does not make a
2098     # remark about this special case. bug?
2099     if (($B == 0) || (($A == 0) && ($C == 1) && ($D == 0)))
2100       {
2101       #print "L4b1: u=$u v=$v\n";
2102       ($u,$v) = ($v,bmod($u,$v));
2103       #$t = $u % $v; $u = $v->copy(); $v = $t;
2104       #print "L4b12: u=$u v=$v\n";
2105       }
2106     else
2107       {
2108       #print "L4b: $u $v $A $B $C $D\n";
2109       $t = $A*$u + $B*$v; $v *= $D; $v += $C*$u; $u = $t;
2110       #print "L4b2: $u $v\n";
2111       }
2112     } # back to L1
2113
2114   return _gcd_old($u,$v) if $v != 1;    # v too small
2115   return $v;                            # 1
2116   }
2117
2118 ###############################################################################
2119 # this method return 0 if the object can be modified, or 1 for not
2120 # We use a fast use constant statement here, to avoid costly calls. Subclasses
2121 # may override it with special code (f.i. Math::BigInt::Constant does so)
2122
2123 use constant modify => 0;
2124
2125 #sub modify
2126 #  {
2127 #  my $self = shift;
2128 #  my $method = shift;
2129 #  print "original $self modify by $method\n";
2130 #  return 0; # $self;
2131 #  }
2132
2133 1;
2134 __END__
2135
2136 =head1 NAME
2137
2138 Math::BigInt - Arbitrary size integer math package
2139
2140 =head1 SYNOPSIS
2141
2142   use Math::BigInt;
2143
2144   # Number creation
2145   $x = Math::BigInt->new($str); # defaults to 0
2146   $nan  = Math::BigInt->bnan(); # create a NotANumber
2147   $zero = Math::BigInt->bzero();# create a "+0"
2148
2149   # Testing
2150   $x->is_zero();                # return whether arg is zero or not
2151   $x->is_nan();                 # return whether arg is NaN or not
2152   $x->is_one();                 # return true if arg is +1
2153   $x->is_one('-');              # return true if arg is -1
2154   $x->is_odd();                 # return true if odd, false for even
2155   $x->is_even();                # return true if even, false for odd
2156   $x->bcmp($y);                 # compare numbers (undef,<0,=0,>0)
2157   $x->bacmp($y);                # compare absolutely (undef,<0,=0,>0)
2158   $x->sign();                   # return the sign, either +,- or NaN
2159   $x->digit($n);                # return the nth digit, counting from right
2160   $x->digit(-$n);               # return the nth digit, counting from left
2161
2162   # The following all modify their first argument:
2163
2164   # set
2165   $x->bzero();                  # set $x to 0
2166   $x->bnan();                   # set $x to NaN
2167
2168   $x->bneg();                   # negation
2169   $x->babs();                   # absolute value
2170   $x->bnorm();                  # normalize (no-op)
2171   $x->bnot();                   # two's complement (bit wise not)
2172   $x->binc();                   # increment x by 1
2173   $x->bdec();                   # decrement x by 1
2174
2175   $x->badd($y);                 # addition (add $y to $x)
2176   $x->bsub($y);                 # subtraction (subtract $y from $x)
2177   $x->bmul($y);                 # multiplication (multiply $x by $y)
2178   $x->bdiv($y);                 # divide, set $x to quotient
2179                                 # return (quo,rem) or quo if scalar
2180
2181   $x->bmod($y);                 # modulus (x % y)
2182   $x->bpow($y);                 # power of arguments (x ** y)
2183   $x->blsft($y);                # left shift
2184   $x->brsft($y);                # right shift
2185   $x->blsft($y,$n);             # left shift, by base $n (like 10)
2186   $x->brsft($y,$n);             # right shift, by base $n (like 10)
2187
2188   $x->band($y);                 # bitwise and
2189   $x->bior($y);                 # bitwise inclusive or
2190   $x->bxor($y);                 # bitwise exclusive or
2191   $x->bnot();                   # bitwise not (two's complement)
2192
2193   $x->bsqrt();                  # calculate square-root
2194
2195   $x->round($A,$P,$round_mode); # round to accuracy or precision using mode $r
2196   $x->bround($N);               # accuracy: preserve $N digits
2197   $x->bfround($N);              # round to $Nth digit, no-op for BigInts
2198
2199   # The following do not modify their arguments in BigInt, but do in BigFloat:
2200   $x->bfloor();                 # return integer less or equal than $x
2201   $x->bceil();                  # return integer greater or equal than $x
2202
2203   # The following do not modify their arguments:
2204
2205   bgcd(@values);                # greatest common divisor
2206   blcm(@values);                # lowest common multiplicator
2207
2208   $x->bstr();                   # normalized string
2209   $x->bsstr();                  # normalized string in scientific notation
2210   $x->length();                 # return number of digits in number
2211   ($x,$f) = $x->length();       # length of number and length of fraction part
2212
2213   $x->exponent();               # return exponent as BigInt
2214   $x->mantissa();               # return mantissa as BigInt
2215   $x->parts();                  # return (mantissa,exponent) as BigInt
2216
2217 =head1 DESCRIPTION
2218
2219 All operators (inlcuding basic math operations) are overloaded if you
2220 declare your big integers as
2221
2222   $i = new Math::BigInt '123_456_789_123_456_789';
2223
2224 Operations with overloaded operators preserve the arguments which is
2225 exactly what you expect.
2226
2227 =over 2
2228
2229 =item Canonical notation
2230
2231 Big integer values are strings of the form C</^[+-]\d+$/> with leading
2232 zeros suppressed.
2233
2234    '-0'                            canonical value '-0', normalized '0'
2235    '   -123_123_123'               canonical value '-123123123'
2236    '1_23_456_7890'                 canonical value '1234567890'
2237
2238 =item Input
2239
2240 Input values to these routines may be either Math::BigInt objects or
2241 strings of the form C</^\s*[+-]?[\d]+\.?[\d]*E?[+-]?[\d]*$/>.
2242
2243 You can include one underscore between any two digits.
2244
2245 This means integer values like 1.01E2 or even 1000E-2 are also accepted.
2246 Non integer values result in NaN.
2247
2248 Math::BigInt::new() defaults to 0, while Math::BigInt::new('') results
2249 in 'NaN'.
2250
2251 bnorm() on a BigInt object is now effectively a no-op, since the numbers
2252 are always stored in normalized form. On a string, it creates a BigInt
2253 object.
2254
2255 =item Output
2256
2257 Output values are BigInt objects (normalized), except for bstr(), which
2258 returns a string in normalized form.
2259 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2260 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2261 return either undef, <0, 0 or >0 and are suited for sort.
2262
2263 =back
2264
2265 =head2 Rounding
2266
2267 Only C<bround()> is defined for BigInts, for further details of rounding see
2268 L<Math::BigFloat>.
2269
2270 =over 2
2271
2272 =item bfround ( +$scale ) rounds to the $scale'th place left from the '.'
2273
2274 =item bround  ( +$scale ) preserves accuracy to $scale sighnificant digits counted from the left and paddes the number with zeros
2275
2276 =item bround  ( -$scale ) preserves accuracy to $scale significant digits counted from the right and paddes the number with zeros.
2277
2278 =back
2279
2280 C<bfround()> does nothing in case of negative C<$scale>. Both C<bround()> and
2281 C<bfround()> are a no-ops for a scale of 0.
2282
2283 All rounding functions take as a second parameter a rounding mode from one of
2284 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2285
2286 The default is 'even'. By using C<< Math::BigInt->round_mode($rnd_mode); >>
2287 you can get and set the default round mode for subsequent rounding.
2288
2289 The second parameter to the round functions than overrides the default
2290 temporarily.
2291
2292 =head2 Internals
2293
2294 Actual math is done in an internal format consisting of an array of
2295 elements of base 100000 digits with the least significant digit first.
2296 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2297 represent the result when input arguments are not numbers, as well as
2298 the result of dividing by zero.
2299
2300 You sould neither care nor depend on the internal represantation, it might
2301 change without notice. Use only method calls like C<< $x->sign(); >> instead
2302 relying on the internal hash keys like in C<< $x->{sign}; >>.
2303
2304 =head2 mantissa(), exponent() and parts()
2305
2306 C<mantissa()> and C<exponent()> return the said parts of the BigInt such
2307 that:
2308
2309         $m = $x->mantissa();
2310         $e = $x->exponent();
2311         $y = $m * ( 10 ** $e );
2312         print "ok\n" if $x == $y;
2313
2314 C<($m,$e) = $x->parts()> is just a shortcut that gives you both of them in one
2315 go. Both the returned mantissa and exponent do have a sign.
2316
2317 Currently, for BigInts C<$e> will be always 0, except for NaN where it will be
2318 NaN and for $x == 0, then it will be 1 (to be compatible with Math::BigFlaot's
2319 internal representation of a zero as C<0E1>).
2320
2321 C<$m> will always be a copy of the original number. The relation between $e
2322 and $m might change in the future, but will be always equivalent in a
2323 numerical sense, e.g. $m might get minized.
2324
2325 =head1 EXAMPLES
2326
2327   use Math::BigInt qw(bstr bint);
2328   $x = bstr("1234")                     # string "1234"
2329   $x = "$x";                            # same as bstr()
2330   $x = bneg("1234")                     # Bigint "-1234"
2331   $x = Math::BigInt->bneg("1234");      # Bigint "-1234"
2332   $x = Math::BigInt->babs("-12345");    # Bigint "12345"
2333   $x = Math::BigInt->bnorm("-0 00");    # BigInt "0"
2334   $x = bint(1) + bint(2);               # BigInt "3"
2335   $x = bint(1) + "2";                   # ditto (auto-BigIntify of "2")
2336   $x = bint(1);                         # BigInt "1"
2337   $x = $x + 5 / 2;                      # BigInt "3"
2338   $x = $x ** 3;                         # BigInt "27"
2339   $x *= 2;                              # BigInt "54"
2340   $x = new Math::BigInt;                # BigInt "0"
2341   $x--;                                 # BigInt "-1"
2342   $x = Math::BigInt->badd(4,5)          # BigInt "9"
2343   $x = Math::BigInt::badd(4,5)          # BigInt "9"
2344   print $x->bsstr();                    # 9e+0
2345
2346 =head1 Autocreating constants
2347
2348 After C<use Math::BigInt ':constant'> all the B<integer> decimal constants
2349 in the given scope are converted to C<Math::BigInt>. This conversion
2350 happens at compile time.
2351
2352 In particular
2353
2354   perl -MMath::BigInt=:constant -e 'print 2**100,"\n"'
2355
2356 prints the integer value of C<2**100>.  Note that without conversion of
2357 constants the expression 2**100 will be calculated as floating point
2358 number.
2359
2360 Please note that strings and floating point constants are not affected,
2361 so that
2362
2363         use Math::BigInt qw/:constant/;
2364
2365         $x = 1234567890123456789012345678901234567890
2366                 + 123456789123456789;
2367         $x = '1234567890123456789012345678901234567890'
2368                 + '123456789123456789';
2369
2370 do both not work. You need a explicit Math::BigInt->new() around one of them.
2371
2372 =head1 PERFORMANCE
2373
2374 Using the form $x += $y; etc over $x = $x + $y is faster, since a copy of $x
2375 must be made in the second case. For long numbers, the copy can eat up to 20%
2376 of the work (in case of addition/subtraction, less for
2377 multiplication/division). If $y is very small compared to $x, the form
2378 $x += $y is MUCH faster than $x = $x + $y since making the copy of $x takes
2379 more time then the actual addition.
2380
2381 With a technic called copy-on-write the cost of copying with overload could
2382 be minimized or even completely avoided. This is currently not implemented.
2383
2384 The new version of this module is slower on new(), bstr() and numify(). Some
2385 operations may be slower for small numbers, but are significantly faster for
2386 big numbers. Other operations are now constant (O(1), like bneg(), babs()
2387 etc), instead of O(N) and thus nearly always take much less time.
2388
2389 For more benchmark results see http://bloodgate.com/perl/benchmarks.html
2390
2391 =head1 BUGS
2392
2393 =over 2
2394
2395 =item :constant and eval()
2396
2397 Under Perl prior to 5.6.0 having an C<use Math::BigInt ':constant';> and
2398 C<eval()> in your code will crash with "Out of memory". This is probably an
2399 overload/exporter bug. You can workaround by not having C<eval()>
2400 and ':constant' at the same time or upgrade your Perl.
2401
2402 =back
2403
2404 =head1 CAVEATS
2405
2406 Some things might not work as you expect them. Below is documented what is
2407 known to be troublesome:
2408
2409 =over 1
2410
2411 =item stringify, bstr(), bsstr() and 'cmp'
2412
2413 Both stringify and bstr() now drop the leading '+'. The old code would return
2414 '+3', the new returns '3'. This is to be consistent with Perl and to make
2415 cmp (especially with overloading) to work as you expect. It also solves
2416 problems with Test.pm, it's ok() uses 'eq' internally.
2417
2418 Mark said, when asked about to drop the '+' altogether, or make only cmp work:
2419
2420         I agree (with the first alternative), don't add the '+' on positive
2421         numbers.  It's not as important anymore with the new internal
2422         form for numbers.  It made doing things like abs and neg easier,
2423         but those have to be done differently now anyway.
2424
2425 So, the following examples will now work all as expected:
2426
2427         use Test;
2428         BEGIN { plan tests => 1 }
2429         use Math::BigInt;
2430
2431         my $x = new Math::BigInt 3*3;
2432         my $y = new Math::BigInt 3*3;
2433
2434         ok ($x,3*3);
2435         print "$x eq 9" if $x eq $y;
2436         print "$x eq 9" if $x eq '9';
2437         print "$x eq 9" if $x eq 3*3;
2438
2439 Additionally, the following still works:
2440
2441         print "$x == 9" if $x == $y;
2442         print "$x == 9" if $x == 9;
2443         print "$x == 9" if $x == 3*3;
2444
2445 There is now a C<bsstr()> method to get the string in scientific notation aka
2446 C<1e+2> instead of C<100>. Be advised that overloaded 'eq' always uses bstr()
2447 for comparisation, but Perl will represent some numbers as 100 and others
2448 as 1e+308. If in doubt, convert both arguments to Math::BigInt before doing eq:
2449
2450         use Test;
2451         BEGIN { plan tests => 3 }
2452         use Math::BigInt;
2453
2454         $x = Math::BigInt->new('1e56'); $y = 1e56;
2455         ok ($x,$y);                     # will fail
2456         ok ($x->bsstr(),$y);            # okay
2457         $y = Math::BigInt->new($y);
2458         ok ($x,$y);                     # okay
2459
2460 =item int()
2461
2462 C<int()> will return (at least for Perl v5.7.1 and up) another BigInt, not a
2463 Perl scalar:
2464
2465         $x = Math::BigInt->new(123);
2466         $y = int($x);                           # BigInt 123
2467         $x = Math::BigFloat->new(123.45);
2468         $y = int($x);                           # BigInt 123
2469
2470 In all Perl versions you can use C<as_number()> for the same effect:
2471
2472         $x = Math::BigFloat->new(123.45);
2473         $y = $x->as_number();                   # BigInt 123
2474
2475 This also works for other subclasses, like Math::String.
2476
2477 =item bdiv
2478
2479 The following will probably not do what you expect:
2480
2481         print $c->bdiv(10000),"\n";
2482
2483 It prints both quotient and reminder since print calls C<bdiv()> in list
2484 context. Also, C<bdiv()> will modify $c, so be carefull. You probably want
2485 to use
2486
2487         print $c / 10000,"\n";
2488         print scalar $c->bdiv(10000),"\n";  # or if you want to modify $c
2489
2490 instead.
2491
2492 The quotient is always the greatest integer less than or equal to the
2493 real-valued quotient of the two operands, and the remainder (when it is
2494 nonzero) always has the same sign as the second operand; so, for
2495 example,
2496
2497         1 / 4   => ( 0, 1)
2498         1 / -4  => (-1,-3)
2499         -3 / 4  => (-1, 1)
2500         -3 / -4 => ( 0,-3)
2501
2502 As a consequence, the behavior of the operator % agrees with the
2503 behavior of Perl's built-in % operator (as documented in the perlop
2504 manpage), and the equation
2505
2506         $x == ($x / $y) * $y + ($x % $y)
2507
2508 holds true for any $x and $y, which justifies calling the two return
2509 values of bdiv() the quotient and remainder.
2510
2511 Perl's 'use integer;' changes the behaviour of % and / for scalars, but will
2512 not change BigInt's way to do things. This is because under 'use integer' Perl
2513 will do what the underlying C thinks is right and this is different for each
2514 system. If you need BigInt's behaving exactly like Perl's 'use integer', bug
2515 the author to implement it ;)
2516
2517 =item Modifying and =
2518
2519 Beware of:
2520
2521         $x = Math::BigFloat->new(5);
2522         $y = $x;
2523
2524 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2525 a second reference to the B<same> object and stores it in $y. Thus anything
2526 that modifies $x will modify $y, and vice versa.
2527
2528         $x->bmul(2);
2529         print "$x, $y\n";       # prints '10, 10'
2530
2531 If you want a true copy of $x, use:
2532
2533         $y = $x->copy();
2534
2535 See also the documentation in L<overload> regarding C<=>.
2536
2537 =item bpow
2538
2539 C<bpow()> (and the rounding functions) now modifies the first argument and
2540 return it, unlike the old code which left it alone and only returned the
2541 result. This is to be consistent with C<badd()> etc. The first three will
2542 modify $x, the last one won't:
2543
2544         print bpow($x,$i),"\n";         # modify $x
2545         print $x->bpow($i),"\n";        # ditto
2546         print $x **= $i,"\n";           # the same
2547         print $x ** $i,"\n";            # leave $x alone
2548
2549 The form C<$x **= $y> is faster than C<$x = $x ** $y;>, though.
2550
2551 =item Overloading -$x
2552
2553 The following:
2554
2555         $x = -$x;
2556
2557 is slower than
2558
2559         $x->bneg();
2560
2561 since overload calls C<sub($x,0,1);> instead of C<neg($x)>. The first variant
2562 needs to preserve $x since it does not know that it later will get overwritten.
2563 This makes a copy of $x and takes O(N). But $x->bneg() is O(1).
2564
2565 With Copy-On-Write, this issue will be gone. Stay tuned...
2566
2567 =item Mixing different object types
2568
2569 In Perl you will get a floating point value if you do one of the following:
2570
2571         $float = 5.0 + 2;
2572         $float = 2 + 5.0;
2573         $float = 5 / 2;
2574
2575 With overloaded math, only the first two variants will result in a BigFloat:
2576
2577         use Math::BigInt;
2578         use Math::BigFloat;
2579
2580         $mbf = Math::BigFloat->new(5);
2581         $mbi2 = Math::BigInteger->new(5);
2582         $mbi = Math::BigInteger->new(2);
2583
2584                                         # what actually gets called:
2585         $float = $mbf + $mbi;           # $mbf->badd()
2586         $float = $mbf / $mbi;           # $mbf->bdiv()
2587         $integer = $mbi + $mbf;         # $mbi->badd()
2588         $integer = $mbi2 / $mbi;        # $mbi2->bdiv()
2589         $integer = $mbi2 / $mbf;        # $mbi2->bdiv()
2590
2591 This is because math with overloaded operators follows the first (dominating)
2592 operand, this one's operation is called and returns thus the result. So,
2593 Math::BigInt::bdiv() will always return a Math::BigInt, regardless whether
2594 the result should be a Math::BigFloat or the second operant is one.
2595
2596 To get a Math::BigFloat you either need to call the operation manually,
2597 make sure the operands are already of the proper type or casted to that type
2598 via Math::BigFloat->new():
2599
2600         $float = Math::BigFloat->new($mbi2) / $mbi;     # = 2.5
2601
2602 Beware of simple "casting" the entire expression, this would only convert
2603 the already computed result:
2604
2605         $float = Math::BigFloat->new($mbi2 / $mbi);     # = 2.0 thus wrong!
2606
2607 Beware of the order of more complicated expressions like:
2608
2609         $integer = ($mbi2 + $mbi) / $mbf;               # int / float => int
2610         $integer = $mbi2 / Math::BigFloat->new($mbi);   # ditto
2611
2612 If in doubt, break the expression into simpler terms, or cast all operands
2613 to the desired resulting type.
2614
2615 Scalar values are a bit different, since:
2616
2617         $float = 2 + $mbf;
2618         $float = $mbf + 2;
2619
2620 will both result in the proper type due to the way the overloaded math works.
2621
2622 This section also applies to other overloaded math packages, like Math::String.
2623
2624 =item bsqrt()
2625
2626 C<bsqrt()> works only good if the result is an big integer, e.g. the square
2627 root of 144 is 12, but from 12 the square root is 3, regardless of rounding
2628 mode.
2629
2630 If you want a better approximation of the square root, then use:
2631
2632         $x = Math::BigFloat->new(12);
2633         $Math::BigFloat::precision = 0;
2634         Math::BigFloat->round_mode('even');
2635         print $x->copy->bsqrt(),"\n";           # 4
2636
2637         $Math::BigFloat::precision = 2;
2638         print $x->bsqrt(),"\n";                 # 3.46
2639         print $x->bsqrt(3),"\n";                # 3.464
2640
2641 =back
2642
2643 =head1 LICENSE
2644
2645 This program is free software; you may redistribute it and/or modify it under
2646 the same terms as Perl itself.
2647
2648 =head1 AUTHORS
2649
2650 Original code by Mark Biggar, overloaded interface by Ilya Zakharevich.
2651 Completely rewritten by Tels http://bloodgate.com in late 2000, 2001.
2652
2653 =cut