This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Clean up POD for Math::BigInt and Math::BigFloat
[perl5.git] / dist / Math-BigInt / lib / Math / BigFloat.pm
CommitLineData
13a12e00
JH
1package Math::BigFloat;
2
3#
d614cd8b 4# Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
13a12e00
JH
5#
6
58cde26e 7# The following hash values are internally used:
9b924220
RGS
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
ada8209b 15$VERSION = '1.999';
0d71d61a 16require 5.006002;
3a427a11
RGS
17
18require Exporter;
fdb4b05f
T
19@ISA = qw/Math::BigInt/;
20@EXPORT_OK = qw/bpi/;
394e6ffb 21
58cde26e 22use strict;
03874afe 23# $_trap_inf/$_trap_nan are internal and should never be accessed from outside
b282a552
T
24use vars qw/$AUTOLOAD $accuracy $precision $div_scale $round_mode $rnd_mode
25 $upgrade $downgrade $_trap_nan $_trap_inf/;
58cde26e 26my $class = "Math::BigFloat";
a0d0e21e 27
a5f75d66 28use overload
a0ac753d 29'<=>' => sub { my $rc = $_[2] ?
bd05a461 30 ref($_[0])->bcmp($_[1],$_[0]) :
a0ac753d
T
31 ref($_[0])->bcmp($_[0],$_[1]);
32 $rc = 1 unless defined $rc;
33 $rc <=> 0;
34 },
35# we need '>=' to get things like "1 >= NaN" right:
36'>=' => sub { my $rc = $_[2] ?
37 ref($_[0])->bcmp($_[1],$_[0]) :
38 ref($_[0])->bcmp($_[0],$_[1]);
39 # if there was a NaN involved, return false
40 return '' unless defined $rc;
41 $rc >= 0;
42 },
0716bf9b 43'int' => sub { $_[0]->as_number() }, # 'trunc' to bigint
a5f75d66 44;
a0d0e21e 45
0716bf9b 46##############################################################################
990fb837 47# global constants, flags and assorted stuff
0716bf9b 48
990fb837
RGS
49# the following are public, but their usage is not recommended. Use the
50# accessor methods instead.
58cde26e 51
ee15d750 52# class constants, use Class->constant_name() to access
20e2035c
T
53# one of 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'
54$round_mode = 'even';
ee15d750
JH
55$accuracy = undef;
56$precision = undef;
57$div_scale = 40;
58cde26e 58
b3abae2a
JH
59$upgrade = undef;
60$downgrade = undef;
9b924220
RGS
61# the package we are using for our private parts, defaults to:
62# Math::BigInt->config()->{lib}
a90064ab 63my $MBI = 'Math::BigInt::Calc';
990fb837
RGS
64
65# are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
66$_trap_nan = 0;
9b924220 67# the same for infinity
990fb837
RGS
68$_trap_inf = 0;
69
70# constant for easier life
71my $nan = 'NaN';
72
9b924220 73my $IMPORT = 0; # was import() called yet? used to make require work
990fb837
RGS
74
75# some digits of accuracy for blog(undef,10); which we use in blog() for speed
76my $LOG_10 =
77 '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
78my $LOG_10_A = length($LOG_10)-1;
79# ditto for log(2)
80my $LOG_2 =
81 '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
82my $LOG_2_A = length($LOG_2)-1;
7b29e1e6 83my $HALF = '0.5'; # made into an object if nec.
990fb837 84
027dc388
JH
85##############################################################################
86# the old code had $rnd_mode, so we need to support it, too
87
027dc388
JH
88sub TIESCALAR { my ($class) = @_; bless \$round_mode, $class; }
89sub FETCH { return $round_mode; }
90sub STORE { $rnd_mode = $_[0]->round_mode($_[1]); }
91
56b9c951 92BEGIN
990fb837 93 {
7d193e39 94 # when someone sets $rnd_mode, we catch this and check the value to see
990fb837 95 # whether it is valid or not.
7b29e1e6
T
96 $rnd_mode = 'even'; tie $rnd_mode, 'Math::BigFloat';
97
98 # we need both of them in this package:
99 *as_int = \&as_number;
56b9c951 100 }
027dc388
JH
101
102##############################################################################
103
58cde26e 104{
ee15d750 105 # valid method aliases for AUTOLOAD
58cde26e
JH
106 my %methods = map { $_ => 1 }
107 qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
7d193e39
T
108 fint facmp fcmp fzero fnan finf finc fdec ffac fneg
109 fceil ffloor frsft flsft fone flog froot fexp
ee15d750 110 /;
7b29e1e6 111 # valid methods that can be handed up (for AUTOLOAD)
ee15d750 112 my %hand_ups = map { $_ => 1 }
ef9466ea 113 qw / is_nan is_inf is_negative is_positive is_pos is_neg
b68b7ab1 114 accuracy precision div_scale round_mode fabs fnot
28df3e88 115 objectify upgrade downgrade
13a12e00 116 bone binf bnan bzero
a0ac753d 117 bsub
58cde26e
JH
118 /;
119
7b29e1e6
T
120 sub _method_alias { exists $methods{$_[0]||''}; }
121 sub _method_hand_up { exists $hand_ups{$_[0]||''}; }
a0d0e21e 122}
0e8b9368 123
58cde26e
JH
124##############################################################################
125# constructors
a0d0e21e 126
58cde26e
JH
127sub new
128 {
129 # create a new BigFloat object from a string or another bigfloat object.
130 # _e: exponent
131 # _m: mantissa
132 # sign => sign (+/-), or "NaN"
a0d0e21e 133
61f5c3f5 134 my ($class,$wanted,@r) = @_;
b3abae2a 135
61f5c3f5
T
136 # avoid numify-calls by not using || on $wanted!
137 return $class->bzero() if !defined $wanted; # default to 0
138 return $wanted->copy() if UNIVERSAL::isa($wanted,'Math::BigFloat');
a0d0e21e 139
990fb837
RGS
140 $class->import() if $IMPORT == 0; # make require work
141
58cde26e 142 my $self = {}; bless $self, $class;
b22b3e31 143 # shortcut for bigints and its subclasses
a0ac753d 144 if ((ref($wanted)) && UNIVERSAL::can( $wanted, "as_number"))
58cde26e 145 {
9b924220
RGS
146 $self->{_m} = $wanted->as_number()->{value}; # get us a bigint copy
147 $self->{_e} = $MBI->_zero();
148 $self->{_es} = '+';
58cde26e 149 $self->{sign} = $wanted->sign();
0716bf9b 150 return $self->bnorm();
58cde26e 151 }
9681bfa6 152 # else: got a string or something masquerading as number (with overload)
2d2b2744 153
58cde26e 154 # handle '+inf', '-inf' first
233f7bc0 155 if ($wanted =~ /^[+-]?inf\z/)
58cde26e 156 {
28df3e88
JH
157 return $downgrade->new($wanted) if $downgrade;
158
233f7bc0
T
159 $self->{sign} = $wanted; # set a default sign for bstr()
160 return $self->binf($wanted);
58cde26e 161 }
b282a552 162
2d2b2744
T
163 # shortcut for simple forms like '12' that neither have trailing nor leading
164 # zeros
165 if ($wanted =~ /^([+-]?)([1-9][0-9]*[1-9])$/)
166 {
167 $self->{_e} = $MBI->_zero();
168 $self->{_es} = '+';
169 $self->{sign} = $1 || '+';
170 $self->{_m} = $MBI->_new($2);
171 return $self->round(@r) if !$downgrade;
172 }
173
9b924220 174 my ($mis,$miv,$mfv,$es,$ev) = Math::BigInt::_split($wanted);
58cde26e
JH
175 if (!ref $mis)
176 {
990fb837
RGS
177 if ($_trap_nan)
178 {
179 require Carp;
180 Carp::croak ("$wanted is not a number initialized to $class");
181 }
28df3e88
JH
182
183 return $downgrade->bnan() if $downgrade;
184
9b924220
RGS
185 $self->{_e} = $MBI->_zero();
186 $self->{_es} = '+';
187 $self->{_m} = $MBI->_zero();
58cde26e
JH
188 $self->{sign} = $nan;
189 }
190 else
191 {
9b924220
RGS
192 # make integer from mantissa by adjusting exp, then convert to int
193 $self->{_e} = $MBI->_new($$ev); # exponent
194 $self->{_es} = $$es || '+';
195 my $mantissa = "$$miv$$mfv"; # create mant.
196 $mantissa =~ s/^0+(\d)/$1/; # strip leading zeros
197 $self->{_m} = $MBI->_new($mantissa); # create mant.
b282a552 198
58cde26e 199 # 3.123E0 = 3123E-3, and 3.123E-2 => 3123E-5
9b924220
RGS
200 if (CORE::length($$mfv) != 0)
201 {
202 my $len = $MBI->_new( CORE::length($$mfv));
203 ($self->{_e}, $self->{_es}) =
204 _e_sub ($self->{_e}, $len, $self->{_es}, '+');
205 }
2d2b2744
T
206 # we can only have trailing zeros on the mantissa if $$mfv eq ''
207 else
b282a552 208 {
2d2b2744
T
209 # Use a regexp to count the trailing zeros in $$miv instead of _zeros()
210 # because that is faster, especially when _m is not stored in base 10.
211 my $zeros = 0; $zeros = CORE::length($1) if $$miv =~ /[1-9](0*)$/;
b282a552
T
212 if ($zeros != 0)
213 {
9b924220 214 my $z = $MBI->_new($zeros);
2d2b2744 215 # turn '120e2' into '12e3'
9b924220 216 $MBI->_rsft ( $self->{_m}, $z, 10);
7596a890
RGS
217 ($self->{_e}, $self->{_es}) =
218 _e_add ( $self->{_e}, $z, $self->{_es}, '+');
b282a552
T
219 }
220 }
2d2b2744
T
221 $self->{sign} = $$mis;
222
9b924220 223 # for something like 0Ey, set y to 1, and -0 => +0
3c4b39be 224 # Check $$miv for being '0' and $$mfv eq '', because otherwise _m could not
2d2b2744 225 # have become 0. That's faster than to call $MBI->_is_zero().
9b924220 226 $self->{sign} = '+', $self->{_e} = $MBI->_one()
2d2b2744
T
227 if $$miv eq '0' and $$mfv eq '';
228
b282a552 229 return $self->round(@r) if !$downgrade;
58cde26e 230 }
28df3e88
JH
231 # if downgrade, inf, NaN or integers go down
232
9b924220 233 if ($downgrade && $self->{_es} eq '+')
28df3e88 234 {
9b924220 235 if ($MBI->_is_zero( $self->{_e} ))
28df3e88 236 {
9b924220 237 return $downgrade->new($$mis . $MBI->_str( $self->{_m} ));
28df3e88 238 }
8df1e0a2 239 return $downgrade->new($self->bsstr());
28df3e88 240 }
990fb837 241 $self->bnorm()->round(@r); # first normalize, then round
58cde26e 242 }
a0d0e21e 243
9b924220
RGS
244sub copy
245 {
86f0d17a 246 # if two arguments, the first one is the class to "swallow" subclasses
9b924220
RGS
247 if (@_ > 1)
248 {
86f0d17a
T
249 my $self = bless {
250 sign => $_[1]->{sign},
251 _es => $_[1]->{_es},
252 _m => $MBI->_copy($_[1]->{_m}),
253 _e => $MBI->_copy($_[1]->{_e}),
254 }, $_[0] if @_ > 1;
255
256 $self->{_a} = $_[1]->{_a} if defined $_[1]->{_a};
257 $self->{_p} = $_[1]->{_p} if defined $_[1]->{_p};
258 return $self;
9b924220 259 }
9b924220 260
86f0d17a
T
261 my $self = bless {
262 sign => $_[0]->{sign},
263 _es => $_[0]->{_es},
264 _m => $MBI->_copy($_[0]->{_m}),
265 _e => $MBI->_copy($_[0]->{_e}),
266 }, ref($_[0]);
9b924220 267
86f0d17a
T
268 $self->{_a} = $_[0]->{_a} if defined $_[0]->{_a};
269 $self->{_p} = $_[0]->{_p} if defined $_[0]->{_p};
9b924220
RGS
270 $self;
271 }
272
13a12e00 273sub _bnan
58cde26e 274 {
990fb837 275 # used by parent class bone() to initialize number to NaN
58cde26e 276 my $self = shift;
990fb837
RGS
277
278 if ($_trap_nan)
279 {
280 require Carp;
281 my $class = ref($self);
282 Carp::croak ("Tried to set $self to NaN in $class\::_bnan()");
283 }
284
285 $IMPORT=1; # call our import only once
9b924220
RGS
286 $self->{_m} = $MBI->_zero();
287 $self->{_e} = $MBI->_zero();
288 $self->{_es} = '+';
58cde26e 289 }
a0d0e21e 290
13a12e00 291sub _binf
58cde26e 292 {
990fb837 293 # used by parent class bone() to initialize number to +-inf
58cde26e 294 my $self = shift;
990fb837
RGS
295
296 if ($_trap_inf)
297 {
298 require Carp;
299 my $class = ref($self);
300 Carp::croak ("Tried to set $self to +-inf in $class\::_binf()");
301 }
302
303 $IMPORT=1; # call our import only once
9b924220
RGS
304 $self->{_m} = $MBI->_zero();
305 $self->{_e} = $MBI->_zero();
306 $self->{_es} = '+';
58cde26e 307 }
a0d0e21e 308
13a12e00 309sub _bone
574bacfe 310 {
13a12e00 311 # used by parent class bone() to initialize number to 1
574bacfe 312 my $self = shift;
990fb837 313 $IMPORT=1; # call our import only once
9b924220
RGS
314 $self->{_m} = $MBI->_one();
315 $self->{_e} = $MBI->_zero();
316 $self->{_es} = '+';
574bacfe
JH
317 }
318
13a12e00 319sub _bzero
58cde26e 320 {
990fb837 321 # used by parent class bone() to initialize number to 0
58cde26e 322 my $self = shift;
990fb837 323 $IMPORT=1; # call our import only once
9b924220
RGS
324 $self->{_m} = $MBI->_zero();
325 $self->{_e} = $MBI->_one();
326 $self->{_es} = '+';
58cde26e
JH
327 }
328
9393ace2
JH
329sub isa
330 {
331 my ($self,$class) = @_;
56b9c951
JH
332 return if $class =~ /^Math::BigInt/; # we aren't one of these
333 UNIVERSAL::isa($self,$class);
9393ace2
JH
334 }
335
8f675a64
JH
336sub config
337 {
338 # return (later set?) configuration data as hash ref
339 my $class = shift || 'Math::BigFloat';
340
2ebb273f
T
341 if (@_ == 1 && ref($_[0]) ne 'HASH')
342 {
343 my $cfg = $class->SUPER::config();
344 return $cfg->{$_[0]};
345 }
346
990fb837 347 my $cfg = $class->SUPER::config(@_);
8f675a64 348
990fb837 349 # now we need only to override the ones that are different from our parent
8f675a64
JH
350 $cfg->{class} = $class;
351 $cfg->{with} = $MBI;
8f675a64
JH
352 $cfg;
353 }
354
58cde26e 355##############################################################################
9681bfa6 356# string conversion
58cde26e
JH
357
358sub bstr
359 {
360 # (ref to BFLOAT or num_str ) return num_str
361 # Convert number from internal format to (non-scientific) string format.
362 # internal format is always normalized (no leading zeros, "-0" => "+0")
b68b7ab1 363 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
58cde26e 364
574bacfe 365 if ($x->{sign} !~ /^[+-]$/)
58cde26e 366 {
574bacfe
JH
367 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
368 return 'inf'; # +inf
58cde26e 369 }
c38b2de2 370
574bacfe
JH
371 my $es = '0'; my $len = 1; my $cad = 0; my $dot = '.';
372
c38b2de2 373 # $x is zero?
9b924220 374 my $not_zero = !($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
574bacfe 375 if ($not_zero)
58cde26e 376 {
9b924220 377 $es = $MBI->_str($x->{_m});
574bacfe 378 $len = CORE::length($es);
9b924220
RGS
379 my $e = $MBI->_num($x->{_e});
380 $e = -$e if $x->{_es} eq '-';
c38b2de2 381 if ($e < 0)
58cde26e 382 {
c38b2de2
JH
383 $dot = '';
384 # if _e is bigger than a scalar, the following will blow your memory
385 if ($e <= -$len)
574bacfe 386 {
c38b2de2
JH
387 my $r = abs($e) - $len;
388 $es = '0.'. ('0' x $r) . $es; $cad = -($len+$r);
574bacfe
JH
389 }
390 else
391 {
9b924220
RGS
392 substr($es,$e,0) = '.'; $cad = $MBI->_num($x->{_e});
393 $cad = -$cad if $x->{_es} eq '-';
574bacfe 394 }
82cf049f 395 }
c38b2de2
JH
396 elsif ($e > 0)
397 {
398 # expand with zeros
399 $es .= '0' x $e; $len += $e; $cad = 0;
400 }
574bacfe 401 } # if not zero
9b924220 402
c38b2de2
JH
403 $es = '-'.$es if $x->{sign} eq '-';
404 # if set accuracy or precision, pad with zeros on the right side
574bacfe
JH
405 if ((defined $x->{_a}) && ($not_zero))
406 {
407 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
408 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
409 $zeros = $x->{_a} - $len if $cad != $len;
574bacfe 410 $es .= $dot.'0' x $zeros if $zeros > 0;
82cf049f 411 }
c38b2de2 412 elsif ((($x->{_p} || 0) < 0))
58cde26e 413 {
574bacfe
JH
414 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
415 my $zeros = -$x->{_p} + $cad;
574bacfe 416 $es .= $dot.'0' x $zeros if $zeros > 0;
58cde26e 417 }
56b9c951 418 $es;
82cf049f 419 }
f216259d 420
58cde26e
JH
421sub bsstr
422 {
423 # (ref to BFLOAT or num_str ) return num_str
424 # Convert number from internal format to scientific string format.
425 # internal format is always normalized (no leading zeros, "-0E0" => "+0E0")
b68b7ab1 426 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
a0d0e21e 427
574bacfe
JH
428 if ($x->{sign} !~ /^[+-]$/)
429 {
430 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
431 return 'inf'; # +inf
432 }
9b924220 433 my $sep = 'e'.$x->{_es};
56d9de68 434 my $sign = $x->{sign}; $sign = '' if $sign eq '+';
9b924220 435 $sign . $MBI->_str($x->{_m}) . $sep . $MBI->_str($x->{_e});
58cde26e
JH
436 }
437
438sub numify
439 {
fc399313
PJA
440 # Convert a Perl scalar number from a BigFloat object.
441 # Create a string and let Perl's atoi()/atof() handle the rest.
b282a552 442 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
fc399313 443 return 0 + $x->bsstr();
58cde26e 444 }
a0d0e21e 445
58cde26e
JH
446##############################################################################
447# public stuff (usually prefixed with "b")
448
b68b7ab1
T
449sub bneg
450 {
451 # (BINT or num_str) return BINT
452 # negate number or make a negated number from string
453 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
454
455 return $x if $x->modify('bneg');
456
ada8209b 457 # for +0 do not negate (to have always normalized +0). Does nothing for 'NaN'
b68b7ab1
T
458 $x->{sign} =~ tr/+-/-+/ unless ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m}));
459 $x;
460 }
461
574bacfe 462# tels 2001-08-04
b282a552 463# XXX TODO this must be overwritten and return NaN for non-integer values
574bacfe 464# band(), bior(), bxor(), too
58cde26e
JH
465#sub bnot
466# {
467# $class->SUPER::bnot($class,@_);
468# }
469
470sub bcmp
471 {
472 # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
f9a08e12
JH
473
474 # set up parameters
475 my ($self,$x,$y) = (ref($_[0]),@_);
aa45dafa 476
f9a08e12
JH
477 # objectify is costly, so avoid it
478 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
479 {
480 ($self,$x,$y) = objectify(2,@_);
481 }
58cde26e 482
56d9de68
T
483 return $upgrade->bcmp($x,$y) if defined $upgrade &&
484 ((!$x->isa($self)) || (!$y->isa($self)));
485
aa45dafa
PJA
486 # Handle all 'nan' cases.
487
488 return undef if ($x->{sign} eq $nan) || ($y->{sign} eq $nan);
489
490 # Handle all '+inf' and '-inf' cases.
491
492 return 0 if ($x->{sign} eq '+inf' && $y->{sign} eq '+inf' ||
493 $x->{sign} eq '-inf' && $y->{sign} eq '-inf');
494 return +1 if $x->{sign} eq '+inf'; # x = +inf and y < +inf
495 return -1 if $x->{sign} eq '-inf'; # x = -inf and y > -inf
496 return -1 if $y->{sign} eq '+inf'; # x < +inf and y = +inf
497 return +1 if $y->{sign} eq '-inf'; # x > -inf and y = -inf
498
499 # Handle all cases with opposite signs.
0716bf9b 500
aa45dafa
PJA
501 return +1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # also does 0 <=> -y
502 return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # also does -x <=> 0
503
504 # Handle all remaining zero cases.
58cde26e 505
574bacfe
JH
506 my $xz = $x->is_zero();
507 my $yz = $y->is_zero();
aa45dafa
PJA
508 return 0 if $xz && $yz; # 0 <=> 0
509 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
510 return +1 if $yz && $x->{sign} eq '+'; # +x <=> 0
511
512 # Both arguments are now finite, non-zero numbers with the same sign.
513
514 my $cmp;
515
516 # The next step is to compare the exponents, but since each mantissa is an
517 # integer of arbitrary value, the exponents must be normalized by the length
518 # of the mantissas before we can compare them.
519
520 my $mxl = $MBI->_len($x->{_m});
521 my $myl = $MBI->_len($y->{_m});
522
523 # If the mantissas have the same length, there is no point in normalizing the
524 # exponents by the length of the mantissas, so treat that as a special case.
525
526 if ($mxl == $myl) {
527
528 # First handle the two cases where the exponents have different signs.
529
530 if ($x->{_es} eq '+' && $y->{_es} eq '-') {
531 $cmp = +1;
532 }
533
534 elsif ($x->{_es} eq '-' && $y->{_es} eq '+') {
535 $cmp = -1;
536 }
537
538 # Then handle the case where the exponents have the same sign.
539
540 else {
541 $cmp = $MBI->_acmp($x->{_e}, $y->{_e});
542 $cmp = -$cmp if $x->{_es} eq '-';
543 }
544
545 # Adjust for the sign, which is the same for x and y, and bail out if
546 # we're done.
547
548 $cmp = -$cmp if $x->{sign} eq '-'; # 124 > 123, but -124 < -123
549 return $cmp if $cmp;
550
551 }
552
553 # We must normalize each exponent by the length of the corresponding
554 # mantissa. Life is a lot easier if we first make both exponents
555 # non-negative. We do this by adding the same positive value to both
556 # exponent. This is safe, because when comparing the exponents, only the
557 # relative difference is important.
558
559 my $ex;
560 my $ey;
561
562 if ($x->{_es} eq '+') {
563
564 # If the exponent of x is >= 0 and the exponent of y is >= 0, there is no
565 # need to do anything special.
566
567 if ($y->{_es} eq '+') {
568 $ex = $MBI->_copy($x->{_e});
569 $ey = $MBI->_copy($y->{_e});
570 }
571
572 # If the exponent of x is >= 0 and the exponent of y is < 0, add the
573 # absolute value of the exponent of y to both.
574
575 else {
576 $ex = $MBI->_copy($x->{_e});
577 $ex = $MBI->_add($ex, $y->{_e}); # ex + |ey|
578 $ey = $MBI->_zero(); # -ex + |ey| = 0
579 }
580
581 } else {
582
583 # If the exponent of x is < 0 and the exponent of y is >= 0, add the
584 # absolute value of the exponent of x to both.
585
586 if ($y->{_es} eq '+') {
587 $ex = $MBI->_zero(); # -ex + |ex| = 0
588 $ey = $MBI->_copy($y->{_e});
589 $ey = $MBI->_add($ey, $x->{_e}); # ey + |ex|
590 }
591
592 # If the exponent of x is < 0 and the exponent of y is < 0, add the
593 # absolute values of both exponents to both exponents.
594
595 else {
596 $ex = $MBI->_copy($y->{_e}); # -ex + |ey| + |ex| = |ey|
597 $ey = $MBI->_copy($x->{_e}); # -ey + |ex| + |ey| = |ex|
598 }
599
600 }
601
602 # Now we can normalize the exponents by adding lengths of the mantissas.
603
604 $MBI->_add($ex, $MBI->_new($mxl));
605 $MBI->_add($ey, $MBI->_new($myl));
606
607 # We're done if the exponents are different.
608
609 $cmp = $MBI->_acmp($ex, $ey);
610 $cmp = -$cmp if $x->{sign} eq '-'; # 124 > 123, but -124 < -123
611 return $cmp if $cmp;
612
613 # Compare the mantissas, but first normalize them by padding the shorter
614 # mantissa with zeros (shift left) until it has the same length as the longer
615 # mantissa.
616
617 my $mx = $x->{_m};
618 my $my = $y->{_m};
619
620 if ($mxl > $myl) {
621 $my = $MBI->_lsft($MBI->_copy($my), $MBI->_new($mxl - $myl), 10);
622 } elsif ($mxl < $myl) {
623 $mx = $MBI->_lsft($MBI->_copy($mx), $MBI->_new($myl - $mxl), 10);
624 }
625
626 $cmp = $MBI->_acmp($mx, $my);
627 $cmp = -$cmp if $x->{sign} eq '-'; # 124 > 123, but -124 < -123
628 return $cmp;
629
58cde26e
JH
630 }
631
632sub bacmp
633 {
634 # Compares 2 values, ignoring their signs.
635 # Returns one of undef, <0, =0, >0. (suitable for sort)
f9a08e12
JH
636
637 # set up parameters
638 my ($self,$x,$y) = (ref($_[0]),@_);
639 # objectify is costly, so avoid it
640 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
641 {
642 ($self,$x,$y) = objectify(2,@_);
643 }
ee15d750 644
56d9de68
T
645 return $upgrade->bacmp($x,$y) if defined $upgrade &&
646 ((!$x->isa($self)) || (!$y->isa($self)));
647
ee15d750 648 # handle +-inf and NaN's
abcfbf51 649 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/)
ee15d750
JH
650 {
651 return undef if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
652 return 0 if ($x->is_inf() && $y->is_inf());
653 return 1 if ($x->is_inf() && !$y->is_inf());
b3abae2a 654 return -1;
ee15d750
JH
655 }
656
657 # shortcut
658 my $xz = $x->is_zero();
659 my $yz = $y->is_zero();
660 return 0 if $xz && $yz; # 0 <=> 0
661 return -1 if $xz && !$yz; # 0 <=> +y
662 return 1 if $yz && !$xz; # +x <=> 0
663
664 # adjust so that exponents are equal
9b924220
RGS
665 my $lxm = $MBI->_len($x->{_m});
666 my $lym = $MBI->_len($y->{_m});
667 my ($xes,$yes) = (1,1);
668 $xes = -1 if $x->{_es} ne '+';
669 $yes = -1 if $y->{_es} ne '+';
28df3e88 670 # the numify somewhat limits our length, but makes it much faster
9b924220
RGS
671 my $lx = $lxm + $xes * $MBI->_num($x->{_e});
672 my $ly = $lym + $yes * $MBI->_num($y->{_e});
394e6ffb 673 my $l = $lx - $ly;
ee15d750 674 return $l <=> 0 if $l != 0;
58cde26e 675
ee15d750 676 # lengths (corrected by exponent) are equal
394e6ffb 677 # so make mantissa equal-length by padding with zero (shift left)
ee15d750
JH
678 my $diff = $lxm - $lym;
679 my $xm = $x->{_m}; # not yet copy it
680 my $ym = $y->{_m};
681 if ($diff > 0)
682 {
9b924220
RGS
683 $ym = $MBI->_copy($y->{_m});
684 $ym = $MBI->_lsft($ym, $MBI->_new($diff), 10);
ee15d750
JH
685 }
686 elsif ($diff < 0)
687 {
9b924220
RGS
688 $xm = $MBI->_copy($x->{_m});
689 $xm = $MBI->_lsft($xm, $MBI->_new(-$diff), 10);
ee15d750 690 }
9b924220 691 $MBI->_acmp($xm,$ym);
58cde26e 692 }
a0d0e21e 693
58cde26e
JH
694sub badd
695 {
696 # add second arg (BFLOAT or string) to first (BFLOAT) (modifies first)
697 # return result as BFLOAT
f9a08e12
JH
698
699 # set up parameters
80365507 700 my ($self,$x,$y,@r) = (ref($_[0]),@_);
f9a08e12
JH
701 # objectify is costly, so avoid it
702 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
703 {
80365507 704 ($self,$x,$y,@r) = objectify(2,@_);
f9a08e12 705 }
50109ad0
RGS
706
707 return $x if $x->modify('badd');
58cde26e 708
574bacfe
JH
709 # inf and NaN handling
710 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
711 {
712 # NaN first
713 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
13a12e00 714 # inf handling
574bacfe
JH
715 if (($x->{sign} =~ /^[+-]inf$/) && ($y->{sign} =~ /^[+-]inf$/))
716 {
13a12e00
JH
717 # +inf++inf or -inf+-inf => same, rest is NaN
718 return $x if $x->{sign} eq $y->{sign};
719 return $x->bnan();
574bacfe 720 }
56b9c951 721 # +-inf + something => +inf; something +-inf => +-inf
574bacfe
JH
722 $x->{sign} = $y->{sign}, return $x if $y->{sign} =~ /^[+-]inf$/;
723 return $x;
724 }
725
80365507 726 return $upgrade->badd($x,$y,@r) if defined $upgrade &&
8f675a64
JH
727 ((!$x->isa($self)) || (!$y->isa($self)));
728
80365507
T
729 $r[3] = $y; # no push!
730
58cde26e 731 # speed: no add for 0+y or x+0
80365507 732 return $x->bround(@r) if $y->is_zero(); # x+0
58cde26e
JH
733 if ($x->is_zero()) # 0+y
734 {
735 # make copy, clobbering up x (modify in place!)
9b924220
RGS
736 $x->{_e} = $MBI->_copy($y->{_e});
737 $x->{_es} = $y->{_es};
738 $x->{_m} = $MBI->_copy($y->{_m});
58cde26e 739 $x->{sign} = $y->{sign} || $nan;
80365507 740 return $x->round(@r);
a0d0e21e 741 }
58cde26e
JH
742
743 # take lower of the two e's and adapt m1 to it to match m2
28df3e88 744 my $e = $y->{_e};
9b924220
RGS
745 $e = $MBI->_zero() if !defined $e; # if no BFLOAT?
746 $e = $MBI->_copy($e); # make copy (didn't do it yet)
747
748 my $es;
749
750 ($e,$es) = _e_sub($e, $x->{_e}, $y->{_es} || '+', $x->{_es});
751
752 my $add = $MBI->_copy($y->{_m});
753
754 if ($es eq '-') # < 0
58cde26e 755 {
9b924220
RGS
756 $MBI->_lsft( $x->{_m}, $e, 10);
757 ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
58cde26e 758 }
9b924220 759 elsif (!$MBI->_is_zero($e)) # > 0
58cde26e 760 {
9b924220 761 $MBI->_lsft($add, $e, 10);
58cde26e 762 }
61f5c3f5 763 # else: both e are the same, so just leave them
9b924220
RGS
764
765 if ($x->{sign} eq $y->{sign})
766 {
767 # add
768 $x->{_m} = $MBI->_add($x->{_m}, $add);
769 }
770 else
771 {
772 ($x->{_m}, $x->{sign}) =
773 _e_add($x->{_m}, $add, $x->{sign}, $y->{sign});
774 }
775
61f5c3f5 776 # delete trailing zeros, then round
80365507 777 $x->bnorm()->round(@r);
58cde26e
JH
778 }
779
03874afe 780# sub bsub is inherited from Math::BigInt!
58cde26e
JH
781
782sub binc
783 {
784 # increment arg by one
b282a552 785 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
e745a66c 786
50109ad0
RGS
787 return $x if $x->modify('binc');
788
9b924220 789 if ($x->{_es} eq '-')
e745a66c 790 {
b282a552 791 return $x->badd($self->bone(),@r); # digits after dot
e745a66c
JH
792 }
793
9b924220 794 if (!$MBI->_is_zero($x->{_e})) # _e == 0 for NaN, inf, -inf
e745a66c 795 {
b282a552 796 # 1e2 => 100, so after the shift below _m has a '0' as last digit
9b924220
RGS
797 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
798 $x->{_e} = $MBI->_zero(); # normalize
799 $x->{_es} = '+';
b282a552
T
800 # we know that the last digit of $x will be '1' or '9', depending on the
801 # sign
e745a66c
JH
802 }
803 # now $x->{_e} == 0
804 if ($x->{sign} eq '+')
805 {
9b924220 806 $MBI->_inc($x->{_m});
b282a552 807 return $x->bnorm()->bround(@r);
e745a66c
JH
808 }
809 elsif ($x->{sign} eq '-')
810 {
9b924220
RGS
811 $MBI->_dec($x->{_m});
812 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
b282a552 813 return $x->bnorm()->bround(@r);
e745a66c
JH
814 }
815 # inf, nan handling etc
b282a552 816 $x->badd($self->bone(),@r); # badd() does round
58cde26e
JH
817 }
818
819sub bdec
820 {
821 # decrement arg by one
b282a552 822 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
e745a66c 823
50109ad0
RGS
824 return $x if $x->modify('bdec');
825
9b924220 826 if ($x->{_es} eq '-')
e745a66c 827 {
b282a552 828 return $x->badd($self->bone('-'),@r); # digits after dot
e745a66c
JH
829 }
830
9b924220 831 if (!$MBI->_is_zero($x->{_e}))
e745a66c 832 {
9b924220
RGS
833 $x->{_m} = $MBI->_lsft($x->{_m}, $x->{_e},10); # 1e2 => 100
834 $x->{_e} = $MBI->_zero(); # normalize
835 $x->{_es} = '+';
e745a66c
JH
836 }
837 # now $x->{_e} == 0
838 my $zero = $x->is_zero();
839 # <= 0
840 if (($x->{sign} eq '-') || $zero)
841 {
9b924220
RGS
842 $MBI->_inc($x->{_m});
843 $x->{sign} = '-' if $zero; # 0 => 1 => -1
844 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # -1 +1 => -0 => +0
b282a552 845 return $x->bnorm()->round(@r);
e745a66c
JH
846 }
847 # > 0
848 elsif ($x->{sign} eq '+')
849 {
9b924220 850 $MBI->_dec($x->{_m});
b282a552 851 return $x->bnorm()->round(@r);
e745a66c
JH
852 }
853 # inf, nan handling etc
9b924220 854 $x->badd($self->bone('-'),@r); # does round
58cde26e
JH
855 }
856
990fb837
RGS
857sub DEBUG () { 0; }
858
61f5c3f5
T
859sub blog
860 {
990fb837 861 my ($self,$x,$base,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
61f5c3f5 862
50109ad0
RGS
863 return $x if $x->modify('blog');
864
990fb837
RGS
865 # $base > 0, $base != 1; if $base == undef default to $base == e
866 # $x >= 0
9393ace2 867
b3abae2a
JH
868 # we need to limit the accuracy to protect against overflow
869 my $fallback = 0;
990fb837
RGS
870 my ($scale,@params);
871 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
61f5c3f5 872
990fb837
RGS
873 # also takes care of the "error in _find_round_parameters?" case
874 return $x->bnan() if $x->{sign} ne '+' || $x->is_zero();
091c87b1 875
b3abae2a 876 # no rounding at all, so must use fallback
990fb837 877 if (scalar @params == 0)
b3abae2a
JH
878 {
879 # simulate old behaviour
990fb837
RGS
880 $params[0] = $self->div_scale(); # and round to it as accuracy
881 $params[1] = undef; # P = undef
882 $scale = $params[0]+4; # at least four more for proper round
883 $params[2] = $r; # round mode by caller or undef
b3abae2a
JH
884 $fallback = 1; # to clear a/p afterwards
885 }
886 else
887 {
888 # the 4 below is empirical, and there might be cases where it is not
889 # enough...
990fb837 890 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
b3abae2a 891 }
61f5c3f5 892
b3abae2a 893 return $x->bzero(@params) if $x->is_one();
50109ad0 894 # base not defined => base == Euler's number e
990fb837
RGS
895 if (defined $base)
896 {
091c87b1 897 # make object, since we don't feed it through objectify() to still get the
990fb837
RGS
898 # case of $base == undef
899 $base = $self->new($base) unless ref($base);
900 # $base > 0; $base != 1
901 return $x->bnan() if $base->is_zero() || $base->is_one() ||
902 $base->{sign} ne '+';
091c87b1 903 # if $x == $base, we know the result must be 1.0
90d1b129
T
904 if ($x->bcmp($base) == 0)
905 {
906 $x->bone('+',@params);
907 if ($fallback)
908 {
909 # clear a/p after round, since user did not request it
910 delete $x->{_a}; delete $x->{_p};
911 }
912 return $x;
913 }
990fb837 914 }
61f5c3f5 915
b3abae2a 916 # when user set globals, they would interfere with our calculation, so
56d9de68 917 # disable them and later re-enable them
b3abae2a
JH
918 no strict 'refs';
919 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
920 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
921 # we also need to disable any set A or P on $x (_find_round_parameters took
922 # them already into account), since these would interfere, too
923 delete $x->{_a}; delete $x->{_p};
9393ace2 924 # need to disable $upgrade in BigInt, to avoid deep recursion
b3abae2a 925 local $Math::BigInt::upgrade = undef;
93c87d9d 926 local $Math::BigFloat::downgrade = undef;
990fb837
RGS
927
928 # upgrade $x if $x is not a BigFloat (handle BigInt input)
7d193e39 929 # XXX TODO: rebless!
990fb837
RGS
930 if (!$x->isa('Math::BigFloat'))
931 {
932 $x = Math::BigFloat->new($x);
933 $self = ref($x);
934 }
b282a552
T
935
936 my $done = 0;
937
938 # If the base is defined and an integer, try to calculate integer result
939 # first. This is very fast, and in case the real result was found, we can
940 # stop right here.
941 if (defined $base && $base->is_int() && $x->is_int())
942 {
9b924220
RGS
943 my $i = $MBI->_copy( $x->{_m} );
944 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
945 my $int = Math::BigInt->bzero();
946 $int->{value} = $i;
b282a552
T
947 $int->blog($base->as_number());
948 # if ($exact)
9b924220 949 if ($base->as_number()->bpow($int) == $x)
b282a552
T
950 {
951 # found result, return it
9b924220
RGS
952 $x->{_m} = $int->{value};
953 $x->{_e} = $MBI->_zero();
954 $x->{_es} = '+';
b282a552
T
955 $x->bnorm();
956 $done = 1;
957 }
958 }
959
960 if ($done == 0)
9393ace2 961 {
7d193e39
T
962 # base is undef, so base should be e (Euler's number), so first calculate the
963 # log to base e (using reduction by 10 (and probably 2)):
b282a552 964 $self->_log_10($x,$scale);
9b924220 965
b282a552
T
966 # and if a different base was requested, convert it
967 if (defined $base)
968 {
969 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
970 # not ln, but some other base (don't modify $base)
971 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
972 }
9393ace2 973 }
990fb837 974
091c87b1 975 # shortcut to not run through _find_round_parameters again
990fb837 976 if (defined $params[0])
b3abae2a 977 {
990fb837 978 $x->bround($params[0],$params[2]); # then round accordingly
b3abae2a
JH
979 }
980 else
981 {
990fb837 982 $x->bfround($params[1],$params[2]); # then round accordingly
b3abae2a
JH
983 }
984 if ($fallback)
985 {
986 # clear a/p after round, since user did not request it
ef9466ea 987 delete $x->{_a}; delete $x->{_p};
b3abae2a
JH
988 }
989 # restore globals
990 $$abr = $ab; $$pbr = $pb;
991
992 $x;
61f5c3f5
T
993 }
994
50109ad0
RGS
995sub _len_to_steps
996 {
997 # Given D (digits in decimal), compute N so that N! (N factorial) is
998 # at least D digits long. D should be at least 50.
999 my $d = shift;
1000
1001 # two constants for the Ramanujan estimate of ln(N!)
1002 my $lg2 = log(2 * 3.14159265) / 2;
1003 my $lg10 = log(10);
1004
1005 # D = 50 => N => 42, so L = 40 and R = 50
1006 my $l = 40; my $r = $d;
1007
1008 # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
1009 $l = $l->numify if ref($l);
1010 $r = $r->numify if ref($r);
1011 $lg2 = $lg2->numify if ref($lg2);
1012 $lg10 = $lg10->numify if ref($lg10);
1013
1014 # binary search for the right value (could this be written as the reverse of lg(n!)?)
1015 while ($r - $l > 1)
1016 {
1017 my $n = int(($r - $l) / 2) + $l;
1018 my $ramanujan =
1019 int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
1020 $ramanujan > $d ? $r = $n : $l = $n;
1021 }
1022 $l;
1023 }
1024
1025sub bnok
1026 {
1027 # Calculate n over k (binomial coefficient or "choose" function) as integer.
1028 # set up parameters
1029 my ($self,$x,$y,@r) = (ref($_[0]),@_);
1030
1031 # objectify is costly, so avoid it
1032 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1033 {
1034 ($self,$x,$y,@r) = objectify(2,@_);
1035 }
1036
1037 return $x if $x->modify('bnok');
1038
1039 return $x->bnan() if $x->is_nan() || $y->is_nan();
1040 return $x->binf() if $x->is_inf();
1041
1042 my $u = $x->as_int();
1043 $u->bnok($y->as_int());
1044
1045 $x->{_m} = $u->{value};
1046 $x->{_e} = $MBI->_zero();
1047 $x->{_es} = '+';
1048 $x->{sign} = '+';
1049 $x->bnorm(@r);
1050 }
1051
7d193e39
T
1052sub bexp
1053 {
50109ad0 1054 # Calculate e ** X (Euler's number to the power of X)
7d193e39
T
1055 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1056
50109ad0
RGS
1057 return $x if $x->modify('bexp');
1058
7d193e39
T
1059 return $x->binf() if $x->{sign} eq '+inf';
1060 return $x->bzero() if $x->{sign} eq '-inf';
1061
1062 # we need to limit the accuracy to protect against overflow
1063 my $fallback = 0;
1064 my ($scale,@params);
1065 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1066
1067 # also takes care of the "error in _find_round_parameters?" case
1068 return $x if $x->{sign} eq 'NaN';
1069
1070 # no rounding at all, so must use fallback
1071 if (scalar @params == 0)
1072 {
1073 # simulate old behaviour
1074 $params[0] = $self->div_scale(); # and round to it as accuracy
1075 $params[1] = undef; # P = undef
1076 $scale = $params[0]+4; # at least four more for proper round
1077 $params[2] = $r; # round mode by caller or undef
1078 $fallback = 1; # to clear a/p afterwards
1079 }
1080 else
1081 {
1082 # the 4 below is empirical, and there might be cases where it's not enough...
1083 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1084 }
1085
1086 return $x->bone(@params) if $x->is_zero();
1087
1088 if (!$x->isa('Math::BigFloat'))
1089 {
1090 $x = Math::BigFloat->new($x);
1091 $self = ref($x);
1092 }
1093
1094 # when user set globals, they would interfere with our calculation, so
1095 # disable them and later re-enable them
1096 no strict 'refs';
1097 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1098 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1099 # we also need to disable any set A or P on $x (_find_round_parameters took
1100 # them already into account), since these would interfere, too
1101 delete $x->{_a}; delete $x->{_p};
1102 # need to disable $upgrade in BigInt, to avoid deep recursion
1103 local $Math::BigInt::upgrade = undef;
1104 local $Math::BigFloat::downgrade = undef;
1105
1106 my $x_org = $x->copy();
7d193e39
T
1107
1108 # We use the following Taylor series:
1109
1110 # x x^2 x^3 x^4
1111 # e = 1 + --- + --- + --- + --- ...
1112 # 1! 2! 3! 4!
1113
1114 # The difference for each term is X and N, which would result in:
1115 # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term
1116
1117 # But it is faster to compute exp(1) and then raising it to the
1118 # given power, esp. if $x is really big and an integer because:
1119
1120 # * The numerator is always 1, making the computation faster
1121 # * the series converges faster in the case of x == 1
1122 # * We can also easily check when we have reached our limit: when the
1123 # term to be added is smaller than "1E$scale", we can stop - f.i.
1124 # scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
1125 # * we can compute the *exact* result by simulating bigrat math:
1126
1127 # 1 1 gcd(3,4) = 1 1*24 + 1*6 5
1128 # - + - = ---------- = --
1129 # 6 24 6*24 24
1130
1131 # We do not compute the gcd() here, but simple do:
1132 # 1 1 1*24 + 1*6 30
1133 # - + - = --------- = --
1134 # 6 24 6*24 144
1135
1136 # In general:
1137 # a c a*d + c*b and note that c is always 1 and d = (b*f)
1138 # - + - = ---------
1139 # b d b*d
1140
1141 # This leads to: which can be reduced by b to:
1142 # a 1 a*b*f + b a*f + 1
1143 # - + - = --------- = -------
1144 # b b*f b*b*f b*f
1145
1146 # The first terms in the series are:
1147
1148 # 1 1 1 1 1 1 1 1 13700
1149 # -- + -- + -- + -- + -- + --- + --- + ---- = -----
1150 # 1 1 2 6 24 120 720 5040 5040
1151
50109ad0 1152 # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!
7d193e39 1153
50109ad0 1154 if ($scale <= 75)
7d193e39 1155 {
50109ad0
RGS
1156 # set $x directly from a cached string form
1157 $x->{_m} = $MBI->_new(
1158 "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
1159 $x->{sign} = '+';
1160 $x->{_es} = '-';
1161 $x->{_e} = $MBI->_new(79);
7d193e39 1162 }
50109ad0
RGS
1163 else
1164 {
1165 # compute A and B so that e = A / B.
1166
1167 # After some terms we end up with this, so we use it as a starting point:
1168 my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
1169 my $F = $MBI->_new(42); my $step = 42;
1170
1171 # Compute how many steps we need to take to get $A and $B sufficiently big
1172 my $steps = _len_to_steps($scale - 4);
1173# print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
1174 while ($step++ <= $steps)
1175 {
1176 # calculate $a * $f + 1
1177 $A = $MBI->_mul($A, $F);
1178 $A = $MBI->_inc($A);
1179 # increment f
1180 $F = $MBI->_inc($F);
1181 }
1182 # compute $B as factorial of $steps (this is faster than doing it manually)
1183 my $B = $MBI->_fac($MBI->_new($steps));
1184
1185# print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";
7d193e39 1186
50109ad0
RGS
1187 # compute A/B with $scale digits in the result (truncate, not round)
1188 $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
1189 $A = $MBI->_div( $A, $B );
7d193e39 1190
50109ad0
RGS
1191 $x->{_m} = $A;
1192 $x->{sign} = '+';
1193 $x->{_es} = '-';
1194 $x->{_e} = $MBI->_new($scale);
1195 }
7d193e39 1196
50109ad0
RGS
1197 # $x contains now an estimate of e, with some surplus digits, so we can round
1198 if (!$x_org->is_one())
7d193e39 1199 {
50109ad0
RGS
1200 # raise $x to the wanted power and round it in one step:
1201 $x->bpow($x_org, @params);
7d193e39
T
1202 }
1203 else
1204 {
50109ad0
RGS
1205 # else just round the already computed result
1206 delete $x->{_a}; delete $x->{_p};
1207 # shortcut to not run through _find_round_parameters again
1208 if (defined $params[0])
1209 {
1210 $x->bround($params[0],$params[2]); # then round accordingly
1211 }
1212 else
1213 {
1214 $x->bfround($params[1],$params[2]); # then round accordingly
1215 }
7d193e39
T
1216 }
1217 if ($fallback)
1218 {
1219 # clear a/p after round, since user did not request it
1220 delete $x->{_a}; delete $x->{_p};
1221 }
1222 # restore globals
1223 $$abr = $ab; $$pbr = $pb;
1224
1225 $x; # return modified $x
1226 }
1227
990fb837
RGS
1228sub _log
1229 {
091c87b1 1230 # internal log function to calculate ln() based on Taylor series.
990fb837
RGS
1231 # Modifies $x in place.
1232 my ($self,$x,$scale) = @_;
1233
091c87b1
T
1234 # in case of $x == 1, result is 0
1235 return $x->bzero() if $x->is_one();
1236
9681bfa6 1237 # XXX TODO: rewrite this in a similar manner to bexp()
7d193e39 1238
990fb837
RGS
1239 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
1240
1241 # u = x-1, v = x+1
1242 # _ _
1243 # Taylor: | u 1 u^3 1 u^5 |
1244 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
1245 # |_ v 3 v^3 5 v^5 _|
1246
1247 # This takes much more steps to calculate the result and is thus not used
1248 # u = x-1
1249 # _ _
1250 # Taylor: | u 1 u^2 1 u^3 |
1251 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
1252 # |_ x 2 x^2 3 x^3 _|
1253
990fb837
RGS
1254 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
1255
1256 $v = $x->copy(); $v->binc(); # v = x+1
1257 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
1258 $x->bdiv($v,$scale); # first term: u/v
1259 $below = $v->copy();
1260 $over = $u->copy();
1261 $u *= $u; $v *= $v; # u^2, v^2
1262 $below->bmul($v); # u^3, v^3
1263 $over->bmul($u);
1264 $factor = $self->new(3); $f = $self->new(2);
1265
1266 my $steps = 0 if DEBUG;
1267 $limit = $self->new("1E-". ($scale-1));
1268 while (3 < 5)
1269 {
1270 # we calculate the next term, and add it to the last
1271 # when the next term is below our limit, it won't affect the outcome
1272 # anymore, so we stop
1273
1274 # calculating the next term simple from over/below will result in quite
1275 # a time hog if the input has many digits, since over and below will
1276 # accumulate more and more digits, and the result will also have many
1277 # digits, but in the end it is rounded to $scale digits anyway. So if we
1278 # round $over and $below first, we save a lot of time for the division
1279 # (not with log(1.2345), but try log (123**123) to see what I mean. This
1280 # can introduce a rounding error if the division result would be f.i.
1281 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
091c87b1
T
1282 # if we truncated $over and $below we might get 0.12345. Does this matter
1283 # for the end result? So we give $over and $below 4 more digits to be
1284 # on the safe side (unscientific error handling as usual... :+D
7d193e39 1285
990fb837
RGS
1286 $next = $over->copy->bround($scale+4)->bdiv(
1287 $below->copy->bmul($factor)->bround($scale+4),
1288 $scale);
1289
1290## old version:
1291## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
1292
1293 last if $next->bacmp($limit) <= 0;
1294
1295 delete $next->{_a}; delete $next->{_p};
1296 $x->badd($next);
990fb837
RGS
1297 # calculate things for the next term
1298 $over *= $u; $below *= $v; $factor->badd($f);
1299 if (DEBUG)
1300 {
1301 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
1302 }
1303 }
990fb837 1304 print "took $steps steps\n" if DEBUG;
7d193e39 1305 $x->bmul($f); # $x *= 2
990fb837
RGS
1306 }
1307
1308sub _log_10
1309 {
091c87b1
T
1310 # Internal log function based on reducing input to the range of 0.1 .. 9.99
1311 # and then "correcting" the result to the proper one. Modifies $x in place.
990fb837
RGS
1312 my ($self,$x,$scale) = @_;
1313
7d193e39 1314 # Taking blog() from numbers greater than 10 takes a *very long* time, so we
990fb837 1315 # break the computation down into parts based on the observation that:
7d193e39
T
1316 # blog(X*Y) = blog(X) + blog(Y)
1317 # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
1318 # $x is the faster it gets. Since 2*$x takes about 10 times as
1319 # long, we make it faster by about a factor of 100 by dividing $x by 10.
1320
1321 # The same observation is valid for numbers smaller than 0.1, e.g. computing
1322 # log(1) is fastest, and the further away we get from 1, the longer it takes.
1323 # So we also 'break' this down by multiplying $x with 10 and subtract the
990fb837
RGS
1324 # log(10) afterwards to get the correct result.
1325
7d193e39
T
1326 # To get $x even closer to 1, we also divide by 2 and then use log(2) to
1327 # correct for this. For instance if $x is 2.4, we use the formula:
1328 # blog(2.4 * 2) == blog (1.2) + blog(2)
1329 # and thus calculate only blog(1.2) and blog(2), which is faster in total
1330 # than calculating blog(2.4).
1331
1332 # In addition, the values for blog(2) and blog(10) are cached.
1333
1334 # Calculate nr of digits before dot:
9b924220
RGS
1335 my $dbd = $MBI->_num($x->{_e});
1336 $dbd = -$dbd if $x->{_es} eq '-';
1337 $dbd += $MBI->_len($x->{_m});
990fb837
RGS
1338
1339 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
1340 # infinite recursion
1341
1342 my $calc = 1; # do some calculation?
1343
1344 # disable the shortcut for 10, since we need log(10) and this would recurse
1345 # infinitely deep
9b924220 1346 if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
990fb837
RGS
1347 {
1348 $dbd = 0; # disable shortcut
1349 # we can use the cached value in these cases
1350 if ($scale <= $LOG_10_A)
1351 {
7d193e39 1352 $x->bzero(); $x->badd($LOG_10); # modify $x in place
990fb837
RGS
1353 $calc = 0; # no need to calc, but round
1354 }
7d193e39 1355 # if we can't use the shortcut, we continue normally
990fb837 1356 }
091c87b1 1357 else
990fb837 1358 {
091c87b1 1359 # disable the shortcut for 2, since we maybe have it cached
9b924220 1360 if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
990fb837 1361 {
091c87b1
T
1362 $dbd = 0; # disable shortcut
1363 # we can use the cached value in these cases
1364 if ($scale <= $LOG_2_A)
1365 {
7d193e39 1366 $x->bzero(); $x->badd($LOG_2); # modify $x in place
091c87b1
T
1367 $calc = 0; # no need to calc, but round
1368 }
7d193e39 1369 # if we can't use the shortcut, we continue normally
990fb837
RGS
1370 }
1371 }
1372
1373 # if $x = 0.1, we know the result must be 0-log(10)
9b924220
RGS
1374 if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
1375 $MBI->_is_one($x->{_m}))
990fb837
RGS
1376 {
1377 $dbd = 0; # disable shortcut
1378 # we can use the cached value in these cases
1379 if ($scale <= $LOG_10_A)
1380 {
1381 $x->bzero(); $x->bsub($LOG_10);
1382 $calc = 0; # no need to calc, but round
1383 }
1384 }
1385
091c87b1
T
1386 return if $calc == 0; # already have the result
1387
990fb837
RGS
1388 # default: these correction factors are undef and thus not used
1389 my $l_10; # value of ln(10) to A of $scale
1390 my $l_2; # value of ln(2) to A of $scale
1391
7d193e39
T
1392 my $two = $self->new(2);
1393
990fb837
RGS
1394 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
1395 # so don't do this shortcut for 1 or 0
1396 if (($dbd > 1) || ($dbd < 0))
1397 {
1398 # convert our cached value to an object if not already (avoid doing this
1399 # at import() time, since not everybody needs this)
1400 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
1401
1402 #print "x = $x, dbd = $dbd, calc = $calc\n";
1403 # got more than one digit before the dot, or more than one zero after the
1404 # dot, so do:
1405 # log(123) == log(1.23) + log(10) * 2
1406 # log(0.0123) == log(1.23) - log(10) * 2
1407
1408 if ($scale <= $LOG_10_A)
1409 {
1410 # use cached value
990fb837
RGS
1411 $l_10 = $LOG_10->copy(); # copy for mul
1412 }
1413 else
1414 {
7d193e39 1415 # else: slower, compute and cache result
990fb837
RGS
1416 # also disable downgrade for this code path
1417 local $Math::BigFloat::downgrade = undef;
7d193e39
T
1418
1419 # shorten the time to calculate log(10) based on the following:
1420 # log(1.25 * 8) = log(1.25) + log(8)
1421 # = log(1.25) + log(2) + log(2) + log(2)
1422
1423 # first get $l_2 (and possible compute and cache log(2))
1424 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1425 if ($scale <= $LOG_2_A)
1426 {
1427 # use cached value
1428 $l_2 = $LOG_2->copy(); # copy() for the mul below
1429 }
1430 else
1431 {
1432 # else: slower, compute and cache result
1433 $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1434 $LOG_2 = $l_2->copy(); # cache the result for later
1435 # the copy() is for mul below
1436 $LOG_2_A = $scale;
1437 }
1438
1439 # now calculate log(1.25):
1440 $l_10 = $self->new('1.25'); $self->_log($l_10, $scale); # scale+4, actually
1441
1442 # log(1.25) + log(2) + log(2) + log(2):
1443 $l_10->badd($l_2);
1444 $l_10->badd($l_2);
1445 $l_10->badd($l_2);
1446 $LOG_10 = $l_10->copy(); # cache the result for later
1447 # the copy() is for mul below
1448 $LOG_10_A = $scale;
990fb837
RGS
1449 }
1450 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
9b924220
RGS
1451 $l_10->bmul( $self->new($dbd)); # log(10) * (digits_before_dot-1)
1452 my $dbd_sign = '+';
1453 if ($dbd < 0)
1454 {
1455 $dbd = -$dbd;
1456 $dbd_sign = '-';
1457 }
1458 ($x->{_e}, $x->{_es}) =
1459 _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
990fb837
RGS
1460
1461 }
1462
1463 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1464
1465 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1466 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1467
27e7b8bb
T
1468 $HALF = $self->new($HALF) unless ref($HALF);
1469
091c87b1 1470 my $twos = 0; # default: none (0 times)
7d193e39 1471 while ($x->bacmp($HALF) <= 0) # X <= 0.5
990fb837 1472 {
091c87b1
T
1473 $twos--; $x->bmul($two);
1474 }
7d193e39 1475 while ($x->bacmp($two) >= 0) # X >= 2
091c87b1
T
1476 {
1477 $twos++; $x->bdiv($two,$scale+4); # keep all digits
1478 }
7d193e39
T
1479 # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
1480 # So calculate correction factor based on ln(2):
091c87b1
T
1481 if ($twos != 0)
1482 {
1483 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1484 if ($scale <= $LOG_2_A)
990fb837 1485 {
091c87b1 1486 # use cached value
7d193e39 1487 $l_2 = $LOG_2->copy(); # copy() for the mul below
990fb837 1488 }
091c87b1 1489 else
990fb837 1490 {
7d193e39 1491 # else: slower, compute and cache result
091c87b1
T
1492 # also disable downgrade for this code path
1493 local $Math::BigFloat::downgrade = undef;
7d193e39
T
1494 $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1495 $LOG_2 = $l_2->copy(); # cache the result for later
1496 # the copy() is for mul below
1497 $LOG_2_A = $scale;
990fb837 1498 }
091c87b1 1499 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
990fb837
RGS
1500 }
1501
091c87b1
T
1502 $self->_log($x,$scale); # need to do the "normal" way
1503 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
1504 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
7d193e39 1505
990fb837 1506 # all done, $x contains now the result
7d193e39 1507 $x;
990fb837
RGS
1508 }
1509
58cde26e
JH
1510sub blcm
1511 {
ee15d750 1512 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
58cde26e
JH
1513 # does not modify arguments, but returns new object
1514 # Lowest Common Multiplicator
58cde26e
JH
1515
1516 my ($self,@arg) = objectify(0,@_);
1517 my $x = $self->new(shift @arg);
b68b7ab1 1518 while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); }
58cde26e
JH
1519 $x;
1520 }
1521
b68b7ab1
T
1522sub bgcd
1523 {
1524 # (BINT or num_str, BINT or num_str) return BINT
58cde26e 1525 # does not modify arguments, but returns new object
b68b7ab1
T
1526
1527 my $y = shift;
1528 $y = __PACKAGE__->new($y) if !ref($y);
1529 my $self = ref($y);
1530 my $x = $y->copy()->babs(); # keep arguments
1531
1532 return $x->bnan() if $x->{sign} !~ /^[+-]$/ # x NaN?
1533 || !$x->is_int(); # only for integers now
1534
1535 while (@_)
1536 {
1537 my $t = shift; $t = $self->new($t) if !ref($t);
1538 $y = $t->copy()->babs();
1539
1540 return $x->bnan() if $y->{sign} !~ /^[+-]$/ # y NaN?
1541 || !$y->is_int(); # only for integers now
1542
1543 # greatest common divisor
1544 while (! $y->is_zero())
1545 {
1546 ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
1547 }
1548
1549 last if $x->is_one();
1550 }
58cde26e
JH
1551 $x;
1552 }
1553
9b924220 1554##############################################################################
b3abae2a 1555
9b924220 1556sub _e_add
091c87b1 1557 {
9b924220
RGS
1558 # Internal helper sub to take two positive integers and their signs and
1559 # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1560 # output ($CALC,('+'|'-'))
1561 my ($x,$y,$xs,$ys) = @_;
1562
1563 # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1564 if ($xs eq $ys)
1565 {
1566 $x = $MBI->_add ($x, $y ); # a+b
1567 # the sign follows $xs
1568 return ($x, $xs);
1569 }
091c87b1 1570
9b924220
RGS
1571 my $a = $MBI->_acmp($x,$y);
1572 if ($a > 0)
1573 {
1574 $x = $MBI->_sub ($x , $y); # abs sub
1575 }
1576 elsif ($a == 0)
1577 {
1578 $x = $MBI->_zero(); # result is 0
1579 $xs = '+';
1580 }
1581 else # a < 0
1582 {
1583 $x = $MBI->_sub ( $y, $x, 1 ); # abs sub
1584 $xs = $ys;
1585 }
1586 ($x,$xs);
091c87b1
T
1587 }
1588
9b924220
RGS
1589sub _e_sub
1590 {
1591 # Internal helper sub to take two positive integers and their signs and
1592 # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1593 # output ($CALC,('+'|'-'))
1594 my ($x,$y,$xs,$ys) = @_;
1595
1596 # flip sign
1597 $ys =~ tr/+-/-+/;
1598 _e_add($x,$y,$xs,$ys); # call add (does subtract now)
1599 }
1600
1601###############################################################################
1602# is_foo methods (is_negative, is_positive are inherited from BigInt)
1603
b3abae2a
JH
1604sub is_int
1605 {
1606 # return true if arg (BFLOAT or num_str) is an integer
091c87b1 1607 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
b3abae2a 1608
80365507
T
1609 (($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1610 ($x->{_es} eq '+')) ? 1 : 0; # 1e-1 => no integer
b3abae2a
JH
1611 }
1612
58cde26e
JH
1613sub is_zero
1614 {
b3abae2a 1615 # return true if arg (BFLOAT or num_str) is zero
091c87b1 1616 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
574bacfe 1617
80365507 1618 ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m})) ? 1 : 0;
58cde26e
JH
1619 }
1620
1621sub is_one
1622 {
b3abae2a 1623 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
091c87b1 1624 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
ee15d750 1625
990fb837 1626 $sign = '+' if !defined $sign || $sign ne '-';
80365507
T
1627
1628 ($x->{sign} eq $sign &&
1629 $MBI->_is_zero($x->{_e}) &&
1630 $MBI->_is_one($x->{_m}) ) ? 1 : 0;
58cde26e
JH
1631 }
1632
1633sub is_odd
1634 {
ee15d750 1635 # return true if arg (BFLOAT or num_str) is odd or false if even
091c87b1 1636 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
0716bf9b 1637
80365507
T
1638 (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1639 ($MBI->_is_zero($x->{_e})) &&
1640 ($MBI->_is_odd($x->{_m}))) ? 1 : 0;
58cde26e
JH
1641 }
1642
1643sub is_even
1644 {
b22b3e31 1645 # return true if arg (BINT or num_str) is even or false if odd
091c87b1 1646 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
0716bf9b 1647
80365507
T
1648 (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1649 ($x->{_es} eq '+') && # 123.45 isn't
1650 ($MBI->_is_even($x->{_m}))) ? 1 : 0; # but 1200 is
58cde26e
JH
1651 }
1652
80365507 1653sub bmul
58cde26e 1654 {
80365507 1655 # multiply two numbers
f9a08e12
JH
1656
1657 # set up parameters
80365507 1658 my ($self,$x,$y,@r) = (ref($_[0]),@_);
f9a08e12
JH
1659 # objectify is costly, so avoid it
1660 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1661 {
80365507 1662 ($self,$x,$y,@r) = objectify(2,@_);
f9a08e12 1663 }
58cde26e 1664
50109ad0
RGS
1665 return $x if $x->modify('bmul');
1666
58cde26e
JH
1667 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1668
574bacfe
JH
1669 # inf handling
1670 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1671 {
13a12e00 1672 return $x->bnan() if $x->is_zero() || $y->is_zero();
574bacfe
JH
1673 # result will always be +-inf:
1674 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1675 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1676 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1677 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1678 return $x->binf('-');
1679 }
8f675a64 1680
80365507
T
1681 return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1682 ((!$x->isa($self)) || (!$y->isa($self)));
1683
1684 # aEb * cEd = (a*c)E(b+d)
1685 $MBI->_mul($x->{_m},$y->{_m});
1686 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1687
1688 $r[3] = $y; # no push!
1689
1690 # adjust sign:
1691 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1692 $x->bnorm->round(@r);
1693 }
1694
1695sub bmuladd
1696 {
1697 # multiply two numbers and add the third to the result
1698
1699 # set up parameters
913a64d5 1700 my ($self,$x,$y,$z,@r) = objectify(3,@_);
80365507
T
1701
1702 return $x if $x->modify('bmuladd');
1703
1704 return $x->bnan() if (($x->{sign} eq $nan) ||
1705 ($y->{sign} eq $nan) ||
1706 ($z->{sign} eq $nan));
1707
1708 # inf handling
1709 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1710 {
1711 return $x->bnan() if $x->is_zero() || $y->is_zero();
1712 # result will always be +-inf:
1713 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1714 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1715 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1716 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1717 return $x->binf('-');
1718 }
1719
1720 return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
8f675a64 1721 ((!$x->isa($self)) || (!$y->isa($self)));
574bacfe 1722
58cde26e 1723 # aEb * cEd = (a*c)E(b+d)
9b924220
RGS
1724 $MBI->_mul($x->{_m},$y->{_m});
1725 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1726
80365507
T
1727 $r[3] = $y; # no push!
1728
58cde26e
JH
1729 # adjust sign:
1730 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
80365507
T
1731
1732 # z=inf handling (z=NaN handled above)
1733 $x->{sign} = $z->{sign}, return $x if $z->{sign} =~ /^[+-]inf$/;
1734
1735 # take lower of the two e's and adapt m1 to it to match m2
1736 my $e = $z->{_e};
1737 $e = $MBI->_zero() if !defined $e; # if no BFLOAT?
1738 $e = $MBI->_copy($e); # make copy (didn't do it yet)
1739
1740 my $es;
1741
1742 ($e,$es) = _e_sub($e, $x->{_e}, $z->{_es} || '+', $x->{_es});
1743
1744 my $add = $MBI->_copy($z->{_m});
1745
1746 if ($es eq '-') # < 0
1747 {
1748 $MBI->_lsft( $x->{_m}, $e, 10);
1749 ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
1750 }
1751 elsif (!$MBI->_is_zero($e)) # > 0
1752 {
1753 $MBI->_lsft($add, $e, 10);
1754 }
1755 # else: both e are the same, so just leave them
1756
1757 if ($x->{sign} eq $z->{sign})
1758 {
1759 # add
1760 $x->{_m} = $MBI->_add($x->{_m}, $add);
1761 }
1762 else
1763 {
1764 ($x->{_m}, $x->{sign}) =
1765 _e_add($x->{_m}, $add, $x->{sign}, $z->{sign});
1766 }
1767
1768 # delete trailing zeros, then round
1769 $x->bnorm()->round(@r);
58cde26e
JH
1770 }
1771
1772sub bdiv
1773 {
1774 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
9393ace2 1775 # (BFLOAT,BFLOAT) (quo,rem) or BFLOAT (only rem)
f9a08e12
JH
1776
1777 # set up parameters
1778 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1779 # objectify is costly, so avoid it
1780 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1781 {
1782 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1783 }
58cde26e 1784
50109ad0
RGS
1785 return $x if $x->modify('bdiv');
1786
13a12e00
JH
1787 return $self->_div_inf($x,$y)
1788 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/) || $y->is_zero());
574bacfe 1789
13a12e00 1790 # x== 0 # also: or y == 1 or y == -1
394e6ffb 1791 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
0716bf9b 1792
9393ace2
JH
1793 # upgrade ?
1794 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
13a12e00 1795
58cde26e 1796 # we need to limit the accuracy to protect against overflow
574bacfe 1797 my $fallback = 0;
990fb837
RGS
1798 my (@params,$scale);
1799 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1800
1801 return $x if $x->is_nan(); # error in _find_round_parameters?
ee15d750
JH
1802
1803 # no rounding at all, so must use fallback
990fb837 1804 if (scalar @params == 0)
58cde26e 1805 {
0716bf9b 1806 # simulate old behaviour
990fb837
RGS
1807 $params[0] = $self->div_scale(); # and round to it as accuracy
1808 $scale = $params[0]+4; # at least four more for proper round
1809 $params[2] = $r; # round mode by caller or undef
ee15d750
JH
1810 $fallback = 1; # to clear a/p afterwards
1811 }
1812 else
1813 {
1814 # the 4 below is empirical, and there might be cases where it is not
1815 # enough...
990fb837 1816 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
a0d0e21e 1817 }
03874afe
T
1818
1819 my $rem; $rem = $self->bzero() if wantarray;
1820
1821 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1822
9b924220 1823 my $lx = $MBI->_len($x->{_m}); my $ly = $MBI->_len($y->{_m});
58cde26e 1824 $scale = $lx if $lx > $scale;
58cde26e 1825 $scale = $ly if $ly > $scale;
0716bf9b
JH
1826 my $diff = $ly - $lx;
1827 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
a87115f0
RGS
1828
1829 # already handled inf/NaN/-inf above:
1830
233f7bc0
T
1831 # check that $y is not 1 nor -1 and cache the result:
1832 my $y_not_one = !($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m}));
1833
7596a890
RGS
1834 # flipping the sign of $y will also flip the sign of $x for the special
1835 # case of $x->bsub($x); so we can catch it below:
1836 my $xsign = $x->{sign};
1837 $y->{sign} =~ tr/+-/-+/;
1838
a87115f0 1839 if ($xsign ne $x->{sign})
b3abae2a 1840 {
a87115f0 1841 # special case of $x /= $x results in 1
233f7bc0 1842 $x->bone(); # "fixes" also sign of $y, since $x is $y
b3abae2a 1843 }
03874afe 1844 else
58cde26e 1845 {
a87115f0
RGS
1846 # correct $y's sign again
1847 $y->{sign} =~ tr/+-/-+/;
1848 # continue with normal div code:
1849
f603091d 1850 # make copy of $x in case of list context for later remainder calculation
a87115f0 1851 if (wantarray && $y_not_one)
03874afe
T
1852 {
1853 $rem = $x->copy();
1854 }
394e6ffb 1855
03874afe 1856 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
9b924220 1857
03874afe 1858 # check for / +-1 ( +/- 1E0)
a87115f0 1859 if ($y_not_one)
03874afe
T
1860 {
1861 # promote BigInts and it's subclasses (except when already a BigFloat)
1862 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1863
1864 # calculate the result to $scale digits and then round it
1865 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1866 $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1867 $MBI->_div ($x->{_m},$y->{_m}); # a/c
1868
1869 # correct exponent of $x
1870 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1871 # correct for 10**scale
1872 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1873 $x->bnorm(); # remove trailing 0's
1874 }
ada8209b 1875 } # end else $x != $y
a5f75d66 1876
091c87b1 1877 # shortcut to not run through _find_round_parameters again
990fb837 1878 if (defined $params[0])
ee15d750 1879 {
ef9466ea 1880 delete $x->{_a}; # clear before round
990fb837 1881 $x->bround($params[0],$params[2]); # then round accordingly
ee15d750
JH
1882 }
1883 else
1884 {
ef9466ea 1885 delete $x->{_p}; # clear before round
990fb837 1886 $x->bfround($params[1],$params[2]); # then round accordingly
ee15d750 1887 }
574bacfe
JH
1888 if ($fallback)
1889 {
1890 # clear a/p after round, since user did not request it
ef9466ea 1891 delete $x->{_a}; delete $x->{_p};
574bacfe 1892 }
03874afe 1893
58cde26e
JH
1894 if (wantarray)
1895 {
a87115f0 1896 if ($y_not_one)
394e6ffb 1897 {
990fb837 1898 $rem->bmod($y,@params); # copy already done
394e6ffb 1899 }
574bacfe
JH
1900 if ($fallback)
1901 {
1902 # clear a/p after round, since user did not request it
ef9466ea 1903 delete $rem->{_a}; delete $rem->{_p};
574bacfe 1904 }
0716bf9b 1905 return ($x,$rem);
58cde26e 1906 }
9393ace2 1907 $x;
58cde26e 1908 }
a0d0e21e 1909
58cde26e
JH
1910sub bmod
1911 {
f603091d 1912 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return remainder
f9a08e12
JH
1913
1914 # set up parameters
1915 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1916 # objectify is costly, so avoid it
1917 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1918 {
1919 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1920 }
a0d0e21e 1921
50109ad0
RGS
1922 return $x if $x->modify('bmod');
1923
9b924220 1924 # handle NaN, inf, -inf
61f5c3f5
T
1925 if (($x->{sign} !~ /^[+-]$/) || ($y->{sign} !~ /^[+-]$/))
1926 {
1927 my ($d,$re) = $self->SUPER::_div_inf($x,$y);
f9a08e12
JH
1928 $x->{sign} = $re->{sign};
1929 $x->{_e} = $re->{_e};
1930 $x->{_m} = $re->{_m};
1931 return $x->round($a,$p,$r,$y);
61f5c3f5 1932 }
9b924220
RGS
1933 if ($y->is_zero())
1934 {
1935 return $x->bnan() if $x->is_zero();
1936 return $x;
1937 }
233f7bc0 1938
7596a890
RGS
1939 return $x->bzero() if $x->is_zero()
1940 || ($x->is_int() &&
1941 # check that $y == +1 or $y == -1:
1942 ($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m})));
58cde26e 1943
61f5c3f5
T
1944 my $cmp = $x->bacmp($y); # equal or $x < $y?
1945 return $x->bzero($a,$p) if $cmp == 0; # $x == $y => result 0
1946
1947 # only $y of the operands negative?
1948 my $neg = 0; $neg = 1 if $x->{sign} ne $y->{sign};
1949
1950 $x->{sign} = $y->{sign}; # calc sign first
1951 return $x->round($a,$p,$r) if $cmp < 0 && $neg == 0; # $x < $y => result $x
1952
9b924220 1953 my $ym = $MBI->_copy($y->{_m});
61f5c3f5
T
1954
1955 # 2e1 => 20
9b924220
RGS
1956 $MBI->_lsft( $ym, $y->{_e}, 10)
1957 if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
61f5c3f5
T
1958
1959 # if $y has digits after dot
1960 my $shifty = 0; # correct _e of $x by this
9b924220 1961 if ($y->{_es} eq '-') # has digits after dot
61f5c3f5
T
1962 {
1963 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
9b924220
RGS
1964 $shifty = $MBI->_num($y->{_e}); # no more digits after dot
1965 $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
61f5c3f5
T
1966 }
1967 # $ym is now mantissa of $y based on exponent 0
b3abae2a 1968
61f5c3f5 1969 my $shiftx = 0; # correct _e of $x by this
9b924220 1970 if ($x->{_es} eq '-') # has digits after dot
61f5c3f5
T
1971 {
1972 # 123.4 % 20 => 1234 % 200
9b924220
RGS
1973 $shiftx = $MBI->_num($x->{_e}); # no more digits after dot
1974 $MBI->_lsft($ym, $x->{_e}, 10); # 123 => 1230
61f5c3f5
T
1975 }
1976 # 123e1 % 20 => 1230 % 20
9b924220 1977 if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
61f5c3f5 1978 {
9b924220 1979 $MBI->_lsft( $x->{_m}, $x->{_e},10); # es => '+' here
61f5c3f5 1980 }
9b924220
RGS
1981
1982 $x->{_e} = $MBI->_new($shiftx);
1983 $x->{_es} = '+';
1984 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
1985 $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
61f5c3f5
T
1986
1987 # now mantissas are equalized, exponent of $x is adjusted, so calc result
b3abae2a 1988
9b924220 1989 $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
61f5c3f5 1990
9b924220 1991 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
61f5c3f5
T
1992 $x->bnorm();
1993
1994 if ($neg != 0) # one of them negative => correct in place
1995 {
1996 my $r = $y - $x;
1997 $x->{_m} = $r->{_m};
1998 $x->{_e} = $r->{_e};
9b924220
RGS
1999 $x->{_es} = $r->{_es};
2000 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
61f5c3f5
T
2001 $x->bnorm();
2002 }
2003
2004 $x->round($a,$p,$r,$y); # round and return
58cde26e
JH
2005 }
2006
990fb837
RGS
2007sub broot
2008 {
2009 # calculate $y'th root of $x
3a427a11
RGS
2010
2011 # set up parameters
2012 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2013 # objectify is costly, so avoid it
2014 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2015 {
2016 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2017 }
990fb837 2018
50109ad0
RGS
2019 return $x if $x->modify('broot');
2020
990fb837
RGS
2021 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
2022 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
2023 $y->{sign} !~ /^\+$/;
2024
2025 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
2026
2027 # we need to limit the accuracy to protect against overflow
2028 my $fallback = 0;
2029 my (@params,$scale);
2030 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2031
2032 return $x if $x->is_nan(); # error in _find_round_parameters?
2033
2034 # no rounding at all, so must use fallback
2035 if (scalar @params == 0)
2036 {
2037 # simulate old behaviour
2038 $params[0] = $self->div_scale(); # and round to it as accuracy
2039 $scale = $params[0]+4; # at least four more for proper round
ada8209b 2040 $params[2] = $r; # round mode by caller or undef
990fb837
RGS
2041 $fallback = 1; # to clear a/p afterwards
2042 }
2043 else
2044 {
2045 # the 4 below is empirical, and there might be cases where it is not
2046 # enough...
2047 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2048 }
2049
2050 # when user set globals, they would interfere with our calculation, so
2051 # disable them and later re-enable them
2052 no strict 'refs';
2053 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2054 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2055 # we also need to disable any set A or P on $x (_find_round_parameters took
2056 # them already into account), since these would interfere, too
2057 delete $x->{_a}; delete $x->{_p};
2058 # need to disable $upgrade in BigInt, to avoid deep recursion
2059 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
2060
2061 # remember sign and make $x positive, since -4 ** (1/2) => -2
27e7b8bb
T
2062 my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
2063
2064 my $is_two = 0;
2065 if ($y->isa('Math::BigFloat'))
2066 {
2067 $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
2068 }
2069 else
2070 {
2071 $is_two = ($y == 2);
2072 }
990fb837 2073
27e7b8bb
T
2074 # normal square root if $y == 2:
2075 if ($is_two)
990fb837
RGS
2076 {
2077 $x->bsqrt($scale+4);
2078 }
2079 elsif ($y->is_one('-'))
2080 {
2081 # $x ** -1 => 1/$x
2082 my $u = $self->bone()->bdiv($x,$scale);
2083 # copy private parts over
2084 $x->{_m} = $u->{_m};
2085 $x->{_e} = $u->{_e};
9b924220 2086 $x->{_es} = $u->{_es};
990fb837
RGS
2087 }
2088 else
2089 {
3a427a11
RGS
2090 # calculate the broot() as integer result first, and if it fits, return
2091 # it rightaway (but only if $x and $y are integer):
2092
2093 my $done = 0; # not yet
2094 if ($y->is_int() && $x->is_int())
2095 {
9b924220
RGS
2096 my $i = $MBI->_copy( $x->{_m} );
2097 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2098 my $int = Math::BigInt->bzero();
2099 $int->{value} = $i;
3a427a11
RGS
2100 $int->broot($y->as_number());
2101 # if ($exact)
2102 if ($int->copy()->bpow($y) == $x)
2103 {
2104 # found result, return it
9b924220
RGS
2105 $x->{_m} = $int->{value};
2106 $x->{_e} = $MBI->_zero();
2107 $x->{_es} = '+';
3a427a11
RGS
2108 $x->bnorm();
2109 $done = 1;
2110 }
2111 }
2112 if ($done == 0)
2113 {
2114 my $u = $self->bone()->bdiv($y,$scale+4);
2115 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
2116 $x->bpow($u,$scale+4); # el cheapo
2117 }
990fb837
RGS
2118 }
2119 $x->bneg() if $sign == 1;
2120
091c87b1 2121 # shortcut to not run through _find_round_parameters again
990fb837
RGS
2122 if (defined $params[0])
2123 {
2124 $x->bround($params[0],$params[2]); # then round accordingly
2125 }
2126 else
2127 {
2128 $x->bfround($params[1],$params[2]); # then round accordingly
2129 }
2130 if ($fallback)
2131 {
2132 # clear a/p after round, since user did not request it
ef9466ea 2133 delete $x->{_a}; delete $x->{_p};
990fb837
RGS
2134 }
2135 # restore globals
2136 $$abr = $ab; $$pbr = $pb;
2137 $x;
2138 }
2139
58cde26e
JH
2140sub bsqrt
2141 {
990fb837 2142 # calculate square root
ee15d750 2143 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 2144
50109ad0
RGS
2145 return $x if $x->modify('bsqrt');
2146
990fb837
RGS
2147 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
2148 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
2149 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
58cde26e 2150
61f5c3f5 2151 # we need to limit the accuracy to protect against overflow
574bacfe 2152 my $fallback = 0;
990fb837
RGS
2153 my (@params,$scale);
2154 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2155
2156 return $x if $x->is_nan(); # error in _find_round_parameters?
61f5c3f5
T
2157
2158 # no rounding at all, so must use fallback
990fb837 2159 if (scalar @params == 0)
0716bf9b
JH
2160 {
2161 # simulate old behaviour
990fb837
RGS
2162 $params[0] = $self->div_scale(); # and round to it as accuracy
2163 $scale = $params[0]+4; # at least four more for proper round
2164 $params[2] = $r; # round mode by caller or undef
ee15d750 2165 $fallback = 1; # to clear a/p afterwards
0716bf9b 2166 }
61f5c3f5
T
2167 else
2168 {
2169 # the 4 below is empirical, and there might be cases where it is not
2170 # enough...
990fb837 2171 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
61f5c3f5
T
2172 }
2173
2174 # when user set globals, they would interfere with our calculation, so
9393ace2 2175 # disable them and later re-enable them
61f5c3f5
T
2176 no strict 'refs';
2177 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
b3abae2a 2178 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
61f5c3f5
T
2179 # we also need to disable any set A or P on $x (_find_round_parameters took
2180 # them already into account), since these would interfere, too
2181 delete $x->{_a}; delete $x->{_p};
9393ace2
JH
2182 # need to disable $upgrade in BigInt, to avoid deep recursion
2183 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
61f5c3f5 2184
9b924220
RGS
2185 my $i = $MBI->_copy( $x->{_m} );
2186 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2187 my $xas = Math::BigInt->bzero();
2188 $xas->{value} = $i;
2189
394e6ffb 2190 my $gs = $xas->copy()->bsqrt(); # some guess
b3abae2a 2191
9b924220 2192 if (($x->{_es} ne '-') # guess can't be accurate if there are
394e6ffb 2193 # digits after the dot
b3abae2a 2194 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
394e6ffb 2195 {
9b924220
RGS
2196 # exact result, copy result over to keep $x
2197 $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
2198 $x->bnorm();
091c87b1 2199 # shortcut to not run through _find_round_parameters again
990fb837 2200 if (defined $params[0])
61f5c3f5 2201 {
990fb837 2202 $x->bround($params[0],$params[2]); # then round accordingly
61f5c3f5
T
2203 }
2204 else
2205 {
990fb837 2206 $x->bfround($params[1],$params[2]); # then round accordingly
61f5c3f5
T
2207 }
2208 if ($fallback)
2209 {
2210 # clear a/p after round, since user did not request it
ef9466ea 2211 delete $x->{_a}; delete $x->{_p};
61f5c3f5 2212 }
9393ace2 2213 # re-enable A and P, upgrade is taken care of by "local"
b3abae2a 2214 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
61f5c3f5 2215 return $x;
394e6ffb 2216 }
2ab5f49d
T
2217
2218 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
9681bfa6 2219 # of the result by multiplying the input by 100 and then divide the integer
2ab5f49d 2220 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
9b924220
RGS
2221
2222 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
2223 my $y1 = $MBI->_copy($x->{_m});
2224
2225 my $length = $MBI->_len($y1);
2226
2227 # Now calculate how many digits the result of sqrt(y1) would have
2228 my $digits = int($length / 2);
2229
2230 # But we need at least $scale digits, so calculate how many are missing
2231 my $shift = $scale - $digits;
2232
0dceeee6
RGS
2233 # This happens if the input had enough digits
2234 # (we take care of integer guesses above)
2235 $shift = 0 if $shift < 0;
9b924220
RGS
2236
2237 # Multiply in steps of 100, by shifting left two times the "missing" digits
2238 my $s2 = $shift * 2;
2239
2ab5f49d
T
2240 # We now make sure that $y1 has the same odd or even number of digits than
2241 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
2242 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
2243 # steps of 10. The length of $x does not count, since an even or odd number
2244 # of digits before the dot is not changed by adding an even number of digits
2245 # after the dot (the result is still odd or even digits long).
9b924220
RGS
2246 $s2++ if $MBI->_is_odd($x->{_e});
2247
2248 $MBI->_lsft( $y1, $MBI->_new($s2), 10);
2249
2ab5f49d 2250 # now take the square root and truncate to integer
9b924220
RGS
2251 $y1 = $MBI->_sqrt($y1);
2252
2ab5f49d
T
2253 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
2254 # result, which is than later rounded to the desired scale.
990fb837
RGS
2255
2256 # calculate how many zeros $x had after the '.' (or before it, depending
9b924220
RGS
2257 # on sign of $dat, the result should have half as many:
2258 my $dat = $MBI->_num($x->{_e});
2259 $dat = -$dat if $x->{_es} eq '-';
2260 $dat += $length;
990fb837
RGS
2261
2262 if ($dat > 0)
2263 {
2264 # no zeros after the dot (e.g. 1.23, 0.49 etc)
2265 # preserve half as many digits before the dot than the input had
2266 # (but round this "up")
2267 $dat = int(($dat+1)/2);
2268 }
2269 else
2270 {
2271 $dat = int(($dat)/2);
2272 }
9b924220
RGS
2273 $dat -= $MBI->_len($y1);
2274 if ($dat < 0)
2275 {
2276 $dat = abs($dat);
2277 $x->{_e} = $MBI->_new( $dat );
2278 $x->{_es} = '-';
2279 }
2280 else
2281 {
2282 $x->{_e} = $MBI->_new( $dat );
2283 $x->{_es} = '+';
2284 }
2ab5f49d 2285 $x->{_m} = $y1;
9b924220 2286 $x->bnorm();
61f5c3f5 2287
091c87b1 2288 # shortcut to not run through _find_round_parameters again
990fb837 2289 if (defined $params[0])
61f5c3f5 2290 {
990fb837 2291 $x->bround($params[0],$params[2]); # then round accordingly
61f5c3f5
T
2292 }
2293 else
2294 {
990fb837 2295 $x->bfround($params[1],$params[2]); # then round accordingly
61f5c3f5 2296 }
574bacfe
JH
2297 if ($fallback)
2298 {
2299 # clear a/p after round, since user did not request it
ef9466ea 2300 delete $x->{_a}; delete $x->{_p};
574bacfe 2301 }
61f5c3f5 2302 # restore globals
b3abae2a 2303 $$abr = $ab; $$pbr = $pb;
574bacfe 2304 $x;
58cde26e
JH
2305 }
2306
b3abae2a
JH
2307sub bfac
2308 {
28df3e88 2309 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
091c87b1 2310 # compute factorial number, modifies first argument
b3abae2a 2311
b282a552
T
2312 # set up parameters
2313 my ($self,$x,@r) = (ref($_[0]),@_);
2314 # objectify is costly, so avoid it
2315 ($self,$x,@r) = objectify(1,@_) if !ref($x);
2316
50109ad0
RGS
2317 # inf => inf
2318 return $x if $x->modify('bfac') || $x->{sign} eq '+inf';
2319
28df3e88
JH
2320 return $x->bnan()
2321 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
9b924220 2322 ($x->{_es} ne '+')); # digits after dot?
b3abae2a 2323
b3abae2a 2324 # use BigInt's bfac() for faster calc
9b924220 2325 if (! $MBI->_is_zero($x->{_e}))
091c87b1 2326 {
9b924220
RGS
2327 $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
2328 $x->{_e} = $MBI->_zero(); # normalize
2329 $x->{_es} = '+';
091c87b1 2330 }
9b924220 2331 $MBI->_fac($x->{_m}); # calculate factorial
091c87b1 2332 $x->bnorm()->round(@r); # norm again and round result
b3abae2a
JH
2333 }
2334
9393ace2
JH
2335sub _pow
2336 {
60a1aa19
T
2337 # Calculate a power where $y is a non-integer, like 2 ** 0.3
2338 my ($x,$y,@r) = @_;
9393ace2
JH
2339 my $self = ref($x);
2340
2341 # if $y == 0.5, it is sqrt($x)
27e7b8bb 2342 $HALF = $self->new($HALF) unless ref($HALF);
60a1aa19 2343 return $x->bsqrt(@r,$y) if $y->bcmp($HALF) == 0;
9393ace2 2344
990fb837
RGS
2345 # Using:
2346 # a ** x == e ** (x * ln a)
2347
9393ace2 2348 # u = y * ln x
990fb837
RGS
2349 # _ _
2350 # Taylor: | u u^2 u^3 |
2351 # x ** y = 1 + | --- + --- + ----- + ... |
2352 # |_ 1 1*2 1*2*3 _|
9393ace2
JH
2353
2354 # we need to limit the accuracy to protect against overflow
2355 my $fallback = 0;
990fb837 2356 my ($scale,@params);
60a1aa19 2357 ($x,@params) = $x->_find_round_parameters(@r);
990fb837
RGS
2358
2359 return $x if $x->is_nan(); # error in _find_round_parameters?
9393ace2
JH
2360
2361 # no rounding at all, so must use fallback
990fb837 2362 if (scalar @params == 0)
9393ace2
JH
2363 {
2364 # simulate old behaviour
990fb837
RGS
2365 $params[0] = $self->div_scale(); # and round to it as accuracy
2366 $params[1] = undef; # disable P
2367 $scale = $params[0]+4; # at least four more for proper round
60a1aa19 2368 $params[2] = $r[2]; # round mode by caller or undef
9393ace2
JH
2369 $fallback = 1; # to clear a/p afterwards
2370 }
2371 else
2372 {
2373 # the 4 below is empirical, and there might be cases where it is not
2374 # enough...
990fb837 2375 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
9393ace2
JH
2376 }
2377
2378 # when user set globals, they would interfere with our calculation, so
56d9de68 2379 # disable them and later re-enable them
9393ace2
JH
2380 no strict 'refs';
2381 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2382 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2383 # we also need to disable any set A or P on $x (_find_round_parameters took
2384 # them already into account), since these would interfere, too
2385 delete $x->{_a}; delete $x->{_p};
2386 # need to disable $upgrade in BigInt, to avoid deep recursion
2387 local $Math::BigInt::upgrade = undef;
2388
2389 my ($limit,$v,$u,$below,$factor,$next,$over);
2390
990fb837 2391 $u = $x->copy()->blog(undef,$scale)->bmul($y);
9393ace2
JH
2392 $v = $self->bone(); # 1
2393 $factor = $self->new(2); # 2
2394 $x->bone(); # first term: 1
2395
2396 $below = $v->copy();
2397 $over = $u->copy();
ae161977 2398
9393ace2
JH
2399 $limit = $self->new("1E-". ($scale-1));
2400 #my $steps = 0;
2401 while (3 < 5)
2402 {
2403 # we calculate the next term, and add it to the last
2404 # when the next term is below our limit, it won't affect the outcome
7d193e39 2405 # anymore, so we stop:
9393ace2 2406 $next = $over->copy()->bdiv($below,$scale);
990fb837 2407 last if $next->bacmp($limit) <= 0;
9393ace2 2408 $x->badd($next);
9393ace2
JH
2409 # calculate things for the next term
2410 $over *= $u; $below *= $factor; $factor->binc();
9b924220
RGS
2411
2412 last if $x->{sign} !~ /^[-+]$/;
2413
9393ace2
JH
2414 #$steps++;
2415 }
2416
091c87b1 2417 # shortcut to not run through _find_round_parameters again
990fb837 2418 if (defined $params[0])
9393ace2 2419 {
990fb837 2420 $x->bround($params[0],$params[2]); # then round accordingly
9393ace2
JH
2421 }
2422 else
2423 {
990fb837 2424 $x->bfround($params[1],$params[2]); # then round accordingly
9393ace2
JH
2425 }
2426 if ($fallback)
2427 {
2428 # clear a/p after round, since user did not request it
ef9466ea 2429 delete $x->{_a}; delete $x->{_p};
9393ace2
JH
2430 }
2431 # restore globals
2432 $$abr = $ab; $$pbr = $pb;
2433 $x;
2434 }
2435
58cde26e
JH
2436sub bpow
2437 {
2438 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2439 # compute power of two numbers, second arg is used as integer
2440 # modifies first argument
2441
f9a08e12
JH
2442 # set up parameters
2443 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2444 # objectify is costly, so avoid it
2445 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2446 {
2447 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2448 }
58cde26e 2449
50109ad0
RGS
2450 return $x if $x->modify('bpow');
2451
58cde26e 2452 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
2d2b2744
T
2453 return $x if $x->{sign} =~ /^[+-]inf$/;
2454
ae161977
RGS
2455 # cache the result of is_zero
2456 my $y_is_zero = $y->is_zero();
2457 return $x->bone() if $y_is_zero;
58cde26e 2458 return $x if $x->is_one() || $y->is_one();
9393ace2 2459
ae161977 2460 my $x_is_zero = $x->is_zero();
2d2b2744 2461 return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int(); # non-integer power
9393ace2 2462
ae161977 2463 my $y1 = $y->as_number()->{value}; # make MBI part
9b924220 2464
394e6ffb 2465 # if ($x == -1)
9b924220 2466 if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
58cde26e
JH
2467 {
2468 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
9b924220 2469 return $MBI->_is_odd($y1) ? $x : $x->babs(1);
288d023a 2470 }
ae161977 2471 if ($x_is_zero)
28df3e88
JH
2472 {
2473 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
9b924220
RGS
2474 # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
2475 return $x->binf();
28df3e88 2476 }
58cde26e 2477
9b924220 2478 my $new_sign = '+';
ae161977 2479 $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
9b924220 2480
58cde26e 2481 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
9b924220 2482 $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
ae161977 2483 $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
9b924220
RGS
2484
2485 $x->{sign} = $new_sign;
58cde26e
JH
2486 $x->bnorm();
2487 if ($y->{sign} eq '-')
2488 {
2489 # modify $x in place!
ae161977 2490 my $z = $x->copy(); $x->bone();
7b29e1e6 2491 return scalar $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
a0d0e21e 2492 }
28df3e88 2493 $x->round($a,$p,$r,$y);
58cde26e
JH
2494 }
2495
80365507
T
2496sub bmodpow
2497 {
2498 # takes a very large number to a very large exponent in a given very
c4a6f826 2499 # large modulus, quickly, thanks to binary exponentiation. Supports
80365507
T
2500 # negative exponents.
2501 my ($self,$num,$exp,$mod,@r) = objectify(3,@_);
2502
2503 return $num if $num->modify('bmodpow');
2504
2505 # check modulus for valid values
2506 return $num->bnan() if ($mod->{sign} ne '+' # NaN, - , -inf, +inf
2507 || $mod->is_zero());
2508
2509 # check exponent for valid values
2510 if ($exp->{sign} =~ /\w/)
2511 {
2512 # i.e., if it's NaN, +inf, or -inf...
2513 return $num->bnan();
2514 }
2515
2516 $num->bmodinv ($mod) if ($exp->{sign} eq '-');
2517
2518 # check num for valid values (also NaN if there was no inverse but $exp < 0)
2519 return $num->bnan() if $num->{sign} !~ /^[+-]$/;
2520
2521 # $mod is positive, sign on $exp is ignored, result also positive
2522
2523 # XXX TODO: speed it up when all three numbers are integers
2524 $num->bpow($exp)->bmod($mod);
2525 }
2526
58cde26e 2527###############################################################################
fdb4b05f
T
2528# trigonometric functions
2529
20e2035c 2530# helper function for bpi() and batan2(), calculates arcus tanges (1/x)
fdb4b05f 2531
20e2035c 2532sub _atan_inv
fdb4b05f 2533 {
20e2035c
T
2534 # return a/b so that a/b approximates atan(1/x) to at least limit digits
2535 my ($self, $x, $limit) = @_;
fdb4b05f 2536
20e2035c
T
2537 # Taylor: x^3 x^5 x^7 x^9
2538 # atan = x - --- + --- - --- + --- - ...
2539 # 3 5 7 9
fdb4b05f 2540
20e2035c
T
2541 # 1 1 1 1
2542 # atan 1/x = - - ------- + ------- - ------- + ...
2543 # x x^3 * 3 x^5 * 5 x^7 * 7
2544
2545 # 1 1 1 1
2546 # atan 1/x = - - --------- + ---------- - ----------- + ...
2547 # 5 3 * 125 5 * 3125 7 * 78125
2548
2549 # Subtraction/addition of a rational:
2550
2551 # 5 7 5*3 +- 7*4
2552 # - +- - = ----------
2553 # 4 3 4*3
2554
2555 # Term: N N+1
2556 #
2557 # a 1 a * d * c +- b
2558 # ----- +- ------------------ = ----------------
2559 # b d * c b * d * c
2560
2561 # since b1 = b0 * (d-2) * c
2562
2563 # a 1 a * d +- b / c
2564 # ----- +- ------------------ = ----------------
2565 # b d * c b * d
2566
2567 # and d = d + 2
2568 # and c = c * x * x
2569
2570 # u = d * c
2571 # stop if length($u) > limit
2572 # a = a * u +- b
2573 # b = b * u
2574 # d = d + 2
2575 # c = c * x * x
2576 # sign = 1 - sign
2577
2578 my $a = $MBI->_one();
30afc38d 2579 my $b = $MBI->_copy($x);
fdb4b05f 2580
30afc38d 2581 my $x2 = $MBI->_mul( $MBI->_copy($x), $b); # x2 = x * x
20e2035c 2582 my $d = $MBI->_new( 3 ); # d = 3
30afc38d 2583 my $c = $MBI->_mul( $MBI->_copy($x), $x2); # c = x ^ 3
20e2035c
T
2584 my $two = $MBI->_new( 2 );
2585
2586 # run the first step unconditionally
2587 my $u = $MBI->_mul( $MBI->_copy($d), $c);
2588 $a = $MBI->_mul($a, $u);
2589 $a = $MBI->_sub($a, $b);
2590 $b = $MBI->_mul($b, $u);
2591 $d = $MBI->_add($d, $two);
2592 $c = $MBI->_mul($c, $x2);
2593
2594 # a is now a * (d-3) * c
2595 # b is now b * (d-2) * c
2596
2597 # run the second step unconditionally
2598 $u = $MBI->_mul( $MBI->_copy($d), $c);
2599 $a = $MBI->_mul($a, $u);
2600 $a = $MBI->_add($a, $b);
2601 $b = $MBI->_mul($b, $u);
2602 $d = $MBI->_add($d, $two);
2603 $c = $MBI->_mul($c, $x2);
2604
2605 # a is now a * (d-3) * (d-5) * c * c
2606 # b is now b * (d-2) * (d-4) * c * c
2607
2608 # so we can remove c * c from both a and b to shorten the numbers involved:
2609 $a = $MBI->_div($a, $x2);
2610 $b = $MBI->_div($b, $x2);
2611 $a = $MBI->_div($a, $x2);
2612 $b = $MBI->_div($b, $x2);
2613
2614# my $step = 0;
2615 my $sign = 0; # 0 => -, 1 => +
2616 while (3 < 5)
fdb4b05f 2617 {
20e2035c
T
2618# $step++;
2619# if (($i++ % 100) == 0)
2620# {
2621# print "a=",$MBI->_str($a),"\n";
2622# print "b=",$MBI->_str($b),"\n";
2623# }
2624# print "d=",$MBI->_str($d),"\n";
2625# print "x2=",$MBI->_str($x2),"\n";
2626# print "c=",$MBI->_str($c),"\n";
2627
2628 my $u = $MBI->_mul( $MBI->_copy($d), $c);
2629 # use _alen() for libs like GMP where _len() would be O(N^2)
2630 last if $MBI->_alen($u) > $limit;
2631 my ($bc,$r) = $MBI->_div( $MBI->_copy($b), $c);
2632 if ($MBI->_is_zero($r))
fdb4b05f 2633 {
20e2035c
T
2634 # b / c is an integer, so we can remove c from all terms
2635 # this happens almost every time:
2636 $a = $MBI->_mul($a, $d);
2637 $a = $MBI->_sub($a, $bc) if $sign == 0;
2638 $a = $MBI->_add($a, $bc) if $sign == 1;
2639 $b = $MBI->_mul($b, $d);
fdb4b05f
T
2640 }
2641 else
2642 {
20e2035c
T
2643 # b / c is not an integer, so we keep c in the terms
2644 # this happens very rarely, for instance for x = 5, this happens only
2645 # at the following steps:
2646 # 1, 5, 14, 32, 72, 157, 340, ...
2647 $a = $MBI->_mul($a, $u);
2648 $a = $MBI->_sub($a, $b) if $sign == 0;
2649 $a = $MBI->_add($a, $b) if $sign == 1;
2650 $b = $MBI->_mul($b, $u);
fdb4b05f 2651 }
20e2035c
T
2652 $d = $MBI->_add($d, $two);
2653 $c = $MBI->_mul($c, $x2);
2654 $sign = 1 - $sign;
2655
fdb4b05f
T
2656 }
2657
30afc38d 2658# print "Took $step steps for ", $MBI->_str($x),"\n";
20e2035c
T
2659# print "a=",$MBI->_str($a),"\n"; print "b=",$MBI->_str($b),"\n";
2660 # return a/b so that a/b approximates atan(1/x)
2661 ($a,$b);
fdb4b05f
T
2662 }
2663
2664sub bpi
2665 {
fdb4b05f 2666 my ($self,$n) = @_;
80365507
T
2667 if (@_ == 0)
2668 {
2669 $self = $class;
2670 }
fdb4b05f
T
2671 if (@_ == 1)
2672 {
2673 # called like Math::BigFloat::bpi(10);
2674 $n = $self; $self = $class;
60a1aa19
T
2675 # called like Math::BigFloat->bpi();
2676 $n = undef if $n eq 'Math::BigFloat';
fdb4b05f
T
2677 }
2678 $self = ref($self) if ref($self);
36ec1dfe 2679 my $fallback = defined $n ? 0 : 1;
fdb4b05f
T
2680 $n = 40 if !defined $n || $n < 1;
2681
20e2035c
T
2682 # after 黃見利 (Hwang Chien-Lih) (1997)
2683 # pi/4 = 183 * atan(1/239) + 32 * atan(1/1023) – 68 * atan(1/5832)
2684 # + 12 * atan(1/110443) - 12 * atan(1/4841182) - 100 * atan(1/6826318)
2685
2686 # a few more to prevent rounding errors
2687 $n += 4;
2688
30afc38d
T
2689 my ($a,$b) = $self->_atan_inv( $MBI->_new(239),$n);
2690 my ($c,$d) = $self->_atan_inv( $MBI->_new(1023),$n);
2691 my ($e,$f) = $self->_atan_inv( $MBI->_new(5832),$n);
2692 my ($g,$h) = $self->_atan_inv( $MBI->_new(110443),$n);
2693 my ($i,$j) = $self->_atan_inv( $MBI->_new(4841182),$n);
2694 my ($k,$l) = $self->_atan_inv( $MBI->_new(6826318),$n);
20e2035c
T
2695
2696 $MBI->_mul($a, $MBI->_new(732));
2697 $MBI->_mul($c, $MBI->_new(128));
2698 $MBI->_mul($e, $MBI->_new(272));
2699 $MBI->_mul($g, $MBI->_new(48));
2700 $MBI->_mul($i, $MBI->_new(48));
2701 $MBI->_mul($k, $MBI->_new(400));
2702
2703 my $x = $self->bone(); $x->{_m} = $a; my $x_d = $self->bone(); $x_d->{_m} = $b;
2704 my $y = $self->bone(); $y->{_m} = $c; my $y_d = $self->bone(); $y_d->{_m} = $d;
2705 my $z = $self->bone(); $z->{_m} = $e; my $z_d = $self->bone(); $z_d->{_m} = $f;
2706 my $u = $self->bone(); $u->{_m} = $g; my $u_d = $self->bone(); $u_d->{_m} = $h;
2707 my $v = $self->bone(); $v->{_m} = $i; my $v_d = $self->bone(); $v_d->{_m} = $j;
2708 my $w = $self->bone(); $w->{_m} = $k; my $w_d = $self->bone(); $w_d->{_m} = $l;
2709 $x->bdiv($x_d, $n);
2710 $y->bdiv($y_d, $n);
2711 $z->bdiv($z_d, $n);
2712 $u->bdiv($u_d, $n);
2713 $v->bdiv($v_d, $n);
2714 $w->bdiv($w_d, $n);
2715
36ec1dfe
T
2716 delete $x->{_a}; delete $y->{_a}; delete $z->{_a};
2717 delete $u->{_a}; delete $v->{_a}; delete $w->{_a};
20e2035c
T
2718 $x->badd($y)->bsub($z)->badd($u)->bsub($v)->bsub($w);
2719
36ec1dfe
T
2720 $x->bround($n-4);
2721 delete $x->{_a} if $fallback == 1;
2722 $x;
fdb4b05f
T
2723 }
2724
60a1aa19
T
2725sub bcos
2726 {
2727 # Calculate a cosinus of x.
2728 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2729
2730 # Taylor: x^2 x^4 x^6 x^8
2731 # cos = 1 - --- + --- - --- + --- ...
2732 # 2! 4! 6! 8!
2733
2734 # we need to limit the accuracy to protect against overflow
2735 my $fallback = 0;
2736 my ($scale,@params);
2737 ($x,@params) = $x->_find_round_parameters(@r);
2738
2739 # constant object or error in _find_round_parameters?
2740 return $x if $x->modify('bcos') || $x->is_nan();
2741
2742 return $x->bone(@r) if $x->is_zero();
2743
2744 # no rounding at all, so must use fallback
2745 if (scalar @params == 0)
2746 {
2747 # simulate old behaviour
2748 $params[0] = $self->div_scale(); # and round to it as accuracy
2749 $params[1] = undef; # disable P
2750 $scale = $params[0]+4; # at least four more for proper round
2751 $params[2] = $r[2]; # round mode by caller or undef
2752 $fallback = 1; # to clear a/p afterwards
2753 }
2754 else
2755 {
2756 # the 4 below is empirical, and there might be cases where it is not
2757 # enough...
2758 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2759 }
2760
2761 # when user set globals, they would interfere with our calculation, so
2762 # disable them and later re-enable them
2763 no strict 'refs';
2764 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2765 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2766 # we also need to disable any set A or P on $x (_find_round_parameters took
2767 # them already into account), since these would interfere, too
2768 delete $x->{_a}; delete $x->{_p};
2769 # need to disable $upgrade in BigInt, to avoid deep recursion
2770 local $Math::BigInt::upgrade = undef;
2771
2772 my $last = 0;
2773 my $over = $x * $x; # X ^ 2
2774 my $x2 = $over->copy(); # X ^ 2; difference between terms
2775 my $sign = 1; # start with -=
2776 my $below = $self->new(2); my $factorial = $self->new(3);
36ec1dfe 2777 $x->bone(); delete $x->{_a}; delete $x->{_p};
60a1aa19
T
2778
2779 my $limit = $self->new("1E-". ($scale-1));
2780 #my $steps = 0;
2781 while (3 < 5)
2782 {
2783 # we calculate the next term, and add it to the last
2784 # when the next term is below our limit, it won't affect the outcome
2785 # anymore, so we stop:
2786 my $next = $over->copy()->bdiv($below,$scale);
2787 last if $next->bacmp($limit) <= 0;
2788
2789 if ($sign == 0)
2790 {
2791 $x->badd($next);
2792 }
2793 else
2794 {
2795 $x->bsub($next);
2796 }
2797 $sign = 1-$sign; # alternate
2798 # calculate things for the next term
2799 $over->bmul($x2); # $x*$x
2800 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2801 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2802 }
2803
2804 # shortcut to not run through _find_round_parameters again
2805 if (defined $params[0])
2806 {
2807 $x->bround($params[0],$params[2]); # then round accordingly
2808 }
2809 else
2810 {
2811 $x->bfround($params[1],$params[2]); # then round accordingly
2812 }
2813 if ($fallback)
2814 {
2815 # clear a/p after round, since user did not request it
2816 delete $x->{_a}; delete $x->{_p};
2817 }
2818 # restore globals
2819 $$abr = $ab; $$pbr = $pb;
2820 $x;
2821 }
2822
2823sub bsin
2824 {
2825 # Calculate a sinus of x.
2826 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2827
2828 # taylor: x^3 x^5 x^7 x^9
2829 # sin = x - --- + --- - --- + --- ...
2830 # 3! 5! 7! 9!
2831
2832 # we need to limit the accuracy to protect against overflow
2833 my $fallback = 0;
2834 my ($scale,@params);
2835 ($x,@params) = $x->_find_round_parameters(@r);
2836
2837 # constant object or error in _find_round_parameters?
2838 return $x if $x->modify('bsin') || $x->is_nan();
2839
2840 return $x->bzero(@r) if $x->is_zero();
2841
2842 # no rounding at all, so must use fallback
2843 if (scalar @params == 0)
2844 {
2845 # simulate old behaviour
2846 $params[0] = $self->div_scale(); # and round to it as accuracy
2847 $params[1] = undef; # disable P
2848 $scale = $params[0]+4; # at least four more for proper round
2849 $params[2] = $r[2]; # round mode by caller or undef
2850 $fallback = 1; # to clear a/p afterwards
2851 }
2852 else
2853 {
2854 # the 4 below is empirical, and there might be cases where it is not
2855 # enough...
2856 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2857 }
2858
2859 # when user set globals, they would interfere with our calculation, so
2860 # disable them and later re-enable them
2861 no strict 'refs';
2862 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2863 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2864 # we also need to disable any set A or P on $x (_find_round_parameters took
2865 # them already into account), since these would interfere, too
2866 delete $x->{_a}; delete $x->{_p};
2867 # need to disable $upgrade in BigInt, to avoid deep recursion
2868 local $Math::BigInt::upgrade = undef;
2869
2870 my $last = 0;
2871 my $over = $x * $x; # X ^ 2
2872 my $x2 = $over->copy(); # X ^ 2; difference between terms
2873 $over->bmul($x); # X ^ 3 as starting value
2874 my $sign = 1; # start with -=
2875 my $below = $self->new(6); my $factorial = $self->new(4);
36ec1dfe 2876 delete $x->{_a}; delete $x->{_p};
60a1aa19
T
2877
2878 my $limit = $self->new("1E-". ($scale-1));
2879 #my $steps = 0;
2880 while (3 < 5)
2881 {
2882 # we calculate the next term, and add it to the last
2883 # when the next term is below our limit, it won't affect the outcome
2884 # anymore, so we stop:
2885 my $next = $over->copy()->bdiv($below,$scale);
2886 last if $next->bacmp($limit) <= 0;
2887
2888 if ($sign == 0)
2889 {
2890 $x->badd($next);
2891 }
2892 else
2893 {
2894 $x->bsub($next);
2895 }
2896 $sign = 1-$sign; # alternate
2897 # calculate things for the next term
2898 $over->bmul($x2); # $x*$x
2899 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2900 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2901 }
2902
2903 # shortcut to not run through _find_round_parameters again
2904 if (defined $params[0])
2905 {
2906 $x->bround($params[0],$params[2]); # then round accordingly
2907 }
2908 else
2909 {
2910 $x->bfround($params[1],$params[2]); # then round accordingly
2911 }
2912 if ($fallback)
2913 {
2914 # clear a/p after round, since user did not request it
2915 delete $x->{_a}; delete $x->{_p};
2916 }
2917 # restore globals
2918 $$abr = $ab; $$pbr = $pb;
2919 $x;
2920 }
2921
20e2035c
T
2922sub batan2
2923 {
30afc38d 2924 # calculate arcus tangens of ($y/$x)
20e2035c
T
2925
2926 # set up parameters
30afc38d 2927 my ($self,$y,$x,@r) = (ref($_[0]),@_);
20e2035c
T
2928 # objectify is costly, so avoid it
2929 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2930 {
30afc38d 2931 ($self,$y,$x,@r) = objectify(2,@_);
20e2035c
T
2932 }
2933
30afc38d 2934 return $y if $y->modify('batan2');
20e2035c 2935
30afc38d 2936 return $y->bnan() if ($y->{sign} eq $nan) || ($x->{sign} eq $nan);
20e2035c 2937
30afc38d
T
2938 # Y X
2939 # 0 0 result is 0
2940 # 0 +x result is 0
0dceeee6
RGS
2941 # ? inf result is 0
2942 return $y->bzero(@r) if ($x->is_inf('+') && !$y->is_inf()) || ($y->is_zero() && $x->{sign} eq '+');
2943
2944 # Y X
2945 # != 0 -inf result is +- pi
2946 if ($x->is_inf() || $y->is_inf())
2947 {
2948 # calculate PI
2949 my $pi = $self->bpi(@r);
2950 if ($y->is_inf())
2951 {
2952 # upgrade to BigRat etc.
2953 return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
2954 if ($x->{sign} eq '-inf')
2955 {
2956 # calculate 3 pi/4
2957 $MBI->_mul($pi->{_m}, $MBI->_new(3));
2958 $MBI->_div($pi->{_m}, $MBI->_new(4));
2959 }
2960 elsif ($x->{sign} eq '+inf')
2961 {
2962 # calculate pi/4
2963 $MBI->_div($pi->{_m}, $MBI->_new(4));
2964 }
2965 else
2966 {
2967 # calculate pi/2
2968 $MBI->_div($pi->{_m}, $MBI->_new(2));
2969 }
2970 $y->{sign} = substr($y->{sign},0,1); # keep +/-
2971 }
2972 # modify $y in place
2973 $y->{_m} = $pi->{_m};
2974 $y->{_e} = $pi->{_e};
2975 $y->{_es} = $pi->{_es};
2976 # keep the sign of $y
2977 return $y;
2978 }
2979
2980 return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
30afc38d
T
2981
2982 # Y X
2983 # 0 -x result is PI
2984 if ($y->is_zero())
20e2035c 2985 {
30afc38d
T
2986 # calculate PI
2987 my $pi = $self->bpi(@r);
0dceeee6 2988 # modify $y in place
30afc38d
T
2989 $y->{_m} = $pi->{_m};
2990 $y->{_e} = $pi->{_e};
2991 $y->{_es} = $pi->{_es};
2992 $y->{sign} = '+';
2993 return $y;
2994 }
2995
2996 # Y X
2997 # +y 0 result is PI/2
2998 # -y 0 result is -PI/2
0dceeee6 2999 if ($x->is_zero())
30afc38d
T
3000 {
3001 # calculate PI/2
3002 my $pi = $self->bpi(@r);
0dceeee6 3003 # modify $y in place
30afc38d
T
3004 $y->{_m} = $pi->{_m};
3005 $y->{_e} = $pi->{_e};
3006 $y->{_es} = $pi->{_es};
3007 # -y => -PI/2, +y => PI/2
30afc38d
T
3008 $MBI->_div($y->{_m}, $MBI->_new(2));
3009 return $y;
20e2035c
T
3010 }
3011
30afc38d
T
3012 # we need to limit the accuracy to protect against overflow
3013 my $fallback = 0;
3014 my ($scale,@params);
3015 ($y,@params) = $y->_find_round_parameters(@r);
3016
3017 # error in _find_round_parameters?
3018 return $y if $y->is_nan();
3019
3020 # no rounding at all, so must use fallback
3021 if (scalar @params == 0)
3022 {
3023 # simulate old behaviour
3024 $params[0] = $self->div_scale(); # and round to it as accuracy
3025 $params[1] = undef; # disable P
3026 $scale = $params[0]+4; # at least four more for proper round
3027 $params[2] = $r[2]; # round mode by caller or undef
3028 $fallback = 1; # to clear a/p afterwards
3029 }
3030 else
3031 {
3032 # the 4 below is empirical, and there might be cases where it is not
3033 # enough...
3034 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3035 }
3036
3037 # inlined is_one() && is_one('-')
3038 if ($MBI->_is_one($y->{_m}) && $MBI->_is_zero($y->{_e}))
3039 {
3040 # shortcut: 1 1 result is PI/4
3041 # inlined is_one() && is_one('-')
3042 if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
3043 {
3044 # 1,1 => PI/4
3045 my $pi_4 = $self->bpi( $scale - 3);
0dceeee6 3046 # modify $y in place
30afc38d
T
3047 $y->{_m} = $pi_4->{_m};
3048 $y->{_e} = $pi_4->{_e};
3049 $y->{_es} = $pi_4->{_es};
3050 # 1 1 => +
3051 # -1 1 => -
3052 # 1 -1 => -
3053 # -1 -1 => +
3054 $y->{sign} = $x->{sign} eq $y->{sign} ? '+' : '-';
3055 $MBI->_div($y->{_m}, $MBI->_new(4));
3056 return $y;
3057 }
3058 # shortcut: 1 int(X) result is _atan_inv(X)
3059
3060 # is integer
3061 if ($x->{_es} eq '+')
3062 {
3063 my $x1 = $MBI->_copy($x->{_m});
3064 $MBI->_lsft($x1, $x->{_e},10) unless $MBI->_is_zero($x->{_e});
3065
3066 my ($a,$b) = $self->_atan_inv($x1, $scale);
3067 my $y_sign = $y->{sign};
3068 # calculate A/B
3069 $y->bone(); $y->{_m} = $a; my $y_d = $self->bone(); $y_d->{_m} = $b;
3070 $y->bdiv($y_d, @r);
3071 $y->{sign} = $y_sign;
3072 return $y;
3073 }
3074 }
3075
3076 # handle all other cases
3077 # X Y
3078 # +x +y 0 to PI/2
3079 # -x +y PI/2 to PI
3080 # +x -y 0 to -PI/2
3081 # -x -y -PI/2 to -PI
3082
3083 my $y_sign = $y->{sign};
20e2035c
T
3084
3085 # divide $x by $y
30afc38d
T
3086 $y->bdiv($x, $scale) unless $x->is_one();
3087 $y->batan(@r);
20e2035c 3088
30afc38d
T
3089 # restore sign
3090 $y->{sign} = $y_sign;
20e2035c 3091
30afc38d 3092 $y;
20e2035c
T
3093 }
3094
60a1aa19
T
3095sub batan
3096 {
3097 # Calculate a arcus tangens of x.
3098 my ($x,@r) = @_;
3099 my $self = ref($x);
3100
3101 # taylor: x^3 x^5 x^7 x^9
3102 # atan = x - --- + --- - --- + --- ...
3103 # 3 5 7 9
3104
60a1aa19
T
3105 # we need to limit the accuracy to protect against overflow
3106 my $fallback = 0;
3107 my ($scale,@params);
3108 ($x,@params) = $x->_find_round_parameters(@r);
3109
3110 # constant object or error in _find_round_parameters?
3111 return $x if $x->modify('batan') || $x->is_nan();
3112
30afc38d
T
3113 if ($x->{sign} =~ /^[+-]inf\z/)
3114 {
3115 # +inf result is PI/2
3116 # -inf result is -PI/2
3117 # calculate PI/2
3118 my $pi = $self->bpi(@r);
3119 # modify $x in place
3120 $x->{_m} = $pi->{_m};
3121 $x->{_e} = $pi->{_e};
3122 $x->{_es} = $pi->{_es};
3123 # -y => -PI/2, +y => PI/2
3124 $x->{sign} = substr($x->{sign},0,1); # +inf => +
3125 $MBI->_div($x->{_m}, $MBI->_new(2));
3126 return $x;
3127 }
3128
3129 return $x->bzero(@r) if $x->is_zero();
3130
60a1aa19
T
3131 # no rounding at all, so must use fallback
3132 if (scalar @params == 0)
3133 {
3134 # simulate old behaviour
3135 $params[0] = $self->div_scale(); # and round to it as accuracy
3136 $params[1] = undef; # disable P
3137 $scale = $params[0]+4; # at least four more for proper round
3138 $params[2] = $r[2]; # round mode by caller or undef
3139 $fallback = 1; # to clear a/p afterwards
3140 }
3141 else
3142 {
3143 # the 4 below is empirical, and there might be cases where it is not
3144 # enough...
3145 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3146 }
3147
30afc38d
T
3148 # 1 or -1 => PI/4
3149 # inlined is_one() && is_one('-')
3150 if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
3151 {
3152 my $pi = $self->bpi($scale - 3);
3153 # modify $x in place
3154 $x->{_m} = $pi->{_m};
3155 $x->{_e} = $pi->{_e};
3156 $x->{_es} = $pi->{_es};
3157 # leave the sign of $x alone (+1 => +PI/4, -1 => -PI/4)
3158 $MBI->_div($x->{_m}, $MBI->_new(4));
3159 return $x;
3160 }
3161
3162 # This series is only valid if -1 < x < 1, so for other x we need to
c97ef841 3163 # to calculate PI/2 - atan(1/x):
30afc38d
T
3164 my $one = $MBI->_new(1);
3165 my $pi = undef;
3166 if ($x->{_es} eq '+' && ($MBI->_acmp($x->{_m},$one) >= 0))
3167 {
3168 # calculate PI/2
3169 $pi = $self->bpi($scale - 3);
3170 $MBI->_div($pi->{_m}, $MBI->_new(2));
3171 # calculate 1/$x:
3172 my $x_copy = $x->copy();
3173 # modify $x in place
3174 $x->bone(); $x->bdiv($x_copy,$scale);
3175 }
3176
60a1aa19
T
3177 # when user set globals, they would interfere with our calculation, so
3178 # disable them and later re-enable them
3179 no strict 'refs';
3180 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
3181 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
3182 # we also need to disable any set A or P on $x (_find_round_parameters took
3183 # them already into account), since these would interfere, too
3184 delete $x->{_a}; delete $x->{_p};
3185 # need to disable $upgrade in BigInt, to avoid deep recursion
3186 local $Math::BigInt::upgrade = undef;
3187
3188 my $last = 0;
3189 my $over = $x * $x; # X ^ 2
3190 my $x2 = $over->copy(); # X ^ 2; difference between terms
3191 $over->bmul($x); # X ^ 3 as starting value
3192 my $sign = 1; # start with -=
3193 my $below = $self->new(3);
3194 my $two = $self->new(2);
36ec1dfe 3195 delete $x->{_a}; delete $x->{_p};
60a1aa19
T
3196
3197 my $limit = $self->new("1E-". ($scale-1));
3198 #my $steps = 0;
3199 while (3 < 5)
3200 {
3201 # we calculate the next term, and add it to the last
3202 # when the next term is below our limit, it won't affect the outcome
3203 # anymore, so we stop:
3204 my $next = $over->copy()->bdiv($below,$scale);
3205 last if $next->bacmp($limit) <= 0;
3206
3207 if ($sign == 0)
3208 {
3209 $x->badd($next);
3210 }
3211 else
3212 {
3213 $x->bsub($next);
3214 }
3215 $sign = 1-$sign; # alternate
3216 # calculate things for the next term
3217 $over->bmul($x2); # $x*$x
3218 $below->badd($two); # n += 2
3219 }
3220
30afc38d
T
3221 if (defined $pi)
3222 {
3223 my $x_copy = $x->copy();
3224 # modify $x in place
3225 $x->{_m} = $pi->{_m};
3226 $x->{_e} = $pi->{_e};
3227 $x->{_es} = $pi->{_es};
3228 # PI/2 - $x
3229 $x->bsub($x_copy);
3230 }
3231
60a1aa19
T
3232 # shortcut to not run through _find_round_parameters again
3233 if (defined $params[0])
3234 {
3235 $x->bround($params[0],$params[2]); # then round accordingly
3236 }
3237 else
3238 {
3239 $x->bfround($params[1],$params[2]); # then round accordingly
3240 }
3241 if ($fallback)
3242 {
3243 # clear a/p after round, since user did not request it
3244 delete $x->{_a}; delete $x->{_p};
3245 }
3246 # restore globals
3247 $$abr = $ab; $$pbr = $pb;
3248 $x;
3249 }
3250
fdb4b05f 3251###############################################################################
58cde26e
JH
3252# rounding functions
3253
3254sub bfround
3255 {
3256 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
3257 # $n == 0 means round to integer
3258 # expects and returns normalized numbers!
ee15d750 3259 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
a0d0e21e 3260
b68b7ab1
T
3261 my ($scale,$mode) = $x->_scale_p(@_);
3262 return $x if !defined $scale || $x->modify('bfround'); # no-op
58cde26e 3263
574bacfe 3264 # never round a 0, +-inf, NaN
61f5c3f5
T
3265 if ($x->is_zero())
3266 {
3267 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
3268 return $x;
3269 }
3270 return $x if $x->{sign} !~ /^[+-]$/;
58cde26e 3271
ee15d750
JH
3272 # don't round if x already has lower precision
3273 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
3274
3275 $x->{_p} = $scale; # remember round in any case
ef9466ea 3276 delete $x->{_a}; # and clear A
58cde26e
JH
3277 if ($scale < 0)
3278 {
58cde26e 3279 # round right from the '.'
f9a08e12 3280
9b924220 3281 return $x if $x->{_es} eq '+'; # e >= 0 => nothing to round
f9a08e12 3282
58cde26e 3283 $scale = -$scale; # positive for simplicity
9b924220 3284 my $len = $MBI->_len($x->{_m}); # length of mantissa
f9a08e12
JH
3285
3286 # the following poses a restriction on _e, but if _e is bigger than a
3287 # scalar, you got other problems (memory etc) anyway
9b924220 3288 my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e}))); # digits after dot
58cde26e 3289 my $zad = 0; # zeros after dot
f9a08e12 3290 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
9b924220 3291
ada8209b 3292 # print "scale $scale dad $dad zad $zad len $len\n";
58cde26e
JH
3293 # number bsstr len zad dad
3294 # 0.123 123e-3 3 0 3
3295 # 0.0123 123e-4 3 1 4
3296 # 0.001 1e-3 1 2 3
3297 # 1.23 123e-2 3 0 2
3298 # 1.2345 12345e-4 5 0 4
3299
3300 # do not round after/right of the $dad
3301 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
3302
ee15d750
JH
3303 # round to zero if rounding inside the $zad, but not for last zero like:
3304 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
3305 return $x->bzero() if $scale < $zad;
3306 if ($scale == $zad) # for 0.006, scale -3 and trunc
58cde26e 3307 {
b3abae2a 3308 $scale = -$len;
58cde26e
JH
3309 }
3310 else
3311 {
3312 # adjust round-point to be inside mantissa
3313 if ($zad != 0)
3314 {
3315 $scale = $scale-$zad;
3316 }
3317 else
3318 {
3319 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
3320 $scale = $dbd+$scale;
3321 }
3322 }
a0d0e21e 3323 }
58cde26e
JH
3324 else
3325 {
f9a08e12
JH
3326 # round left from the '.'
3327
58cde26e 3328 # 123 => 100 means length(123) = 3 - $scale (2) => 1
a5f75d66 3329
9b924220 3330 my $dbt = $MBI->_len($x->{_m});
b3abae2a 3331 # digits before dot
9b924220 3332 my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
b3abae2a
JH
3333 # should be the same, so treat it as this
3334 $scale = 1 if $scale == 0;
3335 # shortcut if already integer
3336 return $x if $scale == 1 && $dbt <= $dbd;
3337 # maximum digits before dot
3338 ++$dbd;
3339
3340 if ($scale > $dbd)
3341 {
3342 # not enough digits before dot, so round to zero
3343 return $x->bzero;
3344 }
3345 elsif ( $scale == $dbd )
3346 {
3347 # maximum
3348 $scale = -$dbt;
3349 }
58cde26e 3350 else
b3abae2a
JH
3351 {
3352 $scale = $dbd - $scale;
3353 }
a0d0e21e 3354 }
574bacfe 3355 # pass sign to bround for rounding modes '+inf' and '-inf'
ae161977 3356 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
9b924220
RGS
3357 $m->bround($scale,$mode);
3358 $x->{_m} = $m->{value}; # get our mantissa back
58cde26e
JH
3359 $x->bnorm();
3360 }
3361
3362sub bround
3363 {
3364 # accuracy: preserve $N digits, and overwrite the rest with 0's
ee15d750 3365 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
9b924220 3366
990fb837
RGS
3367 if (($_[0] || 0) < 0)
3368 {
3369 require Carp; Carp::croak ('bround() needs positive accuracy');
3370 }
58cde26e 3371
b68b7ab1
T
3372 my ($scale,$mode) = $x->_scale_a(@_);
3373 return $x if !defined $scale || $x->modify('bround'); # no-op
61f5c3f5 3374
ee15d750
JH
3375 # scale is now either $x->{_a}, $accuracy, or the user parameter
3376 # test whether $x already has lower accuracy, do nothing in this case
3377 # but do round if the accuracy is the same, since a math operation might
3378 # want to round a number with A=5 to 5 digits afterwards again
b68b7ab1 3379 return $x if defined $x->{_a} && $x->{_a} < $scale;
58cde26e 3380
61f5c3f5 3381 # scale < 0 makes no sense
b68b7ab1 3382 # scale == 0 => keep all digits
61f5c3f5 3383 # never round a +-inf, NaN
b68b7ab1 3384 return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
58cde26e 3385
b68b7ab1
T
3386 # 1: never round a 0
3387 # 2: if we should keep more digits than the mantissa has, do nothing
3388 if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
61f5c3f5
T
3389 {
3390 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
3391 return $x;
3392 }
f216259d 3393
58cde26e 3394 # pass sign to bround for '+inf' and '-inf' rounding modes
ae161977 3395 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
9b924220
RGS
3396
3397 $m->bround($scale,$mode); # round mantissa
3398 $x->{_m} = $m->{value}; # get our mantissa back
ee15d750 3399 $x->{_a} = $scale; # remember rounding
ef9466ea 3400 delete $x->{_p}; # and clear P
574bacfe 3401 $x->bnorm(); # del trailing zeros gen. by bround()
58cde26e
JH
3402 }
3403
3404sub bfloor
3405 {
3406 # return integer less or equal then $x
ee15d750 3407 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e
JH
3408
3409 return $x if $x->modify('bfloor');
3410
3411 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3412
3413 # if $x has digits after dot
9b924220 3414 if ($x->{_es} eq '-')
58cde26e 3415 {
9b924220
RGS
3416 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3417 $x->{_e} = $MBI->_zero(); # trunc/norm
3418 $x->{_es} = '+'; # abs e
3419 $MBI->_inc($x->{_m}) if $x->{sign} eq '-'; # increment if negative
f216259d 3420 }
61f5c3f5 3421 $x->round($a,$p,$r);
58cde26e 3422 }
288d023a 3423
58cde26e
JH
3424sub bceil
3425 {
3426 # return integer greater or equal then $x
ee15d750 3427 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e
JH
3428
3429 return $x if $x->modify('bceil');
3430 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3431
3432 # if $x has digits after dot
9b924220 3433 if ($x->{_es} eq '-')
58cde26e 3434 {
9b924220
RGS
3435 $x->{_m} = $MBI->_rsft($x->{_m},$x->{_e},10); # cut off digits after dot
3436 $x->{_e} = $MBI->_zero(); # trunc/norm
3437 $x->{_es} = '+'; # abs e
3438 $MBI->_inc($x->{_m}) if $x->{sign} eq '+'; # increment if positive
a0d0e21e 3439 }
61f5c3f5 3440 $x->round($a,$p,$r);
58cde26e
JH
3441 }
3442
394e6ffb
JH
3443sub brsft
3444 {
f9a08e12
JH
3445 # shift right by $y (divide by power of $n)
3446
3447 # set up parameters
3448 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3449 # objectify is costly, so avoid it
3450 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3451 {
3452 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3453 }
394e6ffb
JH
3454
3455 return $x if $x->modify('brsft');
3456 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3457
f9a08e12 3458 $n = 2 if !defined $n; $n = $self->new($n);
7b29e1e6
T
3459
3460 # negative amount?
3461 return $x->blsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3462
f603091d 3463 # the following call to bdiv() will return either quo or (quo,remainder):
f9a08e12 3464 $x->bdiv($n->bpow($y),$a,$p,$r,$y);
394e6ffb
JH
3465 }
3466
3467sub blsft
3468 {
f9a08e12
JH
3469 # shift left by $y (multiply by power of $n)
3470
3471 # set up parameters
3472 my ($self,$x,$y,$n,$a,$p,$r) = (ref($_[0]),@_);
3473 # objectify is costly, so avoid it
3474 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3475 {
3476 ($self,$x,$y,$n,$a,$p,$r) = objectify(2,@_);
3477 }
394e6ffb 3478
f9a08e12 3479 return $x if $x->modify('blsft');
394e6ffb
JH
3480 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3481
f9a08e12 3482 $n = 2 if !defined $n; $n = $self->new($n);
7b29e1e6
T
3483
3484 # negative amount?
3485 return $x->brsft($y->copy()->babs(),$n) if $y->{sign} =~ /^-/;
3486
f9a08e12 3487 $x->bmul($n->bpow($y),$a,$p,$r,$y);
394e6ffb
JH
3488 }
3489
58cde26e 3490###############################################################################
a5f75d66 3491
58cde26e
JH
3492sub DESTROY
3493 {
b282a552 3494 # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
58cde26e
JH
3495 }
3496
3497sub AUTOLOAD
3498 {
b3abae2a
JH
3499 # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
3500 # or falling back to MBI::bxxx()
58cde26e
JH
3501 my $name = $AUTOLOAD;
3502
ef9466ea
T
3503 $name =~ s/(.*):://; # split package
3504 my $c = $1 || $class;
ee15d750 3505 no strict 'refs';
ef9466ea 3506 $c->import() if $IMPORT == 0;
7b29e1e6 3507 if (!_method_alias($name))
58cde26e 3508 {
ee15d750
JH
3509 if (!defined $name)
3510 {
3511 # delayed load of Carp and avoid recursion
3512 require Carp;
ef9466ea 3513 Carp::croak ("$c: Can't call a method without name");
ee15d750 3514 }
7b29e1e6 3515 if (!_method_hand_up($name))
ee15d750
JH
3516 {
3517 # delayed load of Carp and avoid recursion
3518 require Carp;
ef9466ea 3519 Carp::croak ("Can't call $c\-\>$name, not a valid method");
ee15d750
JH
3520 }
3521 # try one level up, but subst. bxxx() for fxxx() since MBI only got bxxx()
3522 $name =~ s/^f/b/;
9b924220 3523 return &{"Math::BigInt"."::$name"}(@_);
a0d0e21e 3524 }
58cde26e 3525 my $bname = $name; $bname =~ s/^f/b/;
ef9466ea
T
3526 $c .= "::$name";
3527 *{$c} = \&{$bname};