This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Upgrade Math::BigInt from 1.999701 to 1.999704
[perl5.git] / cpan / 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
06ce15ad 15$VERSION = '1.999704';
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
06ce15ad
SH
863 # If called as $x -> blog() or $x -> blog(undef), don't objectify the
864 # undefined base, since undef signals that the base is Euler's number.
865 #unless (ref($x) && !defined($base)) {
866 # # objectify is costly, so avoid it
867 # if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) {
868 # ($self,$x,$base,$a,$p,$r) = objectify(2,@_);
869 # }
870 #}
871
50109ad0
RGS
872 return $x if $x->modify('blog');
873
06ce15ad 874 return $x -> bnan() if $x -> is_nan();
9393ace2 875
b3abae2a
JH
876 # we need to limit the accuracy to protect against overflow
877 my $fallback = 0;
990fb837
RGS
878 my ($scale,@params);
879 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
61f5c3f5 880
b3abae2a 881 # no rounding at all, so must use fallback
990fb837 882 if (scalar @params == 0)
b3abae2a
JH
883 {
884 # simulate old behaviour
990fb837
RGS
885 $params[0] = $self->div_scale(); # and round to it as accuracy
886 $params[1] = undef; # P = undef
887 $scale = $params[0]+4; # at least four more for proper round
888 $params[2] = $r; # round mode by caller or undef
b3abae2a
JH
889 $fallback = 1; # to clear a/p afterwards
890 }
891 else
892 {
893 # the 4 below is empirical, and there might be cases where it is not
894 # enough...
990fb837 895 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
b3abae2a 896 }
61f5c3f5 897
06ce15ad
SH
898 my $done = 0;
899 if (defined $base) {
900 $base = $self -> new($base) unless ref $base;
901 if ($base -> is_nan() || $base -> is_one()) {
902 $x -> bnan();
903 $done = 1;
904 } elsif ($base -> is_inf() || $base -> is_zero()) {
905 if ($x -> is_inf() || $x -> is_zero()) {
906 $x -> bnan();
907 } else {
908 $x -> bzero(@params);
909 }
910 $done = 1;
911 } elsif ($base -> is_negative()) { # -inf < base < 0
912 if ($x -> is_one()) { # x = 1
913 $x -> bzero(@params);
914 } elsif ($x == $base) {
915 $x -> bone('+', @params); # x = base
916 } else {
917 $x -> bnan(); # otherwise
918 }
919 $done = 1;
920 } elsif ($x == $base) {
921 $x -> bone('+', @params); # 0 < base && 0 < x < inf
922 $done = 1;
923 }
924 }
925
926 # We now know that the base is either undefined or positive and finite.
927
928 unless ($done) {
929 if ($x -> is_inf()) { # x = +/-inf
930 my $sign = defined $base && $base < 1 ? '-' : '+';
931 $x -> binf($sign);
932 $done = 1;
933 } elsif ($x -> is_neg()) { # -inf < x < 0
934 $x -> bnan();
935 $done = 1;
936 } elsif ($x -> is_one()) { # x = 1
937 $x -> bzero(@params);
938 $done = 1;
939 } elsif ($x -> is_zero()) { # x = 0
940 my $sign = defined $base && $base < 1 ? '+' : '-';
941 $x -> binf($sign);
942 $done = 1;
943 }
944 }
945
946 if ($done) {
947 if ($fallback) {
90d1b129 948 # clear a/p after round, since user did not request it
06ce15ad
SH
949 delete $x->{_a};
950 delete $x->{_p};
90d1b129
T
951 }
952 return $x;
953 }
61f5c3f5 954
b3abae2a 955 # when user set globals, they would interfere with our calculation, so
56d9de68 956 # disable them and later re-enable them
b3abae2a
JH
957 no strict 'refs';
958 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
959 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
960 # we also need to disable any set A or P on $x (_find_round_parameters took
961 # them already into account), since these would interfere, too
962 delete $x->{_a}; delete $x->{_p};
9393ace2 963 # need to disable $upgrade in BigInt, to avoid deep recursion
b3abae2a 964 local $Math::BigInt::upgrade = undef;
93c87d9d 965 local $Math::BigFloat::downgrade = undef;
990fb837
RGS
966
967 # upgrade $x if $x is not a BigFloat (handle BigInt input)
7d193e39 968 # XXX TODO: rebless!
990fb837
RGS
969 if (!$x->isa('Math::BigFloat'))
970 {
971 $x = Math::BigFloat->new($x);
972 $self = ref($x);
973 }
b282a552 974
06ce15ad 975 $done = 0;
b282a552
T
976
977 # If the base is defined and an integer, try to calculate integer result
978 # first. This is very fast, and in case the real result was found, we can
979 # stop right here.
980 if (defined $base && $base->is_int() && $x->is_int())
981 {
9b924220
RGS
982 my $i = $MBI->_copy( $x->{_m} );
983 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
984 my $int = Math::BigInt->bzero();
985 $int->{value} = $i;
b282a552
T
986 $int->blog($base->as_number());
987 # if ($exact)
9b924220 988 if ($base->as_number()->bpow($int) == $x)
b282a552
T
989 {
990 # found result, return it
9b924220
RGS
991 $x->{_m} = $int->{value};
992 $x->{_e} = $MBI->_zero();
993 $x->{_es} = '+';
b282a552
T
994 $x->bnorm();
995 $done = 1;
996 }
997 }
998
999 if ($done == 0)
9393ace2 1000 {
7d193e39
T
1001 # base is undef, so base should be e (Euler's number), so first calculate the
1002 # log to base e (using reduction by 10 (and probably 2)):
b282a552 1003 $self->_log_10($x,$scale);
9b924220 1004
b282a552
T
1005 # and if a different base was requested, convert it
1006 if (defined $base)
1007 {
1008 $base = Math::BigFloat->new($base) unless $base->isa('Math::BigFloat');
1009 # not ln, but some other base (don't modify $base)
1010 $x->bdiv( $base->copy()->blog(undef,$scale), $scale );
1011 }
9393ace2 1012 }
990fb837 1013
091c87b1 1014 # shortcut to not run through _find_round_parameters again
990fb837 1015 if (defined $params[0])
b3abae2a 1016 {
990fb837 1017 $x->bround($params[0],$params[2]); # then round accordingly
b3abae2a
JH
1018 }
1019 else
1020 {
990fb837 1021 $x->bfround($params[1],$params[2]); # then round accordingly
b3abae2a
JH
1022 }
1023 if ($fallback)
1024 {
1025 # clear a/p after round, since user did not request it
ef9466ea 1026 delete $x->{_a}; delete $x->{_p};
b3abae2a
JH
1027 }
1028 # restore globals
1029 $$abr = $ab; $$pbr = $pb;
1030
1031 $x;
61f5c3f5
T
1032 }
1033
50109ad0
RGS
1034sub _len_to_steps
1035 {
1036 # Given D (digits in decimal), compute N so that N! (N factorial) is
1037 # at least D digits long. D should be at least 50.
1038 my $d = shift;
1039
1040 # two constants for the Ramanujan estimate of ln(N!)
1041 my $lg2 = log(2 * 3.14159265) / 2;
1042 my $lg10 = log(10);
1043
1044 # D = 50 => N => 42, so L = 40 and R = 50
1045 my $l = 40; my $r = $d;
1046
1047 # Otherwise this does not work under -Mbignum and we do not yet have "no bignum;" :(
1048 $l = $l->numify if ref($l);
1049 $r = $r->numify if ref($r);
1050 $lg2 = $lg2->numify if ref($lg2);
1051 $lg10 = $lg10->numify if ref($lg10);
1052
1053 # binary search for the right value (could this be written as the reverse of lg(n!)?)
1054 while ($r - $l > 1)
1055 {
1056 my $n = int(($r - $l) / 2) + $l;
1057 my $ramanujan =
1058 int(($n * log($n) - $n + log( $n * (1 + 4*$n*(1+2*$n)) ) / 6 + $lg2) / $lg10);
1059 $ramanujan > $d ? $r = $n : $l = $n;
1060 }
1061 $l;
1062 }
1063
1064sub bnok
1065 {
1066 # Calculate n over k (binomial coefficient or "choose" function) as integer.
1067 # set up parameters
1068 my ($self,$x,$y,@r) = (ref($_[0]),@_);
1069
1070 # objectify is costly, so avoid it
1071 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1072 {
1073 ($self,$x,$y,@r) = objectify(2,@_);
1074 }
1075
1076 return $x if $x->modify('bnok');
1077
1078 return $x->bnan() if $x->is_nan() || $y->is_nan();
1079 return $x->binf() if $x->is_inf();
1080
1081 my $u = $x->as_int();
1082 $u->bnok($y->as_int());
1083
1084 $x->{_m} = $u->{value};
1085 $x->{_e} = $MBI->_zero();
1086 $x->{_es} = '+';
1087 $x->{sign} = '+';
1088 $x->bnorm(@r);
1089 }
1090
7d193e39
T
1091sub bexp
1092 {
50109ad0 1093 # Calculate e ** X (Euler's number to the power of X)
7d193e39
T
1094 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
1095
50109ad0
RGS
1096 return $x if $x->modify('bexp');
1097
7d193e39
T
1098 return $x->binf() if $x->{sign} eq '+inf';
1099 return $x->bzero() if $x->{sign} eq '-inf';
1100
1101 # we need to limit the accuracy to protect against overflow
1102 my $fallback = 0;
1103 my ($scale,@params);
1104 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
1105
1106 # also takes care of the "error in _find_round_parameters?" case
1107 return $x if $x->{sign} eq 'NaN';
1108
1109 # no rounding at all, so must use fallback
1110 if (scalar @params == 0)
1111 {
1112 # simulate old behaviour
1113 $params[0] = $self->div_scale(); # and round to it as accuracy
1114 $params[1] = undef; # P = undef
1115 $scale = $params[0]+4; # at least four more for proper round
1116 $params[2] = $r; # round mode by caller or undef
1117 $fallback = 1; # to clear a/p afterwards
1118 }
1119 else
1120 {
1121 # the 4 below is empirical, and there might be cases where it's not enough...
1122 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
1123 }
1124
1125 return $x->bone(@params) if $x->is_zero();
1126
1127 if (!$x->isa('Math::BigFloat'))
1128 {
1129 $x = Math::BigFloat->new($x);
1130 $self = ref($x);
1131 }
1132
1133 # when user set globals, they would interfere with our calculation, so
1134 # disable them and later re-enable them
1135 no strict 'refs';
1136 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
1137 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
1138 # we also need to disable any set A or P on $x (_find_round_parameters took
1139 # them already into account), since these would interfere, too
1140 delete $x->{_a}; delete $x->{_p};
1141 # need to disable $upgrade in BigInt, to avoid deep recursion
1142 local $Math::BigInt::upgrade = undef;
1143 local $Math::BigFloat::downgrade = undef;
1144
1145 my $x_org = $x->copy();
7d193e39
T
1146
1147 # We use the following Taylor series:
1148
1149 # x x^2 x^3 x^4
1150 # e = 1 + --- + --- + --- + --- ...
1151 # 1! 2! 3! 4!
1152
1153 # The difference for each term is X and N, which would result in:
1154 # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term
1155
1156 # But it is faster to compute exp(1) and then raising it to the
1157 # given power, esp. if $x is really big and an integer because:
1158
1159 # * The numerator is always 1, making the computation faster
1160 # * the series converges faster in the case of x == 1
1161 # * We can also easily check when we have reached our limit: when the
1162 # term to be added is smaller than "1E$scale", we can stop - f.i.
1163 # scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
1164 # * we can compute the *exact* result by simulating bigrat math:
1165
1166 # 1 1 gcd(3,4) = 1 1*24 + 1*6 5
1167 # - + - = ---------- = --
1168 # 6 24 6*24 24
1169
1170 # We do not compute the gcd() here, but simple do:
1171 # 1 1 1*24 + 1*6 30
1172 # - + - = --------- = --
1173 # 6 24 6*24 144
1174
1175 # In general:
1176 # a c a*d + c*b and note that c is always 1 and d = (b*f)
1177 # - + - = ---------
1178 # b d b*d
1179
1180 # This leads to: which can be reduced by b to:
1181 # a 1 a*b*f + b a*f + 1
1182 # - + - = --------- = -------
1183 # b b*f b*b*f b*f
1184
1185 # The first terms in the series are:
1186
1187 # 1 1 1 1 1 1 1 1 13700
1188 # -- + -- + -- + -- + -- + --- + --- + ---- = -----
1189 # 1 1 2 6 24 120 720 5040 5040
1190
50109ad0 1191 # Note that we cannot simple reduce 13700/5040 to 685/252, but must keep A and B!
7d193e39 1192
50109ad0 1193 if ($scale <= 75)
7d193e39 1194 {
50109ad0
RGS
1195 # set $x directly from a cached string form
1196 $x->{_m} = $MBI->_new(
1197 "27182818284590452353602874713526624977572470936999595749669676277240766303535476");
1198 $x->{sign} = '+';
1199 $x->{_es} = '-';
1200 $x->{_e} = $MBI->_new(79);
7d193e39 1201 }
50109ad0
RGS
1202 else
1203 {
1204 # compute A and B so that e = A / B.
1205
1206 # After some terms we end up with this, so we use it as a starting point:
1207 my $A = $MBI->_new("90933395208605785401971970164779391644753259799242");
1208 my $F = $MBI->_new(42); my $step = 42;
1209
1210 # Compute how many steps we need to take to get $A and $B sufficiently big
1211 my $steps = _len_to_steps($scale - 4);
1212# print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
1213 while ($step++ <= $steps)
1214 {
1215 # calculate $a * $f + 1
1216 $A = $MBI->_mul($A, $F);
1217 $A = $MBI->_inc($A);
1218 # increment f
1219 $F = $MBI->_inc($F);
1220 }
1221 # compute $B as factorial of $steps (this is faster than doing it manually)
1222 my $B = $MBI->_fac($MBI->_new($steps));
1223
1224# print "A ", $MBI->_str($A), "\nB ", $MBI->_str($B), "\n";
7d193e39 1225
50109ad0
RGS
1226 # compute A/B with $scale digits in the result (truncate, not round)
1227 $A = $MBI->_lsft( $A, $MBI->_new($scale), 10);
1228 $A = $MBI->_div( $A, $B );
7d193e39 1229
50109ad0
RGS
1230 $x->{_m} = $A;
1231 $x->{sign} = '+';
1232 $x->{_es} = '-';
1233 $x->{_e} = $MBI->_new($scale);
1234 }
7d193e39 1235
50109ad0
RGS
1236 # $x contains now an estimate of e, with some surplus digits, so we can round
1237 if (!$x_org->is_one())
7d193e39 1238 {
50109ad0
RGS
1239 # raise $x to the wanted power and round it in one step:
1240 $x->bpow($x_org, @params);
7d193e39
T
1241 }
1242 else
1243 {
50109ad0
RGS
1244 # else just round the already computed result
1245 delete $x->{_a}; delete $x->{_p};
1246 # shortcut to not run through _find_round_parameters again
1247 if (defined $params[0])
1248 {
1249 $x->bround($params[0],$params[2]); # then round accordingly
1250 }
1251 else
1252 {
1253 $x->bfround($params[1],$params[2]); # then round accordingly
1254 }
7d193e39
T
1255 }
1256 if ($fallback)
1257 {
1258 # clear a/p after round, since user did not request it
1259 delete $x->{_a}; delete $x->{_p};
1260 }
1261 # restore globals
1262 $$abr = $ab; $$pbr = $pb;
1263
1264 $x; # return modified $x
1265 }
1266
990fb837
RGS
1267sub _log
1268 {
091c87b1 1269 # internal log function to calculate ln() based on Taylor series.
990fb837
RGS
1270 # Modifies $x in place.
1271 my ($self,$x,$scale) = @_;
1272
091c87b1
T
1273 # in case of $x == 1, result is 0
1274 return $x->bzero() if $x->is_one();
1275
9681bfa6 1276 # XXX TODO: rewrite this in a similar manner to bexp()
7d193e39 1277
990fb837
RGS
1278 # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
1279
1280 # u = x-1, v = x+1
1281 # _ _
1282 # Taylor: | u 1 u^3 1 u^5 |
1283 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
1284 # |_ v 3 v^3 5 v^5 _|
1285
1286 # This takes much more steps to calculate the result and is thus not used
1287 # u = x-1
1288 # _ _
1289 # Taylor: | u 1 u^2 1 u^3 |
1290 # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
1291 # |_ x 2 x^2 3 x^3 _|
1292
990fb837
RGS
1293 my ($limit,$v,$u,$below,$factor,$two,$next,$over,$f);
1294
1295 $v = $x->copy(); $v->binc(); # v = x+1
1296 $x->bdec(); $u = $x->copy(); # u = x-1; x = x-1
1297 $x->bdiv($v,$scale); # first term: u/v
1298 $below = $v->copy();
1299 $over = $u->copy();
1300 $u *= $u; $v *= $v; # u^2, v^2
1301 $below->bmul($v); # u^3, v^3
1302 $over->bmul($u);
1303 $factor = $self->new(3); $f = $self->new(2);
1304
1305 my $steps = 0 if DEBUG;
1306 $limit = $self->new("1E-". ($scale-1));
1307 while (3 < 5)
1308 {
1309 # we calculate the next term, and add it to the last
1310 # when the next term is below our limit, it won't affect the outcome
1311 # anymore, so we stop
1312
1313 # calculating the next term simple from over/below will result in quite
1314 # a time hog if the input has many digits, since over and below will
1315 # accumulate more and more digits, and the result will also have many
1316 # digits, but in the end it is rounded to $scale digits anyway. So if we
1317 # round $over and $below first, we save a lot of time for the division
1318 # (not with log(1.2345), but try log (123**123) to see what I mean. This
1319 # can introduce a rounding error if the division result would be f.i.
1320 # 0.1234500000001 and we round it to 5 digits it would become 0.12346, but
091c87b1
T
1321 # if we truncated $over and $below we might get 0.12345. Does this matter
1322 # for the end result? So we give $over and $below 4 more digits to be
1323 # on the safe side (unscientific error handling as usual... :+D
7d193e39 1324
990fb837
RGS
1325 $next = $over->copy->bround($scale+4)->bdiv(
1326 $below->copy->bmul($factor)->bround($scale+4),
1327 $scale);
1328
1329## old version:
1330## $next = $over->copy()->bdiv($below->copy()->bmul($factor),$scale);
1331
1332 last if $next->bacmp($limit) <= 0;
1333
1334 delete $next->{_a}; delete $next->{_p};
1335 $x->badd($next);
990fb837
RGS
1336 # calculate things for the next term
1337 $over *= $u; $below *= $v; $factor->badd($f);
1338 if (DEBUG)
1339 {
1340 $steps++; print "step $steps = $x\n" if $steps % 10 == 0;
1341 }
1342 }
990fb837 1343 print "took $steps steps\n" if DEBUG;
7d193e39 1344 $x->bmul($f); # $x *= 2
990fb837
RGS
1345 }
1346
1347sub _log_10
1348 {
091c87b1
T
1349 # Internal log function based on reducing input to the range of 0.1 .. 9.99
1350 # and then "correcting" the result to the proper one. Modifies $x in place.
990fb837
RGS
1351 my ($self,$x,$scale) = @_;
1352
7d193e39 1353 # Taking blog() from numbers greater than 10 takes a *very long* time, so we
990fb837 1354 # break the computation down into parts based on the observation that:
7d193e39
T
1355 # blog(X*Y) = blog(X) + blog(Y)
1356 # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
1357 # $x is the faster it gets. Since 2*$x takes about 10 times as
1358 # long, we make it faster by about a factor of 100 by dividing $x by 10.
1359
1360 # The same observation is valid for numbers smaller than 0.1, e.g. computing
1361 # log(1) is fastest, and the further away we get from 1, the longer it takes.
1362 # So we also 'break' this down by multiplying $x with 10 and subtract the
990fb837
RGS
1363 # log(10) afterwards to get the correct result.
1364
7d193e39
T
1365 # To get $x even closer to 1, we also divide by 2 and then use log(2) to
1366 # correct for this. For instance if $x is 2.4, we use the formula:
1367 # blog(2.4 * 2) == blog (1.2) + blog(2)
1368 # and thus calculate only blog(1.2) and blog(2), which is faster in total
1369 # than calculating blog(2.4).
1370
1371 # In addition, the values for blog(2) and blog(10) are cached.
1372
1373 # Calculate nr of digits before dot:
9b924220
RGS
1374 my $dbd = $MBI->_num($x->{_e});
1375 $dbd = -$dbd if $x->{_es} eq '-';
1376 $dbd += $MBI->_len($x->{_m});
990fb837
RGS
1377
1378 # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
1379 # infinite recursion
1380
1381 my $calc = 1; # do some calculation?
1382
1383 # disable the shortcut for 10, since we need log(10) and this would recurse
1384 # infinitely deep
9b924220 1385 if ($x->{_es} eq '+' && $MBI->_is_one($x->{_e}) && $MBI->_is_one($x->{_m}))
990fb837
RGS
1386 {
1387 $dbd = 0; # disable shortcut
1388 # we can use the cached value in these cases
1389 if ($scale <= $LOG_10_A)
1390 {
7d193e39 1391 $x->bzero(); $x->badd($LOG_10); # modify $x in place
990fb837
RGS
1392 $calc = 0; # no need to calc, but round
1393 }
7d193e39 1394 # if we can't use the shortcut, we continue normally
990fb837 1395 }
091c87b1 1396 else
990fb837 1397 {
091c87b1 1398 # disable the shortcut for 2, since we maybe have it cached
9b924220 1399 if (($MBI->_is_zero($x->{_e}) && $MBI->_is_two($x->{_m})))
990fb837 1400 {
091c87b1
T
1401 $dbd = 0; # disable shortcut
1402 # we can use the cached value in these cases
1403 if ($scale <= $LOG_2_A)
1404 {
7d193e39 1405 $x->bzero(); $x->badd($LOG_2); # modify $x in place
091c87b1
T
1406 $calc = 0; # no need to calc, but round
1407 }
7d193e39 1408 # if we can't use the shortcut, we continue normally
990fb837
RGS
1409 }
1410 }
1411
1412 # if $x = 0.1, we know the result must be 0-log(10)
9b924220
RGS
1413 if ($calc != 0 && $x->{_es} eq '-' && $MBI->_is_one($x->{_e}) &&
1414 $MBI->_is_one($x->{_m}))
990fb837
RGS
1415 {
1416 $dbd = 0; # disable shortcut
1417 # we can use the cached value in these cases
1418 if ($scale <= $LOG_10_A)
1419 {
1420 $x->bzero(); $x->bsub($LOG_10);
1421 $calc = 0; # no need to calc, but round
1422 }
1423 }
1424
091c87b1
T
1425 return if $calc == 0; # already have the result
1426
990fb837
RGS
1427 # default: these correction factors are undef and thus not used
1428 my $l_10; # value of ln(10) to A of $scale
1429 my $l_2; # value of ln(2) to A of $scale
1430
7d193e39
T
1431 my $two = $self->new(2);
1432
990fb837
RGS
1433 # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
1434 # so don't do this shortcut for 1 or 0
1435 if (($dbd > 1) || ($dbd < 0))
1436 {
1437 # convert our cached value to an object if not already (avoid doing this
1438 # at import() time, since not everybody needs this)
1439 $LOG_10 = $self->new($LOG_10,undef,undef) unless ref $LOG_10;
1440
1441 #print "x = $x, dbd = $dbd, calc = $calc\n";
1442 # got more than one digit before the dot, or more than one zero after the
1443 # dot, so do:
1444 # log(123) == log(1.23) + log(10) * 2
1445 # log(0.0123) == log(1.23) - log(10) * 2
1446
1447 if ($scale <= $LOG_10_A)
1448 {
1449 # use cached value
990fb837
RGS
1450 $l_10 = $LOG_10->copy(); # copy for mul
1451 }
1452 else
1453 {
7d193e39 1454 # else: slower, compute and cache result
990fb837
RGS
1455 # also disable downgrade for this code path
1456 local $Math::BigFloat::downgrade = undef;
7d193e39
T
1457
1458 # shorten the time to calculate log(10) based on the following:
1459 # log(1.25 * 8) = log(1.25) + log(8)
1460 # = log(1.25) + log(2) + log(2) + log(2)
1461
1462 # first get $l_2 (and possible compute and cache log(2))
1463 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1464 if ($scale <= $LOG_2_A)
1465 {
1466 # use cached value
1467 $l_2 = $LOG_2->copy(); # copy() for the mul below
1468 }
1469 else
1470 {
1471 # else: slower, compute and cache result
1472 $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1473 $LOG_2 = $l_2->copy(); # cache the result for later
1474 # the copy() is for mul below
1475 $LOG_2_A = $scale;
1476 }
1477
1478 # now calculate log(1.25):
1479 $l_10 = $self->new('1.25'); $self->_log($l_10, $scale); # scale+4, actually
1480
1481 # log(1.25) + log(2) + log(2) + log(2):
1482 $l_10->badd($l_2);
1483 $l_10->badd($l_2);
1484 $l_10->badd($l_2);
1485 $LOG_10 = $l_10->copy(); # cache the result for later
1486 # the copy() is for mul below
1487 $LOG_10_A = $scale;
990fb837
RGS
1488 }
1489 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
9b924220
RGS
1490 $l_10->bmul( $self->new($dbd)); # log(10) * (digits_before_dot-1)
1491 my $dbd_sign = '+';
1492 if ($dbd < 0)
1493 {
1494 $dbd = -$dbd;
1495 $dbd_sign = '-';
1496 }
1497 ($x->{_e}, $x->{_es}) =
1498 _e_sub( $x->{_e}, $MBI->_new($dbd), $x->{_es}, $dbd_sign); # 123 => 1.23
990fb837
RGS
1499
1500 }
1501
1502 # Now: 0.1 <= $x < 10 (and possible correction in l_10)
1503
1504 ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
1505 ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
1506
27e7b8bb
T
1507 $HALF = $self->new($HALF) unless ref($HALF);
1508
091c87b1 1509 my $twos = 0; # default: none (0 times)
7d193e39 1510 while ($x->bacmp($HALF) <= 0) # X <= 0.5
990fb837 1511 {
091c87b1
T
1512 $twos--; $x->bmul($two);
1513 }
7d193e39 1514 while ($x->bacmp($two) >= 0) # X >= 2
091c87b1
T
1515 {
1516 $twos++; $x->bdiv($two,$scale+4); # keep all digits
1517 }
0d2463eb 1518 $x->bround($scale+4);
7d193e39
T
1519 # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
1520 # So calculate correction factor based on ln(2):
091c87b1
T
1521 if ($twos != 0)
1522 {
1523 $LOG_2 = $self->new($LOG_2,undef,undef) unless ref $LOG_2;
1524 if ($scale <= $LOG_2_A)
990fb837 1525 {
091c87b1 1526 # use cached value
7d193e39 1527 $l_2 = $LOG_2->copy(); # copy() for the mul below
990fb837 1528 }
091c87b1 1529 else
990fb837 1530 {
7d193e39 1531 # else: slower, compute and cache result
091c87b1
T
1532 # also disable downgrade for this code path
1533 local $Math::BigFloat::downgrade = undef;
7d193e39
T
1534 $l_2 = $two->copy(); $self->_log($l_2, $scale); # scale+4, actually
1535 $LOG_2 = $l_2->copy(); # cache the result for later
1536 # the copy() is for mul below
1537 $LOG_2_A = $scale;
990fb837 1538 }
091c87b1 1539 $l_2->bmul($twos); # * -2 => subtract, * 2 => add
990fb837 1540 }
e16b8d1b 1541 else
1542 {
1543 undef $l_2;
1544 }
990fb837 1545
091c87b1
T
1546 $self->_log($x,$scale); # need to do the "normal" way
1547 $x->badd($l_10) if defined $l_10; # correct it by ln(10)
1548 $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
7d193e39 1549
990fb837 1550 # all done, $x contains now the result
7d193e39 1551 $x;
990fb837
RGS
1552 }
1553
58cde26e
JH
1554sub blcm
1555 {
ee15d750 1556 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
58cde26e
JH
1557 # does not modify arguments, but returns new object
1558 # Lowest Common Multiplicator
58cde26e
JH
1559
1560 my ($self,@arg) = objectify(0,@_);
1561 my $x = $self->new(shift @arg);
b68b7ab1 1562 while (@arg) { $x = Math::BigInt::__lcm($x,shift @arg); }
58cde26e
JH
1563 $x;
1564 }
1565
b68b7ab1
T
1566sub bgcd
1567 {
1568 # (BINT or num_str, BINT or num_str) return BINT
58cde26e 1569 # does not modify arguments, but returns new object
b68b7ab1
T
1570
1571 my $y = shift;
1572 $y = __PACKAGE__->new($y) if !ref($y);
1573 my $self = ref($y);
1574 my $x = $y->copy()->babs(); # keep arguments
1575
1576 return $x->bnan() if $x->{sign} !~ /^[+-]$/ # x NaN?
1577 || !$x->is_int(); # only for integers now
1578
1579 while (@_)
1580 {
1581 my $t = shift; $t = $self->new($t) if !ref($t);
1582 $y = $t->copy()->babs();
1583
1584 return $x->bnan() if $y->{sign} !~ /^[+-]$/ # y NaN?
1585 || !$y->is_int(); # only for integers now
1586
1587 # greatest common divisor
1588 while (! $y->is_zero())
1589 {
1590 ($x,$y) = ($y->copy(), $x->copy()->bmod($y));
1591 }
1592
1593 last if $x->is_one();
1594 }
58cde26e
JH
1595 $x;
1596 }
1597
9b924220 1598##############################################################################
b3abae2a 1599
9b924220 1600sub _e_add
091c87b1 1601 {
9b924220
RGS
1602 # Internal helper sub to take two positive integers and their signs and
1603 # then add them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1604 # output ($CALC,('+'|'-'))
1605 my ($x,$y,$xs,$ys) = @_;
1606
1607 # if the signs are equal we can add them (-5 + -3 => -(5 + 3) => -8)
1608 if ($xs eq $ys)
1609 {
1610 $x = $MBI->_add ($x, $y ); # a+b
1611 # the sign follows $xs
1612 return ($x, $xs);
1613 }
091c87b1 1614
9b924220
RGS
1615 my $a = $MBI->_acmp($x,$y);
1616 if ($a > 0)
1617 {
1618 $x = $MBI->_sub ($x , $y); # abs sub
1619 }
1620 elsif ($a == 0)
1621 {
1622 $x = $MBI->_zero(); # result is 0
1623 $xs = '+';
1624 }
1625 else # a < 0
1626 {
1627 $x = $MBI->_sub ( $y, $x, 1 ); # abs sub
1628 $xs = $ys;
1629 }
1630 ($x,$xs);
091c87b1
T
1631 }
1632
9b924220
RGS
1633sub _e_sub
1634 {
1635 # Internal helper sub to take two positive integers and their signs and
1636 # then subtract them. Input ($CALC,$CALC,('+'|'-'),('+'|'-')),
1637 # output ($CALC,('+'|'-'))
1638 my ($x,$y,$xs,$ys) = @_;
1639
1640 # flip sign
1641 $ys =~ tr/+-/-+/;
1642 _e_add($x,$y,$xs,$ys); # call add (does subtract now)
1643 }
1644
1645###############################################################################
1646# is_foo methods (is_negative, is_positive are inherited from BigInt)
1647
b3abae2a
JH
1648sub is_int
1649 {
1650 # return true if arg (BFLOAT or num_str) is an integer
091c87b1 1651 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
b3abae2a 1652
80365507
T
1653 (($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1654 ($x->{_es} eq '+')) ? 1 : 0; # 1e-1 => no integer
b3abae2a
JH
1655 }
1656
58cde26e
JH
1657sub is_zero
1658 {
b3abae2a 1659 # return true if arg (BFLOAT or num_str) is zero
091c87b1 1660 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
574bacfe 1661
80365507 1662 ($x->{sign} eq '+' && $MBI->_is_zero($x->{_m})) ? 1 : 0;
58cde26e
JH
1663 }
1664
1665sub is_one
1666 {
b3abae2a 1667 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
091c87b1 1668 my ($self,$x,$sign) = ref($_[0]) ? (undef,@_) : objectify(1,@_);
ee15d750 1669
990fb837 1670 $sign = '+' if !defined $sign || $sign ne '-';
80365507
T
1671
1672 ($x->{sign} eq $sign &&
1673 $MBI->_is_zero($x->{_e}) &&
1674 $MBI->_is_one($x->{_m}) ) ? 1 : 0;
58cde26e
JH
1675 }
1676
1677sub is_odd
1678 {
ee15d750 1679 # return true if arg (BFLOAT or num_str) is odd or false if even
091c87b1 1680 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
0716bf9b 1681
80365507
T
1682 (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1683 ($MBI->_is_zero($x->{_e})) &&
1684 ($MBI->_is_odd($x->{_m}))) ? 1 : 0;
58cde26e
JH
1685 }
1686
1687sub is_even
1688 {
b22b3e31 1689 # return true if arg (BINT or num_str) is even or false if odd
091c87b1 1690 my ($self,$x) = ref($_[0]) ? (undef,$_[0]) : objectify(1,@_);
0716bf9b 1691
80365507
T
1692 (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1693 ($x->{_es} eq '+') && # 123.45 isn't
1694 ($MBI->_is_even($x->{_m}))) ? 1 : 0; # but 1200 is
58cde26e
JH
1695 }
1696
80365507 1697sub bmul
58cde26e 1698 {
80365507 1699 # multiply two numbers
f9a08e12
JH
1700
1701 # set up parameters
80365507 1702 my ($self,$x,$y,@r) = (ref($_[0]),@_);
f9a08e12
JH
1703 # objectify is costly, so avoid it
1704 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1705 {
80365507 1706 ($self,$x,$y,@r) = objectify(2,@_);
f9a08e12 1707 }
58cde26e 1708
50109ad0
RGS
1709 return $x if $x->modify('bmul');
1710
58cde26e
JH
1711 return $x->bnan() if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1712
574bacfe
JH
1713 # inf handling
1714 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1715 {
13a12e00 1716 return $x->bnan() if $x->is_zero() || $y->is_zero();
574bacfe
JH
1717 # result will always be +-inf:
1718 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1719 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1720 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1721 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1722 return $x->binf('-');
1723 }
8f675a64 1724
80365507
T
1725 return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
1726 ((!$x->isa($self)) || (!$y->isa($self)));
1727
1728 # aEb * cEd = (a*c)E(b+d)
1729 $MBI->_mul($x->{_m},$y->{_m});
1730 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1731
1732 $r[3] = $y; # no push!
1733
1734 # adjust sign:
1735 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1736 $x->bnorm->round(@r);
1737 }
1738
1739sub bmuladd
1740 {
1741 # multiply two numbers and add the third to the result
1742
1743 # set up parameters
913a64d5 1744 my ($self,$x,$y,$z,@r) = objectify(3,@_);
80365507
T
1745
1746 return $x if $x->modify('bmuladd');
1747
1748 return $x->bnan() if (($x->{sign} eq $nan) ||
1749 ($y->{sign} eq $nan) ||
1750 ($z->{sign} eq $nan));
1751
1752 # inf handling
1753 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/))
1754 {
1755 return $x->bnan() if $x->is_zero() || $y->is_zero();
1756 # result will always be +-inf:
1757 # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1758 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1759 return $x->binf() if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1760 return $x->binf() if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1761 return $x->binf('-');
1762 }
1763
1764 return $upgrade->bmul($x,$y,@r) if defined $upgrade &&
8f675a64 1765 ((!$x->isa($self)) || (!$y->isa($self)));
574bacfe 1766
58cde26e 1767 # aEb * cEd = (a*c)E(b+d)
9b924220
RGS
1768 $MBI->_mul($x->{_m},$y->{_m});
1769 ($x->{_e}, $x->{_es}) = _e_add($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1770
80365507
T
1771 $r[3] = $y; # no push!
1772
58cde26e
JH
1773 # adjust sign:
1774 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
80365507
T
1775
1776 # z=inf handling (z=NaN handled above)
1777 $x->{sign} = $z->{sign}, return $x if $z->{sign} =~ /^[+-]inf$/;
1778
1779 # take lower of the two e's and adapt m1 to it to match m2
1780 my $e = $z->{_e};
1781 $e = $MBI->_zero() if !defined $e; # if no BFLOAT?
1782 $e = $MBI->_copy($e); # make copy (didn't do it yet)
1783
1784 my $es;
1785
1786 ($e,$es) = _e_sub($e, $x->{_e}, $z->{_es} || '+', $x->{_es});
1787
1788 my $add = $MBI->_copy($z->{_m});
1789
1790 if ($es eq '-') # < 0
1791 {
1792 $MBI->_lsft( $x->{_m}, $e, 10);
1793 ($x->{_e},$x->{_es}) = _e_add($x->{_e}, $e, $x->{_es}, $es);
1794 }
1795 elsif (!$MBI->_is_zero($e)) # > 0
1796 {
1797 $MBI->_lsft($add, $e, 10);
1798 }
1799 # else: both e are the same, so just leave them
1800
1801 if ($x->{sign} eq $z->{sign})
1802 {
1803 # add
1804 $x->{_m} = $MBI->_add($x->{_m}, $add);
1805 }
1806 else
1807 {
1808 ($x->{_m}, $x->{sign}) =
1809 _e_add($x->{_m}, $add, $x->{sign}, $z->{sign});
1810 }
1811
1812 # delete trailing zeros, then round
1813 $x->bnorm()->round(@r);
58cde26e
JH
1814 }
1815
1816sub bdiv
1817 {
1818 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
b4a10d33 1819 # (BFLOAT, BFLOAT) (quo, rem) or BFLOAT (only quo)
f9a08e12
JH
1820
1821 # set up parameters
1822 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
1823 # objectify is costly, so avoid it
1824 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
1825 {
1826 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
1827 }
58cde26e 1828
50109ad0
RGS
1829 return $x if $x->modify('bdiv');
1830
b4a10d33 1831 my $wantarray = wantarray; # call only once
574bacfe 1832
b4a10d33
PJA
1833 # At least one argument is NaN. This is handled the same way as in
1834 # Math::BigInt -> bdiv().
1835
1836 if ($x -> is_nan() || $y -> is_nan()) {
1837 return $wantarray ? ($x -> bnan(), $self -> bnan()) : $x -> bnan();
1838 }
1839
1840 # Divide by zero and modulo zero. This is handled the same way as in
1841 # Math::BigInt -> bdiv(). See the comment in the code for Math::BigInt ->
1842 # bdiv() for further details.
1843
1844 if ($y -> is_zero()) {
1845 my ($quo, $rem);
1846 if ($wantarray) {
1847 $rem = $x -> copy();
1848 }
1849 if ($x -> is_zero()) {
1850 $quo = $x -> bnan();
1851 } else {
1852 $quo = $x -> binf($x -> {sign});
1853 }
1854 return $wantarray ? ($quo, $rem) : $quo;
1855 }
1856
1857 # Numerator (dividend) is +/-inf. This is handled the same way as in
1858 # Math::BigInt -> bdiv(). See the comment in the code for Math::BigInt ->
1859 # bdiv() for further details.
1860
1861 if ($x -> is_inf()) {
1862 my ($quo, $rem);
1863 $rem = $self -> bnan() if $wantarray;
1864 if ($y -> is_inf()) {
1865 $quo = $x -> bnan();
1866 } else {
1867 my $sign = $x -> bcmp(0) == $y -> bcmp(0) ? '+' : '-';
1868 $quo = $x -> binf($sign);
1869 }
1870 return $wantarray ? ($quo, $rem) : $quo;
1871 }
1872
1873 # Denominator (divisor) is +/-inf. This is handled the same way as in
1874 # Math::BigInt -> bdiv(), with one exception: In scalar context,
1875 # Math::BigFloat does true division (although rounded), not floored division
1876 # (F-division), so a finite number divided by +/-inf is always zero. See the
1877 # comment in the code for Math::BigInt -> bdiv() for further details.
1878
1879 if ($y -> is_inf()) {
1880 my ($quo, $rem);
1881 if ($wantarray) {
1882 if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) {
1883 $rem = $x -> copy();
1884 $quo = $x -> bzero();
1885 } else {
1886 $rem = $self -> binf($y -> {sign});
1887 $quo = $x -> bone('-');
1888 }
1889 return ($quo, $rem);
1890 } else {
1891 if ($y -> is_inf()) {
1892 if ($x -> is_nan() || $x -> is_inf()) {
1893 return $x -> bnan();
1894 } else {
1895 return $x -> bzero();
1896 }
1897 }
1898 }
1899 }
1900
1901 # At this point, both the numerator and denominator are finite numbers, and
1902 # the denominator (divisor) is non-zero.
1903
1904 # x == 0?
394e6ffb 1905 return wantarray ? ($x,$self->bzero()) : $x if $x->is_zero();
0716bf9b 1906
9393ace2
JH
1907 # upgrade ?
1908 return $upgrade->bdiv($upgrade->new($x),$y,$a,$p,$r) if defined $upgrade;
13a12e00 1909
58cde26e 1910 # we need to limit the accuracy to protect against overflow
574bacfe 1911 my $fallback = 0;
990fb837
RGS
1912 my (@params,$scale);
1913 ($x,@params) = $x->_find_round_parameters($a,$p,$r,$y);
1914
1915 return $x if $x->is_nan(); # error in _find_round_parameters?
ee15d750
JH
1916
1917 # no rounding at all, so must use fallback
b4a10d33 1918 if (scalar @params == 0)
58cde26e 1919 {
0716bf9b 1920 # simulate old behaviour
990fb837
RGS
1921 $params[0] = $self->div_scale(); # and round to it as accuracy
1922 $scale = $params[0]+4; # at least four more for proper round
1923 $params[2] = $r; # round mode by caller or undef
ee15d750 1924 $fallback = 1; # to clear a/p afterwards
b4a10d33 1925 } else {
ee15d750
JH
1926 # the 4 below is empirical, and there might be cases where it is not
1927 # enough...
990fb837 1928 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
a0d0e21e 1929 }
03874afe 1930
b4a10d33
PJA
1931 my $rem;
1932 $rem = $self -> bzero() if wantarray;
03874afe
T
1933
1934 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1935
b4a10d33 1936 my $lx = $MBI -> _len($x->{_m}); my $ly = $MBI -> _len($y->{_m});
58cde26e 1937 $scale = $lx if $lx > $scale;
58cde26e 1938 $scale = $ly if $ly > $scale;
0716bf9b
JH
1939 my $diff = $ly - $lx;
1940 $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
a87115f0 1941
233f7bc0
T
1942 # check that $y is not 1 nor -1 and cache the result:
1943 my $y_not_one = !($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m}));
1944
7596a890
RGS
1945 # flipping the sign of $y will also flip the sign of $x for the special
1946 # case of $x->bsub($x); so we can catch it below:
1947 my $xsign = $x->{sign};
1948 $y->{sign} =~ tr/+-/-+/;
1949
a87115f0 1950 if ($xsign ne $x->{sign})
b3abae2a 1951 {
a87115f0 1952 # special case of $x /= $x results in 1
233f7bc0 1953 $x->bone(); # "fixes" also sign of $y, since $x is $y
b3abae2a 1954 }
03874afe 1955 else
58cde26e 1956 {
a87115f0
RGS
1957 # correct $y's sign again
1958 $y->{sign} =~ tr/+-/-+/;
1959 # continue with normal div code:
1960
f603091d 1961 # make copy of $x in case of list context for later remainder calculation
a87115f0 1962 if (wantarray && $y_not_one)
03874afe
T
1963 {
1964 $rem = $x->copy();
1965 }
394e6ffb 1966
03874afe 1967 $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
9b924220 1968
03874afe 1969 # check for / +-1 ( +/- 1E0)
a87115f0 1970 if ($y_not_one)
03874afe
T
1971 {
1972 # promote BigInts and it's subclasses (except when already a BigFloat)
1973 $y = $self->new($y) unless $y->isa('Math::BigFloat');
1974
1975 # calculate the result to $scale digits and then round it
1976 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
1977 $MBI->_lsft($x->{_m},$MBI->_new($scale),10);
1978 $MBI->_div ($x->{_m},$y->{_m}); # a/c
1979
1980 # correct exponent of $x
1981 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $y->{_e}, $x->{_es}, $y->{_es});
1982 # correct for 10**scale
1983 ($x->{_e},$x->{_es}) = _e_sub($x->{_e}, $MBI->_new($scale), $x->{_es}, '+');
1984 $x->bnorm(); # remove trailing 0's
1985 }
ada8209b 1986 } # end else $x != $y
a5f75d66 1987
091c87b1 1988 # shortcut to not run through _find_round_parameters again
990fb837 1989 if (defined $params[0])
ee15d750 1990 {
ef9466ea 1991 delete $x->{_a}; # clear before round
990fb837 1992 $x->bround($params[0],$params[2]); # then round accordingly
ee15d750
JH
1993 }
1994 else
1995 {
ef9466ea 1996 delete $x->{_p}; # clear before round
990fb837 1997 $x->bfround($params[1],$params[2]); # then round accordingly
ee15d750 1998 }
574bacfe
JH
1999 if ($fallback)
2000 {
2001 # clear a/p after round, since user did not request it
ef9466ea 2002 delete $x->{_a}; delete $x->{_p};
574bacfe 2003 }
03874afe 2004
58cde26e
JH
2005 if (wantarray)
2006 {
a87115f0 2007 if ($y_not_one)
394e6ffb 2008 {
b4a10d33 2009 $x -> bfloor();
990fb837 2010 $rem->bmod($y,@params); # copy already done
394e6ffb 2011 }
574bacfe
JH
2012 if ($fallback)
2013 {
2014 # clear a/p after round, since user did not request it
ef9466ea 2015 delete $rem->{_a}; delete $rem->{_p};
574bacfe 2016 }
0716bf9b 2017 return ($x,$rem);
58cde26e 2018 }
9393ace2 2019 $x;
58cde26e 2020 }
a0d0e21e 2021
58cde26e
JH
2022sub bmod
2023 {
f603091d 2024 # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return remainder
f9a08e12
JH
2025
2026 # set up parameters
2027 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2028 # objectify is costly, so avoid it
2029 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2030 {
2031 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2032 }
a0d0e21e 2033
50109ad0
RGS
2034 return $x if $x->modify('bmod');
2035
b4a10d33
PJA
2036 # At least one argument is NaN. This is handled the same way as in
2037 # Math::BigInt -> bmod().
2038
2039 if ($x -> is_nan() || $y -> is_nan()) {
2040 return $x -> bnan();
61f5c3f5 2041 }
b4a10d33
PJA
2042
2043 # Modulo zero. This is handled the same way as in Math::BigInt -> bmod().
2044
2045 if ($y -> is_zero()) {
9b924220
RGS
2046 return $x;
2047 }
233f7bc0 2048
b4a10d33
PJA
2049 # Numerator (dividend) is +/-inf. This is handled the same way as in
2050 # Math::BigInt -> bmod().
2051
2052 if ($x -> is_inf()) {
2053 return $x -> bnan();
2054 }
2055
2056 # Denominator (divisor) is +/-inf. This is handled the same way as in
2057 # Math::BigInt -> bmod().
2058
2059 if ($y -> is_inf()) {
2060 if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) {
2061 return $x;
2062 } else {
2063 return $x -> binf($y -> sign());
2064 }
2065 }
2066
7596a890
RGS
2067 return $x->bzero() if $x->is_zero()
2068 || ($x->is_int() &&
2069 # check that $y == +1 or $y == -1:
2070 ($MBI->_is_zero($y->{_e}) && $MBI->_is_one($y->{_m})));
58cde26e 2071
61f5c3f5 2072 my $cmp = $x->bacmp($y); # equal or $x < $y?
b4a10d33
PJA
2073 if ($cmp == 0) { # $x == $y => result 0
2074 return $x -> bzero($a, $p);
2075 }
61f5c3f5
T
2076
2077 # only $y of the operands negative?
b4a10d33 2078 my $neg = $x->{sign} ne $y->{sign} ? 1 : 0;
61f5c3f5
T
2079
2080 $x->{sign} = $y->{sign}; # calc sign first
b4a10d33
PJA
2081 if ($cmp < 0 && $neg == 0) { # $x < $y => result $x
2082 return $x -> round($a, $p, $r);
2083 }
61f5c3f5 2084
9b924220 2085 my $ym = $MBI->_copy($y->{_m});
61f5c3f5
T
2086
2087 # 2e1 => 20
9b924220
RGS
2088 $MBI->_lsft( $ym, $y->{_e}, 10)
2089 if $y->{_es} eq '+' && !$MBI->_is_zero($y->{_e});
61f5c3f5
T
2090
2091 # if $y has digits after dot
2092 my $shifty = 0; # correct _e of $x by this
9b924220 2093 if ($y->{_es} eq '-') # has digits after dot
61f5c3f5
T
2094 {
2095 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
9b924220
RGS
2096 $shifty = $MBI->_num($y->{_e}); # no more digits after dot
2097 $MBI->_lsft($x->{_m}, $y->{_e}, 10);# 123 => 1230, $y->{_m} is already 25
61f5c3f5
T
2098 }
2099 # $ym is now mantissa of $y based on exponent 0
b3abae2a 2100
61f5c3f5 2101 my $shiftx = 0; # correct _e of $x by this
9b924220 2102 if ($x->{_es} eq '-') # has digits after dot
61f5c3f5
T
2103 {
2104 # 123.4 % 20 => 1234 % 200
9b924220
RGS
2105 $shiftx = $MBI->_num($x->{_e}); # no more digits after dot
2106 $MBI->_lsft($ym, $x->{_e}, 10); # 123 => 1230
61f5c3f5
T
2107 }
2108 # 123e1 % 20 => 1230 % 20
9b924220 2109 if ($x->{_es} eq '+' && !$MBI->_is_zero($x->{_e}))
61f5c3f5 2110 {
9b924220 2111 $MBI->_lsft( $x->{_m}, $x->{_e},10); # es => '+' here
61f5c3f5 2112 }
9b924220
RGS
2113
2114 $x->{_e} = $MBI->_new($shiftx);
2115 $x->{_es} = '+';
2116 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
2117 $MBI->_add( $x->{_e}, $MBI->_new($shifty)) if $shifty != 0;
61f5c3f5
T
2118
2119 # now mantissas are equalized, exponent of $x is adjusted, so calc result
b3abae2a 2120
9b924220 2121 $x->{_m} = $MBI->_mod( $x->{_m}, $ym);
61f5c3f5 2122
9b924220 2123 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
61f5c3f5
T
2124 $x->bnorm();
2125
b4a10d33 2126 if ($neg != 0 && ! $x -> is_zero()) # one of them negative => correct in place
61f5c3f5
T
2127 {
2128 my $r = $y - $x;
2129 $x->{_m} = $r->{_m};
2130 $x->{_e} = $r->{_e};
9b924220
RGS
2131 $x->{_es} = $r->{_es};
2132 $x->{sign} = '+' if $MBI->_is_zero($x->{_m}); # fix sign for -0
61f5c3f5
T
2133 $x->bnorm();
2134 }
2135
2136 $x->round($a,$p,$r,$y); # round and return
58cde26e
JH
2137 }
2138
990fb837
RGS
2139sub broot
2140 {
2141 # calculate $y'th root of $x
3a427a11
RGS
2142
2143 # set up parameters
2144 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2145 # objectify is costly, so avoid it
2146 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2147 {
2148 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2149 }
990fb837 2150
50109ad0
RGS
2151 return $x if $x->modify('broot');
2152
990fb837
RGS
2153 # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
2154 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
2155 $y->{sign} !~ /^\+$/;
2156
2157 return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
2158
2159 # we need to limit the accuracy to protect against overflow
2160 my $fallback = 0;
2161 my (@params,$scale);
2162 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2163
2164 return $x if $x->is_nan(); # error in _find_round_parameters?
2165
2166 # no rounding at all, so must use fallback
2167 if (scalar @params == 0)
2168 {
2169 # simulate old behaviour
2170 $params[0] = $self->div_scale(); # and round to it as accuracy
2171 $scale = $params[0]+4; # at least four more for proper round
ada8209b 2172 $params[2] = $r; # round mode by caller or undef
990fb837
RGS
2173 $fallback = 1; # to clear a/p afterwards
2174 }
2175 else
2176 {
2177 # the 4 below is empirical, and there might be cases where it is not
2178 # enough...
2179 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2180 }
2181
2182 # when user set globals, they would interfere with our calculation, so
2183 # disable them and later re-enable them
2184 no strict 'refs';
2185 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2186 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2187 # we also need to disable any set A or P on $x (_find_round_parameters took
2188 # them already into account), since these would interfere, too
2189 delete $x->{_a}; delete $x->{_p};
2190 # need to disable $upgrade in BigInt, to avoid deep recursion
2191 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
2192
2193 # remember sign and make $x positive, since -4 ** (1/2) => -2
27e7b8bb
T
2194 my $sign = 0; $sign = 1 if $x->{sign} eq '-'; $x->{sign} = '+';
2195
2196 my $is_two = 0;
2197 if ($y->isa('Math::BigFloat'))
2198 {
2199 $is_two = ($y->{sign} eq '+' && $MBI->_is_two($y->{_m}) && $MBI->_is_zero($y->{_e}));
2200 }
2201 else
2202 {
2203 $is_two = ($y == 2);
2204 }
990fb837 2205
27e7b8bb
T
2206 # normal square root if $y == 2:
2207 if ($is_two)
990fb837
RGS
2208 {
2209 $x->bsqrt($scale+4);
2210 }
2211 elsif ($y->is_one('-'))
2212 {
2213 # $x ** -1 => 1/$x
2214 my $u = $self->bone()->bdiv($x,$scale);
2215 # copy private parts over
2216 $x->{_m} = $u->{_m};
2217 $x->{_e} = $u->{_e};
9b924220 2218 $x->{_es} = $u->{_es};
990fb837
RGS
2219 }
2220 else
2221 {
3a427a11
RGS
2222 # calculate the broot() as integer result first, and if it fits, return
2223 # it rightaway (but only if $x and $y are integer):
2224
2225 my $done = 0; # not yet
2226 if ($y->is_int() && $x->is_int())
2227 {
9b924220
RGS
2228 my $i = $MBI->_copy( $x->{_m} );
2229 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2230 my $int = Math::BigInt->bzero();
2231 $int->{value} = $i;
3a427a11
RGS
2232 $int->broot($y->as_number());
2233 # if ($exact)
2234 if ($int->copy()->bpow($y) == $x)
2235 {
2236 # found result, return it
9b924220
RGS
2237 $x->{_m} = $int->{value};
2238 $x->{_e} = $MBI->_zero();
2239 $x->{_es} = '+';
3a427a11
RGS
2240 $x->bnorm();
2241 $done = 1;
2242 }
2243 }
2244 if ($done == 0)
2245 {
2246 my $u = $self->bone()->bdiv($y,$scale+4);
2247 delete $u->{_a}; delete $u->{_p}; # otherwise it conflicts
2248 $x->bpow($u,$scale+4); # el cheapo
2249 }
990fb837
RGS
2250 }
2251 $x->bneg() if $sign == 1;
2252
091c87b1 2253 # shortcut to not run through _find_round_parameters again
990fb837
RGS
2254 if (defined $params[0])
2255 {
2256 $x->bround($params[0],$params[2]); # then round accordingly
2257 }
2258 else
2259 {
2260 $x->bfround($params[1],$params[2]); # then round accordingly
2261 }
2262 if ($fallback)
2263 {
2264 # clear a/p after round, since user did not request it
ef9466ea 2265 delete $x->{_a}; delete $x->{_p};
990fb837
RGS
2266 }
2267 # restore globals
2268 $$abr = $ab; $$pbr = $pb;
2269 $x;
2270 }
2271
58cde26e
JH
2272sub bsqrt
2273 {
990fb837 2274 # calculate square root
ee15d750 2275 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e 2276
50109ad0
RGS
2277 return $x if $x->modify('bsqrt');
2278
990fb837
RGS
2279 return $x->bnan() if $x->{sign} !~ /^[+]/; # NaN, -inf or < 0
2280 return $x if $x->{sign} eq '+inf'; # sqrt(inf) == inf
2281 return $x->round($a,$p,$r) if $x->is_zero() || $x->is_one();
58cde26e 2282
61f5c3f5 2283 # we need to limit the accuracy to protect against overflow
574bacfe 2284 my $fallback = 0;
990fb837
RGS
2285 my (@params,$scale);
2286 ($x,@params) = $x->_find_round_parameters($a,$p,$r);
2287
2288 return $x if $x->is_nan(); # error in _find_round_parameters?
61f5c3f5
T
2289
2290 # no rounding at all, so must use fallback
990fb837 2291 if (scalar @params == 0)
0716bf9b
JH
2292 {
2293 # simulate old behaviour
990fb837
RGS
2294 $params[0] = $self->div_scale(); # and round to it as accuracy
2295 $scale = $params[0]+4; # at least four more for proper round
2296 $params[2] = $r; # round mode by caller or undef
ee15d750 2297 $fallback = 1; # to clear a/p afterwards
0716bf9b 2298 }
61f5c3f5
T
2299 else
2300 {
2301 # the 4 below is empirical, and there might be cases where it is not
2302 # enough...
990fb837 2303 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
61f5c3f5
T
2304 }
2305
2306 # when user set globals, they would interfere with our calculation, so
9393ace2 2307 # disable them and later re-enable them
61f5c3f5
T
2308 no strict 'refs';
2309 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
b3abae2a 2310 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
61f5c3f5
T
2311 # we also need to disable any set A or P on $x (_find_round_parameters took
2312 # them already into account), since these would interfere, too
2313 delete $x->{_a}; delete $x->{_p};
9393ace2
JH
2314 # need to disable $upgrade in BigInt, to avoid deep recursion
2315 local $Math::BigInt::upgrade = undef; # should be really parent class vs MBI
61f5c3f5 2316
9b924220
RGS
2317 my $i = $MBI->_copy( $x->{_m} );
2318 $MBI->_lsft( $i, $x->{_e}, 10 ) unless $MBI->_is_zero($x->{_e});
2319 my $xas = Math::BigInt->bzero();
2320 $xas->{value} = $i;
2321
394e6ffb 2322 my $gs = $xas->copy()->bsqrt(); # some guess
b3abae2a 2323
9b924220 2324 if (($x->{_es} ne '-') # guess can't be accurate if there are
394e6ffb 2325 # digits after the dot
b3abae2a 2326 && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
394e6ffb 2327 {
9b924220
RGS
2328 # exact result, copy result over to keep $x
2329 $x->{_m} = $gs->{value}; $x->{_e} = $MBI->_zero(); $x->{_es} = '+';
2330 $x->bnorm();
091c87b1 2331 # shortcut to not run through _find_round_parameters again
990fb837 2332 if (defined $params[0])
61f5c3f5 2333 {
990fb837 2334 $x->bround($params[0],$params[2]); # then round accordingly
61f5c3f5
T
2335 }
2336 else
2337 {
990fb837 2338 $x->bfround($params[1],$params[2]); # then round accordingly
61f5c3f5
T
2339 }
2340 if ($fallback)
2341 {
2342 # clear a/p after round, since user did not request it
ef9466ea 2343 delete $x->{_a}; delete $x->{_p};
61f5c3f5 2344 }
9393ace2 2345 # re-enable A and P, upgrade is taken care of by "local"
b3abae2a 2346 ${"$self\::accuracy"} = $ab; ${"$self\::precision"} = $pb;
61f5c3f5 2347 return $x;
394e6ffb 2348 }
2ab5f49d
T
2349
2350 # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy
9681bfa6 2351 # of the result by multiplying the input by 100 and then divide the integer
2ab5f49d 2352 # result of sqrt(input) by 10. Rounding afterwards returns the real result.
9b924220
RGS
2353
2354 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
2355 my $y1 = $MBI->_copy($x->{_m});
2356
2357 my $length = $MBI->_len($y1);
2358
2359 # Now calculate how many digits the result of sqrt(y1) would have
2360 my $digits = int($length / 2);
2361
2362 # But we need at least $scale digits, so calculate how many are missing
2363 my $shift = $scale - $digits;
2364
0dceeee6
RGS
2365 # This happens if the input had enough digits
2366 # (we take care of integer guesses above)
2367 $shift = 0 if $shift < 0;
9b924220
RGS
2368
2369 # Multiply in steps of 100, by shifting left two times the "missing" digits
2370 my $s2 = $shift * 2;
2371
2ab5f49d
T
2372 # We now make sure that $y1 has the same odd or even number of digits than
2373 # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
2374 # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
2375 # steps of 10. The length of $x does not count, since an even or odd number
2376 # of digits before the dot is not changed by adding an even number of digits
2377 # after the dot (the result is still odd or even digits long).
9b924220
RGS
2378 $s2++ if $MBI->_is_odd($x->{_e});
2379
2380 $MBI->_lsft( $y1, $MBI->_new($s2), 10);
2381
2ab5f49d 2382 # now take the square root and truncate to integer
9b924220
RGS
2383 $y1 = $MBI->_sqrt($y1);
2384
2ab5f49d
T
2385 # By "shifting" $y1 right (by creating a negative _e) we calculate the final
2386 # result, which is than later rounded to the desired scale.
990fb837
RGS
2387
2388 # calculate how many zeros $x had after the '.' (or before it, depending
9b924220
RGS
2389 # on sign of $dat, the result should have half as many:
2390 my $dat = $MBI->_num($x->{_e});
2391 $dat = -$dat if $x->{_es} eq '-';
2392 $dat += $length;
990fb837
RGS
2393
2394 if ($dat > 0)
2395 {
2396 # no zeros after the dot (e.g. 1.23, 0.49 etc)
2397 # preserve half as many digits before the dot than the input had
2398 # (but round this "up")
2399 $dat = int(($dat+1)/2);
2400 }
2401 else
2402 {
2403 $dat = int(($dat)/2);
2404 }
9b924220
RGS
2405 $dat -= $MBI->_len($y1);
2406 if ($dat < 0)
2407 {
2408 $dat = abs($dat);
2409 $x->{_e} = $MBI->_new( $dat );
2410 $x->{_es} = '-';
2411 }
2412 else
2413 {
2414 $x->{_e} = $MBI->_new( $dat );
2415 $x->{_es} = '+';
2416 }
2ab5f49d 2417 $x->{_m} = $y1;
9b924220 2418 $x->bnorm();
61f5c3f5 2419
091c87b1 2420 # shortcut to not run through _find_round_parameters again
990fb837 2421 if (defined $params[0])
61f5c3f5 2422 {
990fb837 2423 $x->bround($params[0],$params[2]); # then round accordingly
61f5c3f5
T
2424 }
2425 else
2426 {
990fb837 2427 $x->bfround($params[1],$params[2]); # then round accordingly
61f5c3f5 2428 }
574bacfe
JH
2429 if ($fallback)
2430 {
2431 # clear a/p after round, since user did not request it
ef9466ea 2432 delete $x->{_a}; delete $x->{_p};
574bacfe 2433 }
61f5c3f5 2434 # restore globals
b3abae2a 2435 $$abr = $ab; $$pbr = $pb;
574bacfe 2436 $x;
58cde26e
JH
2437 }
2438
b3abae2a
JH
2439sub bfac
2440 {
28df3e88 2441 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
091c87b1 2442 # compute factorial number, modifies first argument
b3abae2a 2443
b282a552
T
2444 # set up parameters
2445 my ($self,$x,@r) = (ref($_[0]),@_);
2446 # objectify is costly, so avoid it
2447 ($self,$x,@r) = objectify(1,@_) if !ref($x);
2448
50109ad0
RGS
2449 # inf => inf
2450 return $x if $x->modify('bfac') || $x->{sign} eq '+inf';
2451
28df3e88
JH
2452 return $x->bnan()
2453 if (($x->{sign} ne '+') || # inf, NaN, <0 etc => NaN
9b924220 2454 ($x->{_es} ne '+')); # digits after dot?
b3abae2a 2455
b3abae2a 2456 # use BigInt's bfac() for faster calc
9b924220 2457 if (! $MBI->_is_zero($x->{_e}))
091c87b1 2458 {
9b924220
RGS
2459 $MBI->_lsft($x->{_m}, $x->{_e},10); # change 12e1 to 120e0
2460 $x->{_e} = $MBI->_zero(); # normalize
2461 $x->{_es} = '+';
091c87b1 2462 }
9b924220 2463 $MBI->_fac($x->{_m}); # calculate factorial
091c87b1 2464 $x->bnorm()->round(@r); # norm again and round result
b3abae2a
JH
2465 }
2466
9393ace2
JH
2467sub _pow
2468 {
60a1aa19
T
2469 # Calculate a power where $y is a non-integer, like 2 ** 0.3
2470 my ($x,$y,@r) = @_;
9393ace2
JH
2471 my $self = ref($x);
2472
2473 # if $y == 0.5, it is sqrt($x)
27e7b8bb 2474 $HALF = $self->new($HALF) unless ref($HALF);
60a1aa19 2475 return $x->bsqrt(@r,$y) if $y->bcmp($HALF) == 0;
9393ace2 2476
990fb837
RGS
2477 # Using:
2478 # a ** x == e ** (x * ln a)
2479
9393ace2 2480 # u = y * ln x
990fb837
RGS
2481 # _ _
2482 # Taylor: | u u^2 u^3 |
2483 # x ** y = 1 + | --- + --- + ----- + ... |
2484 # |_ 1 1*2 1*2*3 _|
9393ace2
JH
2485
2486 # we need to limit the accuracy to protect against overflow
2487 my $fallback = 0;
990fb837 2488 my ($scale,@params);
60a1aa19 2489 ($x,@params) = $x->_find_round_parameters(@r);
990fb837
RGS
2490
2491 return $x if $x->is_nan(); # error in _find_round_parameters?
9393ace2
JH
2492
2493 # no rounding at all, so must use fallback
990fb837 2494 if (scalar @params == 0)
9393ace2
JH
2495 {
2496 # simulate old behaviour
990fb837
RGS
2497 $params[0] = $self->div_scale(); # and round to it as accuracy
2498 $params[1] = undef; # disable P
2499 $scale = $params[0]+4; # at least four more for proper round
60a1aa19 2500 $params[2] = $r[2]; # round mode by caller or undef
9393ace2
JH
2501 $fallback = 1; # to clear a/p afterwards
2502 }
2503 else
2504 {
2505 # the 4 below is empirical, and there might be cases where it is not
2506 # enough...
990fb837 2507 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
9393ace2
JH
2508 }
2509
2510 # when user set globals, they would interfere with our calculation, so
56d9de68 2511 # disable them and later re-enable them
9393ace2
JH
2512 no strict 'refs';
2513 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2514 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2515 # we also need to disable any set A or P on $x (_find_round_parameters took
2516 # them already into account), since these would interfere, too
2517 delete $x->{_a}; delete $x->{_p};
2518 # need to disable $upgrade in BigInt, to avoid deep recursion
2519 local $Math::BigInt::upgrade = undef;
2520
2521 my ($limit,$v,$u,$below,$factor,$next,$over);
2522
990fb837 2523 $u = $x->copy()->blog(undef,$scale)->bmul($y);
9393ace2
JH
2524 $v = $self->bone(); # 1
2525 $factor = $self->new(2); # 2
2526 $x->bone(); # first term: 1
2527
2528 $below = $v->copy();
2529 $over = $u->copy();
ae161977 2530
9393ace2
JH
2531 $limit = $self->new("1E-". ($scale-1));
2532 #my $steps = 0;
2533 while (3 < 5)
2534 {
2535 # we calculate the next term, and add it to the last
2536 # when the next term is below our limit, it won't affect the outcome
7d193e39 2537 # anymore, so we stop:
9393ace2 2538 $next = $over->copy()->bdiv($below,$scale);
990fb837 2539 last if $next->bacmp($limit) <= 0;
9393ace2 2540 $x->badd($next);
9393ace2
JH
2541 # calculate things for the next term
2542 $over *= $u; $below *= $factor; $factor->binc();
9b924220
RGS
2543
2544 last if $x->{sign} !~ /^[-+]$/;
2545
9393ace2
JH
2546 #$steps++;
2547 }
2548
091c87b1 2549 # shortcut to not run through _find_round_parameters again
990fb837 2550 if (defined $params[0])
9393ace2 2551 {
990fb837 2552 $x->bround($params[0],$params[2]); # then round accordingly
9393ace2
JH
2553 }
2554 else
2555 {
990fb837 2556 $x->bfround($params[1],$params[2]); # then round accordingly
9393ace2
JH
2557 }
2558 if ($fallback)
2559 {
2560 # clear a/p after round, since user did not request it
ef9466ea 2561 delete $x->{_a}; delete $x->{_p};
9393ace2
JH
2562 }
2563 # restore globals
2564 $$abr = $ab; $$pbr = $pb;
2565 $x;
2566 }
2567
58cde26e
JH
2568sub bpow
2569 {
2570 # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2571 # compute power of two numbers, second arg is used as integer
2572 # modifies first argument
2573
f9a08e12
JH
2574 # set up parameters
2575 my ($self,$x,$y,$a,$p,$r) = (ref($_[0]),@_);
2576 # objectify is costly, so avoid it
2577 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
2578 {
2579 ($self,$x,$y,$a,$p,$r) = objectify(2,@_);
2580 }
58cde26e 2581
50109ad0
RGS
2582 return $x if $x->modify('bpow');
2583
58cde26e 2584 return $x->bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
2d2b2744
T
2585 return $x if $x->{sign} =~ /^[+-]inf$/;
2586
ae161977
RGS
2587 # cache the result of is_zero
2588 my $y_is_zero = $y->is_zero();
2589 return $x->bone() if $y_is_zero;
58cde26e 2590 return $x if $x->is_one() || $y->is_one();
9393ace2 2591
ae161977 2592 my $x_is_zero = $x->is_zero();
2d2b2744 2593 return $x->_pow($y,$a,$p,$r) if !$x_is_zero && !$y->is_int(); # non-integer power
9393ace2 2594
ae161977 2595 my $y1 = $y->as_number()->{value}; # make MBI part
9b924220 2596
394e6ffb 2597 # if ($x == -1)
9b924220 2598 if ($x->{sign} eq '-' && $MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
58cde26e
JH
2599 {
2600 # if $x == -1 and odd/even y => +1/-1 because +-1 ^ (+-1) => +-1
9b924220 2601 return $MBI->_is_odd($y1) ? $x : $x->babs(1);
288d023a 2602 }
ae161977 2603 if ($x_is_zero)
28df3e88
JH
2604 {
2605 return $x if $y->{sign} eq '+'; # 0**y => 0 (if not y <= 0)
9b924220
RGS
2606 # 0 ** -y => 1 / (0 ** y) => 1 / 0! (1 / 0 => +inf)
2607 return $x->binf();
28df3e88 2608 }
58cde26e 2609
9b924220 2610 my $new_sign = '+';
ae161977 2611 $new_sign = $MBI->_is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
9b924220 2612
58cde26e 2613 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
9b924220 2614 $x->{_m} = $MBI->_pow( $x->{_m}, $y1);
ae161977 2615 $x->{_e} = $MBI->_mul ($x->{_e}, $y1);
9b924220
RGS
2616
2617 $x->{sign} = $new_sign;
58cde26e
JH
2618 $x->bnorm();
2619 if ($y->{sign} eq '-')
2620 {
2621 # modify $x in place!
ae161977 2622 my $z = $x->copy(); $x->bone();
7b29e1e6 2623 return scalar $x->bdiv($z,$a,$p,$r); # round in one go (might ignore y's A!)
a0d0e21e 2624 }
28df3e88 2625 $x->round($a,$p,$r,$y);
58cde26e
JH
2626 }
2627
80365507
T
2628sub bmodpow
2629 {
2630 # takes a very large number to a very large exponent in a given very
c4a6f826 2631 # large modulus, quickly, thanks to binary exponentiation. Supports
80365507
T
2632 # negative exponents.
2633 my ($self,$num,$exp,$mod,@r) = objectify(3,@_);
2634
2635 return $num if $num->modify('bmodpow');
2636
2637 # check modulus for valid values
2638 return $num->bnan() if ($mod->{sign} ne '+' # NaN, - , -inf, +inf
2639 || $mod->is_zero());
2640
2641 # check exponent for valid values
2642 if ($exp->{sign} =~ /\w/)
2643 {
2644 # i.e., if it's NaN, +inf, or -inf...
2645 return $num->bnan();
2646 }
2647
2648 $num->bmodinv ($mod) if ($exp->{sign} eq '-');
2649
2650 # check num for valid values (also NaN if there was no inverse but $exp < 0)
2651 return $num->bnan() if $num->{sign} !~ /^[+-]$/;
2652
2653 # $mod is positive, sign on $exp is ignored, result also positive
2654
2655 # XXX TODO: speed it up when all three numbers are integers
2656 $num->bpow($exp)->bmod($mod);
2657 }
2658
58cde26e 2659###############################################################################
fdb4b05f
T
2660# trigonometric functions
2661
20e2035c 2662# helper function for bpi() and batan2(), calculates arcus tanges (1/x)
fdb4b05f 2663
20e2035c 2664sub _atan_inv
fdb4b05f 2665 {
20e2035c
T
2666 # return a/b so that a/b approximates atan(1/x) to at least limit digits
2667 my ($self, $x, $limit) = @_;
fdb4b05f 2668
20e2035c
T
2669 # Taylor: x^3 x^5 x^7 x^9
2670 # atan = x - --- + --- - --- + --- - ...
2671 # 3 5 7 9
fdb4b05f 2672
20e2035c
T
2673 # 1 1 1 1
2674 # atan 1/x = - - ------- + ------- - ------- + ...
2675 # x x^3 * 3 x^5 * 5 x^7 * 7
2676
2677 # 1 1 1 1
2678 # atan 1/x = - - --------- + ---------- - ----------- + ...
2679 # 5 3 * 125 5 * 3125 7 * 78125
2680
2681 # Subtraction/addition of a rational:
2682
2683 # 5 7 5*3 +- 7*4
2684 # - +- - = ----------
2685 # 4 3 4*3
2686
2687 # Term: N N+1
2688 #
2689 # a 1 a * d * c +- b
2690 # ----- +- ------------------ = ----------------
2691 # b d * c b * d * c
2692
2693 # since b1 = b0 * (d-2) * c
2694
2695 # a 1 a * d +- b / c
2696 # ----- +- ------------------ = ----------------
2697 # b d * c b * d
2698
2699 # and d = d + 2
2700 # and c = c * x * x
2701
2702 # u = d * c
2703 # stop if length($u) > limit
2704 # a = a * u +- b
2705 # b = b * u
2706 # d = d + 2
2707 # c = c * x * x
2708 # sign = 1 - sign
2709
2710 my $a = $MBI->_one();
30afc38d 2711 my $b = $MBI->_copy($x);
fdb4b05f 2712
30afc38d 2713 my $x2 = $MBI->_mul( $MBI->_copy($x), $b); # x2 = x * x
20e2035c 2714 my $d = $MBI->_new( 3 ); # d = 3
30afc38d 2715 my $c = $MBI->_mul( $MBI->_copy($x), $x2); # c = x ^ 3
20e2035c
T
2716 my $two = $MBI->_new( 2 );
2717
2718 # run the first step unconditionally
2719 my $u = $MBI->_mul( $MBI->_copy($d), $c);
2720 $a = $MBI->_mul($a, $u);
2721 $a = $MBI->_sub($a, $b);
2722 $b = $MBI->_mul($b, $u);
2723 $d = $MBI->_add($d, $two);
2724 $c = $MBI->_mul($c, $x2);
2725
2726 # a is now a * (d-3) * c
2727 # b is now b * (d-2) * c
2728
2729 # run the second step unconditionally
2730 $u = $MBI->_mul( $MBI->_copy($d), $c);
2731 $a = $MBI->_mul($a, $u);
2732 $a = $MBI->_add($a, $b);
2733 $b = $MBI->_mul($b, $u);
2734 $d = $MBI->_add($d, $two);
2735 $c = $MBI->_mul($c, $x2);
2736
2737 # a is now a * (d-3) * (d-5) * c * c
2738 # b is now b * (d-2) * (d-4) * c * c
2739
2740 # so we can remove c * c from both a and b to shorten the numbers involved:
2741 $a = $MBI->_div($a, $x2);
2742 $b = $MBI->_div($b, $x2);
2743 $a = $MBI->_div($a, $x2);
2744 $b = $MBI->_div($b, $x2);
2745
2746# my $step = 0;
2747 my $sign = 0; # 0 => -, 1 => +
2748 while (3 < 5)
fdb4b05f 2749 {
20e2035c
T
2750# $step++;
2751# if (($i++ % 100) == 0)
2752# {
2753# print "a=",$MBI->_str($a),"\n";
2754# print "b=",$MBI->_str($b),"\n";
2755# }
2756# print "d=",$MBI->_str($d),"\n";
2757# print "x2=",$MBI->_str($x2),"\n";
2758# print "c=",$MBI->_str($c),"\n";
2759
2760 my $u = $MBI->_mul( $MBI->_copy($d), $c);
2761 # use _alen() for libs like GMP where _len() would be O(N^2)
2762 last if $MBI->_alen($u) > $limit;
2763 my ($bc,$r) = $MBI->_div( $MBI->_copy($b), $c);
2764 if ($MBI->_is_zero($r))
fdb4b05f 2765 {
20e2035c
T
2766 # b / c is an integer, so we can remove c from all terms
2767 # this happens almost every time:
2768 $a = $MBI->_mul($a, $d);
2769 $a = $MBI->_sub($a, $bc) if $sign == 0;
2770 $a = $MBI->_add($a, $bc) if $sign == 1;
2771 $b = $MBI->_mul($b, $d);
fdb4b05f
T
2772 }
2773 else
2774 {
20e2035c
T
2775 # b / c is not an integer, so we keep c in the terms
2776 # this happens very rarely, for instance for x = 5, this happens only
2777 # at the following steps:
2778 # 1, 5, 14, 32, 72, 157, 340, ...
2779 $a = $MBI->_mul($a, $u);
2780 $a = $MBI->_sub($a, $b) if $sign == 0;
2781 $a = $MBI->_add($a, $b) if $sign == 1;
2782 $b = $MBI->_mul($b, $u);
fdb4b05f 2783 }
20e2035c
T
2784 $d = $MBI->_add($d, $two);
2785 $c = $MBI->_mul($c, $x2);
2786 $sign = 1 - $sign;
2787
fdb4b05f
T
2788 }
2789
30afc38d 2790# print "Took $step steps for ", $MBI->_str($x),"\n";
20e2035c
T
2791# print "a=",$MBI->_str($a),"\n"; print "b=",$MBI->_str($b),"\n";
2792 # return a/b so that a/b approximates atan(1/x)
2793 ($a,$b);
fdb4b05f
T
2794 }
2795
2796sub bpi
2797 {
fdb4b05f 2798 my ($self,$n) = @_;
80365507
T
2799 if (@_ == 0)
2800 {
2801 $self = $class;
2802 }
fdb4b05f
T
2803 if (@_ == 1)
2804 {
2805 # called like Math::BigFloat::bpi(10);
2806 $n = $self; $self = $class;
60a1aa19
T
2807 # called like Math::BigFloat->bpi();
2808 $n = undef if $n eq 'Math::BigFloat';
fdb4b05f
T
2809 }
2810 $self = ref($self) if ref($self);
36ec1dfe 2811 my $fallback = defined $n ? 0 : 1;
fdb4b05f
T
2812 $n = 40 if !defined $n || $n < 1;
2813
20e2035c
T
2814 # after 黃見利 (Hwang Chien-Lih) (1997)
2815 # pi/4 = 183 * atan(1/239) + 32 * atan(1/1023) – 68 * atan(1/5832)
2816 # + 12 * atan(1/110443) - 12 * atan(1/4841182) - 100 * atan(1/6826318)
2817
2818 # a few more to prevent rounding errors
2819 $n += 4;
2820
30afc38d
T
2821 my ($a,$b) = $self->_atan_inv( $MBI->_new(239),$n);
2822 my ($c,$d) = $self->_atan_inv( $MBI->_new(1023),$n);
2823 my ($e,$f) = $self->_atan_inv( $MBI->_new(5832),$n);
2824 my ($g,$h) = $self->_atan_inv( $MBI->_new(110443),$n);
2825 my ($i,$j) = $self->_atan_inv( $MBI->_new(4841182),$n);
2826 my ($k,$l) = $self->_atan_inv( $MBI->_new(6826318),$n);
20e2035c
T
2827
2828 $MBI->_mul($a, $MBI->_new(732));
2829 $MBI->_mul($c, $MBI->_new(128));
2830 $MBI->_mul($e, $MBI->_new(272));
2831 $MBI->_mul($g, $MBI->_new(48));
2832 $MBI->_mul($i, $MBI->_new(48));
2833 $MBI->_mul($k, $MBI->_new(400));
2834
2835 my $x = $self->bone(); $x->{_m} = $a; my $x_d = $self->bone(); $x_d->{_m} = $b;
2836 my $y = $self->bone(); $y->{_m} = $c; my $y_d = $self->bone(); $y_d->{_m} = $d;
2837 my $z = $self->bone(); $z->{_m} = $e; my $z_d = $self->bone(); $z_d->{_m} = $f;
2838 my $u = $self->bone(); $u->{_m} = $g; my $u_d = $self->bone(); $u_d->{_m} = $h;
2839 my $v = $self->bone(); $v->{_m} = $i; my $v_d = $self->bone(); $v_d->{_m} = $j;
2840 my $w = $self->bone(); $w->{_m} = $k; my $w_d = $self->bone(); $w_d->{_m} = $l;
2841 $x->bdiv($x_d, $n);
2842 $y->bdiv($y_d, $n);
2843 $z->bdiv($z_d, $n);
2844 $u->bdiv($u_d, $n);
2845 $v->bdiv($v_d, $n);
2846 $w->bdiv($w_d, $n);
2847
36ec1dfe
T
2848 delete $x->{_a}; delete $y->{_a}; delete $z->{_a};
2849 delete $u->{_a}; delete $v->{_a}; delete $w->{_a};
20e2035c
T
2850 $x->badd($y)->bsub($z)->badd($u)->bsub($v)->bsub($w);
2851
36ec1dfe
T
2852 $x->bround($n-4);
2853 delete $x->{_a} if $fallback == 1;
2854 $x;
fdb4b05f
T
2855 }
2856
60a1aa19
T
2857sub bcos
2858 {
2859 # Calculate a cosinus of x.
2860 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2861
2862 # Taylor: x^2 x^4 x^6 x^8
2863 # cos = 1 - --- + --- - --- + --- ...
2864 # 2! 4! 6! 8!
2865
2866 # we need to limit the accuracy to protect against overflow
2867 my $fallback = 0;
2868 my ($scale,@params);
2869 ($x,@params) = $x->_find_round_parameters(@r);
2870
2871 # constant object or error in _find_round_parameters?
2872 return $x if $x->modify('bcos') || $x->is_nan();
2873
2874 return $x->bone(@r) if $x->is_zero();
2875
2876 # no rounding at all, so must use fallback
2877 if (scalar @params == 0)
2878 {
2879 # simulate old behaviour
2880 $params[0] = $self->div_scale(); # and round to it as accuracy
2881 $params[1] = undef; # disable P
2882 $scale = $params[0]+4; # at least four more for proper round
2883 $params[2] = $r[2]; # round mode by caller or undef
2884 $fallback = 1; # to clear a/p afterwards
2885 }
2886 else
2887 {
2888 # the 4 below is empirical, and there might be cases where it is not
2889 # enough...
2890 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2891 }
2892
2893 # when user set globals, they would interfere with our calculation, so
2894 # disable them and later re-enable them
2895 no strict 'refs';
2896 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2897 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2898 # we also need to disable any set A or P on $x (_find_round_parameters took
2899 # them already into account), since these would interfere, too
2900 delete $x->{_a}; delete $x->{_p};
2901 # need to disable $upgrade in BigInt, to avoid deep recursion
2902 local $Math::BigInt::upgrade = undef;
2903
2904 my $last = 0;
2905 my $over = $x * $x; # X ^ 2
2906 my $x2 = $over->copy(); # X ^ 2; difference between terms
2907 my $sign = 1; # start with -=
2908 my $below = $self->new(2); my $factorial = $self->new(3);
36ec1dfe 2909 $x->bone(); delete $x->{_a}; delete $x->{_p};
60a1aa19
T
2910
2911 my $limit = $self->new("1E-". ($scale-1));
2912 #my $steps = 0;
2913 while (3 < 5)
2914 {
2915 # we calculate the next term, and add it to the last
2916 # when the next term is below our limit, it won't affect the outcome
2917 # anymore, so we stop:
2918 my $next = $over->copy()->bdiv($below,$scale);
2919 last if $next->bacmp($limit) <= 0;
2920
2921 if ($sign == 0)
2922 {
2923 $x->badd($next);
2924 }
2925 else
2926 {
2927 $x->bsub($next);
2928 }
2929 $sign = 1-$sign; # alternate
2930 # calculate things for the next term
2931 $over->bmul($x2); # $x*$x
2932 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2933 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
2934 }
2935
2936 # shortcut to not run through _find_round_parameters again
2937 if (defined $params[0])
2938 {
2939 $x->bround($params[0],$params[2]); # then round accordingly
2940 }
2941 else
2942 {
2943 $x->bfround($params[1],$params[2]); # then round accordingly
2944 }
2945 if ($fallback)
2946 {
2947 # clear a/p after round, since user did not request it
2948 delete $x->{_a}; delete $x->{_p};
2949 }
2950 # restore globals
2951 $$abr = $ab; $$pbr = $pb;
2952 $x;
2953 }
2954
2955sub bsin
2956 {
2957 # Calculate a sinus of x.
2958 my ($self,$x,@r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
2959
2960 # taylor: x^3 x^5 x^7 x^9
2961 # sin = x - --- + --- - --- + --- ...
2962 # 3! 5! 7! 9!
2963
2964 # we need to limit the accuracy to protect against overflow
2965 my $fallback = 0;
2966 my ($scale,@params);
2967 ($x,@params) = $x->_find_round_parameters(@r);
2968
2969 # constant object or error in _find_round_parameters?
2970 return $x if $x->modify('bsin') || $x->is_nan();
2971
2972 return $x->bzero(@r) if $x->is_zero();
2973
2974 # no rounding at all, so must use fallback
2975 if (scalar @params == 0)
2976 {
2977 # simulate old behaviour
2978 $params[0] = $self->div_scale(); # and round to it as accuracy
2979 $params[1] = undef; # disable P
2980 $scale = $params[0]+4; # at least four more for proper round
2981 $params[2] = $r[2]; # round mode by caller or undef
2982 $fallback = 1; # to clear a/p afterwards
2983 }
2984 else
2985 {
2986 # the 4 below is empirical, and there might be cases where it is not
2987 # enough...
2988 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2989 }
2990
2991 # when user set globals, they would interfere with our calculation, so
2992 # disable them and later re-enable them
2993 no strict 'refs';
2994 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
2995 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
2996 # we also need to disable any set A or P on $x (_find_round_parameters took
2997 # them already into account), since these would interfere, too
2998 delete $x->{_a}; delete $x->{_p};
2999 # need to disable $upgrade in BigInt, to avoid deep recursion
3000 local $Math::BigInt::upgrade = undef;
3001
3002 my $last = 0;
3003 my $over = $x * $x; # X ^ 2
3004 my $x2 = $over->copy(); # X ^ 2; difference between terms
3005 $over->bmul($x); # X ^ 3 as starting value
3006 my $sign = 1; # start with -=
3007 my $below = $self->new(6); my $factorial = $self->new(4);
36ec1dfe 3008 delete $x->{_a}; delete $x->{_p};
60a1aa19
T
3009
3010 my $limit = $self->new("1E-". ($scale-1));
3011 #my $steps = 0;
3012 while (3 < 5)
3013 {
3014 # we calculate the next term, and add it to the last
3015 # when the next term is below our limit, it won't affect the outcome
3016 # anymore, so we stop:
3017 my $next = $over->copy()->bdiv($below,$scale);
3018 last if $next->bacmp($limit) <= 0;
3019
3020 if ($sign == 0)
3021 {
3022 $x->badd($next);
3023 }
3024 else
3025 {
3026 $x->bsub($next);
3027 }
3028 $sign = 1-$sign; # alternate
3029 # calculate things for the next term
3030 $over->bmul($x2); # $x*$x
3031 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
3032 $below->bmul($factorial); $factorial->binc(); # n*(n+1)
3033 }
3034
3035 # shortcut to not run through _find_round_parameters again
3036 if (defined $params[0])
3037 {
3038 $x->bround($params[0],$params[2]); # then round accordingly
3039 }
3040 else
3041 {
3042 $x->bfround($params[1],$params[2]); # then round accordingly
3043 }
3044 if ($fallback)
3045 {
3046 # clear a/p after round, since user did not request it
3047 delete $x->{_a}; delete $x->{_p};
3048 }
3049 # restore globals
3050 $$abr = $ab; $$pbr = $pb;
3051 $x;
3052 }
3053
20e2035c
T
3054sub batan2
3055 {
30afc38d 3056 # calculate arcus tangens of ($y/$x)
20e2035c
T
3057
3058 # set up parameters
30afc38d 3059 my ($self,$y,$x,@r) = (ref($_[0]),@_);
20e2035c
T
3060 # objectify is costly, so avoid it
3061 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1])))
3062 {
30afc38d 3063 ($self,$y,$x,@r) = objectify(2,@_);
20e2035c
T
3064 }
3065
30afc38d 3066 return $y if $y->modify('batan2');
20e2035c 3067
30afc38d 3068 return $y->bnan() if ($y->{sign} eq $nan) || ($x->{sign} eq $nan);
20e2035c 3069
30afc38d
T
3070 # Y X
3071 # 0 0 result is 0
3072 # 0 +x result is 0
0dceeee6
RGS
3073 # ? inf result is 0
3074 return $y->bzero(@r) if ($x->is_inf('+') && !$y->is_inf()) || ($y->is_zero() && $x->{sign} eq '+');
3075
3076 # Y X
3077 # != 0 -inf result is +- pi
3078 if ($x->is_inf() || $y->is_inf())
3079 {
3080 # calculate PI
3081 my $pi = $self->bpi(@r);
3082 if ($y->is_inf())
3083 {
3084 # upgrade to BigRat etc.
3085 return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
3086 if ($x->{sign} eq '-inf')
3087 {
3088 # calculate 3 pi/4
3089 $MBI->_mul($pi->{_m}, $MBI->_new(3));
3090 $MBI->_div($pi->{_m}, $MBI->_new(4));
3091 }
3092 elsif ($x->{sign} eq '+inf')
3093 {
3094 # calculate pi/4
3095 $MBI->_div($pi->{_m}, $MBI->_new(4));
3096 }
3097 else
3098 {
3099 # calculate pi/2
3100 $MBI->_div($pi->{_m}, $MBI->_new(2));
3101 }
3102 $y->{sign} = substr($y->{sign},0,1); # keep +/-
3103 }
3104 # modify $y in place
3105 $y->{_m} = $pi->{_m};
3106 $y->{_e} = $pi->{_e};
3107 $y->{_es} = $pi->{_es};
3108 # keep the sign of $y
3109 return $y;
3110 }
3111
3112 return $upgrade->new($y)->batan2($upgrade->new($x),@r) if defined $upgrade;
30afc38d
T
3113
3114 # Y X
3115 # 0 -x result is PI
3116 if ($y->is_zero())
20e2035c 3117 {
30afc38d
T
3118 # calculate PI
3119 my $pi = $self->bpi(@r);
0dceeee6 3120 # modify $y in place
30afc38d
T
3121 $y->{_m} = $pi->{_m};
3122 $y->{_e} = $pi->{_e};
3123 $y->{_es} = $pi->{_es};
3124 $y->{sign} = '+';
3125 return $y;
3126 }
3127
3128 # Y X
3129 # +y 0 result is PI/2
3130 # -y 0 result is -PI/2
0dceeee6 3131 if ($x->is_zero())
30afc38d
T
3132 {
3133 # calculate PI/2
3134 my $pi = $self->bpi(@r);
0dceeee6 3135 # modify $y in place
30afc38d
T
3136 $y->{_m} = $pi->{_m};
3137 $y->{_e} = $pi->{_e};
3138 $y->{_es} = $pi->{_es};
3139 # -y => -PI/2, +y => PI/2
30afc38d
T
3140 $MBI->_div($y->{_m}, $MBI->_new(2));
3141 return $y;
20e2035c
T
3142 }
3143
30afc38d
T
3144 # we need to limit the accuracy to protect against overflow
3145 my $fallback = 0;
3146 my ($scale,@params);
3147 ($y,@params) = $y->_find_round_parameters(@r);
3148
3149 # error in _find_round_parameters?
3150 return $y if $y->is_nan();
3151
3152 # no rounding at all, so must use fallback
3153 if (scalar @params == 0)
3154 {
3155 # simulate old behaviour
3156 $params[0] = $self->div_scale(); # and round to it as accuracy
3157 $params[1] = undef; # disable P
3158 $scale = $params[0]+4; # at least four more for proper round
3159 $params[2] = $r[2]; # round mode by caller or undef
3160 $fallback = 1; # to clear a/p afterwards
3161 }
3162 else
3163 {
3164 # the 4 below is empirical, and there might be cases where it is not
3165 # enough...
3166 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3167 }
3168
3169 # inlined is_one() && is_one('-')
3170 if ($MBI->_is_one($y->{_m}) && $MBI->_is_zero($y->{_e}))
3171 {
3172 # shortcut: 1 1 result is PI/4
3173 # inlined is_one() && is_one('-')
3174 if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
3175 {
3176 # 1,1 => PI/4
3177 my $pi_4 = $self->bpi( $scale - 3);
0dceeee6 3178 # modify $y in place
30afc38d
T
3179 $y->{_m} = $pi_4->{_m};
3180 $y->{_e} = $pi_4->{_e};
3181 $y->{_es} = $pi_4->{_es};
3182 # 1 1 => +
3183 # -1 1 => -
3184 # 1 -1 => -
3185 # -1 -1 => +
3186 $y->{sign} = $x->{sign} eq $y->{sign} ? '+' : '-';
3187 $MBI->_div($y->{_m}, $MBI->_new(4));
3188 return $y;
3189 }
3190 # shortcut: 1 int(X) result is _atan_inv(X)
3191
3192 # is integer
3193 if ($x->{_es} eq '+')
3194 {
3195 my $x1 = $MBI->_copy($x->{_m});
3196 $MBI->_lsft($x1, $x->{_e},10) unless $MBI->_is_zero($x->{_e});
3197
3198 my ($a,$b) = $self->_atan_inv($x1, $scale);
3199 my $y_sign = $y->{sign};
3200 # calculate A/B
3201 $y->bone(); $y->{_m} = $a; my $y_d = $self->bone(); $y_d->{_m} = $b;
3202 $y->bdiv($y_d, @r);
3203 $y->{sign} = $y_sign;
3204 return $y;
3205 }
3206 }
3207
3208 # handle all other cases
3209 # X Y
3210 # +x +y 0 to PI/2
3211 # -x +y PI/2 to PI
3212 # +x -y 0 to -PI/2
3213 # -x -y -PI/2 to -PI
3214
3215 my $y_sign = $y->{sign};
20e2035c
T
3216
3217 # divide $x by $y
30afc38d
T
3218 $y->bdiv($x, $scale) unless $x->is_one();
3219 $y->batan(@r);
20e2035c 3220
30afc38d
T
3221 # restore sign
3222 $y->{sign} = $y_sign;
20e2035c 3223
30afc38d 3224 $y;
20e2035c
T
3225 }
3226
60a1aa19
T
3227sub batan
3228 {
3229 # Calculate a arcus tangens of x.
3230 my ($x,@r) = @_;
3231 my $self = ref($x);
3232
3233 # taylor: x^3 x^5 x^7 x^9
3234 # atan = x - --- + --- - --- + --- ...
3235 # 3 5 7 9
3236
60a1aa19
T
3237 # we need to limit the accuracy to protect against overflow
3238 my $fallback = 0;
3239 my ($scale,@params);
3240 ($x,@params) = $x->_find_round_parameters(@r);
3241
3242 # constant object or error in _find_round_parameters?
3243 return $x if $x->modify('batan') || $x->is_nan();
3244
30afc38d
T
3245 if ($x->{sign} =~ /^[+-]inf\z/)
3246 {
3247 # +inf result is PI/2
3248 # -inf result is -PI/2
3249 # calculate PI/2
3250 my $pi = $self->bpi(@r);
3251 # modify $x in place
3252 $x->{_m} = $pi->{_m};
3253 $x->{_e} = $pi->{_e};
3254 $x->{_es} = $pi->{_es};
3255 # -y => -PI/2, +y => PI/2
3256 $x->{sign} = substr($x->{sign},0,1); # +inf => +
3257 $MBI->_div($x->{_m}, $MBI->_new(2));
3258 return $x;
3259 }
3260
3261 return $x->bzero(@r) if $x->is_zero();
3262
60a1aa19
T
3263 # no rounding at all, so must use fallback
3264 if (scalar @params == 0)
3265 {
3266 # simulate old behaviour
3267 $params[0] = $self->div_scale(); # and round to it as accuracy
3268 $params[1] = undef; # disable P
3269 $scale = $params[0]+4; # at least four more for proper round
3270 $params[2] = $r[2]; # round mode by caller or undef
3271 $fallback = 1; # to clear a/p afterwards
3272 }
3273 else
3274 {
3275 # the 4 below is empirical, and there might be cases where it is not
3276 # enough...
3277 $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3278 }
3279
30afc38d
T
3280 # 1 or -1 => PI/4
3281 # inlined is_one() && is_one('-')
3282 if ($MBI->_is_one($x->{_m}) && $MBI->_is_zero($x->{_e}))
3283 {
3284 my $pi = $self->bpi($scale - 3);
3285 # modify $x in place
3286 $x->{_m} = $pi->{_m};
3287 $x->{_e} = $pi->{_e};
3288 $x->{_es} = $pi->{_es};
3289 # leave the sign of $x alone (+1 => +PI/4, -1 => -PI/4)
3290 $MBI->_div($x->{_m}, $MBI->_new(4));
3291 return $x;
3292 }
3293
3294 # This series is only valid if -1 < x < 1, so for other x we need to
c97ef841 3295 # to calculate PI/2 - atan(1/x):
30afc38d
T
3296 my $one = $MBI->_new(1);
3297 my $pi = undef;
3298 if ($x->{_es} eq '+' && ($MBI->_acmp($x->{_m},$one) >= 0))
3299 {
3300 # calculate PI/2
3301 $pi = $self->bpi($scale - 3);
3302 $MBI->_div($pi->{_m}, $MBI->_new(2));
3303 # calculate 1/$x:
3304 my $x_copy = $x->copy();
3305 # modify $x in place
3306 $x->bone(); $x->bdiv($x_copy,$scale);
3307 }
3308
60a1aa19
T
3309 # when user set globals, they would interfere with our calculation, so
3310 # disable them and later re-enable them
3311 no strict 'refs';
3312 my $abr = "$self\::accuracy"; my $ab = $$abr; $$abr = undef;
3313 my $pbr = "$self\::precision"; my $pb = $$pbr; $$pbr = undef;
3314 # we also need to disable any set A or P on $x (_find_round_parameters took
3315 # them already into account), since these would interfere, too
3316 delete $x->{_a}; delete $x->{_p};
3317 # need to disable $upgrade in BigInt, to avoid deep recursion
3318 local $Math::BigInt::upgrade = undef;
3319
3320 my $last = 0;
3321 my $over = $x * $x; # X ^ 2
3322 my $x2 = $over->copy(); # X ^ 2; difference between terms
3323 $over->bmul($x); # X ^ 3 as starting value
3324 my $sign = 1; # start with -=
3325 my $below = $self->new(3);
3326 my $two = $self->new(2);
36ec1dfe 3327 delete $x->{_a}; delete $x->{_p};
60a1aa19
T
3328
3329 my $limit = $self->new("1E-". ($scale-1));
3330 #my $steps = 0;
3331 while (3 < 5)
3332 {
3333 # we calculate the next term, and add it to the last
3334 # when the next term is below our limit, it won't affect the outcome
3335 # anymore, so we stop:
3336 my $next = $over->copy()->bdiv($below,$scale);
3337 last if $next->bacmp($limit) <= 0;
3338
3339 if ($sign == 0)
3340 {
3341 $x->badd($next);
3342 }
3343 else
3344 {
3345 $x->bsub($next);
3346 }
3347 $sign = 1-$sign; # alternate
3348 # calculate things for the next term
3349 $over->bmul($x2); # $x*$x
3350 $below->badd($two); # n += 2
3351 }
3352
30afc38d
T
3353 if (defined $pi)
3354 {
3355 my $x_copy = $x->copy();
3356 # modify $x in place
3357 $x->{_m} = $pi->{_m};
3358 $x->{_e} = $pi->{_e};
3359 $x->{_es} = $pi->{_es};
3360 # PI/2 - $x
3361 $x->bsub($x_copy);
3362 }
3363
60a1aa19
T
3364 # shortcut to not run through _find_round_parameters again
3365 if (defined $params[0])
3366 {
3367 $x->bround($params[0],$params[2]); # then round accordingly
3368 }
3369 else
3370 {
3371 $x->bfround($params[1],$params[2]); # then round accordingly
3372 }
3373 if ($fallback)
3374 {
3375 # clear a/p after round, since user did not request it
3376 delete $x->{_a}; delete $x->{_p};
3377 }
3378 # restore globals
3379 $$abr = $ab; $$pbr = $pb;
3380 $x;
3381 }
3382
fdb4b05f 3383###############################################################################
58cde26e
JH
3384# rounding functions
3385
3386sub bfround
3387 {
3388 # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
3389 # $n == 0 means round to integer
3390 # expects and returns normalized numbers!
ee15d750 3391 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
a0d0e21e 3392
b68b7ab1
T
3393 my ($scale,$mode) = $x->_scale_p(@_);
3394 return $x if !defined $scale || $x->modify('bfround'); # no-op
58cde26e 3395
574bacfe 3396 # never round a 0, +-inf, NaN
61f5c3f5
T
3397 if ($x->is_zero())
3398 {
3399 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
3400 return $x;
3401 }
3402 return $x if $x->{sign} !~ /^[+-]$/;
58cde26e 3403
ee15d750
JH
3404 # don't round if x already has lower precision
3405 return $x if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p});
3406
3407 $x->{_p} = $scale; # remember round in any case
ef9466ea 3408 delete $x->{_a}; # and clear A
58cde26e
JH
3409 if ($scale < 0)
3410 {
58cde26e 3411 # round right from the '.'
f9a08e12 3412
9b924220 3413 return $x if $x->{_es} eq '+'; # e >= 0 => nothing to round
f9a08e12 3414
58cde26e 3415 $scale = -$scale; # positive for simplicity
9b924220 3416 my $len = $MBI->_len($x->{_m}); # length of mantissa
f9a08e12
JH
3417
3418 # the following poses a restriction on _e, but if _e is bigger than a
3419 # scalar, you got other problems (memory etc) anyway
9b924220 3420 my $dad = -(0+ ($x->{_es}.$MBI->_num($x->{_e}))); # digits after dot
58cde26e 3421 my $zad = 0; # zeros after dot
f9a08e12 3422 $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
9b924220 3423
ada8209b 3424 # print "scale $scale dad $dad zad $zad len $len\n";
58cde26e
JH
3425 # number bsstr len zad dad
3426 # 0.123 123e-3 3 0 3
3427 # 0.0123 123e-4 3 1 4
3428 # 0.001 1e-3 1 2 3
3429 # 1.23 123e-2 3 0 2
3430 # 1.2345 12345e-4 5 0 4
3431
3432 # do not round after/right of the $dad
3433 return $x if $scale > $dad; # 0.123, scale >= 3 => exit
3434
ee15d750
JH
3435 # round to zero if rounding inside the $zad, but not for last zero like:
3436 # 0.0065, scale -2, round last '0' with following '65' (scale == zad case)
3437 return $x->bzero() if $scale < $zad;
3438 if ($scale == $zad) # for 0.006, scale -3 and trunc
58cde26e 3439 {
b3abae2a 3440 $scale = -$len;
58cde26e
JH
3441 }
3442 else
3443 {
3444 # adjust round-point to be inside mantissa
3445 if ($zad != 0)
3446 {
3447 $scale = $scale-$zad;
3448 }
3449 else
3450 {
3451 my $dbd = $len - $dad; $dbd = 0 if $dbd < 0; # digits before dot
3452 $scale = $dbd+$scale;
3453 }
3454 }
a0d0e21e 3455 }
58cde26e
JH
3456 else
3457 {
f9a08e12
JH
3458 # round left from the '.'
3459
58cde26e 3460 # 123 => 100 means length(123) = 3 - $scale (2) => 1
a5f75d66 3461
9b924220 3462 my $dbt = $MBI->_len($x->{_m});
b3abae2a 3463 # digits before dot
9b924220 3464 my $dbd = $dbt + ($x->{_es} . $MBI->_num($x->{_e}));
b3abae2a
JH
3465 # should be the same, so treat it as this
3466 $scale = 1 if $scale == 0;
3467 # shortcut if already integer
3468 return $x if $scale == 1 && $dbt <= $dbd;
3469 # maximum digits before dot
3470 ++$dbd;
3471
3472 if ($scale > $dbd)
3473 {
3474 # not enough digits before dot, so round to zero
3475 return $x->bzero;
3476 }
3477 elsif ( $scale == $dbd )
3478 {
3479 # maximum
3480 $scale = -$dbt;
3481 }
58cde26e 3482 else
b3abae2a
JH
3483 {
3484 $scale = $dbd - $scale;
3485 }
a0d0e21e 3486 }
574bacfe 3487 # pass sign to bround for rounding modes '+inf' and '-inf'
ae161977 3488 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
9b924220
RGS
3489 $m->bround($scale,$mode);
3490 $x->{_m} = $m->{value}; # get our mantissa back
58cde26e
JH
3491 $x->bnorm();
3492 }
3493
3494sub bround
3495 {
3496 # accuracy: preserve $N digits, and overwrite the rest with 0's
ee15d750 3497 my $x = shift; my $self = ref($x) || $x; $x = $self->new(shift) if !ref($x);
9b924220 3498
990fb837
RGS
3499 if (($_[0] || 0) < 0)
3500 {
3501 require Carp; Carp::croak ('bround() needs positive accuracy');
3502 }
58cde26e 3503
b68b7ab1
T
3504 my ($scale,$mode) = $x->_scale_a(@_);
3505 return $x if !defined $scale || $x->modify('bround'); # no-op
61f5c3f5 3506
ee15d750
JH
3507 # scale is now either $x->{_a}, $accuracy, or the user parameter
3508 # test whether $x already has lower accuracy, do nothing in this case
3509 # but do round if the accuracy is the same, since a math operation might
3510 # want to round a number with A=5 to 5 digits afterwards again
b68b7ab1 3511 return $x if defined $x->{_a} && $x->{_a} < $scale;
58cde26e 3512
61f5c3f5 3513 # scale < 0 makes no sense
b68b7ab1 3514 # scale == 0 => keep all digits
61f5c3f5 3515 # never round a +-inf, NaN
b68b7ab1 3516 return $x if ($scale <= 0) || $x->{sign} !~ /^[+-]$/;
58cde26e 3517
b68b7ab1
T
3518 # 1: never round a 0
3519 # 2: if we should keep more digits than the mantissa has, do nothing
3520 if ($x->is_zero() || $MBI->_len($x->{_m}) <= $scale)
61f5c3f5
T
3521 {
3522 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
3523 return $x;
3524 }
f216259d 3525
58cde26e 3526 # pass sign to bround for '+inf' and '-inf' rounding modes
ae161977 3527 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
9b924220
RGS
3528
3529 $m->bround($scale,$mode); # round mantissa
3530 $x->{_m} = $m->{value}; # get our mantissa back
ee15d750 3531 $x->{_a} = $scale; # remember rounding
ef9466ea 3532 delete $x->{_p}; # and clear P
574bacfe 3533 $x->bnorm(); # del trailing zeros gen. by bround()
58cde26e
JH
3534 }
3535
3536sub bfloor
3537 {
6358232b 3538 # round towards minus infinity
ee15d750 3539 my ($self,$x,$a,$p,$r) = ref($_[0]) ? (ref($_[0]),@_) : objectify(1,@_);
58cde26e
JH
3540
3541 return $x if $x->modify('bfloor');
3542
3543 return $x if $x->{sign} !~ /^[+-]$/; # nan, +inf, -inf
3544
3545 # if $x has digits after dot
9b924220 3546 if ($x->{_es} eq '-')
58cde26e