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