projects / p5sagit/p5-mst-13.2.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | history | HEAD
Math::BigInt v1.67 released
[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
091c87b1 15$VERSION = '1.41';
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.
8df1e0a2 175 # print $self->{_e}, " ", $self->{_m},"\n";
58cde26e 176 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
8df1e0a2 177 $self->{_e} -= CORE::length($$mfv) if CORE::length($$mfv) != 0;
027dc388 178 $self->{sign} = $$mis;
58cde26e 179 }
28df3e88 180 # if downgrade, inf, NaN or integers go down
181
182 if ($downgrade && $self->{_e}->{sign} eq '+')
183 {
990fb837 184 #print "downgrading $$miv$$mfv"."E$$es$$ev";
28df3e88 185 if ($self->{_e}->is_zero())
186 {
187 $self->{_m}->{sign} = $$mis; # negative if wanted
188 return $downgrade->new($self->{_m});
189 }
8df1e0a2 190 return $downgrade->new($self->bsstr());
28df3e88 191 }
990fb837 192 #print "mbf new $self->{sign} $self->{_m} e $self->{_e} ",ref($self),"\n";
193 $self->bnorm()->round(@r); # first normalize, then round
58cde26e 194 }
a0d0e21e 195
13a12e00 196sub _bnan
58cde26e 197 {
990fb837 198 # used by parent class bone() to initialize number to NaN
58cde26e 199 my $self = shift;
990fb837 200
201 if ($_trap_nan)
202 {
203 require Carp;
204 my $class = ref($self);
205 Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
206 }
207
208 $IMPORT=1; # call our import only once
56b9c951 209 $self->{_m} = $MBI->bzero();
210 $self->{_e} = $MBI->bzero();
58cde26e 211 }
a0d0e21e 212
13a12e00 213sub _binf
58cde26e 214 {
990fb837 215 # used by parent class bone() to initialize number to +-inf
58cde26e 216 my $self = shift;
990fb837 217
218 if ($_trap_inf)
219 {
220 require Carp;
221 my $class = ref($self);
222 Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
223 }
224
225 $IMPORT=1; # call our import only once
56b9c951 226 $self->{_m} = $MBI->bzero();
227 $self->{_e} = $MBI->bzero();
58cde26e 228 }
a0d0e21e 229
13a12e00 230sub _bone
574bacfe 231 {
13a12e00 232 # used by parent class bone() to initialize number to 1
574bacfe 233 my $self = shift;
990fb837 234 $IMPORT=1; # call our import only once
56b9c951 235 $self->{_m} = $MBI->bone();
236 $self->{_e} = $MBI->bzero();
574bacfe 237 }
238
13a12e00 239sub _bzero
58cde26e 240 {
990fb837 241 # used by parent class bone() to initialize number to 0
58cde26e 242 my $self = shift;
990fb837 243 $IMPORT=1; # call our import only once
56b9c951 244 $self->{_m} = $MBI->bzero();
245 $self->{_e} = $MBI->bone();
58cde26e 246 }
247
9393ace2 248sub isa
249 {
250 my ($self,$class) = @_;
56b9c951 251 return if $class =~ /^Math::BigInt/; # we aren't one of these
252 UNIVERSAL::isa($self,$class);
9393ace2 253 }
254
8f675a64 255sub config
256 {
257 # return (later set?) configuration data as hash ref
258 my $class = shift || 'Math::BigFloat';
259
990fb837 260 my $cfg = $class->SUPER::config(@_);
8f675a64 261
990fb837 262 # now we need only to override the ones that are different from our parent
8f675a64 263 $cfg->{class} = $class;
264 $cfg->{with} = $MBI;
8f675a64 265 $cfg;
266 }
267
58cde26e 268##############################################################################
269# string conversation
270
271sub bstr
272 {
273 # (ref to BFLOAT or num_str ) return num_str
274 # Convert number from internal format to (non-scientific) string format.
275 # internal format is always normalized (no leading zeros, "-0" => "+0")
ee15d750 276 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 277
574bacfe 278 if ($x->{sign} !~ /^[+-]$/)
58cde26e 279 {
574bacfe 280 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
281 return 'inf'; # +inf
58cde26e 282 }
c38b2de2 283
574bacfe 284 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
285
c38b2de2 286 # $x is zero?
287 my $not_zero = !($x->{sign} eq '+' && $x->{_m}->is_zero());
574bacfe 288 if ($not_zero)
58cde26e 289 {
574bacfe 290 $es = $x->{_m}->bstr();
291 $len = CORE::length($es);
c38b2de2 292 my $e = $x->{_e}->numify();
293 if ($e < 0)
58cde26e 294 {
c38b2de2 295 $dot = '';
296 # if _e is bigger than a scalar, the following will blow your memory
297 if ($e <= -$len)
574bacfe 298 {
c38b2de2 299 #print "style: 0.xxxx\n";
300 my $r = abs($e) - $len;
301 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
574bacfe 302 }
303 else
304 {
c38b2de2 305 #print "insert '.' at $e in '$es'\n";
306 substr($es,$e,0) = '.'; $cad = $x->{_e};
574bacfe 307 }
82cf049f 308 }
c38b2de2 309 elsif ($e > 0)
310 {
311 # expand with zeros
312 $es .= '0' x $e; $len += $e; $cad = 0;
313 }
574bacfe 314 } # if not zero
c38b2de2 315 $es = '-'.$es if $x->{sign} eq '-';
316 # if set accuracy or precision, pad with zeros on the right side
574bacfe 317 if ((defined $x->{_a}) && ($not_zero))
318 {
319 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
320 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
321 $zeros = $x->{_a} - $len if $cad != $len;
574bacfe 322 $es .= $dot.'0' x $zeros if $zeros > 0;
82cf049f 323 }
c38b2de2 324 elsif ((($x->{_p} || 0) < 0))
58cde26e 325 {
574bacfe 326 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
327 my $zeros = -$x->{_p} + $cad;
574bacfe 328 $es .= $dot.'0' x $zeros if $zeros > 0;
58cde26e 329 }
56b9c951 330 $es;
82cf049f 331 }
f216259d 332
58cde26e 333sub bsstr
334 {
335 # (ref to BFLOAT or num_str ) return num_str
336 # Convert number from internal format to scientific string format.
337 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
ee15d750 338 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
a0d0e21e 339
574bacfe 340 if ($x->{sign} !~ /^[+-]$/)
341 {
342 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
343 return 'inf'; # +inf
344 }
091c87b1 345 # do $esign, because we need '1e+1', since $x->{_e}->bstr() misses the +
56d9de68 346 my $esign = $x->{_e}->{sign}; $esign = '' if $esign eq '-';
347 my $sep = 'e'.$esign;
348 my $sign = $x->{sign}; $sign = '' if $sign eq '+';
349 $sign . $x->{_m}->bstr() . $sep . $x->{_e}->bstr();
58cde26e 350 }
351
352sub numify
353 {
354 # Make a number from a BigFloat object
574bacfe 355 # simple return string and let Perl's atoi()/atof() handle the rest
ee15d750 356 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
56b9c951 357 $x->bsstr();
58cde26e 358 }
a0d0e21e 359
58cde26e 360##############################################################################
361# public stuff (usually prefixed with "b")
362
574bacfe 363# tels 2001-08-04
364# todo: this must be overwritten and return NaN for non-integer values
365# band(), bior(), bxor(), too
58cde26e 366#sub bnot
367# {
368# $class->SUPER::bnot($class,@_);
369# }
370
371sub bcmp
372 {
373 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
374 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
f9a08e12 375
376 # set up parameters
377 my ($self,$x,$y) = (ref($_[0]),@_);
378 # objectify is costly, so avoid it
379 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
380 {
381 ($self,$x,$y) = objectify(2,@_);
382 }
58cde26e 383
56d9de68 384 return $upgrade->bcmp($x,$y) if defined $upgrade &&
385 ((!$x->isa($self)) || (!$y->isa($self)));
386
0716bf9b 387 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
388 {
389 # handle +-inf and NaN
390 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
391 return 0 if ($x->{sign} eq $y->{sign}) && ($x->{sign} =~ /^[+-]inf$/);
392 return +1 if $x->{sign} eq '+inf';
393 return -1 if $x->{sign} eq '-inf';
394 return -1 if $y->{sign} eq '+inf';
b3abae2a 395 return +1;
0716bf9b 396 }
397
398 # check sign for speed first
574bacfe 399 return 1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # does also 0 <=> -y
58cde26e 400 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # does also -x <=> 0
401
574bacfe 402 # shortcut
403 my $xz = $x->is_zero();
404 my $yz = $y->is_zero();
405 return 0 if $xz && $yz; # 0 <=> 0
406 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
407 return 1 if $yz && $x->{sign} eq '+'; # +x <=> 0
58cde26e 408
409 # adjust so that exponents are equal
bd05a461 410 my $lxm = $x->{_m}->length();
411 my $lym = $y->{_m}->length();
28df3e88 412 # the numify somewhat limits our length, but makes it much faster
413 my $lx = $lxm + $x->{_e}->numify();
414 my $ly = $lym + $y->{_e}->numify();
415 my $l = $lx - $ly; $l = -$l if $x->{sign} eq '-';
bd05a461 416 return $l <=> 0 if $l != 0;
58cde26e 417
bd05a461 418 # lengths (corrected by exponent) are equal
28df3e88 419 # so make mantissa equal length by padding with zero (shift left)
bd05a461 420 my $diff = $lxm - $lym;
421 my $xm = $x->{_m}; # not yet copy it
422 my $ym = $y->{_m};
423 if ($diff > 0)
424 {
425 $ym = $y->{_m}->copy()->blsft($diff,10);
426 }
427 elsif ($diff < 0)
428 {
429 $xm = $x->{_m}->copy()->blsft(-$diff,10);
430 }
28df3e88 431 my $rc = $xm->bacmp($ym);
58cde26e 432 $rc = -$rc if $x->{sign} eq '-'; # -124 < -123
b3abae2a 433 $rc <=> 0;
58cde26e 434 }
435
436sub bacmp
437 {
438 # Compares 2 values, ignoring their signs.
439 # Returns one of undef, <0, =0, >0. (suitable for sort)
440 # (BFLOAT or num_str, BFLOAT or num_str) return cond_code
f9a08e12 441
442 # set up parameters
443 my ($self,$x,$y) = (ref($_[0]),@_);
444 # objectify is costly, so avoid it
445 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
446 {
447 ($self,$x,$y) = objectify(2,@_);
448 }
ee15d750 449
56d9de68 450 return $upgrade->bacmp($x,$y) if defined $upgrade &&
451 ((!$x->isa($self)) || (!$y->isa($self)));
452
ee15d750 453 # handle +-inf and NaN's
abcfbf51 454 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
ee15d750 455 {
456 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
457 return 0 if ($x->is_inf() && $y->is_inf());
458 return 1 if ($x->is_inf() && !$y->is_inf());
b3abae2a 459 return -1;
ee15d750 460 }
461
462 # shortcut
463 my $xz = $x->is_zero();
464 my $yz = $y->is_zero();
465 return 0 if $xz && $yz; # 0 <=> 0
466 return -1 if $xz && !$yz; # 0 <=> +y
467 return 1 if $yz && !$xz; # +x <=> 0
468
469 # adjust so that exponents are equal
470 my $lxm = $x->{_m}->length();
471 my $lym = $y->{_m}->length();
28df3e88 472 # the numify somewhat limits our length, but makes it much faster
473 my $lx = $lxm + $x->{_e}->numify();
474 my $ly = $lym + $y->{_e}->numify();
394e6ffb 475 my $l = $lx - $ly;
ee15d750 476 return $l <=> 0 if $l != 0;
58cde26e 477
ee15d750 478 # lengths (corrected by exponent) are equal
394e6ffb 479 # so make mantissa equal-length by padding with zero (shift left)
ee15d750 480 my $diff = $lxm - $lym;
481 my $xm = $x->{_m}; # not yet copy it
482 my $ym = $y->{_m};
483 if ($diff > 0)
484 {
485 $ym = $y->{_m}->copy()->blsft($diff,10);
486 }
487 elsif ($diff < 0)
488 {
489 $xm = $x->{_m}->copy()->blsft(-$diff,10);
490 }
28df3e88 491 $xm->bacmp($ym) <=> 0;
58cde26e 492 }
a0d0e21e 493
58cde26e 494sub badd
495 {
496 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
497 # return result as BFLOAT
f9a08e12 498
499 # set up parameters
500 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
501 # objectify is costly, so avoid it
502 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
503 {
504 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
505 }
58cde26e 506
574bacfe 507 # inf and NaN handling
508 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
509 {
510 # NaN first
511 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
13a12e00 512 # inf handling
574bacfe 513 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
514 {
13a12e00 515 # +inf++inf or -inf+-inf => same, rest is NaN
516 return $x if $x->{sign} eq $y->{sign};
517 return $x->bnan();
574bacfe 518 }
56b9c951 519 # +-inf + something => +inf; something +-inf => +-inf
574bacfe 520 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
521 return $x;
522 }
523
8f675a64 524 return $upgrade->badd($x,$y,$a,$p,$r) if defined $upgrade &&
525 ((!$x->isa($self)) || (!$y->isa($self)));
526
58cde26e 527 # speed: no add for 0+y or x+0
28df3e88 528 return $x->bround($a,$p,$r) if $y->is_zero(); # x+0
58cde26e 529 if ($x->is_zero()) # 0+y
530 {
531 # make copy, clobbering up x (modify in place!)
532 $x->{_e} = $y->{_e}->copy();
533 $x->{_m} = $y->{_m}->copy();
534 $x->{sign} = $y->{sign} || $nan;
535 return $x->round($a,$p,$r,$y);
a0d0e21e 536 }
58cde26e 537
538 # take lower of the two e's and adapt m1 to it to match m2
28df3e88 539 my $e = $y->{_e};
091c87b1 540 $e = $MBI->bzero() if !defined $e; # if no BFLOAT ?
541 $e = $e->copy(); # make copy (didn't do it yet)
542 $e->bsub($x->{_e}); # Ye - Xe
58cde26e 543 my $add = $y->{_m}->copy();
091c87b1 544 if ($e->{sign} eq '-') # < 0
58cde26e 545 {
091c87b1 546 $x->{_e} += $e; # need the sign of e
547 $x->{_m}->blsft($e->babs(),10); # destroys copy of _e
58cde26e 548 }
091c87b1 549 elsif (!$e->is_zero()) # > 0
58cde26e 550 {
28df3e88 551 $add->blsft($e,10);
58cde26e 552 }
61f5c3f5 553 # else: both e are the same, so just leave them
554 $x->{_m}->{sign} = $x->{sign}; # fiddle with signs
58cde26e 555 $add->{sign} = $y->{sign};
61f5c3f5 556 $x->{_m} += $add; # finally do add/sub
557 $x->{sign} = $x->{_m}->{sign}; # re-adjust signs
558 $x->{_m}->{sign} = '+'; # mantissa always positiv
559 # delete trailing zeros, then round
091c87b1 560 $x->bnorm()->round($a,$p,$r,$y);
58cde26e 561 }
562
563sub bsub
564 {
0716bf9b 565 # (BigFloat or num_str, BigFloat or num_str) return BigFloat
58cde26e 566 # subtract second arg from first, modify first
f9a08e12 567
568 # set up parameters
569 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
570 # objectify is costly, so avoid it
571 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
572 {
573 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
574 }
a0d0e21e 575
091c87b1 576 # XXX TODO: remove?
28df3e88 577 if ($y->is_zero()) # still round for not adding zero
e745a66c 578 {
28df3e88 579 return $x->round($a,$p,$r);
e745a66c 580 }
091c87b1 581
582 # $x - $y = -$x + $y
583 $y->{sign} =~ tr/+-/-+/; # does nothing for NaN
28df3e88 584 $x->badd($y,$a,$p,$r); # badd does not leave internal zeros
091c87b1 585 $y->{sign} =~ tr/+-/-+/; # refix $y (does nothing for NaN)
e745a66c 586 $x; # already rounded by badd()
58cde26e 587 }
588
589sub binc
590 {
591 # increment arg by one
ee15d750 592 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
e745a66c 593
594 if ($x->{_e}->sign() eq '-')
595 {
596 return $x->badd($self->bone(),$a,$p,$r); # digits after dot
597 }
598
599 if (!$x->{_e}->is_zero())
600 {
601 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
602 $x->{_e}->bzero();
603 }
604 # now $x->{_e} == 0
605 if ($x->{sign} eq '+')
606 {
607 $x->{_m}->binc();
608 return $x->bnorm()->bround($a,$p,$r);
609 }
610 elsif ($x->{sign} eq '-')
611 {
612 $x->{_m}->bdec();
613 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
614 return $x->bnorm()->bround($a,$p,$r);
615 }
616 # inf, nan handling etc
091c87b1 617 $x->badd($self->bone(),$a,$p,$r); # does round
58cde26e 618 }
619
620sub bdec
621 {
622 # decrement arg by one
ee15d750 623 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
e745a66c 624
625 if ($x->{_e}->sign() eq '-')
626 {
627 return $x->badd($self->bone('-'),$a,$p,$r); # digits after dot
628 }
629
630 if (!$x->{_e}->is_zero())
631 {
632 $x->{_m}->blsft($x->{_e},10); # 1e2 => 100
633 $x->{_e}->bzero();
634 }
635 # now $x->{_e} == 0
636 my $zero = $x->is_zero();
637 # <= 0
638 if (($x->{sign} eq '-') || $zero)
639 {
640 $x->{_m}->binc();
641 $x->{sign} = '-' if $zero; # 0 => 1 => -1
642 $x->{sign} = '+' if $x->{_m}->is_zero(); # -1 +1 => -0 => +0
643 return $x->bnorm()->round($a,$p,$r);
644 }
645 # > 0
646 elsif ($x->{sign} eq '+')
647 {
648 $x->{_m}->bdec();
649 return $x->bnorm()->round($a,$p,$r);
650 }
651 # inf, nan handling etc
652 $x->badd($self->bone('-'),$a,$p,$r); # does round
58cde26e 653 }
654
990fb837 655sub DEBUG () { 0; }
656
61f5c3f5 657sub blog
658 {
990fb837 659 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
61f5c3f5 660
990fb837 661 # $base > 0, $base != 1; if $base == undef default to $base == e
662 # $x >= 0
9393ace2 663
b3abae2a 664 # we need to limit the accuracy to protect against overflow
665 my $fallback = 0;
990fb837 666 my ($scale,@params);
667 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
61f5c3f5 668
990fb837 669 # also takes care of the "error in _find_round_parameters?" case
670 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
091c87b1 671
b3abae2a 672 # no rounding at all, so must use fallback
990fb837 673 if (scalar @params == 0)
b3abae2a 674 {
675 # simulate old behaviour
990fb837 676 $params[0] = $self->div_scale(); # and round to it as accuracy
677 $params[1] = undef; # P = undef
678 $scale = $params[0]+4; # at least four more for proper round
679 $params[2] = $r; # round mode by caller or undef
b3abae2a 680 $fallback = 1; # to clear a/p afterwards
681 }
682 else
683 {
684 # the 4 below is empirical, and there might be cases where it is not
685 # enough...
990fb837 686 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
b3abae2a 687 }
61f5c3f5 688
b3abae2a 689 return $x->bzero(@params) if $x->is_one();
990fb837 690 # base not defined => base == Euler's constant e
691 if (defined $base)
692 {
091c87b1 693 # make object, since we don't feed it through objectify() to still get the
990fb837 694 # case of $base == undef
695 $base = $self->new($base) unless ref($base);
696 # $base > 0; $base != 1
697 return $x->bnan() if $base->is_zero() || $base->is_one() ||
698 $base->{sign} ne '+';
091c87b1 699 # if $x == $base, we know the result must be 1.0
990fb837 700 return $x->bone('+',@params) if $x->bcmp($base) == 0;
701 }
61f5c3f5 702
b3abae2a 703 # when user set globals, they would interfere with our calculation, so
56d9de68 704 # disable them and later re-enable them
b3abae2a 705 no strict 'refs';
706 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
707 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
708 # we also need to disable any set A or P on $x (_find_round_parameters took
709 # them already into account), since these would interfere, too
710 delete $x->{_a}; delete $x->{_p};
9393ace2 711 # need to disable $upgrade in BigInt, to avoid deep recursion
b3abae2a 712 local $Math::BigInt::upgrade = undef;
93c87d9d 713 local $Math::BigFloat::downgrade = undef;
990fb837 714
715 # upgrade $x if $x is not a BigFloat (handle BigInt input)
716 if (!$x->isa('Math::BigFloat'))
717 {
718 $x = Math::BigFloat->new($x);
719 $self = ref($x);
720 }
721 # first calculate the log to base e (using reduction by 10 (and probably 2))
722 $self->_log_10($x,$scale);
9393ace2 723
990fb837 724 # and if a different base was requested, convert it
725 if (defined $base)
9393ace2 726 {
93c87d9d 727 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
091c87b1 728 # not ln, but some other base (don't modify $base)
990fb837 729 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
9393ace2 730 }
990fb837 731
091c87b1 732 # shortcut to not run through _find_round_parameters again
990fb837 733 if (defined $params[0])
b3abae2a 734 {
990fb837 735 $x->bround($params[0],$params[2]); # then round accordingly
b3abae2a 736 }
737 else
738 {
990fb837 739 $x->bfround($params[1],$params[2]); # then round accordingly
b3abae2a 740 }
741 if ($fallback)
742 {
743 # clear a/p after round, since user did not request it
744 $x->{_a} = undef; $x->{_p} = undef;
745 }
746 # restore globals
747 $$abr = $ab; $$pbr = $pb;
748
749 $x;
61f5c3f5 750 }
751
990fb837 752sub _log
753 {
091c87b1 754 # internal log function to calculate ln() based on Taylor series.
990fb837 755 # Modifies $x in place.
756 my ($self,$x,$scale) = @_;
757
091c87b1 758 # in case of $x == 1, result is 0
759 return $x->bzero() if $x->is_one();
760
990fb837 761 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
762
763 # u = x-1, v = x+1
764 # _ _
765 # Taylor: | u 1 u^3 1 u^5 |
766 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
767 # |_ v 3 v^3 5 v^5 _|
768
769 # This takes much more steps to calculate the result and is thus not used
770 # u = x-1
771 # _ _
772 # Taylor: | u 1 u^2 1 u^3 |
773 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
774 # |_ x 2 x^2 3 x^3 _|
775
990fb837 776 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
777
778 $v = $x->copy(); $v->binc(); # v = x+1
779 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
780 $x->bdiv($v,$scale); # first term: u/v
781 $below = $v->copy();
782 $over = $u->copy();
783 $u *= $u; $v *= $v; # u^2, v^2
784 $below->bmul($v); # u^3, v^3
785 $over->bmul($u);
786 $factor = $self->new(3); $f = $self->new(2);
787
788 my $steps = 0 if DEBUG;
789 $limit = $self->new("1E-". ($scale-1));
790 while (3 < 5)
791 {
792 # we calculate the next term, and add it to the last
793 # when the next term is below our limit, it won't affect the outcome
794 # anymore, so we stop
795
796 # calculating the next term simple from over/below will result in quite
797 # a time hog if the input has many digits, since over and below will
798 # accumulate more and more digits, and the result will also have many
799 # digits, but in the end it is rounded to $scale digits anyway. So if we
800 # round $over and $below first, we save a lot of time for the division
801 # (not with log(1.2345), but try log (123**123) to see what I mean. This
802 # can introduce a rounding error if the division result would be f.i.
803 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
091c87b1 804 # if we truncated $over and $below we might get 0.12345. Does this matter
805 # for the end result? So we give $over and $below 4 more digits to be
806 # on the safe side (unscientific error handling as usual... :+D
990fb837 807
808 $next = $over->copy->bround($scale+4)->bdiv(
809 $below->copy->bmul($factor)->bround($scale+4),
810 $scale);
811
812## old version:
813## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
814
815 last if $next->bacmp($limit) <= 0;
816
817 delete $next->{_a}; delete $next->{_p};
818 $x->badd($next);
819 #print "step $x\n ($next - $limit = ",$next - $limit,")\n";
820 # calculate things for the next term
821 $over *= $u; $below *= $v; $factor->badd($f);
822 if (DEBUG)
823 {
824 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
825 }
826 }
827 $x->bmul($f); # $x *= 2
828 print "took $steps steps\n" if DEBUG;
829 }
830
831sub _log_10
832 {
091c87b1 833 # Internal log function based on reducing input to the range of 0.1 .. 9.99
834 # and then "correcting" the result to the proper one. Modifies $x in place.
990fb837 835 my ($self,$x,$scale) = @_;
836
837 # taking blog() from numbers greater than 10 takes a *very long* time, so we
838 # break the computation down into parts based on the observation that:
839 # blog(x*y) = blog(x) + blog(y)
840 # We set $y here to multiples of 10 so that $x is below 1 (the smaller $x is
841 # the faster it get's, especially because 2*$x takes about 10 times as long,
842 # so by dividing $x by 10 we make it at least factor 100 faster...)
843
844 # The same observation is valid for numbers smaller than 0.1 (e.g. computing
845 # log(1) is fastest, and the farther away we get from 1, the longer it takes)
846 # so we also 'break' this down by multiplying $x with 10 and subtract the
847 # log(10) afterwards to get the correct result.
848
849 # calculate nr of digits before dot
850 my $dbd = $x->{_m}->length() + $x->{_e}->numify();
851
852 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
853 # infinite recursion
854
855 my $calc = 1; # do some calculation?
856
857 # disable the shortcut for 10, since we need log(10) and this would recurse
858 # infinitely deep
859 if ($x->{_e}->is_one() && $x->{_m}->is_one())
860 {
861 $dbd = 0; # disable shortcut
862 # we can use the cached value in these cases
863 if ($scale <= $LOG_10_A)
864 {
865 $x->bzero(); $x->badd($LOG_10);
866 $calc = 0; # no need to calc, but round
867 }
868 }
091c87b1 869 else
990fb837 870 {
091c87b1 871 # disable the shortcut for 2, since we maybe have it cached
872 if ($x->{_e}->is_zero() && $x->{_m}->bcmp(2) == 0)
990fb837 873 {
091c87b1 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 }
990fb837 881 }
882 }
883
884 # if $x = 0.1, we know the result must be 0-log(10)
091c87b1 885 if ($calc != 0 && $x->{_e}->is_one('-') && $x->{_m}->is_one())
990fb837 886 {
887 $dbd = 0; # disable shortcut
888 # we can use the cached value in these cases
889 if ($scale <= $LOG_10_A)
890 {
891 $x->bzero(); $x->bsub($LOG_10);
892 $calc = 0; # no need to calc, but round
893 }
894 }
895
091c87b1 896 return if $calc == 0; # already have the result
897
990fb837 898 # default: these correction factors are undef and thus not used
899 my $l_10; # value of ln(10) to A of $scale
900 my $l_2; # value of ln(2) to A of $scale
901
902 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
903 # so don't do this shortcut for 1 or 0
904 if (($dbd > 1) || ($dbd < 0))
905 {
906 # convert our cached value to an object if not already (avoid doing this
907 # at import() time, since not everybody needs this)
908 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
909
910 #print "x = $x, dbd = $dbd, calc = $calc\n";
911 # got more than one digit before the dot, or more than one zero after the
912 # dot, so do:
913 # log(123) == log(1.23) + log(10) * 2
914 # log(0.0123) == log(1.23) - log(10) * 2
915
916 if ($scale <= $LOG_10_A)
917 {
918 # use cached value
919 #print "using cached value for l_10\n";
920 $l_10 = $LOG_10->copy(); # copy for mul
921 }
922 else
923 {
924 # else: slower, compute it (but don't cache it, because it could be big)
925 # also disable downgrade for this code path
926 local $Math::BigFloat::downgrade = undef;
927 #print "l_10 = $l_10 (self = $self',
928 # ", ref(l_10) = ",ref($l_10)," scale $scale)\n";
929 #print "calculating value for l_10, scale $scale\n";
930 $l_10 = $self->new(10)->blog(undef,$scale); # scale+4, actually
931 }
932 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
933 # make object
934 $dbd = $self->new($dbd);
935 #print "dbd $dbd\n";
936 $l_10->bmul($dbd); # log(10) * (digits_before_dot-1)
937 #print "l_10 = $l_10\n";
938 #print "x = $x";
939 $x->{_e}->bsub($dbd); # 123 => 1.23
940 #print " => $x\n";
941 #print "calculating log($x) with scale=$scale\n";
942
943 }
944
945 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
946
947 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
948 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
949
091c87b1 950 my $half = $self->new('0.5');
951 my $twos = 0; # default: none (0 times)
952 my $two = $self->new(2);
953 while ($x->bacmp($half) <= 0)
990fb837 954 {
091c87b1 955 $twos--; $x->bmul($two);
956 }
957 while ($x->bacmp($two) >= 0)
958 {
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)
990fb837 968 {
091c87b1 969 # use cached value
970 #print "using cached value for l_10\n";
971 $l_2 = $LOG_2->copy(); # copy for mul
990fb837 972 }
091c87b1 973 else
990fb837 974 {
091c87b1 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
990fb837 980 }
091c87b1 981 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
990fb837 982 }
983
091c87b1 984 $self->_log($x,$scale); # need to do the "normal" way
985 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
986 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
990fb837 987 # all done, $x contains now the result
988 }
989
58cde26e 990sub blcm
991 {
ee15d750 992 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
58cde26e 993 # does not modify arguments, but returns new object
994 # Lowest Common Multiplicator
58cde26e 995
996 my ($self,@arg) = objectify(0,@_);
997 my $x = $self->new(shift @arg);
998 while (@arg) { $x = _lcm($x,shift @arg); }
999 $x;
1000 }
1001
1002sub bgcd
1003 {
ee15d750 1004 # (BFLOAT or num_str, BFLOAT or num_str) return BINT
58cde26e 1005 # does not modify arguments, but returns new object
1006 # GCD -- Euclids algorithm Knuth Vol 2 pg 296
58cde26e 1007
1008 my ($self,@arg) = objectify(0,@_);
1009 my $x = $self->new(shift @arg);
1010 while (@arg) { $x = _gcd($x,shift @arg); }
1011 $x;
1012 }
1013
b3abae2a 1014###############################################################################
1015# is_foo methods (is_negative, is_positive are inherited from BigInt)
1016
091c87b1 1017sub _is_zero_or_one
1018 {
1019 # internal, return true if BigInt arg is zero or one, saving the
1020 # two calls to is_zero() and is_one()
1021 my $x = $_[0];
1022
1023 $x->{sign} eq '+' && ($x->is_zero() || $x->is_one());
1024 }
1025
b3abae2a 1026sub is_int
1027 {
1028 # return true if arg (BFLOAT or num_str) is an integer
091c87b1 1029 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
b3abae2a 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
091c87b1 1039 my ($self,$x) = ref($_[0]) ? (undef,$_[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
091c87b1 1048 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : 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
091c87b1 1059 my ($self,$x) = ref($_[0]) ? (undef,$_[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
091c87b1 1069 my ($self,$x) = ref($_[0]) ? (undef,$_[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
091c87b1 1194 # shortcut to not run through _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
091c87b1 1386 # shortcut to not run through _find_round_parameters again
990fb837 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();
091c87b1 1457 # shortcut to not run through _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
091c87b1 1521 # shortcut to not run through _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
091c87b1 1543 # compute factorial number, modifies first argument
b3abae2a 1544 my ($self,$x,@r) = objectify(1,@_);
1545
28df3e88 1546 return $x->bnan()
1547 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
1548 ($x->{_e}->{sign} ne '+')); # digits after dot?
b3abae2a 1549
b3abae2a 1550 # use BigInt's bfac() for faster calc
091c87b1 1551 if (! _is_zero_or_one($x->{_e}))
1552 {
1553 $x->{_m}->blsft($x->{_e},10); # unnorm
1554 $x->{_e}->bzero(); # norm again
1555 }
b3abae2a 1556 $x->{_m}->blsft($x->{_e},10); # un-norm m
091c87b1 1557 $x->{_e}->bzero(); # norm again
1558 $x->{_m}->bfac(); # calculate factorial
1559 $x->bnorm()->round(@r); # norm again and round result
b3abae2a 1560 }
1561
9393ace2 1562sub _pow
1563 {
1564 # Calculate a power where $y is a non-integer, like 2 ** 0.5
1565 my ($x,$y,$a,$p,$r) = @_;
1566 my $self = ref($x);
1567
1568 # if $y == 0.5, it is sqrt($x)
1569 return $x->bsqrt($a,$p,$r,$y) if $y->bcmp('0.5') == 0;
1570
990fb837 1571 # Using:
1572 # a ** x == e ** (x * ln a)
1573
9393ace2 1574 # u = y * ln x
990fb837 1575 # _ _
1576 # Taylor: | u u^2 u^3 |
1577 # x ** y = 1 + | --- + --- + ----- + ... |
1578 # |_ 1 1*2 1*2*3 _|
9393ace2 1579
1580 # we need to limit the accuracy to protect against overflow
1581 my $fallback = 0;
990fb837 1582 my ($scale,@params);
1583 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1584
1585 return $x if $x->is_nan(); # error in _find_round_parameters?
9393ace2 1586
1587 # no rounding at all, so must use fallback
990fb837 1588 if (scalar @params == 0)
9393ace2 1589 {
1590 # simulate old behaviour
990fb837 1591 $params[0] = $self->div_scale(); # and round to it as accuracy
1592 $params[1] = undef; # disable P
1593 $scale = $params[0]+4; # at least four more for proper round
1594 $params[2] = $r; # round mode by caller or undef
9393ace2 1595 $fallback = 1; # to clear a/p afterwards
1596 }
1597 else
1598 {
1599 # the 4 below is empirical, and there might be cases where it is not
1600 # enough...
990fb837 1601 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
9393ace2 1602 }
1603
1604 # when user set globals, they would interfere with our calculation, so
56d9de68 1605 # disable them and later re-enable them
9393ace2 1606 no strict 'refs';
1607 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1608 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1609 # we also need to disable any set A or P on $x (_find_round_parameters took
1610 # them already into account), since these would interfere, too
1611 delete $x->{_a}; delete $x->{_p};
1612 # need to disable $upgrade in BigInt, to avoid deep recursion
1613 local $Math::BigInt::upgrade = undef;
1614
1615 my ($limit,$v,$u,$below,$factor,$next,$over);
1616
990fb837 1617 $u = $x->copy()->blog(undef,$scale)->bmul($y);
9393ace2 1618 $v = $self->bone(); # 1
1619 $factor = $self->new(2); # 2
1620 $x->bone(); # first term: 1
1621
1622 $below = $v->copy();
1623 $over = $u->copy();
1624
1625 $limit = $self->new("1E-". ($scale-1));
1626 #my $steps = 0;
1627 while (3 < 5)
1628 {
1629 # we calculate the next term, and add it to the last
1630 # when the next term is below our limit, it won't affect the outcome
1631 # anymore, so we stop
1632 $next = $over->copy()->bdiv($below,$scale);
990fb837 1633 last if $next->bacmp($limit) <= 0;
9393ace2 1634 $x->badd($next);
9393ace2 1635 # calculate things for the next term
1636 $over *= $u; $below *= $factor; $factor->binc();
1637 #$steps++;
1638 }
1639
091c87b1 1640 # shortcut to not run through _find_round_parameters again
990fb837 1641 if (defined $params[0])
9393ace2 1642 {
990fb837 1643 $x->bround($params[0],$params[2]); # then round accordingly
9393ace2 1644 }
1645 else
1646 {
990fb837 1647 $x->bfround($params[1],$params[2]); # then round accordingly
9393ace2 1648 }
1649 if ($fallback)
1650 {
1651 # clear a/p after round, since user did not request it
1652 $x->{_a} = undef; $x->{_p} = undef;
1653 }
1654 # restore globals
1655 $$abr = $ab; $$pbr = $pb;
1656 $x;
1657 }
1658
58cde26e 1659sub bpow
1660 {
1661 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
1662 # compute power of two numbers, second arg is used as integer
1663 # modifies first argument
1664
f9a08e12 1665 # set up parameters
1666 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1667 # objectify is costly, so avoid it
1668 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1669 {
1670 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1671 }
58cde26e 1672
0716bf9b 1673 return $x if $x->{sign} =~ /^[+-]inf$/;
58cde26e 1674 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
574bacfe 1675 return $x->bone() if $y->is_zero();
58cde26e 1676 return $x if $x->is_one() || $y->is_one();
9393ace2 1677
d614cd8b 1678 return $x->_pow($y,$a,$p,$r) if !$y->is_int(); # non-integer power
9393ace2 1679
1680 my $y1 = $y->as_number(); # make bigint
394e6ffb 1681 # if ($x == -1)
1682 if ($x->{sign} eq '-' && $x->{_m}->is_one() && $x->{_e}->is_zero())
58cde26e 1683 {
1684 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
0716bf9b 1685 return $y1->is_odd() ? $x : $x->babs(1);
288d023a 1686 }
28df3e88 1687 if ($x->is_zero())
1688 {
1689 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
1690 # 0 ** -y => 1 / (0 ** y) => / 0! (1 / 0 => +inf)
1691 $x->binf();
1692 }
58cde26e 1693
1694 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
1695 $y1->babs();
1696 $x->{_m}->bpow($y1);
1697 $x->{_e}->bmul($y1);
1698 $x->{sign} = $nan if $x->{_m}->{sign} eq $nan || $x->{_e}->{sign} eq $nan;
1699 $x->bnorm();
1700 if ($y->{sign} eq '-')
1701 {
1702 # modify $x in place!
0716bf9b 1703 my $z = $x->copy(); $x->bzero()->binc();
58cde26e 1704 return $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
a0d0e21e 1705 }
28df3e88 1706 $x->round($a,$p,$r,$y);
58cde26e 1707 }
1708
1709###############################################################################
1710# rounding functions
1711
1712sub bfround
1713 {
1714 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
1715 # $n == 0 means round to integer
1716 # expects and returns normalized numbers!
ee15d750 1717 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
a0d0e21e 1718
58cde26e 1719 return $x if $x->modify('bfround');
1720
ee15d750 1721 my ($scale,$mode) = $x->_scale_p($self->precision(),$self->round_mode(),@_);
58cde26e 1722 return $x if !defined $scale; # no-op
1723
574bacfe 1724 # never round a 0, +-inf, NaN
61f5c3f5 1725 if ($x->is_zero())
1726 {
1727 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
1728 return $x;
1729 }
1730 return $x if $x->{sign} !~ /^[+-]$/;
58cde26e 1731
ee15d750 1732 # don't round if x already has lower precision
1733 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
1734
1735 $x->{_p} = $scale; # remember round in any case
1736 $x->{_a} = undef; # and clear A
58cde26e 1737 if ($scale < 0)
1738 {
58cde26e 1739 # round right from the '.'
f9a08e12 1740
1741 return $x if $x->{_e}->{sign} eq '+'; # e >= 0 => nothing to round
1742
58cde26e 1743 $scale = -$scale; # positive for simplicity
1744 my $len = $x->{_m}->length(); # length of mantissa
f9a08e12 1745
1746 # the following poses a restriction on _e, but if _e is bigger than a
1747 # scalar, you got other problems (memory etc) anyway
1748 my $dad = -($x->{_e}->numify()); # digits after dot
58cde26e 1749 my $zad = 0; # zeros after dot
f9a08e12 1750 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
1751
ee15d750 1752 #print "scale $scale dad $dad zad $zad len $len\n";
58cde26e 1753 # number bsstr len zad dad
1754 # 0.123 123e-3 3 0 3
1755 # 0.0123 123e-4 3 1 4
1756 # 0.001 1e-3 1 2 3
1757 # 1.23 123e-2 3 0 2
1758 # 1.2345 12345e-4 5 0 4
1759
1760 # do not round after/right of the $dad
1761 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
1762
ee15d750 1763 # round to zero if rounding inside the $zad, but not for last zero like:
1764 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
1765 return $x->bzero() if $scale < $zad;
1766 if ($scale == $zad) # for 0.006, scale -3 and trunc
58cde26e 1767 {
b3abae2a 1768 $scale = -$len;
58cde26e 1769 }
1770 else
1771 {
1772 # adjust round-point to be inside mantissa
1773 if ($zad != 0)
1774 {
1775 $scale = $scale-$zad;
1776 }
1777 else
1778 {
1779 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
1780 $scale = $dbd+$scale;
1781 }
1782 }
a0d0e21e 1783 }
58cde26e 1784 else
1785 {
f9a08e12 1786 # round left from the '.'
1787
58cde26e 1788 # 123 => 100 means length(123) = 3 - $scale (2) => 1
a5f75d66 1789
b3abae2a 1790 my $dbt = $x->{_m}->length();
1791 # digits before dot
f9a08e12 1792 my $dbd = $dbt + $x->{_e}->numify();
b3abae2a 1793 # should be the same, so treat it as this
1794 $scale = 1 if $scale == 0;
1795 # shortcut if already integer
1796 return $x if $scale == 1 && $dbt <= $dbd;
1797 # maximum digits before dot
1798 ++$dbd;
1799
1800 if ($scale > $dbd)
1801 {
1802 # not enough digits before dot, so round to zero
1803 return $x->bzero;
1804 }
1805 elsif ( $scale == $dbd )
1806 {
1807 # maximum
1808 $scale = -$dbt;
1809 }
58cde26e 1810 else
b3abae2a 1811 {
1812 $scale = $dbd - $scale;
1813 }
a0d0e21e 1814 }
574bacfe 1815 # pass sign to bround for rounding modes '+inf' and '-inf'
58cde26e 1816 $x->{_m}->{sign} = $x->{sign};
1817 $x->{_m}->bround($scale,$mode);
1818 $x->{_m}->{sign} = '+'; # fix sign back
1819 $x->bnorm();
1820 }
1821
1822sub bround
1823 {
1824 # accuracy: preserve $N digits, and overwrite the rest with 0's
ee15d750 1825 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
1826
990fb837 1827 if (($_[0] || 0) < 0)
1828 {
1829 require Carp; Carp::croak ('bround() needs positive accuracy');
1830 }
58cde26e 1831
ee15d750 1832 my ($scale,$mode) = $x->_scale_a($self->accuracy(),$self->round_mode(),@_);
1833 return $x if !defined $scale; # no-op
61f5c3f5 1834
58cde26e 1835 return $x if $x->modify('bround');
61f5c3f5 1836
ee15d750 1837 # scale is now either $x->{_a}, $accuracy, or the user parameter
1838 # test whether $x already has lower accuracy, do nothing in this case
1839 # but do round if the accuracy is the same, since a math operation might
1840 # want to round a number with A=5 to 5 digits afterwards again
1841 return $x if defined $_[0] && defined $x->{_a} && $x->{_a} < $_[0];
58cde26e 1842
61f5c3f5 1843 # scale < 0 makes no sense
1844 # never round a +-inf, NaN
1845 return $x if ($scale < 0) || $x->{sign} !~ /^[+-]$/;
58cde26e 1846
61f5c3f5 1847 # 1: $scale == 0 => keep all digits
1848 # 2: never round a 0
1849 # 3: if we should keep more digits than the mantissa has, do nothing
1850 if ($scale == 0 || $x->is_zero() || $x->{_m}->length() <= $scale)
1851 {
1852 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
1853 return $x;
1854 }
f216259d 1855
58cde26e 1856 # pass sign to bround for '+inf' and '-inf' rounding modes
1857 $x->{_m}->{sign} = $x->{sign};
1858 $x->{_m}->bround($scale,$mode); # round mantissa
1859 $x->{_m}->{sign} = '+'; # fix sign back
61f5c3f5 1860 # $x->{_m}->{_a} = undef; $x->{_m}->{_p} = undef;
ee15d750 1861 $x->{_a} = $scale; # remember rounding
1862 $x->{_p} = undef; # and clear P
574bacfe 1863 $x->bnorm(); # del trailing zeros gen. by bround()
58cde26e 1864 }
1865
1866sub bfloor
1867 {
1868 # return integer less or equal then $x
ee15d750 1869 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 1870
1871 return $x if $x->modify('bfloor');
1872
1873 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1874
1875 # if $x has digits after dot
1876 if ($x->{_e}->{sign} eq '-')
1877 {
28df3e88 1878 $x->{_e}->{sign} = '+'; # negate e
1879 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1880 $x->{_e}->bzero(); # trunc/norm
1881 $x->{_m}->binc() if $x->{sign} eq '-'; # decrement if negative
f216259d 1882 }
61f5c3f5 1883 $x->round($a,$p,$r);
58cde26e 1884 }
288d023a 1885
58cde26e 1886sub bceil
1887 {
1888 # return integer greater or equal then $x
ee15d750 1889 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 1890
1891 return $x if $x->modify('bceil');
1892 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1893
1894 # if $x has digits after dot
1895 if ($x->{_e}->{sign} eq '-')
1896 {
28df3e88 1897 #$x->{_m}->brsft(-$x->{_e},10);
1898 #$x->{_e}->bzero();
1899 #$x++ if $x->{sign} eq '+';
1900
1901 $x->{_e}->{sign} = '+'; # negate e
1902 $x->{_m}->brsft($x->{_e},10); # cut off digits after dot
1903 $x->{_e}->bzero(); # trunc/norm
1904 $x->{_m}->binc() if $x->{sign} eq '+'; # decrement if negative
a0d0e21e 1905 }
61f5c3f5 1906 $x->round($a,$p,$r);
58cde26e 1907 }
1908
394e6ffb 1909sub brsft
1910 {
f9a08e12 1911 # shift right by $y (divide by power of $n)
1912
1913 # set up parameters
1914 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1915 # objectify is costly, so avoid it
1916 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1917 {
1918 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1919 }
394e6ffb 1920
1921 return $x if $x->modify('brsft');
1922 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1923
f9a08e12 1924 $n = 2 if !defined $n; $n = $self->new($n);
1925 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
394e6ffb 1926 }
1927
1928sub blsft
1929 {
f9a08e12 1930 # shift left by $y (multiply by power of $n)
1931
1932 # set up parameters
1933 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
1934 # objectify is costly, so avoid it
1935 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1936 {
1937 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
1938 }
394e6ffb 1939
f9a08e12 1940 return $x if $x->modify('blsft');
394e6ffb 1941 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
1942
f9a08e12 1943 $n = 2 if !defined $n; $n = $self->new($n);
1944 $x->bmul($n->bpow($y),$a,$p,$r,$y);
394e6ffb 1945 }
1946
58cde26e 1947###############################################################################
a5f75d66 1948
58cde26e 1949sub DESTROY
1950 {
ee15d750 1951 # going through AUTOLOAD for every DESTROY is costly, so avoid it by empty sub
58cde26e 1952 }
1953
1954sub AUTOLOAD
1955 {
b3abae2a 1956 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
1957 # or falling back to MBI::bxxx()
58cde26e 1958 my $name = $AUTOLOAD;
1959
1960 $name =~ s/.*:://; # split package
ee15d750 1961 no strict 'refs';
990fb837 1962 $class->import() if $IMPORT == 0;
ee15d750 1963 if (!method_alias($name))
58cde26e 1964 {
ee15d750 1965 if (!defined $name)
1966 {
1967 # delayed load of Carp and avoid recursion
1968 require Carp;
1969 Carp::croak ("Can't call a method without name");
1970 }
ee15d750 1971 if (!method_hand_up($name))
1972 {
1973 # delayed load of Carp and avoid recursion
1974 require Carp;
1975 Carp::croak ("Can't call $class\-\>$name, not a valid method");
1976 }
1977 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
1978 $name =~ s/^f/b/;
56b9c951 1979 return &{"$MBI"."::$name"}(@_);
a0d0e21e 1980 }
58cde26e 1981 my $bname = $name; $bname =~ s/^f/b/;
b3abae2a 1982 *{$class."::$name"} = \&$bname;
58cde26e 1983 &$bname; # uses @_
1984 }
1985
1986sub exponent
1987 {
1988 # return a copy of the exponent
ee15d750 1989 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 1990
ee15d750 1991 if ($x->{sign} !~ /^[+-]$/)
1992 {
1993 my $s = $x->{sign}; $s =~ s/^[+-]//;
1994 return $self->new($s); # -inf, +inf => +inf
1995 }
1996 return $x->{_e}->copy();
58cde26e 1997 }
1998
1999sub mantissa
2000 {
2001 # return a copy of the mantissa
ee15d750 2002 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2003
ee15d750 2004 if ($x->{sign} !~ /^[+-]$/)
2005 {
2006 my $s = $x->{sign}; $s =~ s/^[+]//;
2007 return $self->new($s); # -inf, +inf => +inf
2008 }
2009 my $m = $x->{_m}->copy(); # faster than going via bstr()
2010 $m->bneg() if $x->{sign} eq '-';
58cde26e 2011
61f5c3f5 2012 $m;
58cde26e 2013 }
2014
2015sub parts
2016 {
2017 # return a copy of both the exponent and the mantissa
ee15d750 2018 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2019
ee15d750 2020 if ($x->{sign} !~ /^[+-]$/)
2021 {
2022 my $s = $x->{sign}; $s =~ s/^[+]//; my $se = $s; $se =~ s/^[-]//;
2023 return ($self->new($s),$self->new($se)); # +inf => inf and -inf,+inf => inf
2024 }
2025 my $m = $x->{_m}->copy(); # faster than going via bstr()
2026 $m->bneg() if $x->{sign} eq '-';
2027 return ($m,$x->{_e}->copy());
58cde26e 2028 }
2029
2030##############################################################################
2031# private stuff (internal use only)
2032
58cde26e 2033sub import
2034 {
2035 my $self = shift;
8f675a64 2036 my $l = scalar @_;
2037 my $lib = ''; my @a;
990fb837 2038 $IMPORT=1;
8f675a64 2039 for ( my $i = 0; $i < $l ; $i++)
58cde26e 2040 {
2041 if ( $_[$i] eq ':constant' )
2042 {
091c87b1 2043 # This causes overlord er load to step in. 'binary' and 'integer'
2044 # are handled by BigInt.
58cde26e 2045 overload::constant float => sub { $self->new(shift); };
b3abae2a 2046 }
2047 elsif ($_[$i] eq 'upgrade')
2048 {
2049 # this causes upgrading
28df3e88 2050 $upgrade = $_[$i+1]; # or undef to disable
8f675a64 2051 $i++;
28df3e88 2052 }
2053 elsif ($_[$i] eq 'downgrade')
2054 {
2055 # this causes downgrading
2056 $downgrade = $_[$i+1]; # or undef to disable
8f675a64 2057 $i++;
58cde26e 2058 }
56b9c951 2059 elsif ($_[$i] eq 'lib')
2060 {
990fb837 2061 # alternative library
56b9c951 2062 $lib = $_[$i+1] || ''; # default Calc
8f675a64 2063 $i++;
56b9c951 2064 }
2065 elsif ($_[$i] eq 'with')
2066 {
990fb837 2067 # alternative class for our private parts()
56b9c951 2068 $MBI = $_[$i+1] || 'Math::BigInt'; # default Math::BigInt
8f675a64 2069 $i++;
2070 }
2071 else
2072 {
2073 push @a, $_[$i];
56b9c951 2074 }
58cde26e 2075 }
8f675a64 2076
56b9c951 2077 # let use Math::BigInt lib => 'GMP'; use Math::BigFloat; still work
2078 my $mbilib = eval { Math::BigInt->config()->{lib} };
8f675a64 2079 if ((defined $mbilib) && ($MBI eq 'Math::BigInt'))
2080 {
2081 # MBI already loaded
2082 $MBI->import('lib',"$lib,$mbilib", 'objectify');
2083 }
2084 else
2085 {
2086 # MBI not loaded, or with ne "Math::BigInt"
2087 $lib .= ",$mbilib" if defined $mbilib;
07d34614 2088 $lib =~ s/^,//; # don't leave empty
990fb837 2089 # replacement library can handle lib statement, but also could ignore it
8f675a64 2090 if ($] < 5.006)
2091 {
2092 # Perl < 5.6.0 dies with "out of memory!" when eval() and ':constant' is
2093 # used in the same script, or eval inside import().
2094 my @parts = split /::/, $MBI; # Math::BigInt => Math BigInt
2095 my $file = pop @parts; $file .= '.pm'; # BigInt => BigInt.pm
07d34614 2096 require File::Spec;
8f675a64 2097 $file = File::Spec->catfile (@parts, $file);
07d34614 2098 eval { require "$file"; };
2099 $MBI->import( lib => $lib, 'objectify' );
8f675a64 2100 }
2101 else
2102 {
2103 my $rc = "use $MBI lib => '$lib', 'objectify';";
2104 eval $rc;
2105 }
2106 }
990fb837 2107 if ($@)
2108 {
2109 require Carp; Carp::croak ("Couldn't load $MBI: $! $@");
2110 }
56b9c951 2111
58cde26e 2112 # any non :constant stuff is handled by our parent, Exporter
2113 # even if @_ is empty, to give it a chance
b3abae2a 2114 $self->SUPER::import(@a); # for subclasses
2115 $self->export_to_level(1,$self,@a); # need this, too
58cde26e 2116 }
2117
2118sub bnorm
2119 {
2120 # adjust m and e so that m is smallest possible
2121 # round number according to accuracy and precision settings
ee15d750 2122 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2123
0716bf9b 2124 return $x if $x->{sign} !~ /^[+-]$/; # inf, nan etc
58cde26e 2125
28df3e88 2126# if (!$x->{_m}->is_odd())
2127# {
2128 my $zeros = $x->{_m}->_trailing_zeros(); # correct for trailing zeros
2129 if ($zeros != 0)
2130 {
2131 $x->{_m}->brsft($zeros,10); $x->{_e}->badd($zeros);
2132 }
2133 # for something like 0Ey, set y to 1, and -0 => +0
2134 $x->{sign} = '+', $x->{_e}->bone() if $x->{_m}->is_zero();
2135# }
ee15d750 2136 # this is to prevent automatically rounding when MBI's globals are set
0716bf9b 2137 $x->{_m}->{_f} = MB_NEVER_ROUND;
2138 $x->{_e}->{_f} = MB_NEVER_ROUND;
ee15d750 2139 # 'forget' that mantissa was rounded via MBI::bround() in MBF's bfround()
2140 $x->{_m}->{_a} = undef; $x->{_e}->{_a} = undef;
2141 $x->{_m}->{_p} = undef; $x->{_e}->{_p} = undef;
61f5c3f5 2142 $x; # MBI bnorm is no-op, so dont call it
2143 }
58cde26e 2144
2145##############################################################################
56d9de68 2146
2147sub as_hex
2148 {
2149 # return number as hexadecimal string (only for integers defined)
2150 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2151
2152 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2153 return '0x0' if $x->is_zero();
2154
990fb837 2155 return $nan if $x->{_e}->{sign} ne '+'; # how to do 1e-1 in hex!?
56d9de68 2156
2157 my $z = $x->{_m}->copy();
2158 if (!$x->{_e}->is_zero()) # > 0
2159 {
2160 $z->blsft($x->{_e},10);
2161 }
2162 $z->{sign} = $x->{sign};
2163 $z->as_hex();
2164 }
2165
2166sub as_bin
2167 {
2168 # return number as binary digit string (only for integers defined)
2169 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
2170
2171 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
2172 return '0b0' if $x->is_zero();
2173
990fb837 2174 return $nan if $x->{_e}->{sign} ne '+'; # how to do 1e-1 in hex!?
56d9de68 2175
2176 my $z = $x->{_m}->copy();
2177 if (!$x->{_e}->is_zero()) # > 0
2178 {
2179 $z->blsft($x->{_e},10);
2180 }
2181 $z->{sign} = $x->{sign};
2182 $z->as_bin();
2183 }
58cde26e 2184
2185sub as_number
2186 {
394e6ffb 2187 # return copy as a bigint representation of this BigFloat number
2188 my ($self,$x) = ref($_[0]) ? (ref($_[0]),$_[0]) : objectify(1,@_);
58cde26e 2189
28df3e88 2190 my $z = $x->{_m}->copy();
2191 if ($x->{_e}->{sign} eq '-') # < 0
58cde26e 2192 {
28df3e88 2193 $x->{_e}->{sign} = '+'; # flip
2194 $z->brsft($x->{_e},10);
2195 $x->{_e}->{sign} = '-'; # flip back
0716bf9b 2196 }
28df3e88 2197 elsif (!$x->{_e}->is_zero()) # > 0
0716bf9b 2198 {
2199 $z->blsft($x->{_e},10);
58cde26e 2200 }
58cde26e 2201 $z->{sign} = $x->{sign};
61f5c3f5 2202 $z;
58cde26e 2203 }
2204
2205sub length
2206 {
ee15d750 2207 my $x = shift;
2208 my $class = ref($x) || $x;
2209 $x = $class->new(shift) unless ref($x);
58cde26e 2210
ee15d750 2211 return 1 if $x->{_m}->is_zero();
58cde26e 2212 my $len = $x->{_m}->length();
2213 $len += $x->{_e} if $x->{_e}->sign() eq '+';
2214 if (wantarray())
2215 {
56b9c951 2216 my $t = $MBI->bzero();
58cde26e 2217 $t = $x->{_e}->copy()->babs() if $x->{_e}->sign() eq '-';
2218 return ($len,$t);
2219 }
61f5c3f5 2220 $len;
58cde26e 2221 }
a0d0e21e 2222
22231;
a5f75d66 2224__END__
2225
2226=head1 NAME
2227
58cde26e 2228Math::BigFloat - Arbitrary size floating point math package
a5f75d66 2229
2230=head1 SYNOPSIS
2231
a2008d6d 2232 use Math::BigFloat;
58cde26e 2233
b3abae2a 2234 # Number creation
2235 $x = Math::BigFloat->new($str); # defaults to 0
2236 $nan = Math::BigFloat->bnan(); # create a NotANumber
2237 $zero = Math::BigFloat->bzero(); # create a +0
2238 $inf = Math::BigFloat->binf(); # create a +inf
2239 $inf = Math::BigFloat->binf('-'); # create a -inf
2240 $one = Math::BigFloat->bone(); # create a +1
2241 $one = Math::BigFloat->bone('-'); # create a -1
58cde26e 2242
2243 # Testing
b3abae2a 2244 $x->is_zero(); # true if arg is +0
2245 $x->is_nan(); # true if arg is NaN
0716bf9b 2246 $x->is_one(); # true if arg is +1
2247 $x->is_one('-'); # true if arg is -1
2248 $x->is_odd(); # true if odd, false for even
2249 $x->is_even(); # true if even, false for odd
2250 $x->is_positive(); # true if >= 0
2251 $x->is_negative(); # true if < 0
b3abae2a 2252 $x->is_inf(sign); # true if +inf, or -inf (default is '+')
2253
58cde26e 2254 $x->bcmp($y); # compare numbers (undef,<0,=0,>0)
2255 $x->bacmp($y); # compare absolutely (undef,<0,=0,>0)
2256 $x->sign(); # return the sign, either +,- or NaN
b3abae2a 2257 $x->digit($n); # return the nth digit, counting from right
2258 $x->digit(-$n); # return the nth digit, counting from left
58cde26e 2259
990fb837 2260 # The following all modify their first argument. If you want to preserve
2261 # $x, use $z = $x->copy()->bXXX($y); See under L<CAVEATS> for why this is
2262 # neccessary when mixing $a = $b assigments with non-overloaded math.
2263
58cde26e 2264 # set
2265 $x->bzero(); # set $i to 0
2266 $x->bnan(); # set $i to NaN
b3abae2a 2267 $x->bone(); # set $x to +1
2268 $x->bone('-'); # set $x to -1
2269 $x->binf(); # set $x to inf
2270 $x->binf('-'); # set $x to -inf
58cde26e 2271
2272 $x->bneg(); # negation
2273 $x->babs(); # absolute value
2274 $x->bnorm(); # normalize (no-op)
2275 $x->bnot(); # two's complement (bit wise not)
2276 $x->binc(); # increment x by 1
2277 $x->bdec(); # decrement x by 1
2278
2279 $x->badd($y); # addition (add $y to $x)
2280 $x->bsub($y); # subtraction (subtract $y from $x)
2281 $x->bmul($y); # multiplication (multiply $x by $y)
990fb837 2282 $x->bdiv($y); # divide, set $x to quotient
58cde26e 2283 # return (quo,rem) or quo if scalar
2284
990fb837 2285 $x->bmod($y); # modulus ($x % $y)
2286 $x->bpow($y); # power of arguments ($x ** $y)
58cde26e 2287 $x->blsft($y); # left shift
2288 $x->brsft($y); # right shift
2289 # return (quo,rem) or quo if scalar
2290
990fb837 2291 $x->blog(); # logarithm of $x to base e (Euler's number)
2292 $x->blog($base); # logarithm of $x to base $base (f.i. 2)
61f5c3f5 2293
58cde26e 2294 $x->band($y); # bit-wise and
2295 $x->bior($y); # bit-wise inclusive or
2296 $x->bxor($y); # bit-wise exclusive or
2297 $x->bnot(); # bit-wise not (two's complement)
b3abae2a 2298
2299 $x->bsqrt(); # calculate square-root
990fb837 2300 $x->broot($y); # $y'th root of $x (e.g. $y == 3 => cubic root)
b3abae2a 2301 $x->bfac(); # factorial of $x (1*2*3*4*..$x)
2302
990fb837 2303 $x->bround($N); # accuracy: preserve $N digits
58cde26e 2304 $x->bfround($N); # precision: round to the $Nth digit
2305
990fb837 2306 $x->bfloor(); # return integer less or equal than $x
2307 $x->bceil(); # return integer greater or equal than $x
2308
58cde26e 2309 # The following do not modify their arguments:
990fb837 2310
58cde26e 2311 bgcd(@values); # greatest common divisor
2312 blcm(@values); # lowest common multiplicator
2313
2314 $x->bstr(); # return string
2315 $x->bsstr(); # return string in scientific notation
b3abae2a 2316
58cde26e 2317 $x->exponent(); # return exponent as BigInt
2318 $x->mantissa(); # return mantissa as BigInt
2319 $x->parts(); # return (mantissa,exponent) as BigInt
2320
2321 $x->length(); # number of digits (w/o sign and '.')
2322 ($l,$f) = $x->length(); # number of digits, and length of fraction
a5f75d66 2323
f9a08e12 2324 $x->precision(); # return P of $x (or global, if P of $x undef)
2325 $x->precision($n); # set P of $x to $n
2326 $x->accuracy(); # return A of $x (or global, if A of $x undef)
723d369b 2327 $x->accuracy($n); # set A $x to $n
f9a08e12 2328
990fb837 2329 # these get/set the appropriate global value for all BigFloat objects
2330 Math::BigFloat->precision(); # Precision
2331 Math::BigFloat->accuracy(); # Accuracy
2332 Math::BigFloat->round_mode(); # rounding mode
f9a08e12 2333
a5f75d66 2334=head1 DESCRIPTION
2335
58cde26e 2336All operators (inlcuding basic math operations) are overloaded if you
2337declare your big floating point numbers as
a5f75d66 2338
58cde26e 2339 $i = new Math::BigFloat '12_3.456_789_123_456_789E-2';
2340
2341Operations with overloaded operators preserve the arguments, which is
2342exactly what you expect.
2343
2344=head2 Canonical notation
2345
2346Input to these routines are either BigFloat objects, or strings of the
2347following four forms:
a5f75d66 2348
2349=over 2
2350
58cde26e 2351=item *
2352
2353C</^[+-]\d+$/>
a5f75d66 2354
58cde26e 2355=item *
a5f75d66 2356
58cde26e 2357C</^[+-]\d+\.\d*$/>
a5f75d66 2358
58cde26e 2359=item *
a5f75d66 2360
58cde26e 2361C</^[+-]\d+E[+-]?\d+$/>
a5f75d66 2362
58cde26e 2363=item *
a5f75d66 2364
58cde26e 2365C</^[+-]\d*\.\d+E[+-]?\d+$/>
5d7098d5 2366
58cde26e 2367=back
2368
2369all with optional leading and trailing zeros and/or spaces. Additonally,
2370numbers are allowed to have an underscore between any two digits.
2371
2372Empty strings as well as other illegal numbers results in 'NaN'.
2373
2374bnorm() on a BigFloat object is now effectively a no-op, since the numbers
2375are always stored in normalized form. On a string, it creates a BigFloat
2376object.
2377
2378=head2 Output
2379
2380Output values are BigFloat objects (normalized), except for bstr() and bsstr().
2381
2382The string output will always have leading and trailing zeros stripped and drop
2383a plus sign. C<bstr()> will give you always the form with a decimal point,
990fb837 2384while C<bsstr()> (s for scientific) gives you the scientific notation.
58cde26e 2385
2386 Input bstr() bsstr()
2387 '-0' '0' '0E1'
2388 ' -123 123 123' '-123123123' '-123123123E0'
2389 '00.0123' '0.0123' '123E-4'
2390 '123.45E-2' '1.2345' '12345E-4'
2391 '10E+3' '10000' '1E4'
2392
2393Some routines (C<is_odd()>, C<is_even()>, C<is_zero()>, C<is_one()>,
2394C<is_nan()>) return true or false, while others (C<bcmp()>, C<bacmp()>)
2395return either undef, <0, 0 or >0 and are suited for sort.
2396
990fb837 2397Actual math is done by using the class defined with C<with => Class;> (which
2398defaults to BigInts) to represent the mantissa and exponent.
2399
58cde26e 2400The sign C</^[+-]$/> is stored separately. The string 'NaN' is used to
2401represent the result when input arguments are not numbers, as well as
2402the result of dividing by zero.
2403
2404=head2 C<mantissa()>, C<exponent()> and C<parts()>
2405
2406C<mantissa()> and C<exponent()> return the said parts of the BigFloat
2407as BigInts such that:
2408
2409 $m = $x->mantissa();
2410 $e = $x->exponent();
2411 $y = $m * ( 10 ** $e );
2412 print "ok\n" if $x == $y;
2413
2414C<< ($m,$e) = $x->parts(); >> is just a shortcut giving you both of them.
2415
2416A zero is represented and returned as C<0E1>, B<not> C<0E0> (after Knuth).
2417
2418Currently the mantissa is reduced as much as possible, favouring higher
2419exponents over lower ones (e.g. returning 1e7 instead of 10e6 or 10000000e0).
2420This might change in the future, so do not depend on it.
2421
2422=head2 Accuracy vs. Precision
2423
2424See also: L<Rounding|Rounding>.
2425
027dc388 2426Math::BigFloat supports both precision and accuracy. For a full documentation,
2427examples and tips on these topics please see the large section in
2428L<Math::BigInt>.
5d7098d5 2429
58cde26e 2430Since things like sqrt(2) or 1/3 must presented with a limited precision lest
2431a operation consumes all resources, each operation produces no more than
990fb837 2432the requested number of digits.
2433
2434Please refer to BigInt's documentation for the precedence rules of which
2435accuracy/precision setting will be used.
2436
2437If there is no gloabl precision set, B<and> the operation inquestion was not
2438called with a requested precision or accuracy, B<and> the input $x has no
2439accuracy or precision set, then a fallback parameter will be used. For
2440historical reasons, it is called C<div_scale> and can be accessed via:
2441
2442 $d = Math::BigFloat->div_scale(); # query
2443 Math::BigFloat->div_scale($n); # set to $n digits
2444
2445The default value is 40 digits.
58cde26e 2446
2447In case the result of one operation has more precision than specified,
2448it is rounded. The rounding mode taken is either the default mode, or the one
2449supplied to the operation after the I<scale>:
2450
2451 $x = Math::BigFloat->new(2);
990fb837 2452 Math::BigFloat->precision(5); # 5 digits max
58cde26e 2453 $y = $x->copy()->bdiv(3); # will give 0.66666
2454 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2455 $y = $x->copy()->bdiv(3,6,'odd'); # will give 0.666667
990fb837 2456 Math::BigFloat->round_mode('zero');
58cde26e 2457 $y = $x->copy()->bdiv(3,6); # will give 0.666666
2458
2459=head2 Rounding
2460
2461=over 2
2462
5dc6f178 2463=item ffround ( +$scale )
58cde26e 2464
0716bf9b 2465Rounds to the $scale'th place left from the '.', counting from the dot.
2466The first digit is numbered 1.
58cde26e 2467
5dc6f178 2468=item ffround ( -$scale )
58cde26e 2469
0716bf9b 2470Rounds to the $scale'th place right from the '.', counting from the dot.
58cde26e 2471
5dc6f178 2472=item ffround ( 0 )
2473
0716bf9b 2474Rounds to an integer.
5dc6f178 2475
2476=item fround ( +$scale )
2477
0716bf9b 2478Preserves accuracy to $scale digits from the left (aka significant digits)
2479and pads the rest with zeros. If the number is between 1 and -1, the
2480significant digits count from the first non-zero after the '.'
5dc6f178 2481
2482=item fround ( -$scale ) and fround ( 0 )
2483
990fb837 2484These are effectively no-ops.
5d7098d5 2485
a5f75d66 2486=back
2487
0716bf9b 2488All rounding functions take as a second parameter a rounding mode from one of
2489the following: 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
58cde26e 2490
2491The default rounding mode is 'even'. By using
990fb837 2492C<< Math::BigFloat->round_mode($round_mode); >> you can get and set the default
ee15d750 2493mode for subsequent rounding. The usage of C<$Math::BigFloat::$round_mode> is
0716bf9b 2494no longer supported.
b22b3e31 2495The second parameter to the round functions then overrides the default
0716bf9b 2496temporarily.
58cde26e 2497
990fb837 2498The C<as_number()> function returns a BigInt from a Math::BigFloat. It uses
58cde26e 2499'trunc' as rounding mode to make it equivalent to:
2500
2501 $x = 2.5;
2502 $y = int($x) + 2;
2503
2504You can override this by passing the desired rounding mode as parameter to
2505C<as_number()>:
2506
2507 $x = Math::BigFloat->new(2.5);
2508 $y = $x->as_number('odd'); # $y = 3
2509
2510=head1 EXAMPLES
2511
58cde26e 2512 # not ready yet
58cde26e 2513
2514=head1 Autocreating constants
2515
2516After C<use Math::BigFloat ':constant'> all the floating point constants
2517in the given scope are converted to C<Math::BigFloat>. This conversion
2518happens at compile time.
2519
2520In particular
2521
2522 perl -MMath::BigFloat=:constant -e 'print 2E-100,"\n"'
2523
56b9c951 2524prints the value of C<2E-100>. Note that without conversion of
58cde26e 2525constants the expression 2E-100 will be calculated as normal floating point
2526number.
2527
56b9c951 2528Please note that ':constant' does not affect integer constants, nor binary
2529nor hexadecimal constants. Use L<bignum> or L<Math::BigInt> to get this to
2530work.
2531
2532=head2 Math library
2533
2534Math with the numbers is done (by default) by a module called
2535Math::BigInt::Calc. This is equivalent to saying:
2536
2537 use Math::BigFloat lib => 'Calc';
2538
2539You can change this by using:
2540
2541 use Math::BigFloat lib => 'BitVect';
2542
2543The following would first try to find Math::BigInt::Foo, then
2544Math::BigInt::Bar, and when this also fails, revert to Math::BigInt::Calc:
2545
2546 use Math::BigFloat lib => 'Foo,Math::BigInt::Bar';
2547
2548Calc.pm uses as internal format an array of elements of some decimal base
2549(usually 1e7, but this might be differen for some systems) with the least
2550significant digit first, while BitVect.pm uses a bit vector of base 2, most
2551significant bit first. Other modules might use even different means of
2552representing the numbers. See the respective module documentation for further
2553details.
2554
2555Please note that Math::BigFloat does B<not> use the denoted library itself,
2556but it merely passes the lib argument to Math::BigInt. So, instead of the need
2557to do:
2558
2559 use Math::BigInt lib => 'GMP';
2560 use Math::BigFloat;
2561
2562you can roll it all into one line:
2563
2564 use Math::BigFloat lib => 'GMP';
2565
990fb837 2566It is also possible to just require Math::BigFloat:
2567
2568 require Math::BigFloat;
2569
2570This will load the neccessary things (like BigInt) when they are needed, and
2571automatically.
2572
2573Use the lib, Luke! And see L<Using Math::BigInt::Lite> for more details than
2574you ever wanted to know about loading a different library.
56b9c951 2575
2576=head2 Using Math::BigInt::Lite
2577
2578It is possible to use L<Math::BigInt::Lite> with Math::BigFloat:
2579
2580 # 1
2581 use Math::BigFloat with => 'Math::BigInt::Lite';
2582
2583There is no need to "use Math::BigInt" or "use Math::BigInt::Lite", but you
2584can combine these if you want. For instance, you may want to use
2585Math::BigInt objects in your main script, too.
2586
2587 # 2
2588 use Math::BigInt;
2589 use Math::BigFloat with => 'Math::BigInt::Lite';
2590
2591Of course, you can combine this with the C<lib> parameter.
2592
2593 # 3
2594 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2595
990fb837 2596There is no need for a "use Math::BigInt;" statement, even if you want to
2597use Math::BigInt's, since Math::BigFloat will needs Math::BigInt and thus
2598always loads it. But if you add it, add it B<before>:
56b9c951 2599
2600 # 4
2601 use Math::BigInt;
2602 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'GMP,Pari';
2603
2604Notice that the module with the last C<lib> will "win" and thus
2605it's lib will be used if the lib is available:
2606
2607 # 5
2608 use Math::BigInt lib => 'Bar,Baz';
2609 use Math::BigFloat with => 'Math::BigInt::Lite', lib => 'Foo';
2610
2611That would try to load Foo, Bar, Baz and Calc (in that order). Or in other
2612words, Math::BigFloat will try to retain previously loaded libs when you
990fb837 2613don't specify it onem but if you specify one, it will try to load them.
56b9c951 2614
2615Actually, the lib loading order would be "Bar,Baz,Calc", and then
2616"Foo,Bar,Baz,Calc", but independend of which lib exists, the result is the
990fb837 2617same as trying the latter load alone, except for the fact that one of Bar or
2618Baz might be loaded needlessly in an intermidiate step (and thus hang around
2619and waste memory). If neither Bar nor Baz exist (or don't work/compile), they
2620will still be tried to be loaded, but this is not as time/memory consuming as
2621actually loading one of them. Still, this type of usage is not recommended due
2622to these issues.
56b9c951 2623
990fb837 2624The old way (loading the lib only in BigInt) still works though:
56b9c951 2625
2626 # 6
2627 use Math::BigInt lib => 'Bar,Baz';
2628 use Math::BigFloat;
2629
990fb837 2630You can even load Math::BigInt afterwards:
56b9c951 2631
990fb837 2632 # 7
2633 use Math::BigFloat;
2634 use Math::BigInt lib => 'Bar,Baz';
a5f75d66 2635
990fb837 2636But this has the same problems like #5, it will first load Calc
2637(Math::BigFloat needs Math::BigInt and thus loads it) and then later Bar or
2638Baz, depending on which of them works and is usable/loadable. Since this
2639loads Calc unnecc., it is not recommended.
58cde26e 2640
990fb837 2641Since it also possible to just require Math::BigFloat, this poses the question
2642about what libary this will use:
58cde26e 2643
990fb837 2644 require Math::BigFloat;
2645 my $x = Math::BigFloat->new(123); $x += 123;
58cde26e 2646
990fb837 2647It will use Calc. Please note that the call to import() is still done, but
2648only when you use for the first time some Math::BigFloat math (it is triggered
2649via any constructor, so the first time you create a Math::BigFloat, the load
2650will happen in the background). This means:
58cde26e 2651
990fb837 2652 require Math::BigFloat;
2653 Math::BigFloat->import ( lib => 'Foo,Bar' );
58cde26e 2654
990fb837 2655would be the same as:
58cde26e 2656
990fb837 2657 use Math::BigFloat lib => 'Foo, Bar';
2658
2659But don't try to be clever to insert some operations in between:
2660
2661 require Math::BigFloat;
2662 my $x = Math::BigFloat->bone() + 4; # load BigInt and Calc
2663 Math::BigFloat->import( lib => 'Pari' ); # load Pari, too
2664 $x = Math::BigFloat->bone()+4; # now use Pari
2665
2666While this works, it loads Calc needlessly. But maybe you just wanted that?
2667
2668B<Examples #3 is highly recommended> for daily usage.
2669
2670=head1 BUGS
2671
2672Please see the file BUGS in the CPAN distribution Math::BigInt for known bugs.
58cde26e 2673
990fb837 2674=head1 CAVEATS
58cde26e 2675
2676=over 1
2677
2678=item stringify, bstr()
2679
2680Both stringify and bstr() now drop the leading '+'. The old code would return
2681'+1.23', the new returns '1.23'. See the documentation in L<Math::BigInt> for
2682reasoning and details.
2683
2684=item bdiv
2685
2686The following will probably not do what you expect:
2687
2688 print $c->bdiv(123.456),"\n";
2689
2690It prints both quotient and reminder since print works in list context. Also,
2691bdiv() will modify $c, so be carefull. You probably want to use
2692
2693 print $c / 123.456,"\n";
2694 print scalar $c->bdiv(123.456),"\n"; # or if you want to modify $c
2695
2696instead.
2697
2698=item Modifying and =
2699
2700Beware of:
2701
2702 $x = Math::BigFloat->new(5);
2703 $y = $x;
2704
2705It will not do what you think, e.g. making a copy of $x. Instead it just makes
2706a second reference to the B<same> object and stores it in $y. Thus anything
990fb837 2707that modifies $x will modify $y (except overloaded math operators), and vice
2708versa. See L<Math::BigInt> for details and how to avoid that.
58cde26e 2709
2710=item bpow
2711
2712C<bpow()> now modifies the first argument, unlike the old code which left
2713it alone and only returned the result. This is to be consistent with
2714C<badd()> etc. The first will modify $x, the second one won't:
2715
2716 print bpow($x,$i),"\n"; # modify $x
2717 print $x->bpow($i),"\n"; # ditto
2718 print $x ** $i,"\n"; # leave $x alone
2719
2720=back
2721
990fb837 2722=head1 SEE ALSO
2723
2724L<Math::BigInt>, L<Math::BigRat> and L<Math::Big> as well as
2725L<Math::BigInt::BitVect>, L<Math::BigInt::Pari> and L<Math::BigInt::GMP>.
2726
2727The pragmas L<bignum>, L<bigint> and L<bigrat> might also be of interest
2728because they solve the autoupgrading/downgrading issue, at least partly.
2729
2730The package at
2731L<http://search.cpan.org/search?mode=module&query=Math%3A%3ABigInt> contains
2732more documentation including a full version history, testcases, empty
2733subclass files and benchmarks.
2734
58cde26e 2735=head1 LICENSE
a5f75d66 2736
58cde26e 2737This program is free software; you may redistribute it and/or modify it under
2738the same terms as Perl itself.
5d7098d5 2739
58cde26e 2740=head1 AUTHORS
5d7098d5 2741
58cde26e 2742Mark Biggar, overloaded interface by Ilya Zakharevich.
990fb837 2743Completely rewritten by Tels http://bloodgate.com in 2001, 2002, and still
2744at it in 2003.
a5f75d66 2745
a5f75d66 2746=cut
Unnamed repository; edit this file to name it for gitweb.