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