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