1package Math::BigFloat;
2
3#
4# Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
5#
6
7# The following hash values are internally used:
8#   _e	: exponent (ref to $CALC object)
9#   _m	: mantissa (ref to $CALC object)
10#   _es	: sign of _e
11# sign	: +,-,+inf,-inf, or "NaN" if not a number
12#   _a	: accuracy
13#   _p	: precision
14
15$VERSION = '1.51';
16require 5.005;
17
18require Exporter;
19@ISA =       qw(Exporter Math::BigInt);
20
21use strict;
22# $_trap_inf/$_trap_nan are internal and should never be accessed from outside
23use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode
24	    $upgrade $downgrade $_trap_nan $_trap_inf/;
25my $class = "Math::BigFloat";
26
27use overload
28'<=>'	=>	sub { $_[2] ?
29                      ref($_[0])->bcmp($_[1],$_[0]) :
30                      ref($_[0])->bcmp($_[0],$_[1])},
31'int'	=>	sub { $_[0]->as_number() },		# 'trunc' to bigint
32;
33
34##############################################################################
35# global constants, flags and assorted stuff
36
37# the following are public, but their usage is not recommended. Use the
38# accessor methods instead.
39
40# class constants, use Class->constant_name() to access
41$round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
42$accuracy   = undef;
43$precision  = undef;
44$div_scale  = 40;
45
46$upgrade = undef;
47$downgrade = undef;
48# the package we are using for our private parts, defaults to:
49# Math::BigInt->config()->{lib}
50my $MBI = 'Math::BigInt::FastCalc';
51
52# are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
53$_trap_nan = 0;
54# the same for infinity
55$_trap_inf = 0;
56
57# constant for easier life
58my $nan = 'NaN';
59
60my $IMPORT = 0;	# was import() called yet? used to make require work
61
62# some digits of accuracy for blog(undef,10); which we use in blog() for speed
63my $LOG_10 =
64 '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
65my $LOG_10_A = length($LOG_10)-1;
66# ditto for log(2)
67my $LOG_2 =
68 '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
69my $LOG_2_A = length($LOG_2)-1;
70my $HALF = '0.5';			# made into an object if necc.
71
72##############################################################################
73# the old code had $rnd_mode, so we need to support it, too
74
75sub TIESCALAR   { my ($class) = @_; bless \$round_mode, $class; }
76sub FETCH       { return $round_mode; }
77sub STORE       { $rnd_mode = $_[0]->round_mode($_[1]); }
78
79BEGIN
80  {
81  # when someone set's $rnd_mode, we catch this and check the value to see
82  # whether it is valid or not.
83  $rnd_mode   = 'even'; tie $rnd_mode, 'Math::BigFloat';
84  }
85
86##############################################################################
87
88{
89  # valid method aliases for AUTOLOAD
90  my %methods = map { $_ => 1 }
91   qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
92        fint facmp fcmp fzero fnan finf finc fdec flog ffac fneg
93	fceil ffloor frsft flsft fone flog froot
94      /;
95  # valid method's that can be hand-ed up (for AUTOLOAD)
96  my %hand_ups = map { $_ => 1 }
97   qw / is_nan is_inf is_negative is_positive is_pos is_neg
98        accuracy precision div_scale round_mode fabs fnot
99        objectify upgrade downgrade
100	bone binf bnan bzero
101      /;
102
103  sub method_alias { exists $methods{$_[0]||''}; }
104  sub method_hand_up { exists $hand_ups{$_[0]||''}; }
105}
106
107##############################################################################
108# constructors
109
110sub new
111  {
112  # create a new BigFloat object from a string or another bigfloat object.
113  # _e: exponent
114  # _m: mantissa
115  # sign  => sign (+/-), or "NaN"
116
117  my ($class,$wanted,@r) = @_;
118
119  # avoid numify-calls by not using || on $wanted!
120  return $class->bzero() if !defined $wanted;	# default to 0
121  return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
122
123  $class->import() if $IMPORT == 0;             # make require work
124
125  my $self = {}; bless $self, $class;
126  # shortcut for bigints and its subclasses
127  if ((ref($wanted)) && (ref($wanted) ne $class))
128    {
129    $self->{_m} = $wanted->as_number()->{value}; # get us a bigint copy
130    $self->{_e} = $MBI->_zero();
131    $self->{_es} = '+';
132    $self->{sign} = $wanted->sign();
133    return $self->bnorm();
134    }
135  # else: got a string
136
137  # handle '+inf', '-inf' first
138  if ($wanted =~ /^[+-]?inf\z/)
139    {
140    return $downgrade->new($wanted) if $downgrade;
141
142    $self->{sign} = $wanted;		# set a default sign for bstr()
143    return $self->binf($wanted);
144    }
145
146  # shortcut for simple forms like '12' that neither have trailing nor leading
147  # zeros
148  if ($wanted =~ /^([+-]?)([1-9][0-9]*[1-9])$/)
149    {
150    $self->{_e} = $MBI->_zero();
151    $self->{_es} = '+';
152    $self->{sign} = $1 || '+';
153    $self->{_m} = $MBI->_new($2);
154    return $self->round(@r) if !$downgrade;
155    }
156
157  my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
158  if (!ref $mis)
159    {
160    if ($_trap_nan)
161      {
162      require Carp;
163      Carp::croak ("$wanted is not a number initialized to $class");
164      }
165
166    return $downgrade->bnan() if $downgrade;
167
168    $self->{_e} = $MBI->_zero();
169    $self->{_es} = '+';
170    $self->{_m} = $MBI->_zero();
171    $self->{sign} = $nan;
172    }
173  else
174    {
175    # make integer from mantissa by adjusting exp, then convert to int
176    $self->{_e} = $MBI->_new($$ev);		# exponent
177    $self->{_es} = $$es || '+';
178    my $mantissa = "$$miv$$mfv"; 		# create mant.
179    $mantissa =~ s/^0+(\d)/$1/;			# strip leading zeros
180    $self->{_m} = $MBI->_new($mantissa); 	# create mant.
181
182    # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
183    if (CORE::length($$mfv) != 0)
184      {
185      my $len = $MBI->_new( CORE::length($$mfv));
186      ($self->{_e}, $self->{_es}) =
187	_e_sub ($self->{_e}, $len, $self->{_es}, '+');
188      }
189    # we can only have trailing zeros on the mantissa if $$mfv eq ''
190    else
191      {
192      # Use a regexp to count the trailing zeros in $$miv instead of _zeros()
193      # because that is faster, especially when _m is not stored in base 10.
194      my $zeros = 0; $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/;
195      if ($zeros != 0)
196        {
197        my $z = $MBI->_new($zeros);
198        # turn '120e2' into '12e3'
199        $MBI->_rsft ( $self->{_m}, $z, 10);
200        ($self->{_e}, $self->{_es}) =
201	  _e_add ( $self->{_e}, $z, $self->{_es}, '+');
202        }
203      }
204    $self->{sign} = $$mis;
205
206    # for something like 0Ey, set y to 1, and -0 => +0
207    # Check $$miv for beeing '0' and $$mfv eq '', because otherwise _m could not
208    # have become 0. That's faster than to call $MBI->_is_zero().
209    $self->{sign} = '+', $self->{_e} = $MBI->_one()
210     if $$miv eq '0' and $$mfv eq '';
211
212    return $self->round(@r) if !$downgrade;
213    }
214  # if downgrade, inf, NaN or integers go down
215
216  if ($downgrade && $self->{_es} eq '+')
217    {
218    if ($MBI->_is_zero( $self->{_e} ))
219      {
220      return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
221      }
222    return $downgrade->new($self->bsstr());
223    }
224  $self->bnorm()->round(@r);			# first normalize, then round
225  }
226
227sub copy
228  {
229  my ($c,$x);
230  if (@_ > 1)
231    {
232    # if two arguments, the first one is the class to "swallow" subclasses
233    ($c,$x) = @_;
234    }
235  else
236    {
237    $x = shift;
238    $c = ref($x);
239    }
240  return unless ref($x); # only for objects
241
242  my $self = {}; bless $self,$c;
243
244  $self->{sign} = $x->{sign};
245  $self->{_es} = $x->{_es};
246  $self->{_m} = $MBI->_copy($x->{_m});
247  $self->{_e} = $MBI->_copy($x->{_e});
248  $self->{_a} = $x->{_a} if defined $x->{_a};
249  $self->{_p} = $x->{_p} if defined $x->{_p};
250  $self;
251  }
252
253sub _bnan
254  {
255  # used by parent class bone() to initialize number to NaN
256  my $self = shift;
257
258  if ($_trap_nan)
259    {
260    require Carp;
261    my $class = ref($self);
262    Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
263    }
264
265  $IMPORT=1;					# call our import only once
266  $self->{_m} = $MBI->_zero();
267  $self->{_e} = $MBI->_zero();
268  $self->{_es} = '+';
269  }
270
271sub _binf
272  {
273  # used by parent class bone() to initialize number to +-inf
274  my $self = shift;
275
276  if ($_trap_inf)
277    {
278    require Carp;
279    my $class = ref($self);
280    Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
281    }
282
283  $IMPORT=1;					# call our import only once
284  $self->{_m} = $MBI->_zero();
285  $self->{_e} = $MBI->_zero();
286  $self->{_es} = '+';
287  }
288
289sub _bone
290  {
291  # used by parent class bone() to initialize number to 1
292  my $self = shift;
293  $IMPORT=1;					# call our import only once
294  $self->{_m} = $MBI->_one();
295  $self->{_e} = $MBI->_zero();
296  $self->{_es} = '+';
297  }
298
299sub _bzero
300  {
301  # used by parent class bone() to initialize number to 0
302  my $self = shift;
303  $IMPORT=1;					# call our import only once
304  $self->{_m} = $MBI->_zero();
305  $self->{_e} = $MBI->_one();
306  $self->{_es} = '+';
307  }
308
309sub isa
310  {
311  my ($self,$class) = @_;
312  return if $class =~ /^Math::BigInt/;		# we aren't one of these
313  UNIVERSAL::isa($self,$class);
314  }
315
316sub config
317  {
318  # return (later set?) configuration data as hash ref
319  my $class = shift || 'Math::BigFloat';
320
321  my $cfg = $class->SUPER::config(@_);
322
323  # now we need only to override the ones that are different from our parent
324  $cfg->{class} = $class;
325  $cfg->{with} = $MBI;
326  $cfg;
327  }
328
329##############################################################################
330# string conversation
331
332sub bstr
333  {
334  # (ref to BFLOAT or num_str ) return num_str
335  # Convert number from internal format to (non-scientific) string format.
336  # internal format is always normalized (no leading zeros, "-0" => "+0")
337  my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
338
339  if ($x->{sign} !~ /^[+-]$/)
340    {
341    return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
342    return 'inf';                                       # +inf
343    }
344
345  my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
346
347  # $x is zero?
348  my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
349  if ($not_zero)
350    {
351    $es = $MBI->_str($x->{_m});
352    $len = CORE::length($es);
353    my $e = $MBI->_num($x->{_e});
354    $e = -$e if $x->{_es} eq '-';
355    if ($e < 0)
356      {
357      $dot = '';
358      # if _e is bigger than a scalar, the following will blow your memory
359      if ($e <= -$len)
360        {
361        my $r = abs($e) - $len;
362        $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
363        }
364      else
365        {
366        substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
367        $cad = -$cad if $x->{_es} eq '-';
368        }
369      }
370    elsif ($e > 0)
371      {
372      # expand with zeros
373      $es .= '0' x $e; $len += $e; $cad = 0;
374      }
375    } # if not zero
376
377  $es = '-'.$es if $x->{sign} eq '-';
378  # if set accuracy or precision, pad with zeros on the right side
379  if ((defined $x->{_a}) && ($not_zero))
380    {
381    # 123400 => 6, 0.1234 => 4, 0.001234 => 4
382    my $zeros = $x->{_a} - $cad;		# cad == 0 => 12340
383    $zeros = $x->{_a} - $len if $cad != $len;
384    $es .= $dot.'0' x $zeros if $zeros > 0;
385    }
386  elsif ((($x->{_p} || 0) < 0))
387    {
388    # 123400 => 6, 0.1234 => 4, 0.001234 => 6
389    my $zeros = -$x->{_p} + $cad;
390    $es .= $dot.'0' x $zeros if $zeros > 0;
391    }
392  $es;
393  }
394
395sub bsstr
396  {
397  # (ref to BFLOAT or num_str ) return num_str
398  # Convert number from internal format to scientific string format.
399  # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
400  my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
401
402  if ($x->{sign} !~ /^[+-]$/)
403    {
404    return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
405    return 'inf';                                       # +inf
406    }
407  my $sep = 'e'.$x->{_es};
408  my $sign = $x->{sign}; $sign = '' if $sign eq '+';
409  $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
410  }
411
412sub numify
413  {
414  # Make a number from a BigFloat object
415  # simple return a string and let Perl's atoi()/atof() handle the rest
416  my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
417  $x->bsstr();
418  }
419
420##############################################################################
421# public stuff (usually prefixed with "b")
422
423sub bneg
424  {
425  # (BINT or num_str) return BINT
426  # negate number or make a negated number from string
427  my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
428
429  return $x if $x->modify('bneg');
430
431  # for +0 dont negate (to have always normalized +0). Does nothing for 'NaN'
432  $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
433  $x;
434  }
435
436# tels 2001-08-04
437# XXX TODO this must be overwritten and return NaN for non-integer values
438# band(), bior(), bxor(), too
439#sub bnot
440#  {
441#  $class->SUPER::bnot($class,@_);
442#  }
443
444sub bcmp
445  {
446  # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
447
448  # set up parameters
449  my ($self,$x,$y) = (ref($_[0]),@_);
450  # objectify is costly, so avoid it
451  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
452    {
453    ($self,$x,$y) = objectify(2,@_);
454    }
455
456  return $upgrade->bcmp($x,$y) if defined $upgrade &&
457    ((!$x->isa($self)) || (!$y->isa($self)));
458
459  if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
460    {
461    # handle +-inf and NaN
462    return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
463    return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
464    return +1 if $x->{sign} eq '+inf';
465    return -1 if $x->{sign} eq '-inf';
466    return -1 if $y->{sign} eq '+inf';
467    return +1;
468    }
469
470  # check sign for speed first
471  return 1 if $x->{sign} eq '+' && $y->{sign} eq '-';	# does also 0 <=> -y
472  return -1 if $x->{sign} eq '-' && $y->{sign} eq '+';	# does also -x <=> 0
473
474  # shortcut
475  my $xz = $x->is_zero();
476  my $yz = $y->is_zero();
477  return 0 if $xz && $yz;				# 0 <=> 0
478  return -1 if $xz && $y->{sign} eq '+';		# 0 <=> +y
479  return 1 if $yz && $x->{sign} eq '+';			# +x <=> 0
480
481  # adjust so that exponents are equal
482  my $lxm = $MBI->_len($x->{_m});
483  my $lym = $MBI->_len($y->{_m});
484  # the numify somewhat limits our length, but makes it much faster
485  my ($xes,$yes) = (1,1);
486  $xes = -1 if $x->{_es} ne '+';
487  $yes = -1 if $y->{_es} ne '+';
488  my $lx = $lxm + $xes * $MBI->_num($x->{_e});
489  my $ly = $lym + $yes * $MBI->_num($y->{_e});
490  my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
491  return $l <=> 0 if $l != 0;
492
493  # lengths (corrected by exponent) are equal
494  # so make mantissa equal length by padding with zero (shift left)
495  my $diff = $lxm - $lym;
496  my $xm = $x->{_m};		# not yet copy it
497  my $ym = $y->{_m};
498  if ($diff > 0)
499    {
500    $ym = $MBI->_copy($y->{_m});
501    $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
502    }
503  elsif ($diff < 0)
504    {
505    $xm = $MBI->_copy($x->{_m});
506    $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
507    }
508  my $rc = $MBI->_acmp($xm,$ym);
509  $rc = -$rc if $x->{sign} eq '-';		# -124 < -123
510  $rc <=> 0;
511  }
512
513sub bacmp
514  {
515  # Compares 2 values, ignoring their signs.
516  # Returns one of undef, <0, =0, >0. (suitable for sort)
517
518  # set up parameters
519  my ($self,$x,$y) = (ref($_[0]),@_);
520  # objectify is costly, so avoid it
521  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
522    {
523    ($self,$x,$y) = objectify(2,@_);
524    }
525
526  return $upgrade->bacmp($x,$y) if defined $upgrade &&
527    ((!$x->isa($self)) || (!$y->isa($self)));
528
529  # handle +-inf and NaN's
530  if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
531    {
532    return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
533    return 0 if ($x->is_inf() && $y->is_inf());
534    return 1 if ($x->is_inf() && !$y->is_inf());
535    return -1;
536    }
537
538  # shortcut
539  my $xz = $x->is_zero();
540  my $yz = $y->is_zero();
541  return 0 if $xz && $yz;				# 0 <=> 0
542  return -1 if $xz && !$yz;				# 0 <=> +y
543  return 1 if $yz && !$xz;				# +x <=> 0
544
545  # adjust so that exponents are equal
546  my $lxm = $MBI->_len($x->{_m});
547  my $lym = $MBI->_len($y->{_m});
548  my ($xes,$yes) = (1,1);
549  $xes = -1 if $x->{_es} ne '+';
550  $yes = -1 if $y->{_es} ne '+';
551  # the numify somewhat limits our length, but makes it much faster
552  my $lx = $lxm + $xes * $MBI->_num($x->{_e});
553  my $ly = $lym + $yes * $MBI->_num($y->{_e});
554  my $l = $lx - $ly;
555  return $l <=> 0 if $l != 0;
556
557  # lengths (corrected by exponent) are equal
558  # so make mantissa equal-length by padding with zero (shift left)
559  my $diff = $lxm - $lym;
560  my $xm = $x->{_m};		# not yet copy it
561  my $ym = $y->{_m};
562  if ($diff > 0)
563    {
564    $ym = $MBI->_copy($y->{_m});
565    $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
566    }
567  elsif ($diff < 0)
568    {
569    $xm = $MBI->_copy($x->{_m});
570    $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
571    }
572  $MBI->_acmp($xm,$ym);
573  }
574
575sub badd
576  {
577  # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
578  # return result as BFLOAT
579
580  # set up parameters
581  my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
582  # objectify is costly, so avoid it
583  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
584    {
585    ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
586    }
587
588  # inf and NaN handling
589  if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
590    {
591    # NaN first
592    return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
593    # inf handling
594    if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
595      {
596      # +inf++inf or -inf+-inf => same, rest is NaN
597      return $x if $x->{sign} eq $y->{sign};
598      return $x->bnan();
599      }
600    # +-inf + something => +inf; something +-inf => +-inf
601    $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
602    return $x;
603    }
604
605  return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
606   ((!$x->isa($self)) || (!$y->isa($self)));
607
608  # speed: no add for 0+y or x+0
609  return $x->bround($a,$p,$r) if $y->is_zero();		# x+0
610  if ($x->is_zero())					# 0+y
611    {
612    # make copy, clobbering up x (modify in place!)
613    $x->{_e} = $MBI->_copy($y->{_e});
614    $x->{_es} = $y->{_es};
615    $x->{_m} = $MBI->_copy($y->{_m});
616    $x->{sign} = $y->{sign} || $nan;
617    return $x->round($a,$p,$r,$y);
618    }
619
620  # take lower of the two e's and adapt m1 to it to match m2
621  my $e = $y->{_e};
622  $e = $MBI->_zero() if !defined $e;		# if no BFLOAT?
623  $e = $MBI->_copy($e);				# make copy (didn't do it yet)
624
625  my $es;
626
627  ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});
628
629  my $add = $MBI->_copy($y->{_m});
630
631  if ($es eq '-')				# < 0
632    {
633    $MBI->_lsft( $x->{_m}, $e, 10);
634    ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
635    }
636  elsif (!$MBI->_is_zero($e))			# > 0
637    {
638    $MBI->_lsft($add, $e, 10);
639    }
640  # else: both e are the same, so just leave them
641
642  if ($x->{sign} eq $y->{sign})
643    {
644    # add
645    $x->{_m} = $MBI->_add($x->{_m}, $add);
646    }
647  else
648    {
649    ($x->{_m}, $x->{sign}) =
650     _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
651    }
652
653  # delete trailing zeros, then round
654  $x->bnorm()->round($a,$p,$r,$y);
655  }
656
657# sub bsub is inherited from Math::BigInt!
658
659sub binc
660  {
661  # increment arg by one
662  my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
663
664  if ($x->{_es} eq '-')
665    {
666    return $x->badd($self->bone(),@r);	#  digits after dot
667    }
668
669  if (!$MBI->_is_zero($x->{_e}))		# _e == 0 for NaN, inf, -inf
670    {
671    # 1e2 => 100, so after the shift below _m has a '0' as last digit
672    $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10);	# 1e2 => 100
673    $x->{_e} = $MBI->_zero();				# normalize
674    $x->{_es} = '+';
675    # we know that the last digit of $x will be '1' or '9', depending on the
676    # sign
677    }
678  # now $x->{_e} == 0
679  if ($x->{sign} eq '+')
680    {
681    $MBI->_inc($x->{_m});
682    return $x->bnorm()->bround(@r);
683    }
684  elsif ($x->{sign} eq '-')
685    {
686    $MBI->_dec($x->{_m});
687    $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
688    return $x->bnorm()->bround(@r);
689    }
690  # inf, nan handling etc
691  $x->badd($self->bone(),@r);			# badd() does round
692  }
693
694sub bdec
695  {
696  # decrement arg by one
697  my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
698
699  if ($x->{_es} eq '-')
700    {
701    return $x->badd($self->bone('-'),@r);	#  digits after dot
702    }
703
704  if (!$MBI->_is_zero($x->{_e}))
705    {
706    $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10);	# 1e2 => 100
707    $x->{_e} = $MBI->_zero();				# normalize
708    $x->{_es} = '+';
709    }
710  # now $x->{_e} == 0
711  my $zero = $x->is_zero();
712  # <= 0
713  if (($x->{sign} eq '-') || $zero)
714    {
715    $MBI->_inc($x->{_m});
716    $x->{sign} = '-' if $zero;				# 0 => 1 => -1
717    $x->{sign} = '+' if $MBI->_is_zero($x->{_m});	# -1 +1 => -0 => +0
718    return $x->bnorm()->round(@r);
719    }
720  # > 0
721  elsif ($x->{sign} eq '+')
722    {
723    $MBI->_dec($x->{_m});
724    return $x->bnorm()->round(@r);
725    }
726  # inf, nan handling etc
727  $x->badd($self->bone('-'),@r);		# does round
728  }
729
730sub DEBUG () { 0; }
731
732sub blog
733  {
734  my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
735
736  # $base > 0, $base != 1; if $base == undef default to $base == e
737  # $x >= 0
738
739  # we need to limit the accuracy to protect against overflow
740  my $fallback = 0;
741  my ($scale,@params);
742  ($x,@params) = $x->_find_round_parameters($a,$p,$r);
743
744  # also takes care of the "error in _find_round_parameters?" case
745  return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
746
747
748  # no rounding at all, so must use fallback
749  if (scalar @params == 0)
750    {
751    # simulate old behaviour
752    $params[0] = $self->div_scale();	# and round to it as accuracy
753    $params[1] = undef;			# P = undef
754    $scale = $params[0]+4; 		# at least four more for proper round
755    $params[2] = $r;			# round mode by caller or undef
756    $fallback = 1;			# to clear a/p afterwards
757    }
758  else
759    {
760    # the 4 below is empirical, and there might be cases where it is not
761    # enough...
762    $scale = abs($params[0] || $params[1]) + 4;	# take whatever is defined
763    }
764
765  return $x->bzero(@params) if $x->is_one();
766  # base not defined => base == Euler's constant e
767  if (defined $base)
768    {
769    # make object, since we don't feed it through objectify() to still get the
770    # case of $base == undef
771    $base = $self->new($base) unless ref($base);
772    # $base > 0; $base != 1
773    return $x->bnan() if $base->is_zero() || $base->is_one() ||
774      $base->{sign} ne '+';
775    # if $x == $base, we know the result must be 1.0
776    if ($x->bcmp($base) == 0)
777      {
778      $x->bone('+',@params);
779      if ($fallback)
780        {
781        # clear a/p after round, since user did not request it
782        delete $x->{_a}; delete $x->{_p};
783        }
784      return $x;
785      }
786    }
787
788  # when user set globals, they would interfere with our calculation, so
789  # disable them and later re-enable them
790  no strict 'refs';
791  my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
792  my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
793  # we also need to disable any set A or P on $x (_find_round_parameters took
794  # them already into account), since these would interfere, too
795  delete $x->{_a}; delete $x->{_p};
796  # need to disable $upgrade in BigInt, to avoid deep recursion
797  local $Math::BigInt::upgrade = undef;
798  local $Math::BigFloat::downgrade = undef;
799
800  # upgrade $x if $x is not a BigFloat (handle BigInt input)
801  if (!$x->isa('Math::BigFloat'))
802    {
803    $x = Math::BigFloat->new($x);
804    $self = ref($x);
805    }
806
807  my $done = 0;
808
809  # If the base is defined and an integer, try to calculate integer result
810  # first. This is very fast, and in case the real result was found, we can
811  # stop right here.
812  if (defined $base && $base->is_int() && $x->is_int())
813    {
814    my $i = $MBI->_copy( $x->{_m} );
815    $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
816    my $int = Math::BigInt->bzero();
817    $int->{value} = $i;
818    $int->blog($base->as_number());
819    # if ($exact)
820    if ($base->as_number()->bpow($int) == $x)
821      {
822      # found result, return it
823      $x->{_m} = $int->{value};
824      $x->{_e} = $MBI->_zero();
825      $x->{_es} = '+';
826      $x->bnorm();
827      $done = 1;
828      }
829    }
830
831  if ($done == 0)
832    {
833    # first calculate the log to base e (using reduction by 10 (and probably 2))
834    $self->_log_10($x,$scale);
835
836    # and if a different base was requested, convert it
837    if (defined $base)
838      {
839      $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
840      # not ln, but some other base (don't modify $base)
841      $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
842      }
843    }
844
845  # shortcut to not run through _find_round_parameters again
846  if (defined $params[0])
847    {
848    $x->bround($params[0],$params[2]);		# then round accordingly
849    }
850  else
851    {
852    $x->bfround($params[1],$params[2]);		# then round accordingly
853    }
854  if ($fallback)
855    {
856    # clear a/p after round, since user did not request it
857    delete $x->{_a}; delete $x->{_p};
858    }
859  # restore globals
860  $$abr = $ab; $$pbr = $pb;
861
862  $x;
863  }
864
865sub _log
866  {
867  # internal log function to calculate ln() based on Taylor series.
868  # Modifies $x in place.
869  my ($self,$x,$scale) = @_;
870
871  # in case of $x == 1, result is 0
872  return $x->bzero() if $x->is_one();
873
874  # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
875
876  # u = x-1, v = x+1
877  #              _                               _
878  # Taylor:     |    u    1   u^3   1   u^5       |
879  # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 0
880  #             |_   v    3   v^3   5   v^5      _|
881
882  # This takes much more steps to calculate the result and is thus not used
883  # u = x-1
884  #              _                               _
885  # Taylor:     |    u    1   u^2   1   u^3       |
886  # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 1/2
887  #             |_   x    2   x^2   3   x^3      _|
888
889  my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
890
891  $v = $x->copy(); $v->binc();		# v = x+1
892  $x->bdec(); $u = $x->copy();		# u = x-1; x = x-1
893  $x->bdiv($v,$scale);			# first term: u/v
894  $below = $v->copy();
895  $over = $u->copy();
896  $u *= $u; $v *= $v;				# u^2, v^2
897  $below->bmul($v);				# u^3, v^3
898  $over->bmul($u);
899  $factor = $self->new(3); $f = $self->new(2);
900
901  my $steps = 0 if DEBUG;
902  $limit = $self->new("1E-". ($scale-1));
903  while (3 < 5)
904    {
905    # we calculate the next term, and add it to the last
906    # when the next term is below our limit, it won't affect the outcome
907    # anymore, so we stop
908
909    # calculating the next term simple from over/below will result in quite
910    # a time hog if the input has many digits, since over and below will
911    # accumulate more and more digits, and the result will also have many
912    # digits, but in the end it is rounded to $scale digits anyway. So if we
913    # round $over and $below first, we save a lot of time for the division
914    # (not with log(1.2345), but try log (123**123) to see what I mean. This
915    # can introduce a rounding error if the division result would be f.i.
916    # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
917    # if we truncated $over and $below we might get 0.12345. Does this matter
918    # for the end result? So we give $over and $below 4 more digits to be
919    # on the safe side (unscientific error handling as usual... :+D
920
921    $next = $over->copy->bround($scale+4)->bdiv(
922      $below->copy->bmul($factor)->bround($scale+4),
923      $scale);
924
925## old version:
926##    $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
927
928    last if $next->bacmp($limit) <= 0;
929
930    delete $next->{_a}; delete $next->{_p};
931    $x->badd($next);
932    # calculate things for the next term
933    $over *= $u; $below *= $v; $factor->badd($f);
934    if (DEBUG)
935      {
936      $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
937      }
938    }
939  $x->bmul($f);					# $x *= 2
940  print "took $steps steps\n" if DEBUG;
941  }
942
943sub _log_10
944  {
945  # Internal log function based on reducing input to the range of 0.1 .. 9.99
946  # and then "correcting" the result to the proper one. Modifies $x in place.
947  my ($self,$x,$scale) = @_;
948
949  # taking blog() from numbers greater than 10 takes a *very long* time, so we
950  # break the computation down into parts based on the observation that:
951  #  blog(x*y) = blog(x) + blog(y)
952  # We set $y here to multiples of 10 so that $x is below 1 (the smaller $x is
953  # the faster it get's, especially because 2*$x takes about 10 times as long,
954  # so by dividing $x by 10 we make it at least factor 100 faster...)
955
956  # The same observation is valid for numbers smaller than 0.1 (e.g. computing
957  # log(1) is fastest, and the farther away we get from 1, the longer it takes)
958  # so we also 'break' this down by multiplying $x with 10 and subtract the
959  # log(10) afterwards to get the correct result.
960
961  # calculate nr of digits before dot
962  my $dbd = $MBI->_num($x->{_e});
963  $dbd = -$dbd if $x->{_es} eq '-';
964  $dbd += $MBI->_len($x->{_m});
965
966  # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
967  # infinite recursion
968
969  my $calc = 1;					# do some calculation?
970
971  # disable the shortcut for 10, since we need log(10) and this would recurse
972  # infinitely deep
973  if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
974    {
975    $dbd = 0;					# disable shortcut
976    # we can use the cached value in these cases
977    if ($scale <= $LOG_10_A)
978      {
979      $x->bzero(); $x->badd($LOG_10);
980      $calc = 0; 				# no need to calc, but round
981      }
982    }
983  else
984    {
985    # disable the shortcut for 2, since we maybe have it cached
986    if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
987      {
988      $dbd = 0;					# disable shortcut
989      # we can use the cached value in these cases
990      if ($scale <= $LOG_2_A)
991        {
992        $x->bzero(); $x->badd($LOG_2);
993        $calc = 0; 				# no need to calc, but round
994        }
995      }
996    }
997
998  # if $x = 0.1, we know the result must be 0-log(10)
999  if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
1000      $MBI->_is_one($x->{_m}))
1001    {
1002    $dbd = 0;					# disable shortcut
1003    # we can use the cached value in these cases
1004    if ($scale <= $LOG_10_A)
1005      {
1006      $x->bzero(); $x->bsub($LOG_10);
1007      $calc = 0; 				# no need to calc, but round
1008      }
1009    }
1010
1011  return if $calc == 0;				# already have the result
1012
1013  # default: these correction factors are undef and thus not used
1014  my $l_10;				# value of ln(10) to A of $scale
1015  my $l_2;				# value of ln(2) to A of $scale
1016
1017  # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
1018  # so don't do this shortcut for 1 or 0
1019  if (($dbd > 1) || ($dbd < 0))
1020    {
1021    # convert our cached value to an object if not already (avoid doing this
1022    # at import() time, since not everybody needs this)
1023    $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
1024
1025    #print "x = $x, dbd = $dbd, calc = $calc\n";
1026    # got more than one digit before the dot, or more than one zero after the
1027    # dot, so do:
1028    #  log(123)    == log(1.23) + log(10) * 2
1029    #  log(0.0123) == log(1.23) - log(10) * 2
1030
1031    if ($scale <= $LOG_10_A)
1032      {
1033      # use cached value
1034      $l_10 = $LOG_10->copy();		# copy for mul
1035      }
1036    else
1037      {
1038      # else: slower, compute it (but don't cache it, because it could be big)
1039      # also disable downgrade for this code path
1040      local $Math::BigFloat::downgrade = undef;
1041      $l_10 = $self->new(10)->blog(undef,$scale);	# scale+4, actually
1042      }
1043    $dbd-- if ($dbd > 1); 		# 20 => dbd=2, so make it dbd=1
1044    $l_10->bmul( $self->new($dbd));	# log(10) * (digits_before_dot-1)
1045    my $dbd_sign = '+';
1046    if ($dbd < 0)
1047      {
1048      $dbd = -$dbd;
1049      $dbd_sign = '-';
1050      }
1051    ($x->{_e}, $x->{_es}) =
1052	_e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
1053
1054    }
1055
1056  # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1057
1058  ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1059  ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1060
1061  $HALF = $self->new($HALF) unless ref($HALF);
1062
1063  my $twos = 0;				# default: none (0 times)
1064  my $two = $self->new(2);
1065  while ($x->bacmp($HALF) <= 0)
1066    {
1067    $twos--; $x->bmul($two);
1068    }
1069  while ($x->bacmp($two) >= 0)
1070    {
1071    $twos++; $x->bdiv($two,$scale+4);		# keep all digits
1072    }
1073  # $twos > 0 => did mul 2, < 0 => did div 2 (never both)
1074  # calculate correction factor based on ln(2)
1075  if ($twos != 0)
1076    {
1077    $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1078    if ($scale <= $LOG_2_A)
1079      {
1080      # use cached value
1081      $l_2 = $LOG_2->copy();			# copy for mul
1082      }
1083    else
1084      {
1085      # else: slower, compute it (but don't cache it, because it could be big)
1086      # also disable downgrade for this code path
1087      local $Math::BigFloat::downgrade = undef;
1088      $l_2 = $two->blog(undef,$scale);	# scale+4, actually
1089      }
1090    $l_2->bmul($twos);		# * -2 => subtract, * 2 => add
1091    }
1092
1093  $self->_log($x,$scale);			# need to do the "normal" way
1094  $x->badd($l_10) if defined $l_10; 		# correct it by ln(10)
1095  $x->badd($l_2) if defined $l_2;		# and maybe by ln(2)
1096  # all done, $x contains now the result
1097  }
1098
1099sub blcm
1100  {
1101  # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1102  # does not modify arguments, but returns new object
1103  # Lowest Common Multiplicator
1104
1105  my ($self,@arg) = objectify(0,@_);
1106  my $x = $self->new(shift @arg);
1107  while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); }
1108  $x;
1109  }
1110
1111sub bgcd
1112  {
1113  # (BINT or num_str, BINT or num_str) return BINT
1114  # does not modify arguments, but returns new object
1115
1116  my $y = shift;
1117  $y = __PACKAGE__->new($y) if !ref($y);
1118  my $self = ref($y);
1119  my $x = $y->copy()->babs();			# keep arguments
1120
1121  return $x->bnan() if $x->{sign} !~ /^[+-]$/	# x NaN?
1122	|| !$x->is_int();			# only for integers now
1123
1124  while (@_)
1125    {
1126    my $t = shift; $t = $self->new($t) if !ref($t);
1127    $y = $t->copy()->babs();
1128
1129    return $x->bnan() if $y->{sign} !~ /^[+-]$/	# y NaN?
1130     	|| !$y->is_int();			# only for integers now
1131
1132    # greatest common divisor
1133    while (! $y->is_zero())
1134      {
1135      ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
1136      }
1137
1138    last if $x->is_one();
1139    }
1140  $x;
1141  }
1142
1143##############################################################################
1144
1145sub _e_add
1146  {
1147  # Internal helper sub to take two positive integers and their signs and
1148  # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1149  # output ($CALC,('+'|'-'))
1150  my ($x,$y,$xs,$ys) = @_;
1151
1152  # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1153  if ($xs eq $ys)
1154    {
1155    $x = $MBI->_add ($x, $y );		# a+b
1156    # the sign follows $xs
1157    return ($x, $xs);
1158    }
1159
1160  my $a = $MBI->_acmp($x,$y);
1161  if ($a > 0)
1162    {
1163    $x = $MBI->_sub ($x , $y);				# abs sub
1164    }
1165  elsif ($a == 0)
1166    {
1167    $x = $MBI->_zero();					# result is 0
1168    $xs = '+';
1169    }
1170  else # a < 0
1171    {
1172    $x = $MBI->_sub ( $y, $x, 1 );			# abs sub
1173    $xs = $ys;
1174    }
1175  ($x,$xs);
1176  }
1177
1178sub _e_sub
1179  {
1180  # Internal helper sub to take two positive integers and their signs and
1181  # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1182  # output ($CALC,('+'|'-'))
1183  my ($x,$y,$xs,$ys) = @_;
1184
1185  # flip sign
1186  $ys =~ tr/+-/-+/;
1187  _e_add($x,$y,$xs,$ys);		# call add (does subtract now)
1188  }
1189
1190###############################################################################
1191# is_foo methods (is_negative, is_positive are inherited from BigInt)
1192
1193sub is_int
1194  {
1195  # return true if arg (BFLOAT or num_str) is an integer
1196  my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1197
1198  return 1 if ($x->{sign} =~ /^[+-]$/) &&	# NaN and +-inf aren't
1199    $x->{_es} eq '+';				# 1e-1 => no integer
1200  0;
1201  }
1202
1203sub is_zero
1204  {
1205  # return true if arg (BFLOAT or num_str) is zero
1206  my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1207
1208  return 1 if $x->{sign} eq '+' && $MBI->_is_zero($x->{_m});
1209  0;
1210  }
1211
1212sub is_one
1213  {
1214  # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1215  my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1216
1217  $sign = '+' if !defined $sign || $sign ne '-';
1218  return 1
1219   if ($x->{sign} eq $sign &&
1220    $MBI->_is_zero($x->{_e}) && $MBI->_is_one($x->{_m}));
1221  0;
1222  }
1223
1224sub is_odd
1225  {
1226  # return true if arg (BFLOAT or num_str) is odd or false if even
1227  my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1228
1229  return 1 if ($x->{sign} =~ /^[+-]$/) &&		# NaN & +-inf aren't
1230    ($MBI->_is_zero($x->{_e}) && $MBI->_is_odd($x->{_m}));
1231  0;
1232  }
1233
1234sub is_even
1235  {
1236  # return true if arg (BINT or num_str) is even or false if odd
1237  my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1238
1239  return 0 if $x->{sign} !~ /^[+-]$/;			# NaN & +-inf aren't
1240  return 1 if ($x->{_es} eq '+'	 			# 123.45 is never
1241     && $MBI->_is_even($x->{_m}));			# but 1200 is
1242  0;
1243  }
1244
1245sub bmul
1246  {
1247  # multiply two numbers -- stolen from Knuth Vol 2 pg 233
1248  # (BINT or num_str, BINT or num_str) return BINT
1249
1250  # set up parameters
1251  my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1252  # objectify is costly, so avoid it
1253  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1254    {
1255    ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1256    }
1257
1258  return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1259
1260  # inf handling
1261  if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1262    {
1263    return $x->bnan() if $x->is_zero() || $y->is_zero();
1264    # result will always be +-inf:
1265    # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1266    # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1267    return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1268    return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1269    return $x->binf('-');
1270    }
1271  # handle result = 0
1272  return $x->bzero() if $x->is_zero() || $y->is_zero();
1273
1274  return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
1275   ((!$x->isa($self)) || (!$y->isa($self)));
1276
1277  # aEb * cEd = (a*c)E(b+d)
1278  $MBI->_mul($x->{_m},$y->{_m});
1279  ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1280
1281  # adjust sign:
1282  $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1283  return $x->bnorm()->round($a,$p,$r,$y);
1284  }
1285
1286sub bdiv
1287  {
1288  # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
1289  # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
1290
1291  # set up parameters
1292  my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1293  # objectify is costly, so avoid it
1294  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1295    {
1296    ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1297    }
1298
1299  return $self->_div_inf($x,$y)
1300   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
1301
1302  # x== 0 # also: or y == 1 or y == -1
1303  return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1304
1305  # upgrade ?
1306  return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1307
1308  # we need to limit the accuracy to protect against overflow
1309  my $fallback = 0;
1310  my (@params,$scale);
1311  ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1312
1313  return $x if $x->is_nan();		# error in _find_round_parameters?
1314
1315  # no rounding at all, so must use fallback
1316  if (scalar @params == 0)
1317    {
1318    # simulate old behaviour
1319    $params[0] = $self->div_scale();	# and round to it as accuracy
1320    $scale = $params[0]+4; 		# at least four more for proper round
1321    $params[2] = $r;			# round mode by caller or undef
1322    $fallback = 1;			# to clear a/p afterwards
1323    }
1324  else
1325    {
1326    # the 4 below is empirical, and there might be cases where it is not
1327    # enough...
1328    $scale = abs($params[0] || $params[1]) + 4;	# take whatever is defined
1329    }
1330
1331  my $rem; $rem = $self->bzero() if wantarray;
1332
1333  $y = $self->new($y) unless $y->isa('Math::BigFloat');
1334
1335  my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
1336  $scale = $lx if $lx > $scale;
1337  $scale = $ly if $ly > $scale;
1338  my $diff = $ly - $lx;
1339  $scale += $diff if $diff > 0;		# if lx << ly, but not if ly << lx!
1340
1341  # already handled inf/NaN/-inf above:
1342
1343  # check that $y is not 1 nor -1 and cache the result:
1344  my $y_not_one = !($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m}));
1345
1346  # flipping the sign of $y will also flip the sign of $x for the special
1347  # case of $x->bsub($x); so we can catch it below:
1348  my $xsign = $x->{sign};
1349  $y->{sign} =~ tr/+-/-+/;
1350
1351  if ($xsign ne $x->{sign})
1352    {
1353    # special case of $x /= $x results in 1
1354    $x->bone();			# "fixes" also sign of $y, since $x is $y
1355    }
1356  else
1357    {
1358    # correct $y's sign again
1359    $y->{sign} =~ tr/+-/-+/;
1360    # continue with normal div code:
1361
1362    # make copy of $x in case of list context for later reminder calculation
1363    if (wantarray && $y_not_one)
1364      {
1365      $rem = $x->copy();
1366      }
1367
1368    $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
1369
1370    # check for / +-1 ( +/- 1E0)
1371    if ($y_not_one)
1372      {
1373      # promote BigInts and it's subclasses (except when already a BigFloat)
1374      $y = $self->new($y) unless $y->isa('Math::BigFloat');
1375
1376      # calculate the result to $scale digits and then round it
1377      # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1378      $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1379      $MBI->_div ($x->{_m},$y->{_m});	# a/c
1380
1381      # correct exponent of $x
1382      ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1383      # correct for 10**scale
1384      ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1385      $x->bnorm();		# remove trailing 0's
1386      }
1387    } # ende else $x != $y
1388
1389  # shortcut to not run through _find_round_parameters again
1390  if (defined $params[0])
1391    {
1392    delete $x->{_a}; 				# clear before round
1393    $x->bround($params[0],$params[2]);		# then round accordingly
1394    }
1395  else
1396    {
1397    delete $x->{_p}; 				# clear before round
1398    $x->bfround($params[1],$params[2]);		# then round accordingly
1399    }
1400  if ($fallback)
1401    {
1402    # clear a/p after round, since user did not request it
1403    delete $x->{_a}; delete $x->{_p};
1404    }
1405
1406  if (wantarray)
1407    {
1408    if ($y_not_one)
1409      {
1410      $rem->bmod($y,@params);			# copy already done
1411      }
1412    if ($fallback)
1413      {
1414      # clear a/p after round, since user did not request it
1415      delete $rem->{_a}; delete $rem->{_p};
1416      }
1417    return ($x,$rem);
1418    }
1419  $x;
1420  }
1421
1422sub bmod
1423  {
1424  # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder
1425
1426  # set up parameters
1427  my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1428  # objectify is costly, so avoid it
1429  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1430    {
1431    ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1432    }
1433
1434  # handle NaN, inf, -inf
1435  if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1436    {
1437    my ($d,$re) = $self->SUPER::_div_inf($x,$y);
1438    $x->{sign} = $re->{sign};
1439    $x->{_e} = $re->{_e};
1440    $x->{_m} = $re->{_m};
1441    return $x->round($a,$p,$r,$y);
1442    }
1443  if ($y->is_zero())
1444    {
1445    return $x->bnan() if $x->is_zero();
1446    return $x;
1447    }
1448
1449  return $x->bzero() if $x->is_zero()
1450 || ($x->is_int() &&
1451  # check that $y == +1 or $y == -1:
1452    ($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m})));
1453
1454  my $cmp = $x->bacmp($y);			# equal or $x < $y?
1455  return $x->bzero($a,$p) if $cmp == 0;		# $x == $y => result 0
1456
1457  # only $y of the operands negative?
1458  my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1459
1460  $x->{sign} = $y->{sign};				# calc sign first
1461  return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0;	# $x < $y => result $x
1462
1463  my $ym = $MBI->_copy($y->{_m});
1464
1465  # 2e1 => 20
1466  $MBI->_lsft( $ym, $y->{_e}, 10)
1467   if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
1468
1469  # if $y has digits after dot
1470  my $shifty = 0;			# correct _e of $x by this
1471  if ($y->{_es} eq '-')			# has digits after dot
1472    {
1473    # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1474    $shifty = $MBI->_num($y->{_e}); 	# no more digits after dot
1475    $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
1476    }
1477  # $ym is now mantissa of $y based on exponent 0
1478
1479  my $shiftx = 0;			# correct _e of $x by this
1480  if ($x->{_es} eq '-')			# has digits after dot
1481    {
1482    # 123.4 % 20 => 1234 % 200
1483    $shiftx = $MBI->_num($x->{_e});	# no more digits after dot
1484    $MBI->_lsft($ym, $x->{_e}, 10);	# 123 => 1230
1485    }
1486  # 123e1 % 20 => 1230 % 20
1487  if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
1488    {
1489    $MBI->_lsft( $x->{_m}, $x->{_e},10);	# es => '+' here
1490    }
1491
1492  $x->{_e} = $MBI->_new($shiftx);
1493  $x->{_es} = '+';
1494  $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1495  $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
1496
1497  # now mantissas are equalized, exponent of $x is adjusted, so calc result
1498
1499  $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
1500
1501  $x->{sign} = '+' if $MBI->_is_zero($x->{_m});		# fix sign for -0
1502  $x->bnorm();
1503
1504  if ($neg != 0)	# one of them negative => correct in place
1505    {
1506    my $r = $y - $x;
1507    $x->{_m} = $r->{_m};
1508    $x->{_e} = $r->{_e};
1509    $x->{_es} = $r->{_es};
1510    $x->{sign} = '+' if $MBI->_is_zero($x->{_m});	# fix sign for -0
1511    $x->bnorm();
1512    }
1513
1514  $x->round($a,$p,$r,$y);	# round and return
1515  }
1516
1517sub broot
1518  {
1519  # calculate $y'th root of $x
1520
1521  # set up parameters
1522  my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1523  # objectify is costly, so avoid it
1524  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1525    {
1526    ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1527    }
1528
1529  # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
1530  return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
1531         $y->{sign} !~ /^\+$/;
1532
1533  return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
1534
1535  # we need to limit the accuracy to protect against overflow
1536  my $fallback = 0;
1537  my (@params,$scale);
1538  ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1539
1540  return $x if $x->is_nan();		# error in _find_round_parameters?
1541
1542  # no rounding at all, so must use fallback
1543  if (scalar @params == 0)
1544    {
1545    # simulate old behaviour
1546    $params[0] = $self->div_scale();	# and round to it as accuracy
1547    $scale = $params[0]+4; 		# at least four more for proper round
1548    $params[2] = $r;			# iound mode by caller or undef
1549    $fallback = 1;			# to clear a/p afterwards
1550    }
1551  else
1552    {
1553    # the 4 below is empirical, and there might be cases where it is not
1554    # enough...
1555    $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1556    }
1557
1558  # when user set globals, they would interfere with our calculation, so
1559  # disable them and later re-enable them
1560  no strict 'refs';
1561  my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1562  my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1563  # we also need to disable any set A or P on $x (_find_round_parameters took
1564  # them already into account), since these would interfere, too
1565  delete $x->{_a}; delete $x->{_p};
1566  # need to disable $upgrade in BigInt, to avoid deep recursion
1567  local $Math::BigInt::upgrade = undef;	# should be really parent class vs MBI
1568
1569  # remember sign and make $x positive, since -4 ** (1/2) => -2
1570  my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
1571
1572  my $is_two = 0;
1573  if ($y->isa('Math::BigFloat'))
1574    {
1575    $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
1576    }
1577  else
1578    {
1579    $is_two = ($y == 2);
1580    }
1581
1582  # normal square root if $y == 2:
1583  if ($is_two)
1584    {
1585    $x->bsqrt($scale+4);
1586    }
1587  elsif ($y->is_one('-'))
1588    {
1589    # $x ** -1 => 1/$x
1590    my $u = $self->bone()->bdiv($x,$scale);
1591    # copy private parts over
1592    $x->{_m} = $u->{_m};
1593    $x->{_e} = $u->{_e};
1594    $x->{_es} = $u->{_es};
1595    }
1596  else
1597    {
1598    # calculate the broot() as integer result first, and if it fits, return
1599    # it rightaway (but only if $x and $y are integer):
1600
1601    my $done = 0;				# not yet
1602    if ($y->is_int() && $x->is_int())
1603      {
1604      my $i = $MBI->_copy( $x->{_m} );
1605      $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1606      my $int = Math::BigInt->bzero();
1607      $int->{value} = $i;
1608      $int->broot($y->as_number());
1609      # if ($exact)
1610      if ($int->copy()->bpow($y) == $x)
1611        {
1612        # found result, return it
1613        $x->{_m} = $int->{value};
1614        $x->{_e} = $MBI->_zero();
1615        $x->{_es} = '+';
1616        $x->bnorm();
1617        $done = 1;
1618        }
1619      }
1620    if ($done == 0)
1621      {
1622      my $u = $self->bone()->bdiv($y,$scale+4);
1623      delete $u->{_a}; delete $u->{_p};         # otherwise it conflicts
1624      $x->bpow($u,$scale+4);                    # el cheapo
1625      }
1626    }
1627  $x->bneg() if $sign == 1;
1628
1629  # shortcut to not run through _find_round_parameters again
1630  if (defined $params[0])
1631    {
1632    $x->bround($params[0],$params[2]);		# then round accordingly
1633    }
1634  else
1635    {
1636    $x->bfround($params[1],$params[2]);		# then round accordingly
1637    }
1638  if ($fallback)
1639    {
1640    # clear a/p after round, since user did not request it
1641    delete $x->{_a}; delete $x->{_p};
1642    }
1643  # restore globals
1644  $$abr = $ab; $$pbr = $pb;
1645  $x;
1646  }
1647
1648sub bsqrt
1649  {
1650  # calculate square root
1651  my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1652
1653  return $x->bnan() if $x->{sign} !~ /^[+]/;	# NaN, -inf or < 0
1654  return $x if $x->{sign} eq '+inf';		# sqrt(inf) == inf
1655  return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
1656
1657  # we need to limit the accuracy to protect against overflow
1658  my $fallback = 0;
1659  my (@params,$scale);
1660  ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1661
1662  return $x if $x->is_nan();		# error in _find_round_parameters?
1663
1664  # no rounding at all, so must use fallback
1665  if (scalar @params == 0)
1666    {
1667    # simulate old behaviour
1668    $params[0] = $self->div_scale();	# and round to it as accuracy
1669    $scale = $params[0]+4; 		# at least four more for proper round
1670    $params[2] = $r;			# round mode by caller or undef
1671    $fallback = 1;			# to clear a/p afterwards
1672    }
1673  else
1674    {
1675    # the 4 below is empirical, and there might be cases where it is not
1676    # enough...
1677    $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1678    }
1679
1680  # when user set globals, they would interfere with our calculation, so
1681  # disable them and later re-enable them
1682  no strict 'refs';
1683  my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1684  my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1685  # we also need to disable any set A or P on $x (_find_round_parameters took
1686  # them already into account), since these would interfere, too
1687  delete $x->{_a}; delete $x->{_p};
1688  # need to disable $upgrade in BigInt, to avoid deep recursion
1689  local $Math::BigInt::upgrade = undef;	# should be really parent class vs MBI
1690
1691  my $i = $MBI->_copy( $x->{_m} );
1692  $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
1693  my $xas = Math::BigInt->bzero();
1694  $xas->{value} = $i;
1695
1696  my $gs = $xas->copy()->bsqrt();	# some guess
1697
1698  if (($x->{_es} ne '-')		# guess can't be accurate if there are
1699					# digits after the dot
1700   && ($xas->bacmp($gs * $gs) == 0))	# guess hit the nail on the head?
1701    {
1702    # exact result, copy result over to keep $x
1703    $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
1704    $x->bnorm();
1705    # shortcut to not run through _find_round_parameters again
1706    if (defined $params[0])
1707      {
1708      $x->bround($params[0],$params[2]);	# then round accordingly
1709      }
1710    else
1711      {
1712      $x->bfround($params[1],$params[2]);	# then round accordingly
1713      }
1714    if ($fallback)
1715      {
1716      # clear a/p after round, since user did not request it
1717      delete $x->{_a}; delete $x->{_p};
1718      }
1719    # re-enable A and P, upgrade is taken care of by "local"
1720    ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1721    return $x;
1722    }
1723
1724  # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
1725  # of the result by multipyling the input by 100 and then divide the integer
1726  # result of sqrt(input) by 10. Rounding afterwards returns the real result.
1727
1728  # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
1729  my $y1 = $MBI->_copy($x->{_m});
1730
1731  my $length = $MBI->_len($y1);
1732
1733  # Now calculate how many digits the result of sqrt(y1) would have
1734  my $digits = int($length / 2);
1735
1736  # But we need at least $scale digits, so calculate how many are missing
1737  my $shift = $scale - $digits;
1738
1739  # That should never happen (we take care of integer guesses above)
1740  # $shift = 0 if $shift < 0;
1741
1742  # Multiply in steps of 100, by shifting left two times the "missing" digits
1743  my $s2 = $shift * 2;
1744
1745  # We now make sure that $y1 has the same odd or even number of digits than
1746  # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
1747  # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
1748  # steps of 10. The length of $x does not count, since an even or odd number
1749  # of digits before the dot is not changed by adding an even number of digits
1750  # after the dot (the result is still odd or even digits long).
1751  $s2++ if $MBI->_is_odd($x->{_e});
1752
1753  $MBI->_lsft( $y1, $MBI->_new($s2), 10);
1754
1755  # now take the square root and truncate to integer
1756  $y1 = $MBI->_sqrt($y1);
1757
1758  # By "shifting" $y1 right (by creating a negative _e) we calculate the final
1759  # result, which is than later rounded to the desired scale.
1760
1761  # calculate how many zeros $x had after the '.' (or before it, depending
1762  # on sign of $dat, the result should have half as many:
1763  my $dat = $MBI->_num($x->{_e});
1764  $dat = -$dat if $x->{_es} eq '-';
1765  $dat += $length;
1766
1767  if ($dat > 0)
1768    {
1769    # no zeros after the dot (e.g. 1.23, 0.49 etc)
1770    # preserve half as many digits before the dot than the input had
1771    # (but round this "up")
1772    $dat = int(($dat+1)/2);
1773    }
1774  else
1775    {
1776    $dat = int(($dat)/2);
1777    }
1778  $dat -= $MBI->_len($y1);
1779  if ($dat < 0)
1780    {
1781    $dat = abs($dat);
1782    $x->{_e} = $MBI->_new( $dat );
1783    $x->{_es} = '-';
1784    }
1785  else
1786    {
1787    $x->{_e} = $MBI->_new( $dat );
1788    $x->{_es} = '+';
1789    }
1790  $x->{_m} = $y1;
1791  $x->bnorm();
1792
1793  # shortcut to not run through _find_round_parameters again
1794  if (defined $params[0])
1795    {
1796    $x->bround($params[0],$params[2]);		# then round accordingly
1797    }
1798  else
1799    {
1800    $x->bfround($params[1],$params[2]);		# then round accordingly
1801    }
1802  if ($fallback)
1803    {
1804    # clear a/p after round, since user did not request it
1805    delete $x->{_a}; delete $x->{_p};
1806    }
1807  # restore globals
1808  $$abr = $ab; $$pbr = $pb;
1809  $x;
1810  }
1811
1812sub bfac
1813  {
1814  # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1815  # compute factorial number, modifies first argument
1816
1817  # set up parameters
1818  my ($self,$x,@r) = (ref($_[0]),@_);
1819  # objectify is costly, so avoid it
1820  ($self,$x,@r) = objectify(1,@_) if !ref($x);
1821
1822 return $x if $x->{sign} eq '+inf';	# inf => inf
1823  return $x->bnan()
1824    if (($x->{sign} ne '+') ||		# inf, NaN, <0 etc => NaN
1825     ($x->{_es} ne '+'));		# digits after dot?
1826
1827  # use BigInt's bfac() for faster calc
1828  if (! $MBI->_is_zero($x->{_e}))
1829    {
1830    $MBI->_lsft($x->{_m}, $x->{_e},10);	# change 12e1 to 120e0
1831    $x->{_e} = $MBI->_zero();		# normalize
1832    $x->{_es} = '+';
1833    }
1834  $MBI->_fac($x->{_m});			# calculate factorial
1835  $x->bnorm()->round(@r); 		# norm again and round result
1836  }
1837
1838sub _pow
1839  {
1840  # Calculate a power where $y is a non-integer, like 2 ** 0.5
1841  my ($x,$y,$a,$p,$r) = @_;
1842  my $self = ref($x);
1843
1844  # if $y == 0.5, it is sqrt($x)
1845  $HALF = $self->new($HALF) unless ref($HALF);
1846  return $x->bsqrt($a,$p,$r,$y) if $y->bcmp($HALF) == 0;
1847
1848  # Using:
1849  # a ** x == e ** (x * ln a)
1850
1851  # u = y * ln x
1852  #                _                         _
1853  # Taylor:       |   u    u^2    u^3         |
1854  # x ** y  = 1 + |  --- + --- + ----- + ...  |
1855  #               |_  1    1*2   1*2*3       _|
1856
1857  # we need to limit the accuracy to protect against overflow
1858  my $fallback = 0;
1859  my ($scale,@params);
1860  ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1861
1862  return $x if $x->is_nan();		# error in _find_round_parameters?
1863
1864  # no rounding at all, so must use fallback
1865  if (scalar @params == 0)
1866    {
1867    # simulate old behaviour
1868    $params[0] = $self->div_scale();	# and round to it as accuracy
1869    $params[1] = undef;			# disable P
1870    $scale = $params[0]+4; 		# at least four more for proper round
1871    $params[2] = $r;			# round mode by caller or undef
1872    $fallback = 1;			# to clear a/p afterwards
1873    }
1874  else
1875    {
1876    # the 4 below is empirical, and there might be cases where it is not
1877    # enough...
1878    $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1879    }
1880
1881  # when user set globals, they would interfere with our calculation, so
1882  # disable them and later re-enable them
1883  no strict 'refs';
1884  my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1885  my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1886  # we also need to disable any set A or P on $x (_find_round_parameters took
1887  # them already into account), since these would interfere, too
1888  delete $x->{_a}; delete $x->{_p};
1889  # need to disable $upgrade in BigInt, to avoid deep recursion
1890  local $Math::BigInt::upgrade = undef;
1891
1892  my ($limit,$v,$u,$below,$factor,$next,$over);
1893
1894  $u = $x->copy()->blog(undef,$scale)->bmul($y);
1895  $v = $self->bone();				# 1
1896  $factor = $self->new(2);			# 2
1897  $x->bone();					# first term: 1
1898
1899  $below = $v->copy();
1900  $over = $u->copy();
1901
1902  $limit = $self->new("1E-". ($scale-1));
1903  #my $steps = 0;
1904  while (3 < 5)
1905    {
1906    # we calculate the next term, and add it to the last
1907    # when the next term is below our limit, it won't affect the outcome
1908    # anymore, so we stop
1909    $next = $over->copy()->bdiv($below,$scale);
1910    last if $next->bacmp($limit) <= 0;
1911    $x->badd($next);
1912    # calculate things for the next term
1913    $over *= $u; $below *= $factor; $factor->binc();
1914
1915    last if $x->{sign} !~ /^[-+]$/;
1916
1917    #$steps++;
1918    }
1919
1920  # shortcut to not run through _find_round_parameters again
1921  if (defined $params[0])
1922    {
1923    $x->bround($params[0],$params[2]);		# then round accordingly
1924    }
1925  else
1926    {
1927    $x->bfround($params[1],$params[2]);		# then round accordingly
1928    }
1929  if ($fallback)
1930    {
1931    # clear a/p after round, since user did not request it
1932    delete $x->{_a}; delete $x->{_p};
1933    }
1934  # restore globals
1935  $$abr = $ab; $$pbr = $pb;
1936  $x;
1937  }
1938
1939sub bpow
1940  {
1941  # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1942  # compute power of two numbers, second arg is used as integer
1943  # modifies first argument
1944
1945  # set up parameters
1946  my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1947  # objectify is costly, so avoid it
1948  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1949    {
1950    ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1951    }
1952
1953  return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1954  return $x if $x->{sign} =~ /^[+-]inf$/;
1955
1956  # -2 ** -2 => NaN
1957  return $x->bnan() if $x->{sign} eq '-' && $y->{sign} eq '-';
1958
1959  # cache the result of is_zero
1960  my $y_is_zero = $y->is_zero();
1961  return $x->bone() if $y_is_zero;
1962  return $x         if $x->is_one() || $y->is_one();
1963
1964  my $x_is_zero = $x->is_zero();
1965  return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int();		# non-integer power
1966
1967  my $y1 = $y->as_number()->{value};			# make MBI part
1968
1969  # if ($x == -1)
1970  if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
1971    {
1972    # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
1973    return $MBI->_is_odd($y1) ? $x : $x->babs(1);
1974    }
1975  if ($x_is_zero)
1976    {
1977    return $x->bone() if $y_is_zero;
1978    return $x if $y->{sign} eq '+'; 	# 0**y => 0 (if not y <= 0)
1979    # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
1980    return $x->binf();
1981    }
1982
1983  my $new_sign = '+';
1984  $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
1985
1986  # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1987  $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
1988  $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
1989
1990  $x->{sign} = $new_sign;
1991  $x->bnorm();
1992  if ($y->{sign} eq '-')
1993    {
1994    # modify $x in place!
1995    my $z = $x->copy(); $x->bone();
1996    return $x->bdiv($z,$a,$p,$r);	# round in one go (might ignore y's A!)
1997    }
1998  $x->round($a,$p,$r,$y);
1999  }
2000
2001###############################################################################
2002# rounding functions
2003
2004sub bfround
2005  {
2006  # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
2007  # $n == 0 means round to integer
2008  # expects and returns normalized numbers!
2009  my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2010
2011  my ($scale,$mode) = $x->_scale_p(@_);
2012  return $x if !defined $scale || $x->modify('bfround'); # no-op
2013
2014  # never round a 0, +-inf, NaN
2015  if ($x->is_zero())
2016    {
2017    $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
2018    return $x;
2019    }
2020  return $x if $x->{sign} !~ /^[+-]$/;
2021
2022  # don't round if x already has lower precision
2023  return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
2024
2025  $x->{_p} = $scale;			# remember round in any case
2026  delete $x->{_a};			# and clear A
2027  if ($scale < 0)
2028    {
2029    # round right from the '.'
2030
2031    return $x if $x->{_es} eq '+';		# e >= 0 => nothing to round
2032
2033    $scale = -$scale;				# positive for simplicity
2034    my $len = $MBI->_len($x->{_m});		# length of mantissa
2035
2036    # the following poses a restriction on _e, but if _e is bigger than a
2037    # scalar, you got other problems (memory etc) anyway
2038    my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e})));	# digits after dot
2039    my $zad = 0;				# zeros after dot
2040    $zad = $dad - $len if (-$dad < -$len);	# for 0.00..00xxx style
2041
2042    # p rint "scale $scale dad $dad zad $zad len $len\n";
2043    # number  bsstr   len zad dad
2044    # 0.123   123e-3	3   0 3
2045    # 0.0123  123e-4	3   1 4
2046    # 0.001   1e-3      1   2 3
2047    # 1.23    123e-2	3   0 2
2048    # 1.2345  12345e-4	5   0 4
2049
2050    # do not round after/right of the $dad
2051    return $x if $scale > $dad;			# 0.123, scale >= 3 => exit
2052
2053    # round to zero if rounding inside the $zad, but not for last zero like:
2054    # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
2055    return $x->bzero() if $scale < $zad;
2056    if ($scale == $zad)			# for 0.006, scale -3 and trunc
2057      {
2058      $scale = -$len;
2059      }
2060    else
2061      {
2062      # adjust round-point to be inside mantissa
2063      if ($zad != 0)
2064        {
2065	$scale = $scale-$zad;
2066        }
2067      else
2068        {
2069        my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;	# digits before dot
2070	$scale = $dbd+$scale;
2071        }
2072      }
2073    }
2074  else
2075    {
2076    # round left from the '.'
2077
2078    # 123 => 100 means length(123) = 3 - $scale (2) => 1
2079
2080    my $dbt = $MBI->_len($x->{_m});
2081    # digits before dot
2082    my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
2083    # should be the same, so treat it as this
2084    $scale = 1 if $scale == 0;
2085    # shortcut if already integer
2086    return $x if $scale == 1 && $dbt <= $dbd;
2087    # maximum digits before dot
2088    ++$dbd;
2089
2090    if ($scale > $dbd)
2091       {
2092       # not enough digits before dot, so round to zero
2093       return $x->bzero;
2094       }
2095    elsif ( $scale == $dbd )
2096       {
2097       # maximum
2098       $scale = -$dbt;
2099       }
2100    else
2101       {
2102       $scale = $dbd - $scale;
2103       }
2104    }
2105  # pass sign to bround for rounding modes '+inf' and '-inf'
2106  my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2107  $m->bround($scale,$mode);
2108  $x->{_m} = $m->{value};			# get our mantissa back
2109  $x->bnorm();
2110  }
2111
2112sub bround
2113  {
2114  # accuracy: preserve $N digits, and overwrite the rest with 0's
2115  my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
2116
2117  if (($_[0] || 0) < 0)
2118    {
2119    require Carp; Carp::croak ('bround() needs positive accuracy');
2120    }
2121
2122  my ($scale,$mode) = $x->_scale_a(@_);
2123  return $x if !defined $scale || $x->modify('bround');	# no-op
2124
2125  # scale is now either $x->{_a}, $accuracy, or the user parameter
2126  # test whether $x already has lower accuracy, do nothing in this case
2127  # but do round if the accuracy is the same, since a math operation might
2128  # want to round a number with A=5 to 5 digits afterwards again
2129  return $x if defined $x->{_a} && $x->{_a} < $scale;
2130
2131  # scale < 0 makes no sense
2132  # scale == 0 => keep all digits
2133  # never round a +-inf, NaN
2134  return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
2135
2136  # 1: never round a 0
2137  # 2: if we should keep more digits than the mantissa has, do nothing
2138  if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
2139    {
2140    $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
2141    return $x;
2142    }
2143
2144  # pass sign to bround for '+inf' and '-inf' rounding modes
2145  my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
2146
2147  $m->bround($scale,$mode);		# round mantissa
2148  $x->{_m} = $m->{value};		# get our mantissa back
2149  $x->{_a} = $scale;			# remember rounding
2150  delete $x->{_p};			# and clear P
2151  $x->bnorm();				# del trailing zeros gen. by bround()
2152  }
2153
2154sub bfloor
2155  {
2156  # return integer less or equal then $x
2157  my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2158
2159  return $x if $x->modify('bfloor');
2160
2161  return $x if $x->{sign} !~ /^[+-]$/;	# nan, +inf, -inf
2162
2163  # if $x has digits after dot
2164  if ($x->{_es} eq '-')
2165    {
2166    $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2167    $x->{_e} = $MBI->_zero();			# trunc/norm
2168    $x->{_es} = '+';				# abs e
2169    $MBI->_inc($x->{_m}) if $x->{sign} eq '-';	# increment if negative
2170    }
2171  $x->round($a,$p,$r);
2172  }
2173
2174sub bceil
2175  {
2176  # return integer greater or equal then $x
2177  my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2178
2179  return $x if $x->modify('bceil');
2180  return $x if $x->{sign} !~ /^[+-]$/;	# nan, +inf, -inf
2181
2182  # if $x has digits after dot
2183  if ($x->{_es} eq '-')
2184    {
2185    $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
2186    $x->{_e} = $MBI->_zero();			# trunc/norm
2187    $x->{_es} = '+';				# abs e
2188    $MBI->_inc($x->{_m}) if $x->{sign} eq '+';	# increment if positive
2189    }
2190  $x->round($a,$p,$r);
2191  }
2192
2193sub brsft
2194  {
2195  # shift right by $y (divide by power of $n)
2196
2197  # set up parameters
2198  my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2199  # objectify is costly, so avoid it
2200  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2201    {
2202    ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2203    }
2204
2205  return $x if $x->modify('brsft');
2206  return $x if $x->{sign} !~ /^[+-]$/;	# nan, +inf, -inf
2207
2208  $n = 2 if !defined $n; $n = $self->new($n);
2209  $x->bdiv($n->bpow($y),$a,$p,$r,$y);
2210  }
2211
2212sub blsft
2213  {
2214  # shift left by $y (multiply by power of $n)
2215
2216  # set up parameters
2217  my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
2218  # objectify is costly, so avoid it
2219  if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2220    {
2221    ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
2222    }
2223
2224  return $x if $x->modify('blsft');
2225  return $x if $x->{sign} !~ /^[+-]$/;	# nan, +inf, -inf
2226
2227  $n = 2 if !defined $n; $n = $self->new($n);
2228  $x->bmul($n->bpow($y),$a,$p,$r,$y);
2229  }
2230
2231###############################################################################
2232
2233sub DESTROY
2234  {
2235  # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
2236  }
2237
2238sub AUTOLOAD
2239  {
2240  # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
2241  # or falling back to MBI::bxxx()
2242  my $name = $AUTOLOAD;
2243
2244  $name =~ s/(.*):://;	# split package
2245  my $c = $1 || $class;
2246  no strict 'refs';
2247  $c->import() if $IMPORT == 0;
2248  if (!method_alias($name))
2249    {
2250    if (!defined $name)
2251      {
2252      # delayed load of Carp and avoid recursion
2253      require Carp;
2254      Carp::croak ("$c: Can't call a method without name");
2255      }
2256    if (!method_hand_up($name))
2257      {
2258      # delayed load of Carp and avoid recursion
2259      require Carp;
2260      Carp::croak ("Can't call $c\-\>$name, not a valid method");
2261      }
2262    # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
2263    $name =~ s/^f/b/;
2264    return &{"Math::BigInt"."::$name"}(@_);
2265    }
2266  my $bname = $name; $bname =~ s/^f/b/;
2267  $c .= "::$name";
2268  *{$c} = \&{$bname};
2269  &{$c};	# uses @_
2270  }
2271
2272sub exponent
2273  {
2274  # return a copy of the exponent
2275  my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2276
2277  if ($x->{sign} !~ /^[+-]$/)
2278    {
2279    my $s = $x->{sign}; $s =~ s/^[+-]//;
2280    return Math::BigInt->new($s); 		# -inf, +inf => +inf
2281    }
2282  Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
2283  }
2284
2285sub mantissa
2286  {
2287  # return a copy of the mantissa
2288  my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2289
2290  if ($x->{sign} !~ /^[+-]$/)
2291    {
2292    my $s = $x->{sign}; $s =~ s/^[+]//;
2293    return Math::BigInt->new($s);		# -inf, +inf => +inf
2294    }
2295  my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
2296  $m->bneg() if $x->{sign} eq '-';
2297
2298  $m;
2299  }
2300
2301sub parts
2302  {
2303  # return a copy of both the exponent and the mantissa
2304  my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2305
2306  if ($x->{sign} !~ /^[+-]$/)
2307    {
2308    my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2309    return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2310    }
2311  my $m = Math::BigInt->bzero();
2312  $m->{value} = $MBI->_copy($x->{_m});
2313  $m->bneg() if $x->{sign} eq '-';
2314  ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
2315  }
2316
2317##############################################################################
2318# private stuff (internal use only)
2319
2320sub import
2321  {
2322  my $self = shift;
2323  my $l = scalar @_;
2324  my $lib = ''; my @a;
2325  $IMPORT=1;
2326  for ( my $i = 0; $i < $l ; $i++)
2327    {
2328    if ( $_[$i] eq ':constant' )
2329      {
2330      # This causes overlord er load to step in. 'binary' and 'integer'
2331      # are handled by BigInt.
2332      overload::constant float => sub { $self->new(shift); };
2333      }
2334    elsif ($_[$i] eq 'upgrade')
2335      {
2336      # this causes upgrading
2337      $upgrade = $_[$i+1];		# or undef to disable
2338      $i++;
2339      }
2340    elsif ($_[$i] eq 'downgrade')
2341      {
2342      # this causes downgrading
2343      $downgrade = $_[$i+1];		# or undef to disable
2344      $i++;
2345      }
2346    elsif ($_[$i] eq 'lib')
2347      {
2348      # alternative library
2349      $lib = $_[$i+1] || '';		# default Calc
2350      $i++;
2351      }
2352    elsif ($_[$i] eq 'with')
2353      {
2354      # alternative class for our private parts()
2355      # XXX: no longer supported
2356      # $MBI = $_[$i+1] || 'Math::BigInt';
2357      $i++;
2358      }
2359    else
2360      {
2361      push @a, $_[$i];
2362      }
2363    }
2364
2365  $lib =~ tr/a-zA-Z0-9,://cd;		# restrict to sane characters
2366  # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2367  my $mbilib = eval { Math::BigInt->config()->{lib} };
2368  if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
2369    {
2370    # MBI already loaded
2371    Math::BigInt->import('lib',"$lib,$mbilib", 'objectify');
2372    }
2373  else
2374    {
2375    # MBI not loaded, or with ne "Math::BigInt::Calc"
2376    $lib .= ",$mbilib" if defined $mbilib;
2377    $lib =~ s/^,//;				# don't leave empty
2378
2379    # replacement library can handle lib statement, but also could ignore it
2380
2381    # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2382    # used in the same script, or eval inside import(). So we require MBI:
2383    require Math::BigInt;
2384    Math::BigInt->import( lib => $lib, 'objectify' );
2385    }
2386  if ($@)
2387    {
2388    require Carp; Carp::croak ("Couldn't load $lib: $! $@");
2389    }
2390  # find out which one was actually loaded
2391  $MBI = Math::BigInt->config()->{lib};
2392
2393  # register us with MBI to get notified of future lib changes
2394  Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
2395
2396  # any non :constant stuff is handled by our parent, Exporter
2397  # even if @_ is empty, to give it a chance
2398  $self->SUPER::import(@a);      	# for subclasses
2399  $self->export_to_level(1,$self,@a);	# need this, too
2400  }
2401
2402sub bnorm
2403  {
2404  # adjust m and e so that m is smallest possible
2405  # round number according to accuracy and precision settings
2406  my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
2407
2408  return $x if $x->{sign} !~ /^[+-]$/;		# inf, nan etc
2409
2410  my $zeros = $MBI->_zeros($x->{_m});		# correct for trailing zeros
2411  if ($zeros != 0)
2412    {
2413    my $z = $MBI->_new($zeros);
2414    $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
2415    if ($x->{_es} eq '-')
2416      {
2417      if ($MBI->_acmp($x->{_e},$z) >= 0)
2418        {
2419        $x->{_e} = $MBI->_sub  ($x->{_e}, $z);
2420        $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
2421        }
2422      else
2423        {
2424        $x->{_e} = $MBI->_sub  ( $MBI->_copy($z), $x->{_e});
2425        $x->{_es} = '+';
2426        }
2427      }
2428    else
2429      {
2430      $x->{_e} = $MBI->_add  ($x->{_e}, $z);
2431      }
2432    }
2433  else
2434    {
2435    # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
2436    # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
2437    $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
2438     if $MBI->_is_zero($x->{_m});
2439    }
2440
2441  $x;					# MBI bnorm is no-op, so dont call it
2442  }
2443
2444##############################################################################
2445
2446sub as_hex
2447  {
2448  # return number as hexadecimal string (only for integers defined)
2449  my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2450
2451  return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
2452  return '0x0' if $x->is_zero();
2453
2454  return $nan if $x->{_es} ne '+';		# how to do 1e-1 in hex!?
2455
2456  my $z = $MBI->_copy($x->{_m});
2457  if (! $MBI->_is_zero($x->{_e}))		# > 0
2458    {
2459    $MBI->_lsft( $z, $x->{_e},10);
2460    }
2461  $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2462  $z->as_hex();
2463  }
2464
2465sub as_bin
2466  {
2467  # return number as binary digit string (only for integers defined)
2468  my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2469
2470  return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
2471  return '0b0' if $x->is_zero();
2472
2473  return $nan if $x->{_es} ne '+';		# how to do 1e-1 in hex!?
2474
2475  my $z = $MBI->_copy($x->{_m});
2476  if (! $MBI->_is_zero($x->{_e}))		# > 0
2477    {
2478    $MBI->_lsft( $z, $x->{_e},10);
2479    }
2480  $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2481  $z->as_bin();
2482  }
2483
2484sub as_number
2485  {
2486  # return copy as a bigint representation of this BigFloat number
2487  my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2488
2489  my $z = $MBI->_copy($x->{_m});
2490  if ($x->{_es} eq '-')			# < 0
2491    {
2492    $MBI->_rsft( $z, $x->{_e},10);
2493    }
2494  elsif (! $MBI->_is_zero($x->{_e}))	# > 0
2495    {
2496    $MBI->_lsft( $z, $x->{_e},10);
2497    }
2498  $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
2499  $z;
2500  }
2501
2502sub length
2503  {
2504  my $x = shift;
2505  my $class = ref($x) || $x;
2506  $x = $class->new(shift) unless ref($x);
2507
2508  return 1 if $MBI->_is_zero($x->{_m});
2509
2510  my $len = $MBI->_len($x->{_m});
2511  $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
2512  if (wantarray())
2513    {
2514    my $t = 0;
2515    $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
2516    return ($len, $t);
2517    }
2518  $len;
2519  }
2520
25211;
2522__END__
2523
2524=head1 NAME
2525
2526Math::BigFloat - Arbitrary size floating point math package
2527
2528=head1 SYNOPSIS
2529
2530  use Math::BigFloat;
2531
2532  # Number creation
2533  $x = Math::BigFloat->new($str);	# defaults to 0
2534  $nan  = Math::BigFloat->bnan();	# create a NotANumber
2535  $zero = Math::BigFloat->bzero();	# create a +0
2536  $inf = Math::BigFloat->binf();	# create a +inf
2537  $inf = Math::BigFloat->binf('-');	# create a -inf
2538  $one = Math::BigFloat->bone();	# create a +1
2539  $one = Math::BigFloat->bone('-');	# create a -1
2540
2541  # Testing
2542  $x->is_zero();		# true if arg is +0
2543  $x->is_nan();			# true if arg is NaN
2544  $x->is_one();			# true if arg is +1
2545  $x->is_one('-');		# true if arg is -1
2546  $x->is_odd();			# true if odd, false for even
2547  $x->is_even();		# true if even, false for odd
2548  $x->is_pos();			# true if >= 0
2549  $x->is_neg();			# true if <  0
2550  $x->is_inf(sign);		# true if +inf, or -inf (default is '+')
2551
2552  $x->bcmp($y);			# compare numbers (undef,<0,=0,>0)
2553  $x->bacmp($y);		# compare absolutely (undef,<0,=0,>0)
2554  $x->sign();			# return the sign, either +,- or NaN
2555  $x->digit($n);		# return the nth digit, counting from right
2556  $x->digit(-$n);		# return the nth digit, counting from left
2557
2558  # The following all modify their first argument. If you want to preserve
2559  # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2560  # neccessary when mixing $a = $b assigments with non-overloaded math.
2561
2562  # set
2563  $x->bzero();			# set $i to 0
2564  $x->bnan();			# set $i to NaN
2565  $x->bone();                   # set $x to +1
2566  $x->bone('-');                # set $x to -1
2567  $x->binf();                   # set $x to inf
2568  $x->binf('-');                # set $x to -inf
2569
2570  $x->bneg();			# negation
2571  $x->babs();			# absolute value
2572  $x->bnorm();			# normalize (no-op)
2573  $x->bnot();			# two's complement (bit wise not)
2574  $x->binc();			# increment x by 1
2575  $x->bdec();			# decrement x by 1
2576
2577  $x->badd($y);			# addition (add $y to $x)
2578  $x->bsub($y);			# subtraction (subtract $y from $x)
2579  $x->bmul($y);			# multiplication (multiply $x by $y)
2580  $x->bdiv($y);			# divide, set $x to quotient
2581				# return (quo,rem) or quo if scalar
2582
2583  $x->bmod($y);			# modulus ($x % $y)
2584  $x->bpow($y);			# power of arguments ($x ** $y)
2585  $x->blsft($y);		# left shift
2586  $x->brsft($y);		# right shift
2587				# return (quo,rem) or quo if scalar
2588
2589  $x->blog();			# logarithm of $x to base e (Euler's number)
2590  $x->blog($base);		# logarithm of $x to base $base (f.i. 2)
2591
2592  $x->band($y);			# bit-wise and
2593  $x->bior($y);			# bit-wise inclusive or
2594  $x->bxor($y);			# bit-wise exclusive or
2595  $x->bnot();			# bit-wise not (two's complement)
2596
2597  $x->bsqrt();			# calculate square-root
2598  $x->broot($y);		# $y'th root of $x (e.g. $y == 3 => cubic root)
2599  $x->bfac();			# factorial of $x (1*2*3*4*..$x)
2600
2601  $x->bround($N); 		# accuracy: preserve $N digits
2602  $x->bfround($N);		# precision: round to the $Nth digit
2603
2604  $x->bfloor();			# return integer less or equal than $x
2605  $x->bceil();			# return integer greater or equal than $x
2606
2607  # The following do not modify their arguments:
2608
2609  bgcd(@values);		# greatest common divisor
2610  blcm(@values);		# lowest common multiplicator
2611
2612  $x->bstr();			# return string
2613  $x->bsstr();			# return string in scientific notation
2614
2615  $x->as_int();			# return $x as BigInt
2616  $x->exponent();		# return exponent as BigInt
2617  $x->mantissa();		# return mantissa as BigInt
2618  $x->parts();			# return (mantissa,exponent) as BigInt
2619
2620  $x->length();			# number of digits (w/o sign and '.')
2621  ($l,$f) = $x->length();	# number of digits, and length of fraction
2622
2623  $x->precision();		# return P of $x (or global, if P of $x undef)
2624  $x->precision($n);		# set P of $x to $n
2625  $x->accuracy();		# return A of $x (or global, if A of $x undef)
2626  $x->accuracy($n);		# set A $x to $n
2627
2628  # these get/set the appropriate global value for all BigFloat objects
2629  Math::BigFloat->precision();	# Precision
2630  Math::BigFloat->accuracy();	# Accuracy
2631  Math::BigFloat->round_mode();	# rounding mode
2632
2633=head1 DESCRIPTION
2634
2635All operators (inlcuding basic math operations) are overloaded if you
2636declare your big floating point numbers as
2637
2638  $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2639
2640Operations with overloaded operators preserve the arguments, which is
2641exactly what you expect.
2642
2643=head2 Canonical notation
2644
2645Input to these routines are either BigFloat objects, or strings of the
2646following four forms:
2647
2648=over 2
2649
2650=item *
2651
2652C</^[+-]\d+$/>
2653
2654=item *
2655
2656C</^[+-]\d+\.\d*$/>
2657
2658=item *
2659
2660C</^[+-]\d+E[+-]?\d+$/>
2661
2662=item *
2663
2664C</^[+-]\d*\.\d+E[+-]?\d+$/>
2665
2666=back
2667
2668all with optional leading and trailing zeros and/or spaces. Additonally,
2669numbers are allowed to have an underscore between any two digits.
2670
2671Empty strings as well as other illegal numbers results in 'NaN'.
2672
2673bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2674are always stored in normalized form. On a string, it creates a BigFloat
2675object.
2676
2677=head2 Output
2678
2679Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2680
2681The string output will always have leading and trailing zeros stripped and drop
2682a plus sign. C<bstr()> will give you always the form with a decimal point,
2683while C<bsstr()> (s for scientific) gives you the scientific notation.
2684
2685	Input			bstr()		bsstr()
2686	'-0'			'0'		'0E1'
2687   	'  -123 123 123'	'-123123123'	'-123123123E0'
2688	'00.0123'		'0.0123'	'123E-4'
2689	'123.45E-2'		'1.2345'	'12345E-4'
2690	'10E+3'			'10000'		'1E4'
2691
2692Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2693C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2694return either undef, <0, 0 or >0 and are suited for sort.
2695
2696Actual math is done by using the class defined with C<with => Class;> (which
2697defaults to BigInts) to represent the mantissa and exponent.
2698
2699The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2700represent the result when input arguments are not numbers, as well as
2701the result of dividing by zero.
2702
2703=head2 C<mantissa()>, C<exponent()> and C<parts()>
2704
2705C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2706as BigInts such that:
2707
2708	$m = $x->mantissa();
2709	$e = $x->exponent();
2710	$y = $m * ( 10 ** $e );
2711	print "ok\n" if $x == $y;
2712
2713C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2714
2715A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2716
2717Currently the mantissa is reduced as much as possible, favouring higher
2718exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2719This might change in the future, so do not depend on it.
2720
2721=head2 Accuracy vs. Precision
2722
2723See also: L<Rounding|Rounding>.
2724
2725Math::BigFloat supports both precision (rounding to a certain place before or
2726after the dot) and accuracy (rounding to a certain number of digits). For a
2727full documentation, examples and tips on these topics please see the large
2728section about rounding in L<Math::BigInt>.
2729
2730Since things like C<sqrt(2)> or C<1 / 3> must presented with a limited
2731accuracy lest a operation consumes all resources, each operation produces
2732no more than the requested number of digits.
2733
2734If there is no gloabl precision or accuracy set, B<and> the operation in
2735question was not called with a requested precision or accuracy, B<and> the
2736input $x has no accuracy or precision set, then a fallback parameter will
2737be used. For historical reasons, it is called C<div_scale> and can be accessed
2738via:
2739
2740	$d = Math::BigFloat->div_scale();		# query
2741	Math::BigFloat->div_scale($n);			# set to $n digits
2742
2743The default value for C<div_scale> is 40.
2744
2745In case the result of one operation has more digits than specified,
2746it is rounded. The rounding mode taken is either the default mode, or the one
2747supplied to the operation after the I<scale>:
2748
2749	$x = Math::BigFloat->new(2);
2750	Math::BigFloat->accuracy(5);		# 5 digits max
2751	$y = $x->copy()->bdiv(3);		# will give 0.66667
2752	$y = $x->copy()->bdiv(3,6);		# will give 0.666667
2753	$y = $x->copy()->bdiv(3,6,undef,'odd');	# will give 0.666667
2754	Math::BigFloat->round_mode('zero');
2755	$y = $x->copy()->bdiv(3,6);		# will also give 0.666667
2756
2757Note that C<< Math::BigFloat->accuracy() >> and C<< Math::BigFloat->precision() >>
2758set the global variables, and thus B<any> newly created number will be subject
2759to the global rounding B<immidiately>. This means that in the examples above, the
2760C<3> as argument to C<bdiv()> will also get an accuracy of B<5>.
2761
2762It is less confusing to either calculate the result fully, and afterwards
2763round it explicitely, or use the additional parameters to the math
2764functions like so:
2765
2766	use Math::BigFloat;
2767	$x = Math::BigFloat->new(2);
2768	$y = $x->copy()->bdiv(3);
2769	print $y->bround(5),"\n";		# will give 0.66667
2770
2771	or
2772
2773	use Math::BigFloat;
2774	$x = Math::BigFloat->new(2);
2775	$y = $x->copy()->bdiv(3,5);		# will give 0.66667
2776	print "$y\n";
2777
2778=head2 Rounding
2779
2780=over 2
2781
2782=item ffround ( +$scale )
2783
2784Rounds to the $scale'th place left from the '.', counting from the dot.
2785The first digit is numbered 1.
2786
2787=item ffround ( -$scale )
2788
2789Rounds to the $scale'th place right from the '.', counting from the dot.
2790
2791=item ffround ( 0 )
2792
2793Rounds to an integer.
2794
2795=item fround  ( +$scale )
2796
2797Preserves accuracy to $scale digits from the left (aka significant digits)
2798and pads the rest with zeros. If the number is between 1 and -1, the
2799significant digits count from the first non-zero after the '.'
2800
2801=item fround  ( -$scale ) and fround ( 0 )
2802
2803These are effectively no-ops.
2804
2805=back
2806
2807All rounding functions take as a second parameter a rounding mode from one of
2808the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2809
2810The default rounding mode is 'even'. By using
2811C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
2812mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2813no longer supported.
2814The second parameter to the round functions then overrides the default
2815temporarily.
2816
2817The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
2818'trunc' as rounding mode to make it equivalent to:
2819
2820	$x = 2.5;
2821	$y = int($x) + 2;
2822
2823You can override this by passing the desired rounding mode as parameter to
2824C<as_number()>:
2825
2826	$x = Math::BigFloat->new(2.5);
2827	$y = $x->as_number('odd');	# $y = 3
2828
2829=head1 METHODS
2830
2831=head2 accuracy
2832
2833        $x->accuracy(5);                # local for $x
2834        CLASS->accuracy(5);             # global for all members of CLASS
2835                                        # Note: This also applies to new()!
2836
2837        $A = $x->accuracy();            # read out accuracy that affects $x
2838        $A = CLASS->accuracy();         # read out global accuracy
2839
2840Set or get the global or local accuracy, aka how many significant digits the
2841results have. If you set a global accuracy, then this also applies to new()!
2842
2843Warning! The accuracy I<sticks>, e.g. once you created a number under the
2844influence of C<< CLASS->accuracy($A) >>, all results from math operations with
2845that number will also be rounded.
2846
2847In most cases, you should probably round the results explicitely using one of
2848L<round()>, L<bround()> or L<bfround()> or by passing the desired accuracy
2849to the math operation as additional parameter:
2850
2851        my $x = Math::BigInt->new(30000);
2852        my $y = Math::BigInt->new(7);
2853        print scalar $x->copy()->bdiv($y, 2);           # print 4300
2854        print scalar $x->copy()->bdiv($y)->bround(2);   # print 4300
2855
2856=head2 precision()
2857
2858        $x->precision(-2);      # local for $x, round at the second digit right of the dot
2859        $x->precision(2);       # ditto, round at the second digit left of the dot
2860
2861        CLASS->precision(5);    # Global for all members of CLASS
2862                                # This also applies to new()!
2863        CLASS->precision(-5);   # ditto
2864
2865        $P = CLASS->precision();        # read out global precision
2866        $P = $x->precision();           # read out precision that affects $x
2867
2868Note: You probably want to use L<accuracy()> instead. With L<accuracy> you
2869set the number of digits each result should have, with L<precision> you
2870set the place where to round!
2871
2872=head1 Autocreating constants
2873
2874After C<use Math::BigFloat ':constant'> all the floating point constants
2875in the given scope are converted to C<Math::BigFloat>. This conversion
2876happens at compile time.
2877
2878In particular
2879
2880  perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2881
2882prints the value of C<2E-100>. Note that without conversion of
2883constants the expression 2E-100 will be calculated as normal floating point
2884number.
2885
2886Please note that ':constant' does not affect integer constants, nor binary
2887nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2888work.
2889
2890=head2 Math library
2891
2892Math with the numbers is done (by default) by a module called
2893Math::BigInt::Calc. This is equivalent to saying:
2894
2895	use Math::BigFloat lib => 'Calc';
2896
2897You can change this by using:
2898
2899	use Math::BigFloat lib => 'BitVect';
2900
2901The following would first try to find Math::BigInt::Foo, then
2902Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2903
2904	use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2905
2906Calc.pm uses as internal format an array of elements of some decimal base
2907(usually 1e7, but this might be differen for some systems) with the least
2908significant digit first, while BitVect.pm uses a bit vector of base 2, most
2909significant bit first. Other modules might use even different means of
2910representing the numbers. See the respective module documentation for further
2911details.
2912
2913Please note that Math::BigFloat does B<not> use the denoted library itself,
2914but it merely passes the lib argument to Math::BigInt. So, instead of the need
2915to do:
2916
2917	use Math::BigInt lib => 'GMP';
2918	use Math::BigFloat;
2919
2920you can roll it all into one line:
2921
2922	use Math::BigFloat lib => 'GMP';
2923
2924It is also possible to just require Math::BigFloat:
2925
2926	require Math::BigFloat;
2927
2928This will load the neccessary things (like BigInt) when they are needed, and
2929automatically.
2930
2931Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
2932you ever wanted to know about loading a different library.
2933
2934=head2 Using Math::BigInt::Lite
2935
2936It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2937
2938        # 1
2939        use Math::BigFloat with => 'Math::BigInt::Lite';
2940
2941There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2942can combine these if you want. For instance, you may want to use
2943Math::BigInt objects in your main script, too.
2944
2945        # 2
2946        use Math::BigInt;
2947        use Math::BigFloat with => 'Math::BigInt::Lite';
2948
2949Of course, you can combine this with the C<lib> parameter.
2950
2951        # 3
2952        use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2953
2954There is no need for a "use Math::BigInt;" statement, even if you want to
2955use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
2956always loads it. But if you add it, add it B<before>:
2957
2958        # 4
2959        use Math::BigInt;
2960        use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2961
2962Notice that the module with the last C<lib> will "win" and thus
2963it's lib will be used if the lib is available:
2964
2965        # 5
2966        use Math::BigInt lib => 'Bar,Baz';
2967        use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2968
2969That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2970words, Math::BigFloat will try to retain previously loaded libs when you
2971don't specify it onem but if you specify one, it will try to load them.
2972
2973Actually, the lib loading order would be "Bar,Baz,Calc", and then
2974"Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2975same as trying the latter load alone, except for the fact that one of Bar or
2976Baz might be loaded needlessly in an intermidiate step (and thus hang around
2977and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
2978will still be tried to be loaded, but this is not as time/memory consuming as
2979actually loading one of them. Still, this type of usage is not recommended due
2980to these issues.
2981
2982The old way (loading the lib only in BigInt) still works though:
2983
2984        # 6
2985        use Math::BigInt lib => 'Bar,Baz';
2986        use Math::BigFloat;
2987
2988You can even load Math::BigInt afterwards:
2989
2990        # 7
2991        use Math::BigFloat;
2992        use Math::BigInt lib => 'Bar,Baz';
2993
2994But this has the same problems like #5, it will first load Calc
2995(Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
2996Baz, depending on which of them works and is usable/loadable. Since this
2997loads Calc unnecc., it is not recommended.
2998
2999Since it also possible to just require Math::BigFloat, this poses the question
3000about what libary this will use:
3001
3002	require Math::BigFloat;
3003	my $x = Math::BigFloat->new(123); $x += 123;
3004
3005It will use Calc. Please note that the call to import() is still done, but
3006only when you use for the first time some Math::BigFloat math (it is triggered
3007via any constructor, so the first time you create a Math::BigFloat, the load
3008will happen in the background). This means:
3009
3010	require Math::BigFloat;
3011	Math::BigFloat->import ( lib => 'Foo,Bar' );
3012
3013would be the same as:
3014
3015	use Math::BigFloat lib => 'Foo, Bar';
3016
3017But don't try to be clever to insert some operations in between:
3018
3019	require Math::BigFloat;
3020	my $x = Math::BigFloat->bone() + 4;		# load BigInt and Calc
3021	Math::BigFloat->import( lib => 'Pari' );	# load Pari, too
3022	$x = Math::BigFloat->bone()+4;			# now use Pari
3023
3024While this works, it loads Calc needlessly. But maybe you just wanted that?
3025
3026B<Examples #3 is highly recommended> for daily usage.
3027
3028=head1 BUGS
3029
3030Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
3031
3032=head1 CAVEATS
3033
3034=over 1
3035
3036=item stringify, bstr()
3037
3038Both stringify and bstr() now drop the leading '+'. The old code would return
3039'+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
3040reasoning and details.
3041
3042=item bdiv
3043
3044The following will probably not do what you expect:
3045
3046	print $c->bdiv(123.456),"\n";
3047
3048It prints both quotient and reminder since print works in list context. Also,
3049bdiv() will modify $c, so be carefull. You probably want to use
3050
3051	print $c / 123.456,"\n";
3052	print scalar $c->bdiv(123.456),"\n";  # or if you want to modify $c
3053
3054instead.
3055
3056=item Modifying and =
3057
3058Beware of:
3059
3060	$x = Math::BigFloat->new(5);
3061	$y = $x;
3062
3063It will not do what you think, e.g. making a copy of $x. Instead it just makes
3064a second reference to the B<same> object and stores it in $y. Thus anything
3065that modifies $x will modify $y (except overloaded math operators), and vice
3066versa. See L<Math::BigInt> for details and how to avoid that.
3067
3068=item bpow
3069
3070C<bpow()> now modifies the first argument, unlike the old code which left
3071it alone and only returned the result. This is to be consistent with
3072C<badd()> etc. The first will modify $x, the second one won't:
3073
3074	print bpow($x,$i),"\n"; 	# modify $x
3075	print $x->bpow($i),"\n"; 	# ditto
3076	print $x ** $i,"\n";		# leave $x alone
3077
3078=item precision() vs. accuracy()
3079
3080A common pitfall is to use L<precision()> when you want to round a result to
3081a certain number of digits:
3082
3083	use Math::BigFloat;
3084
3085	Math::BigFloat->precision(4);		# does not do what you think it does
3086	my $x = Math::BigFloat->new(12345);	# rounds $x to "12000"!
3087	print "$x\n";				# print "12000"
3088	my $y = Math::BigFloat->new(3);		# rounds $y to "0"!
3089	print "$y\n";				# print "0"
3090	$z = $x / $y;				# 12000 / 0 => NaN!
3091	print "$z\n";
3092	print $z->precision(),"\n";		# 4
3093
3094Replacing L<precision> with L<accuracy> is probably not what you want, either:
3095
3096	use Math::BigFloat;
3097
3098	Math::BigFloat->accuracy(4);		# enables global rounding:
3099	my $x = Math::BigFloat->new(123456);	# rounded immidiately to "12350"
3100	print "$x\n";				# print "123500"
3101	my $y = Math::BigFloat->new(3);		# rounded to "3
3102	print "$y\n";				# print "3"
3103	print $z = $x->copy()->bdiv($y),"\n";	# 41170
3104	print $z->accuracy(),"\n";		# 4
3105
3106What you want to use instead is:
3107
3108	use Math::BigFloat;
3109
3110	my $x = Math::BigFloat->new(123456);	# no rounding
3111	print "$x\n";				# print "123456"
3112	my $y = Math::BigFloat->new(3);		# no rounding
3113	print "$y\n";				# print "3"
3114	print $z = $x->copy()->bdiv($y,4),"\n";	# 41150
3115	print $z->accuracy(),"\n";		# undef
3116
3117In addition to computing what you expected, the last example also does B<not>
3118"taint" the result with an accuracy or precision setting, which would
3119influence any further operation.
3120
3121=back
3122
3123=head1 SEE ALSO
3124
3125L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
3126L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and  L<Math::BigInt::GMP>.
3127
3128The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
3129because they solve the autoupgrading/downgrading issue, at least partly.
3130
3131The package at
3132L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
3133more documentation including a full version history, testcases, empty
3134subclass files and benchmarks.
3135
3136=head1 LICENSE
3137
3138This program is free software; you may redistribute it and/or modify it under
3139the same terms as Perl itself.
3140
3141=head1 AUTHORS
3142
3143Mark Biggar, overloaded interface by Ilya Zakharevich.
3144Completely rewritten by Tels L<http://bloodgate.com> in 2001 - 2004, and still
3145at it in 2005.
3146
3147=cut
3148