source: for-distributions/trunk/bin/windows/perl/lib/Math/BigFloat.pm@ 14489

Last change on this file since 14489 was 14489, checked in by oranfry, 17 years ago

upgrading to perl 5.8

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