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