This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
a12cf28a84d1c6f89d21292751524d5dd0ace5c3
[perl5.git] / lib / Math / BigFloat.pm
1 package Math::BigFloat;
2
3
4 # Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
5 #
6
7 # The following hash values are internally used:
8 #   _e: exponent (BigInt)
9 #   _m: mantissa (absolute BigInt)
10 # sign: +,-,"NaN" if not a number
11 #   _a: accuracy
12 #   _p: precision
13 #   _f: flags, used to signal MBI not to touch our private parts
14
15 $VERSION = '1.35';
16 require 5.005;
17 use Exporter;
18 use File::Spec;
19 # use Math::BigInt;
20 @ISA =       qw( Exporter Math::BigInt);
21
22 use strict;
23 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode/;
24 use vars qw/$upgrade $downgrade/;
25 my $class = "Math::BigFloat";
26
27 use overload
28 '<=>'   =>      sub { $_[2] ?
29                       ref($_[0])->bcmp($_[1],$_[0]) : 
30                       ref($_[0])->bcmp($_[0],$_[1])},
31 'int'   =>      sub { $_[0]->as_number() },             # 'trunc' to bigint
32 ;
33
34 ##############################################################################
35 # global constants, flags and accessory
36
37 use constant MB_NEVER_ROUND => 0x0001;
38
39 # are NaNs ok?
40 my $NaNOK=1;
41 # constant for easier life
42 my $nan = 'NaN'; 
43
44 # class constants, use Class->constant_name() to access
45 $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
46 $accuracy   = undef;
47 $precision  = undef;
48 $div_scale  = 40;
49
50 $upgrade = undef;
51 $downgrade = undef;
52 my $MBI = 'Math::BigInt'; # the package we are using for our private parts
53                           # changable by use Math::BigFloat with => 'package'
54
55 ##############################################################################
56 # the old code had $rnd_mode, so we need to support it, too
57
58 sub TIESCALAR   { my ($class) = @_; bless \$round_mode, $class; }
59 sub FETCH       { return $round_mode; }
60 sub STORE       { $rnd_mode = $_[0]->round_mode($_[1]); }
61
62 BEGIN
63   { 
64   $rnd_mode   = 'even';
65   tie $rnd_mode, 'Math::BigFloat'; 
66   }
67  
68 ##############################################################################
69
70 # in case we call SUPER::->foo() and this wants to call modify()
71 # sub modify () { 0; }
72
73 {
74   # valid method aliases for AUTOLOAD
75   my %methods = map { $_ => 1 }  
76    qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
77         fint facmp fcmp fzero fnan finf finc fdec flog ffac
78         fceil ffloor frsft flsft fone flog
79       /;
80   # valid method's that can be hand-ed up (for AUTOLOAD)
81   my %hand_ups = map { $_ => 1 }  
82    qw / is_nan is_inf is_negative is_positive
83         accuracy precision div_scale round_mode fneg fabs babs fnot
84         objectify upgrade downgrade
85         bone binf bnan bzero
86       /;
87
88   sub method_alias { return exists $methods{$_[0]||''}; } 
89   sub method_hand_up { return exists $hand_ups{$_[0]||''}; } 
90 }
91
92 ##############################################################################
93 # constructors
94
95 sub new 
96   {
97   # create a new BigFloat object from a string or another bigfloat object. 
98   # _e: exponent
99   # _m: mantissa
100   # sign  => sign (+/-), or "NaN"
101
102   my ($class,$wanted,@r) = @_;
103
104   # avoid numify-calls by not using || on $wanted!
105   return $class->bzero() if !defined $wanted;   # default to 0
106   return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
107
108   my $self = {}; bless $self, $class;
109   # shortcut for bigints and its subclasses
110   if ((ref($wanted)) && (ref($wanted) ne $class))
111     {
112     $self->{_m} = $wanted->as_number();         # get us a bigint copy
113     $self->{_e} = $MBI->bzero();
114     $self->{_m}->babs();
115     $self->{sign} = $wanted->sign();
116     return $self->bnorm();
117     }
118   # got string
119   # handle '+inf', '-inf' first
120   if ($wanted =~ /^[+-]?inf$/)
121     {
122     return $downgrade->new($wanted) if $downgrade;
123
124     $self->{_e} = $MBI->bzero();
125     $self->{_m} = $MBI->bzero();
126     $self->{sign} = $wanted;
127     $self->{sign} = '+inf' if $self->{sign} eq 'inf';
128     return $self->bnorm();
129     }
130   #print "new string '$wanted'\n";
131   my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
132   if (!ref $mis)
133     {
134     die "$wanted is not a number initialized to $class" if !$NaNOK;
135     
136     return $downgrade->bnan() if $downgrade;
137     
138     $self->{_e} = $MBI->bzero();
139     $self->{_m} = $MBI->bzero();
140     $self->{sign} = $nan;
141     }
142   else
143     {
144     # make integer from mantissa by adjusting exp, then convert to bigint
145     # undef,undef to signal MBI that we don't need no bloody rounding
146     $self->{_e} = $MBI->new("$$es$$ev",undef,undef);    # exponent
147     $self->{_m} = $MBI->new("$$miv$$mfv",undef,undef);  # create mant.
148     # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
149     $self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0;             
150     $self->{sign} = $$mis;
151     }
152   # if downgrade, inf, NaN or integers go down
153
154   if ($downgrade && $self->{_e}->{sign} eq '+')
155     {
156 #   print "downgrading $$miv$$mfv"."E$$es$$ev";
157     if ($self->{_e}->is_zero())
158       {
159       $self->{_m}->{sign} = $$mis;              # negative if wanted
160       return $downgrade->new($self->{_m});
161       }
162     return $downgrade->new("$$mis$$miv$$mfv"."E$$es$$ev");
163     }
164   # print "mbf new $self->{sign} $self->{_m} e $self->{_e} ",ref($self),"\n";
165   $self->bnorm()->round(@r);            # first normalize, then round
166   }
167
168 sub _bnan
169   {
170   # used by parent class bone() to initialize number to 1
171   my $self = shift;
172   $self->{_m} = $MBI->bzero();
173   $self->{_e} = $MBI->bzero();
174   }
175
176 sub _binf
177   {
178   # used by parent class bone() to initialize number to 1
179   my $self = shift;
180   $self->{_m} = $MBI->bzero();
181   $self->{_e} = $MBI->bzero();
182   }
183
184 sub _bone
185   {
186   # used by parent class bone() to initialize number to 1
187   my $self = shift;
188   $self->{_m} = $MBI->bone();
189   $self->{_e} = $MBI->bzero();
190   }
191
192 sub _bzero
193   {
194   # used by parent class bone() to initialize number to 1
195   my $self = shift;
196   $self->{_m} = $MBI->bzero();
197   $self->{_e} = $MBI->bone();
198   }
199
200 sub isa
201   {
202   my ($self,$class) = @_;
203   return if $class =~ /^Math::BigInt/;          # we aren't one of these
204   UNIVERSAL::isa($self,$class);
205   }
206
207 sub config
208   {
209   # return (later set?) configuration data as hash ref
210   my $class = shift || 'Math::BigFloat';
211
212   my $cfg = $MBI->config();
213
214   no strict 'refs';
215   $cfg->{class} = $class;
216   $cfg->{with} = $MBI;
217   foreach (
218    qw/upgrade downgrade precision accuracy round_mode VERSION div_scale/)
219     {
220     $cfg->{lc($_)} = ${"${class}::$_"};
221     };
222   $cfg;
223   }
224
225 ##############################################################################
226 # string conversation
227
228 sub bstr 
229   {
230   # (ref to BFLOAT or num_str ) return num_str
231   # Convert number from internal format to (non-scientific) string format.
232   # internal format is always normalized (no leading zeros, "-0" => "+0")
233   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
234   #my $x = shift; my $class = ref($x) || $x;
235   #$x = $class->new(shift) unless ref($x);
236
237   #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
238   #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
239   if ($x->{sign} !~ /^[+-]$/)
240     {
241     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
242     return 'inf';                                       # +inf
243     }
244  
245   my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
246
247   my $not_zero = ! $x->is_zero();
248   if ($not_zero)
249     {
250     $es = $x->{_m}->bstr();
251     $len = CORE::length($es);
252     if (!$x->{_e}->is_zero())
253       {
254       if ($x->{_e}->sign() eq '-')
255         {
256         $dot = '';
257         if ($x->{_e} <= -$len)
258           {
259           # print "style: 0.xxxx\n";
260           my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
261           $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
262           }
263         else
264           {
265           # print "insert '.' at $x->{_e} in '$es'\n";
266           substr($es,$x->{_e},0) = '.'; $cad = $x->{_e};
267           }
268         }
269       else
270         {
271         # expand with zeros
272         $es .= '0' x $x->{_e}; $len += $x->{_e}; $cad = 0;
273         }
274       }
275     } # if not zero
276   $es = $x->{sign}.$es if $x->{sign} eq '-';
277   # if set accuracy or precision, pad with zeros
278   if ((defined $x->{_a}) && ($not_zero))
279     {
280     # 123400 => 6, 0.1234 => 4, 0.001234 => 4
281     my $zeros = $x->{_a} - $cad;                # cad == 0 => 12340
282     $zeros = $x->{_a} - $len if $cad != $len;
283     $es .= $dot.'0' x $zeros if $zeros > 0;
284     }
285   elsif ($x->{_p} || 0 < 0)
286     {
287     # 123400 => 6, 0.1234 => 4, 0.001234 => 6
288     my $zeros = -$x->{_p} + $cad;
289     $es .= $dot.'0' x $zeros if $zeros > 0;
290     }
291   $es;
292   }
293
294 sub bsstr
295   {
296   # (ref to BFLOAT or num_str ) return num_str
297   # Convert number from internal format to scientific string format.
298   # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
299   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
300   #my $x = shift; my $class = ref($x) || $x;
301   #$x = $class->new(shift) unless ref($x);
302
303   #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
304   #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
305   if ($x->{sign} !~ /^[+-]$/)
306     {
307     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
308     return 'inf';                                       # +inf
309     }
310   my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
311   my $sep = 'e'.$sign;
312   $x->{_m}->bstr().$sep.$x->{_e}->bstr();
313   }
314     
315 sub numify 
316   {
317   # Make a number from a BigFloat object
318   # simple return string and let Perl's atoi()/atof() handle the rest
319   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
320   $x->bsstr(); 
321   }
322
323 ##############################################################################
324 # public stuff (usually prefixed with "b")
325
326 # tels 2001-08-04 
327 # todo: this must be overwritten and return NaN for non-integer values
328 # band(), bior(), bxor(), too
329 #sub bnot
330 #  {
331 #  $class->SUPER::bnot($class,@_);
332 #  }
333
334 sub bcmp 
335   {
336   # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
337   # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
338
339   # set up parameters
340   my ($self,$x,$y) = (ref($_[0]),@_);
341   # objectify is costly, so avoid it
342   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
343     {
344     ($self,$x,$y) = objectify(2,@_);
345     }
346
347   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
348     {
349     # handle +-inf and NaN
350     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
351     return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
352     return +1 if $x->{sign} eq '+inf';
353     return -1 if $x->{sign} eq '-inf';
354     return -1 if $y->{sign} eq '+inf';
355     return +1;
356     }
357
358   # check sign for speed first
359   return 1 if $x->{sign} eq '+' && $y->{sign} eq '-';   # does also 0 <=> -y
360   return -1 if $x->{sign} eq '-' && $y->{sign} eq '+';  # does also -x <=> 0
361
362   # shortcut 
363   my $xz = $x->is_zero();
364   my $yz = $y->is_zero();
365   return 0 if $xz && $yz;                               # 0 <=> 0
366   return -1 if $xz && $y->{sign} eq '+';                # 0 <=> +y
367   return 1 if $yz && $x->{sign} eq '+';                 # +x <=> 0
368
369   # adjust so that exponents are equal
370   my $lxm = $x->{_m}->length();
371   my $lym = $y->{_m}->length();
372   # the numify somewhat limits our length, but makes it much faster
373   my $lx = $lxm + $x->{_e}->numify();
374   my $ly = $lym + $y->{_e}->numify();
375   my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
376   return $l <=> 0 if $l != 0;
377   
378   # lengths (corrected by exponent) are equal
379   # so make mantissa equal length by padding with zero (shift left)
380   my $diff = $lxm - $lym;
381   my $xm = $x->{_m};            # not yet copy it
382   my $ym = $y->{_m};
383   if ($diff > 0)
384     {
385     $ym = $y->{_m}->copy()->blsft($diff,10);
386     }
387   elsif ($diff < 0)
388     {
389     $xm = $x->{_m}->copy()->blsft(-$diff,10);
390     }
391   my $rc = $xm->bacmp($ym);
392   $rc = -$rc if $x->{sign} eq '-';              # -124 < -123
393   $rc <=> 0;
394   }
395
396 sub bacmp 
397   {
398   # Compares 2 values, ignoring their signs. 
399   # Returns one of undef, <0, =0, >0. (suitable for sort)
400   # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
401   
402   # set up parameters
403   my ($self,$x,$y) = (ref($_[0]),@_);
404   # objectify is costly, so avoid it
405   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
406     {
407     ($self,$x,$y) = objectify(2,@_);
408     }
409
410   # handle +-inf and NaN's
411   if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
412     {
413     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
414     return 0 if ($x->is_inf() && $y->is_inf());
415     return 1 if ($x->is_inf() && !$y->is_inf());
416     return -1;
417     }
418
419   # shortcut 
420   my $xz = $x->is_zero();
421   my $yz = $y->is_zero();
422   return 0 if $xz && $yz;                               # 0 <=> 0
423   return -1 if $xz && !$yz;                             # 0 <=> +y
424   return 1 if $yz && !$xz;                              # +x <=> 0
425
426   # adjust so that exponents are equal
427   my $lxm = $x->{_m}->length();
428   my $lym = $y->{_m}->length();
429   # the numify somewhat limits our length, but makes it much faster
430   my $lx = $lxm + $x->{_e}->numify();
431   my $ly = $lym + $y->{_e}->numify();
432   my $l = $lx - $ly;
433   return $l <=> 0 if $l != 0;
434   
435   # lengths (corrected by exponent) are equal
436   # so make mantissa equal-length by padding with zero (shift left)
437   my $diff = $lxm - $lym;
438   my $xm = $x->{_m};            # not yet copy it
439   my $ym = $y->{_m};
440   if ($diff > 0)
441     {
442     $ym = $y->{_m}->copy()->blsft($diff,10);
443     }
444   elsif ($diff < 0)
445     {
446     $xm = $x->{_m}->copy()->blsft(-$diff,10);
447     }
448   $xm->bacmp($ym) <=> 0;
449   }
450
451 sub badd 
452   {
453   # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
454   # return result as BFLOAT
455
456   # set up parameters
457   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
458   # objectify is costly, so avoid it
459   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
460     {
461     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
462     }
463
464   # inf and NaN handling
465   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
466     {
467     # NaN first
468     return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
469     # inf handling
470     if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
471       {
472       # +inf++inf or -inf+-inf => same, rest is NaN
473       return $x if $x->{sign} eq $y->{sign};
474       return $x->bnan();
475       }
476     # +-inf + something => +inf; something +-inf => +-inf
477     $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
478     return $x;
479     }
480
481   return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
482    ((!$x->isa($self)) || (!$y->isa($self)));
483
484   # speed: no add for 0+y or x+0
485   return $x->bround($a,$p,$r) if $y->is_zero();         # x+0
486   if ($x->is_zero())                                    # 0+y
487     {
488     # make copy, clobbering up x (modify in place!)
489     $x->{_e} = $y->{_e}->copy();
490     $x->{_m} = $y->{_m}->copy();
491     $x->{sign} = $y->{sign} || $nan;
492     return $x->round($a,$p,$r,$y);
493     }
494  
495   # take lower of the two e's and adapt m1 to it to match m2
496   my $e = $y->{_e};
497   $e = $MBI->bzero() if !defined $e;    # if no BFLOAT ?
498   $e = $e->copy();                      # make copy (didn't do it yet)
499   $e->bsub($x->{_e});
500   my $add = $y->{_m}->copy();
501   if ($e->{sign} eq '-')                # < 0
502     {
503     my $e1 = $e->copy()->babs();
504     #$x->{_m} *= (10 ** $e1);
505     $x->{_m}->blsft($e1,10);
506     $x->{_e} += $e;                     # need the sign of e
507     }
508   elsif (!$e->is_zero())                # > 0
509     {
510     #$add *= (10 ** $e);
511     $add->blsft($e,10);
512     }
513   # else: both e are the same, so just leave them
514   $x->{_m}->{sign} = $x->{sign};                # fiddle with signs
515   $add->{sign} = $y->{sign};
516   $x->{_m} += $add;                             # finally do add/sub
517   $x->{sign} = $x->{_m}->{sign};                # re-adjust signs
518   $x->{_m}->{sign} = '+';                       # mantissa always positiv
519   # delete trailing zeros, then round
520   return $x->bnorm()->round($a,$p,$r,$y);
521   }
522
523 sub bsub 
524   {
525   # (BigFloat or num_str, BigFloat or num_str) return BigFloat
526   # subtract second arg from first, modify first
527
528   # set up parameters
529   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
530   # objectify is costly, so avoid it
531   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
532     {
533     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
534     }
535
536   if ($y->is_zero())            # still round for not adding zero
537     {
538     return $x->round($a,$p,$r);
539     }
540   
541   $y->{sign} =~ tr/+\-/-+/;     # does nothing for NaN
542   $x->badd($y,$a,$p,$r);        # badd does not leave internal zeros
543   $y->{sign} =~ tr/+\-/-+/;     # refix $y (does nothing for NaN)
544   $x;                           # already rounded by badd()
545   }
546
547 sub binc
548   {
549   # increment arg by one
550   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
551
552   if ($x->{_e}->sign() eq '-')
553     {
554     return $x->badd($self->bone(),$a,$p,$r);    #  digits after dot
555     }
556
557   if (!$x->{_e}->is_zero())
558     {
559     $x->{_m}->blsft($x->{_e},10);               # 1e2 => 100
560     $x->{_e}->bzero();
561     }
562   # now $x->{_e} == 0
563   if ($x->{sign} eq '+')
564     {
565     $x->{_m}->binc();
566     return $x->bnorm()->bround($a,$p,$r);
567     }
568   elsif ($x->{sign} eq '-')
569     {
570     $x->{_m}->bdec();
571     $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
572     return $x->bnorm()->bround($a,$p,$r);
573     }
574   # inf, nan handling etc
575   $x->badd($self->__one(),$a,$p,$r);            # does round 
576   }
577
578 sub bdec
579   {
580   # decrement arg by one
581   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
582
583   if ($x->{_e}->sign() eq '-')
584     {
585     return $x->badd($self->bone('-'),$a,$p,$r); #  digits after dot
586     }
587
588   if (!$x->{_e}->is_zero())
589     {
590     $x->{_m}->blsft($x->{_e},10);               # 1e2 => 100
591     $x->{_e}->bzero();
592     }
593   # now $x->{_e} == 0
594   my $zero = $x->is_zero();
595   # <= 0
596   if (($x->{sign} eq '-') || $zero)
597     {
598     $x->{_m}->binc();
599     $x->{sign} = '-' if $zero;                  # 0 => 1 => -1
600     $x->{sign} = '+' if $x->{_m}->is_zero();    # -1 +1 => -0 => +0
601     return $x->bnorm()->round($a,$p,$r);
602     }
603   # > 0
604   elsif ($x->{sign} eq '+')
605     {
606     $x->{_m}->bdec();
607     return $x->bnorm()->round($a,$p,$r);
608     }
609   # inf, nan handling etc
610   $x->badd($self->bone('-'),$a,$p,$r);          # does round 
611   } 
612
613 sub blog
614   {
615   my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(2,@_);
616
617   # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
618
619   # u = x-1, v = x+1
620   #              _                               _
621   # Taylor:     |    u    1   u^3   1   u^5       |
622   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 0
623   #             |_   v    3   v^3   5   v^5      _|
624
625   # This takes much more steps to calculate the result: 
626   # u = x-1
627   #              _                               _
628   # Taylor:     |    u    1   u^2   1   u^3       |
629   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 1/2
630   #             |_   x    2   x^2   3   x^3      _|
631
632   # we need to limit the accuracy to protect against overflow
633   my $fallback = 0;
634   my $scale = 0;
635   my @params = $x->_find_round_parameters($a,$p,$r);
636
637   # no rounding at all, so must use fallback
638   if (scalar @params == 1)
639     {
640     # simulate old behaviour
641     $params[1] = $self->div_scale();    # and round to it as accuracy
642     $params[0] = undef;
643     $scale = $params[1]+4;              # at least four more for proper round
644     $params[3] = $r;                    # round mode by caller or undef
645     $fallback = 1;                      # to clear a/p afterwards
646     }
647   else
648     {
649     # the 4 below is empirical, and there might be cases where it is not
650     # enough...
651     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
652     }
653
654   return $x->bzero(@params) if $x->is_one();
655   return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
656   return $x->bone('+',@params) if $x->bcmp($base) == 0;
657
658   # when user set globals, they would interfere with our calculation, so
659   # disable then and later re-enable them
660   no strict 'refs';
661   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
662   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
663   # we also need to disable any set A or P on $x (_find_round_parameters took
664   # them already into account), since these would interfere, too
665   delete $x->{_a}; delete $x->{_p};
666   # need to disable $upgrade in BigInt, to avoid deep recursion
667   local $Math::BigInt::upgrade = undef;
668  
669   my ($case,$limit,$v,$u,$below,$factor,$two,$next,$over,$f);
670
671   if (3 < 5)
672   #if ($x <= Math::BigFloat->new("0.5"))
673     {
674     $case = 0;
675   #  print "case $case $x < 0.5\n";
676     $v = $x->copy(); $v->binc();                # v = x+1
677     $x->bdec(); $u = $x->copy();                # u = x-1; x = x-1
678     $x->bdiv($v,$scale);                        # first term: u/v
679     $below = $v->copy();
680     $over = $u->copy();
681     $u *= $u; $v *= $v;                         # u^2, v^2
682     $below->bmul($v);                           # u^3, v^3
683     $over->bmul($u);
684     $factor = $self->new(3); $f = $self->new(2);
685     }
686   #else
687   #  {
688   #  $case = 1;
689   #  print "case 1 $x > 0.5\n";
690   #  $v = $x->copy();                           # v = x
691   #  $u = $x->copy(); $u->bdec();               # u = x-1;
692   #  $x->bdec(); $x->bdiv($v,$scale);           # first term: x-1/x
693   #  $below = $v->copy();
694   #  $over = $u->copy();
695   #  $below->bmul($v);                          # u^2, v^2
696   #  $over->bmul($u);
697   #  $factor = $self->new(2); $f = $self->bone();
698   #  }
699   $limit = $self->new("1E-". ($scale-1));
700   #my $steps = 0;
701   while (3 < 5)
702     {
703     # we calculate the next term, and add it to the last
704     # when the next term is below our limit, it won't affect the outcome
705     # anymore, so we stop
706     $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
707     last if $next->bcmp($limit) <= 0;
708     $x->badd($next);
709     # print "step  $x\n";
710     # calculate things for the next term
711     $over *= $u; $below *= $v; $factor->badd($f);
712     #$steps++;
713     }
714   $x->bmul(2) if $case == 0;
715   #print "took $steps steps\n";
716   
717   # shortcut to not run trough _find_round_parameters again
718   if (defined $params[1])
719     {
720     $x->bround($params[1],$params[3]);          # then round accordingly
721     }
722   else
723     {
724     $x->bfround($params[2],$params[3]);         # then round accordingly
725     }
726   if ($fallback)
727     {
728     # clear a/p after round, since user did not request it
729     $x->{_a} = undef; $x->{_p} = undef;
730     }
731   # restore globals
732   $$abr = $ab; $$pbr = $pb;
733
734   $x;
735   }
736
737 sub blcm 
738   { 
739   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
740   # does not modify arguments, but returns new object
741   # Lowest Common Multiplicator
742
743   my ($self,@arg) = objectify(0,@_);
744   my $x = $self->new(shift @arg);
745   while (@arg) { $x = _lcm($x,shift @arg); } 
746   $x;
747   }
748
749 sub bgcd 
750   { 
751   # (BFLOAT or num_str, BFLOAT or num_str) return BINT
752   # does not modify arguments, but returns new object
753   # GCD -- Euclids algorithm Knuth Vol 2 pg 296
754    
755   my ($self,@arg) = objectify(0,@_);
756   my $x = $self->new(shift @arg);
757   while (@arg) { $x = _gcd($x,shift @arg); } 
758   $x;
759   }
760
761 ###############################################################################
762 # is_foo methods (is_negative, is_positive are inherited from BigInt)
763
764 sub is_int
765   {
766   # return true if arg (BFLOAT or num_str) is an integer
767   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
768
769   return 1 if ($x->{sign} =~ /^[+-]$/) &&       # NaN and +-inf aren't
770     $x->{_e}->{sign} eq '+';                    # 1e-1 => no integer
771   0;
772   }
773
774 sub is_zero
775   {
776   # return true if arg (BFLOAT or num_str) is zero
777   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
778
779   return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
780   0;
781   }
782
783 sub is_one
784   {
785   # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
786   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
787
788   my $sign = shift || ''; $sign = '+' if $sign ne '-';
789   return 1
790    if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one()); 
791   0;
792   }
793
794 sub is_odd
795   {
796   # return true if arg (BFLOAT or num_str) is odd or false if even
797   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
798   
799   return 1 if ($x->{sign} =~ /^[+-]$/) &&               # NaN & +-inf aren't
800     ($x->{_e}->is_zero() && $x->{_m}->is_odd()); 
801   0;
802   }
803
804 sub is_even
805   {
806   # return true if arg (BINT or num_str) is even or false if odd
807   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
808
809   return 0 if $x->{sign} !~ /^[+-]$/;                   # NaN & +-inf aren't
810   return 1 if ($x->{_e}->{sign} eq '+'                  # 123.45 is never
811      && $x->{_m}->is_even());                           # but 1200 is
812   0;
813   }
814
815 sub bmul 
816   { 
817   # multiply two numbers -- stolen from Knuth Vol 2 pg 233
818   # (BINT or num_str, BINT or num_str) return BINT
819   
820   # set up parameters
821   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
822   # objectify is costly, so avoid it
823   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
824     {
825     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
826     }
827
828   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
829
830   # inf handling
831   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
832     {
833     return $x->bnan() if $x->is_zero() || $y->is_zero(); 
834     # result will always be +-inf:
835     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
836     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
837     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
838     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
839     return $x->binf('-');
840     }
841   # handle result = 0
842   return $x->bzero() if $x->is_zero() || $y->is_zero();
843   
844   return $upgrade->bmul($x,$y,$a,$p,$r) if defined $upgrade &&
845    ((!$x->isa($self)) || (!$y->isa($self)));
846
847   # aEb * cEd = (a*c)E(b+d)
848   $x->{_m}->bmul($y->{_m});
849   $x->{_e}->badd($y->{_e});
850   # adjust sign:
851   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
852   return $x->bnorm()->round($a,$p,$r,$y);
853   }
854
855 sub bdiv 
856   {
857   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return 
858   # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
859
860   # set up parameters
861   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
862   # objectify is costly, so avoid it
863   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
864     {
865     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
866     }
867
868   return $self->_div_inf($x,$y)
869    if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
870
871   # x== 0 # also: or y == 1 or y == -1
872   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
873
874   # upgrade ?
875   return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
876
877   # we need to limit the accuracy to protect against overflow
878   my $fallback = 0;
879   my $scale = 0;
880   my @params = $x->_find_round_parameters($a,$p,$r,$y);
881
882   # no rounding at all, so must use fallback
883   if (scalar @params == 1)
884     {
885     # simulate old behaviour
886     $params[1] = $self->div_scale();    # and round to it as accuracy
887     $scale = $params[1]+4;              # at least four more for proper round
888     $params[3] = $r;                    # round mode by caller or undef
889     $fallback = 1;                      # to clear a/p afterwards
890     }
891   else
892     {
893     # the 4 below is empirical, and there might be cases where it is not
894     # enough...
895     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
896     }
897   my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
898   $scale = $lx if $lx > $scale;
899   $scale = $ly if $ly > $scale;
900   my $diff = $ly - $lx;
901   $scale += $diff if $diff > 0;         # if lx << ly, but not if ly << lx!
902     
903   # make copy of $x in case of list context for later reminder calculation
904   my $rem;
905   if (wantarray && !$y->is_one())
906     {
907     $rem = $x->copy();
908     }
909
910   $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+'; 
911
912   # check for / +-1 ( +/- 1E0)
913   if (!$y->is_one())
914     {
915     # promote BigInts and it's subclasses (except when already a BigFloat)
916     $y = $self->new($y) unless $y->isa('Math::BigFloat'); 
917
918     #print "bdiv $y ",ref($y),"\n";
919     # need to disable $upgrade in BigInt, to avoid deep recursion
920     local $Math::BigInt::upgrade = undef;       # should be parent class vs MBI
921
922     # calculate the result to $scale digits and then round it
923     # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
924     $x->{_m}->blsft($scale,10);
925     $x->{_m}->bdiv( $y->{_m} ); # a/c
926     $x->{_e}->bsub( $y->{_e} ); # b-d
927     $x->{_e}->bsub($scale);     # correct for 10**scale
928     $x->bnorm();                # remove trailing 0's
929     }
930
931   # shortcut to not run trough _find_round_parameters again
932   if (defined $params[1])
933     {
934     $x->bround($params[1],$params[3]);          # then round accordingly
935     }
936   else
937     {
938     $x->bfround($params[2],$params[3]);         # then round accordingly
939     }
940   if ($fallback)
941     {
942     # clear a/p after round, since user did not request it
943     $x->{_a} = undef; $x->{_p} = undef;
944     }
945   
946   if (wantarray)
947     {
948     if (!$y->is_one())
949       {
950       $rem->bmod($y,$params[1],$params[2],$params[3]);  # copy already done
951       }
952     else
953       {
954       $rem = $self->bzero();
955       }
956     if ($fallback)
957       {
958       # clear a/p after round, since user did not request it
959       $rem->{_a} = undef; $rem->{_p} = undef;
960       }
961     return ($x,$rem);
962     }
963   $x;
964   }
965
966 sub bmod 
967   {
968   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder 
969
970   # set up parameters
971   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
972   # objectify is costly, so avoid it
973   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
974     {
975     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
976     }
977
978   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
979     {
980     my ($d,$re) = $self->SUPER::_div_inf($x,$y);
981     $x->{sign} = $re->{sign};
982     $x->{_e} = $re->{_e};
983     $x->{_m} = $re->{_m};
984     return $x->round($a,$p,$r,$y);
985     } 
986   return $x->bnan() if $x->is_zero() && $y->is_zero();
987   return $x if $y->is_zero();
988   return $x->bnan() if $x->is_nan() || $y->is_nan();
989   return $x->bzero() if $y->is_one() || $x->is_zero();
990
991   # inf handling is missing here
992  
993   my $cmp = $x->bacmp($y);                      # equal or $x < $y?
994   return $x->bzero($a,$p) if $cmp == 0;         # $x == $y => result 0
995
996   # only $y of the operands negative? 
997   my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
998
999   $x->{sign} = $y->{sign};                              # calc sign first
1000   return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0;  # $x < $y => result $x
1001   
1002   my $ym = $y->{_m}->copy();
1003   
1004   # 2e1 => 20
1005   $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
1006  
1007   # if $y has digits after dot
1008   my $shifty = 0;                       # correct _e of $x by this
1009   if ($y->{_e}->{sign} eq '-')          # has digits after dot
1010     {
1011     # 123 % 2.5 => 1230 % 25 => 5 => 0.5
1012     $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
1013     $x->blsft($shifty,10);              # 123 => 1230, $y->{_m} is already 25
1014     }
1015   # $ym is now mantissa of $y based on exponent 0
1016
1017   my $shiftx = 0;                       # correct _e of $x by this
1018   if ($x->{_e}->{sign} eq '-')          # has digits after dot
1019     {
1020     # 123.4 % 20 => 1234 % 200
1021     $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
1022     $ym->blsft($shiftx,10);
1023     }
1024   # 123e1 % 20 => 1230 % 20
1025   if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
1026     {
1027     $x->{_m}->blsft($x->{_e},10);
1028     }
1029   $x->{_e} = $MBI->bzero() unless $x->{_e}->is_zero();
1030   
1031   $x->{_e}->bsub($shiftx) if $shiftx != 0;
1032   $x->{_e}->bsub($shifty) if $shifty != 0;
1033   
1034   # now mantissas are equalized, exponent of $x is adjusted, so calc result
1035
1036   $x->{_m}->bmod($ym);
1037
1038   $x->{sign} = '+' if $x->{_m}->is_zero();              # fix sign for -0
1039   $x->bnorm();
1040
1041   if ($neg != 0)        # one of them negative => correct in place
1042     {
1043     my $r = $y - $x;
1044     $x->{_m} = $r->{_m};
1045     $x->{_e} = $r->{_e};
1046     $x->{sign} = '+' if $x->{_m}->is_zero();            # fix sign for -0
1047     $x->bnorm();
1048     }
1049
1050   $x->round($a,$p,$r,$y);       # round and return
1051   }
1052
1053 sub bsqrt
1054   { 
1055   # calculate square root; this should probably
1056   # use a different test to see whether the accuracy we want is...
1057   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1058
1059   return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
1060   return $x if $x->{sign} eq '+inf';                              # +inf
1061   return $x if $x->is_zero() || $x->is_one();
1062
1063   # we need to limit the accuracy to protect against overflow
1064   my $fallback = 0;
1065   my $scale = 0;
1066   my @params = $x->_find_round_parameters($a,$p,$r);
1067
1068   # no rounding at all, so must use fallback
1069   if (scalar @params == 1)
1070     {
1071     # simulate old behaviour
1072     $params[1] = $self->div_scale();    # and round to it as accuracy
1073     $scale = $params[1]+4;              # at least four more for proper round
1074     $params[3] = $r;                    # round mode by caller or undef
1075     $fallback = 1;                      # to clear a/p afterwards
1076     }
1077   else
1078     {
1079     # the 4 below is empirical, and there might be cases where it is not
1080     # enough...
1081     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1082     }
1083
1084   # when user set globals, they would interfere with our calculation, so
1085   # disable them and later re-enable them
1086   no strict 'refs';
1087   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1088   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1089   # we also need to disable any set A or P on $x (_find_round_parameters took
1090   # them already into account), since these would interfere, too
1091   delete $x->{_a}; delete $x->{_p};
1092   # need to disable $upgrade in BigInt, to avoid deep recursion
1093   local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
1094
1095   my $xas = $x->as_number();
1096   my $gs = $xas->copy()->bsqrt();       # some guess
1097
1098 #  print "guess $gs\n";
1099   if (($x->{_e}->{sign} ne '-')         # guess can't be accurate if there are
1100                                         # digits after the dot
1101    && ($xas->bacmp($gs * $gs) == 0))    # guess hit the nail on the head?
1102     {
1103     # exact result
1104     $x->{_m} = $gs; $x->{_e} = $MBI->bzero(); $x->bnorm();
1105     # shortcut to not run trough _find_round_parameters again
1106     if (defined $params[1])
1107       {
1108       $x->bround($params[1],$params[3]);        # then round accordingly
1109       }
1110     else
1111       {
1112       $x->bfround($params[2],$params[3]);       # then round accordingly
1113       }
1114     if ($fallback)
1115       {
1116       # clear a/p after round, since user did not request it
1117       $x->{_a} = undef; $x->{_p} = undef;
1118       }
1119     # re-enable A and P, upgrade is taken care of by "local"
1120     ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
1121     return $x;
1122     }
1123   $gs = $self->new( $gs );              # BigInt to BigFloat
1124
1125   my $lx = $x->{_m}->length();
1126   $scale = $lx if $scale < $lx;
1127   my $e = $self->new("1E-$scale");      # make test variable
1128
1129   my $y = $x->copy();
1130   my $two = $self->new(2);
1131   my $diff = $e;
1132   # promote BigInts and it's subclasses (except when already a BigFloat)
1133   $y = $self->new($y) unless $y->isa('Math::BigFloat'); 
1134
1135   my $rem;
1136   while ($diff->bacmp($e) >= 0)
1137     {
1138     $rem = $y->copy()->bdiv($gs,$scale);
1139     $rem = $y->copy()->bdiv($gs,$scale)->badd($gs)->bdiv($two,$scale);
1140     $diff = $rem->copy()->bsub($gs);
1141     $gs = $rem->copy();
1142     }
1143   # copy over to modify $x
1144   $x->{_m} = $rem->{_m}; $x->{_e} = $rem->{_e};
1145   
1146   # shortcut to not run trough _find_round_parameters again
1147   if (defined $params[1])
1148     {
1149     $x->bround($params[1],$params[3]);          # then round accordingly
1150     }
1151   else
1152     {
1153     $x->bfround($params[2],$params[3]);         # then round accordingly
1154     }
1155   if ($fallback)
1156     {
1157     # clear a/p after round, since user did not request it
1158     $x->{_a} = undef; $x->{_p} = undef;
1159     }
1160   # restore globals
1161   $$abr = $ab; $$pbr = $pb;
1162   $x;
1163   }
1164
1165 sub bfac
1166   {
1167   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1168   # compute factorial numbers
1169   # modifies first argument
1170   my ($self,$x,@r) = objectify(1,@_);
1171
1172   return $x->bnan() 
1173     if (($x->{sign} ne '+') ||          # inf, NaN, <0 etc => NaN
1174      ($x->{_e}->{sign} ne '+'));        # digits after dot?
1175
1176   return $x->bone('+',@r) if $x->is_zero() || $x->is_one();     # 0 or 1 => 1
1177   
1178   # use BigInt's bfac() for faster calc
1179   $x->{_m}->blsft($x->{_e},10);         # un-norm m
1180   $x->{_e}->bzero();                    # norm $x again
1181   $x->{_m}->bfac();                     # factorial
1182   $x->bnorm()->round(@r);
1183   }
1184
1185 sub _pow2
1186   {
1187   # Calculate a power where $y is a non-integer, like 2 ** 0.5
1188   my ($x,$y,$a,$p,$r) = @_;
1189   my $self = ref($x);
1190   
1191   # we need to limit the accuracy to protect against overflow
1192   my $fallback = 0;
1193   my $scale = 0;
1194   my @params = $x->_find_round_parameters($a,$p,$r);
1195
1196   # no rounding at all, so must use fallback
1197   if (scalar @params == 1)
1198     {
1199     # simulate old behaviour
1200     $params[1] = $self->div_scale();    # and round to it as accuracy
1201     $scale = $params[1]+4;              # at least four more for proper round
1202     $params[3] = $r;                    # round mode by caller or undef
1203     $fallback = 1;                      # to clear a/p afterwards
1204     }
1205   else
1206     {
1207     # the 4 below is empirical, and there might be cases where it is not
1208     # enough...
1209     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1210     }
1211
1212   # when user set globals, they would interfere with our calculation, so
1213   # disable then and later re-enable them
1214   no strict 'refs';
1215   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1216   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1217   # we also need to disable any set A or P on $x (_find_round_parameters took
1218   # them already into account), since these would interfere, too
1219   delete $x->{_a}; delete $x->{_p};
1220   # need to disable $upgrade in BigInt, to avoid deep recursion
1221   local $Math::BigInt::upgrade = undef;
1222  
1223   # split the second argument into its integer and fraction part
1224   # we calculate the result then from these two parts, like in
1225   # 2 ** 2.4 == (2 ** 2) * (2 ** 0.4)
1226   my $c = $self->new($y->as_number());  # integer part
1227   my $d = $y-$c;                        # fractional part
1228   my $xc = $x->copy();                  # a temp. copy
1229   
1230   # now calculate binary fraction from the decimal fraction on the fly
1231   # f.i. 0.654:
1232   # 0.654 * 2 = 1.308 > 1 => 0.1        ( 1.308 - 1 = 0.308)
1233   # 0.308 * 2 = 0.616 < 1 => 0.10
1234   # 0.616 * 2 = 1.232 > 1 => 0.101      ( 1.232 - 1 = 0.232)
1235   # and so on...
1236   # The process stops when the result is exactly one, or when we have
1237   # enough accuracy
1238
1239   # From the binary fraction we calculate the result as follows:
1240   # we assume the fraction ends in 1, and we remove this one first.
1241   # For each digit after the dot, assume 1 eq R and 0 eq XR, where R means
1242   # take square root and X multiply with the original X. 
1243   
1244   my $i = 0;
1245   while ($i++ < 50)
1246     {
1247     $d->badd($d);                                               # * 2
1248     last if $d->is_one();                                       # == 1
1249     $x->bsqrt();                                                # 0
1250     if ($d > 1)
1251       {
1252       $x->bsqrt(); $x->bmul($xc); $d->bdec();                   # 1
1253       }
1254     }
1255   # assume fraction ends in 1
1256   $x->bsqrt();                                                  # 1
1257   if (!$c->is_one())
1258     {
1259     $x->bmul( $xc->bpow($c) );
1260     }
1261   elsif (!$c->is_zero())
1262     {
1263     $x->bmul( $xc );
1264     }
1265   # done
1266
1267   # shortcut to not run trough _find_round_parameters again
1268   if (defined $params[1])
1269     {
1270     $x->bround($params[1],$params[3]);          # then round accordingly
1271     }
1272   else
1273     {
1274     $x->bfround($params[2],$params[3]);         # then round accordingly
1275     }
1276   if ($fallback)
1277     {
1278     # clear a/p after round, since user did not request it
1279     $x->{_a} = undef; $x->{_p} = undef;
1280     }
1281   # restore globals
1282   $$abr = $ab; $$pbr = $pb;
1283   $x;
1284   }
1285
1286 sub _pow
1287   {
1288   # Calculate a power where $y is a non-integer, like 2 ** 0.5
1289   my ($x,$y,$a,$p,$r) = @_;
1290   my $self = ref($x);
1291
1292   # if $y == 0.5, it is sqrt($x)
1293   return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
1294
1295   # u = y * ln x
1296   #                _                             _
1297   # Taylor:       |    u     u^2      u^3         |
1298   # x ** y  = 1 + |   --- +  --- + * ----- + ...  |
1299   #               |_   1     1*2     1*2*3       _|
1300
1301   # we need to limit the accuracy to protect against overflow
1302   my $fallback = 0;
1303   my $scale = 0;
1304   my @params = $x->_find_round_parameters($a,$p,$r);
1305
1306   # no rounding at all, so must use fallback
1307   if (scalar @params == 1)
1308     {
1309     # simulate old behaviour
1310     $params[1] = $self->div_scale();    # and round to it as accuracy
1311     $scale = $params[1]+4;              # at least four more for proper round
1312     $params[3] = $r;                    # round mode by caller or undef
1313     $fallback = 1;                      # to clear a/p afterwards
1314     }
1315   else
1316     {
1317     # the 4 below is empirical, and there might be cases where it is not
1318     # enough...
1319     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
1320     }
1321
1322   # when user set globals, they would interfere with our calculation, so
1323   # disable then and later re-enable them
1324   no strict 'refs';
1325   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1326   my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1327   # we also need to disable any set A or P on $x (_find_round_parameters took
1328   # them already into account), since these would interfere, too
1329   delete $x->{_a}; delete $x->{_p};
1330   # need to disable $upgrade in BigInt, to avoid deep recursion
1331   local $Math::BigInt::upgrade = undef;
1332  
1333   my ($limit,$v,$u,$below,$factor,$next,$over);
1334
1335   $u = $x->copy()->blog($scale)->bmul($y);
1336   $v = $self->bone();                           # 1
1337   $factor = $self->new(2);                      # 2
1338   $x->bone();                                   # first term: 1
1339
1340   $below = $v->copy();
1341   $over = $u->copy();
1342  
1343   $limit = $self->new("1E-". ($scale-1));
1344   #my $steps = 0;
1345   while (3 < 5)
1346     {
1347     # we calculate the next term, and add it to the last
1348     # when the next term is below our limit, it won't affect the outcome
1349     # anymore, so we stop
1350     $next = $over->copy()->bdiv($below,$scale);
1351     last if $next->bcmp($limit) <= 0;
1352     $x->badd($next);
1353 #    print "at $x\n";
1354     # calculate things for the next term
1355     $over *= $u; $below *= $factor; $factor->binc();
1356     #$steps++;
1357     }
1358   
1359   # shortcut to not run trough _find_round_parameters again
1360   if (defined $params[1])
1361     {
1362     $x->bround($params[1],$params[3]);          # then round accordingly
1363     }
1364   else
1365     {
1366     $x->bfround($params[2],$params[3]);         # then round accordingly
1367     }
1368   if ($fallback)
1369     {
1370     # clear a/p after round, since user did not request it
1371     $x->{_a} = undef; $x->{_p} = undef;
1372     }
1373   # restore globals
1374   $$abr = $ab; $$pbr = $pb;
1375   $x;
1376   }
1377
1378 sub bpow 
1379   {
1380   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1381   # compute power of two numbers, second arg is used as integer
1382   # modifies first argument
1383
1384   # set up parameters
1385   my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1386   # objectify is costly, so avoid it
1387   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1388     {
1389     ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1390     }
1391
1392   return $x if $x->{sign} =~ /^[+-]inf$/;
1393   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1394   return $x->bone() if $y->is_zero();
1395   return $x         if $x->is_one() || $y->is_one();
1396
1397   return $x->_pow($y,$a,$p,$r) if !$y->is_int();        # non-integer power
1398
1399   my $y1 = $y->as_number();             # make bigint
1400   # if ($x == -1)
1401   if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
1402     {
1403     # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
1404     return $y1->is_odd() ? $x : $x->babs(1);
1405     }
1406   if ($x->is_zero())
1407     {
1408     return $x if $y->{sign} eq '+';     # 0**y => 0 (if not y <= 0)
1409     # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1410     $x->binf();
1411     }
1412
1413   # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1414   $y1->babs();
1415   $x->{_m}->bpow($y1);
1416   $x->{_e}->bmul($y1);
1417   $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1418   $x->bnorm();
1419   if ($y->{sign} eq '-')
1420     {
1421     # modify $x in place!
1422     my $z = $x->copy(); $x->bzero()->binc();
1423     return $x->bdiv($z,$a,$p,$r);       # round in one go (might ignore y's A!)
1424     }
1425   $x->round($a,$p,$r,$y);
1426   }
1427
1428 ###############################################################################
1429 # rounding functions
1430
1431 sub bfround
1432   {
1433   # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1434   # $n == 0 means round to integer
1435   # expects and returns normalized numbers!
1436   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1437
1438   return $x if $x->modify('bfround');
1439   
1440   my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1441   return $x if !defined $scale;                 # no-op
1442
1443   # never round a 0, +-inf, NaN
1444   if ($x->is_zero())
1445     {
1446     $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1447     return $x; 
1448     }
1449   return $x if $x->{sign} !~ /^[+-]$/;
1450
1451   # don't round if x already has lower precision
1452   return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1453
1454   $x->{_p} = $scale;                    # remember round in any case
1455   $x->{_a} = undef;                     # and clear A
1456   if ($scale < 0)
1457     {
1458     # round right from the '.'
1459
1460     return $x if $x->{_e}->{sign} eq '+';       # e >= 0 => nothing to round
1461
1462     $scale = -$scale;                           # positive for simplicity
1463     my $len = $x->{_m}->length();               # length of mantissa
1464
1465     # the following poses a restriction on _e, but if _e is bigger than a
1466     # scalar, you got other problems (memory etc) anyway
1467     my $dad = -($x->{_e}->numify());            # digits after dot
1468     my $zad = 0;                                # zeros after dot
1469     $zad = $dad - $len if (-$dad < -$len);      # for 0.00..00xxx style
1470     
1471     #print "scale $scale dad $dad zad $zad len $len\n";
1472     # number  bsstr   len zad dad       
1473     # 0.123   123e-3    3   0 3
1474     # 0.0123  123e-4    3   1 4
1475     # 0.001   1e-3      1   2 3
1476     # 1.23    123e-2    3   0 2
1477     # 1.2345  12345e-4  5   0 4
1478
1479     # do not round after/right of the $dad
1480     return $x if $scale > $dad;                 # 0.123, scale >= 3 => exit
1481
1482     # round to zero if rounding inside the $zad, but not for last zero like:
1483     # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1484     return $x->bzero() if $scale < $zad;
1485     if ($scale == $zad)                 # for 0.006, scale -3 and trunc
1486       {
1487       $scale = -$len;
1488       }
1489     else
1490       {
1491       # adjust round-point to be inside mantissa
1492       if ($zad != 0)
1493         {
1494         $scale = $scale-$zad;
1495         }
1496       else
1497         {
1498         my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;    # digits before dot
1499         $scale = $dbd+$scale;
1500         }
1501       }
1502     }
1503   else
1504     {
1505     # round left from the '.'
1506
1507     # 123 => 100 means length(123) = 3 - $scale (2) => 1
1508
1509     my $dbt = $x->{_m}->length(); 
1510     # digits before dot 
1511     my $dbd = $dbt + $x->{_e}->numify(); 
1512     # should be the same, so treat it as this 
1513     $scale = 1 if $scale == 0; 
1514     # shortcut if already integer 
1515     return $x if $scale == 1 && $dbt <= $dbd; 
1516     # maximum digits before dot 
1517     ++$dbd;
1518
1519     if ($scale > $dbd) 
1520        { 
1521        # not enough digits before dot, so round to zero 
1522        return $x->bzero; 
1523        }
1524     elsif ( $scale == $dbd )
1525        { 
1526        # maximum 
1527        $scale = -$dbt; 
1528        } 
1529     else
1530        { 
1531        $scale = $dbd - $scale; 
1532        }
1533     }
1534   # pass sign to bround for rounding modes '+inf' and '-inf'
1535   $x->{_m}->{sign} = $x->{sign};
1536   $x->{_m}->bround($scale,$mode);
1537   $x->{_m}->{sign} = '+';               # fix sign back
1538   $x->bnorm();
1539   }
1540
1541 sub bround
1542   {
1543   # accuracy: preserve $N digits, and overwrite the rest with 0's
1544   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1545   
1546   die ('bround() needs positive accuracy') if ($_[0] || 0) < 0;
1547
1548   my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1549   return $x if !defined $scale;                         # no-op
1550
1551   return $x if $x->modify('bround');
1552
1553   # scale is now either $x->{_a}, $accuracy, or the user parameter
1554   # test whether $x already has lower accuracy, do nothing in this case 
1555   # but do round if the accuracy is the same, since a math operation might
1556   # want to round a number with A=5 to 5 digits afterwards again
1557   return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
1558
1559   # scale < 0 makes no sense
1560   # never round a +-inf, NaN
1561   return $x if ($scale < 0) ||  $x->{sign} !~ /^[+-]$/;
1562
1563   # 1: $scale == 0 => keep all digits
1564   # 2: never round a 0
1565   # 3: if we should keep more digits than the mantissa has, do nothing
1566   if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1567     {
1568     $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1569     return $x; 
1570     }
1571
1572   # pass sign to bround for '+inf' and '-inf' rounding modes
1573   $x->{_m}->{sign} = $x->{sign};
1574   $x->{_m}->bround($scale,$mode);       # round mantissa
1575   $x->{_m}->{sign} = '+';               # fix sign back
1576   # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
1577   $x->{_a} = $scale;                    # remember rounding
1578   $x->{_p} = undef;                     # and clear P
1579   $x->bnorm();                          # del trailing zeros gen. by bround()
1580   }
1581
1582 sub bfloor
1583   {
1584   # return integer less or equal then $x
1585   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1586
1587   return $x if $x->modify('bfloor');
1588    
1589   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1590
1591   # if $x has digits after dot
1592   if ($x->{_e}->{sign} eq '-')
1593     {
1594     $x->{_e}->{sign} = '+';                     # negate e
1595     $x->{_m}->brsft($x->{_e},10);               # cut off digits after dot
1596     $x->{_e}->bzero();                          # trunc/norm    
1597     $x->{_m}->binc() if $x->{sign} eq '-';      # decrement if negative
1598     }
1599   $x->round($a,$p,$r);
1600   }
1601
1602 sub bceil
1603   {
1604   # return integer greater or equal then $x
1605   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1606
1607   return $x if $x->modify('bceil');
1608   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1609
1610   # if $x has digits after dot
1611   if ($x->{_e}->{sign} eq '-')
1612     {
1613     #$x->{_m}->brsft(-$x->{_e},10);
1614     #$x->{_e}->bzero();
1615     #$x++ if $x->{sign} eq '+';
1616
1617     $x->{_e}->{sign} = '+';                     # negate e
1618     $x->{_m}->brsft($x->{_e},10);               # cut off digits after dot
1619     $x->{_e}->bzero();                          # trunc/norm    
1620     $x->{_m}->binc() if $x->{sign} eq '+';      # decrement if negative
1621     }
1622   $x->round($a,$p,$r);
1623   }
1624
1625 sub brsft
1626   {
1627   # shift right by $y (divide by power of $n)
1628   
1629   # set up parameters
1630   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1631   # objectify is costly, so avoid it
1632   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1633     {
1634     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1635     }
1636
1637   return $x if $x->modify('brsft');
1638   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1639
1640   $n = 2 if !defined $n; $n = $self->new($n);
1641   $x->bdiv($n->bpow($y),$a,$p,$r,$y);
1642   }
1643
1644 sub blsft
1645   {
1646   # shift left by $y (multiply by power of $n)
1647   
1648   # set up parameters
1649   my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1650   # objectify is costly, so avoid it
1651   if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1652     {
1653     ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1654     }
1655
1656   return $x if $x->modify('blsft');
1657   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1658
1659   $n = 2 if !defined $n; $n = $self->new($n);
1660   $x->bmul($n->bpow($y),$a,$p,$r,$y);
1661   }
1662
1663 ###############################################################################
1664
1665 sub DESTROY
1666   {
1667   # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
1668   }
1669
1670 sub AUTOLOAD
1671   {
1672   # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
1673   # or falling back to MBI::bxxx()
1674   my $name = $AUTOLOAD;
1675
1676   $name =~ s/.*:://;    # split package
1677   no strict 'refs';
1678   if (!method_alias($name))
1679     {
1680     if (!defined $name)
1681       {
1682       # delayed load of Carp and avoid recursion        
1683       require Carp;
1684       Carp::croak ("Can't call a method without name");
1685       }
1686     if (!method_hand_up($name))
1687       {
1688       # delayed load of Carp and avoid recursion        
1689       require Carp;
1690       Carp::croak ("Can't call $class\-\>$name, not a valid method");
1691       }
1692     # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1693     $name =~ s/^f/b/;
1694     return &{"$MBI"."::$name"}(@_);
1695     }
1696   my $bname = $name; $bname =~ s/^f/b/;
1697   *{$class."::$name"} = \&$bname;
1698   &$bname;      # uses @_
1699   }
1700
1701 sub exponent
1702   {
1703   # return a copy of the exponent
1704   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1705
1706   if ($x->{sign} !~ /^[+-]$/)
1707     {
1708     my $s = $x->{sign}; $s =~ s/^[+-]//;
1709     return $self->new($s);                      # -inf, +inf => +inf
1710     }
1711   return $x->{_e}->copy();
1712   }
1713
1714 sub mantissa
1715   {
1716   # return a copy of the mantissa
1717   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1718  
1719   if ($x->{sign} !~ /^[+-]$/)
1720     {
1721     my $s = $x->{sign}; $s =~ s/^[+]//;
1722     return $self->new($s);                      # -inf, +inf => +inf
1723     }
1724   my $m = $x->{_m}->copy();             # faster than going via bstr()
1725   $m->bneg() if $x->{sign} eq '-';
1726
1727   $m;
1728   }
1729
1730 sub parts
1731   {
1732   # return a copy of both the exponent and the mantissa
1733   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1734
1735   if ($x->{sign} !~ /^[+-]$/)
1736     {
1737     my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
1738     return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
1739     }
1740   my $m = $x->{_m}->copy();     # faster than going via bstr()
1741   $m->bneg() if $x->{sign} eq '-';
1742   return ($m,$x->{_e}->copy());
1743   }
1744
1745 ##############################################################################
1746 # private stuff (internal use only)
1747
1748 sub import
1749   {
1750   my $self = shift;
1751   my $l = scalar @_;
1752   my $lib = ''; my @a;
1753   for ( my $i = 0; $i < $l ; $i++)
1754     {
1755 #    print "at $_[$i] (",$_[$i+1]||'undef',")\n";
1756     if ( $_[$i] eq ':constant' )
1757       {
1758       # this rest causes overlord er load to step in
1759       # print "overload @_\n";
1760       overload::constant float => sub { $self->new(shift); }; 
1761       }
1762     elsif ($_[$i] eq 'upgrade')
1763       {
1764       # this causes upgrading
1765       $upgrade = $_[$i+1];              # or undef to disable
1766       $i++;
1767       }
1768     elsif ($_[$i] eq 'downgrade')
1769       {
1770       # this causes downgrading
1771       $downgrade = $_[$i+1];            # or undef to disable
1772       $i++;
1773       }
1774     elsif ($_[$i] eq 'lib')
1775       {
1776       $lib = $_[$i+1] || '';            # default Calc
1777       $i++;
1778       }
1779     elsif ($_[$i] eq 'with')
1780       {
1781       $MBI = $_[$i+1] || 'Math::BigInt';        # default Math::BigInt
1782       $i++;
1783       }
1784     else
1785       {
1786       push @a, $_[$i];
1787       }
1788     }
1789
1790   # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
1791   my $mbilib = eval { Math::BigInt->config()->{lib} };
1792   if ((defined $mbilib) && ($MBI eq 'Math::BigInt'))
1793     {
1794     # MBI already loaded
1795     $MBI->import('lib',"$lib,$mbilib", 'objectify');
1796     }
1797   else
1798     {
1799     # MBI not loaded, or with ne "Math::BigInt"
1800     $lib .= ",$mbilib" if defined $mbilib;
1801     $lib =~ s/^,//;                             # don't leave empty 
1802     if ($] < 5.006)
1803       {
1804       # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
1805       # used in the same script, or eval inside import().
1806       my @parts = split /::/, $MBI;             # Math::BigInt => Math BigInt
1807       my $file = pop @parts; $file .= '.pm';    # BigInt => BigInt.pm
1808       require File::Spec;
1809       $file = File::Spec->catfile (@parts, $file);
1810       eval { require "$file"; };
1811       $MBI->import( lib => $lib, 'objectify' );
1812       }
1813     else
1814       {
1815       my $rc = "use $MBI lib => '$lib', 'objectify';";
1816       eval $rc;
1817       }
1818     }
1819   die ("Couldn't load $MBI: $! $@") if $@;
1820
1821   # any non :constant stuff is handled by our parent, Exporter
1822   # even if @_ is empty, to give it a chance
1823   $self->SUPER::import(@a);             # for subclasses
1824   $self->export_to_level(1,$self,@a);   # need this, too
1825   }
1826
1827 sub bnorm
1828   {
1829   # adjust m and e so that m is smallest possible
1830   # round number according to accuracy and precision settings
1831   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1832
1833   return $x if $x->{sign} !~ /^[+-]$/;          # inf, nan etc
1834
1835 #  if (!$x->{_m}->is_odd())
1836 #    {
1837     my $zeros = $x->{_m}->_trailing_zeros();    # correct for trailing zeros 
1838     if ($zeros != 0)
1839       {
1840       $x->{_m}->brsft($zeros,10); $x->{_e}->badd($zeros);
1841       }
1842     # for something like 0Ey, set y to 1, and -0 => +0
1843     $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
1844 #    }
1845   # this is to prevent automatically rounding when MBI's globals are set
1846   $x->{_m}->{_f} = MB_NEVER_ROUND;
1847   $x->{_e}->{_f} = MB_NEVER_ROUND;
1848   # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
1849   $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
1850   $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
1851   $x;                                   # MBI bnorm is no-op, so dont call it
1852   } 
1853  
1854 ##############################################################################
1855 # internal calculation routines
1856
1857 sub as_number
1858   {
1859   # return copy as a bigint representation of this BigFloat number
1860   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1861
1862   my $z = $x->{_m}->copy();
1863   if ($x->{_e}->{sign} eq '-')          # < 0
1864     {
1865     $x->{_e}->{sign} = '+';             # flip
1866     $z->brsft($x->{_e},10);
1867     $x->{_e}->{sign} = '-';             # flip back
1868     } 
1869   elsif (!$x->{_e}->is_zero())          # > 0 
1870     {
1871     $z->blsft($x->{_e},10);
1872     }
1873   $z->{sign} = $x->{sign};
1874   $z;
1875   }
1876
1877 sub length
1878   {
1879   my $x = shift;
1880   my $class = ref($x) || $x;
1881   $x = $class->new(shift) unless ref($x);
1882
1883   return 1 if $x->{_m}->is_zero();
1884   my $len = $x->{_m}->length();
1885   $len += $x->{_e} if $x->{_e}->sign() eq '+';
1886   if (wantarray())
1887     {
1888     my $t = $MBI->bzero();
1889     $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1890     return ($len,$t);
1891     }
1892   $len;
1893   }
1894
1895 1;
1896 __END__
1897
1898 =head1 NAME
1899
1900 Math::BigFloat - Arbitrary size floating point math package
1901
1902 =head1 SYNOPSIS
1903
1904   use Math::BigFloat;
1905
1906   # Number creation
1907   $x = Math::BigFloat->new($str);       # defaults to 0
1908   $nan  = Math::BigFloat->bnan();       # create a NotANumber
1909   $zero = Math::BigFloat->bzero();      # create a +0
1910   $inf = Math::BigFloat->binf();        # create a +inf
1911   $inf = Math::BigFloat->binf('-');     # create a -inf
1912   $one = Math::BigFloat->bone();        # create a +1
1913   $one = Math::BigFloat->bone('-');     # create a -1
1914
1915   # Testing
1916   $x->is_zero();                # true if arg is +0
1917   $x->is_nan();                 # true if arg is NaN
1918   $x->is_one();                 # true if arg is +1
1919   $x->is_one('-');              # true if arg is -1
1920   $x->is_odd();                 # true if odd, false for even
1921   $x->is_even();                # true if even, false for odd
1922   $x->is_positive();            # true if >= 0
1923   $x->is_negative();            # true if <  0
1924   $x->is_inf(sign);             # true if +inf, or -inf (default is '+')
1925
1926   $x->bcmp($y);                 # compare numbers (undef,<0,=0,>0)
1927   $x->bacmp($y);                # compare absolutely (undef,<0,=0,>0)
1928   $x->sign();                   # return the sign, either +,- or NaN
1929   $x->digit($n);                # return the nth digit, counting from right
1930   $x->digit(-$n);               # return the nth digit, counting from left 
1931
1932   # The following all modify their first argument:
1933   
1934   # set 
1935   $x->bzero();                  # set $i to 0
1936   $x->bnan();                   # set $i to NaN
1937   $x->bone();                   # set $x to +1
1938   $x->bone('-');                # set $x to -1
1939   $x->binf();                   # set $x to inf
1940   $x->binf('-');                # set $x to -inf
1941
1942   $x->bneg();                   # negation
1943   $x->babs();                   # absolute value
1944   $x->bnorm();                  # normalize (no-op)
1945   $x->bnot();                   # two's complement (bit wise not)
1946   $x->binc();                   # increment x by 1
1947   $x->bdec();                   # decrement x by 1
1948   
1949   $x->badd($y);                 # addition (add $y to $x)
1950   $x->bsub($y);                 # subtraction (subtract $y from $x)
1951   $x->bmul($y);                 # multiplication (multiply $x by $y)
1952   $x->bdiv($y);                 # divide, set $i to quotient
1953                                 # return (quo,rem) or quo if scalar
1954
1955   $x->bmod($y);                 # modulus
1956   $x->bpow($y);                 # power of arguments (a**b)
1957   $x->blsft($y);                # left shift
1958   $x->brsft($y);                # right shift 
1959                                 # return (quo,rem) or quo if scalar
1960   
1961   $x->blog($base);              # logarithm of $x, base defaults to e
1962                                 # (other bases than e not supported yet)
1963   
1964   $x->band($y);                 # bit-wise and
1965   $x->bior($y);                 # bit-wise inclusive or
1966   $x->bxor($y);                 # bit-wise exclusive or
1967   $x->bnot();                   # bit-wise not (two's complement)
1968  
1969   $x->bsqrt();                  # calculate square-root
1970   $x->bfac();                   # factorial of $x (1*2*3*4*..$x)
1971  
1972   $x->bround($N);               # accuracy: preserver $N digits
1973   $x->bfround($N);              # precision: round to the $Nth digit
1974
1975   # The following do not modify their arguments:
1976   bgcd(@values);                # greatest common divisor
1977   blcm(@values);                # lowest common multiplicator
1978   
1979   $x->bstr();                   # return string
1980   $x->bsstr();                  # return string in scientific notation
1981  
1982   $x->bfloor();                 # return integer less or equal than $x
1983   $x->bceil();                  # return integer greater or equal than $x
1984  
1985   $x->exponent();               # return exponent as BigInt
1986   $x->mantissa();               # return mantissa as BigInt
1987   $x->parts();                  # return (mantissa,exponent) as BigInt
1988
1989   $x->length();                 # number of digits (w/o sign and '.')
1990   ($l,$f) = $x->length();       # number of digits, and length of fraction      
1991
1992   $x->precision();              # return P of $x (or global, if P of $x undef)
1993   $x->precision($n);            # set P of $x to $n
1994   $x->accuracy();               # return A of $x (or global, if A of $x undef)
1995   $x->accuracy($n);             # set A $x to $n
1996
1997   Math::BigFloat->precision();  # get/set global P for all BigFloat objects
1998   Math::BigFloat->accuracy();   # get/set global A for all BigFloat objects
1999
2000 =head1 DESCRIPTION
2001
2002 All operators (inlcuding basic math operations) are overloaded if you
2003 declare your big floating point numbers as
2004
2005   $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2006
2007 Operations with overloaded operators preserve the arguments, which is
2008 exactly what you expect.
2009
2010 =head2 Canonical notation
2011
2012 Input to these routines are either BigFloat objects, or strings of the
2013 following four forms:
2014
2015 =over 2
2016
2017 =item *
2018
2019 C</^[+-]\d+$/>
2020
2021 =item *
2022
2023 C</^[+-]\d+\.\d*$/>
2024
2025 =item *
2026
2027 C</^[+-]\d+E[+-]?\d+$/>
2028
2029 =item *
2030
2031 C</^[+-]\d*\.\d+E[+-]?\d+$/>
2032
2033 =back
2034
2035 all with optional leading and trailing zeros and/or spaces. Additonally,
2036 numbers are allowed to have an underscore between any two digits.
2037
2038 Empty strings as well as other illegal numbers results in 'NaN'.
2039
2040 bnorm() on a BigFloat object is now effectively a no-op, since the numbers 
2041 are always stored in normalized form. On a string, it creates a BigFloat 
2042 object.
2043
2044 =head2 Output
2045
2046 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2047
2048 The string output will always have leading and trailing zeros stripped and drop
2049 a plus sign. C<bstr()> will give you always the form with a decimal point,
2050 while C<bsstr()> (for scientific) gives you the scientific notation.
2051
2052         Input                   bstr()          bsstr()
2053         '-0'                    '0'             '0E1'
2054         '  -123 123 123'        '-123123123'    '-123123123E0'
2055         '00.0123'               '0.0123'        '123E-4'
2056         '123.45E-2'             '1.2345'        '12345E-4'
2057         '10E+3'                 '10000'         '1E4'
2058
2059 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2060 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2061 return either undef, <0, 0 or >0 and are suited for sort.
2062
2063 Actual math is done by using BigInts to represent the mantissa and exponent.
2064 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to 
2065 represent the result when input arguments are not numbers, as well as 
2066 the result of dividing by zero.
2067
2068 =head2 C<mantissa()>, C<exponent()> and C<parts()>
2069
2070 C<mantissa()> and C<exponent()> return the said parts of the BigFloat 
2071 as BigInts such that:
2072
2073         $m = $x->mantissa();
2074         $e = $x->exponent();
2075         $y = $m * ( 10 ** $e );
2076         print "ok\n" if $x == $y;
2077
2078 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2079
2080 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2081
2082 Currently the mantissa is reduced as much as possible, favouring higher
2083 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2084 This might change in the future, so do not depend on it.
2085
2086 =head2 Accuracy vs. Precision
2087
2088 See also: L<Rounding|Rounding>.
2089
2090 Math::BigFloat supports both precision and accuracy. For a full documentation,
2091 examples and tips on these topics please see the large section in
2092 L<Math::BigInt>.
2093
2094 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2095 a operation consumes all resources, each operation produces no more than
2096 C<Math::BigFloat::precision()> digits.
2097
2098 In case the result of one operation has more precision than specified,
2099 it is rounded. The rounding mode taken is either the default mode, or the one
2100 supplied to the operation after the I<scale>:
2101
2102         $x = Math::BigFloat->new(2);
2103         Math::BigFloat::precision(5);           # 5 digits max
2104         $y = $x->copy()->bdiv(3);               # will give 0.66666
2105         $y = $x->copy()->bdiv(3,6);             # will give 0.666666
2106         $y = $x->copy()->bdiv(3,6,'odd');       # will give 0.666667
2107         Math::BigFloat::round_mode('zero');
2108         $y = $x->copy()->bdiv(3,6);             # will give 0.666666
2109
2110 =head2 Rounding
2111
2112 =over 2
2113
2114 =item ffround ( +$scale )
2115
2116 Rounds to the $scale'th place left from the '.', counting from the dot.
2117 The first digit is numbered 1. 
2118
2119 =item ffround ( -$scale )
2120
2121 Rounds to the $scale'th place right from the '.', counting from the dot.
2122
2123 =item ffround ( 0 )
2124
2125 Rounds to an integer.
2126
2127 =item fround  ( +$scale )
2128
2129 Preserves accuracy to $scale digits from the left (aka significant digits)
2130 and pads the rest with zeros. If the number is between 1 and -1, the
2131 significant digits count from the first non-zero after the '.'
2132
2133 =item fround  ( -$scale ) and fround ( 0 )
2134
2135 These are effetively no-ops.
2136
2137 =back
2138
2139 All rounding functions take as a second parameter a rounding mode from one of
2140 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
2141
2142 The default rounding mode is 'even'. By using
2143 C<< Math::BigFloat::round_mode($round_mode); >> you can get and set the default
2144 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
2145 no longer supported.
2146 The second parameter to the round functions then overrides the default
2147 temporarily. 
2148
2149 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
2150 'trunc' as rounding mode to make it equivalent to:
2151
2152         $x = 2.5;
2153         $y = int($x) + 2;
2154
2155 You can override this by passing the desired rounding mode as parameter to
2156 C<as_number()>:
2157
2158         $x = Math::BigFloat->new(2.5);
2159         $y = $x->as_number('odd');      # $y = 3
2160
2161 =head1 EXAMPLES
2162  
2163   # not ready yet
2164
2165 =head1 Autocreating constants
2166
2167 After C<use Math::BigFloat ':constant'> all the floating point constants
2168 in the given scope are converted to C<Math::BigFloat>. This conversion
2169 happens at compile time.
2170
2171 In particular
2172
2173   perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2174
2175 prints the value of C<2E-100>. Note that without conversion of 
2176 constants the expression 2E-100 will be calculated as normal floating point 
2177 number.
2178
2179 Please note that ':constant' does not affect integer constants, nor binary 
2180 nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2181 work.
2182
2183 =head2 Math library
2184
2185 Math with the numbers is done (by default) by a module called
2186 Math::BigInt::Calc. This is equivalent to saying:
2187
2188         use Math::BigFloat lib => 'Calc';
2189
2190 You can change this by using:
2191
2192         use Math::BigFloat lib => 'BitVect';
2193
2194 The following would first try to find Math::BigInt::Foo, then
2195 Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2196
2197         use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2198
2199 Calc.pm uses as internal format an array of elements of some decimal base
2200 (usually 1e7, but this might be differen for some systems) with the least
2201 significant digit first, while BitVect.pm uses a bit vector of base 2, most
2202 significant bit first. Other modules might use even different means of
2203 representing the numbers. See the respective module documentation for further
2204 details.
2205
2206 Please note that Math::BigFloat does B<not> use the denoted library itself,
2207 but it merely passes the lib argument to Math::BigInt. So, instead of the need
2208 to do:
2209
2210         use Math::BigInt lib => 'GMP';
2211         use Math::BigFloat;
2212
2213 you can roll it all into one line:
2214
2215         use Math::BigFloat lib => 'GMP';
2216
2217 Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details.
2218
2219 =head2 Using Math::BigInt::Lite
2220
2221 It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2222
2223         # 1
2224         use Math::BigFloat with => 'Math::BigInt::Lite';
2225
2226 There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2227 can combine these if you want. For instance, you may want to use
2228 Math::BigInt objects in your main script, too.
2229
2230         # 2
2231         use Math::BigInt;
2232         use Math::BigFloat with => 'Math::BigInt::Lite';
2233
2234 Of course, you can combine this with the C<lib> parameter.
2235
2236         # 3
2237         use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2238
2239 If you want to use Math::BigInt's, too, simple add a Math::BigInt B<before>:
2240
2241         # 4
2242         use Math::BigInt;
2243         use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2244
2245 Notice that the module with the last C<lib> will "win" and thus
2246 it's lib will be used if the lib is available:
2247
2248         # 5
2249         use Math::BigInt lib => 'Bar,Baz';
2250         use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2251
2252 That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2253 words, Math::BigFloat will try to retain previously loaded libs when you
2254 don't specify it one.
2255
2256 Actually, the lib loading order would be "Bar,Baz,Calc", and then
2257 "Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
2258 same as trying the latter load alone, except for the fact that Bar or Baz
2259 might be loaded needlessly in an intermidiate step
2260
2261 The old way still works though:
2262
2263         # 6
2264         use Math::BigInt lib => 'Bar,Baz';
2265         use Math::BigFloat;
2266
2267 But B<examples #3 and #4 are recommended> for usage.
2268
2269 =head1 BUGS
2270
2271 =over 2
2272
2273 =item *
2274
2275 The following does not work yet:
2276
2277         $m = $x->mantissa();
2278         $e = $x->exponent();
2279         $y = $m * ( 10 ** $e );
2280         print "ok\n" if $x == $y;
2281
2282 =item *
2283
2284 There is no fmod() function yet.
2285
2286 =back
2287
2288 =head1 CAVEAT
2289
2290 =over 1
2291
2292 =item stringify, bstr()
2293
2294 Both stringify and bstr() now drop the leading '+'. The old code would return
2295 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2296 reasoning and details.
2297
2298 =item bdiv
2299
2300 The following will probably not do what you expect:
2301
2302         print $c->bdiv(123.456),"\n";
2303
2304 It prints both quotient and reminder since print works in list context. Also,
2305 bdiv() will modify $c, so be carefull. You probably want to use
2306         
2307         print $c / 123.456,"\n";
2308         print scalar $c->bdiv(123.456),"\n";  # or if you want to modify $c
2309
2310 instead.
2311
2312 =item Modifying and =
2313
2314 Beware of:
2315
2316         $x = Math::BigFloat->new(5);
2317         $y = $x;
2318
2319 It will not do what you think, e.g. making a copy of $x. Instead it just makes
2320 a second reference to the B<same> object and stores it in $y. Thus anything
2321 that modifies $x will modify $y, and vice versa.
2322
2323         $x->bmul(2);
2324         print "$x, $y\n";       # prints '10, 10'
2325
2326 If you want a true copy of $x, use:
2327         
2328         $y = $x->copy();
2329
2330 See also the documentation in L<overload> regarding C<=>.
2331
2332 =item bpow
2333
2334 C<bpow()> now modifies the first argument, unlike the old code which left
2335 it alone and only returned the result. This is to be consistent with
2336 C<badd()> etc. The first will modify $x, the second one won't:
2337
2338         print bpow($x,$i),"\n";         # modify $x
2339         print $x->bpow($i),"\n";        # ditto
2340         print $x ** $i,"\n";            # leave $x alone 
2341
2342 =back
2343
2344 =head1 LICENSE
2345
2346 This program is free software; you may redistribute it and/or modify it under
2347 the same terms as Perl itself.
2348
2349 =head1 AUTHORS
2350
2351 Mark Biggar, overloaded interface by Ilya Zakharevich.
2352 Completely rewritten by Tels http://bloodgate.com in 2001.
2353
2354 =cut