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