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