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