This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Math-BigInt v1.49 released
[perl5.git] / lib / Math / BigFloat.pm
1 #!/usr/bin/perl -w
2
3 # The following hash values are internally used:
4 #   _e: exponent (BigInt)
5 #   _m: mantissa (absolute BigInt)
6 # sign: +,-,"NaN" if not a number
7 #   _a: accuracy
8 #   _p: precision
9 #   _f: flags, used to signal MBI not to touch our private parts
10
11 package Math::BigFloat;
12
13 $VERSION = '1.27';
14 require 5.005;
15 use Exporter;
16 use Math::BigInt qw/objectify/;
17 @ISA =       qw( Exporter Math::BigInt);
18
19 use strict;
20 use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode/;
21 my $class = "Math::BigFloat";
22
23 use overload
24 '<=>'   =>      sub { $_[2] ?
25                       ref($_[0])->bcmp($_[1],$_[0]) : 
26                       ref($_[0])->bcmp($_[0],$_[1])},
27 'int'   =>      sub { $_[0]->as_number() },             # 'trunc' to bigint
28 ;
29
30 ##############################################################################
31 # global constants, flags and accessory
32
33 use constant MB_NEVER_ROUND => 0x0001;
34
35 # are NaNs ok?
36 my $NaNOK=1;
37 # constant for easier life
38 my $nan = 'NaN'; 
39
40 # class constants, use Class->constant_name() to access
41 $round_mode = 'even'; # one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'
42 $accuracy   = undef;
43 $precision  = undef;
44 $div_scale  = 40;
45
46 ##############################################################################
47 # the old code had $rnd_mode, so we need to support it, too
48
49 $rnd_mode   = 'even';
50 sub TIESCALAR   { my ($class) = @_; bless \$round_mode, $class; }
51 sub FETCH       { return $round_mode; }
52 sub STORE       { $rnd_mode = $_[0]->round_mode($_[1]); }
53
54 BEGIN { tie $rnd_mode, 'Math::BigFloat'; }
55  
56 ##############################################################################
57
58 # in case we call SUPER::->foo() and this wants to call modify()
59 # sub modify () { 0; }
60
61 {
62   # valid method aliases for AUTOLOAD
63   my %methods = map { $_ => 1 }  
64    qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
65         fint facmp fcmp fzero fnan finf finc fdec flog
66         fceil ffloor frsft flsft fone flog
67       /;
68   # valid method's that can be hand-ed up (for AUTOLOAD)
69   my %hand_ups = map { $_ => 1 }  
70    qw / is_nan is_inf is_negative is_positive
71         accuracy precision div_scale round_mode fneg fabs babs fnot
72       /;
73
74   sub method_alias { return exists $methods{$_[0]||''}; } 
75   sub method_hand_up { return exists $hand_ups{$_[0]||''}; } 
76 }
77
78 ##############################################################################
79 # constructors
80
81 sub new 
82   {
83   # create a new BigFloat object from a string or another bigfloat object. 
84   # _e: exponent
85   # _m: mantissa
86   # sign  => sign (+/-), or "NaN"
87
88   my ($class,$wanted,@r) = @_;
89  
90   # avoid numify-calls by not using || on $wanted!
91   return $class->bzero() if !defined $wanted;   # default to 0
92   return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
93
94   my $self = {}; bless $self, $class;
95   # shortcut for bigints and its subclasses
96   if ((ref($wanted)) && (ref($wanted) ne $class))
97     {
98     $self->{_m} = $wanted->as_number();         # get us a bigint copy
99     $self->{_e} = Math::BigInt->bzero();
100     $self->{_m}->babs();
101     $self->{sign} = $wanted->sign();
102     return $self->bnorm();
103     }
104   # got string
105   # handle '+inf', '-inf' first
106   if ($wanted =~ /^[+-]?inf$/)
107     {
108     $self->{_e} = Math::BigInt->bzero();
109     $self->{_m} = Math::BigInt->bzero();
110     $self->{sign} = $wanted;
111     $self->{sign} = '+inf' if $self->{sign} eq 'inf';
112     return $self->bnorm();
113     }
114   #print "new string '$wanted'\n";
115   my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split(\$wanted);
116   if (!ref $mis)
117     {
118     die "$wanted is not a number initialized to $class" if !$NaNOK;
119     $self->{_e} = Math::BigInt->bzero();
120     $self->{_m} = Math::BigInt->bzero();
121     $self->{sign} = $nan;
122     }
123   else
124     {
125     # make integer from mantissa by adjusting exp, then convert to bigint
126     # undef,undef to signal MBI that we don't need no bloody rounding
127     $self->{_e} = Math::BigInt->new("$$es$$ev",undef,undef);    # exponent
128     $self->{_m} = Math::BigInt->new("$$miv$$mfv",undef,undef);  # create mant.
129     # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
130     $self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0;             
131     $self->{sign} = $$mis;
132     }
133   # print "mbf new ",join(' ',@r),"\n";
134   $self->bnorm()->round(@r);            # first normalize, then round
135   }
136
137 sub bnan
138   {
139   # create a bigfloat 'NaN', if given a BigFloat, set it to 'NaN'
140   my $self = shift;
141   $self = $class if !defined $self;
142   if (!ref($self))
143     {
144     my $c = $self; $self = {}; bless $self, $c;
145     }
146   $self->{_m} = Math::BigInt->bzero();
147   $self->{_e} = Math::BigInt->bzero();
148   $self->{sign} = $nan;
149   $self->{_a} = undef; $self->{_p} = undef;
150   $self;
151   }
152
153 sub binf
154   {
155   # create a bigfloat '+-inf', if given a BigFloat, set it to '+-inf'
156   my $self = shift;
157   my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
158
159   $self = $class if !defined $self;
160   if (!ref($self))
161     {
162     my $c = $self; $self = {}; bless $self, $c;
163     }
164   $self->{_m} = Math::BigInt->bzero();
165   $self->{_e} = Math::BigInt->bzero();
166   $self->{sign} = $sign.'inf';
167   $self->{_a} = undef; $self->{_p} = undef;
168   $self;
169   }
170
171 sub bone
172   {
173   # create a bigfloat '+-1', if given a BigFloat, set it to '+-1'
174   my $self = shift;
175   my $sign = shift; $sign = '+' if !defined $sign || $sign ne '-';
176
177   $self = $class if !defined $self;
178   if (!ref($self))
179     {
180     my $c = $self; $self = {}; bless $self, $c;
181     }
182   $self->{_m} = Math::BigInt->bone();
183   $self->{_e} = Math::BigInt->bzero();
184   $self->{sign} = $sign;
185   if (@_ > 0)
186     {
187     $self->{_a} = $_[0]
188      if (defined $self->{_a} && defined $_[0] && $_[0] > $self->{_a});
189     $self->{_p} = $_[1]
190      if (defined $self->{_p} && defined $_[1] && $_[1] < $self->{_p});
191     }
192   return $self;
193   }
194
195 sub bzero
196   {
197   # create a bigfloat '+0', if given a BigFloat, set it to 0
198   my $self = shift;
199   $self = $class if !defined $self;
200   if (!ref($self))
201     {
202     my $c = $self; $self = {}; bless $self, $c;
203     }
204   $self->{_m} = Math::BigInt->bzero();
205   $self->{_e} = Math::BigInt->bone();
206   $self->{sign} = '+';
207   if (@_ > 0)
208     {
209     $self->{_a} = $_[0]
210      if (defined $self->{_a} && defined $_[0] && $_[0] > $self->{_a});
211     $self->{_p} = $_[1]
212      if (defined $self->{_p} && defined $_[1] && $_[1] < $self->{_p});
213     }
214   return $self;
215   }
216
217 ##############################################################################
218 # string conversation
219
220 sub bstr 
221   {
222   # (ref to BFLOAT or num_str ) return num_str
223   # Convert number from internal format to (non-scientific) string format.
224   # internal format is always normalized (no leading zeros, "-0" => "+0")
225   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
226   #my $x = shift; my $class = ref($x) || $x;
227   #$x = $class->new(shift) unless ref($x);
228
229   #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
230   #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
231   if ($x->{sign} !~ /^[+-]$/)
232     {
233     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
234     return 'inf';                                       # +inf
235     }
236  
237   my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
238
239   my $not_zero = !$x->is_zero();
240   if ($not_zero)
241     {
242     $es = $x->{_m}->bstr();
243     $len = CORE::length($es);
244     if (!$x->{_e}->is_zero())
245       {
246       if ($x->{_e}->sign() eq '-')
247         {
248         $dot = '';
249         if ($x->{_e} <= -$len)
250           {
251           # print "style: 0.xxxx\n";
252           my $r = $x->{_e}->copy(); $r->babs()->bsub( CORE::length($es) );
253           $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
254           }
255         else
256           {
257           # print "insert '.' at $x->{_e} in '$es'\n";
258           substr($es,$x->{_e},0) = '.'; $cad = $x->{_e};
259           }
260         }
261       else
262         {
263         # expand with zeros
264         $es .= '0' x $x->{_e}; $len += $x->{_e}; $cad = 0;
265         }
266       }
267     } # if not zero
268   $es = $x->{sign}.$es if $x->{sign} eq '-';
269   # if set accuracy or precision, pad with zeros
270   if ((defined $x->{_a}) && ($not_zero))
271     {
272     # 123400 => 6, 0.1234 => 4, 0.001234 => 4
273     my $zeros = $x->{_a} - $cad;                # cad == 0 => 12340
274     $zeros = $x->{_a} - $len if $cad != $len;
275     $es .= $dot.'0' x $zeros if $zeros > 0;
276     }
277   elsif ($x->{_p} || 0 < 0)
278     {
279     # 123400 => 6, 0.1234 => 4, 0.001234 => 6
280     my $zeros = -$x->{_p} + $cad;
281     $es .= $dot.'0' x $zeros if $zeros > 0;
282     }
283   return $es;
284   }
285
286 sub bsstr
287   {
288   # (ref to BFLOAT or num_str ) return num_str
289   # Convert number from internal format to scientific string format.
290   # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
291   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
292   #my $x = shift; my $class = ref($x) || $x;
293   #$x = $class->new(shift) unless ref($x);
294
295   #die "Oups! e was $nan" if $x->{_e}->{sign} eq $nan;
296   #die "Oups! m was $nan" if $x->{_m}->{sign} eq $nan;
297   if ($x->{sign} !~ /^[+-]$/)
298     {
299     return $x->{sign} unless $x->{sign} eq '+inf';      # -inf, NaN
300     return 'inf';                                       # +inf
301     }
302   my $sign = $x->{_e}->{sign}; $sign = '' if $sign eq '-';
303   my $sep = 'e'.$sign;
304   return $x->{_m}->bstr().$sep.$x->{_e}->bstr();
305   }
306     
307 sub numify 
308   {
309   # Make a number from a BigFloat object
310   # simple return string and let Perl's atoi()/atof() handle the rest
311   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
312   return $x->bsstr(); 
313   }
314
315 ##############################################################################
316 # public stuff (usually prefixed with "b")
317
318 # tels 2001-08-04 
319 # todo: this must be overwritten and return NaN for non-integer values
320 # band(), bior(), bxor(), too
321 #sub bnot
322 #  {
323 #  $class->SUPER::bnot($class,@_);
324 #  }
325
326 sub bcmp 
327   {
328   # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
329   # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
330   my ($self,$x,$y) = objectify(2,@_);
331
332   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
333     {
334     # handle +-inf and NaN
335     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
336     return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
337     return +1 if $x->{sign} eq '+inf';
338     return -1 if $x->{sign} eq '-inf';
339     return -1 if $y->{sign} eq '+inf';
340     return +1 if $y->{sign} eq '-inf';
341     }
342
343   # check sign for speed first
344   return 1 if $x->{sign} eq '+' && $y->{sign} eq '-';   # does also 0 <=> -y
345   return -1 if $x->{sign} eq '-' && $y->{sign} eq '+';  # does also -x <=> 0
346
347   # shortcut 
348   my $xz = $x->is_zero();
349   my $yz = $y->is_zero();
350   return 0 if $xz && $yz;                               # 0 <=> 0
351   return -1 if $xz && $y->{sign} eq '+';                # 0 <=> +y
352   return 1 if $yz && $x->{sign} eq '+';                 # +x <=> 0
353
354   # adjust so that exponents are equal
355   my $lxm = $x->{_m}->length();
356   my $lym = $y->{_m}->length();
357   my $lx = $lxm + $x->{_e};
358   my $ly = $lym + $y->{_e};
359   # print "x $x y $y lx $lx ly $ly\n";
360   my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
361   # print "$l $x->{sign}\n";
362   return $l <=> 0 if $l != 0;
363   
364   # lengths (corrected by exponent) are equal
365   # so make mantissa euqal length by padding with zero (shift left)
366   my $diff = $lxm - $lym;
367   my $xm = $x->{_m};            # not yet copy it
368   my $ym = $y->{_m};
369   if ($diff > 0)
370     {
371     $ym = $y->{_m}->copy()->blsft($diff,10);
372     }
373   elsif ($diff < 0)
374     {
375     $xm = $x->{_m}->copy()->blsft(-$diff,10);
376     }
377   my $rc = $xm->bcmp($ym);
378   $rc = -$rc if $x->{sign} eq '-';              # -124 < -123
379   return $rc <=> 0;
380   }
381
382 sub bacmp 
383   {
384   # Compares 2 values, ignoring their signs. 
385   # Returns one of undef, <0, =0, >0. (suitable for sort)
386   # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
387   my ($self,$x,$y) = objectify(2,@_);
388
389   # handle +-inf and NaN's
390   if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
391     {
392     return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
393     return 0 if ($x->is_inf() && $y->is_inf());
394     return 1 if ($x->is_inf() && !$y->is_inf());
395     return -1 if (!$x->is_inf() && $y->is_inf());
396     }
397
398   # shortcut 
399   my $xz = $x->is_zero();
400   my $yz = $y->is_zero();
401   return 0 if $xz && $yz;                               # 0 <=> 0
402   return -1 if $xz && !$yz;                             # 0 <=> +y
403   return 1 if $yz && !$xz;                              # +x <=> 0
404
405   # adjust so that exponents are equal
406   my $lxm = $x->{_m}->length();
407   my $lym = $y->{_m}->length();
408   my $lx = $lxm + $x->{_e};
409   my $ly = $lym + $y->{_e};
410   # print "x $x y $y lx $lx ly $ly\n";
411   my $l = $lx - $ly;
412   # print "$l $x->{sign}\n";
413   return $l <=> 0 if $l != 0;
414   
415   # lengths (corrected by exponent) are equal
416   # so make mantissa equal-length by padding with zero (shift left)
417   my $diff = $lxm - $lym;
418   my $xm = $x->{_m};            # not yet copy it
419   my $ym = $y->{_m};
420   if ($diff > 0)
421     {
422     $ym = $y->{_m}->copy()->blsft($diff,10);
423     }
424   elsif ($diff < 0)
425     {
426     $xm = $x->{_m}->copy()->blsft(-$diff,10);
427     }
428   my $rc = $xm->bcmp($ym);
429   return $rc <=> 0;
430   }
431
432 sub badd 
433   {
434   # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
435   # return result as BFLOAT
436   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
437
438   # inf and NaN handling
439   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
440     {
441     # NaN first
442     return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
443     # inf handline
444     if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
445       {
446       # + and + => +, - and - => -, + and - => 0, - and + => 0
447       return $x->bzero() if $x->{sign} ne $y->{sign};
448       return $x;
449       }
450     # +-inf + something => +inf
451     # something +-inf => +-inf
452     $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
453     return $x;
454     }
455
456   # speed: no add for 0+y or x+0
457   return $x if $y->is_zero();                           # x+0
458   if ($x->is_zero())                                    # 0+y
459     {
460     # make copy, clobbering up x (modify in place!)
461     $x->{_e} = $y->{_e}->copy();
462     $x->{_m} = $y->{_m}->copy();
463     $x->{sign} = $y->{sign} || $nan;
464     return $x->round($a,$p,$r,$y);
465     }
466  
467   # take lower of the two e's and adapt m1 to it to match m2
468   my $e = $y->{_e}; $e = Math::BigInt::bzero() if !defined $e;  # if no BFLOAT
469   $e = $e - $x->{_e};
470   my $add = $y->{_m}->copy();
471   if ($e < 0)
472     {
473     my $e1 = $e->copy()->babs();
474     $x->{_m} *= (10 ** $e1);
475     $x->{_e} += $e;                     # need the sign of e
476     }
477   elsif ($e > 0)
478     {
479     $add *= (10 ** $e);
480     }
481   # else: both e are the same, so just leave them
482   $x->{_m}->{sign} = $x->{sign};                # fiddle with signs
483   $add->{sign} = $y->{sign};
484   $x->{_m} += $add;                             # finally do add/sub
485   $x->{sign} = $x->{_m}->{sign};                # re-adjust signs
486   $x->{_m}->{sign} = '+';                       # mantissa always positiv
487   # delete trailing zeros, then round
488   return $x->bnorm()->round($a,$p,$r,$y);
489   }
490
491 sub bsub 
492   {
493   # (BigFloat or num_str, BigFloat or num_str) return BigFloat
494   # subtract second arg from first, modify first
495   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
496
497   if (!$y->is_zero())           # don't need to do anything if $y is 0
498     {
499     $y->{sign} =~ tr/+\-/-+/;   # does nothing for NaN
500     $x->badd($y,$a,$p,$r);      # badd does not leave internal zeros
501     $y->{sign} =~ tr/+\-/-+/;   # refix $y (does nothing for NaN)
502     }
503   $x;                           # already rounded by badd()
504   }
505
506 sub binc
507   {
508   # increment arg by one
509   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
510
511   if ($x->{_e}->sign() eq '-')
512     {
513     return $x->badd($self->bone(),$a,$p,$r);    #  digits after dot
514     }
515
516   if (!$x->{_e}->is_zero())
517     {
518     $x->{_m}->blsft($x->{_e},10);               # 1e2 => 100
519     $x->{_e}->bzero();
520     }
521   # now $x->{_e} == 0
522   if ($x->{sign} eq '+')
523     {
524     $x->{_m}->binc();
525     return $x->bnorm()->bround($a,$p,$r);
526     }
527   elsif ($x->{sign} eq '-')
528     {
529     $x->{_m}->bdec();
530     $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
531     return $x->bnorm()->bround($a,$p,$r);
532     }
533   # inf, nan handling etc
534   $x->badd($self->__one(),$a,$p,$r);            # does round 
535   }
536
537 sub bdec
538   {
539   # decrement arg by one
540   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
541
542   if ($x->{_e}->sign() eq '-')
543     {
544     return $x->badd($self->bone('-'),$a,$p,$r); #  digits after dot
545     }
546
547   if (!$x->{_e}->is_zero())
548     {
549     $x->{_m}->blsft($x->{_e},10);               # 1e2 => 100
550     $x->{_e}->bzero();
551     }
552   # now $x->{_e} == 0
553   my $zero = $x->is_zero();
554   # <= 0
555   if (($x->{sign} eq '-') || $zero)
556     {
557     $x->{_m}->binc();
558     $x->{sign} = '-' if $zero;                  # 0 => 1 => -1
559     $x->{sign} = '+' if $x->{_m}->is_zero();    # -1 +1 => -0 => +0
560     return $x->bnorm()->round($a,$p,$r);
561     }
562   # > 0
563   elsif ($x->{sign} eq '+')
564     {
565     $x->{_m}->bdec();
566     return $x->bnorm()->round($a,$p,$r);
567     }
568   # inf, nan handling etc
569   $x->badd($self->bone('-'),$a,$p,$r);          # does round 
570   } 
571
572 sub blog
573   {
574   my ($self,$x,$base,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(2,@_);
575
576   # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
577
578   # u = x-1, v = x +1
579   #              _                               _
580   # taylor:     |    u    1   u^3   1   u^5       |
581   # ln (x)  = 2 |   --- + - * --- + - * --- + ... |  x > 0
582   #             |_   v    3    v    5   v        _|
583
584   return $x->bzero(@r) if $x->is_one();
585   return $x->bone(@r) if $x->bcmp($base) == 0;
586
587   my $d = $r[0] || $self->accuracy() || 40;
588   $d += 2;                                      # 2 more for rounding
589  
590   my $u = $x->copy(); $u->bdec();
591   my $v = $x->copy(); $v->binc();
592
593   $x->bdec()->bdiv($v,$d);                      # first term: u/v
594
595   $u *= $u; $v *= $v;
596   my $below = $v->copy()->bmul($v);
597   my $over = $u->copy()->bmul($u);
598   my $factor = $self->new(3); my $two = $self->new(2);
599
600   my $diff = $self->bone();
601   my $limit = $self->new("1E-". ($d-1)); my $last;
602   # print "diff $diff limit $limit\n";
603   while ($diff > $limit)
604     {
605     print "$x $over $below $factor\n";
606     $diff = $x->copy()->bsub($last)->babs();
607     print "diff $diff $limit\n";
608     $last = $x->copy();
609     $x += $over->copy()->bdiv($below->copy()->bmul($factor),$d);
610     $over *= $u; $below *= $v; $factor->badd($two);
611     }
612   $x->bmul($two);
613   return $x->round(@r);
614   }
615
616 sub blcm 
617   { 
618   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
619   # does not modify arguments, but returns new object
620   # Lowest Common Multiplicator
621
622   my ($self,@arg) = objectify(0,@_);
623   my $x = $self->new(shift @arg);
624   while (@arg) { $x = _lcm($x,shift @arg); } 
625   $x;
626   }
627
628 sub bgcd 
629   { 
630   # (BFLOAT or num_str, BFLOAT or num_str) return BINT
631   # does not modify arguments, but returns new object
632   # GCD -- Euclids algorithm Knuth Vol 2 pg 296
633    
634   my ($self,@arg) = objectify(0,@_);
635   my $x = $self->new(shift @arg);
636   while (@arg) { $x = _gcd($x,shift @arg); } 
637   $x;
638   }
639
640 sub is_zero
641   {
642   # return true if arg (BFLOAT or num_str) is zero (array '+', '0')
643   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
644
645   return 1 if $x->{sign} eq '+' && $x->{_m}->is_zero();
646   return 0;
647   }
648
649 sub is_one
650   {
651   # return true if arg (BFLOAT or num_str) is +1 (array '+', '1')
652   # or -1 if signis given
653   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
654
655   my $sign = shift || ''; $sign = '+' if $sign ne '-';
656   return 1
657    if ($x->{sign} eq $sign && $x->{_e}->is_zero() && $x->{_m}->is_one()); 
658   return 0;
659   }
660
661 sub is_odd
662   {
663   # return true if arg (BFLOAT or num_str) is odd or false if even
664   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
665   
666   return 0 if $x->{sign} !~ /^[+-]$/;                   # NaN & +-inf aren't
667   return 1 if ($x->{_e}->is_zero() && $x->{_m}->is_odd()); 
668   return 0;
669   }
670
671 sub is_even
672   {
673   # return true if arg (BINT or num_str) is even or false if odd
674   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
675
676   return 0 if $x->{sign} !~ /^[+-]$/;                   # NaN & +-inf aren't
677   return 1 if $x->{_m}->is_zero();                      # 0e1 is even
678   return 1 if ($x->{_e}->is_zero() && $x->{_m}->is_even()); # 123.45 is never
679   return 0;
680   }
681
682 sub bmul 
683   { 
684   # multiply two numbers -- stolen from Knuth Vol 2 pg 233
685   # (BINT or num_str, BINT or num_str) return BINT
686   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
687
688   # print "mbf bmul $x->{_m}e$x->{_e} $y->{_m}e$y->{_e}\n";
689   return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
690
691   # handle result = 0
692   return $x->bzero() if $x->is_zero() || $y->is_zero();
693   # inf handling
694   if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
695     {
696     # result will always be +-inf:
697     # +inf * +/+inf => +inf, -inf * -/-inf => +inf
698     # +inf * -/-inf => -inf, -inf * +/+inf => -inf
699     return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
700     return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
701     return $x->binf('-');
702     }
703
704   # aEb * cEd = (a*c)E(b+d)
705   $x->{_m}->bmul($y->{_m});
706   $x->{_e}->badd($y->{_e});
707   # adjust sign:
708   $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
709   return $x->bnorm()->round($a,$p,$r,$y);
710   }
711
712 sub bdiv 
713   {
714   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return 
715   # (BFLOAT,BFLOAT) (quo,rem) or BINT (only rem)
716   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
717
718   # x / +-inf => 0, reminder x
719   return wantarray ? ($x->bzero(),$x->copy()) : $x->bzero()
720    if $y->{sign} =~ /^[+-]inf$/;
721
722   # NaN if x == NaN or y == NaN or x==y==0
723   return wantarray ? ($x->bnan(),bnan()) : $x->bnan()
724    if (($x->is_nan() || $y->is_nan()) ||
725       ($x->is_zero() && $y->is_zero()));
726
727   # 5 / 0 => +inf, -6 / 0 => -inf
728   return wantarray
729    ? ($x->binf($x->{sign}),$self->bnan()) : $x->binf($x->{sign})
730    if ($x->{sign} =~ /^[+-]$/ && $y->is_zero());
731
732   # x== 0 or y == 1 or y == -1
733   return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
734
735   # we need to limit the accuracy to protect against overflow
736   my $fallback = 0;
737   my $scale = 0;
738   my @params = $x->_find_round_parameters($a,$p,$r,$y);
739
740   # no rounding at all, so must use fallback
741   if (scalar @params == 1)
742     {
743     # simulate old behaviour
744     $params[1] = $self->div_scale();    # and round to it as accuracy
745     $scale = $params[1]+4;              # at least four more for proper round
746     $params[3] = $r;                    # round mode by caller or undef
747     $fallback = 1;                      # to clear a/p afterwards
748     }
749   else
750     {
751     # the 4 below is empirical, and there might be cases where it is not
752     # enough...
753     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
754     }
755   my $lx = $x->{_m}->length(); my $ly = $y->{_m}->length();
756   $scale = $lx if $lx > $scale;
757   $scale = $ly if $ly > $scale;
758   my $diff = $ly - $lx;
759   $scale += $diff if $diff > 0;         # if lx << ly, but not if ly << lx!
760
761   $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+'; 
762
763   # check for / +-1 ( +/- 1E0)
764   if (!$y->is_one())
765     {
766     # promote BigInts and it's subclasses (except when already a BigFloat)
767     $y = $self->new($y) unless $y->isa('Math::BigFloat'); 
768
769     # calculate the result to $scale digits and then round it
770     # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
771     $x->{_m}->blsft($scale,10);
772     $x->{_m}->bdiv( $y->{_m} ); # a/c
773     $x->{_e}->bsub( $y->{_e} ); # b-d
774     $x->{_e}->bsub($scale);     # correct for 10**scale
775     $x->bnorm();                # remove trailing 0's
776     }
777
778   # shortcut to not run trough _find_round_parameters again
779   if (defined $params[1])
780     {
781     $x->bround($params[1],$params[3]);          # then round accordingly
782     }
783   else
784     {
785     $x->bfround($params[2],$params[3]);         # then round accordingly
786     }
787   if ($fallback)
788     {
789     # clear a/p after round, since user did not request it
790     $x->{_a} = undef; $x->{_p} = undef;
791     }
792   
793   if (wantarray)
794     {
795     my $rem;
796     if (!$y->is_one())
797       {
798       $rem = $x->copy();
799       $rem->bmod($y,$params[1],$params[2],$params[3]);
800       }
801     else
802       {
803       $rem = $self->bzero();
804       }
805     if ($fallback)
806       {
807       # clear a/p after round, since user did not request it
808       $rem->{_a} = undef; $rem->{_p} = undef;
809       }
810     return ($x,$rem);
811     }
812   return $x;
813   }
814
815 sub bmod 
816   {
817   # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return reminder 
818   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
819
820   if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
821     {
822     my ($d,$re) = $self->SUPER::_div_inf($x,$y);
823     return $re->round($a,$p,$r,$y);
824     } 
825   return $x->bnan() if $x->is_zero() && $y->is_zero();
826   return $x if $y->is_zero();
827   return $x->bnan() if $x->is_nan() || $y->is_nan();
828   return $x->bzero() if $y->is_one() || $x->is_zero();
829
830   # inf handling is missing here
831  
832   my $cmp = $x->bacmp($y);                      # equal or $x < $y?
833   return $x->bzero($a,$p) if $cmp == 0;         # $x == $y => result 0
834
835   # only $y of the operands negative? 
836   my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
837
838   $x->{sign} = $y->{sign};                              # calc sign first
839   return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0;  # $x < $y => result $x
840   
841   my $ym = $y->{_m}->copy();
842   
843   # 2e1 => 20
844   $ym->blsft($y->{_e},10) if $y->{_e}->{sign} eq '+' && !$y->{_e}->is_zero();
845  
846   # if $y has digits after dot
847   my $shifty = 0;                       # correct _e of $x by this
848   if ($y->{_e}->{sign} eq '-')          # has digits after dot
849     {
850     # 123 % 2.5 => 1230 % 25 => 5 => 0.5
851     $shifty = $y->{_e}->copy()->babs(); # no more digits after dot
852     $x->blsft($shifty,10);              # 123 => 1230, $y->{_m} is already 25
853     }
854   # $ym is now mantissa of $y based on exponent 0
855   
856   my $shiftx = 0;                       # correct _e of $x by this
857   if ($x->{_e}->{sign} eq '-')          # has digits after dot
858     {
859     # 123.4 % 20 => 1234 % 200
860     $shiftx = $x->{_e}->copy()->babs(); # no more digits after dot
861     $ym->blsft($shiftx,10);
862     }
863   # 123e1 % 20 => 1230 % 20
864   if ($x->{_e}->{sign} eq '+' && !$x->{_e}->is_zero())
865     {
866     $x->{_m}->blsft($x->{_e},10);
867     }
868   $x->{_e} = Math::BigInt->bzero() unless $x->{_e}->is_zero();
869   
870   $x->{_e}->bsub($shiftx) if $shiftx != 0;
871   $x->{_e}->bsub($shifty) if $shifty != 0;
872   
873   # now mantissas are equalized, exponent of $x is adjusted, so calc result
874   $x->{_m}->bmod($ym);
875
876   $x->{sign} = '+' if $x->{_m}->is_zero();              # fix sign for -0
877   $x->bnorm();
878
879   if ($neg != 0)        # one of them negative => correct in place
880     {
881     my $r = $y - $x;
882     $x->{_m} = $r->{_m};
883     $x->{_e} = $r->{_e};
884     $x->{sign} = '+' if $x->{_m}->is_zero();            # fix sign for -0
885     $x->bnorm();
886     }
887
888   $x->round($a,$p,$r,$y);       # round and return
889   }
890
891 sub bsqrt
892   { 
893   # calculate square root; this should probably
894   # use a different test to see whether the accuracy we want is...
895   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
896
897   return $x->bnan() if $x->{sign} eq 'NaN' || $x->{sign} =~ /^-/; # <0, NaN
898   return $x if $x->{sign} eq '+inf';                              # +inf
899   return $x if $x->is_zero() || $x->is_one();
900
901   # we need to limit the accuracy to protect against overflow
902   my $fallback = 0;
903   my $scale = 0;
904   my @params = $x->_find_round_parameters($a,$p,$r);
905
906   # no rounding at all, so must use fallback
907   if (scalar @params == 1)
908     {
909     # simulate old behaviour
910     $params[1] = $self->div_scale();    # and round to it as accuracy
911     $scale = $params[1]+4;              # at least four more for proper round
912     $params[3] = $r;                    # round mode by caller or undef
913     $fallback = 1;                      # to clear a/p afterwards
914     }
915   else
916     {
917     # the 4 below is empirical, and there might be cases where it is not
918     # enough...
919     $scale = abs($params[1] || $params[2]) + 4; # take whatever is defined
920     }
921
922   # when user set globals, they would interfere with our calculation, so
923   # disable then and later re-enable them
924   no strict 'refs';
925   my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
926   $abr = "$self\::precision"; my $pb = $$abr; $$abr = undef;
927   # we also need to disable any set A or P on $x (_find_round_parameters took
928   # them already into account), since these would interfere, too
929   delete $x->{_a}; delete $x->{_p};
930
931   my $xas = $x->as_number();
932   my $gs = $xas->copy()->bsqrt();       # some guess
933   if (($x->{_e}->{sign} ne '-')         # guess can't be accurate if there are
934                                         # digits after the dot
935    && ($xas->bcmp($gs * $gs) == 0))     # guess hit the nail on the head?
936     {
937     # exact result
938     $x->{_m} = $gs; $x->{_e} = Math::BigInt->bzero(); $x->bnorm();
939     # shortcut to not run trough _find_round_parameters again
940     if (defined $params[1])
941       {
942       $x->bround($params[1],$params[3]);        # then round accordingly
943       }
944     else
945       {
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     return $x;
954     }
955   $gs = $self->new( $gs );              # BigInt to BigFloat
956
957   my $lx = $x->{_m}->length();
958   $scale = $lx if $scale < $lx;
959   my $e = $self->new("1E-$scale");      # make test variable
960   return $x->bnan() if $e->sign() eq 'NaN';
961
962   my $y = $x->copy();
963   my $two = $self->new(2);
964   my $diff = $e;
965   # promote BigInts and it's subclasses (except when already a BigFloat)
966   $y = $self->new($y) unless $y->isa('Math::BigFloat'); 
967
968   my $rem;
969 #  my $steps = 0;
970   while ($diff >= $e)
971     {
972 #    return $x->bnan() if $gs->is_zero();
973
974     $rem = $y->copy()->bdiv($gs,$scale)->badd($gs)->bdiv($two,$scale);
975     $diff = $rem->copy()->bsub($gs)->babs();
976     $gs = $rem->copy();
977 #    $steps++;
978     }
979 #  print "steps $steps\n";
980   # copy over to modify $x
981   $x->{_m} = $rem->{_m}; $x->{_e} = $rem->{_e};
982   
983   # shortcut to not run trough _find_round_parameters again
984   if (defined $params[1])
985     {
986     $x->bround($params[1],$params[3]);          # then round accordingly
987     }
988   else
989     {
990     $x->bfround($params[2],$params[3]);         # then round accordingly
991     }
992   if ($fallback)
993     {
994     # clear a/p after round, since user did not request it
995     $x->{_a} = undef; $x->{_p} = undef;
996     }
997   # restore globals
998   ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
999   $x;
1000   }
1001
1002 sub bpow 
1003   {
1004   # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1005   # compute power of two numbers, second arg is used as integer
1006   # modifies first argument
1007
1008   my ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1009
1010   return $x if $x->{sign} =~ /^[+-]inf$/;
1011   return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
1012   return $x->bone() if $y->is_zero();
1013   return $x         if $x->is_one() || $y->is_one();
1014   my $y1 = $y->as_number();             # make bigint (trunc)
1015   # if ($x == -1)
1016   if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
1017     {
1018     # if $x == -1 and odd/even y => +1/-1  because +-1 ^ (+-1) => +-1
1019     return $y1->is_odd() ? $x : $x->babs(1);
1020     }
1021   return $x if $x->is_zero() && $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1022   # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1023   return $x->binf() if $x->is_zero() && $y->{sign} eq '-';
1024
1025   # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1026   $y1->babs();
1027   $x->{_m}->bpow($y1);
1028   $x->{_e}->bmul($y1);
1029   $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1030   $x->bnorm();
1031   if ($y->{sign} eq '-')
1032     {
1033     # modify $x in place!
1034     my $z = $x->copy(); $x->bzero()->binc();
1035     return $x->bdiv($z,$a,$p,$r);       # round in one go (might ignore y's A!)
1036     }
1037   return $x->round($a,$p,$r,$y);
1038   }
1039
1040 ###############################################################################
1041 # rounding functions
1042
1043 sub bfround
1044   {
1045   # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1046   # $n == 0 means round to integer
1047   # expects and returns normalized numbers!
1048   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1049
1050   return $x if $x->modify('bfround');
1051   
1052   my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
1053   return $x if !defined $scale;                 # no-op
1054
1055   # never round a 0, +-inf, NaN
1056   if ($x->is_zero())
1057     {
1058     $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1059     return $x; 
1060     }
1061   return $x if $x->{sign} !~ /^[+-]$/;
1062   # print "MBF bfround $x to scale $scale mode $mode\n";
1063
1064   # don't round if x already has lower precision
1065   return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1066
1067   $x->{_p} = $scale;                    # remember round in any case
1068   $x->{_a} = undef;                     # and clear A
1069   if ($scale < 0)
1070     {
1071     # print "bfround scale $scale e $x->{_e}\n";
1072     # round right from the '.'
1073     return $x if $x->{_e} >= 0;                 # nothing to round
1074     $scale = -$scale;                           # positive for simplicity
1075     my $len = $x->{_m}->length();               # length of mantissa
1076     my $dad = -$x->{_e};                        # digits after dot
1077     my $zad = 0;                                # zeros after dot
1078     $zad = -$len-$x->{_e} if ($x->{_e} < -$len);# for 0.00..00xxx style
1079     #print "scale $scale dad $dad zad $zad len $len\n";
1080
1081     # number  bsstr   len zad dad       
1082     # 0.123   123e-3    3   0 3
1083     # 0.0123  123e-4    3   1 4
1084     # 0.001   1e-3      1   2 3
1085     # 1.23    123e-2    3   0 2
1086     # 1.2345  12345e-4  5   0 4
1087
1088     # do not round after/right of the $dad
1089     return $x if $scale > $dad;                 # 0.123, scale >= 3 => exit
1090
1091     # round to zero if rounding inside the $zad, but not for last zero like:
1092     # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1093     return $x->bzero() if $scale < $zad;
1094     if ($scale == $zad)                 # for 0.006, scale -3 and trunc
1095       {
1096       $scale = -$len-1;
1097       }
1098     else
1099       {
1100       # adjust round-point to be inside mantissa
1101       if ($zad != 0)
1102         {
1103         $scale = $scale-$zad;
1104         }
1105       else
1106         {
1107         my $dbd = $len - $dad; $dbd = 0 if $dbd < 0;    # digits before dot
1108         $scale = $dbd+$scale;
1109         }
1110       }
1111     # print "round to $x->{_m} to $scale\n";
1112     }
1113   else
1114     {
1115     # 123 => 100 means length(123) = 3 - $scale (2) => 1
1116
1117     # calculate digits before dot
1118     my $dbt = $x->{_m}->length(); $dbt += $x->{_e} if $x->{_e}->sign() eq '-';
1119     # if not enough digits before dot, round to zero
1120     return $x->bzero() if ($scale > $dbt) && ($dbt < 0);
1121     # scale always >= 0 here
1122     if ($dbt == 0)
1123       {
1124       # 0.49->bfround(1): scale == 1, dbt == 0: => 0.0
1125       # 0.51->bfround(0): scale == 0, dbt == 0: => 1.0
1126       # 0.5->bfround(0):  scale == 0, dbt == 0: => 0
1127       # 0.05->bfround(0): scale == 0, dbt == 0: => 0
1128       # print "$scale $dbt $x->{_m}\n";
1129       $scale = -$x->{_m}->length();
1130       }
1131     elsif ($dbt > 0)
1132       {
1133       # correct by subtracting scale
1134       $scale = $dbt - $scale;
1135       }
1136     else
1137       {
1138       $scale = $x->{_m}->length() - $scale;
1139       }
1140     }
1141   # print "using $scale for $x->{_m} with '$mode'\n";
1142   # pass sign to bround for rounding modes '+inf' and '-inf'
1143   $x->{_m}->{sign} = $x->{sign};
1144   $x->{_m}->bround($scale,$mode);
1145   $x->{_m}->{sign} = '+';               # fix sign back
1146   $x->bnorm();
1147   }
1148
1149 sub bround
1150   {
1151   # accuracy: preserve $N digits, and overwrite the rest with 0's
1152   my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1153   
1154   die ('bround() needs positive accuracy') if ($_[0] || 0) < 0;
1155
1156   my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1157   return $x if !defined $scale;                         # no-op
1158
1159   return $x if $x->modify('bround');
1160
1161   # scale is now either $x->{_a}, $accuracy, or the user parameter
1162   # test whether $x already has lower accuracy, do nothing in this case 
1163   # but do round if the accuracy is the same, since a math operation might
1164   # want to round a number with A=5 to 5 digits afterwards again
1165   return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
1166
1167   # scale < 0 makes no sense
1168   # never round a +-inf, NaN
1169   return $x if ($scale < 0) ||  $x->{sign} !~ /^[+-]$/;
1170
1171   # 1: $scale == 0 => keep all digits
1172   # 2: never round a 0
1173   # 3: if we should keep more digits than the mantissa has, do nothing
1174   if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1175     {
1176     $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1177     return $x; 
1178     }
1179
1180   # pass sign to bround for '+inf' and '-inf' rounding modes
1181   $x->{_m}->{sign} = $x->{sign};
1182   $x->{_m}->bround($scale,$mode);       # round mantissa
1183   $x->{_m}->{sign} = '+';               # fix sign back
1184   # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
1185   $x->{_a} = $scale;                    # remember rounding
1186   $x->{_p} = undef;                     # and clear P
1187   $x->bnorm();                          # del trailing zeros gen. by bround()
1188   }
1189
1190 sub bfloor
1191   {
1192   # return integer less or equal then $x
1193   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1194
1195   return $x if $x->modify('bfloor');
1196    
1197   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1198
1199   # if $x has digits after dot
1200   if ($x->{_e}->{sign} eq '-')
1201     {
1202     $x->{_m}->brsft(-$x->{_e},10);
1203     $x->{_e}->bzero();
1204     $x-- if $x->{sign} eq '-';
1205     }
1206   $x->round($a,$p,$r);
1207   }
1208
1209 sub bceil
1210   {
1211   # return integer greater or equal then $x
1212   my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1213
1214   return $x if $x->modify('bceil');
1215   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1216
1217   # if $x has digits after dot
1218   if ($x->{_e}->{sign} eq '-')
1219     {
1220     $x->{_m}->brsft(-$x->{_e},10);
1221     $x->{_e}->bzero();
1222     $x++ if $x->{sign} eq '+';
1223     }
1224   $x->round($a,$p,$r);
1225   }
1226
1227 sub brsft
1228   {
1229   # shift right by $y (divide by power of 2)
1230   my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1231
1232   return $x if $x->modify('brsft');
1233   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1234
1235   $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1236   $x->bdiv($n ** $y,$a,$p,$r,$y);
1237   }
1238
1239 sub blsft
1240   {
1241   # shift right by $y (divide by power of 2)
1242   my ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1243
1244   return $x if $x->modify('brsft');
1245   return $x if $x->{sign} !~ /^[+-]$/;  # nan, +inf, -inf
1246
1247   $n = 2 if !defined $n; $n = Math::BigFloat->new($n);
1248   $x->bmul($n ** $y,$a,$p,$r,$y);
1249   }
1250
1251 ###############################################################################
1252
1253 sub DESTROY
1254   {
1255   # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
1256   }
1257
1258 sub AUTOLOAD
1259   {
1260   # make fxxx and bxxx work
1261   # my $self = $_[0];
1262   my $name = $AUTOLOAD;
1263
1264   $name =~ s/.*:://;    # split package
1265   #print "$name\n";
1266   no strict 'refs';
1267   if (!method_alias($name))
1268     {
1269     if (!defined $name)
1270       {
1271       # delayed load of Carp and avoid recursion        
1272       require Carp;
1273       Carp::croak ("Can't call a method without name");
1274       }
1275     if (!method_hand_up($name))
1276       {
1277       # delayed load of Carp and avoid recursion        
1278       require Carp;
1279       Carp::croak ("Can't call $class\-\>$name, not a valid method");
1280       }
1281     # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1282     $name =~ s/^f/b/;
1283     return &{'Math::BigInt'."::$name"}(@_);
1284     }
1285   my $bname = $name; $bname =~ s/^f/b/;
1286   *{$class."\:\:$name"} = \&$bname;
1287   &$bname;      # uses @_
1288   }
1289
1290 sub exponent
1291   {
1292   # return a copy of the exponent
1293   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1294
1295   if ($x->{sign} !~ /^[+-]$/)
1296     {
1297     my $s = $x->{sign}; $s =~ s/^[+-]//;
1298     return $self->new($s);                      # -inf, +inf => +inf
1299     }
1300   return $x->{_e}->copy();
1301   }
1302
1303 sub mantissa
1304   {
1305   # return a copy of the mantissa
1306   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1307  
1308   if ($x->{sign} !~ /^[+-]$/)
1309     {
1310     my $s = $x->{sign}; $s =~ s/^[+]//;
1311     return $self->new($s);                      # -inf, +inf => +inf
1312     }
1313   my $m = $x->{_m}->copy();             # faster than going via bstr()
1314   $m->bneg() if $x->{sign} eq '-';
1315
1316   $m;
1317   }
1318
1319 sub parts
1320   {
1321   # return a copy of both the exponent and the mantissa
1322   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1323
1324   if ($x->{sign} !~ /^[+-]$/)
1325     {
1326     my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
1327     return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
1328     }
1329   my $m = $x->{_m}->copy();     # faster than going via bstr()
1330   $m->bneg() if $x->{sign} eq '-';
1331   return ($m,$x->{_e}->copy());
1332   }
1333
1334 ##############################################################################
1335 # private stuff (internal use only)
1336
1337 sub import
1338   {
1339   my $self = shift;
1340   for ( my $i = 0; $i < @_ ; $i++ )
1341     {
1342     if ( $_[$i] eq ':constant' )
1343       {
1344       # this rest causes overlord er load to step in
1345       # print "overload @_\n";
1346       overload::constant float => sub { $self->new(shift); }; 
1347       splice @_, $i, 1; last;
1348       }
1349     }
1350   # any non :constant stuff is handled by our parent, Exporter
1351   # even if @_ is empty, to give it a chance
1352   $self->SUPER::import(@_);             # for subclasses
1353   $self->export_to_level(1,$self,@_);   # need this, too
1354   }
1355
1356 sub bnorm
1357   {
1358   # adjust m and e so that m is smallest possible
1359   # round number according to accuracy and precision settings
1360   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1361
1362   return $x if $x->{sign} !~ /^[+-]$/;          # inf, nan etc
1363
1364   my $zeros = $x->{_m}->_trailing_zeros();      # correct for trailing zeros 
1365   if ($zeros != 0)
1366     {
1367     $x->{_m}->brsft($zeros,10); $x->{_e} += $zeros;
1368     }
1369   # for something like 0Ey, set y to 1, and -0 => +0
1370   $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
1371   # this is to prevent automatically rounding when MBI's globals are set
1372   $x->{_m}->{_f} = MB_NEVER_ROUND;
1373   $x->{_e}->{_f} = MB_NEVER_ROUND;
1374   # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
1375   $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
1376   $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
1377   $x;                                   # MBI bnorm is no-op, so dont call it
1378   } 
1379  
1380 ##############################################################################
1381 # internal calculation routines
1382
1383 sub as_number
1384   {
1385   # return copy as a bigint representation of this BigFloat number
1386   my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
1387
1388   my $z;
1389   if ($x->{_e}->is_zero())
1390     {
1391     $z = $x->{_m}->copy();
1392     $z->{sign} = $x->{sign};
1393     return $z;
1394     }
1395   $z = $x->{_m}->copy();
1396   if ($x->{_e} < 0)
1397     {
1398     $z->brsft(-$x->{_e},10);
1399     } 
1400   else
1401     {
1402     $z->blsft($x->{_e},10);
1403     }
1404   $z->{sign} = $x->{sign};
1405   $z;
1406   }
1407
1408 sub length
1409   {
1410   my $x = shift;
1411   my $class = ref($x) || $x;
1412   $x = $class->new(shift) unless ref($x);
1413
1414   return 1 if $x->{_m}->is_zero();
1415   my $len = $x->{_m}->length();
1416   $len += $x->{_e} if $x->{_e}->sign() eq '+';
1417   if (wantarray())
1418     {
1419     my $t = Math::BigInt::bzero();
1420     $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
1421     return ($len,$t);
1422     }
1423   $len;
1424   }
1425
1426 1;
1427 __END__
1428
1429 =head1 NAME
1430
1431 Math::BigFloat - Arbitrary size floating point math package
1432
1433 =head1 SYNOPSIS
1434
1435   use Math::BigFloat;
1436
1437   # Number creation     
1438   $x = Math::BigInt->new($str); # defaults to 0
1439   $nan  = Math::BigInt->bnan(); # create a NotANumber
1440   $zero = Math::BigInt->bzero();# create a "+0"
1441
1442   # Testing
1443   $x->is_zero();                # return whether arg is zero or not
1444   $x->is_nan();                 # return whether arg is NaN or not
1445   $x->is_one();                 # true if arg is +1
1446   $x->is_one('-');              # true if arg is -1
1447   $x->is_odd();                 # true if odd, false for even
1448   $x->is_even();                # true if even, false for odd
1449   $x->is_positive();            # true if >= 0
1450   $x->is_negative();            # true if <  0
1451   $x->is_inf(sign)              # true if +inf or -inf (sign default '+')
1452   $x->bcmp($y);                 # compare numbers (undef,<0,=0,>0)
1453   $x->bacmp($y);                # compare absolutely (undef,<0,=0,>0)
1454   $x->sign();                   # return the sign, either +,- or NaN
1455
1456   # The following all modify their first argument:
1457
1458   # set 
1459   $x->bzero();                  # set $i to 0
1460   $x->bnan();                   # set $i to NaN
1461
1462   $x->bneg();                   # negation
1463   $x->babs();                   # absolute value
1464   $x->bnorm();                  # normalize (no-op)
1465   $x->bnot();                   # two's complement (bit wise not)
1466   $x->binc();                   # increment x by 1
1467   $x->bdec();                   # decrement x by 1
1468   
1469   $x->badd($y);                 # addition (add $y to $x)
1470   $x->bsub($y);                 # subtraction (subtract $y from $x)
1471   $x->bmul($y);                 # multiplication (multiply $x by $y)
1472   $x->bdiv($y);                 # divide, set $i to quotient
1473                                 # return (quo,rem) or quo if scalar
1474
1475   $x->bmod($y);                 # modulus
1476   $x->bpow($y);                 # power of arguments (a**b)
1477   $x->blsft($y);                # left shift
1478   $x->brsft($y);                # right shift 
1479                                 # return (quo,rem) or quo if scalar
1480   
1481   $x->blog($base);              # logarithm of $x, base defaults to e
1482                                 # (other bases than e not supported yet)
1483   
1484   $x->band($y);                 # bit-wise and
1485   $x->bior($y);                 # bit-wise inclusive or
1486   $x->bxor($y);                 # bit-wise exclusive or
1487   $x->bnot();                   # bit-wise not (two's complement)
1488   
1489   $x->bround($N);               # accuracy: preserver $N digits
1490   $x->bfround($N);              # precision: round to the $Nth digit
1491
1492   # The following do not modify their arguments:
1493
1494   bgcd(@values);                # greatest common divisor
1495   blcm(@values);                # lowest common multiplicator
1496   
1497   $x->bstr();                   # return string
1498   $x->bsstr();                  # return string in scientific notation
1499   
1500   $x->exponent();               # return exponent as BigInt
1501   $x->mantissa();               # return mantissa as BigInt
1502   $x->parts();                  # return (mantissa,exponent) as BigInt
1503
1504   $x->length();                 # number of digits (w/o sign and '.')
1505   ($l,$f) = $x->length();       # number of digits, and length of fraction      
1506
1507 =head1 DESCRIPTION
1508
1509 All operators (inlcuding basic math operations) are overloaded if you
1510 declare your big floating point numbers as
1511
1512   $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
1513
1514 Operations with overloaded operators preserve the arguments, which is
1515 exactly what you expect.
1516
1517 =head2 Canonical notation
1518
1519 Input to these routines are either BigFloat objects, or strings of the
1520 following four forms:
1521
1522 =over 2
1523
1524 =item *
1525
1526 C</^[+-]\d+$/>
1527
1528 =item *
1529
1530 C</^[+-]\d+\.\d*$/>
1531
1532 =item *
1533
1534 C</^[+-]\d+E[+-]?\d+$/>
1535
1536 =item *
1537
1538 C</^[+-]\d*\.\d+E[+-]?\d+$/>
1539
1540 =back
1541
1542 all with optional leading and trailing zeros and/or spaces. Additonally,
1543 numbers are allowed to have an underscore between any two digits.
1544
1545 Empty strings as well as other illegal numbers results in 'NaN'.
1546
1547 bnorm() on a BigFloat object is now effectively a no-op, since the numbers 
1548 are always stored in normalized form. On a string, it creates a BigFloat 
1549 object.
1550
1551 =head2 Output
1552
1553 Output values are BigFloat objects (normalized), except for bstr() and bsstr().
1554
1555 The string output will always have leading and trailing zeros stripped and drop
1556 a plus sign. C<bstr()> will give you always the form with a decimal point,
1557 while C<bsstr()> (for scientific) gives you the scientific notation.
1558
1559         Input                   bstr()          bsstr()
1560         '-0'                    '0'             '0E1'
1561         '  -123 123 123'        '-123123123'    '-123123123E0'
1562         '00.0123'               '0.0123'        '123E-4'
1563         '123.45E-2'             '1.2345'        '12345E-4'
1564         '10E+3'                 '10000'         '1E4'
1565
1566 Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
1567 C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
1568 return either undef, <0, 0 or >0 and are suited for sort.
1569
1570 Actual math is done by using BigInts to represent the mantissa and exponent.
1571 The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to 
1572 represent the result when input arguments are not numbers, as well as 
1573 the result of dividing by zero.
1574
1575 =head2 C<mantissa()>, C<exponent()> and C<parts()>
1576
1577 C<mantissa()> and C<exponent()> return the said parts of the BigFloat 
1578 as BigInts such that:
1579
1580         $m = $x->mantissa();
1581         $e = $x->exponent();
1582         $y = $m * ( 10 ** $e );
1583         print "ok\n" if $x == $y;
1584
1585 C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
1586
1587 A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
1588
1589 Currently the mantissa is reduced as much as possible, favouring higher
1590 exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
1591 This might change in the future, so do not depend on it.
1592
1593 =head2 Accuracy vs. Precision
1594
1595 See also: L<Rounding|Rounding>.
1596
1597 Math::BigFloat supports both precision and accuracy. For a full documentation,
1598 examples and tips on these topics please see the large section in
1599 L<Math::BigInt>.
1600
1601 Since things like sqrt(2) or 1/3 must presented with a limited precision lest
1602 a operation consumes all resources, each operation produces no more than
1603 C<Math::BigFloat::precision()> digits.
1604
1605 In case the result of one operation has more precision than specified,
1606 it is rounded. The rounding mode taken is either the default mode, or the one
1607 supplied to the operation after the I<scale>:
1608
1609         $x = Math::BigFloat->new(2);
1610         Math::BigFloat::precision(5);           # 5 digits max
1611         $y = $x->copy()->bdiv(3);               # will give 0.66666
1612         $y = $x->copy()->bdiv(3,6);             # will give 0.666666
1613         $y = $x->copy()->bdiv(3,6,'odd');       # will give 0.666667
1614         Math::BigFloat::round_mode('zero');
1615         $y = $x->copy()->bdiv(3,6);             # will give 0.666666
1616
1617 =head2 Rounding
1618
1619 =over 2
1620
1621 =item ffround ( +$scale )
1622
1623 Rounds to the $scale'th place left from the '.', counting from the dot.
1624 The first digit is numbered 1. 
1625
1626 =item ffround ( -$scale )
1627
1628 Rounds to the $scale'th place right from the '.', counting from the dot.
1629
1630 =item ffround ( 0 )
1631
1632 Rounds to an integer.
1633
1634 =item fround  ( +$scale )
1635
1636 Preserves accuracy to $scale digits from the left (aka significant digits)
1637 and pads the rest with zeros. If the number is between 1 and -1, the
1638 significant digits count from the first non-zero after the '.'
1639
1640 =item fround  ( -$scale ) and fround ( 0 )
1641
1642 These are effetively no-ops.
1643
1644 =back
1645
1646 All rounding functions take as a second parameter a rounding mode from one of
1647 the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
1648
1649 The default rounding mode is 'even'. By using
1650 C<< Math::BigFloat::round_mode($round_mode); >> you can get and set the default
1651 mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
1652 no longer supported.
1653 The second parameter to the round functions then overrides the default
1654 temporarily. 
1655
1656 The C<< as_number() >> function returns a BigInt from a Math::BigFloat. It uses
1657 'trunc' as rounding mode to make it equivalent to:
1658
1659         $x = 2.5;
1660         $y = int($x) + 2;
1661
1662 You can override this by passing the desired rounding mode as parameter to
1663 C<as_number()>:
1664
1665         $x = Math::BigFloat->new(2.5);
1666         $y = $x->as_number('odd');      # $y = 3
1667
1668 =head1 EXAMPLES
1669  
1670   # not ready yet
1671
1672 =head1 Autocreating constants
1673
1674 After C<use Math::BigFloat ':constant'> all the floating point constants
1675 in the given scope are converted to C<Math::BigFloat>. This conversion
1676 happens at compile time.
1677
1678 In particular
1679
1680   perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
1681
1682 prints the value of C<2E-100>.  Note that without conversion of 
1683 constants the expression 2E-100 will be calculated as normal floating point 
1684 number.
1685
1686 =head1 BUGS
1687
1688 =over 2
1689
1690 =item *
1691
1692 The following does not work yet:
1693
1694         $m = $x->mantissa();
1695         $e = $x->exponent();
1696         $y = $m * ( 10 ** $e );
1697         print "ok\n" if $x == $y;
1698
1699 =item *
1700
1701 There is no fmod() function yet.
1702
1703 =back
1704
1705 =head1 CAVEAT
1706
1707 =over 1
1708
1709 =item stringify, bstr()
1710
1711 Both stringify and bstr() now drop the leading '+'. The old code would return
1712 '+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
1713 reasoning and details.
1714
1715 =item bdiv
1716
1717 The following will probably not do what you expect:
1718
1719         print $c->bdiv(123.456),"\n";
1720
1721 It prints both quotient and reminder since print works in list context. Also,
1722 bdiv() will modify $c, so be carefull. You probably want to use
1723         
1724         print $c / 123.456,"\n";
1725         print scalar $c->bdiv(123.456),"\n";  # or if you want to modify $c
1726
1727 instead.
1728
1729 =item Modifying and =
1730
1731 Beware of:
1732
1733         $x = Math::BigFloat->new(5);
1734         $y = $x;
1735
1736 It will not do what you think, e.g. making a copy of $x. Instead it just makes
1737 a second reference to the B<same> object and stores it in $y. Thus anything
1738 that modifies $x will modify $y, and vice versa.
1739
1740         $x->bmul(2);
1741         print "$x, $y\n";       # prints '10, 10'
1742
1743 If you want a true copy of $x, use:
1744         
1745         $y = $x->copy();
1746
1747 See also the documentation in L<overload> regarding C<=>.
1748
1749 =item bpow
1750
1751 C<bpow()> now modifies the first argument, unlike the old code which left
1752 it alone and only returned the result. This is to be consistent with
1753 C<badd()> etc. The first will modify $x, the second one won't:
1754
1755         print bpow($x,$i),"\n";         # modify $x
1756         print $x->bpow($i),"\n";        # ditto
1757         print $x ** $i,"\n";            # leave $x alone 
1758
1759 =back
1760
1761 =head1 LICENSE
1762
1763 This program is free software; you may redistribute it and/or modify it under
1764 the same terms as Perl itself.
1765
1766 =head1 AUTHORS
1767
1768 Mark Biggar, overloaded interface by Ilya Zakharevich.
1769 Completely rewritten by Tels http://bloodgate.com in 2001.
1770
1771 =cut