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