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