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