This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Upgrade Math::BigInt from 1.999701 to 1.999704
[perl5.git] / cpan / Math-BigInt / 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.999704';
16 require 5.006002;
17
18 require Exporter;
19 @ISA            = qw/Math::BigInt/;
20 @EXPORT_OK      = qw/bpi/;
21
22 use strict;
23 # $_trap_inf/$_trap_nan are internal and should never be accessed from outside
24 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode
25             $upgrade $downgrade $_trap_nan $_trap_inf/;
26 my $class = "Math::BigFloat";
27
28 use overload
29 '<=>'   =>      sub { my $rc = $_[2] ?
30                       ref($_[0])->bcmp($_[1],$_[0]) : 
31                       ref($_[0])->bcmp($_[0],$_[1]); 
32                       $rc = 1 unless defined $rc;
33                       $rc <=> 0;
34                 },
35 # we need '>=' to get things like "1 >= NaN" right:
36 '>='    =>      sub { my $rc = $_[2] ?
37                       ref($_[0])->bcmp($_[1],$_[0]) : 
38                       ref($_[0])->bcmp($_[0],$_[1]);
39                       # if there was a NaN involved, return false
40                       return '' unless defined $rc;
41                       $rc >= 0;
42                 },
43 'int'   =>      sub { $_[0]->as_number() },             # 'trunc' to bigint
44 ;
45
46 ##############################################################################
47 # global constants, flags and assorted stuff
48
49 # the following are public, but their usage is not recommended. Use the
50 # accessor methods instead.
51
52 # class constants, use Class->constant_name() to access
53 # one of 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'
54 $round_mode = 'even';
55 $accuracy   = undef;
56 $precision  = undef;
57 $div_scale  = 40;
58
59 $upgrade = undef;
60 $downgrade = undef;
61 # the package we are using for our private parts, defaults to:
62 # Math::BigInt->config()->{lib}
63 my $MBI = 'Math::BigInt::Calc';
64
65 # are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
66 $_trap_nan = 0;
67 # the same for infinity
68 $_trap_inf = 0;
69
70 # constant for easier life
71 my $nan = 'NaN'; 
72
73 my $IMPORT = 0; # was import() called yet? used to make require work
74
75 # some digits of accuracy for blog(undef,10); which we use in blog() for speed
76 my $LOG_10 = 
77  '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
78 my $LOG_10_A = length($LOG_10)-1;
79 # ditto for log(2)
80 my $LOG_2 = 
81  '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
82 my $LOG_2_A = length($LOG_2)-1;
83 my $HALF = '0.5';                       # made into an object if nec.
84
85 ##############################################################################
86 # the old code had $rnd_mode, so we need to support it, too
87
88 sub TIESCALAR   { my ($class) = @_; bless \$round_mode, $class; }
89 sub FETCH       { return $round_mode; }
90 sub STORE       { $rnd_mode = $_[0]->round_mode($_[1]); }
91
92 BEGIN
93   {
94   # when someone sets $rnd_mode, we catch this and check the value to see
95   # whether it is valid or not. 
96   $rnd_mode   = 'even'; tie $rnd_mode, 'Math::BigFloat';
97
98   # we need both of them in this package:
99   *as_int = \&as_number;
100   }
101  
102 ##############################################################################
103
104 {
105   # valid method aliases for AUTOLOAD
106   my %methods = map { $_ => 1 }  
107    qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
108         fint facmp fcmp fzero fnan finf finc fdec ffac fneg
109         fceil ffloor frsft flsft fone flog froot fexp
110       /;
111   # valid methods that can be handed up (for AUTOLOAD)
112   my %hand_ups = map { $_ => 1 }  
113    qw / is_nan is_inf is_negative is_positive is_pos is_neg
114         accuracy precision div_scale round_mode fabs fnot
115         objectify upgrade downgrade
116         bone binf bnan bzero
117         bsub
118       /;
119
120   sub _method_alias { exists $methods{$_[0]||''}; } 
121   sub _method_hand_up { exists $hand_ups{$_[0]||''}; } 
122 }
123
124 ##############################################################################
125 # constructors
126
127 sub new 
128   {
129   # create a new BigFloat object from a string or another bigfloat object. 
130   # _e: exponent
131   # _m: mantissa
132   # sign  => sign (+/-), or "NaN"
133
134   my ($class,$wanted,@r) = @_;
135
136   # avoid numify-calls by not using || on $wanted!
137   return $class->bzero() if !defined $wanted;   # default to 0
138   return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
139
140   $class->import() if $IMPORT == 0;             # make require work
141
142   my $self = {}; bless $self, $class;
143   # shortcut for bigints and its subclasses
144   if ((ref($wanted)) && UNIVERSAL::can( $wanted, "as_number"))
145     {
146     $self->{_m} = $wanted->as_number()->{value}; # get us a bigint copy
147     $self->{_e} = $MBI->_zero();
148     $self->{_es} = '+';
149     $self->{sign} = $wanted->sign();
150     return $self->bnorm();
151     }
152   # else: got a string or something masquerading as number (with overload)
153
154   # handle '+inf', '-inf' first
155   if ($wanted =~ /^[+-]?inf\z/)
156     {
157     return $downgrade->new($wanted) if $downgrade;
158
159     $self->{sign} = $wanted;            # set a default sign for bstr()
160     return $self->binf($wanted);
161     }
162
163   # shortcut for simple forms like '12' that neither have trailing nor leading
164   # zeros
165   if ($wanted =~ /^([+-]?)([1-9][0-9]*[1-9])$/)
166     {
167     $self->{_e} = $MBI->_zero();
168     $self->{_es} = '+';
169     $self->{sign} = $1 || '+';
170     $self->{_m} = $MBI->_new($2);
171     return $self->round(@r) if !$downgrade;
172     }
173
174   my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
175   if (!ref $mis)
176     {
177     if ($_trap_nan)
178       {
179       require Carp;
180       Carp::croak ("$wanted is not a number initialized to $class");
181       }
182     
183     return $downgrade->bnan() if $downgrade;
184     
185     $self->{_e} = $MBI->_zero();
186     $self->{_es} = '+';
187     $self->{_m} = $MBI->_zero();
188     $self->{sign} = $nan;
189     }
190   else
191     {
192     # make integer from mantissa by adjusting exp, then convert to int
193     $self->{_e} = $MBI->_new($$ev);             # exponent
194     $self->{_es} = $$es || '+';
195     my $mantissa = "$$miv$$mfv";                # create mant.
196     $mantissa =~ s/^0+(\d)/$1/;                 # strip leading zeros
197     $self->{_m} = $MBI->_new($mantissa);        # create mant.
198
199     # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
200     if (CORE::length($$mfv) != 0)
201       {
202       my $len = $MBI->_new( CORE::length($$mfv));
203       ($self->{_e}, $self->{_es}) =
204         _e_sub ($self->{_e}, $len, $self->{_es}, '+');
205       }
206     # we can only have trailing zeros on the mantissa if $$mfv eq ''
207     else
208       {
209       # Use a regexp to count the trailing zeros in $$miv instead of _zeros()
210       # because that is faster, especially when _m is not stored in base 10.
211       my $zeros = 0; $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/; 
212       if ($zeros != 0)
213         {
214         my $z = $MBI->_new($zeros);
215         # turn '120e2' into '12e3'
216         $MBI->_rsft ( $self->{_m}, $z, 10);
217         ($self->{_e}, $self->{_es}) =
218           _e_add ( $self->{_e}, $z, $self->{_es}, '+');
219         }
220       }
221     $self->{sign} = $$mis;
222
223     # for something like 0Ey, set y to 1, and -0 => +0
224     # Check $$miv for being '0' and $$mfv eq '', because otherwise _m could not
225     # have become 0. That's faster than to call $MBI->_is_zero().
226     $self->{sign} = '+', $self->{_e} = $MBI->_one()
227      if $$miv eq '0' and $$mfv eq '';
228
229     return $self->round(@r) if !$downgrade;
230     }
231   # if downgrade, inf, NaN or integers go down
232
233   if ($downgrade && $self->{_es} eq '+')
234     {
235     if ($MBI->_is_zero( $self->{_e} ))
236       {
237       return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
238       }
239     return $downgrade->new($self->bsstr()); 
240     }
241   $self->bnorm()->round(@r);                    # first normalize, then round
242   }
243
244 sub copy
245   {
246   # if two arguments, the first one is the class to "swallow" subclasses
247   if (@_ > 1)
248     {
249     my  $self = bless {
250         sign => $_[1]->{sign}, 
251         _es => $_[1]->{_es}, 
252         _m => $MBI->_copy($_[1]->{_m}),
253         _e => $MBI->_copy($_[1]->{_e}),
254     }, $_[0] if @_ > 1;
255
256     $self->{_a} = $_[1]->{_a} if defined $_[1]->{_a};
257     $self->{_p} = $_[1]->{_p} if defined $_[1]->{_p};
258     return $self;
259     }
260
261   my $self = bless {
262         sign => $_[0]->{sign}, 
263         _es => $_[0]->{_es}, 
264         _m => $MBI->_copy($_[0]->{_m}),
265         _e => $MBI->_copy($_[0]->{_e}),
266         }, ref($_[0]);
267
268   $self->{_a} = $_[0]->{_a} if defined $_[0]->{_a};
269   $self->{_p} = $_[0]->{_p} if defined $_[0]->{_p};
270   $self;
271   }
272
273 sub _bnan
274   {
275   # used by parent class bone() to initialize number to NaN
276   my $self = shift;
277   
278   if ($_trap_nan)
279     {
280     require Carp;
281     my $class = ref($self);
282     Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
283     }
284
285   $IMPORT=1;                                    # call our import only once
286   $self->{_m} = $MBI->_zero();
287   $self->{_e} = $MBI->_zero();
288   $self->{_es} = '+';
289   }
290
291 sub _binf
292   {
293   # used by parent class bone() to initialize number to +-inf
294   my $self = shift;
295   
296   if ($_trap_inf)
297     {
298     require Carp;
299     my $class = ref($self);
300     Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
301     }
302
303   $IMPORT=1;                                    # call our import only once
304   $self->{_m} = $MBI->_zero();
305   $self->{_e} = $MBI->_zero();
306   $self->{_es} = '+';
307   }
308
309 sub _bone
310   {
311   # used by parent class bone() to initialize number to 1
312   my $self = shift;
313   $IMPORT=1;                                    # call our import only once
314   $self->{_m} = $MBI->_one();
315   $self->{_e} = $MBI->_zero();
316   $self->{_es} = '+';
317   }
318
319 sub _bzero
320   {
321   # used by parent class bone() to initialize number to 0
322   my $self = shift;
323   $IMPORT=1;                                    # call our import only once
324   $self->{_m} = $MBI->_zero();
325   $self->{_e} = $MBI->_one();
326   $self->{_es} = '+';
327   }
328
329 sub isa
330   {
331   my ($self,$class) = @_;
332   return if $class =~ /^Math::BigInt/;          # we aren't one of these
333   UNIVERSAL::isa($self,$class);
334   }
335
336 sub config
337   {
338   # return (later set?) configuration data as hash ref
339   my $class = shift || 'Math::BigFloat';
340
341   if (@_ == 1 && ref($_[0]) ne 'HASH')
342     {
343     my $cfg = $class->SUPER::config();
344     return $cfg->{$_[0]};
345     }
346
347   my $cfg = $class->SUPER::config(@_);
348
349   # now we need only to override the ones that are different from our parent
350   $cfg->{class} = $class;
351   $cfg->{with} = $MBI;
352   $cfg;
353   }
354
355 ##############################################################################
356 # string conversion
357
358 sub bstr 
359   {
360   # (ref to BFLOAT or num_str ) return num_str
361   # Convert number from internal format to (non-scientific) string format.
362   # internal format is always normalized (no leading zeros, "-0" => "+0")
363   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
364
365   if ($x->{sign} !~ /^[+-]$/)
366     {
367     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
368     return 'inf';                                       # +inf
369     }
370
371   my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
372
373   # $x is zero?
374   my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
375   if ($not_zero)
376     {
377     $es = $MBI->_str($x->{_m});
378     $len = CORE::length($es);
379     my $e = $MBI->_num($x->{_e});       
380     $e = -$e if $x->{_es} eq '-';
381     if ($e < 0)
382       {
383       $dot = '';
384       # if _e is bigger than a scalar, the following will blow your memory
385       if ($e <= -$len)
386         {
387         my $r = abs($e) - $len;
388         $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
389         }
390       else
391         {
392         substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
393         $cad = -$cad if $x->{_es} eq '-';
394         }
395       }
396     elsif ($e > 0)
397       {
398       # expand with zeros
399       $es .= '0' x $e; $len += $e; $cad = 0;
400       }
401     } # if not zero
402
403   $es = '-'.$es if $x->{sign} eq '-';
404   # if set accuracy or precision, pad with zeros on the right side
405   if ((defined $x->{_a}) && ($not_zero))
406     {
407     # 123400 => 6, 0.1234 => 4, 0.001234 => 4
408     my $zeros = $x->{_a} - $cad;                # cad == 0 => 12340
409     $zeros = $x->{_a} - $len if $cad != $len;
410     $es .= $dot.'0' x $zeros if $zeros > 0;
411     }
412   elsif ((($x->{_p} || 0) < 0))
413     {
414     # 123400 => 6, 0.1234 => 4, 0.001234 => 6
415     my $zeros = -$x->{_p} + $cad;
416     $es .= $dot.'0' x $zeros if $zeros > 0;
417     }
418   $es;
419   }
420
421 sub bsstr
422   {
423   # (ref to BFLOAT or num_str ) return num_str
424   # Convert number from internal format to scientific string format.
425   # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
426   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
427
428   if ($x->{sign} !~ /^[+-]$/)
429     {
430     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
431     return 'inf';                                       # +inf
432     }
433   my $sep = 'e'.$x->{_es};
434   my $sign = $x->{sign}; $sign = '' if $sign eq '+';
435   $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
436   }
437     
438 sub numify 
439   {
440   # Convert a Perl scalar number from a BigFloat object.
441   # Create a string and let Perl's atoi()/atof() handle the rest.
442   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
443   return 0 + $x->bsstr(); 
444   }
445
446 ##############################################################################
447 # public stuff (usually prefixed with "b")
448
449 sub bneg
450   {
451   # (BINT or num_str) return BINT
452   # negate number or make a negated number from string
453   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
454
455   return $x if $x->modify('bneg');
456
457   # for +0 do not negate (to have always normalized +0). Does nothing for 'NaN'
458   $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
459   $x;
460   }
461
462 # tels 2001-08-04 
463 # XXX TODO this must be overwritten and return NaN for non-integer values
464 # band(), bior(), bxor(), too
465 #sub bnot
466 #  {
467 #  $class->SUPER::bnot($class,@_);
468 #  }
469
470 sub bcmp 
471   {
472   # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
473
474   # set up parameters
475   my ($self,$x,$y) = (ref($_[0]),@_);
476
477   # objectify is costly, so avoid it
478   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
479     {
480     ($self,$x,$y) = objectify(2,@_);
481     }
482
483   return $upgrade->bcmp($x,$y) if defined $upgrade &&
484     ((!$x->isa($self)) || (!$y->isa($self)));
485
486   # Handle all 'nan' cases.
487
488   return undef if ($x->{sign} eq $nan) || ($y->{sign} eq $nan);
489
490   # Handle all '+inf' and '-inf' cases.
491
492   return  0 if ($x->{sign} eq '+inf' && $y->{sign} eq '+inf' ||
493                 $x->{sign} eq '-inf' && $y->{sign} eq '-inf');
494   return +1 if $x->{sign} eq '+inf';    # x = +inf and y < +inf
495   return -1 if $x->{sign} eq '-inf';    # x = -inf and y > -inf
496   return -1 if $y->{sign} eq '+inf';    # x < +inf and y = +inf
497   return +1 if $y->{sign} eq '-inf';    # x > -inf and y = -inf
498
499   # Handle all cases with opposite signs.
500
501   return +1 if $x->{sign} eq '+' && $y->{sign} eq '-';  # also does 0 <=> -y
502   return -1 if $x->{sign} eq '-' && $y->{sign} eq '+';  # also does -x <=> 0
503
504   # Handle all remaining zero cases.
505
506   my $xz = $x->is_zero();
507   my $yz = $y->is_zero();
508   return  0 if $xz && $yz;                              # 0 <=> 0
509   return -1 if $xz && $y->{sign} eq '+';                # 0 <=> +y
510   return +1 if $yz && $x->{sign} eq '+';                # +x <=> 0
511
512   # Both arguments are now finite, non-zero numbers with the same sign.
513
514   my $cmp;
515
516   # The next step is to compare the exponents, but since each mantissa is an
517   # integer of arbitrary value, the exponents must be normalized by the length
518   # of the mantissas before we can compare them.
519
520   my $mxl = $MBI->_len($x->{_m});
521   my $myl = $MBI->_len($y->{_m});
522
523   # If the mantissas have the same length, there is no point in normalizing the
524   # exponents by the length of the mantissas, so treat that as a special case.
525
526   if ($mxl == $myl) {
527
528       # First handle the two cases where the exponents have different signs.
529
530       if ($x->{_es} eq '+' && $y->{_es} eq '-') {
531           $cmp = +1;
532       }
533
534       elsif ($x->{_es} eq '-' && $y->{_es} eq '+') {
535           $cmp = -1;
536       }
537
538       # Then handle the case where the exponents have the same sign.
539
540       else {
541           $cmp = $MBI->_acmp($x->{_e}, $y->{_e});
542           $cmp = -$cmp if $x->{_es} eq '-';
543       }
544
545       # Adjust for the sign, which is the same for x and y, and bail out if
546       # we're done.
547
548       $cmp = -$cmp if $x->{sign} eq '-';        # 124 > 123, but -124 < -123
549       return $cmp if $cmp;
550
551   }
552
553   # We must normalize each exponent by the length of the corresponding
554   # mantissa. Life is a lot easier if we first make both exponents
555   # non-negative. We do this by adding the same positive value to both
556   # exponent. This is safe, because when comparing the exponents, only the
557   # relative difference is important.
558
559   my $ex;
560   my $ey;
561
562   if ($x->{_es} eq '+') {
563
564       # If the exponent of x is >= 0 and the exponent of y is >= 0, there is no
565       # need to do anything special.
566
567       if ($y->{_es} eq '+') {
568           $ex = $MBI->_copy($x->{_e});
569           $ey = $MBI->_copy($y->{_e});
570       }
571
572       # If the exponent of x is >= 0 and the exponent of y is < 0, add the
573       # absolute value of the exponent of y to both.
574
575       else {
576           $ex = $MBI->_copy($x->{_e});
577           $ex = $MBI->_add($ex, $y->{_e});      # ex + |ey|
578           $ey = $MBI->_zero();                  # -ex + |ey| = 0
579       }
580
581   } else {
582
583       # If the exponent of x is < 0 and the exponent of y is >= 0, add the
584       # absolute value of the exponent of x to both.
585
586       if ($y->{_es} eq '+') {
587           $ex = $MBI->_zero();                  # -ex + |ex| = 0
588           $ey = $MBI->_copy($y->{_e});
589           $ey = $MBI->_add($ey, $x->{_e});      # ey + |ex|
590       }
591
592       # If the exponent of x is < 0 and the exponent of y is < 0, add the
593       # absolute values of both exponents to both exponents.
594
595       else {
596           $ex = $MBI->_copy($y->{_e});          # -ex + |ey| + |ex| = |ey|
597           $ey = $MBI->_copy($x->{_e});          # -ey + |ex| + |ey| = |ex|
598       }
599
600   }
601
602   # Now we can normalize the exponents by adding lengths of the mantissas.
603
604   $MBI->_add($ex, $MBI->_new($mxl));
605   $MBI->_add($ey, $MBI->_new($myl));
606
607   # We're done if the exponents are different.
608
609   $cmp = $MBI->_acmp($ex, $ey);
610   $cmp = -$cmp if $x->{sign} eq '-';            # 124 > 123, but -124 < -123
611   return $cmp if $cmp;
612
613   # Compare the mantissas, but first normalize them by padding the shorter
614   # mantissa with zeros (shift left) until it has the same length as the longer
615   # mantissa.
616
617   my $mx = $x->{_m};
618   my $my = $y->{_m};
619
620   if ($mxl > $myl) {
621       $my = $MBI->_lsft($MBI->_copy($my), $MBI->_new($mxl - $myl), 10);
622   } elsif ($mxl < $myl) {
623       $mx = $MBI->_lsft($MBI->_copy($mx), $MBI->_new($myl - $mxl), 10);
624   }
625
626   $cmp = $MBI->_acmp($mx, $my);
627   $cmp = -$cmp if $x->{sign} eq '-';            # 124 > 123, but -124 < -123
628   return $cmp;
629
630   }
631
632 sub bacmp 
633   {
634   # Compares 2 values, ignoring their signs. 
635   # Returns one of undef, <0, =0, >0. (suitable for sort)
636   
637   # set up parameters
638   my ($self,$x,$y) = (ref($_[0]),@_);
639   # objectify is costly, so avoid it
640   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
641     {
642     ($self,$x,$y) = objectify(2,@_);
643     }
644
645   return $upgrade->bacmp($x,$y) if defined $upgrade &&
646     ((!$x->isa($self)) || (!$y->isa($self)));
647
648   # handle +-inf and NaN's
649   if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
650     {
651     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
652     return 0 if ($x->is_inf() && $y->is_inf());
653     return 1 if ($x->is_inf() && !$y->is_inf());
654     return -1;
655     }
656
657   # shortcut 
658   my $xz = $x->is_zero();
659   my $yz = $y->is_zero();
660   return 0 if $xz && $yz;                               # 0 <=> 0
661   return -1 if $xz && !$yz;                             # 0 <=> +y
662   return 1 if $yz && !$xz;                              # +x <=> 0
663
664   # adjust so that exponents are equal
665   my $lxm = $MBI->_len($x->{_m});
666   my $lym = $MBI->_len($y->{_m});
667   my ($xes,$yes) = (1,1);
668   $xes = -1 if $x->{_es} ne '+';
669   $yes = -1 if $y->{_es} ne '+';
670   # the numify somewhat limits our length, but makes it much faster
671   my $lx = $lxm + $xes * $MBI->_num($x->{_e});
672   my $ly = $lym + $yes * $MBI->_num($y->{_e});
673   my $l = $lx - $ly;
674   return $l <=> 0 if $l != 0;
675   
676   # lengths (corrected by exponent) are equal
677   # so make mantissa equal-length by padding with zero (shift left)
678   my $diff = $lxm - $lym;
679   my $xm = $x->{_m};            # not yet copy it
680   my $ym = $y->{_m};
681   if ($diff > 0)
682     {
683     $ym = $MBI->_copy($y->{_m});
684     $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
685     }
686   elsif ($diff < 0)
687     {
688     $xm = $MBI->_copy($x->{_m});
689     $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
690     }
691   $MBI->_acmp($xm,$ym);
692   }
693
694 sub badd 
695   {
696   # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
697   # return result as BFLOAT
698
699   # set up parameters
700   my ($self,$x,$y,@r) = (ref($_[0]),@_);
701   # objectify is costly, so avoid it
702   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
703     {
704     ($self,$x,$y,@r) = objectify(2,@_);
705     }
706  
707   return $x if $x->modify('badd');
708
709   # inf and NaN handling
710   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
711     {
712     # NaN first
713     return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
714     # inf handling
715     if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
716       {
717       # +inf++inf or -inf+-inf => same, rest is NaN
718       return $x if $x->{sign} eq $y->{sign};
719       return $x->bnan();
720       }
721     # +-inf + something => +inf; something +-inf => +-inf
722     $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
723     return $x;
724     }
725
726   return $upgrade->badd($x,$y,@r) if defined $upgrade &&
727    ((!$x->isa($self)) || (!$y->isa($self)));
728
729   $r[3] = $y;                                           # no push!
730
731   # speed: no add for 0+y or x+0
732   return $x->bround(@r) if $y->is_zero();               # x+0
733   if ($x->is_zero())                                    # 0+y
734     {
735     # make copy, clobbering up x (modify in place!)
736     $x->{_e} = $MBI->_copy($y->{_e});
737     $x->{_es} = $y->{_es};
738     $x->{_m} = $MBI->_copy($y->{_m});
739     $x->{sign} = $y->{sign} || $nan;
740     return $x->round(@r);
741     }
742  
743   # take lower of the two e's and adapt m1 to it to match m2
744   my $e = $y->{_e};
745   $e = $MBI->_zero() if !defined $e;            # if no BFLOAT?
746   $e = $MBI->_copy($e);                         # make copy (didn't do it yet)
747
748   my $es;
749
750   ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});
751
752   my $add = $MBI->_copy($y->{_m});
753
754   if ($es eq '-')                               # < 0
755     {
756     $MBI->_lsft( $x->{_m}, $e, 10);
757     ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
758     }
759   elsif (!$MBI->_is_zero($e))                   # > 0
760     {
761     $MBI->_lsft($add, $e, 10);
762     }
763   # else: both e are the same, so just leave them
764
765   if ($x->{sign} eq $y->{sign})
766     {
767     # add
768     $x->{_m} = $MBI->_add($x->{_m}, $add);
769     }
770   else
771     {
772     ($x->{_m}, $x->{sign}) = 
773      _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
774     }
775
776   # delete trailing zeros, then round
777   $x->bnorm()->round(@r);
778   }
779
780 # sub bsub is inherited from Math::BigInt!
781
782 sub binc
783   {
784   # increment arg by one
785   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
786
787   return $x if $x->modify('binc');
788
789   if ($x->{_es} eq '-')
790     {
791     return $x->badd($self->bone(),@r);  #  digits after dot
792     }
793
794   if (!$MBI->_is_zero($x->{_e}))                # _e == 0 for NaN, inf, -inf
795     {
796     # 1e2 => 100, so after the shift below _m has a '0' as last digit
797     $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10);      # 1e2 => 100
798     $x->{_e} = $MBI->_zero();                           # normalize
799     $x->{_es} = '+';
800     # we know that the last digit of $x will be '1' or '9', depending on the
801     # sign
802     }
803   # now $x->{_e} == 0
804   if ($x->{sign} eq '+')
805     {
806     $MBI->_inc($x->{_m});
807     return $x->bnorm()->bround(@r);
808     }
809   elsif ($x->{sign} eq '-')
810     {
811     $MBI->_dec($x->{_m});
812     $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
813     return $x->bnorm()->bround(@r);
814     }
815   # inf, nan handling etc
816   $x->badd($self->bone(),@r);                   # badd() does round 
817   }
818
819 sub bdec
820   {
821   # decrement arg by one
822   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
823
824   return $x if $x->modify('bdec');
825
826   if ($x->{_es} eq '-')
827     {
828     return $x->badd($self->bone('-'),@r);       #  digits after dot
829     }
830
831   if (!$MBI->_is_zero($x->{_e}))
832     {
833     $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10);      # 1e2 => 100
834     $x->{_e} = $MBI->_zero();                           # normalize
835     $x->{_es} = '+';
836     }
837   # now $x->{_e} == 0
838   my $zero = $x->is_zero();
839   # <= 0
840   if (($x->{sign} eq '-') || $zero)
841     {
842     $MBI->_inc($x->{_m});
843     $x->{sign} = '-' if $zero;                          # 0 => 1 => -1
844     $x->{sign} = '+' if $MBI->_is_zero($x->{_m});       # -1 +1 => -0 => +0
845     return $x->bnorm()->round(@r);
846     }
847   # > 0
848   elsif ($x->{sign} eq '+')
849     {
850     $MBI->_dec($x->{_m});
851     return $x->bnorm()->round(@r);
852     }
853   # inf, nan handling etc
854   $x->badd($self->bone('-'),@r);                # does round
855   } 
856
857 sub DEBUG () { 0; }
858
859 sub blog
860   {
861   my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
862
863   # If called as $x -> blog() or $x -> blog(undef), don't objectify the
864   # undefined base, since undef signals that the base is Euler's number.
865   #unless (ref($x) && !defined($base)) {
866   #    # objectify is costly, so avoid it
867   #    if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) {
868   #        ($self,$x,$base,$a,$p,$r) = objectify(2,@_);
869   #    }
870   #}
871
872   return $x if $x->modify('blog');
873
874   return $x -> bnan() if $x -> is_nan();
875
876   # we need to limit the accuracy to protect against overflow
877   my $fallback = 0;
878   my ($scale,@params);
879   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
880
881   # no rounding at all, so must use fallback
882   if (scalar @params == 0)
883     {
884     # simulate old behaviour
885     $params[0] = $self->div_scale();    # and round to it as accuracy
886     $params[1] = undef;                 # P = undef
887     $scale = $params[0]+4;              # at least four more for proper round
888     $params[2] = $r;                    # round mode by caller or undef
889     $fallback = 1;                      # to clear a/p afterwards
890     }
891   else
892     {
893     # the 4 below is empirical, and there might be cases where it is not
894     # enough...
895     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
896     }
897
898   my $done = 0;
899   if (defined $base) {
900       $base = $self -> new($base) unless ref $base;
901       if ($base -> is_nan() || $base -> is_one()) {
902           $x -> bnan();
903           $done = 1;
904       } elsif ($base -> is_inf() || $base -> is_zero()) {
905           if ($x -> is_inf() || $x -> is_zero()) {
906               $x -> bnan();
907           } else {
908               $x -> bzero(@params);
909           }
910           $done = 1;
911       } elsif ($base -> is_negative()) {        # -inf < base < 0
912           if ($x -> is_one()) {                 #     x = 1
913               $x -> bzero(@params);
914           } elsif ($x == $base) {
915               $x -> bone('+', @params);         #     x = base
916           } else {
917               $x -> bnan();                     #     otherwise
918           }
919           $done = 1;
920       } elsif ($x == $base) {
921           $x -> bone('+', @params);             # 0 < base && 0 < x < inf
922           $done = 1;
923       }
924   }
925
926   # We now know that the base is either undefined or positive and finite.
927
928   unless ($done) {
929       if ($x -> is_inf()) {             #   x = +/-inf
930           my $sign = defined $base && $base < 1 ? '-' : '+';
931           $x -> binf($sign);
932           $done = 1;
933       } elsif ($x -> is_neg()) {        #   -inf < x < 0
934           $x -> bnan();
935           $done = 1;
936       } elsif ($x -> is_one()) {        #   x = 1
937           $x -> bzero(@params);
938           $done = 1;
939       } elsif ($x -> is_zero()) {       #   x = 0
940           my $sign = defined $base && $base < 1 ? '+' : '-';
941           $x -> binf($sign);
942           $done = 1;
943       }
944   }
945
946   if ($done) {
947       if ($fallback) {
948         # clear a/p after round, since user did not request it
949           delete $x->{_a};
950           delete $x->{_p};
951         }
952       return $x;
953       }
954
955   # when user set globals, they would interfere with our calculation, so
956   # disable them and later re-enable them
957   no strict 'refs';
958   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
959   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
960   # we also need to disable any set A or P on $x (_find_round_parameters took
961   # them already into account), since these would interfere, too
962   delete $x->{_a}; delete $x->{_p};
963   # need to disable $upgrade in BigInt, to avoid deep recursion
964   local $Math::BigInt::upgrade = undef;
965   local $Math::BigFloat::downgrade = undef;
966
967   # upgrade $x if $x is not a BigFloat (handle BigInt input)
968   # XXX TODO: rebless!
969   if (!$x->isa('Math::BigFloat'))
970     {
971     $x = Math::BigFloat->new($x);
972     $self = ref($x);
973     }
974   
975   $done = 0;
976
977   # If the base is defined and an integer, try to calculate integer result
978   # first. This is very fast, and in case the real result was found, we can
979   # stop right here.
980   if (defined $base && $base->is_int() && $x->is_int())
981     {
982     my $i = $MBI->_copy( $x->{_m} );
983     $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
984     my $int = Math::BigInt->bzero();
985     $int->{value} = $i;
986     $int->blog($base->as_number());
987     # if ($exact)
988     if ($base->as_number()->bpow($int) == $x)
989       {
990       # found result, return it
991       $x->{_m} = $int->{value};
992       $x->{_e} = $MBI->_zero();
993       $x->{_es} = '+';
994       $x->bnorm();
995       $done = 1;
996       }
997     }
998
999   if ($done == 0)
1000     {
1001     # base is undef, so base should be e (Euler's number), so first calculate the
1002     # log to base e (using reduction by 10 (and probably 2)):
1003     $self->_log_10($x,$scale);
1004
1005     # and if a different base was requested, convert it
1006     if (defined $base)
1007       {
1008       $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
1009       # not ln, but some other base (don't modify $base)
1010       $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
1011       }
1012     }
1013  
1014   # shortcut to not run through _find_round_parameters again
1015   if (defined $params[0])
1016     {
1017     $x->bround($params[0],$params[2]);          # then round accordingly
1018     }
1019   else
1020     {
1021     $x->bfround($params[1],$params[2]);         # then round accordingly
1022     }
1023   if ($fallback)
1024     {
1025     # clear a/p after round, since user did not request it
1026     delete $x->{_a}; delete $x->{_p};
1027     }
1028   # restore globals
1029   $$abr = $ab; $$pbr = $pb;
1030
1031   $x;
1032   }
1033
1034 sub _len_to_steps
1035   {
1036   # Given D (digits in decimal), compute N so that N! (N factorial) is
1037   # at least D digits long. D should be at least 50.
1038   my $d = shift;
1039
1040   # two constants for the Ramanujan estimate of ln(N!)
1041   my $lg2 = log(2 * 3.14159265) / 2;
1042   my $lg10 = log(10);
1043
1044   # D = 50 => N => 42, so L = 40 and R = 50
1045   my $l = 40; my $r = $d;
1046
1047   # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
1048   $l = $l->numify if ref($l);
1049   $r = $r->numify if ref($r);
1050   $lg2 = $lg2->numify if ref($lg2);
1051   $lg10 = $lg10->numify if ref($lg10);
1052
1053   # binary search for the right value (could this be written as the reverse of lg(n!)?)
1054   while ($r - $l > 1)
1055     {
1056     my $n = int(($r - $l) / 2) + $l;
1057     my $ramanujan = 
1058       int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
1059     $ramanujan > $d ? $r = $n : $l = $n;
1060     }
1061   $l;
1062   }
1063
1064 sub bnok
1065   {
1066   # Calculate n over k (binomial coefficient or "choose" function) as integer.
1067   # set up parameters
1068   my ($self,$x,$y,@r) = (ref($_[0]),@_);
1069
1070   # objectify is costly, so avoid it
1071   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1072     {
1073     ($self,$x,$y,@r) = objectify(2,@_);
1074     }
1075
1076   return $x if $x->modify('bnok');
1077
1078   return $x->bnan() if $x->is_nan() || $y->is_nan();
1079   return $x->binf() if $x->is_inf();
1080
1081   my $u = $x->as_int();
1082   $u->bnok($y->as_int());
1083
1084   $x->{_m} = $u->{value};
1085   $x->{_e} = $MBI->_zero();
1086   $x->{_es} = '+';
1087   $x->{sign} = '+';
1088   $x->bnorm(@r);
1089   }
1090
1091 sub bexp
1092   {
1093   # Calculate e ** X (Euler's number to the power of X)
1094   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1095
1096   return $x if $x->modify('bexp');
1097
1098   return $x->binf() if $x->{sign} eq '+inf';
1099   return $x->bzero() if $x->{sign} eq '-inf';
1100
1101   # we need to limit the accuracy to protect against overflow
1102   my $fallback = 0;
1103   my ($scale,@params);
1104   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1105
1106   # also takes care of the "error in _find_round_parameters?" case
1107   return $x if $x->{sign} eq 'NaN';
1108
1109   # no rounding at all, so must use fallback
1110   if (scalar @params == 0)
1111     {
1112     # simulate old behaviour
1113     $params[0] = $self->div_scale();    # and round to it as accuracy
1114     $params[1] = undef;                 # P = undef
1115     $scale = $params[0]+4;              # at least four more for proper round
1116     $params[2] = $r;                    # round mode by caller or undef
1117     $fallback = 1;                      # to clear a/p afterwards
1118     }
1119   else
1120     {
1121     # the 4 below is empirical, and there might be cases where it's not enough...
1122     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1123     }
1124
1125   return $x->bone(@params) if $x->is_zero();
1126
1127   if (!$x->isa('Math::BigFloat'))
1128     {
1129     $x = Math::BigFloat->new($x);
1130     $self = ref($x);
1131     }
1132   
1133   # when user set globals, they would interfere with our calculation, so
1134   # disable them and later re-enable them
1135   no strict 'refs';
1136   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1137   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1138   # we also need to disable any set A or P on $x (_find_round_parameters took
1139   # them already into account), since these would interfere, too
1140   delete $x->{_a}; delete $x->{_p};
1141   # need to disable $upgrade in BigInt, to avoid deep recursion
1142   local $Math::BigInt::upgrade = undef;
1143   local $Math::BigFloat::downgrade = undef;
1144
1145   my $x_org = $x->copy();
1146
1147   # We use the following Taylor series:
1148
1149   #           x    x^2   x^3   x^4
1150   #  e = 1 + --- + --- + --- + --- ...
1151   #           1!    2!    3!    4!
1152
1153   # The difference for each term is X and N, which would result in:
1154   # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term
1155
1156   # But it is faster to compute exp(1) and then raising it to the
1157   # given power, esp. if $x is really big and an integer because:
1158
1159   #  * The numerator is always 1, making the computation faster
1160   #  * the series converges faster in the case of x == 1
1161   #  * We can also easily check when we have reached our limit: when the
1162   #    term to be added is smaller than "1E$scale", we can stop - f.i.
1163   #    scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
1164   #  * we can compute the *exact* result by simulating bigrat math:
1165
1166   #  1   1    gcd(3,4) = 1    1*24 + 1*6    5
1167   #  - + -                  = ---------- =  --                 
1168   #  6   24                      6*24       24
1169
1170   # We do not compute the gcd() here, but simple do:
1171   #  1   1    1*24 + 1*6   30
1172   #  - + -  = --------- =  --                 
1173   #  6   24       6*24     144
1174
1175   # In general:
1176   #  a   c    a*d + c*b         and note that c is always 1 and d = (b*f)
1177   #  - + -  = ---------
1178   #  b   d       b*d
1179
1180   # This leads to:         which can be reduced by b to:
1181   #  a   1     a*b*f + b    a*f + 1
1182   #  - + -   = --------- =  -------
1183   #  b   b*f     b*b*f        b*f
1184
1185   # The first terms in the series are:
1186
1187   # 1     1    1    1    1    1     1     1     13700
1188   # -- + -- + -- + -- + -- + --- + --- + ---- = -----
1189   # 1     1    2    6   24   120   720   5040   5040
1190
1191   # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!
1192
1193   if ($scale <= 75)
1194     {
1195     # set $x directly from a cached string form
1196     $x->{_m} = $MBI->_new(
1197     "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
1198     $x->{sign} = '+';
1199     $x->{_es} = '-';
1200     $x->{_e} = $MBI->_new(79);
1201     }
1202   else
1203     {
1204     # compute A and B so that e = A / B.
1205  
1206     # After some terms we end up with this, so we use it as a starting point:
1207     my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
1208     my $F = $MBI->_new(42); my $step = 42;
1209
1210     # Compute how many steps we need to take to get $A and $B sufficiently big
1211     my $steps = _len_to_steps($scale - 4);
1212 #    print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
1213     while ($step++ <= $steps)
1214       {
1215       # calculate $a * $f + 1
1216       $A = $MBI->_mul($A, $F);
1217       $A = $MBI->_inc($A);
1218       # increment f
1219       $F = $MBI->_inc($F);
1220       }
1221     # compute $B as factorial of $steps (this is faster than doing it manually)
1222     my $B = $MBI->_fac($MBI->_new($steps));
1223     
1224 #  print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";
1225
1226     # compute A/B with $scale digits in the result (truncate, not round)
1227     $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
1228     $A = $MBI->_div( $A, $B );
1229
1230     $x->{_m} = $A;
1231     $x->{sign} = '+';
1232     $x->{_es} = '-';
1233     $x->{_e} = $MBI->_new($scale);
1234     }
1235
1236   # $x contains now an estimate of e, with some surplus digits, so we can round
1237   if (!$x_org->is_one())
1238     {
1239     # raise $x to the wanted power and round it in one step:
1240     $x->bpow($x_org, @params);
1241     }
1242   else
1243     {
1244     # else just round the already computed result
1245     delete $x->{_a}; delete $x->{_p};
1246     # shortcut to not run through _find_round_parameters again
1247     if (defined $params[0])
1248       {
1249       $x->bround($params[0],$params[2]);                # then round accordingly
1250       }
1251     else
1252       {
1253       $x->bfround($params[1],$params[2]);               # then round accordingly
1254       }
1255     }
1256   if ($fallback)
1257     {
1258     # clear a/p after round, since user did not request it
1259     delete $x->{_a}; delete $x->{_p};
1260     }
1261   # restore globals
1262   $$abr = $ab; $$pbr = $pb;
1263
1264   $x;                                           # return modified $x
1265   }
1266
1267 sub _log
1268   {
1269   # internal log function to calculate ln() based on Taylor series.
1270   # Modifies $x in place.
1271   my ($self,$x,$scale) = @_;
1272
1273   # in case of $x == 1, result is 0
1274   return $x->bzero() if $x->is_one();
1275
1276   # XXX TODO: rewrite this in a similar manner to bexp()
1277
1278   # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
1279
1280   # u = x-1, v = x+1
1281   #              _                               _
1282   # Taylor:     |    u    1   u^3   1   u^5       |
1283   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 0
1284   #             |_   v    3   v^3   5   v^5      _|
1285
1286   # This takes much more steps to calculate the result and is thus not used
1287   # u = x-1
1288   #              _                               _
1289   # Taylor:     |    u    1   u^2   1   u^3       |
1290   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 1/2
1291   #             |_   x    2   x^2   3   x^3      _|
1292
1293   my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
1294
1295   $v = $x->copy(); $v->binc();          # v = x+1
1296   $x->bdec(); $u = $x->copy();          # u = x-1; x = x-1
1297   $x->bdiv($v,$scale);                  # first term: u/v
1298   $below = $v->copy();
1299   $over = $u->copy();
1300   $u *= $u; $v *= $v;                           # u^2, v^2
1301   $below->bmul($v);                             # u^3, v^3
1302   $over->bmul($u);
1303   $factor = $self->new(3); $f = $self->new(2);
1304
1305   my $steps = 0 if DEBUG;  
1306   $limit = $self->new("1E-". ($scale-1));
1307   while (3 < 5)
1308     {
1309     # we calculate the next term, and add it to the last
1310     # when the next term is below our limit, it won't affect the outcome
1311     # anymore, so we stop
1312
1313     # calculating the next term simple from over/below will result in quite
1314     # a time hog if the input has many digits, since over and below will
1315     # accumulate more and more digits, and the result will also have many
1316     # digits, but in the end it is rounded to $scale digits anyway. So if we
1317     # round $over and $below first, we save a lot of time for the division
1318     # (not with log(1.2345), but try log (123**123) to see what I mean. This
1319     # can introduce a rounding error if the division result would be f.i.
1320     # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
1321     # if we truncated $over and $below we might get 0.12345. Does this matter
1322     # for the end result? So we give $over and $below 4 more digits to be
1323     # on the safe side (unscientific error handling as usual... :+D
1324
1325     $next = $over->copy->bround($scale+4)->bdiv(
1326       $below->copy->bmul($factor)->bround($scale+4), 
1327       $scale);
1328
1329 ## old version:    
1330 ##    $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
1331
1332     last if $next->bacmp($limit) <= 0;
1333
1334     delete $next->{_a}; delete $next->{_p};
1335     $x->badd($next);
1336     # calculate things for the next term
1337     $over *= $u; $below *= $v; $factor->badd($f);
1338     if (DEBUG)
1339       {
1340       $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
1341       }
1342     }
1343   print "took $steps steps\n" if DEBUG;
1344   $x->bmul($f);                                 # $x *= 2
1345   }
1346
1347 sub _log_10
1348   {
1349   # Internal log function based on reducing input to the range of 0.1 .. 9.99
1350   # and then "correcting" the result to the proper one. Modifies $x in place.
1351   my ($self,$x,$scale) = @_;
1352
1353   # Taking blog() from numbers greater than 10 takes a *very long* time, so we
1354   # break the computation down into parts based on the observation that:
1355   #  blog(X*Y) = blog(X) + blog(Y)
1356   # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
1357   # $x is the faster it gets. Since 2*$x takes about 10 times as
1358   # long, we make it faster by about a factor of 100 by dividing $x by 10.
1359
1360   # The same observation is valid for numbers smaller than 0.1, e.g. computing
1361   # log(1) is fastest, and the further away we get from 1, the longer it takes.
1362   # So we also 'break' this down by multiplying $x with 10 and subtract the
1363   # log(10) afterwards to get the correct result.
1364
1365   # To get $x even closer to 1, we also divide by 2 and then use log(2) to
1366   # correct for this. For instance if $x is 2.4, we use the formula:
1367   #  blog(2.4 * 2) == blog (1.2) + blog(2)
1368   # and thus calculate only blog(1.2) and blog(2), which is faster in total
1369   # than calculating blog(2.4).
1370
1371   # In addition, the values for blog(2) and blog(10) are cached.
1372
1373   # Calculate nr of digits before dot:
1374   my $dbd = $MBI->_num($x->{_e});
1375   $dbd = -$dbd if $x->{_es} eq '-';
1376   $dbd += $MBI->_len($x->{_m});
1377
1378   # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
1379   # infinite recursion
1380
1381   my $calc = 1;                                 # do some calculation?
1382
1383   # disable the shortcut for 10, since we need log(10) and this would recurse
1384   # infinitely deep
1385   if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
1386     {
1387     $dbd = 0;                                   # disable shortcut
1388     # we can use the cached value in these cases
1389     if ($scale <= $LOG_10_A)
1390       {
1391       $x->bzero(); $x->badd($LOG_10);           # modify $x in place
1392       $calc = 0;                                # no need to calc, but round
1393       }
1394     # if we can't use the shortcut, we continue normally
1395     }
1396   else
1397     {
1398     # disable the shortcut for 2, since we maybe have it cached
1399     if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
1400       {
1401       $dbd = 0;                                 # disable shortcut
1402       # we can use the cached value in these cases
1403       if ($scale <= $LOG_2_A)
1404         {
1405         $x->bzero(); $x->badd($LOG_2);          # modify $x in place
1406         $calc = 0;                              # no need to calc, but round
1407         }
1408       # if we can't use the shortcut, we continue normally
1409       }
1410     }
1411
1412   # if $x = 0.1, we know the result must be 0-log(10)
1413   if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
1414       $MBI->_is_one($x->{_m}))
1415     {
1416     $dbd = 0;                                   # disable shortcut
1417     # we can use the cached value in these cases
1418     if ($scale <= $LOG_10_A)
1419       {
1420       $x->bzero(); $x->bsub($LOG_10);
1421       $calc = 0;                                # no need to calc, but round
1422       }
1423     }
1424
1425   return if $calc == 0;                         # already have the result
1426
1427   # default: these correction factors are undef and thus not used
1428   my $l_10;                             # value of ln(10) to A of $scale
1429   my $l_2;                              # value of ln(2) to A of $scale
1430
1431   my $two = $self->new(2);
1432
1433   # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
1434   # so don't do this shortcut for 1 or 0
1435   if (($dbd > 1) || ($dbd < 0))
1436     {
1437     # convert our cached value to an object if not already (avoid doing this
1438     # at import() time, since not everybody needs this)
1439     $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
1440
1441     #print "x = $x, dbd = $dbd, calc = $calc\n";
1442     # got more than one digit before the dot, or more than one zero after the
1443     # dot, so do:
1444     #  log(123)    == log(1.23) + log(10) * 2
1445     #  log(0.0123) == log(1.23) - log(10) * 2
1446   
1447     if ($scale <= $LOG_10_A)
1448       {
1449       # use cached value
1450       $l_10 = $LOG_10->copy();          # copy for mul
1451       }
1452     else
1453       {
1454       # else: slower, compute and cache result
1455       # also disable downgrade for this code path
1456       local $Math::BigFloat::downgrade = undef;
1457
1458       # shorten the time to calculate log(10) based on the following:
1459       # log(1.25 * 8) = log(1.25) + log(8)
1460       #               = log(1.25) + log(2) + log(2) + log(2)
1461
1462       # first get $l_2 (and possible compute and cache log(2))
1463       $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1464       if ($scale <= $LOG_2_A)
1465         {
1466         # use cached value
1467         $l_2 = $LOG_2->copy();                  # copy() for the mul below
1468         }
1469       else
1470         {
1471         # else: slower, compute and cache result
1472         $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1473         $LOG_2 = $l_2->copy();                  # cache the result for later
1474                                                 # the copy() is for mul below
1475         $LOG_2_A = $scale;
1476         }
1477
1478       # now calculate log(1.25):
1479       $l_10 = $self->new('1.25'); $self->_log($l_10, $scale); # scale+4, actually
1480
1481       # log(1.25) + log(2) + log(2) + log(2):
1482       $l_10->badd($l_2);
1483       $l_10->badd($l_2);
1484       $l_10->badd($l_2);
1485       $LOG_10 = $l_10->copy();          # cache the result for later
1486                                         # the copy() is for mul below
1487       $LOG_10_A = $scale;
1488       }
1489     $dbd-- if ($dbd > 1);               # 20 => dbd=2, so make it dbd=1 
1490     $l_10->bmul( $self->new($dbd));     # log(10) * (digits_before_dot-1)
1491     my $dbd_sign = '+';
1492     if ($dbd < 0)
1493       {
1494       $dbd = -$dbd;
1495       $dbd_sign = '-';
1496       }
1497     ($x->{_e}, $x->{_es}) = 
1498         _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
1499  
1500     }
1501
1502   # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1503
1504   ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1505   ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1506
1507   $HALF = $self->new($HALF) unless ref($HALF);
1508
1509   my $twos = 0;                         # default: none (0 times)       
1510   while ($x->bacmp($HALF) <= 0)         # X <= 0.5
1511     {
1512     $twos--; $x->bmul($two);
1513     }
1514   while ($x->bacmp($two) >= 0)          # X >= 2
1515     {
1516     $twos++; $x->bdiv($two,$scale+4);           # keep all digits
1517     }
1518   $x->bround($scale+4);
1519   # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
1520   # So calculate correction factor based on ln(2):
1521   if ($twos != 0)
1522     {
1523     $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1524     if ($scale <= $LOG_2_A)
1525       {
1526       # use cached value
1527       $l_2 = $LOG_2->copy();                    # copy() for the mul below
1528       }
1529     else
1530       {
1531       # else: slower, compute and cache result
1532       # also disable downgrade for this code path
1533       local $Math::BigFloat::downgrade = undef;
1534       $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1535       $LOG_2 = $l_2->copy();                    # cache the result for later
1536                                                 # the copy() is for mul below
1537       $LOG_2_A = $scale;
1538       }
1539     $l_2->bmul($twos);          # * -2 => subtract, * 2 => add
1540     }
1541   else
1542     {
1543       undef $l_2;
1544     }
1545   
1546   $self->_log($x,$scale);                       # need to do the "normal" way
1547   $x->badd($l_10) if defined $l_10;             # correct it by ln(10)
1548   $x->badd($l_2) if defined $l_2;               # and maybe by ln(2)
1549
1550   # all done, $x contains now the result
1551   $x;
1552   }
1553
1554 sub blcm 
1555   { 
1556   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1557   # does not modify arguments, but returns new object
1558   # Lowest Common Multiplicator
1559
1560   my ($self,@arg) = objectify(0,@_);
1561   my $x = $self->new(shift @arg);
1562   while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); } 
1563   $x;
1564   }
1565
1566 sub bgcd
1567   {
1568   # (BINT or num_str, BINT or num_str) return BINT
1569   # does not modify arguments, but returns new object
1570
1571   my $y = shift;
1572   $y = __PACKAGE__->new($y) if !ref($y);
1573   my $self = ref($y);
1574   my $x = $y->copy()->babs();                   # keep arguments
1575
1576   return $x->bnan() if $x->{sign} !~ /^[+-]$/   # x NaN?
1577         || !$x->is_int();                       # only for integers now
1578
1579   while (@_)
1580     {
1581     my $t = shift; $t = $self->new($t) if !ref($t);
1582     $y = $t->copy()->babs();
1583     
1584     return $x->bnan() if $y->{sign} !~ /^[+-]$/ # y NaN?
1585         || !$y->is_int();                       # only for integers now
1586
1587     # greatest common divisor
1588     while (! $y->is_zero())
1589       {
1590       ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
1591       }
1592
1593     last if $x->is_one();
1594     }
1595   $x;
1596   }
1597
1598 ##############################################################################
1599
1600 sub _e_add
1601   {
1602   # Internal helper sub to take two positive integers and their signs and
1603   # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')), 
1604   # output ($CALC,('+'|'-'))
1605   my ($x,$y,$xs,$ys) = @_;
1606
1607   # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1608   if ($xs eq $ys)
1609     {
1610     $x = $MBI->_add ($x, $y );          # a+b
1611     # the sign follows $xs
1612     return ($x, $xs);
1613     }
1614
1615   my $a = $MBI->_acmp($x,$y);
1616   if ($a > 0)
1617     {
1618     $x = $MBI->_sub ($x , $y);                          # abs sub
1619     }
1620   elsif ($a == 0)
1621     {
1622     $x = $MBI->_zero();                                 # result is 0
1623     $xs = '+';
1624     }
1625   else # a < 0
1626     {
1627     $x = $MBI->_sub ( $y, $x, 1 );                      # abs sub
1628     $xs = $ys;
1629     }
1630   ($x,$xs);
1631   }
1632
1633 sub _e_sub
1634   {
1635   # Internal helper sub to take two positive integers and their signs and
1636   # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')), 
1637   # output ($CALC,('+'|'-'))
1638   my ($x,$y,$xs,$ys) = @_;
1639
1640   # flip sign
1641   $ys =~ tr/+-/-+/;
1642   _e_add($x,$y,$xs,$ys);                # call add (does subtract now)
1643   }
1644
1645 ###############################################################################
1646 # is_foo methods (is_negative, is_positive are inherited from BigInt)
1647
1648 sub is_int
1649   {
1650   # return true if arg (BFLOAT or num_str) is an integer
1651   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1652
1653   (($x->{sign} =~ /^[+-]$/) &&                  # NaN and +-inf aren't
1654    ($x->{_es} eq '+')) ? 1 : 0;                 # 1e-1 => no integer
1655   }
1656
1657 sub is_zero
1658   {
1659   # return true if arg (BFLOAT or num_str) is zero
1660   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1661
1662   ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m})) ? 1 : 0;
1663   }
1664
1665 sub is_one
1666   {
1667   # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1668   my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
1669
1670   $sign = '+' if !defined $sign || $sign ne '-';
1671
1672   ($x->{sign} eq $sign && 
1673    $MBI->_is_zero($x->{_e}) &&
1674    $MBI->_is_one($x->{_m}) ) ? 1 : 0; 
1675   }
1676
1677 sub is_odd
1678   {
1679   # return true if arg (BFLOAT or num_str) is odd or false if even
1680   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1681   
1682   (($x->{sign} =~ /^[+-]$/) &&          # NaN & +-inf aren't
1683    ($MBI->_is_zero($x->{_e})) &&
1684    ($MBI->_is_odd($x->{_m}))) ? 1 : 0; 
1685   }
1686
1687 sub is_even
1688   {
1689   # return true if arg (BINT or num_str) is even or false if odd
1690   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
1691
1692   (($x->{sign} =~ /^[+-]$/) &&                  # NaN & +-inf aren't
1693    ($x->{_es} eq '+') &&                        # 123.45 isn't
1694    ($MBI->_is_even($x->{_m}))) ? 1 : 0;         # but 1200 is
1695   }
1696
1697 sub bmul
1698   { 
1699   # multiply two numbers
1700   
1701   # set up parameters
1702   my ($self,$x,$y,@r) = (ref($_[0]),@_);
1703   # objectify is costly, so avoid it
1704   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1705     {
1706     ($self,$x,$y,@r) = objectify(2,@_);
1707     }
1708
1709   return $x if $x->modify('bmul');
1710
1711   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1712
1713   # inf handling
1714   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1715     {
1716     return $x->bnan() if $x->is_zero() || $y->is_zero(); 
1717     # result will always be +-inf:
1718     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1719     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1720     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1721     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1722     return $x->binf('-');
1723     }
1724   
1725   return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1726    ((!$x->isa($self)) || (!$y->isa($self)));
1727
1728   # aEb * cEd = (a*c)E(b+d)
1729   $MBI->_mul($x->{_m},$y->{_m});
1730   ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1731
1732   $r[3] = $y;                           # no push!
1733
1734   # adjust sign:
1735   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1736   $x->bnorm->round(@r);
1737   }
1738
1739 sub bmuladd
1740   { 
1741   # multiply two numbers and add the third to the result
1742   
1743   # set up parameters
1744   my ($self,$x,$y,$z,@r) = objectify(3,@_);
1745
1746   return $x if $x->modify('bmuladd');
1747
1748   return $x->bnan() if (($x->{sign} eq $nan) ||
1749                         ($y->{sign} eq $nan) ||
1750                         ($z->{sign} eq $nan));
1751
1752   # inf handling
1753   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1754     {
1755     return $x->bnan() if $x->is_zero() || $y->is_zero(); 
1756     # result will always be +-inf:
1757     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1758     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1759     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1760     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1761     return $x->binf('-');
1762     }
1763
1764   return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1765    ((!$x->isa($self)) || (!$y->isa($self)));
1766
1767   # aEb * cEd = (a*c)E(b+d)
1768   $MBI->_mul($x->{_m},$y->{_m});
1769   ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1770
1771   $r[3] = $y;                           # no push!
1772
1773   # adjust sign:
1774   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1775
1776   # z=inf handling (z=NaN handled above)
1777   $x->{sign} = $z->{sign}, return $x if $z->{sign} =~ /^[+-]inf$/;
1778
1779   # take lower of the two e's and adapt m1 to it to match m2
1780   my $e = $z->{_e};
1781   $e = $MBI->_zero() if !defined $e;            # if no BFLOAT?
1782   $e = $MBI->_copy($e);                         # make copy (didn't do it yet)
1783
1784   my $es;
1785
1786   ($e,$es) = _e_sub($e, $x->{_e}, $z->{_es} || '+', $x->{_es});
1787
1788   my $add = $MBI->_copy($z->{_m});
1789
1790   if ($es eq '-')                               # < 0
1791     {
1792     $MBI->_lsft( $x->{_m}, $e, 10);
1793     ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
1794     }
1795   elsif (!$MBI->_is_zero($e))                   # > 0
1796     {
1797     $MBI->_lsft($add, $e, 10);
1798     }
1799   # else: both e are the same, so just leave them
1800
1801   if ($x->{sign} eq $z->{sign})
1802     {
1803     # add
1804     $x->{_m} = $MBI->_add($x->{_m}, $add);
1805     }
1806   else
1807     {
1808     ($x->{_m}, $x->{sign}) = 
1809      _e_add($x->{_m}, $add, $x->{sign}, $z->{sign});
1810     }
1811
1812   # delete trailing zeros, then round
1813   $x->bnorm()->round(@r);
1814   }
1815
1816 sub bdiv 
1817   {
1818   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return 
1819     # (BFLOAT, BFLOAT) (quo, rem) or BFLOAT (only quo)
1820
1821   # set up parameters
1822   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1823   # objectify is costly, so avoid it
1824   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1825     {
1826     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1827     }
1828
1829   return $x if $x->modify('bdiv');
1830
1831     my $wantarray = wantarray;          # call only once
1832
1833     # At least one argument is NaN. This is handled the same way as in
1834     # Math::BigInt -> bdiv().
1835
1836     if ($x -> is_nan() || $y -> is_nan()) {
1837         return $wantarray ? ($x -> bnan(), $self -> bnan()) : $x -> bnan();
1838     }
1839
1840     # Divide by zero and modulo zero. This is handled the same way as in
1841     # Math::BigInt -> bdiv(). See the comment in the code for Math::BigInt ->
1842     # bdiv() for further details.
1843
1844     if ($y -> is_zero()) {
1845         my ($quo, $rem);
1846         if ($wantarray) {
1847             $rem = $x -> copy();
1848         }
1849         if ($x -> is_zero()) {
1850             $quo = $x -> bnan();
1851         } else {
1852             $quo = $x -> binf($x -> {sign});
1853         }
1854         return $wantarray ? ($quo, $rem) : $quo;
1855     }
1856
1857     # Numerator (dividend) is +/-inf. This is handled the same way as in
1858     # Math::BigInt -> bdiv(). See the comment in the code for Math::BigInt ->
1859     # bdiv() for further details.
1860
1861     if ($x -> is_inf()) {
1862         my ($quo, $rem);
1863         $rem = $self -> bnan() if $wantarray;
1864         if ($y -> is_inf()) {
1865             $quo = $x -> bnan();
1866         } else {
1867             my $sign = $x -> bcmp(0) == $y -> bcmp(0) ? '+' : '-';
1868             $quo = $x -> binf($sign);
1869         }
1870         return $wantarray ? ($quo, $rem) : $quo;
1871     }
1872
1873   # Denominator (divisor) is +/-inf. This is handled the same way as in
1874   # Math::BigInt -> bdiv(), with one exception: In scalar context,
1875   # Math::BigFloat does true division (although rounded), not floored division
1876   # (F-division), so a finite number divided by +/-inf is always zero. See the
1877   # comment in the code for Math::BigInt -> bdiv() for further details.
1878
1879   if ($y -> is_inf()) {
1880       my ($quo, $rem);
1881       if ($wantarray) {
1882           if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) {
1883               $rem = $x -> copy();
1884               $quo = $x -> bzero();
1885           } else {
1886               $rem = $self -> binf($y -> {sign});
1887               $quo = $x -> bone('-');
1888           }
1889           return ($quo, $rem);
1890       } else {
1891         if ($y -> is_inf()) {
1892             if ($x -> is_nan() || $x -> is_inf()) {
1893                 return $x -> bnan();
1894             } else {
1895                 return $x -> bzero();
1896             }
1897         }
1898     }
1899     }
1900
1901   # At this point, both the numerator and denominator are finite numbers, and
1902   # the denominator (divisor) is non-zero.
1903
1904   # x == 0?
1905   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
1906
1907   # upgrade ?
1908   return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
1909
1910   # we need to limit the accuracy to protect against overflow
1911   my $fallback = 0;
1912   my (@params,$scale);
1913   ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1914
1915   return $x if $x->is_nan();            # error in _find_round_parameters?
1916
1917   # no rounding at all, so must use fallback
1918     if (scalar @params == 0)
1919     {
1920     # simulate old behaviour
1921     $params[0] = $self->div_scale();    # and round to it as accuracy
1922     $scale = $params[0]+4;              # at least four more for proper round
1923     $params[2] = $r;                    # round mode by caller or undef
1924     $fallback = 1;                      # to clear a/p afterwards
1925     } else {
1926     # the 4 below is empirical, and there might be cases where it is not
1927     # enough...
1928     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1929     }
1930
1931     my $rem;
1932     $rem = $self -> bzero() if wantarray;
1933
1934   $y = $self->new($y) unless $y->isa('Math::BigFloat');
1935
1936   my $lx = $MBI -> _len($x->{_m}); my $ly = $MBI -> _len($y->{_m});
1937   $scale = $lx if $lx > $scale;
1938   $scale = $ly if $ly > $scale;
1939   my $diff = $ly - $lx;
1940   $scale += $diff if $diff > 0;         # if lx << ly, but not if ly << lx!
1941
1942   # check that $y is not 1 nor -1 and cache the result:
1943   my $y_not_one = !($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m}));
1944
1945   # flipping the sign of $y will also flip the sign of $x for the special
1946   # case of $x->bsub($x); so we can catch it below:
1947   my $xsign = $x->{sign};
1948   $y->{sign} =~ tr/+-/-+/;
1949
1950   if ($xsign ne $x->{sign})
1951     {
1952     # special case of $x /= $x results in 1
1953     $x->bone();                 # "fixes" also sign of $y, since $x is $y
1954     }
1955   else
1956     {
1957     # correct $y's sign again
1958     $y->{sign} =~ tr/+-/-+/;
1959     # continue with normal div code:
1960
1961     # make copy of $x in case of list context for later remainder calculation
1962     if (wantarray && $y_not_one)
1963       {
1964       $rem = $x->copy();
1965       }
1966
1967     $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+'; 
1968
1969     # check for / +-1 ( +/- 1E0)
1970     if ($y_not_one)
1971       {
1972       # promote BigInts and it's subclasses (except when already a BigFloat)
1973       $y = $self->new($y) unless $y->isa('Math::BigFloat'); 
1974
1975       # calculate the result to $scale digits and then round it
1976       # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1977       $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1978       $MBI->_div ($x->{_m},$y->{_m});   # a/c
1979
1980       # correct exponent of $x
1981       ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1982       # correct for 10**scale
1983       ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1984       $x->bnorm();              # remove trailing 0's
1985       }
1986     } # end else $x != $y
1987
1988   # shortcut to not run through _find_round_parameters again
1989   if (defined $params[0])
1990     {
1991     delete $x->{_a};                            # clear before round
1992     $x->bround($params[0],$params[2]);          # then round accordingly
1993     }
1994   else
1995     {
1996     delete $x->{_p};                            # clear before round
1997     $x->bfround($params[1],$params[2]);         # then round accordingly
1998     }
1999   if ($fallback)
2000     {
2001     # clear a/p after round, since user did not request it
2002     delete $x->{_a}; delete $x->{_p};
2003     }
2004
2005   if (wantarray)
2006     {
2007     if ($y_not_one)
2008       {
2009       $x -> bfloor();
2010       $rem->bmod($y,@params);                   # copy already done
2011       }
2012     if ($fallback)
2013       {
2014       # clear a/p after round, since user did not request it
2015       delete $rem->{_a}; delete $rem->{_p};
2016       }
2017     return ($x,$rem);
2018     }
2019   $x;
2020   }
2021
2022 sub bmod 
2023   {
2024   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return remainder
2025
2026   # set up parameters
2027   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2028   # objectify is costly, so avoid it
2029   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2030     {
2031     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2032     }
2033
2034   return $x if $x->modify('bmod');
2035
2036     # At least one argument is NaN. This is handled the same way as in
2037     # Math::BigInt -> bmod().
2038
2039     if ($x -> is_nan() || $y -> is_nan()) {
2040         return $x -> bnan();
2041     } 
2042
2043     # Modulo zero. This is handled the same way as in Math::BigInt -> bmod().
2044
2045     if ($y -> is_zero()) {
2046     return $x;
2047     }
2048
2049     # Numerator (dividend) is +/-inf. This is handled the same way as in
2050     # Math::BigInt -> bmod().
2051
2052     if ($x -> is_inf()) {
2053         return $x -> bnan();
2054     }
2055
2056     # Denominator (divisor) is +/-inf. This is handled the same way as in
2057     # Math::BigInt -> bmod().
2058
2059     if ($y -> is_inf()) {
2060         if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) {
2061             return $x;
2062         } else {
2063             return $x -> binf($y -> sign());
2064         }
2065     }
2066
2067   return $x->bzero() if $x->is_zero()
2068  || ($x->is_int() &&
2069   # check that $y == +1 or $y == -1:
2070     ($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m})));
2071
2072   my $cmp = $x->bacmp($y);                      # equal or $x < $y?
2073     if ($cmp == 0) {                            # $x == $y => result 0
2074         return $x -> bzero($a, $p);
2075     }
2076
2077   # only $y of the operands negative? 
2078     my $neg = $x->{sign} ne $y->{sign} ? 1 : 0;
2079
2080   $x->{sign} = $y->{sign};                              # calc sign first
2081     if ($cmp < 0 && $neg == 0) {                # $x < $y => result $x
2082         return $x -> round($a, $p, $r);
2083     }
2084   
2085   my $ym = $MBI->_copy($y->{_m});
2086   
2087   # 2e1 => 20
2088   $MBI->_lsft( $ym, $y->{_e}, 10) 
2089    if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
2090  
2091   # if $y has digits after dot
2092   my $shifty = 0;                       # correct _e of $x by this
2093   if ($y->{_es} eq '-')                 # has digits after dot
2094     {
2095     # 123 % 2.5 => 1230 % 25 => 5 => 0.5
2096     $shifty = $MBI->_num($y->{_e});     # no more digits after dot
2097     $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
2098     }
2099   # $ym is now mantissa of $y based on exponent 0
2100
2101   my $shiftx = 0;                       # correct _e of $x by this
2102   if ($x->{_es} eq '-')                 # has digits after dot
2103     {
2104     # 123.4 % 20 => 1234 % 200
2105     $shiftx = $MBI->_num($x->{_e});     # no more digits after dot
2106     $MBI->_lsft($ym, $x->{_e}, 10);     # 123 => 1230
2107     }
2108   # 123e1 % 20 => 1230 % 20
2109   if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
2110     {
2111     $MBI->_lsft( $x->{_m}, $x->{_e},10);        # es => '+' here
2112     }
2113
2114   $x->{_e} = $MBI->_new($shiftx);
2115   $x->{_es} = '+'; 
2116   $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
2117   $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
2118   
2119   # now mantissas are equalized, exponent of $x is adjusted, so calc result
2120
2121   $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
2122
2123   $x->{sign} = '+' if $MBI->_is_zero($x->{_m});         # fix sign for -0
2124   $x->bnorm();
2125
2126     if ($neg != 0 && ! $x -> is_zero())   # one of them negative => correct in place
2127     {
2128     my $r = $y - $x;
2129     $x->{_m} = $r->{_m};
2130     $x->{_e} = $r->{_e};
2131     $x->{_es} = $r->{_es};
2132     $x->{sign} = '+' if $MBI->_is_zero($x->{_m});       # fix sign for -0
2133     $x->bnorm();
2134     }
2135
2136   $x->round($a,$p,$r,$y);       # round and return
2137   }
2138
2139 sub broot
2140   {
2141   # calculate $y'th root of $x
2142   
2143   # set up parameters
2144   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2145   # objectify is costly, so avoid it
2146   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2147     {
2148     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2149     }
2150
2151   return $x if $x->modify('broot');
2152
2153   # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
2154   return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
2155          $y->{sign} !~ /^\+$/;
2156
2157   return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
2158   
2159   # we need to limit the accuracy to protect against overflow
2160   my $fallback = 0;
2161   my (@params,$scale);
2162   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2163
2164   return $x if $x->is_nan();            # error in _find_round_parameters?
2165
2166   # no rounding at all, so must use fallback
2167   if (scalar @params == 0) 
2168     {
2169     # simulate old behaviour
2170     $params[0] = $self->div_scale();    # and round to it as accuracy
2171     $scale = $params[0]+4;              # at least four more for proper round
2172     $params[2] = $r;                    # round mode by caller or undef
2173     $fallback = 1;                      # to clear a/p afterwards
2174     }
2175   else
2176     {
2177     # the 4 below is empirical, and there might be cases where it is not
2178     # enough...
2179     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2180     }
2181
2182   # when user set globals, they would interfere with our calculation, so
2183   # disable them and later re-enable them
2184   no strict 'refs';
2185   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2186   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2187   # we also need to disable any set A or P on $x (_find_round_parameters took
2188   # them already into account), since these would interfere, too
2189   delete $x->{_a}; delete $x->{_p};
2190   # need to disable $upgrade in BigInt, to avoid deep recursion
2191   local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
2192
2193   # remember sign and make $x positive, since -4 ** (1/2) => -2
2194   my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
2195
2196   my $is_two = 0;
2197   if ($y->isa('Math::BigFloat'))
2198     {
2199     $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
2200     }
2201   else
2202     {
2203     $is_two = ($y == 2);
2204     }
2205
2206   # normal square root if $y == 2:
2207   if ($is_two)
2208     {
2209     $x->bsqrt($scale+4);
2210     }
2211   elsif ($y->is_one('-'))
2212     {
2213     # $x ** -1 => 1/$x
2214     my $u = $self->bone()->bdiv($x,$scale);
2215     # copy private parts over
2216     $x->{_m} = $u->{_m};
2217     $x->{_e} = $u->{_e};
2218     $x->{_es} = $u->{_es};
2219     }
2220   else
2221     {
2222     # calculate the broot() as integer result first, and if it fits, return
2223     # it rightaway (but only if $x and $y are integer):
2224
2225     my $done = 0;                               # not yet
2226     if ($y->is_int() && $x->is_int())
2227       {
2228       my $i = $MBI->_copy( $x->{_m} );
2229       $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2230       my $int = Math::BigInt->bzero();
2231       $int->{value} = $i;
2232       $int->broot($y->as_number());
2233       # if ($exact)
2234       if ($int->copy()->bpow($y) == $x)
2235         {
2236         # found result, return it
2237         $x->{_m} = $int->{value};
2238         $x->{_e} = $MBI->_zero();
2239         $x->{_es} = '+';
2240         $x->bnorm();
2241         $done = 1;
2242         }
2243       }
2244     if ($done == 0)
2245       {
2246       my $u = $self->bone()->bdiv($y,$scale+4);
2247       delete $u->{_a}; delete $u->{_p};         # otherwise it conflicts
2248       $x->bpow($u,$scale+4);                    # el cheapo
2249       }
2250     }
2251   $x->bneg() if $sign == 1;
2252   
2253   # shortcut to not run through _find_round_parameters again
2254   if (defined $params[0])
2255     {
2256     $x->bround($params[0],$params[2]);          # then round accordingly
2257     }
2258   else
2259     {
2260     $x->bfround($params[1],$params[2]);         # then round accordingly
2261     }
2262   if ($fallback)
2263     {
2264     # clear a/p after round, since user did not request it
2265     delete $x->{_a}; delete $x->{_p};
2266     }
2267   # restore globals
2268   $$abr = $ab; $$pbr = $pb;
2269   $x;
2270   }
2271
2272 sub bsqrt
2273   { 
2274   # calculate square root
2275   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2276
2277   return $x if $x->modify('bsqrt');
2278
2279   return $x->bnan() if $x->{sign} !~ /^[+]/;    # NaN, -inf or < 0
2280   return $x if $x->{sign} eq '+inf';            # sqrt(inf) == inf
2281   return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
2282
2283   # we need to limit the accuracy to protect against overflow
2284   my $fallback = 0;
2285   my (@params,$scale);
2286   ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2287
2288   return $x if $x->is_nan();            # error in _find_round_parameters?
2289
2290   # no rounding at all, so must use fallback
2291   if (scalar @params == 0) 
2292     {
2293     # simulate old behaviour
2294     $params[0] = $self->div_scale();    # and round to it as accuracy
2295     $scale = $params[0]+4;              # at least four more for proper round
2296     $params[2] = $r;                    # round mode by caller or undef
2297     $fallback = 1;                      # to clear a/p afterwards
2298     }
2299   else
2300     {
2301     # the 4 below is empirical, and there might be cases where it is not
2302     # enough...
2303     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2304     }
2305
2306   # when user set globals, they would interfere with our calculation, so
2307   # disable them and later re-enable them
2308   no strict 'refs';
2309   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2310   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2311   # we also need to disable any set A or P on $x (_find_round_parameters took
2312   # them already into account), since these would interfere, too
2313   delete $x->{_a}; delete $x->{_p};
2314   # need to disable $upgrade in BigInt, to avoid deep recursion
2315   local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
2316
2317   my $i = $MBI->_copy( $x->{_m} );
2318   $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2319   my $xas = Math::BigInt->bzero();
2320   $xas->{value} = $i;
2321
2322   my $gs = $xas->copy()->bsqrt();       # some guess
2323
2324   if (($x->{_es} ne '-')                # guess can't be accurate if there are
2325                                         # digits after the dot
2326    && ($xas->bacmp($gs * $gs) == 0))    # guess hit the nail on the head?
2327     {
2328     # exact result, copy result over to keep $x
2329     $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
2330     $x->bnorm();
2331     # shortcut to not run through _find_round_parameters again
2332     if (defined $params[0])
2333       {
2334       $x->bround($params[0],$params[2]);        # then round accordingly
2335       }
2336     else
2337       {
2338       $x->bfround($params[1],$params[2]);       # then round accordingly
2339       }
2340     if ($fallback)
2341       {
2342       # clear a/p after round, since user did not request it
2343       delete $x->{_a}; delete $x->{_p};
2344       }
2345     # re-enable A and P, upgrade is taken care of by "local"
2346     ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
2347     return $x;
2348     }
2349  
2350   # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
2351   # of the result by multiplying the input by 100 and then divide the integer
2352   # result of sqrt(input) by 10. Rounding afterwards returns the real result.
2353
2354   # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
2355   my $y1 = $MBI->_copy($x->{_m});
2356
2357   my $length = $MBI->_len($y1);
2358   
2359   # Now calculate how many digits the result of sqrt(y1) would have
2360   my $digits = int($length / 2);
2361
2362   # But we need at least $scale digits, so calculate how many are missing
2363   my $shift = $scale - $digits;
2364
2365   # This happens if the input had enough digits
2366   # (we take care of integer guesses above)
2367   $shift = 0 if $shift < 0; 
2368
2369   # Multiply in steps of 100, by shifting left two times the "missing" digits
2370   my $s2 = $shift * 2;
2371
2372   # We now make sure that $y1 has the same odd or even number of digits than
2373   # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
2374   # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
2375   # steps of 10. The length of $x does not count, since an even or odd number
2376   # of digits before the dot is not changed by adding an even number of digits
2377   # after the dot (the result is still odd or even digits long).
2378   $s2++ if $MBI->_is_odd($x->{_e});
2379
2380   $MBI->_lsft( $y1, $MBI->_new($s2), 10);
2381
2382   # now take the square root and truncate to integer
2383   $y1 = $MBI->_sqrt($y1);
2384
2385   # By "shifting" $y1 right (by creating a negative _e) we calculate the final
2386   # result, which is than later rounded to the desired scale.
2387
2388   # calculate how many zeros $x had after the '.' (or before it, depending
2389   # on sign of $dat, the result should have half as many:
2390   my $dat = $MBI->_num($x->{_e});
2391   $dat = -$dat if $x->{_es} eq '-';
2392   $dat += $length;
2393
2394   if ($dat > 0)
2395     {
2396     # no zeros after the dot (e.g. 1.23, 0.49 etc)
2397     # preserve half as many digits before the dot than the input had 
2398     # (but round this "up")
2399     $dat = int(($dat+1)/2);
2400     }
2401   else
2402     {
2403     $dat = int(($dat)/2);
2404     }
2405   $dat -= $MBI->_len($y1);
2406   if ($dat < 0)
2407     {
2408     $dat = abs($dat);
2409     $x->{_e} = $MBI->_new( $dat );
2410     $x->{_es} = '-';
2411     }
2412   else
2413     {    
2414     $x->{_e} = $MBI->_new( $dat );
2415     $x->{_es} = '+';
2416     }
2417   $x->{_m} = $y1;
2418   $x->bnorm();
2419
2420   # shortcut to not run through _find_round_parameters again
2421   if (defined $params[0])
2422     {
2423     $x->bround($params[0],$params[2]);          # then round accordingly
2424     }
2425   else
2426     {
2427     $x->bfround($params[1],$params[2]);         # then round accordingly
2428     }
2429   if ($fallback)
2430     {
2431     # clear a/p after round, since user did not request it
2432     delete $x->{_a}; delete $x->{_p};
2433     }
2434   # restore globals
2435   $$abr = $ab; $$pbr = $pb;
2436   $x;
2437   }
2438
2439 sub bfac
2440   {
2441   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2442   # compute factorial number, modifies first argument
2443
2444   # set up parameters
2445   my ($self,$x,@r) = (ref($_[0]),@_);
2446   # objectify is costly, so avoid it
2447   ($self,$x,@r) = objectify(1,@_) if !ref($x);
2448
2449   # inf => inf
2450   return $x if $x->modify('bfac') || $x->{sign} eq '+inf';      
2451
2452   return $x->bnan() 
2453     if (($x->{sign} ne '+') ||          # inf, NaN, <0 etc => NaN
2454      ($x->{_es} ne '+'));               # digits after dot?
2455
2456   # use BigInt's bfac() for faster calc
2457   if (! $MBI->_is_zero($x->{_e}))
2458     {
2459     $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
2460     $x->{_e} = $MBI->_zero();           # normalize
2461     $x->{_es} = '+';
2462     }
2463   $MBI->_fac($x->{_m});                 # calculate factorial
2464   $x->bnorm()->round(@r);               # norm again and round result
2465   }
2466
2467 sub _pow
2468   {
2469   # Calculate a power where $y is a non-integer, like 2 ** 0.3
2470   my ($x,$y,@r) = @_;
2471   my $self = ref($x);
2472
2473   # if $y == 0.5, it is sqrt($x)
2474   $HALF = $self->new($HALF) unless ref($HALF);
2475   return $x->bsqrt(@r,$y) if $y->bcmp($HALF) == 0;
2476
2477   # Using:
2478   # a ** x == e ** (x * ln a)
2479
2480   # u = y * ln x
2481   #                _                         _
2482   # Taylor:       |   u    u^2    u^3         |
2483   # x ** y  = 1 + |  --- + --- + ----- + ...  |
2484   #               |_  1    1*2   1*2*3       _|
2485
2486   # we need to limit the accuracy to protect against overflow
2487   my $fallback = 0;
2488   my ($scale,@params);
2489   ($x,@params) = $x->_find_round_parameters(@r);
2490     
2491   return $x if $x->is_nan();            # error in _find_round_parameters?
2492
2493   # no rounding at all, so must use fallback
2494   if (scalar @params == 0)
2495     {
2496     # simulate old behaviour
2497     $params[0] = $self->div_scale();    # and round to it as accuracy
2498     $params[1] = undef;                 # disable P
2499     $scale = $params[0]+4;              # at least four more for proper round
2500     $params[2] = $r[2];                 # round mode by caller or undef
2501     $fallback = 1;                      # to clear a/p afterwards
2502     }
2503   else
2504     {
2505     # the 4 below is empirical, and there might be cases where it is not
2506     # enough...
2507     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2508     }
2509
2510   # when user set globals, they would interfere with our calculation, so
2511   # disable them and later re-enable them
2512   no strict 'refs';
2513   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2514   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2515   # we also need to disable any set A or P on $x (_find_round_parameters took
2516   # them already into account), since these would interfere, too
2517   delete $x->{_a}; delete $x->{_p};
2518   # need to disable $upgrade in BigInt, to avoid deep recursion
2519   local $Math::BigInt::upgrade = undef;
2520  
2521   my ($limit,$v,$u,$below,$factor,$next,$over);
2522
2523   $u = $x->copy()->blog(undef,$scale)->bmul($y);
2524   $v = $self->bone();                           # 1
2525   $factor = $self->new(2);                      # 2
2526   $x->bone();                                   # first term: 1
2527
2528   $below = $v->copy();
2529   $over = $u->copy();
2530
2531   $limit = $self->new("1E-". ($scale-1));
2532   #my $steps = 0;
2533   while (3 < 5)
2534     {
2535     # we calculate the next term, and add it to the last
2536     # when the next term is below our limit, it won't affect the outcome
2537     # anymore, so we stop:
2538     $next = $over->copy()->bdiv($below,$scale);
2539     last if $next->bacmp($limit) <= 0;
2540     $x->badd($next);
2541     # calculate things for the next term
2542     $over *= $u; $below *= $factor; $factor->binc();
2543
2544     last if $x->{sign} !~ /^[-+]$/;
2545
2546     #$steps++;
2547     }
2548   
2549   # shortcut to not run through _find_round_parameters again
2550   if (defined $params[0])
2551     {
2552     $x->bround($params[0],$params[2]);          # then round accordingly
2553     }
2554   else
2555     {
2556     $x->bfround($params[1],$params[2]);         # then round accordingly
2557     }
2558   if ($fallback)
2559     {
2560     # clear a/p after round, since user did not request it
2561     delete $x->{_a}; delete $x->{_p};
2562     }
2563   # restore globals
2564   $$abr = $ab; $$pbr = $pb;
2565   $x;
2566   }
2567
2568 sub bpow 
2569   {
2570   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2571   # compute power of two numbers, second arg is used as integer
2572   # modifies first argument
2573
2574   # set up parameters
2575   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2576   # objectify is costly, so avoid it
2577   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2578     {
2579     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2580     }
2581
2582   return $x if $x->modify('bpow');
2583
2584   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
2585   return $x if $x->{sign} =~ /^[+-]inf$/;
2586   
2587   # cache the result of is_zero
2588   my $y_is_zero = $y->is_zero();
2589   return $x->bone() if $y_is_zero;
2590   return $x         if $x->is_one() || $y->is_one();
2591
2592   my $x_is_zero = $x->is_zero();
2593   return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int();         # non-integer power
2594
2595   my $y1 = $y->as_number()->{value};                    # make MBI part
2596
2597   # if ($x == -1)
2598   if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
2599     {
2600     # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
2601     return $MBI->_is_odd($y1) ? $x : $x->babs(1);
2602     }
2603   if ($x_is_zero)
2604     {
2605     return $x if $y->{sign} eq '+';     # 0**y => 0 (if not y <= 0)
2606     # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
2607     return $x->binf();
2608     }
2609
2610   my $new_sign = '+';
2611   $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
2612
2613   # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
2614   $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
2615   $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
2616
2617   $x->{sign} = $new_sign;
2618   $x->bnorm();
2619   if ($y->{sign} eq '-')
2620     {
2621     # modify $x in place!
2622     my $z = $x->copy(); $x->bone();
2623     return scalar $x->bdiv($z,$a,$p,$r);        # round in one go (might ignore y's A!)
2624     }
2625   $x->round($a,$p,$r,$y);
2626   }
2627
2628 sub bmodpow
2629   {
2630   # takes a very large number to a very large exponent in a given very
2631   # large modulus, quickly, thanks to binary exponentiation. Supports
2632   # negative exponents.
2633   my ($self,$num,$exp,$mod,@r) = objectify(3,@_);
2634
2635   return $num if $num->modify('bmodpow');
2636
2637   # check modulus for valid values
2638   return $num->bnan() if ($mod->{sign} ne '+'           # NaN, - , -inf, +inf
2639                        || $mod->is_zero());
2640
2641   # check exponent for valid values
2642   if ($exp->{sign} =~ /\w/)
2643     {
2644     # i.e., if it's NaN, +inf, or -inf...
2645     return $num->bnan();
2646     }
2647
2648   $num->bmodinv ($mod) if ($exp->{sign} eq '-');
2649
2650   # check num for valid values (also NaN if there was no inverse but $exp < 0)
2651   return $num->bnan() if $num->{sign} !~ /^[+-]$/;
2652
2653   # $mod is positive, sign on $exp is ignored, result also positive
2654
2655   # XXX TODO: speed it up when all three numbers are integers
2656   $num->bpow($exp)->bmod($mod);
2657   }
2658
2659 ###############################################################################
2660 # trigonometric functions
2661
2662 # helper function for bpi() and batan2(), calculates arcus tanges (1/x)
2663
2664 sub _atan_inv
2665   {
2666   # return a/b so that a/b approximates atan(1/x) to at least limit digits
2667   my ($self, $x, $limit) = @_;
2668
2669   # Taylor:       x^3   x^5   x^7   x^9
2670   #    atan = x - --- + --- - --- + --- - ...
2671   #                3     5     7     9 
2672
2673   #               1      1         1        1
2674   #    atan 1/x = - - ------- + ------- - ------- + ...
2675   #               x   x^3 * 3   x^5 * 5   x^7 * 7 
2676
2677   #               1      1         1            1
2678   #    atan 1/x = - - --------- + ---------- - ----------- + ... 
2679   #               5    3 * 125     5 * 3125     7 * 78125
2680
2681   # Subtraction/addition of a rational:
2682
2683   #  5    7    5*3 +- 7*4
2684   #  - +- -  = ----------
2685   #  4    3       4*3
2686
2687   # Term:  N        N+1
2688   #
2689   #        a             1                  a * d * c +- b
2690   #        ----- +- ------------------  =  ----------------
2691   #        b           d * c                b * d * c
2692
2693   #  since b1 = b0 * (d-2) * c
2694
2695   #        a             1                  a * d +- b / c
2696   #        ----- +- ------------------  =  ----------------
2697   #        b           d * c                b * d 
2698
2699   # and  d = d + 2
2700   # and  c = c * x * x
2701
2702   #        u = d * c
2703   #        stop if length($u) > limit 
2704   #        a = a * u +- b
2705   #        b = b * u
2706   #        d = d + 2
2707   #        c = c * x * x
2708   #        sign = 1 - sign
2709
2710   my $a = $MBI->_one();
2711   my $b = $MBI->_copy($x);
2712  
2713   my $x2  = $MBI->_mul( $MBI->_copy($x), $b);           # x2 = x * x
2714   my $d   = $MBI->_new( 3 );                            # d = 3
2715   my $c   = $MBI->_mul( $MBI->_copy($x), $x2);          # c = x ^ 3
2716   my $two = $MBI->_new( 2 );
2717
2718   # run the first step unconditionally
2719   my $u = $MBI->_mul( $MBI->_copy($d), $c);
2720   $a = $MBI->_mul($a, $u);
2721   $a = $MBI->_sub($a, $b);
2722   $b = $MBI->_mul($b, $u);
2723   $d = $MBI->_add($d, $two);
2724   $c = $MBI->_mul($c, $x2);
2725
2726   # a is now a * (d-3) * c
2727   # b is now b * (d-2) * c
2728
2729   # run the second step unconditionally
2730   $u = $MBI->_mul( $MBI->_copy($d), $c);
2731   $a = $MBI->_mul($a, $u);
2732   $a = $MBI->_add($a, $b);
2733   $b = $MBI->_mul($b, $u);
2734   $d = $MBI->_add($d, $two);
2735   $c = $MBI->_mul($c, $x2);
2736
2737   # a is now a * (d-3) * (d-5) * c * c  
2738   # b is now b * (d-2) * (d-4) * c * c
2739
2740   # so we can remove c * c from both a and b to shorten the numbers involved:
2741   $a = $MBI->_div($a, $x2);
2742   $b = $MBI->_div($b, $x2);
2743   $a = $MBI->_div($a, $x2);
2744   $b = $MBI->_div($b, $x2);
2745
2746 #  my $step = 0; 
2747   my $sign = 0;                                         # 0 => -, 1 => +
2748   while (3 < 5)
2749     {
2750 #    $step++;
2751 #    if (($i++ % 100) == 0)
2752 #      {
2753 #    print "a=",$MBI->_str($a),"\n";
2754 #    print "b=",$MBI->_str($b),"\n";
2755 #      }
2756 #    print "d=",$MBI->_str($d),"\n";
2757 #    print "x2=",$MBI->_str($x2),"\n";
2758 #    print "c=",$MBI->_str($c),"\n";
2759
2760     my $u = $MBI->_mul( $MBI->_copy($d), $c);
2761     # use _alen() for libs like GMP where _len() would be O(N^2)
2762     last if $MBI->_alen($u) > $limit;
2763     my ($bc,$r) = $MBI->_div( $MBI->_copy($b), $c);
2764     if ($MBI->_is_zero($r))
2765       {
2766       # b / c is an integer, so we can remove c from all terms
2767       # this happens almost every time:
2768       $a = $MBI->_mul($a, $d);
2769       $a = $MBI->_sub($a, $bc) if $sign == 0;
2770       $a = $MBI->_add($a, $bc) if $sign == 1;
2771       $b = $MBI->_mul($b, $d);
2772       }
2773     else
2774       {
2775       # b / c is not an integer, so we keep c in the terms
2776       # this happens very rarely, for instance for x = 5, this happens only
2777       # at the following steps:
2778       # 1, 5, 14, 32, 72, 157, 340, ...
2779       $a = $MBI->_mul($a, $u);
2780       $a = $MBI->_sub($a, $b) if $sign == 0;
2781       $a = $MBI->_add($a, $b) if $sign == 1;
2782       $b = $MBI->_mul($b, $u);
2783       }
2784     $d = $MBI->_add($d, $two);
2785     $c = $MBI->_mul($c, $x2);
2786     $sign = 1 - $sign;
2787
2788     }
2789
2790 #  print "Took $step steps for ", $MBI->_str($x),"\n";
2791 #  print "a=",$MBI->_str($a),"\n"; print "b=",$MBI->_str($b),"\n";
2792   # return a/b so that a/b approximates atan(1/x)
2793   ($a,$b);
2794   }
2795
2796 sub bpi
2797   {
2798   my ($self,$n) = @_;
2799   if (@_ == 0)
2800     {
2801     $self = $class;
2802     }
2803   if (@_ == 1)
2804     {
2805     # called like Math::BigFloat::bpi(10);
2806     $n = $self; $self = $class;
2807     # called like Math::BigFloat->bpi();
2808     $n = undef if $n eq 'Math::BigFloat';
2809     }
2810   $self = ref($self) if ref($self);
2811   my $fallback = defined $n ? 0 : 1;
2812   $n = 40 if !defined $n || $n < 1;
2813
2814   # after 黃見利 (Hwang Chien-Lih) (1997)
2815   # pi/4 = 183 * atan(1/239) + 32 * atan(1/1023) – 68 * atan(1/5832)
2816   #      + 12 * atan(1/110443) - 12 * atan(1/4841182) - 100 * atan(1/6826318)
2817
2818   # a few more to prevent rounding errors
2819   $n += 4;
2820
2821   my ($a,$b) = $self->_atan_inv( $MBI->_new(239),$n);
2822   my ($c,$d) = $self->_atan_inv( $MBI->_new(1023),$n);
2823   my ($e,$f) = $self->_atan_inv( $MBI->_new(5832),$n);
2824   my ($g,$h) = $self->_atan_inv( $MBI->_new(110443),$n);
2825   my ($i,$j) = $self->_atan_inv( $MBI->_new(4841182),$n);
2826   my ($k,$l) = $self->_atan_inv( $MBI->_new(6826318),$n);
2827
2828   $MBI->_mul($a, $MBI->_new(732));
2829   $MBI->_mul($c, $MBI->_new(128));
2830   $MBI->_mul($e, $MBI->_new(272));
2831   $MBI->_mul($g, $MBI->_new(48));
2832   $MBI->_mul($i, $MBI->_new(48));
2833   $MBI->_mul($k, $MBI->_new(400));
2834
2835   my $x = $self->bone(); $x->{_m} = $a; my $x_d = $self->bone(); $x_d->{_m} = $b;
2836   my $y = $self->bone(); $y->{_m} = $c; my $y_d = $self->bone(); $y_d->{_m} = $d;
2837   my $z = $self->bone(); $z->{_m} = $e; my $z_d = $self->bone(); $z_d->{_m} = $f;
2838   my $u = $self->bone(); $u->{_m} = $g; my $u_d = $self->bone(); $u_d->{_m} = $h;
2839   my $v = $self->bone(); $v->{_m} = $i; my $v_d = $self->bone(); $v_d->{_m} = $j;
2840   my $w = $self->bone(); $w->{_m} = $k; my $w_d = $self->bone(); $w_d->{_m} = $l;
2841   $x->bdiv($x_d, $n);
2842   $y->bdiv($y_d, $n);
2843   $z->bdiv($z_d, $n);
2844   $u->bdiv($u_d, $n);
2845   $v->bdiv($v_d, $n);
2846   $w->bdiv($w_d, $n);
2847
2848   delete $x->{_a}; delete $y->{_a}; delete $z->{_a};
2849   delete $u->{_a}; delete $v->{_a}; delete $w->{_a};
2850   $x->badd($y)->bsub($z)->badd($u)->bsub($v)->bsub($w);
2851
2852   $x->bround($n-4);
2853   delete $x->{_a} if $fallback == 1;
2854   $x;
2855   }
2856
2857 sub bcos
2858   {
2859   # Calculate a cosinus of x.
2860   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2861
2862   # Taylor:      x^2   x^4   x^6   x^8
2863   #    cos = 1 - --- + --- - --- + --- ...
2864   #               2!    4!    6!    8!
2865
2866   # we need to limit the accuracy to protect against overflow
2867   my $fallback = 0;
2868   my ($scale,@params);
2869   ($x,@params) = $x->_find_round_parameters(@r);
2870     
2871   #         constant object       or error in _find_round_parameters?
2872   return $x if $x->modify('bcos') || $x->is_nan();
2873
2874   return $x->bone(@r) if $x->is_zero();
2875
2876   # no rounding at all, so must use fallback
2877   if (scalar @params == 0)
2878     {
2879     # simulate old behaviour
2880     $params[0] = $self->div_scale();    # and round to it as accuracy
2881     $params[1] = undef;                 # disable P
2882     $scale = $params[0]+4;              # at least four more for proper round
2883     $params[2] = $r[2];                 # round mode by caller or undef
2884     $fallback = 1;                      # to clear a/p afterwards
2885     }
2886   else
2887     {
2888     # the 4 below is empirical, and there might be cases where it is not
2889     # enough...
2890     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2891     }
2892
2893   # when user set globals, they would interfere with our calculation, so
2894   # disable them and later re-enable them
2895   no strict 'refs';
2896   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2897   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2898   # we also need to disable any set A or P on $x (_find_round_parameters took
2899   # them already into account), since these would interfere, too
2900   delete $x->{_a}; delete $x->{_p};
2901   # need to disable $upgrade in BigInt, to avoid deep recursion
2902   local $Math::BigInt::upgrade = undef;
2903  
2904   my $last = 0;
2905   my $over = $x * $x;                   # X ^ 2
2906   my $x2 = $over->copy();               # X ^ 2; difference between terms
2907   my $sign = 1;                         # start with -=
2908   my $below = $self->new(2); my $factorial = $self->new(3);
2909   $x->bone(); delete $x->{_a}; delete $x->{_p};
2910
2911   my $limit = $self->new("1E-". ($scale-1));
2912   #my $steps = 0;
2913   while (3 < 5)
2914     {
2915     # we calculate the next term, and add it to the last
2916     # when the next term is below our limit, it won't affect the outcome
2917     # anymore, so we stop:
2918     my $next = $over->copy()->bdiv($below,$scale);
2919     last if $next->bacmp($limit) <= 0;
2920
2921     if ($sign == 0)
2922       {
2923       $x->badd($next);
2924       }
2925     else
2926       {
2927       $x->bsub($next);
2928       }
2929     $sign = 1-$sign;                                    # alternate
2930     # calculate things for the next term
2931     $over->bmul($x2);                                   # $x*$x
2932     $below->bmul($factorial); $factorial->binc();       # n*(n+1)
2933     $below->bmul($factorial); $factorial->binc();       # n*(n+1)
2934     }
2935
2936   # shortcut to not run through _find_round_parameters again
2937   if (defined $params[0])
2938     {
2939     $x->bround($params[0],$params[2]);          # then round accordingly
2940     }
2941   else
2942     {
2943     $x->bfround($params[1],$params[2]);         # then round accordingly
2944     }
2945   if ($fallback)
2946     {
2947     # clear a/p after round, since user did not request it
2948     delete $x->{_a}; delete $x->{_p};
2949     }
2950   # restore globals
2951   $$abr = $ab; $$pbr = $pb;
2952   $x;
2953   }
2954
2955 sub bsin
2956   {
2957   # Calculate a sinus of x.
2958   my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2959
2960   # taylor:      x^3   x^5   x^7   x^9
2961   #    sin = x - --- + --- - --- + --- ...
2962   #               3!    5!    7!    9!
2963
2964   # we need to limit the accuracy to protect against overflow
2965   my $fallback = 0;
2966   my ($scale,@params);
2967   ($x,@params) = $x->_find_round_parameters(@r);
2968     
2969   #         constant object       or error in _find_round_parameters?
2970   return $x if $x->modify('bsin') || $x->is_nan();
2971
2972   return $x->bzero(@r) if $x->is_zero();
2973
2974   # no rounding at all, so must use fallback
2975   if (scalar @params == 0)
2976     {
2977     # simulate old behaviour
2978     $params[0] = $self->div_scale();    # and round to it as accuracy
2979     $params[1] = undef;                 # disable P
2980     $scale = $params[0]+4;              # at least four more for proper round
2981     $params[2] = $r[2];                 # round mode by caller or undef
2982     $fallback = 1;                      # to clear a/p afterwards
2983     }
2984   else
2985     {
2986     # the 4 below is empirical, and there might be cases where it is not
2987     # enough...
2988     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2989     }
2990
2991   # when user set globals, they would interfere with our calculation, so
2992   # disable them and later re-enable them
2993   no strict 'refs';
2994   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2995   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2996   # we also need to disable any set A or P on $x (_find_round_parameters took
2997   # them already into account), since these would interfere, too
2998   delete $x->{_a}; delete $x->{_p};
2999   # need to disable $upgrade in BigInt, to avoid deep recursion
3000   local $Math::BigInt::upgrade = undef;
3001  
3002   my $last = 0;
3003   my $over = $x * $x;                   # X ^ 2
3004   my $x2 = $over->copy();               # X ^ 2; difference between terms
3005   $over->bmul($x);                      # X ^ 3 as starting value
3006   my $sign = 1;                         # start with -=
3007   my $below = $self->new(6); my $factorial = $self->new(4);
3008   delete $x->{_a}; delete $x->{_p};
3009
3010   my $limit = $self->new("1E-". ($scale-1));
3011   #my $steps = 0;
3012   while (3 < 5)
3013     {
3014     # we calculate the next term, and add it to the last
3015     # when the next term is below our limit, it won't affect the outcome
3016     # anymore, so we stop:
3017     my $next = $over->copy()->bdiv($below,$scale);
3018     last if $next->bacmp($limit) <= 0;
3019
3020     if ($sign == 0)
3021       {
3022       $x->badd($next);
3023       }
3024     else
3025       {
3026       $x->bsub($next);
3027       }
3028     $sign = 1-$sign;                                    # alternate
3029     # calculate things for the next term
3030     $over->bmul($x2);                                   # $x*$x
3031     $below->bmul($factorial); $factorial->binc();       # n*(n+1)
3032     $below->bmul($factorial); $factorial->binc();       # n*(n+1)
3033     }
3034
3035   # shortcut to not run through _find_round_parameters again
3036   if (defined $params[0])
3037     {
3038     $x->bround($params[0],$params[2]);          # then round accordingly
3039     }
3040   else
3041     {
3042     $x->bfround($params[1],$params[2]);         # then round accordingly
3043     }
3044   if ($fallback)
3045     {
3046     # clear a/p after round, since user did not request it
3047     delete $x->{_a}; delete $x->{_p};
3048     }
3049   # restore globals
3050   $$abr = $ab; $$pbr = $pb;
3051   $x;
3052   }
3053
3054 sub batan2
3055   { 
3056   # calculate arcus tangens of ($y/$x)
3057   
3058   # set up parameters
3059   my ($self,$y,$x,@r) = (ref($_[0]),@_);
3060   # objectify is costly, so avoid it
3061   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3062     {
3063     ($self,$y,$x,@r) = objectify(2,@_);
3064     }
3065
3066   return $y if $y->modify('batan2');
3067
3068   return $y->bnan() if ($y->{sign} eq $nan) || ($x->{sign} eq $nan);
3069
3070   # Y X
3071   # 0 0 result is 0
3072   # 0 +x result is 0
3073   # ? inf result is 0
3074   return $y->bzero(@r) if ($x->is_inf('+') && !$y->is_inf()) || ($y->is_zero() && $x->{sign} eq '+');
3075
3076   # Y    X
3077   # != 0 -inf result is +- pi
3078   if ($x->is_inf() || $y->is_inf())
3079     {
3080     # calculate PI
3081     my $pi = $self->bpi(@r);
3082     if ($y->is_inf())
3083       {
3084       # upgrade to BigRat etc. 
3085       return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
3086       if ($x->{sign} eq '-inf')
3087         {
3088         # calculate 3 pi/4
3089         $MBI->_mul($pi->{_m}, $MBI->_new(3));
3090         $MBI->_div($pi->{_m}, $MBI->_new(4));
3091         }
3092       elsif ($x->{sign} eq '+inf')
3093         {
3094         # calculate pi/4
3095         $MBI->_div($pi->{_m}, $MBI->_new(4));
3096         }
3097       else
3098         {
3099         # calculate pi/2
3100         $MBI->_div($pi->{_m}, $MBI->_new(2));
3101         }
3102       $y->{sign} = substr($y->{sign},0,1); # keep +/-
3103       }
3104     # modify $y in place
3105     $y->{_m} = $pi->{_m};
3106     $y->{_e} = $pi->{_e};
3107     $y->{_es} = $pi->{_es};
3108     # keep the sign of $y
3109     return $y;
3110     }
3111
3112   return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
3113
3114   # Y X
3115   # 0 -x result is PI
3116   if ($y->is_zero())
3117     {
3118     # calculate PI
3119     my $pi = $self->bpi(@r);
3120     # modify $y in place
3121     $y->{_m} = $pi->{_m};
3122     $y->{_e} = $pi->{_e};
3123     $y->{_es} = $pi->{_es};
3124     $y->{sign} = '+';
3125     return $y;
3126     }
3127
3128   # Y X
3129   # +y 0 result is PI/2
3130   # -y 0 result is -PI/2
3131   if ($x->is_zero())
3132     {
3133     # calculate PI/2
3134     my $pi = $self->bpi(@r);
3135     # modify $y in place
3136     $y->{_m} = $pi->{_m};
3137     $y->{_e} = $pi->{_e};
3138     $y->{_es} = $pi->{_es};
3139     # -y => -PI/2, +y => PI/2
3140     $MBI->_div($y->{_m}, $MBI->_new(2));
3141     return $y;
3142     }
3143
3144   # we need to limit the accuracy to protect against overflow
3145   my $fallback = 0;
3146   my ($scale,@params);
3147   ($y,@params) = $y->_find_round_parameters(@r);
3148     
3149   # error in _find_round_parameters?
3150   return $y if $y->is_nan();
3151
3152   # no rounding at all, so must use fallback
3153   if (scalar @params == 0)
3154     {
3155     # simulate old behaviour
3156     $params[0] = $self->div_scale();    # and round to it as accuracy
3157     $params[1] = undef;                 # disable P
3158     $scale = $params[0]+4;              # at least four more for proper round
3159     $params[2] = $r[2];                 # round mode by caller or undef
3160     $fallback = 1;                      # to clear a/p afterwards
3161     }
3162   else
3163     {
3164     # the 4 below is empirical, and there might be cases where it is not
3165     # enough...
3166     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3167     }
3168
3169   # inlined is_one() && is_one('-')
3170   if ($MBI->_is_one($y->{_m}) && $MBI->_is_zero($y->{_e}))
3171     {
3172     # shortcut: 1 1 result is PI/4
3173     # inlined is_one() && is_one('-')
3174     if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
3175       {
3176       # 1,1 => PI/4
3177       my $pi_4 = $self->bpi( $scale - 3);
3178       # modify $y in place
3179       $y->{_m} = $pi_4->{_m};
3180       $y->{_e} = $pi_4->{_e};
3181       $y->{_es} = $pi_4->{_es};
3182       # 1 1 => +
3183       # -1 1 => -
3184       # 1 -1 => -
3185       # -1 -1 => +
3186       $y->{sign} = $x->{sign} eq $y->{sign} ? '+' : '-';
3187       $MBI->_div($y->{_m}, $MBI->_new(4));
3188       return $y;
3189       }
3190     # shortcut: 1 int(X) result is _atan_inv(X)
3191
3192     # is integer
3193     if ($x->{_es} eq '+')
3194       {
3195       my $x1 = $MBI->_copy($x->{_m});
3196       $MBI->_lsft($x1, $x->{_e},10) unless $MBI->_is_zero($x->{_e});
3197
3198       my ($a,$b) = $self->_atan_inv($x1, $scale);
3199       my $y_sign = $y->{sign};
3200       # calculate A/B
3201       $y->bone(); $y->{_m} = $a; my $y_d = $self->bone(); $y_d->{_m} = $b;
3202       $y->bdiv($y_d, @r);
3203       $y->{sign} = $y_sign;
3204       return $y;
3205       }
3206     }
3207
3208   # handle all other cases
3209   #  X  Y
3210   # +x +y 0 to PI/2
3211   # -x +y PI/2 to PI
3212   # +x -y 0 to -PI/2
3213   # -x -y -PI/2 to -PI 
3214
3215   my $y_sign = $y->{sign};
3216
3217   # divide $x by $y
3218   $y->bdiv($x, $scale) unless $x->is_one();
3219   $y->batan(@r);
3220
3221   # restore sign
3222   $y->{sign} = $y_sign;
3223
3224   $y;
3225   }
3226
3227 sub batan
3228   {
3229   # Calculate a arcus tangens of x.
3230   my ($x,@r) = @_;
3231   my $self = ref($x);
3232
3233   # taylor:       x^3   x^5   x^7   x^9
3234   #    atan = x - --- + --- - --- + --- ...
3235   #                3     5     7     9 
3236
3237   # we need to limit the accuracy to protect against overflow
3238   my $fallback = 0;
3239   my ($scale,@params);
3240   ($x,@params) = $x->_find_round_parameters(@r);
3241     
3242   #         constant object       or error in _find_round_parameters?
3243   return $x if $x->modify('batan') || $x->is_nan();
3244
3245   if ($x->{sign} =~ /^[+-]inf\z/)
3246     {
3247     # +inf result is PI/2
3248     # -inf result is -PI/2
3249     # calculate PI/2
3250     my $pi = $self->bpi(@r);
3251     # modify $x in place
3252     $x->{_m} = $pi->{_m};
3253     $x->{_e} = $pi->{_e};
3254     $x->{_es} = $pi->{_es};
3255     # -y => -PI/2, +y => PI/2
3256     $x->{sign} = substr($x->{sign},0,1);                # +inf => +
3257     $MBI->_div($x->{_m}, $MBI->_new(2));
3258     return $x;
3259     }
3260
3261   return $x->bzero(@r) if $x->is_zero();
3262
3263   # no rounding at all, so must use fallback
3264   if (scalar @params == 0)
3265     {
3266     # simulate old behaviour
3267     $params[0] = $self->div_scale();    # and round to it as accuracy
3268     $params[1] = undef;                 # disable P
3269     $scale = $params[0]+4;              # at least four more for proper round
3270     $params[2] = $r[2];                 # round mode by caller or undef
3271     $fallback = 1;                      # to clear a/p afterwards
3272     }
3273   else
3274     {
3275     # the 4 below is empirical, and there might be cases where it is not
3276     # enough...
3277     $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3278     }
3279
3280   # 1 or -1 => PI/4
3281   # inlined is_one() && is_one('-')
3282   if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
3283     {
3284     my $pi = $self->bpi($scale - 3);
3285     # modify $x in place
3286     $x->{_m} = $pi->{_m};
3287     $x->{_e} = $pi->{_e};
3288     $x->{_es} = $pi->{_es};
3289     # leave the sign of $x alone (+1 => +PI/4, -1 => -PI/4)
3290     $MBI->_div($x->{_m}, $MBI->_new(4));
3291     return $x;
3292     }
3293   
3294   # This series is only valid if -1 < x < 1, so for other x we need to
3295   # to calculate PI/2 - atan(1/x):
3296   my $one = $MBI->_new(1);
3297   my $pi = undef;
3298   if ($x->{_es} eq '+' && ($MBI->_acmp($x->{_m},$one) >= 0))
3299     {
3300     # calculate PI/2
3301     $pi = $self->bpi($scale - 3);
3302     $MBI->_div($pi->{_m}, $MBI->_new(2));
3303     # calculate 1/$x:
3304     my $x_copy = $x->copy();
3305     # modify $x in place
3306     $x->bone(); $x->bdiv($x_copy,$scale);
3307     }
3308
3309   # when user set globals, they would interfere with our calculation, so
3310   # disable them and later re-enable them
3311   no strict 'refs';
3312   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
3313   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
3314   # we also need to disable any set A or P on $x (_find_round_parameters took
3315   # them already into account), since these would interfere, too
3316   delete $x->{_a}; delete $x->{_p};
3317   # need to disable $upgrade in BigInt, to avoid deep recursion
3318   local $Math::BigInt::upgrade = undef;
3319  
3320   my $last = 0;
3321   my $over = $x * $x;                   # X ^ 2
3322   my $x2 = $over->copy();               # X ^ 2; difference between terms
3323   $over->bmul($x);                      # X ^ 3 as starting value
3324   my $sign = 1;                         # start with -=
3325   my $below = $self->new(3);
3326   my $two = $self->new(2);
3327   delete $x->{_a}; delete $x->{_p};
3328
3329   my $limit = $self->new("1E-". ($scale-1));
3330   #my $steps = 0;
3331   while (3 < 5)
3332     {
3333     # we calculate the next term, and add it to the last
3334     # when the next term is below our limit, it won't affect the outcome
3335     # anymore, so we stop:
3336     my $next = $over->copy()->bdiv($below,$scale);
3337     last if $next->bacmp($limit) <= 0;
3338
3339     if ($sign == 0)
3340       {
3341       $x->badd($next);
3342       }
3343     else
3344       {
3345       $x->bsub($next);
3346       }
3347     $sign = 1-$sign;                                    # alternate
3348     # calculate things for the next term
3349     $over->bmul($x2);                                   # $x*$x
3350     $below->badd($two);                                 # n += 2
3351     }
3352
3353   if (defined $pi)
3354     {
3355     my $x_copy = $x->copy();
3356     # modify $x in place
3357     $x->{_m} = $pi->{_m};
3358     $x->{_e} = $pi->{_e};
3359     $x->{_es} = $pi->{_es};
3360     # PI/2 - $x
3361     $x->bsub($x_copy);
3362     }
3363
3364   # shortcut to not run through _find_round_parameters again
3365   if (defined $params[0])
3366     {
3367     $x->bround($params[0],$params[2]);          # then round accordingly
3368     }
3369   else
3370     {
3371     $x->bfround($params[1],$params[2]);         # then round accordingly
3372     }
3373   if ($fallback)
3374     {
3375     # clear a/p after round, since user did not request it
3376     delete $x->{_a}; delete $x->{_p};
3377     }
3378   # restore globals
3379   $$abr = $ab; $$pbr = $pb;
3380   $x;
3381   }
3382
3383 ###############################################################################
3384 # rounding functions
3385
3386 sub bfround
3387   {
3388   # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
3389   # $n == 0 means round to integer
3390   # expects and returns normalized numbers!
3391   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
3392
3393   my ($scale,$mode) = $x->_scale_p(@_);
3394   return $x if !defined $scale || $x->modify('bfround'); # no-op
3395
3396   # never round a 0, +-inf, NaN
3397   if ($x->is_zero())
3398     {
3399     $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
3400     return $x; 
3401     }
3402   return $x if $x->{sign} !~ /^[+-]$/;
3403
3404   # don't round if x already has lower precision
3405   return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
3406
3407   $x->{_p} = $scale;                    # remember round in any case
3408   delete $x->{_a};                      # and clear A
3409   if ($scale < 0)
3410     {
3411     # round right from the '.'
3412
3413     return $x if $x->{_es} eq '+';              # e >= 0 => nothing to round
3414
3415     $scale = -$scale;                           # positive for simplicity
3416     my $len = $MBI->_len($x->{_m});             # length of mantissa
3417
3418     # the following poses a restriction on _e, but if _e is bigger than a
3419     # scalar, you got other problems (memory etc) anyway
3420     my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e})));   # digits after dot
3421     my $zad = 0;                                # zeros after dot
3422     $zad = $dad - $len if (-$dad < -$len);      # for 0.00..00xxx style
3423    
3424     # print "scale $scale dad $dad zad $zad len $len\n";
3425     # number  bsstr   len zad dad       
3426     # 0.123   123e-3    3   0 3
3427     # 0.0123  123e-4    3   1 4
3428     # 0.001   1e-3      1   2 3
3429     # 1.23    123e-2    3   0 2
3430     # 1.2345  12345e-4  5   0 4
3431
3432     # do not round after/right of the $dad
3433     return $x if $scale > $dad;                 # 0.123, scale >= 3 => exit
3434
3435     # round to zero if rounding inside the $zad, but not for last zero like:
3436     # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
3437     return $x->bzero() if $scale < $zad;
3438     if ($scale == $zad)                 # for 0.006, scale -3 and trunc
3439       {
3440       $scale = -$len;
3441       }
3442     else
3443       {
3444       # adjust round-point to be inside mantissa
3445       if ($zad != 0)
3446         {
3447         $scale = $scale-$zad;
3448         }
3449       else
3450         {
3451         my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;    # digits before dot
3452         $scale = $dbd+$scale;
3453         }
3454       }
3455     }
3456   else
3457     {
3458     # round left from the '.'
3459
3460     # 123 => 100 means length(123) = 3 - $scale (2) => 1
3461
3462     my $dbt = $MBI->_len($x->{_m}); 
3463     # digits before dot 
3464     my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
3465     # should be the same, so treat it as this 
3466     $scale = 1 if $scale == 0; 
3467     # shortcut if already integer 
3468     return $x if $scale == 1 && $dbt <= $dbd; 
3469     # maximum digits before dot 
3470     ++$dbd;
3471
3472     if ($scale > $dbd) 
3473        { 
3474        # not enough digits before dot, so round to zero 
3475        return $x->bzero; 
3476        }
3477     elsif ( $scale == $dbd )
3478        { 
3479        # maximum 
3480        $scale = -$dbt; 
3481        } 
3482     else
3483        { 
3484        $scale = $dbd - $scale; 
3485        }
3486     }
3487   # pass sign to bround for rounding modes '+inf' and '-inf'
3488   my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
3489   $m->bround($scale,$mode);
3490   $x->{_m} = $m->{value};                       # get our mantissa back
3491   $x->bnorm();
3492   }
3493
3494 sub bround
3495   {
3496   # accuracy: preserve $N digits, and overwrite the rest with 0's
3497   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
3498
3499   if (($_[0] || 0) < 0)
3500     {
3501     require Carp; Carp::croak ('bround() needs positive accuracy');
3502     }
3503
3504   my ($scale,$mode) = $x->_scale_a(@_);
3505   return $x if !defined $scale || $x->modify('bround'); # no-op
3506
3507   # scale is now either $x->{_a}, $accuracy, or the user parameter
3508   # test whether $x already has lower accuracy, do nothing in this case 
3509   # but do round if the accuracy is the same, since a math operation might
3510   # want to round a number with A=5 to 5 digits afterwards again
3511   return $x if defined $x->{_a} && $x->{_a} < $scale;
3512
3513   # scale < 0 makes no sense
3514   # scale == 0 => keep all digits
3515   # never round a +-inf, NaN
3516   return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
3517
3518   # 1: never round a 0
3519   # 2: if we should keep more digits than the mantissa has, do nothing
3520   if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
3521     {
3522     $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
3523     return $x; 
3524     }
3525
3526   # pass sign to bround for '+inf' and '-inf' rounding modes
3527   my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
3528
3529   $m->bround($scale,$mode);             # round mantissa
3530   $x->{_m} = $m->{value};               # get our mantissa back
3531   $x->{_a} = $scale;                    # remember rounding
3532   delete $x->{_p};                      # and clear P
3533   $x->bnorm();                          # del trailing zeros gen. by bround()
3534   }
3535
3536 sub bfloor
3537   {
3538   # round towards minus infinity
3539   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3540
3541   return $x if $x->modify('bfloor');
3542    
3543   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3544
3545   # if $x has digits after dot
3546   if ($x->{_es} eq '-')
3547     {
3548     $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3549     $x->{_e} = $MBI->_zero();                   # trunc/norm    
3550     $x->{_es} = '+';                            # abs e
3551     $MBI->_inc($x->{_m}) if $x->{sign} eq '-';  # increment if negative
3552     }
3553   $x->round($a,$p,$r);
3554   }
3555
3556 sub bceil
3557   {
3558   # round towards plus infinity
3559   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3560
3561   return $x if $x->modify('bceil');
3562   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3563
3564   # if $x has digits after dot
3565   if ($x->{_es} eq '-')
3566     {
3567     $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3568     $x->{_e} = $MBI->_zero();                   # trunc/norm    
3569     $x->{_es} = '+';                            # abs e
3570     if ($x->{sign} eq '+') {
3571         $MBI->_inc($x->{_m});                   # increment if positive
3572     } else {
3573         $x->{sign} = '+' if $MBI->_is_zero($x->{_m});   # avoid -0
3574     }
3575     }
3576   $x->round($a,$p,$r);
3577   }
3578
3579 sub bint
3580   {
3581   # round towards zero
3582   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
3583
3584   return $x if $x->modify('bint');
3585   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3586
3587   # if $x has digits after the decimal point
3588   if ($x->{_es} eq '-')
3589     {
3590     $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3591     $x->{_e} = $MBI->_zero();                     # truncate/normalize
3592     $x->{_es} = '+';                              # abs e
3593     $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # avoid -0
3594     }
3595   $x->round($a,$p,$r);
3596   }
3597
3598 sub brsft
3599   {
3600   # shift right by $y (divide by power of $n)
3601   
3602   # set up parameters
3603   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3604   # objectify is costly, so avoid it
3605   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3606     {
3607     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3608     }
3609
3610   return $x if $x->modify('brsft');
3611   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3612
3613   $n = 2 if !defined $n; $n = $self->new($n);
3614
3615   # negative amount?
3616   return $x->blsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3617
3618   # the following call to bdiv() will return either quo or (quo,remainder):
3619   $x->bdiv($n->bpow($y),$a,$p,$r,$y);
3620   }
3621
3622 sub blsft
3623   {
3624   # shift left by $y (multiply by power of $n)
3625   
3626   # set up parameters
3627   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3628   # objectify is costly, so avoid it
3629   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3630     {
3631     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3632     }
3633
3634   return $x if $x->modify('blsft');
3635   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
3636
3637   $n = 2 if !defined $n; $n = $self->new($n);
3638
3639   # negative amount?
3640   return $x->brsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3641
3642   $x->bmul($n->bpow($y),$a,$p,$r,$y);
3643   }
3644
3645 ###############################################################################
3646
3647 sub DESTROY
3648   {
3649   # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
3650   }
3651
3652 sub AUTOLOAD
3653   {
3654   # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
3655   # or falling back to MBI::bxxx()
3656   my $name = $AUTOLOAD;
3657
3658   $name =~ s/(.*):://;  # split package
3659   my $c = $1 || $class;
3660   no strict 'refs';
3661   $c->import() if $IMPORT == 0;
3662   if (!_method_alias($name))
3663     {
3664     if (!defined $name)
3665       {
3666       # delayed load of Carp and avoid recursion        
3667       require Carp;
3668       Carp::croak ("$c: Can't call a method without name");
3669       }
3670     if (!_method_hand_up($name))
3671       {
3672       # delayed load of Carp and avoid recursion        
3673       require Carp;
3674       Carp::croak ("Can't call $c\-\>$name, not a valid method");
3675       }
3676     # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
3677     $name =~ s/^f/b/;
3678     return &{"Math::BigInt"."::$name"}(@_);
3679     }
3680   my $bname = $name; $bname =~ s/^f/b/;
3681   $c .= "::$name";
3682   *{$c} = \&{$bname};
3683   &{$c};        # uses @_
3684   }
3685
3686 sub exponent
3687   {
3688   # return a copy of the exponent
3689   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3690
3691   if ($x->{sign} !~ /^[+-]$/)
3692     {
3693     my $s = $x->{sign}; $s =~ s/^[+-]//;
3694     return Math::BigInt->new($s);               # -inf, +inf => +inf
3695     }
3696   Math::BigInt->new( $x->{_es} . $MBI->_str($x->{_e}));
3697   }
3698
3699 sub mantissa
3700   {
3701   # return a copy of the mantissa
3702   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3703  
3704   if ($x->{sign} !~ /^[+-]$/)
3705     {
3706     my $s = $x->{sign}; $s =~ s/^[+]//;
3707     return Math::BigInt->new($s);               # -inf, +inf => +inf
3708     }
3709   my $m = Math::BigInt->new( $MBI->_str($x->{_m}));
3710   $m->bneg() if $x->{sign} eq '-';
3711
3712   $m;
3713   }
3714
3715 sub parts
3716   {
3717   # return a copy of both the exponent and the mantissa
3718   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3719
3720   if ($x->{sign} !~ /^[+-]$/)
3721     {
3722     my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
3723     return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
3724     }
3725   my $m = Math::BigInt->bzero();
3726   $m->{value} = $MBI->_copy($x->{_m});
3727   $m->bneg() if $x->{sign} eq '-';
3728   ($m, Math::BigInt->new( $x->{_es} . $MBI->_num($x->{_e}) ));
3729   }
3730
3731 ##############################################################################
3732 # private stuff (internal use only)
3733
3734 sub import
3735   {
3736   my $self = shift;
3737   my $l = scalar @_;
3738   my $lib = ''; my @a;
3739   my $lib_kind = 'try';
3740   $IMPORT=1;
3741   for ( my $i = 0; $i < $l ; $i++)
3742     {
3743     if ( $_[$i] eq ':constant' )
3744       {
3745       # This causes overlord er load to step in. 'binary' and 'integer'
3746       # are handled by BigInt.
3747       overload::constant float => sub { $self->new(shift); }; 
3748       }
3749     elsif ($_[$i] eq 'upgrade')
3750       {
3751       # this causes upgrading
3752       $upgrade = $_[$i+1];              # or undef to disable
3753       $i++;
3754       }
3755     elsif ($_[$i] eq 'downgrade')
3756       {
3757       # this causes downgrading
3758       $downgrade = $_[$i+1];            # or undef to disable
3759       $i++;
3760       }
3761     elsif ($_[$i] =~ /^(lib|try|only)\z/)
3762       {
3763       # alternative library
3764       $lib = $_[$i+1] || '';            # default Calc
3765       $lib_kind = $1;                   # lib, try or only
3766       $i++;
3767       }
3768     elsif ($_[$i] eq 'with')
3769       {
3770       # alternative class for our private parts()
3771       # XXX: no longer supported
3772       # $MBI = $_[$i+1] || 'Math::BigInt';
3773       $i++;
3774       }
3775     else
3776       {
3777       push @a, $_[$i];
3778       }
3779     }
3780
3781   $lib =~ tr/a-zA-Z0-9,://cd;           # restrict to sane characters
3782   # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
3783   my $mbilib = eval { Math::BigInt->config()->{lib} };
3784   if ((defined $mbilib) && ($MBI eq 'Math::BigInt::Calc'))
3785     {
3786     # MBI already loaded
3787     Math::BigInt->import( $lib_kind, "$lib,$mbilib", 'objectify');
3788     }
3789   else
3790     {
3791     # MBI not loaded, or with ne "Math::BigInt::Calc"
3792     $lib .= ",$mbilib" if defined $mbilib;
3793     $lib =~ s/^,//;                             # don't leave empty 
3794     
3795     # replacement library can handle lib statement, but also could ignore it
3796     
3797     # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
3798     # used in the same script, or eval inside import(). So we require MBI:
3799     require Math::BigInt;
3800     Math::BigInt->import( $lib_kind => $lib, 'objectify' );
3801     }
3802   if ($@)
3803     {
3804     require Carp; Carp::croak ("Couldn't load $lib: $! $@");
3805     }
3806   # find out which one was actually loaded
3807   $MBI = Math::BigInt->config()->{lib};
3808
3809   # register us with MBI to get notified of future lib changes
3810   Math::BigInt::_register_callback( $self, sub { $MBI = $_[0]; } );
3811
3812   $self->export_to_level(1,$self,@a);           # export wanted functions
3813   }
3814
3815 sub bnorm
3816   {
3817   # adjust m and e so that m is smallest possible
3818   my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
3819
3820   return $x if $x->{sign} !~ /^[+-]$/;          # inf, nan etc
3821
3822   my $zeros = $MBI->_zeros($x->{_m});           # correct for trailing zeros
3823   if ($zeros != 0)
3824     {
3825     my $z = $MBI->_new($zeros);
3826     $x->{_m} = $MBI->_rsft ($x->{_m}, $z, 10);
3827     if ($x->{_es} eq '-')
3828       {
3829       if ($MBI->_acmp($x->{_e},$z) >= 0)
3830         {
3831         $x->{_e} = $MBI->_sub ($x->{_e}, $z);
3832         $x->{_es} = '+' if $MBI->_is_zero($x->{_e});
3833         }
3834       else
3835         {
3836         $x->{_e} = $MBI->_sub ( $MBI->_copy($z), $x->{_e});
3837         $x->{_es} = '+';
3838         }
3839       }
3840     else
3841       {
3842       $x->{_e} = $MBI->_add ($x->{_e}, $z);
3843       }
3844     }
3845   else
3846     {
3847     # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
3848     # zeros). So, for something like 0Ey, set y to 1, and -0 => +0
3849     $x->{sign} = '+', $x->{_es} = '+', $x->{_e} = $MBI->_one()
3850      if $MBI->_is_zero($x->{_m});
3851     }
3852
3853   $x;                                   # MBI bnorm is no-op, so do not call it
3854   } 
3855  
3856 ##############################################################################
3857
3858 sub as_hex
3859   {
3860   # return number as hexadecimal string (only for integers defined)
3861   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3862
3863   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
3864   return '0x0' if $x->is_zero();
3865
3866   return $nan if $x->{_es} ne '+';              # how to do 1e-1 in hex!?
3867
3868   my $z = $MBI->_copy($x->{_m});
3869   if (! $MBI->_is_zero($x->{_e}))               # > 0 
3870     {
3871     $MBI->_lsft( $z, $x->{_e},10);
3872     }
3873   $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3874   $z->as_hex();
3875   }
3876
3877 sub as_bin
3878   {
3879   # return number as binary digit string (only for integers defined)
3880   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3881
3882   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
3883   return '0b0' if $x->is_zero();
3884
3885   return $nan if $x->{_es} ne '+';              # how to do 1e-1 in hex!?
3886
3887   my $z = $MBI->_copy($x->{_m});
3888   if (! $MBI->_is_zero($x->{_e}))               # > 0 
3889     {
3890     $MBI->_lsft( $z, $x->{_e},10);
3891     }
3892   $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3893   $z->as_bin();
3894   }
3895
3896 sub as_oct
3897   {
3898   # return number as octal digit string (only for integers defined)
3899   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3900
3901   return $x->bstr() if $x->{sign} !~ /^[+-]$/;  # inf, nan etc
3902   return '0' if $x->is_zero();
3903
3904   return $nan if $x->{_es} ne '+';              # how to do 1e-1 in hex!?
3905
3906   my $z = $MBI->_copy($x->{_m});
3907   if (! $MBI->_is_zero($x->{_e}))               # > 0 
3908     {
3909     $MBI->_lsft( $z, $x->{_e},10);
3910     }
3911   $z = Math::BigInt->new( $x->{sign} . $MBI->_num($z));
3912   $z->as_oct();
3913   }
3914
3915 sub as_number
3916   {
3917   # return copy as a bigint representation of this BigFloat number
3918   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
3919
3920   return $x if $x->modify('as_number');
3921
3922   if (!$x->isa('Math::BigFloat'))
3923     {
3924     # if the object can as_number(), use it
3925     return $x->as_number() if $x->can('as_number');
3926     # otherwise, get us a float and then a number
3927     $x = $x->can('as_float') ? $x->as_float() : $self->new(0+"$x");
3928     }
3929
3930   return Math::BigInt->binf($x->sign()) if $x->is_inf();
3931   return Math::BigInt->bnan()           if $x->is_nan();
3932
3933   my $z = $MBI->_copy($x->{_m});
3934   if ($x->{_es} eq '-')                 # < 0
3935     {
3936     $MBI->_rsft( $z, $x->{_e},10);
3937     } 
3938   elsif (! $MBI->_is_zero($x->{_e}))    # > 0 
3939     {
3940     $MBI->_lsft( $z, $x->{_e},10);
3941     }
3942   $z = Math::BigInt->new( $x->{sign} . $MBI->_str($z));
3943   $z;
3944   }
3945
3946 sub length
3947   {
3948   my $x = shift;
3949   my $class = ref($x) || $x;
3950   $x = $class->new(shift) unless ref($x);
3951
3952   return 1 if $MBI->_is_zero($x->{_m});
3953
3954   my $len = $MBI->_len($x->{_m});
3955   $len += $MBI->_num($x->{_e}) if $x->{_es} eq '+';
3956   if (wantarray())
3957     {
3958     my $t = 0;
3959     $t = $MBI->_num($x->{_e}) if $x->{_es} eq '-';
3960     return ($len, $t);
3961     }
3962   $len;
3963   }
3964
3965 1;
3966
3967 __END__
3968
3969 =pod
3970
3971 =head1 NAME
3972
3973 Math::BigFloat - Arbitrary size floating point math package
3974
3975 =head1 SYNOPSIS
3976
3977  use Math::BigFloat;
3978
3979  # Number creation
3980  my $x = Math::BigFloat->new($str);     # defaults to 0
3981  my $y = $x->copy();                    # make a true copy
3982  my $nan  = Math::BigFloat->bnan();     # create a NotANumber
3983  my $zero = Math::BigFloat->bzero();    # create a +0
3984  my $inf = Math::BigFloat->binf();      # create a +inf
3985  my $inf = Math::BigFloat->binf('-');   # create a -inf
3986  my $one = Math::BigFloat->bone();      # create a +1
3987  my $mone = Math::BigFloat->bone('-');  # create a -1
3988
3989  my $pi = Math::BigFloat->bpi(100);     # PI to 100 digits
3990
3991  # the following examples compute their result to 100 digits accuracy:
3992  my $cos  = Math::BigFloat->new(1)->bcos(100);        # cosinus(1)
3993  my $sin  = Math::BigFloat->new(1)->bsin(100);        # sinus(1)
3994  my $atan = Math::BigFloat->new(1)->batan(100);       # arcus tangens(1)
3995
3996  my $atan2 = Math::BigFloat->new(  1 )->batan2( 1 ,100); # batan(1)
3997  my $atan2 = Math::BigFloat->new(  1 )->batan2( 8 ,100); # batan(1/8)
3998  my $atan2 = Math::BigFloat->new( -2 )->batan2( 1 ,100); # batan(-2)
3999
4000  # Testing
4001  $x->is_zero();          # true if arg is +0
4002  $x->is_nan();           # true if arg is NaN
4003  $x->is_one();           # true if arg is +1
4004  $x->is_one('-');        # true if arg is -1
4005  $x->is_odd();           # true if odd, false for even
4006  $x->is_even();          # true if even, false for odd
4007  $x->is_pos();           # true if >= 0
4008  $x->is_neg();           # true if <  0
4009  $x->is_inf(sign);       # true if +inf, or -inf (default is '+')
4010
4011  $x->bcmp($y);           # compare numbers (undef,<0,=0,>0)
4012  $x->bacmp($y);          # compare absolutely (undef,<0,=0,>0)
4013  $x->sign();             # return the sign, either +,- or NaN
4014  $x->digit($n);          # return the nth digit, counting from right
4015  $x->digit(-$n);         # return the nth digit, counting from left 
4016
4017  # The following all modify their first argument. If you want to pre-
4018  # serve $x, use $z = $x->copy()->bXXX($y); See under L</CAVEATS> for
4019  # necessary when mixing $a = $b assignments with non-overloaded math.
4020
4021  # set 
4022  $x->bzero();            # set $i to 0
4023  $x->bnan();             # set $i to NaN
4024  $x->bone();             # set $x to +1
4025  $x->bone('-');          # set $x to -1
4026  $x->binf();             # set $x to inf
4027  $x->binf('-');          # set $x to -inf
4028
4029  $x->bneg();             # negation
4030  $x->babs();             # absolute value
4031  $x->bnorm();            # normalize (no-op)
4032  $x->bnot();             # two's complement (bit wise not)
4033  $x->binc();             # increment x by 1
4034  $x->bdec();             # decrement x by 1
4035
4036  $x->badd($y);           # addition (add $y to $x)
4037  $x->bsub($y);           # subtraction (subtract $y from $x)
4038  $x->bmul($y);           # multiplication (multiply $x by $y)
4039  $x->bdiv($y);           # divide, set $x to quotient
4040                          # return (quo,rem) or quo if scalar
4041
4042  $x->bmod($y);           # modulus ($x % $y)
4043  $x->bpow($y);           # power of arguments ($x ** $y)
4044  $x->bmodpow($exp,$mod); # modular exponentiation (($num**$exp) % $mod))
4045  $x->blsft($y, $n);      # left shift by $y places in base $n
4046  $x->brsft($y, $n);      # right shift by $y places in base $n
4047                          # returns (quo,rem) or quo if in scalar context
4048
4049  $x->blog();             # logarithm of $x to base e (Euler's number)
4050  $x->blog($base);        # logarithm of $x to base $base (f.i. 2)
4051  $x->bexp();             # calculate e ** $x where e is Euler's number
4052
4053  $x->band($y);           # bit-wise and
4054  $x->bior($y);           # bit-wise inclusive or
4055  $x->bxor($y);           # bit-wise exclusive or
4056  $x->bnot();             # bit-wise not (two's complement)
4057
4058  $x->bsqrt();            # calculate square-root
4059  $x->broot($y);          # $y'th root of $x (e.g. $y == 3 => cubic root)
4060  $x->bfac();             # factorial of $x (1*2*3*4*..$x)
4061
4062  $x->bround($N);         # accuracy: preserve $N digits
4063  $x->bfround($N);        # precision: round to the $Nth digit
4064
4065  $x->bfloor();           # return integer less or equal than $x
4066  $x->bceil();            # return integer greater or equal than $x
4067  $x->bint();             # round towards zero
4068
4069   # The following do not modify their arguments:
4070
4071  bgcd(@values);          # greatest common divisor
4072  blcm(@values);          # lowest common multiplicator
4073
4074  $x->bstr();             # return string
4075  $x->bsstr();            # return string in scientific notation
4076
4077  $x->as_int();           # return $x as BigInt 
4078  $x->exponent();         # return exponent as BigInt
4079  $x->mantissa();         # return mantissa as BigInt
4080  $x->parts();            # return (mantissa,exponent) as BigInt
4081
4082  $x->length();           # number of digits (w/o sign and '.')
4083  ($l,$f) = $x->length(); # number of digits, and length of fraction
4084
4085  $x->precision();        # return P of $x (or global, if P of $x undef)
4086  $x->precision($n);      # set P of $x to $n
4087  $x->accuracy();         # return A of $x (or global, if A of $x undef)
4088  $x->accuracy($n);       # set A $x to $n
4089
4090  # these get/set the appropriate global value for all BigFloat objects
4091  Math::BigFloat->precision();   # Precision
4092  Math::BigFloat->accuracy();    # Accuracy
4093  Math::BigFloat->round_mode();  # rounding mode
4094
4095 =head1 DESCRIPTION
4096
4097 All operators (including basic math operations) are overloaded if you
4098 declare your big floating point numbers as
4099
4100   $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
4101
4102 Operations with overloaded operators preserve the arguments, which is
4103 exactly what you expect.
4104
4105 =head2 Input
4106
4107 Input to these routines are either BigFloat objects, or strings of the
4108 following four forms:
4109
4110 =over
4111
4112 =item *
4113
4114 C</^[+-]\d+$/>
4115
4116 =item *
4117
4118 C</^[+-]\d+\.\d*$/>
4119
4120 =item *
4121
4122 C</^[+-]\d+E[+-]?\d+$/>
4123
4124 =item *
4125
4126 C</^[+-]\d*\.\d+E[+-]?\d+$/>
4127
4128 =back
4129
4130 all with optional leading and trailing zeros and/or spaces. Additionally,
4131 numbers are allowed to have an underscore between any two digits.
4132
4133 Empty strings as well as other illegal numbers results in 'NaN'.
4134
4135 bnorm() on a BigFloat object is now effectively a no-op, since the numbers 
4136 are always stored in normalized form. On a string, it creates a BigFloat 
4137 object.
4138
4139 =head2 Output
4140
4141 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
4142
4143 The string output will always have leading and trailing zeros stripped and drop
4144 a plus sign. C<bstr()> will give you always the form with a decimal point,
4145 while C<bsstr()> (s for scientific) gives you the scientific notation.
4146
4147         Input                   bstr()          bsstr()
4148         '-0'                    '0'             '0E1'
4149         '  -123 123 123'        '-123123123'    '-123123123E0'
4150         '00.0123'               '0.0123'        '123E-4'
4151         '123.45E-2'             '1.2345'        '12345E-4'
4152         '10E+3'                 '10000'         '1E4'
4153
4154 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
4155 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
4156 return either undef, <0, 0 or >0 and are suited for sort.
4157
4158 Actual math is done by using the class defined with C<< with => Class; >>
4159 (which defaults to BigInts) to represent the mantissa and exponent.
4160
4161 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to 
4162 represent the result when input arguments are not numbers, and 'inf' and
4163 '-inf' are used to represent positive and negative infinity, respectively.
4164
4165 =head2 mantissa(), exponent() and parts()
4166
4167 mantissa() and exponent() return the said parts of the BigFloat
4168 as BigInts such that:
4169
4170         $m = $x->mantissa();
4171         $e = $x->exponent();
4172         $y = $m * ( 10 ** $e );
4173         print "ok\n" if $x == $y;
4174
4175 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
4176
4177 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
4178
4179 Currently the mantissa is reduced as much as possible, favouring higher
4180 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
4181 This might change in the future, so do not depend on it.
4182
4183 =head2 Accuracy vs. Precision
4184
4185 See also: L<Rounding|/Rounding>.
4186
4187 Math::BigFloat supports both precision (rounding to a certain place before or
4188 after the dot) and accuracy (rounding to a certain number of digits). For a
4189 full documentation, examples and tips on these topics please see the large
4190 section about rounding in L<Math::BigInt>.
4191
4192 Since things like C<sqrt(2)> or C<1 / 3> must presented with a limited
4193 accuracy lest a operation consumes all resources, each operation produces
4194 no more than the requested number of digits.
4195
4196 If there is no global precision or accuracy set, B<and> the operation in
4197 question was not called with a requested precision or accuracy, B<and> the
4198 input $x has no accuracy or precision set, then a fallback parameter will
4199 be used. For historical reasons, it is called C<div_scale> and can be accessed
4200 via:
4201
4202         $d = Math::BigFloat->div_scale();       # query
4203         Math::BigFloat->div_scale($n);          # set to $n digits
4204
4205 The default value for C<div_scale> is 40.
4206
4207 In case the result of one operation has more digits than specified,
4208 it is rounded. The rounding mode taken is either the default mode, or the one
4209 supplied to the operation after the I<scale>:
4210
4211     $x = Math::BigFloat->new(2);
4212     Math::BigFloat->accuracy(5);              # 5 digits max
4213     $y = $x->copy()->bdiv(3);                 # will give 0.66667
4214     $y = $x->copy()->bdiv(3,6);               # will give 0.666667
4215     $y = $x->copy()->bdiv(3,6,undef,'odd');   # will give 0.666667
4216     Math::BigFloat->round_mode('zero');
4217     $y = $x->copy()->bdiv(3,6);               # will also give 0.666667
4218
4219 Note that C<< Math::BigFloat->accuracy() >> and C<< Math::BigFloat->precision() >>
4220 set the global variables, and thus B<any> newly created number will be subject
4221 to the global rounding B<immediately>. This means that in the examples above, the
4222 C<3> as argument to C<bdiv()> will also get an accuracy of B<5>.
4223
4224 It is less confusing to either calculate the result fully, and afterwards
4225 round it explicitly, or use the additional parameters to the math
4226 functions like so:
4227
4228         use Math::BigFloat;
4229         $x = Math::BigFloat->new(2);
4230         $y = $x->copy()->bdiv(3);
4231         print $y->bround(5),"\n";               # will give 0.66667
4232
4233         or
4234
4235         use Math::BigFloat;
4236         $x = Math::BigFloat->new(2);
4237         $y = $x->copy()->bdiv(3,5);             # will give 0.66667
4238         print "$y\n";
4239
4240 =head2 Rounding
4241
4242 =over
4243
4244 =item ffround ( +$scale )
4245
4246 Rounds to the $scale'th place left from the '.', counting from the dot.
4247 The first digit is numbered 1. 
4248
4249 =item ffround ( -$scale )
4250
4251 Rounds to the $scale'th place right from the '.', counting from the dot.
4252
4253 =item ffround ( 0 )
4254
4255 Rounds to an integer.
4256
4257 =item fround  ( +$scale )
4258
4259 Preserves accuracy to $scale digits from the left (aka significant digits)
4260 and pads the rest with zeros. If the number is between 1 and -1, the
4261 significant digits count from the first non-zero after the '.'
4262
4263 =item fround  ( -$scale ) and fround ( 0 )
4264
4265 These are effectively no-ops.
4266
4267 =back
4268
4269 All rounding functions take as a second parameter a rounding mode from one of
4270 the following: 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'.
4271
4272 The default rounding mode is 'even'. By using
4273 C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
4274 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
4275 no longer supported.
4276 The second parameter to the round functions then overrides the default
4277 temporarily. 
4278
4279 The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
4280 'trunc' as rounding mode to make it equivalent to:
4281
4282         $x = 2.5;
4283         $y = int($x) + 2;
4284
4285 You can override this by passing the desired rounding mode as parameter to
4286 C<as_number()>:
4287
4288         $x = Math::BigFloat->new(2.5);
4289  &nbs