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