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