File Coverage

blib/lib/Math/BigFloat.pm
Criterion Covered Total %
statement 2054 2528 81.2
branch 1260 1902 66.2
condition 629 1032 60.9
subroutine 135 170 79.4
pod 87 88 98.8
total 4165 5720 72.8


line stmt bran cond sub pod time code
1             package Math::BigFloat;
2              
3             #
4             # Mike grinned. 'Two down, infinity to go' - Mike Nostrus in 'Before and After'
5             #
6              
7             # The following hash values are used internally:
8             # sign : "+", "-", "+inf", "-inf", or "NaN" if not a number
9             # _m : mantissa ($LIB thingy)
10             # _es : sign of _e
11             # _e : exponent ($LIB thingy)
12             # _a : accuracy
13             # _p : precision
14              
15 43     43   1078645 use 5.006001;
  43         304  
16 43     43   274 use strict;
  43         79  
  43         1097  
17 43     43   227 use warnings;
  43         109  
  43         1618  
18              
19 43     43   285 use Carp qw< carp croak >;
  43         93  
  43         3151  
20 43     43   358 use Scalar::Util qw< blessed >;
  43         89  
  43         2428  
21 43     43   29976 use Math::BigInt qw< >;
  43         187  
  43         75487  
22              
23             our $VERSION = '1.999842';
24             $VERSION =~ tr/_//d;
25              
26             require Exporter;
27             our @ISA = qw/Math::BigInt/;
28             our @EXPORT_OK = qw/bpi/;
29              
30             # $_trap_inf/$_trap_nan are internal and should never be accessed from outside
31             our ($AUTOLOAD, $accuracy, $precision, $div_scale, $round_mode, $rnd_mode,
32             $upgrade, $downgrade, $_trap_nan, $_trap_inf);
33              
34             use overload
35              
36             # overload key: with_assign
37              
38 13     13   183 '+' => sub { $_[0] -> copy() -> badd($_[1]); },
39              
40 73     73   1324 '-' => sub { my $c = $_[0] -> copy();
41 73 50       337 $_[2] ? $c -> bneg() -> badd($_[1])
42             : $c -> bsub($_[1]); },
43              
44 153     153   5156 '*' => sub { $_[0] -> copy() -> bmul($_[1]); },
45              
46 12 50   12   2835 '/' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bdiv($_[0])
47             : $_[0] -> copy() -> bdiv($_[1]); },
48              
49 0 0   0   0 '%' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bmod($_[0])
50             : $_[0] -> copy() -> bmod($_[1]); },
51              
52 7 50   7   1139 '**' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bpow($_[0])
53             : $_[0] -> copy() -> bpow($_[1]); },
54              
55 0 0   0   0 '<<' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bblsft($_[0])
56             : $_[0] -> copy() -> bblsft($_[1]); },
57              
58 0 0   0   0 '>>' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bbrsft($_[0])
59             : $_[0] -> copy() -> bbrsft($_[1]); },
60              
61             # overload key: assign
62              
63 29     29   320 '+=' => sub { $_[0] -> badd($_[1]); },
64              
65 28     28   310 '-=' => sub { $_[0] -> bsub($_[1]); },
66              
67 6424     6424   17674 '*=' => sub { $_[0] -> bmul($_[1]); },
68              
69 8     8   105 '/=' => sub { scalar $_[0] -> bdiv($_[1]); },
70              
71 8     8   111 '%=' => sub { $_[0] -> bmod($_[1]); },
72              
73 4     4   75 '**=' => sub { $_[0] -> bpow($_[1]); },
74              
75 8     8   111 '<<=' => sub { $_[0] -> bblsft($_[1]); },
76              
77 8     8   1515 '>>=' => sub { $_[0] -> bbrsft($_[1]); },
78              
79             # 'x=' => sub { },
80              
81             # '.=' => sub { },
82              
83             # overload key: num_comparison
84              
85 301 50   301   1460 '<' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> blt($_[0])
86             : $_[0] -> blt($_[1]); },
87              
88 857 50   857   3059 '<=' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> ble($_[0])
89             : $_[0] -> ble($_[1]); },
90              
91 879 50   879   3910 '>' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bgt($_[0])
92             : $_[0] -> bgt($_[1]); },
93              
94 177 100   177   834 '>=' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bge($_[0])
95             : $_[0] -> bge($_[1]); },
96              
97 186     186   5938 '==' => sub { $_[0] -> beq($_[1]); },
98              
99 9     9   366 '!=' => sub { $_[0] -> bne($_[1]); },
100              
101             # overload key: 3way_comparison
102              
103 0     0   0 '<=>' => sub { my $cmp = $_[0] -> bcmp($_[1]);
104 0 0 0     0 defined($cmp) && $_[2] ? -$cmp : $cmp; },
105              
106 5923 50   5923   1501146 'cmp' => sub { $_[2] ? "$_[1]" cmp $_[0] -> bstr()
107             : $_[0] -> bstr() cmp "$_[1]"; },
108              
109             # overload key: str_comparison
110              
111             # 'lt' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrlt($_[0])
112             # : $_[0] -> bstrlt($_[1]); },
113             #
114             # 'le' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrle($_[0])
115             # : $_[0] -> bstrle($_[1]); },
116             #
117             # 'gt' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrgt($_[0])
118             # : $_[0] -> bstrgt($_[1]); },
119             #
120             # 'ge' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bstrge($_[0])
121             # : $_[0] -> bstrge($_[1]); },
122             #
123             # 'eq' => sub { $_[0] -> bstreq($_[1]); },
124             #
125             # 'ne' => sub { $_[0] -> bstrne($_[1]); },
126              
127             # overload key: binary
128              
129 0 0   0   0 '&' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> band($_[0])
130             : $_[0] -> copy() -> band($_[1]); },
131              
132 0     0   0 '&=' => sub { $_[0] -> band($_[1]); },
133              
134 0 0   0   0 '|' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bior($_[0])
135             : $_[0] -> copy() -> bior($_[1]); },
136              
137 0     0   0 '|=' => sub { $_[0] -> bior($_[1]); },
138              
139 0 0   0   0 '^' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bxor($_[0])
140             : $_[0] -> copy() -> bxor($_[1]); },
141              
142 0     0   0 '^=' => sub { $_[0] -> bxor($_[1]); },
143              
144             # '&.' => sub { },
145              
146             # '&.=' => sub { },
147              
148             # '|.' => sub { },
149              
150             # '|.=' => sub { },
151              
152             # '^.' => sub { },
153              
154             # '^.=' => sub { },
155              
156             # overload key: unary
157              
158 352     352   1005 'neg' => sub { $_[0] -> copy() -> bneg(); },
159              
160             # '!' => sub { },
161              
162 0     0   0 '~' => sub { $_[0] -> copy() -> bnot(); },
163              
164             # '~.' => sub { },
165              
166             # overload key: mutators
167              
168 10     10   126 '++' => sub { $_[0] -> binc() },
169              
170 0     0   0 '--' => sub { $_[0] -> bdec() },
171              
172             # overload key: func
173              
174 0 0   0   0 'atan2' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> batan2($_[0])
175             : $_[0] -> copy() -> batan2($_[1]); },
176              
177 0     0   0 'cos' => sub { $_[0] -> copy() -> bcos(); },
178              
179 0     0   0 'sin' => sub { $_[0] -> copy() -> bsin(); },
180              
181 0     0   0 'exp' => sub { $_[0] -> copy() -> bexp($_[1]); },
182              
183 0     0   0 'abs' => sub { $_[0] -> copy() -> babs(); },
184              
185 40     40   549 'log' => sub { $_[0] -> copy() -> blog(); },
186              
187 0     0   0 'sqrt' => sub { $_[0] -> copy() -> bsqrt(); },
188              
189 140     140   412 'int' => sub { $_[0] -> copy() -> bint(); },
190              
191             # overload key: conversion
192              
193 280 50   280   625 'bool' => sub { $_[0] -> is_zero() ? '' : 1; },
194              
195 742     742   80743 '""' => sub { $_[0] -> bstr(); },
196              
197 8     8   31 '0+' => sub { $_[0] -> numify(); },
198              
199 0     0   0 '=' => sub { $_[0] -> copy(); },
200              
201 43     43   422 ;
  43         101  
  43         3863  
202              
203             ##############################################################################
204             # global constants, flags and assorted stuff
205              
206             # the following are public, but their usage is not recommended. Use the
207             # accessor methods instead.
208              
209             # class constants, use Class->constant_name() to access
210             # one of 'even', 'odd', '+inf', '-inf', 'zero', 'trunc' or 'common'
211             $round_mode = 'even';
212             $accuracy = undef;
213             $precision = undef;
214             $div_scale = 40;
215              
216             $upgrade = undef;
217             $downgrade = undef;
218             # the package we are using for our private parts, defaults to:
219             # Math::BigInt->config('lib')
220             my $LIB = 'Math::BigInt::Calc';
221              
222             # are NaNs ok? (otherwise it dies when encountering an NaN) set w/ config()
223             $_trap_nan = 0;
224             # the same for infinity
225             $_trap_inf = 0;
226              
227             # constant for easier life
228             my $nan = 'NaN';
229              
230             my $IMPORT = 0; # was import() called yet? used to make require work
231              
232             # some digits of accuracy for blog(undef, 10); which we use in blog() for speed
233             my $LOG_10 =
234             '2.3025850929940456840179914546843642076011014886287729760333279009675726097';
235             my $LOG_10_A = length($LOG_10)-1;
236             # ditto for log(2)
237             my $LOG_2 =
238             '0.6931471805599453094172321214581765680755001343602552541206800094933936220';
239             my $LOG_2_A = length($LOG_2)-1;
240             my $HALF = '0.5'; # made into an object if nec.
241              
242             ##############################################################################
243             # the old code had $rnd_mode, so we need to support it, too
244              
245             sub TIESCALAR {
246 43     43   146 my ($class) = @_;
247 43         239 bless \$round_mode, $class;
248             }
249              
250             sub FETCH {
251 1     1   607 return $round_mode;
252             }
253              
254             sub STORE {
255 1     1   651 $rnd_mode = $_[0]->round_mode($_[1]);
256             }
257              
258             BEGIN {
259             # when someone sets $rnd_mode, we catch this and check the value to see
260             # whether it is valid or not.
261 43     43   33632 $rnd_mode = 'even';
262 43         306 tie $rnd_mode, 'Math::BigFloat';
263              
264 43         4288 *as_number = \&as_int;
265             }
266              
267       0     sub DESTROY {
268             # going through AUTOLOAD for every DESTROY is costly, avoid it by empty sub
269             }
270              
271             sub AUTOLOAD {
272             # make fxxx and bxxx both work by selectively mapping fxxx() to MBF::bxxx()
273 3     3   9 my $name = $AUTOLOAD;
274 3         42 $name =~ s/(.*):://; # split package
275 3   50     27 my $c = $1 || __PACKAGE__;
276 43     43   338 no strict 'refs';
  43         98  
  43         175009  
277 3 50       23 $c->import() if $IMPORT == 0;
278 3 50       18 if (!_method_alias($name)) {
279 3 50       11 if (!defined $name) {
280             # delayed load of Carp and avoid recursion
281 0         0 croak("$c: Can't call a method without name");
282             }
283 3 50       11 if (!_method_hand_up($name)) {
284             # delayed load of Carp and avoid recursion
285 0         0 croak("Can't call $c\-\>$name, not a valid method");
286             }
287             # try one level up, but subst. bxxx() for fxxx() since MBI only got
288             # bxxx()
289 3         9 $name =~ s/^f/b/;
290 3         6 return &{"Math::BigInt"."::$name"}(@_);
  3         25  
291             }
292 0         0 my $bname = $name;
293 0         0 $bname =~ s/^f/b/;
294 0         0 $c .= "::$name";
295 0         0 *{$c} = \&{$bname};
  0         0  
  0         0  
296 0         0 &{$c}; # uses @_
  0         0  
297             }
298              
299             ##############################################################################
300              
301             {
302             # valid method aliases for AUTOLOAD
303             my %methods = map { $_ => 1 }
304             qw / fadd fsub fmul fdiv fround ffround fsqrt fmod fstr fsstr fpow fnorm
305             fint facmp fcmp fzero fnan finf finc fdec ffac fneg
306             fceil ffloor frsft flsft fone flog froot fexp
307             /;
308             # valid methods that can be handed up (for AUTOLOAD)
309             my %hand_ups = map { $_ => 1 }
310             qw / is_nan is_inf is_negative is_positive is_pos is_neg
311             accuracy precision div_scale round_mode fabs fnot
312             objectify upgrade downgrade
313             bone binf bnan bzero
314             bsub
315             /;
316              
317 3   50 3   23 sub _method_alias { exists $methods{$_[0]||''}; }
318 3   50 3   17 sub _method_hand_up { exists $hand_ups{$_[0]||''}; }
319             }
320              
321             sub isa {
322 24688     24688 0 49623 my ($self, $class) = @_;
323 24688 100       59363 return if $class =~ /^Math::BigInt/; # we aren't one of these
324 23557         93255 UNIVERSAL::isa($self, $class);
325             }
326              
327             sub config {
328             # return (later set?) configuration data as hash ref
329 55   50 55 1 42881 my $class = shift || 'Math::BigFloat';
330              
331             # Getter/accessor.
332              
333 55 100 100     317 if (@_ == 1 && ref($_[0]) ne 'HASH') {
334 23         49 my $param = shift;
335 23 100       70 return $class if $param eq 'class';
336 21 100       95 return $LIB if $param eq 'with';
337 16         175 return $class->SUPER::config($param);
338             }
339              
340             # Setter.
341              
342 32         112 my $cfg = $class->SUPER::config(@_);
343              
344             # now we need only to override the ones that are different from our parent
345 31         54 $cfg->{class} = $class;
346 31         51 $cfg->{with} = $LIB;
347 31         108 $cfg;
348             }
349              
350             ###############################################################################
351             # Constructor methods
352             ###############################################################################
353              
354             sub new {
355             # Create a new Math::BigFloat object from a string or another bigfloat
356             # object.
357             # _e: exponent
358             # _m: mantissa
359             # sign => ("+", "-", "+inf", "-inf", or "NaN")
360              
361 17120     17120 1 6041553 my $self = shift;
362 17120         30419 my $selfref = ref $self;
363 17120   33     58717 my $class = $selfref || $self;
364              
365             # Make "require" work.
366              
367 17120 100       38070 $class -> import() if $IMPORT == 0;
368              
369             # Although this use has been discouraged for more than 10 years, people
370             # apparently still use it, so we still support it.
371              
372 17120 100       37263 return $class -> bzero() unless @_;
373              
374 17112         35994 my ($wanted, @r) = @_;
375              
376 17112 50       34966 if (!defined($wanted)) {
377             #if (warnings::enabled("uninitialized")) {
378             # warnings::warn("uninitialized",
379             # "Use of uninitialized value in new()");
380             #}
381 0         0 return $class -> bzero(@r);
382             }
383              
384 17112 50 66     60728 if (!ref($wanted) && $wanted eq "") {
385             #if (warnings::enabled("numeric")) {
386             # warnings::warn("numeric",
387             # q|Argument "" isn't numeric in new()|);
388             #}
389             #return $class -> bzero(@r);
390 0         0 return $class -> bnan(@r);
391             }
392              
393             # Initialize a new object.
394              
395 17112 50       49018 $self = bless {}, $class unless $selfref;
396              
397             # Math::BigFloat or subclass
398              
399 17112 100 100     62280 if (defined(blessed($wanted)) && $wanted -> isa(__PACKAGE__)) {
400              
401             # Don't copy the accuracy and precision, because a new object should get
402             # them from the global configuration.
403              
404 344         974 $self -> {sign} = $wanted -> {sign};
405 344         1054 $self -> {_m} = $LIB -> _copy($wanted -> {_m});
406 344         741 $self -> {_es} = $wanted -> {_es};
407 344         824 $self -> {_e} = $LIB -> _copy($wanted -> {_e});
408 344 50 66     1586 $self = $self->round(@r)
      66        
409             unless @r >= 2 && !defined($r[0]) && !defined($r[1]);
410 344         1922 return $self;
411             }
412              
413             # Shortcut for Math::BigInt and its subclasses. This should be improved.
414              
415 16768 100       41487 if (defined(blessed($wanted))) {
416 523 100       1625 if ($wanted -> isa('Math::BigInt')) {
417 522         1687 $self->{sign} = $wanted -> {sign};
418 522         1656 $self->{_m} = $LIB -> _copy($wanted -> {value});
419 522         1062 $self->{_es} = '+';
420 522         1282 $self->{_e} = $LIB -> _zero();
421 522         1491 return $self -> bnorm();
422             }
423              
424 1 50       9 if ($wanted -> can("as_number")) {
425 0         0 $self->{sign} = $wanted -> sign();
426 0         0 $self->{_m} = $wanted -> as_number() -> {value};
427 0         0 $self->{_es} = '+';
428 0         0 $self->{_e} = $LIB -> _zero();
429 0         0 return $self -> bnorm();
430             }
431             }
432              
433             # Shortcut for simple forms like '123' that have no trailing zeros. Trailing
434             # zeros would require a non-zero exponent.
435              
436 16246 100       86460 if ($wanted =~
437             / ^
438             \s* # optional leading whitespace
439             ( [+-]? ) # optional sign
440             0* # optional leading zeros
441             ( [1-9] (?: [0-9]* [1-9] )? ) # significand
442             \s* # optional trailing whitespace
443             $
444             /x)
445             {
446 7885 100       17534 return $downgrade -> new($1 . $2) if defined $downgrade;
447 7877   100     35870 $self->{sign} = $1 || '+';
448 7877         27666 $self->{_m} = $LIB -> _new($2);
449 7877         15330 $self->{_es} = '+';
450 7877         19146 $self->{_e} = $LIB -> _zero();
451 7877 100 100     35569 $self = $self->round(@r)
      100        
452             unless @r >= 2 && !defined $r[0] && !defined $r[1];
453 7877         72319 return $self;
454             }
455              
456             # Handle Infs.
457              
458 8361 100       25382 if ($wanted =~ / ^
459             \s*
460             ( [+-]? )
461             inf (?: inity )?
462             \s*
463             \z
464             /ix)
465             {
466 1727   100     6775 my $sgn = $1 || '+';
467 1727         5641 return $class -> binf($sgn, @r);
468             }
469              
470             # Handle explicit NaNs (not the ones returned due to invalid input).
471              
472 6634 100       16852 if ($wanted =~ / ^
473             \s*
474             ( [+-]? )
475             nan
476             \s*
477             \z
478             /ix)
479             {
480 475         1868 return $class -> bnan(@r);
481             }
482              
483 6159         10232 my @parts;
484              
485 6159 100 66     63565 if (
      33        
      66        
      66        
      66        
      100        
      33        
      66        
486             # Handle hexadecimal numbers. We auto-detect hexadecimal numbers if they
487             # have a "0x", "0X", "x", or "X" prefix, cf. CORE::oct().
488              
489             $wanted =~ /^\s*[+-]?0?[Xx]/ and
490             @parts = $class -> _hex_str_to_flt_lib_parts($wanted)
491              
492             or
493              
494             # Handle octal numbers. We auto-detect octal numbers if they have a
495             # "0o", "0O", "o", "O" prefix, cf. CORE::oct().
496              
497             $wanted =~ /^\s*[+-]?0?[Oo]/ and
498             @parts = $class -> _oct_str_to_flt_lib_parts($wanted)
499              
500             or
501              
502             # Handle binary numbers. We auto-detect binary numbers if they have a
503             # "0b", "0B", "b", or "B" prefix, cf. CORE::oct().
504              
505             $wanted =~ /^\s*[+-]?0?[Bb]/ and
506             @parts = $class -> _bin_str_to_flt_lib_parts($wanted)
507              
508             or
509              
510             # At this point, what is left are decimal numbers that aren't handled
511             # above and octal floating point numbers that don't have any of the
512             # "0o", "0O", "o", or "O" prefixes. First see if it is a decimal number.
513              
514             @parts = $class -> _dec_str_to_flt_lib_parts($wanted)
515             or
516              
517             # See if it is an octal floating point number. The extra check is
518             # included because _oct_str_to_flt_lib_parts() accepts octal numbers
519             # that don't have a prefix (this is needed to make it work with, e.g.,
520             # from_oct() that don't require a prefix). However, Perl requires a
521             # prefix for octal floating point literals. For example, "1p+0" is not
522             # valid, but "01p+0" and "0__1p+0" are.
523              
524             $wanted =~ /^\s*[+-]?0_*\d/ and
525             @parts = $class -> _oct_str_to_flt_lib_parts($wanted))
526             {
527 5708         23260 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts;
528              
529 5708 100 100     25211 $self = $self->round(@r)
      100        
530             unless @r >= 2 && !defined($r[0]) && !defined($r[1]);
531              
532 5708 100 100     15019 return $downgrade -> new($self -> bdstr(), @r)
533             if defined($downgrade) && $self -> is_int();
534 5706         60451 return $self;
535             }
536              
537             # If we get here, the value is neither a valid decimal, binary, octal, or
538             # hexadecimal number. It is not an explicit Inf or a NaN either.
539              
540 451         1390 return $class -> bnan(@r);
541             }
542              
543             sub from_dec {
544 1     1 1 2549 my $self = shift;
545 1         3 my $selfref = ref $self;
546 1   33     7 my $class = $selfref || $self;
547              
548             # Don't modify constant (read-only) objects.
549              
550 1 50 33     4 return $self if $selfref && $self->modify('from_dec');
551              
552 1         2 my $str = shift;
553 1         3 my @r = @_;
554              
555             # If called as a class method, initialize a new object.
556              
557 1 50       6 $self = bless {}, $class unless $selfref;
558              
559 1 50       5 if (my @parts = $class -> _dec_str_to_flt_lib_parts($str)) {
560 1         7 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts;
561              
562 1 0 33     13 $self = $self->round(@r)
      33        
563             unless @r >= 2 && !defined($r[0]) && !defined($r[1]);
564              
565 1 50 33     6 return $downgrade -> new($self -> bdstr(), @r)
566             if defined($downgrade) && $self -> is_int();
567 0         0 return $self;
568             }
569              
570 0         0 return $self -> bnan(@r);
571             }
572              
573             sub from_hex {
574 1     1 1 2540 my $self = shift;
575 1         2 my $selfref = ref $self;
576 1   33     8 my $class = $selfref || $self;
577              
578             # Don't modify constant (read-only) objects.
579              
580 1 50 33     5 return $self if $selfref && $self->modify('from_hex');
581              
582 1         3 my $str = shift;
583 1         2 my @r = @_;
584              
585             # If called as a class method, initialize a new object.
586              
587 1 50       5 $self = bless {}, $class unless $selfref;
588              
589 1 50       12 if (my @parts = $class -> _hex_str_to_flt_lib_parts($str)) {
590 1         14 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts;
591              
592 1 0 33     10 $self = $self->round(@r)
      33        
593             unless @r >= 2 && !defined($r[0]) && !defined($r[1]);
594              
595 1 50 33     6 return $downgrade -> new($self -> bdstr(), @r)
596             if defined($downgrade) && $self -> is_int();
597 0         0 return $self;
598             }
599              
600 0         0 return $self -> bnan(@r);
601             }
602              
603             sub from_oct {
604 1     1 1 2585 my $self = shift;
605 1         5 my $selfref = ref $self;
606 1   33     7 my $class = $selfref || $self;
607              
608             # Don't modify constant (read-only) objects.
609              
610 1 50 33     5 return $self if $selfref && $self->modify('from_oct');
611              
612 1         2 my $str = shift;
613 1         3 my @r = @_;
614              
615             # If called as a class method, initialize a new object.
616              
617 1 50       5 $self = bless {}, $class unless $selfref;
618              
619 1 50       11 if (my @parts = $class -> _oct_str_to_flt_lib_parts($str)) {
620 1         6 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts;
621              
622 1 0 33     8 $self = $self->round(@r)
      33        
623             unless @r >= 2 && !defined($r[0]) && !defined($r[1]);
624              
625 1 50 33     10 return $downgrade -> new($self -> bdstr(), @r)
626             if defined($downgrade) && $self -> is_int();
627 0         0 return $self;
628             }
629              
630 0         0 return $self -> bnan(@r);
631             }
632              
633             sub from_bin {
634 3     3 1 2629 my $self = shift;
635 3         7 my $selfref = ref $self;
636 3   33     13 my $class = $selfref || $self;
637              
638             # Don't modify constant (read-only) objects.
639              
640 3 50 33     10 return $self if $selfref && $self->modify('from_bin');
641              
642 3         7 my $str = shift;
643 3         7 my @r = @_;
644              
645             # If called as a class method, initialize a new object.
646              
647 3 50       11 $self = bless {}, $class unless $selfref;
648              
649 3 50       14 if (my @parts = $class -> _bin_str_to_flt_lib_parts($str)) {
650 3         16 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts;
651              
652 3 0 33     16 $self = $self->round(@r)
      33        
653             unless @r >= 2 && !defined($r[0]) && !defined($r[1]);
654              
655 3 50 33     13 return $downgrade -> new($self -> bdstr(), @r)
656             if defined($downgrade) && $self -> is_int();
657 0         0 return $self;
658             }
659              
660 0         0 return $self -> bnan(@r);
661             }
662              
663             sub from_ieee754 {
664 1     1 1 2573 my $self = shift;
665 1         6 my $selfref = ref $self;
666 1   33     6 my $class = $selfref || $self;
667              
668             # Don't modify constant (read-only) objects.
669              
670 1 50 33     7 return $self if $selfref && $self->modify('from_ieee754');
671              
672 1         3 my $in = shift; # input string (or raw bytes)
673 1         2 my $format = shift; # format ("binary32", "decimal64" etc.)
674 1         4 my $enc; # significand encoding (applies only to decimal)
675             my $k; # storage width in bits
676 1         0 my $b; # base
677 1         4 my @r = @_; # rounding parameters, if any
678              
679 1 50       9 if ($format =~ /^binary(\d+)\z/) {
    0          
    0          
    0          
    0          
    0          
    0          
    0          
680 1         3 $k = $1;
681 1         5 $b = 2;
682             } elsif ($format =~ /^decimal(\d+)(dpd|bcd)?\z/) {
683 0         0 $k = $1;
684 0         0 $b = 10;
685 0   0     0 $enc = $2 || 'dpd'; # default is dencely-packed decimals (DPD)
686             } elsif ($format eq 'half') {
687 0         0 $k = 16;
688 0         0 $b = 2;
689             } elsif ($format eq 'single') {
690 0         0 $k = 32;
691 0         0 $b = 2;
692             } elsif ($format eq 'double') {
693 0         0 $k = 64;
694 0         0 $b = 2;
695             } elsif ($format eq 'quadruple') {
696 0         0 $k = 128;
697 0         0 $b = 2;
698             } elsif ($format eq 'octuple') {
699 0         0 $k = 256;
700 0         0 $b = 2;
701             } elsif ($format eq 'sexdecuple') {
702 0         0 $k = 512;
703 0         0 $b = 2;
704             }
705              
706 1 50       4 if ($b == 2) {
707              
708             # Get the parameters for this format.
709              
710 1         12 my $p; # precision (in bits)
711             my $t; # number of bits in significand
712 1         0 my $w; # number of bits in exponent
713              
714 1 50       17 if ($k == 16) { # binary16 (half-precision)
    50          
    0          
715 0         0 $p = 11;
716 0         0 $t = 10;
717 0         0 $w = 5;
718             } elsif ($k == 32) { # binary32 (single-precision)
719 1         2 $p = 24;
720 1         2 $t = 23;
721 1         2 $w = 8;
722             } elsif ($k == 64) { # binary64 (double-precision)
723 0         0 $p = 53;
724 0         0 $t = 52;
725 0         0 $w = 11;
726             } else { # binaryN (quadruple-precision and above)
727 0 0 0     0 if ($k < 128 || $k != 32 * sprintf('%.0f', $k / 32)) {
728 0         0 croak "Number of bits must be 16, 32, 64, or >= 128 and",
729             " a multiple of 32";
730             }
731 0         0 $p = $k - sprintf('%.0f', 4 * log($k) / log(2)) + 13;
732 0         0 $t = $p - 1;
733 0         0 $w = $k - $t - 1;
734             }
735              
736             # The maximum exponent, minimum exponent, and exponent bias.
737              
738 1         5 my $emax = Math::BigFloat -> new(2) -> bpow($w - 1) -> bdec();
739 1         24 my $emin = 1 - $emax;
740 1         196 my $bias = $emax;
741              
742             # Undefined input.
743              
744 1 50       11 unless (defined $in) {
745 0         0 carp("Input is undefined");
746 0         0 return $self -> bzero(@r);
747             }
748              
749             # Make sure input string is a string of zeros and ones.
750              
751 1         2 my $len = CORE::length $in;
752 1 50       5 if (8 * $len == $k) { # bytes
    0          
    0          
753 1         6 $in = unpack "B*", $in;
754             } elsif (4 * $len == $k) { # hexadecimal
755 0 0       0 if ($in =~ /([^\da-f])/i) {
756 0         0 croak "Illegal hexadecimal digit '$1'";
757             }
758 0         0 $in = unpack "B*", pack "H*", $in;
759             } elsif ($len == $k) { # bits
760 0 0       0 if ($in =~ /([^01])/) {
761 0         0 croak "Illegal binary digit '$1'";
762             }
763             } else {
764 0         0 croak "Unknown input -- $in";
765             }
766              
767             # Split bit string into sign, exponent, and mantissa/significand.
768              
769 1 50       5 my $sign = substr($in, 0, 1) eq '1' ? '-' : '+';
770 1         6 my $expo = $class -> from_bin(substr($in, 1, $w));
771 1         19 my $mant = $class -> from_bin(substr($in, $w + 1));
772              
773 1         4 my $x;
774              
775 1         17 $expo = $expo -> bsub($bias); # subtract bias
776              
777 1 50       5 if ($expo < $emin) { # zero and subnormals
    50          
778 0 0       0 if ($mant == 0) { # zero
779 0         0 $x = $class -> bzero();
780             } else { # subnormals
781             # compute (1/$b)**(N) rather than ($b)**(-N)
782 0         0 $x = $class -> new("0.5"); # 1/$b
783 0         0 $x = $x -> bpow($bias + $t - 1) -> bmul($mant);
784 0 0       0 $x = $x -> bneg() if $sign eq '-';
785             }
786             }
787              
788             elsif ($expo > $emax) { # inf and nan
789 0 0       0 if ($mant == 0) { # inf
790 0         0 $x = $class -> binf($sign);
791             } else { # nan
792 0         0 $x = $class -> bnan(@r);
793             }
794             }
795              
796             else { # normals
797 1         15 $mant = $class -> new(2) -> bpow($t) -> badd($mant);
798 1 50       6 if ($expo < $t) {
799             # compute (1/$b)**(N) rather than ($b)**(-N)
800 1         4 $x = $class -> new("0.5"); # 1/$b
801 1         8 $x = $x -> bpow($t - $expo) -> bmul($mant);
802             } else {
803 0         0 $x = $class -> new(2);
804 0         0 $x = $x -> bpow($expo - $t) -> bmul($mant);
805             }
806 1 50       22 $x = $x -> bneg() if $sign eq '-';
807             }
808              
809 1 50       6 if ($selfref) {
810 0         0 $self -> {sign} = $x -> {sign};
811 0         0 $self -> {_m} = $x -> {_m};
812 0         0 $self -> {_es} = $x -> {_es};
813 0         0 $self -> {_e} = $x -> {_e};
814             } else {
815 1         3 $self = $x;
816             }
817              
818 1 50 33     23 return $downgrade -> new($self -> bdstr(), @r)
819             if defined($downgrade) && $self -> is_int();
820 0         0 return $self -> round(@r);
821             }
822              
823 0         0 croak("The format '$format' is not yet supported.");
824             }
825              
826             sub bzero {
827             # create/assign '+0'
828              
829             # Class::method(...) -> Class->method(...)
830 603 50 66 603 1 20197 unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) ||
      66        
831             $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i))
832             {
833             #carp "Using ", (caller(0))[3], "() as a function is deprecated;",
834             # " use is as a method instead";
835 0         0 unshift @_, __PACKAGE__;
836             }
837              
838 603         1538 my $self = shift;
839 603         1150 my $selfref = ref $self;
840 603   66     1539 my $class = $selfref || $self;
841              
842 603 100       1819 $self->import() if $IMPORT == 0; # make require work
843              
844             # Don't modify constant (read-only) objects.
845              
846 603 50 66     2431 return $self if $selfref && $self->modify('bzero');
847              
848             # Get the rounding parameters, if any.
849              
850 603         1246 my @r = @_;
851              
852 603 100       1391 return $downgrade -> bzero(@r) if defined $downgrade;
853              
854             # If called as a class method, initialize a new object.
855              
856 602 100       1390 $self = bless {}, $class unless $selfref;
857              
858 602         1323 $self -> {sign} = '+';
859 602         1695 $self -> {_m} = $LIB -> _zero();
860 602         1289 $self -> {_es} = '+';
861 602         1465 $self -> {_e} = $LIB -> _zero();
862              
863             # If rounding parameters are given as arguments, use them. If no rounding
864             # parameters are given, and if called as a class method initialize the new
865             # instance with the class variables.
866              
867             #return $self -> round(@r); # this should work, but doesnt; fixme!
868              
869 602 100       1689 if (@r) {
870 64 50 100     287 croak "can't specify both accuracy and precision"
      66        
871             if @r >= 2 && defined($r[0]) && defined($r[1]);
872 64         167 $self->{_a} = $r[0];
873 64         141 $self->{_p} = $r[1];
874             } else {
875 538 100       1344 unless($selfref) {
876 101         320 $self->{_a} = $class -> accuracy();
877 101         357 $self->{_p} = $class -> precision();
878             }
879             }
880              
881 602         4467 return $self;
882             }
883              
884             sub bone {
885             # Create or assign '+1' (or -1 if given sign '-').
886              
887             # Class::method(...) -> Class->method(...)
888 1589 50 66 1589 1 24050 unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) ||
      66        
889             $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i))
890             {
891             #carp "Using ", (caller(0))[3], "() as a function is deprecated;",
892             # " use is as a method instead";
893 0         0 unshift @_, __PACKAGE__;
894             }
895              
896 1589         3561 my $self = shift;
897 1589         2855 my $selfref = ref $self;
898 1589   66     4356 my $class = $selfref || $self;
899              
900 1589 100       3465 $self->import() if $IMPORT == 0; # make require work
901              
902             # Don't modify constant (read-only) objects.
903              
904 1589 50 66     5096 return $self if $selfref && $self->modify('bone');
905              
906 1589 100       3238 return $downgrade -> bone(@_) if defined $downgrade;
907              
908             # Get the sign.
909              
910 1588         2670 my $sign = '+'; # default is to return +1
911 1588 100 100     5327 if (defined($_[0]) && $_[0] =~ /^\s*([+-])\s*$/) {
912 163         484 $sign = $1;
913 163         282 shift;
914             }
915              
916             # Get the rounding parameters, if any.
917              
918 1588         2932 my @r = @_;
919              
920             # If called as a class method, initialize a new object.
921              
922 1588 100       4130 $self = bless {}, $class unless $selfref;
923              
924 1588         3989 $self -> {sign} = $sign;
925 1588         4869 $self -> {_m} = $LIB -> _one();
926 1588         3194 $self -> {_es} = '+';
927 1588         3944 $self -> {_e} = $LIB -> _zero();
928              
929             # If rounding parameters are given as arguments, use them. If no rounding
930             # parameters are given, and if called as a class method initialize the new
931             # instance with the class variables.
932              
933             #return $self -> round(@r); # this should work, but doesnt; fixme!
934              
935 1588 100       4061 if (@r) {
936 29 50 100     175 croak "can't specify both accuracy and precision"
      66        
937             if @r >= 2 && defined($r[0]) && defined($r[1]);
938 29         73 $self->{_a} = $_[0];
939 29         60 $self->{_p} = $_[1];
940             } else {
941 1559 100       3427 unless($selfref) {
942 1000         3057 $self->{_a} = $class -> accuracy();
943 1000         2996 $self->{_p} = $class -> precision();
944             }
945             }
946              
947 1588         6787 return $self;
948             }
949              
950             sub binf {
951             # create/assign a '+inf' or '-inf'
952              
953             # Class::method(...) -> Class->method(...)
954 2071 50 66 2071 1 23288 unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) ||
      66        
955             $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i))
956             {
957             #carp "Using ", (caller(0))[3], "() as a function is deprecated;",
958             # " use is as a method instead";
959 0         0 unshift @_, __PACKAGE__;
960             }
961              
962 2071         4147 my $self = shift;
963 2071         3583 my $selfref = ref $self;
964 2071   66     5781 my $class = $selfref || $self;
965              
966             {
967 43     43   455 no strict 'refs';
  43         121  
  43         20988  
  2071         3192  
968 2071 100       2865 if (${"${class}::_trap_inf"}) {
  2071         8723  
969 5         555 croak("Tried to create +-inf in $class->binf()");
970             }
971             }
972              
973 2066 100       4439 $self->import() if $IMPORT == 0; # make require work
974              
975             # Don't modify constant (read-only) objects.
976              
977 2066 50 66     5546 return $self if $selfref && $self->modify('binf');
978              
979 2066 100       4623 return $downgrade -> binf(@_) if $downgrade;
980              
981             # Get the sign.
982              
983 2059         3369 my $sign = '+'; # default is to return positive infinity
984 2059 100 100     10215 if (defined($_[0]) && $_[0] =~ /^\s*([+-])(inf|$)/i) {
985 1968         4279 $sign = $1;
986 1968         2869 shift;
987             }
988              
989             # Get the rounding parameters, if any.
990              
991 2059         4061 my @r = @_;
992              
993             # If called as a class method, initialize a new object.
994              
995 2059 100       5699 $self = bless {}, $class unless $selfref;
996              
997 2059         6574 $self -> {sign} = $sign . 'inf';
998 2059         6710 $self -> {_m} = $LIB -> _zero();
999 2059         4239 $self -> {_es} = '+';
1000 2059         4686 $self -> {_e} = $LIB -> _zero();
1001              
1002             # If rounding parameters are given as arguments, use them. If no rounding
1003             # parameters are given, and if called as a class method initialize the new
1004             # instance with the class variables.
1005              
1006             #return $self -> round(@r); # this should work, but doesnt; fixme!
1007              
1008 2059 100       5288 if (@r) {
1009 530 50 66     2646 croak "can't specify both accuracy and precision"
      33        
1010             if @r >= 2 && defined($r[0]) && defined($r[1]);
1011 530         1063 $self->{_a} = $r[0];
1012 530         1049 $self->{_p} = $r[1];
1013             } else {
1014 1529 100       3429 unless($selfref) {
1015 1234         3839 $self->{_a} = $class -> accuracy();
1016 1234         3412 $self->{_p} = $class -> precision();
1017             }
1018             }
1019              
1020 2059         22707 return $self;
1021             }
1022              
1023             sub bnan {
1024             # create/assign a 'NaN'
1025              
1026             # Class::method(...) -> Class->method(...)
1027 2115 50 66 2115 1 23423 unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) ||
      66        
1028             $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i))
1029             {
1030             #carp "Using ", (caller(0))[3], "() as a function is deprecated;",
1031             # " use is as a method instead";
1032 0         0 unshift @_, __PACKAGE__;
1033             }
1034              
1035 2115         4415 my $self = shift;
1036 2115         4086 my $selfref = ref $self;
1037 2115   66     5247 my $class = $selfref || $self;
1038              
1039             {
1040 43     43   422 no strict 'refs';
  43         126  
  43         406369  
  2115         3148  
1041 2115 100       2962 if (${"${class}::_trap_nan"}) {
  2115         8669  
1042 3         342 croak("Tried to create NaN in $class->bnan()");
1043             }
1044             }
1045              
1046 2112 100       4480 $self->import() if $IMPORT == 0; # make require work
1047              
1048             # Don't modify constant (read-only) objects.
1049              
1050 2112 50 66     7353 return $self if $selfref && $self->modify('bnan');
1051              
1052 2112 100       4254 return $downgrade -> bnan(@_) if defined $downgrade;
1053              
1054             # Get the rounding parameters, if any.
1055              
1056 2098         3961 my @r = @_;
1057              
1058             # If called as a class method, initialize a new object.
1059              
1060 2098 100       4764 $self = bless {}, $class unless $selfref;
1061              
1062 2098         4872 $self -> {sign} = $nan;
1063 2098         6287 $self -> {_m} = $LIB -> _zero();
1064 2098         4068 $self -> {_es} = '+';
1065 2098         4747 $self -> {_e} = $LIB -> _zero();
1066              
1067             # If rounding parameters are given as arguments, use them. If no rounding
1068             # parameters are given, and if called as a class method initialize the new
1069             # instance with the class variables.
1070              
1071             #return $self -> round(@r); # this should work, but doesnt; fixme!
1072              
1073 2098 100       5037 if (@r) {
1074 296 50 66     1367 croak "can't specify both accuracy and precision"
      33        
1075             if @r >= 2 && defined($r[0]) && defined($r[1]);
1076 296         606 $self->{_a} = $r[0];
1077 296         565 $self->{_p} = $r[1];
1078             } else {
1079 1802 100       3921 unless($selfref) {
1080 751         2231 $self->{_a} = $class -> accuracy();
1081 751         2223 $self->{_p} = $class -> precision();
1082             }
1083             }
1084              
1085 2098         20572 return $self;
1086             }
1087              
1088             sub bpi {
1089              
1090             # Class::method(...) -> Class->method(...)
1091 205 100 66 205 1 3944 unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) ||
      66        
1092             $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i))
1093             {
1094             #carp "Using ", (caller(0))[3], "() as a function is deprecated;",
1095             # " use is as a method instead";
1096 1         5 unshift @_, __PACKAGE__;
1097             }
1098              
1099             # Called as Argument list
1100             # --------- -------------
1101             # Math::BigFloat->bpi() ("Math::BigFloat")
1102             # Math::BigFloat->bpi(10) ("Math::BigFloat", 10)
1103             # $x->bpi() ($x)
1104             # $x->bpi(10) ($x, 10)
1105             # Math::BigFloat::bpi() ()
1106             # Math::BigFloat::bpi(10) (10)
1107             #
1108             # In ambiguous cases, we favour the OO-style, so the following case
1109             #
1110             # $n = Math::BigFloat->new("10");
1111             # $x = Math::BigFloat->bpi($n);
1112             #
1113             # which gives an argument list with the single element $n, is resolved as
1114             #
1115             # $n->bpi();
1116              
1117 205         479 my $self = shift;
1118 205         361 my $selfref = ref $self;
1119 205   66     597 my $class = $selfref || $self;
1120 205         444 my @r = @_; # rounding paramters
1121              
1122 205 100       399 if ($selfref) { # bpi() called as an instance method
1123 83 50       291 return $self if $self -> modify('bpi');
1124             } else { # bpi() called as a class method
1125 122         283 $self = bless {}, $class; # initialize new instance
1126             }
1127              
1128 205         634 ($self, @r) = $self -> _find_round_parameters(@r);
1129              
1130             # The accuracy, i.e., the number of digits. Pi has one digit before the
1131             # dot, so a precision of 4 digits is equivalent to an accuracy of 5 digits.
1132              
1133 205 50       804 my $n = defined $r[0] ? $r[0]
    100          
1134             : defined $r[1] ? 1 - $r[1]
1135             : $self -> div_scale();
1136              
1137 205 100       483 my $rmode = defined $r[2] ? $r[2] : $self -> round_mode();
1138              
1139 205         316 my $pi;
1140              
1141 205 50       430 if ($n <= 1000) {
1142              
1143             # 75 x 14 = 1050 digits
1144              
1145 205         383 my $all_digits = <
1146             314159265358979323846264338327950288419716939937510582097494459230781640628
1147             620899862803482534211706798214808651328230664709384460955058223172535940812
1148             848111745028410270193852110555964462294895493038196442881097566593344612847
1149             564823378678316527120190914564856692346034861045432664821339360726024914127
1150             372458700660631558817488152092096282925409171536436789259036001133053054882
1151             046652138414695194151160943305727036575959195309218611738193261179310511854
1152             807446237996274956735188575272489122793818301194912983367336244065664308602
1153             139494639522473719070217986094370277053921717629317675238467481846766940513
1154             200056812714526356082778577134275778960917363717872146844090122495343014654
1155             958537105079227968925892354201995611212902196086403441815981362977477130996
1156             051870721134999999837297804995105973173281609631859502445945534690830264252
1157             230825334468503526193118817101000313783875288658753320838142061717766914730
1158             359825349042875546873115956286388235378759375195778185778053217122680661300
1159             192787661119590921642019893809525720106548586327886593615338182796823030195
1160             EOF
1161              
1162             # Should we round up?
1163              
1164 205         313 my $round_up;
1165              
1166             # From the string above, we need to extract the number of digits we
1167             # want plus extra characters for the newlines.
1168              
1169 205         578 my $nchrs = $n + int($n / 75);
1170              
1171             # Extract the digits we want.
1172              
1173 205         507 my $digits = substr($all_digits, 0, $nchrs);
1174              
1175             # Find out whether we should round up or down. Rounding is easy, since
1176             # pi is trancendental. With directed rounding, it doesn't matter what
1177             # the following digits are. With rounding to nearest, we only have to
1178             # look at one extra digit.
1179              
1180 205 50       452 if ($rmode eq 'trunc') {
1181 0         0 $round_up = 0;
1182             } else {
1183 205         373 my $next_digit = substr($all_digits, $nchrs, 1);
1184 205 100       509 $round_up = $next_digit lt '5' ? 0 : 1;
1185             }
1186              
1187             # Remove the newlines.
1188              
1189 205         454 $digits =~ tr/0-9//cd;
1190              
1191             # Now do the rounding. We could easily make the regex substitution
1192             # handle all cases, but we avoid using the regex engine when it is
1193             # simple to avoid it.
1194              
1195 205 100       474 if ($round_up) {
1196 119         241 my $last_digit = substr($digits, -1, 1);
1197 119 100       293 if ($last_digit lt '9') {
1198 108         297 substr($digits, -1, 1) = ++$last_digit;
1199             } else {
1200 11         158 $digits =~ s{([0-8])(9+)$}
  11         90  
1201             { ($1 + 1) . ("0" x CORE::length($2)) }e;
1202             }
1203             }
1204              
1205             # Convert to an object.
1206 205         757  
1207             $pi = bless {
1208             sign => '+',
1209             _m => $LIB -> _new($digits),
1210             _es => '-',
1211             _e => $LIB -> _new($n - 1),
1212             }, $class;
1213              
1214             } else {
1215              
1216             # For large accuracy, the arctan formulas become very inefficient with
1217             # Math::BigFloat, so use Brent-Salamin (aka AGM or Gauss-Legendre).
1218              
1219 0         0 # Use a few more digits in the intermediate computations.
1220             $n += 8;
1221 0 0       0  
1222 0         0 $HALF = $class -> new($HALF) unless ref($HALF);
1223             my ($an, $bn, $tn, $pn)
1224             = ($class -> bone, $HALF -> copy() -> bsqrt($n),
1225 0         0 $HALF -> copy() -> bmul($HALF), $class -> bone);
1226 0         0 while ($pn < $n) {
1227 0         0 my $prev_an = $an -> copy();
1228 0         0 $an = $an -> badd($bn) -> bmul($HALF, $n);
1229 0         0 $bn = $bn -> bmul($prev_an) -> bsqrt($n);
1230 0         0 $prev_an = $prev_an -> bsub($an);
1231 0         0 $tn = $tn -> bsub($pn * $prev_an * $prev_an);
1232             $pn = $pn -> badd($pn);
1233 0         0 }
1234 0         0 $an = $an -> badd($bn);
1235             $an = $an -> bmul($an, $n) -> bdiv(4 * $tn, $n);
1236 0         0  
1237 0         0 $an = $an -> round(@r);
1238             $pi = $an;
1239             }
1240 205 100       644  
    50          
1241 191         690 if (defined $r[0]) {
1242             $pi -> accuracy($r[0]);
1243 0         0 } elsif (defined $r[1]) {
1244             $pi -> precision($r[1]);
1245             }
1246 205         461  
1247 1230         2464 for my $key (qw/ sign _m _es _e _a _p /) {
1248             $self -> {$key} = $pi -> {$key};
1249             }
1250 205 50 33     539  
1251             return $downgrade -> new($self -> bdstr(), @r)
1252 205         1261 if defined($downgrade) && $self->is_int();
1253             return $self;
1254             }
1255              
1256 17365     17365 1 193558 sub copy {
1257 17365 50       34358 my ($x, $class);
1258 17365         25698 if (ref($_[0])) { # $y = $x -> copy()
1259 17365         27731 $x = shift;
1260             $class = ref($x);
1261 0         0 } else { # $y = Math::BigInt -> copy($y)
1262 0         0 $class = shift;
1263             $x = shift;
1264             }
1265 17365 50       34133  
1266             carp "Rounding is not supported for ", (caller(0))[3], "()" if @_;
1267 17365         35051  
1268             my $copy = bless {}, $class;
1269 17365         42374  
1270 17365         33154 $copy->{sign} = $x->{sign};
1271 17365         46920 $copy->{_es} = $x->{_es};
1272 17365         40348 $copy->{_m} = $LIB->_copy($x->{_m});
1273 17365 100       44401 $copy->{_e} = $LIB->_copy($x->{_e});
1274 17365 100       38075 $copy->{_a} = $x->{_a} if exists $x->{_a};
1275             $copy->{_p} = $x->{_p} if exists $x->{_p};
1276 17365         45636  
1277             return $copy;
1278             }
1279              
1280             sub as_int {
1281 650 50   650 1 3872 # return copy as a bigint representation of this Math::BigFloat number
1282 650 50       1689 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
1283             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
1284 650 50       1620  
1285             return $x -> copy() if $x -> isa("Math::BigInt");
1286              
1287             # disable upgrading and downgrading
1288 650         4183  
1289 650         2312 require Math::BigInt;
1290 650         1967 my $upg = Math::BigInt -> upgrade();
1291 650         1886 my $dng = Math::BigInt -> downgrade();
1292 650         1786 Math::BigInt -> upgrade(undef);
1293             Math::BigInt -> downgrade(undef);
1294 650         973  
1295 650 100       1671 my $y;
    100          
1296 8         39 if ($x -> is_inf()) {
1297             $y = Math::BigInt -> binf($x->sign());
1298 4         55 } elsif ($x -> is_nan()) {
1299             $y = Math::BigInt -> bnan();
1300 638         2231 } else {
1301 638 100       2550 $y = $LIB->_copy($x->{_m});
    100          
1302 183         715 if ($x->{_es} eq '-') { # < 0
1303             $y = $LIB->_rsft($y, $x->{_e}, 10);
1304 14         94 } elsif (! $LIB->_is_zero($x->{_e})) { # > 0
1305             $y = $LIB->_lsft($y, $x->{_e}, 10);
1306 638         2252 }
1307             $y = Math::BigInt->new($x->{sign} . $LIB->_str($y));
1308             }
1309              
1310             # reset upgrading and downgrading
1311 650         3142  
1312 650         2181 Math::BigInt -> upgrade($upg);
1313             Math::BigInt -> downgrade($dng);
1314 650         3154  
1315             return $y;
1316             }
1317              
1318 2 50   2 1 23 sub as_float {
1319 2 50       7 my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
1320             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
1321 2 50       9  
1322             return $x -> copy() if $x -> isa("Math::BigFloat");
1323              
1324             # disable upgrading and downgrading
1325 0         0  
1326 0         0 require Math::BigFloat;
1327 0         0 my $upg = Math::BigFloat -> upgrade();
1328 0         0 my $dng = Math::BigFloat -> downgrade();
1329 0         0 Math::BigFloat -> upgrade(undef);
1330             Math::BigFloat -> downgrade(undef);
1331 0         0  
1332             my $y = Math::BigFloat -> new($x);
1333              
1334             # reset upgrading and downgrading
1335 0         0  
1336 0         0 Math::BigFloat -> upgrade($upg);
1337             Math::BigFloat -> downgrade($dng);
1338 0         0  
1339             return $y;
1340             }
1341              
1342             ###############################################################################
1343             # Boolean methods
1344             ###############################################################################
1345              
1346             sub is_zero {
1347 110809 100   110809 1 247660 # return true if arg (BFLOAT or num_str) is zero
1348             my (undef, $x) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
1349 110809 100 100     365170  
1350             ($x->{sign} eq '+' && $LIB->_is_zero($x->{_m})) ? 1 : 0;
1351             }
1352              
1353             sub is_one {
1354 2958 100   2958 1 10960 # return true if arg (BFLOAT or num_str) is +1 or -1 if signis given
1355             my (undef, $x, $sign) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
1356 2958 100 100     10366  
1357             $sign = '+' if !defined $sign || $sign ne '-';
1358              
1359             ($x->{sign} eq $sign &&
1360 2958 100 100     11105 $LIB->_is_zero($x->{_e}) &&
1361             $LIB->_is_one($x->{_m})) ? 1 : 0;
1362             }
1363              
1364             sub is_odd {
1365 104 50   104 1 841 # return true if arg (BFLOAT or num_str) is odd or false if even
1366             my (undef, $x) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
1367              
1368             (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1369 104 100 100     702 ($LIB->_is_zero($x->{_e})) &&
1370             ($LIB->_is_odd($x->{_m}))) ? 1 : 0;
1371             }
1372              
1373             sub is_even {
1374 72 50   72 1 958 # return true if arg (BINT or num_str) is even or false if odd
1375             my (undef, $x) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
1376              
1377             (($x->{sign} =~ /^[+-]$/) && # NaN & +-inf aren't
1378 72 100 100     747 ($x->{_es} eq '+') && # 123.45 isn't
1379             ($LIB->_is_even($x->{_m}))) ? 1 : 0; # but 1200 is
1380             }
1381              
1382             sub is_int {
1383 2548 50   2548 1 6901 # return true if arg (BFLOAT or num_str) is an integer
1384             my (undef, $x) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
1385              
1386 2548 100 100     17077 (($x->{sign} =~ /^[+-]$/) && # NaN and +-inf aren't
1387             ($x->{_es} eq '+')) ? 1 : 0; # 1e-1 => no integer
1388             }
1389              
1390             ###############################################################################
1391             # Comparison methods
1392             ###############################################################################
1393              
1394             sub bcmp {
1395             # Compares 2 values. Returns one of undef, <0, =0, >0. (suitable for sort)
1396              
1397 2913 100 100 2913 1 16640 # set up parameters
1398             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
1399             ? (ref($_[0]), @_)
1400             : objectify(2, @_);
1401 2913 50       6136  
1402             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
1403              
1404             # Handle all 'nan' cases.
1405 2913 100 100     11331  
1406             return if ($x->{sign} eq $nan) || ($y->{sign} eq $nan);
1407              
1408             # Handle all '+inf' and '-inf' cases.
1409              
1410 2855 100 100     11975 return 0 if ($x->{sign} eq '+inf' && $y->{sign} eq '+inf' ||
      100        
      100        
1411 2847 100       6041 $x->{sign} eq '-inf' && $y->{sign} eq '-inf');
1412 2791 100       5759 return +1 if $x->{sign} eq '+inf'; # x = +inf and y < +inf
1413 2727 50       5262 return -1 if $x->{sign} eq '-inf'; # x = -inf and y > -inf
1414 2727 100       4966 return -1 if $y->{sign} eq '+inf'; # x < +inf and y = +inf
1415             return +1 if $y->{sign} eq '-inf'; # x > -inf and y = -inf
1416              
1417             # Handle all cases with opposite signs.
1418 2723 100 100     9219  
1419 2267 100 100     6334 return +1 if $x->{sign} eq '+' && $y->{sign} eq '-'; # also does 0 <=> -y
1420             return -1 if $x->{sign} eq '-' && $y->{sign} eq '+'; # also does -x <=> 0
1421              
1422             # Handle all remaining zero cases.
1423 1961         4402  
1424 1961         4066 my $xz = $x->is_zero();
1425 1961 100 100     5631 my $yz = $y->is_zero();
1426 1827 100 66     4739 return 0 if $xz && $yz; # 0 <=> 0
1427 1707 100 66     6174 return -1 if $xz && $y->{sign} eq '+'; # 0 <=> +y
1428             return +1 if $yz && $x->{sign} eq '+'; # +x <=> 0
1429              
1430             # Both arguments are now finite, non-zero numbers with the same sign.
1431 1255         1989  
1432             my $cmp;
1433              
1434             # The next step is to compare the exponents, but since each mantissa is an
1435             # integer of arbitrary value, the exponents must be normalized by the length
1436             # of the mantissas before we can compare them.
1437 1255         3472  
1438 1255         3207 my $mxl = $LIB->_len($x->{_m});
1439             my $myl = $LIB->_len($y->{_m});
1440              
1441             # If the mantissas have the same length, there is no point in normalizing
1442             # the exponents by the length of the mantissas, so treat that as a special
1443             # case.
1444 1255 100       2961  
1445             if ($mxl == $myl) {
1446              
1447             # First handle the two cases where the exponents have different signs.
1448 1134 50 66     5495  
    100 100        
1449 0         0 if ($x->{_es} eq '+' && $y->{_es} eq '-') {
1450             $cmp = +1;
1451 16         45 } elsif ($x->{_es} eq '-' && $y->{_es} eq '+') {
1452             $cmp = -1;
1453             }
1454              
1455             # Then handle the case where the exponents have the same sign.
1456              
1457 1118         3208 else {
1458 1118 100       2849 $cmp = $LIB->_acmp($x->{_e}, $y->{_e});
1459             $cmp = -$cmp if $x->{_es} eq '-';
1460             }
1461              
1462             # Adjust for the sign, which is the same for x and y, and bail out if
1463             # we're done.
1464 1134 100       2319  
1465 1134 100       2618 $cmp = -$cmp if $x->{sign} eq '-'; # 124 > 123, but -124 < -123
1466             return $cmp if $cmp;
1467              
1468             }
1469              
1470             # We must normalize each exponent by the length of the corresponding
1471             # mantissa. Life is a lot easier if we first make both exponents
1472             # non-negative. We do this by adding the same positive value to both
1473             # exponent. This is safe, because when comparing the exponents, only the
1474             # relative difference is important.
1475 1187         2038  
1476             my $ex;
1477             my $ey;
1478 1187 100       2399  
1479             if ($x->{_es} eq '+') {
1480              
1481             # If the exponent of x is >= 0 and the exponent of y is >= 0, there is
1482             # no need to do anything special.
1483 1051 100       1886  
1484 1042         2774 if ($y->{_es} eq '+') {
1485 1042         2449 $ex = $LIB->_copy($x->{_e});
1486             $ey = $LIB->_copy($y->{_e});
1487             }
1488              
1489             # If the exponent of x is >= 0 and the exponent of y is < 0, add the
1490             # absolute value of the exponent of y to both.
1491              
1492 9         48 else {
1493 9         51 $ex = $LIB->_copy($x->{_e});
1494 9         42 $ex = $LIB->_add($ex, $y->{_e}); # ex + |ey|
1495             $ey = $LIB->_zero(); # -ex + |ey| = 0
1496             }
1497              
1498             } else {
1499              
1500             # If the exponent of x is < 0 and the exponent of y is >= 0, add the
1501             # absolute value of the exponent of x to both.
1502 136 100       408  
1503 25         117 if ($y->{_es} eq '+') {
1504 25         131 $ex = $LIB->_zero(); # -ex + |ex| = 0
1505 25         111 $ey = $LIB->_copy($y->{_e});
1506             $ey = $LIB->_add($ey, $x->{_e}); # ey + |ex|
1507             }
1508              
1509             # If the exponent of x is < 0 and the exponent of y is < 0, add the
1510             # absolute values of both exponents to both exponents.
1511              
1512 111         378 else {
1513 111         380 $ex = $LIB->_copy($y->{_e}); # -ex + |ey| + |ex| = |ey|
1514             $ey = $LIB->_copy($x->{_e}); # -ey + |ex| + |ey| = |ex|
1515             }
1516              
1517             }
1518              
1519             # Now we can normalize the exponents by adding lengths of the mantissas.
1520 1187         3365  
1521 1187         3964 $ex = $LIB->_add($ex, $LIB->_new($mxl));
1522             $ey = $LIB->_add($ey, $LIB->_new($myl));
1523              
1524             # We're done if the exponents are different.
1525 1187         3826  
1526 1187 100       2973 $cmp = $LIB->_acmp($ex, $ey);
1527 1187 100       2816 $cmp = -$cmp if $x->{sign} eq '-'; # 124 > 123, but -124 < -123
1528             return $cmp if $cmp;
1529              
1530             # Compare the mantissas, but first normalize them by padding the shorter
1531             # mantissa with zeros (shift left) until it has the same length as the
1532             # longer mantissa.
1533 1097         1731  
1534 1097         1679 my $mx = $x->{_m};
1535             my $my = $y->{_m};
1536 1097 100       2754  
    100          
1537 26         110 if ($mxl > $myl) {
1538             $my = $LIB->_lsft($LIB->_copy($my), $LIB->_new($mxl - $myl), 10);
1539 5         37 } elsif ($mxl < $myl) {
1540             $mx = $LIB->_lsft($LIB->_copy($mx), $LIB->_new($myl - $mxl), 10);
1541             }
1542 1097         2402  
1543 1097 100       2789 $cmp = $LIB->_acmp($mx, $my);
1544 1097         3951 $cmp = -$cmp if $x->{sign} eq '-'; # 124 > 123, but -124 < -123
1545             return $cmp;
1546              
1547             }
1548              
1549             sub bacmp {
1550             # Compares 2 values, ignoring their signs.
1551             # Returns one of undef, <0, =0, >0. (suitable for sort)
1552              
1553 8979 100 66 8979 1 46994 # set up parameters
1554             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
1555             ? (ref($_[0]), @_)
1556             : objectify(2, @_);
1557 8979 50       19943  
1558             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
1559              
1560 8979 100 100     49027 # handle +-inf and NaN's
1561 92 100 100     634 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/) {
1562 64 100 100     195 return if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1563 48 100 66     115 return 0 if ($x->is_inf() && $y->is_inf());
1564 16         160 return 1 if ($x->is_inf() && !$y->is_inf());
1565             return -1;
1566             }
1567              
1568 8887         20486 # shortcut
1569 8887         17734 my $xz = $x->is_zero();
1570 8887 100 100     20693 my $yz = $y->is_zero();
1571 8883 100 66     19380 return 0 if $xz && $yz; # 0 <=> 0
1572 8843 100 66     18046 return -1 if $xz && !$yz; # 0 <=> +y
1573             return 1 if $yz && !$xz; # +x <=> 0
1574              
1575 8803         22090 # adjust so that exponents are equal
1576 8803         20429 my $lxm = $LIB->_len($x->{_m});
1577 8803         16217 my $lym = $LIB->_len($y->{_m});
1578 8803 100       19919 my ($xes, $yes) = (1, 1);
1579 8803 100       17767 $xes = -1 if $x->{_es} ne '+';
1580             $yes = -1 if $y->{_es} ne '+';
1581 8803         21970 # the numify somewhat limits our length, but makes it much faster
1582 8803         20499 my $lx = $lxm + $xes * $LIB->_num($x->{_e});
1583 8803         15049 my $ly = $lym + $yes * $LIB->_num($y->{_e});
1584 8803 100       26216 my $l = $lx - $ly;
1585             return $l <=> 0 if $l != 0;
1586              
1587             # lengths (corrected by exponent) are equal
1588 3821         5429 # so make mantissa equal-length by padding with zero (shift left)
1589 3821         6090 my $diff = $lxm - $lym;
1590 3821         5923 my $xm = $x->{_m}; # not yet copy it
1591 3821 100       9700 my $ym = $y->{_m};
    100          
1592 578         1996 if ($diff > 0) {
1593 578         1964 $ym = $LIB->_copy($y->{_m});
1594             $ym = $LIB->_lsft($ym, $LIB->_new($diff), 10);
1595 258         1153 } elsif ($diff < 0) {
1596 258         1030 $xm = $LIB->_copy($x->{_m});
1597             $xm = $LIB->_lsft($xm, $LIB->_new(-$diff), 10);
1598 3821         10505 }
1599             $LIB->_acmp($xm, $ym);
1600             }
1601              
1602             ###############################################################################
1603             # Arithmetic methods
1604             ###############################################################################
1605              
1606             sub bneg {
1607             # (BINT or num_str) return BINT
1608 475 50   475 1 2006 # negate number or make a negated number from string
1609             my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
1610 475 50       1469  
1611             return $x if $x->modify('bneg');
1612 475 100       1224  
1613             return $x -> bnan(@r) if $x -> is_nan();
1614              
1615             # For +0 do not negate (to have always normalized +0).
1616 466 100 100     2052 $x->{sign} =~ tr/+-/-+/
1617             unless $x->{sign} eq '+' && $LIB->_is_zero($x->{_m});
1618 466 100 66     1236  
      66        
1619             return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade)
1620 463         1322 && ($x -> is_int() || $x -> is_inf() || $x -> is_nan());
1621             return $x -> round(@r);
1622             }
1623              
1624             sub bnorm {
1625             # bnorm() can't support rounding, because bround() and bfround() call
1626             # bnorm(), which would recurse indefinitely.
1627              
1628 72362 50   72362 1 194025 # adjust m and e so that m is smallest possible
1629             my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
1630 72362 50       145826  
1631             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
1632              
1633 72362 100       260832 # inf, nan etc
1634 6 100       44 if ($x->{sign} !~ /^[+-]$/) {
1635 4         24 return $downgrade -> new($x) if defined $downgrade;
1636             return $x;
1637             }
1638 72356         193965  
1639 72356 100       138550 my $zeros = $LIB->_zeros($x->{_m}); # correct for trailing zeros
1640 29435         72550 if ($zeros != 0) {
1641 29435         91500 my $z = $LIB->_new($zeros);
1642 29435 100       64996 $x->{_m} = $LIB->_rsft($x->{_m}, $z, 10);
1643 27727 100       75583 if ($x->{_es} eq '-') {
1644 27387         75363 if ($LIB->_acmp($x->{_e}, $z) >= 0) {
1645 27387 100       69323 $x->{_e} = $LIB->_sub($x->{_e}, $z);
1646             $x->{_es} = '+' if $LIB->_is_zero($x->{_e});
1647 340         1097 } else {
1648 340         1031 $x->{_e} = $LIB->_sub($LIB->_copy($z), $x->{_e});
1649             $x->{_es} = '+';
1650             }
1651 1708         5923 } else {
1652             $x->{_e} = $LIB->_add($x->{_e}, $z);
1653             }
1654             } else {
1655             # $x can only be 0Ey if there are no trailing zeros ('0' has 0 trailing
1656 42921 100       110738 # zeros). So, for something like 0Ey, set y to 0, and -0 => +0
1657 399         997 if ($LIB->_is_zero($x->{_m})) {
1658 399         718 $x->{sign} = '+';
1659 399         941 $x->{_es} = '+';
1660             $x->{_e} = $LIB->_zero();
1661             }
1662             }
1663 72356 100 100     172138  
1664             return $downgrade -> new($x)
1665 72336         238153 if defined($downgrade) && $x->is_int();
1666             return $x;
1667             }
1668              
1669             sub binc {
1670 4362 50   4362 1 12758 # increment arg by one
1671             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
1672 4362 50       12630  
1673             return $x if $x->modify('binc');
1674              
1675             # Inf and NaN
1676 4362 100       11414  
1677 4357 100       10603 return $x -> bnan(@r) if $x -> is_nan();
1678             return $x -> binf($x->{sign}, @r) if $x -> is_inf();
1679              
1680             # Non-integer
1681 4348 100       10526  
1682 461         1549 if ($x->{_es} eq '-') {
1683             return $x->badd($class->bone(), @r);
1684             }
1685              
1686             # If the exponent is non-zero, convert the internal representation, so that,
1687             # e.g., 12e+3 becomes 12000e+0 and we can easily increment the mantissa.
1688 3887 100       10301  
1689 350         1633 if (!$LIB->_is_zero($x->{_e})) {
1690 350         1568 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # 1e2 => 100
1691 350         842 $x->{_e} = $LIB->_zero(); # normalize
1692             $x->{_es} = '+';
1693             # we know that the last digit of $x will be '1' or '9', depending on the
1694             # sign
1695             }
1696              
1697 3887 100       8359 # now $x->{_e} == 0
    50          
1698 3871         9915 if ($x->{sign} eq '+') {
1699 3871         8651 $x->{_m} = $LIB->_inc($x->{_m});
1700             return $x->bnorm()->bround(@r);
1701 16         51 } elsif ($x->{sign} eq '-') {
1702 16 100       58 $x->{_m} = $LIB->_dec($x->{_m});
1703 16         50 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # -1 +1 => -0 => +0
1704             return $x->bnorm()->bround(@r);
1705             }
1706 0 0 0     0  
1707             return $downgrade -> new($x -> bdstr(), @r)
1708 0         0 if defined($downgrade) && $x -> is_int();
1709             return $x;
1710             }
1711              
1712             sub bdec {
1713 168 50   168 1 1289 # decrement arg by one
1714             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
1715 168 50       688  
1716             return $x if $x->modify('bdec');
1717              
1718             # Inf and NaN
1719 168 100       585  
1720 163 100       563 return $x -> bnan(@r) if $x -> is_nan();
1721             return $x -> binf($x->{sign}, @r) if $x -> is_inf();
1722              
1723             # Non-integer
1724 154 100       728  
1725 113         430 if ($x->{_es} eq '-') {
1726             return $x->badd($class->bone('-'), @r);
1727             }
1728              
1729             # If the exponent is non-zero, convert the internal representation, so that,
1730             # e.g., 12e+3 becomes 12000e+0 and we can easily increment the mantissa.
1731 41 100       146  
1732 8         42 if (!$LIB->_is_zero($x->{_e})) {
1733 8         43 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # 1e2 => 100
1734 8         26 $x->{_e} = $LIB->_zero(); # normalize
1735             $x->{_es} = '+';
1736             }
1737              
1738 41         125 # now $x->{_e} == 0
1739 41 100 100     221 my $zero = $x->is_zero();
    50          
1740 21         78 if (($x->{sign} eq '-') || $zero) { # x <= 0
1741 21 100       61 $x->{_m} = $LIB->_inc($x->{_m});
1742 21 50       58 $x->{sign} = '-' if $zero; # 0 => 1 => -1
1743 21         70 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # -1 +1 => -0 => +0
1744             return $x->bnorm()->round(@r);
1745             }
1746 20         95 elsif ($x->{sign} eq '+') { # x > 0
1747 20         57 $x->{_m} = $LIB->_dec($x->{_m});
1748             return $x->bnorm()->round(@r);
1749             }
1750 0 0 0     0  
1751             return $downgrade -> new($x -> bdstr(), @r)
1752 0         0 if defined($downgrade) && $x -> is_int();
1753             return $x -> round(@r);
1754             }
1755              
1756             sub badd {
1757 13458 100 100 13458 1 62496 # set up parameters
1758             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
1759             ? (ref($_[0]), @_)
1760             : objectify(2, @_);
1761 13458 50       39180  
1762             return $x if $x->modify('badd');
1763              
1764 13458 100 100     72060 # inf and NaN handling
1765             if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/) {
1766              
1767 226 100 100     1974 # $x is NaN and/or $y is NaN
    100 100        
    100          
1768 110         244 if ($x->{sign} eq $nan || $y->{sign} eq $nan) {
1769             $x = $x->bnan();
1770             }
1771              
1772             # $x is Inf and $y is Inf
1773             elsif ($x->{sign} =~ /^[+-]inf$/ && $y->{sign} =~ /^[+-]inf$/) {
1774 56 100       231 # +Inf + +Inf or -Inf + -Inf => same, rest is NaN
1775             $x = $x->bnan() if $x->{sign} ne $y->{sign};
1776             }
1777              
1778             # +-inf + something => +-inf; something +-inf => +-inf
1779 38         84 elsif ($y->{sign} =~ /^[+-]inf$/) {
1780             $x->{sign} = $y->{sign};
1781             }
1782 226 100       574  
1783 218         693 return $downgrade -> new($x -> bdstr(), @r) if defined $downgrade;
1784             return $x -> round(@r);
1785             }
1786 13232 100       27302  
1787             return $upgrade->badd($x, $y, @r) if defined $upgrade;
1788 13231         21989  
1789             $r[3] = $y; # no push!
1790              
1791 13231 100       26874 # for speed: no add for $x + 0
    100          
1792 24         121 if ($y->is_zero()) {
1793             $x = $x->round(@r);
1794             }
1795              
1796             # for speed: no add for 0 + $y
1797             elsif ($x->is_zero()) {
1798 102         503 # make copy, clobbering up x (modify in place!)
1799 102         338 $x->{_e} = $LIB->_copy($y->{_e});
1800 102         398 $x->{_es} = $y->{_es};
1801 102   33     450 $x->{_m} = $LIB->_copy($y->{_m});
1802 102         339 $x->{sign} = $y->{sign} || $nan;
1803             $x = $x->round(@r);
1804             }
1805              
1806             # both $x and $y are non-zero
1807             else {
1808              
1809 13105         23011 # take lower of the two e's and adapt m1 to it to match m2
1810 13105 50       26338 my $e = $y->{_e};
1811 13105         31448 $e = $LIB->_zero() if !defined $e; # if no BFLOAT?
1812             $e = $LIB->_copy($e); # make copy (didn't do it yet)
1813 13105         20147  
1814             my $es;
1815 13105   50     49892  
1816             ($e, $es) = $LIB -> _ssub($e, $y->{_es} || '+', $x->{_e}, $x->{_es});
1817 13105         36304  
1818             my $add = $LIB->_copy($y->{_m});
1819 13105 100       31909  
    100          
1820 7475         21096 if ($es eq '-') { # < 0
1821 7475         24278 $x->{_m} = $LIB->_lsft($x->{_m}, $e, 10);
1822             ($x->{_e}, $x->{_es}) = $LIB -> _sadd($x->{_e}, $x->{_es}, $e, $es);
1823 789         2592 } elsif (!$LIB->_is_zero($e)) { # > 0
1824             $add = $LIB->_lsft($add, $e, 10);
1825             }
1826              
1827             # else: both e are the same, so just leave them
1828 13105 100       31498  
1829 11511         29034 if ($x->{sign} eq $y->{sign}) {
1830             $x->{_m} = $LIB->_add($x->{_m}, $add);
1831             } else {
1832 1594         4981 ($x->{_m}, $x->{sign}) =
1833             $LIB -> _sadd($x->{_m}, $x->{sign}, $add, $y->{sign});
1834             }
1835              
1836 13105         32818 # delete trailing zeros, then round
1837             $x = $x->bnorm()->round(@r);
1838             }
1839 13231 100 66     33458  
1840             return $downgrade -> new($x -> bdstr(), @r)
1841 13224         44236 if defined($downgrade) && $x -> is_int();
1842             return $x; # rounding already done above
1843             }
1844              
1845             sub bsub {
1846 1655 100 100 1655 1 10298 # set up parameters
1847             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
1848             ? (ref($_[0]), @_)
1849             : objectify(2, @_);
1850 1655 50       5237  
1851             return $x if $x -> modify('bsub');
1852 1655 100       3407  
1853 42         189 if ($y -> is_zero()) {
1854             $x = $x -> round(@r);
1855             } else {
1856              
1857             # To correctly handle the special case $x -> bsub($x), we note the sign
1858             # of $x, then flip the sign of $y, and if the sign of $x changed too,
1859             # then we know that $x and $y are the same object.
1860 1613         3313  
1861 1613         3297 my $xsign = $x -> {sign};
1862 1613 100       3986 $y -> {sign} =~ tr/+-/-+/; # does nothing for NaN
1863             if ($xsign ne $x -> {sign}) {
1864 24 100       125 # special case of $x -> bsub($x) results in 0
1865 16         50 if ($xsign =~ /^[+-]$/) {
1866             $x = $x -> bzero(@r);
1867 8         34 } else {
1868             $x = $x -> bnan(); # NaN, -inf, +inf
1869 24 50       66 }
1870 24         73 return $downgrade -> new($x -> bdstr(), @r) if defined $downgrade;
1871             return $x -> round(@r);
1872 1589         4117 }
1873 1589         4023 $x = $x -> badd($y, @r); # badd does not leave internal zeros
1874             $y -> {sign} =~ tr/+-/-+/; # reset $y (does nothing for NaN)
1875             }
1876 1631 50 66     3502  
      66        
1877             return $downgrade -> new($x -> bdstr(), @r)
1878 1623         5678 if defined($downgrade) && ($x->is_int() || $x->is_inf() || $x->is_nan());
1879             $x; # already rounded by badd() or no rounding
1880             }
1881              
1882             sub bmul {
1883             # multiply two numbers
1884              
1885 19512 100 100 19512 1 93861 # set up parameters
1886             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
1887             ? (ref($_[0]), @_)
1888             : objectify(2, @_);
1889 19512 50       55592  
1890             return $x if $x->modify('bmul');
1891 19512 100 100     70665  
1892             return $x->bnan(@r) if ($x->{sign} eq $nan) || ($y->{sign} eq $nan);
1893              
1894 19455 100 100     64484 # inf handling
1895 85 100 100     297 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/)) {
1896             return $x->bnan(@r) if $x->is_zero() || $y->is_zero();
1897             # result will always be +-inf:
1898             # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1899 73 100 100     482 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1900 50 100 100     314 return $x->binf(@r) if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1901 36         124 return $x->binf(@r) if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1902             return $x->binf('-', @r);
1903             }
1904 19370 100       35840  
1905             return $upgrade->bmul($x, $y, @r) if defined $upgrade;
1906              
1907 19369         54852 # aEb * cEd = (a*c)E(b+d)
1908             $x->{_m} = $LIB->_mul($x->{_m}, $y->{_m});
1909 19369         65654 ($x->{_e}, $x->{_es})
1910             = $LIB -> _sadd($x->{_e}, $x->{_es}, $y->{_e}, $y->{_es});
1911 19369         36121  
1912             $r[3] = $y; # no push!
1913              
1914 19369 100       45161 # adjust sign:
1915 19369         41603 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1916             $x = $x->bnorm->round(@r);
1917 19369 0 33     39991  
      66        
1918             return $downgrade -> new($x -> bdstr(), @r)
1919 19365         48845 if defined($downgrade) && ($x->is_int() || $x->is_inf() || $x->is_nan());
1920             return $x;
1921             }
1922              
1923             sub bmuladd {
1924             # multiply two numbers and add the third to the result
1925              
1926 247 100 66 247 1 4019 # set up parameters
1927             my ($class, $x, $y, $z, @r)
1928             = ref($_[0]) && ref($_[0]) eq ref($_[1]) && ref($_[1]) eq ref($_[2])
1929             ? (ref($_[0]), @_)
1930             : objectify(3, @_);
1931 247 50       822  
1932             return $x if $x->modify('bmuladd');
1933              
1934             return $x->bnan(@r) if (($x->{sign} eq $nan) ||
1935 247 100 100     1335 ($y->{sign} eq $nan) ||
      100        
1936             ($z->{sign} eq $nan));
1937              
1938 214 100 66     832 # inf handling
1939 18 50 33     50 if (($x->{sign} =~ /^[+-]inf$/) || ($y->{sign} =~ /^[+-]inf$/)) {
1940             return $x->bnan(@r) if $x->is_zero() || $y->is_zero();
1941             # result will always be +-inf:
1942             # +inf * +/+inf => +inf, -inf * -/-inf => +inf
1943 18 100 100     243 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1944 12 100 100     98 return $x->binf(@r) if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1945 8         30 return $x->binf(@r) if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1946             return $x->binf('-', @r);
1947             }
1948              
1949 196         703 # aEb * cEd = (a*c)E(b+d)
1950             $x->{_m} = $LIB->_mul($x->{_m}, $y->{_m});
1951 196         759 ($x->{_e}, $x->{_es})
1952             = $LIB -> _sadd($x->{_e}, $x->{_es}, $y->{_e}, $y->{_es});
1953 196         373  
1954             $r[3] = $y; # no push!
1955              
1956 196 100       490 # adjust sign:
1957             $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1958              
1959 196 100       424 # z=inf handling (z=NaN handled above)
1960 1         4 if ($z->{sign} =~ /^[+-]inf$/) {
1961 1 50       9 $x->{sign} = $z->{sign};
1962 0         0 return $downgrade -> new($x -> bdstr(), @r) if defined $downgrade;
1963             return $x -> round(@r);
1964             }
1965              
1966 195         315 # take lower of the two e's and adapt m1 to it to match m2
1967 195 50       388 my $e = $z->{_e};
1968 195         515 $e = $LIB->_zero() if !defined $e; # if no BFLOAT?
1969             $e = $LIB->_copy($e); # make copy (didn't do it yet)
1970 195         288  
1971             my $es;
1972 195   50     664  
1973             ($e, $es) = $LIB -> _ssub($e, $z->{_es} || '+', $x->{_e}, $x->{_es});
1974 195         588  
1975             my $add = $LIB->_copy($z->{_m});
1976 195 100       603  
    100          
1977             if ($es eq '-') # < 0
1978 4         29 {
1979 4         20 $x->{_m} = $LIB->_lsft($x->{_m}, $e, 10);
1980             ($x->{_e}, $x->{_es}) = $LIB -> _sadd($x->{_e}, $x->{_es}, $e, $es);
1981             } elsif (!$LIB->_is_zero($e)) # > 0
1982 9         40 {
1983             $add = $LIB->_lsft($add, $e, 10);
1984             }
1985             # else: both e are the same, so just leave them
1986 195 100       553  
1987             if ($x->{sign} eq $z->{sign}) {
1988 151         396 # add
1989             $x->{_m} = $LIB->_add($x->{_m}, $add);
1990             } else {
1991 44         139 ($x->{_m}, $x->{sign}) =
1992             $LIB -> _sadd($x->{_m}, $x->{sign}, $add, $z->{sign});
1993             }
1994              
1995 195         560 # delete trailing zeros, then round
1996             $x = $x->bnorm()->round(@r);
1997 195 0 33     486  
      66        
1998             return $downgrade -> new($x -> bdstr(), @r)
1999 192         2391 if defined($downgrade) && ($x->is_int() || $x->is_inf() || $x->is_nan());
2000             return $x;
2001             }
2002              
2003             sub bdiv {
2004             # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return
2005             # (BFLOAT, BFLOAT) (quo, rem) or BFLOAT (only quo)
2006              
2007 9074     9074 1 29070 # set up parameters
2008             my ($class, $x, $y, @r) = (ref($_[0]), @_);
2009 9074 100 100     36016 # objectify is costly, so avoid it
2010 137         452 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) {
2011             ($class, $x, $y, @r) = objectify(2, @_);
2012             }
2013 9074 50       23776  
2014             return $x if $x->modify('bdiv');
2015 9074         15749  
2016             my $wantarray = wantarray; # call only once
2017              
2018             # At least one argument is NaN. This is handled the same way as in
2019             # Math::BigInt -> bdiv().
2020 9074 100 100     21917  
2021 68 100       265 if ($x -> is_nan() || $y -> is_nan()) {
2022             return $wantarray ? ($x -> bnan(@r), $class -> bnan(@r))
2023             : $x -> bnan(@r);
2024             }
2025              
2026             # Divide by zero and modulo zero. This is handled the same way as in
2027             # Math::BigInt -> bdiv(). See the comment in the code for Math::BigInt ->
2028             # bdiv() for further details.
2029 9006 100       19783  
2030 55         122 if ($y -> is_zero()) {
2031 55 100       143 my ($quo, $rem);
2032 16         73 if ($wantarray) {
2033 16 50 33     77 $rem = $x -> copy() -> round(@r);
2034             $rem = $downgrade -> new($rem, @r)
2035             if defined($downgrade) && $rem -> is_int();
2036 55 100       132 }
2037 21         85 if ($x -> is_zero()) {
2038             $quo = $x -> bnan(@r);
2039 34         112 } else {
2040             $quo = $x -> binf($x -> {sign}, @r);
2041 52 100       506 }
2042             return $wantarray ? ($quo, $rem) : $quo;
2043             }
2044              
2045             # Numerator (dividend) is +/-inf. This is handled the same way as in
2046             # Math::BigInt -> bdiv(). See the comment in the code for Math::BigInt ->
2047             # bdiv() for further details.
2048 8951 100       22050  
2049 40         156 if ($x -> is_inf()) {
2050 40 100       161 my ($quo, $rem);
2051 40 100       130 $rem = $class -> bnan(@r) if $wantarray;
2052 16         48 if ($y -> is_inf()) {
2053             $quo = $x -> bnan(@r);
2054 24 100       148 } else {
2055 24         96 my $sign = $x -> bcmp(0) == $y -> bcmp(0) ? '+' : '-';
2056             $quo = $x -> binf($sign, @r);
2057 40 100       226 }
2058             return $wantarray ? ($quo, $rem) : $quo;
2059             }
2060              
2061             # Denominator (divisor) is +/-inf. This is handled the same way as in
2062             # Math::BigInt -> bdiv(), with one exception: In scalar context,
2063             # Math::BigFloat does true division (although rounded), not floored division
2064             # (F-division), so a finite number divided by +/-inf is always zero. See the
2065             # comment in the code for Math::BigInt -> bdiv() for further details.
2066 8911 100       19601  
2067 40         110 if ($y -> is_inf()) {
2068 40 100       104 my ($quo, $rem);
2069 16 100 100     39 if ($wantarray) {
2070 12         42 if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) {
2071 12 50 33     53 $rem = $x -> copy() -> round(@r);
2072             $rem = $downgrade -> new($rem, @r)
2073 12         40 if defined($downgrade) && $rem -> is_int();
2074             $quo = $x -> bzero(@r);
2075 4         32 } else {
2076 4         39 $rem = $class -> binf($y -> {sign}, @r);
2077             $quo = $x -> bone('-', @r);
2078 16         77 }
2079             return ($quo, $rem);
2080 24 50       71 } else {
2081 24 50 33     77 if ($y -> is_inf()) {
2082 0         0 if ($x -> is_nan() || $x -> is_inf()) {
2083             return $x -> bnan(@r);
2084 24         81 } else {
2085             return $x -> bzero(@r);
2086             }
2087             }
2088             }
2089             }
2090              
2091             # At this point, both the numerator and denominator are finite numbers, and
2092             # the denominator (divisor) is non-zero.
2093              
2094 8871 100       16552 # x == 0?
2095 58         166 if ($x->is_zero()) {
2096 58         238 my ($quo, $rem);
2097 58 100 66     347 $quo = $x->round(@r);
2098             $quo = $downgrade -> new($quo, @r)
2099 58 100       172 if defined($downgrade) && $quo -> is_int();
2100 13         60 if ($wantarray) {
2101 13         78 $rem = $class -> bzero(@r);
2102             return $quo, $rem;
2103 45         286 }
2104             return $quo;
2105             }
2106              
2107             # Division might return a value that we can not represent exactly, so
2108             # upgrade, if upgrading is enabled.
2109              
2110 8813 0 33     20680 return $upgrade -> bdiv($x, $y, @r)
      33        
2111             if defined($upgrade) && !wantarray && !$LIB -> _is_one($y -> {_m});
2112              
2113 8813         12522 # we need to limit the accuracy to protect against overflow
2114 8813         13884 my $fallback = 0;
2115 8813         33505 my (@params, $scale);
2116             ($x, @params) = $x->_find_round_parameters($r[0], $r[1], $r[2], $y);
2117 8813 50       28432  
2118             return $x -> round(@r) if $x->is_nan(); # error in _find_round_parameters?
2119              
2120 8813 100       18863 # no rounding at all, so must use fallback
2121             if (scalar @params == 0) {
2122 506         1640 # simulate old behaviour
2123 506         857 $params[0] = $class->div_scale(); # and round to it as accuracy
2124 506         831 $scale = $params[0]+4; # at least four more for proper round
2125 506         736 $params[2] = $r[2]; # round mode by caller or undef
2126             $fallback = 1; # to clear a/p afterwards
2127             } else {
2128             # the 4 below is empirical, and there might be cases where it is not
2129 8307   66     20053 # enough...
2130             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2131             }
2132 8813         12702  
2133 8813 100       16452 my $rem;
2134             $rem = $class -> bzero() if wantarray;
2135 8813 50       20323  
2136             $y = $class->new($y) unless $y->isa('Math::BigFloat');
2137 8813         28996  
2138 8813         20741 my $lx = $LIB -> _len($x->{_m});
2139 8813 100       19497 my $ly = $LIB -> _len($y->{_m});
2140 8813 100       17562 $scale = $lx if $lx > $scale;
2141 8813         12750 $scale = $ly if $ly > $scale;
2142 8813 100       16844 my $diff = $ly - $lx;
2143             $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
2144              
2145 8813   100     21106 # check that $y is not 1 nor -1 and cache the result:
2146             my $y_not_one = !($LIB->_is_zero($y->{_e}) && $LIB->_is_one($y->{_m}));
2147              
2148             # flipping the sign of $y will also flip the sign of $x for the special
2149 8813         16959 # case of $x->bsub($x); so we can catch it below:
2150 8813         18905 my $xsign = $x->{sign};
2151             $y->{sign} =~ tr/+-/-+/;
2152 8813 100       18352  
2153             if ($xsign ne $x->{sign}) {
2154 8         32 # special case of $x /= $x results in 1
2155             $x = $x->bone(); # "fixes" also sign of $y, since $x is $y
2156             } else {
2157 8805         13615 # correct $y's sign again
2158             $y->{sign} =~ tr/+-/-+/;
2159             # continue with normal div code:
2160              
2161             # make copy of $x in case of list context for later remainder
2162 8805 100 100     20680 # calculation
2163 25         96 if (wantarray && $y_not_one) {
2164             $rem = $x->copy();
2165             }
2166 8805 100       23297  
2167             $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
2168              
2169 8805 100       20015 # check for / +-1 (+/- 1E0)
2170             if ($y_not_one) {
2171             # promote Math::BigInt and its subclasses (except when already a
2172 8636 50       16969 # Math::BigFloat)
2173             $y = $class->new($y) unless $y->isa('Math::BigFloat');
2174              
2175             # calculate the result to $scale digits and then round it
2176 8636         27541 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
2177 8636         30758 $x->{_m} = $LIB->_lsft($x->{_m}, $LIB->_new($scale), 10);
2178             $x->{_m} = $LIB->_div($x->{_m}, $y->{_m}); # a/c
2179              
2180             # correct exponent of $x
2181 8636         32701 ($x->{_e}, $x->{_es})
2182             = $LIB -> _ssub($x->{_e}, $x->{_es}, $y->{_e}, $y->{_es});
2183             # correct for 10**scale
2184 8636         27817 ($x->{_e}, $x->{_es})
2185 8636         26987 = $LIB -> _ssub($x->{_e}, $x->{_es}, $LIB->_new($scale), '+');
2186             $x = $x->bnorm(); # remove trailing 0's
2187             }
2188             } # end else $x != $y
2189              
2190 8813 100       18046 # shortcut to not run through _find_round_parameters again
2191 8780         14636 if (defined $params[0]) {
2192 8780         22147 $x->{_a} = undef; # clear before round
2193             $x = $x->bround($params[0], $params[2]); # then round accordingly
2194 33         89 } else {
2195 33         99 $x->{_p} = undef; # clear before round
2196             $x = $x->bfround($params[1], $params[2]); # then round accordingly
2197 8813 100       19682 }
2198             if ($fallback) {
2199 506         933 # clear a/p after round, since user did not request it
2200 506         905 $x->{_a} = undef;
2201             $x->{_p} = undef;
2202             }
2203 8813 100       16671  
2204 49 100       126 if (wantarray) {
2205 25         97 if ($y_not_one) {
2206 25         102 $x = $x -> bfloor();
2207             $rem = $rem->bmod($y, @params); # copy already done
2208 49 50       128 }
2209             if ($fallback) {
2210 49         93 # clear a/p after round, since user did not request it
2211 49         81 $rem->{_a} = undef;
2212             $rem->{_p} = undef;
2213 49 50 33     182 }
2214             $x = $downgrade -> new($x -> bdstr(), @r)
2215 49 50 33     134 if defined($downgrade) && $x -> is_int();
2216             $rem = $downgrade -> new($rem -> bdstr(), @r)
2217 49         264 if defined($downgrade) && $rem -> is_int();
2218             return ($x, $rem);
2219             }
2220 8764 50 66     18325  
2221             $x = $downgrade -> new($x, @r)
2222 8764         29787 if defined($downgrade) && $x -> is_int();
2223             $x; # rounding already done above
2224             }
2225              
2226             sub bmod {
2227             # (dividend: BFLOAT or num_str, divisor: BFLOAT or num_str) return remainder
2228              
2229 743 100 100 743 1 8087 # set up parameters
2230             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
2231             ? (ref($_[0]), @_)
2232             : objectify(2, @_);
2233 743 50       2347  
2234             return $x if $x->modify('bmod');
2235              
2236             # At least one argument is NaN. This is handled the same way as in
2237             # Math::BigInt -> bmod().
2238 743 100 100     2112  
2239             return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan();
2240              
2241             # Modulo zero. This is handled the same way as in Math::BigInt -> bmod().
2242 707 100       1701  
2243 40         131 if ($y -> is_zero()) {
2244             return $x -> round(@r);
2245             }
2246              
2247             # Numerator (dividend) is +/-inf. This is handled the same way as in
2248             # Math::BigInt -> bmod().
2249 667 100       1781  
2250 48         168 if ($x -> is_inf()) {
2251             return $x -> bnan(@r);
2252             }
2253              
2254             # Denominator (divisor) is +/-inf. This is handled the same way as in
2255             # Math::BigInt -> bmod().
2256 619 100       1332  
2257 40 100 100     109 if ($y -> is_inf()) {
2258 28         102 if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) {
2259             return $x -> round(@r);
2260 12         64 } else {
2261             return $x -> binf($y -> sign(), @r);
2262             }
2263             }
2264              
2265             return $x->bzero(@r) if $x->is_zero()
2266             || ($x->is_int() &&
2267 579 100 100     1146 # check that $y == +1 or $y == -1:
      100        
      100        
2268             ($LIB->_is_zero($y->{_e}) && $LIB->_is_one($y->{_m})));
2269 483         1444  
2270 483 100       1099 my $cmp = $x->bacmp($y); # equal or $x < $y?
2271 30         131 if ($cmp == 0) { # $x == $y => result 0
2272             return $x -> bzero(@r);
2273             }
2274              
2275 453 100       1081 # only $y of the operands negative?
2276             my $neg = $x->{sign} ne $y->{sign} ? 1 : 0;
2277 453         790  
2278 453 100 100     1142 $x->{sign} = $y->{sign}; # calc sign first
2279 48         196 if ($cmp < 0 && $neg == 0) { # $x < $y => result $x
2280             return $x -> round(@r);
2281             }
2282 405         1083  
2283             my $ym = $LIB->_copy($y->{_m});
2284              
2285             # 2e1 => 20
2286 405 100 100     1508 $ym = $LIB->_lsft($ym, $y->{_e}, 10)
2287             if $y->{_es} eq '+' && !$LIB->_is_zero($y->{_e});
2288              
2289 405         689 # if $y has digits after dot
2290 405 100       938 my $shifty = 0; # correct _e of $x by this
2291             if ($y->{_es} eq '-') # has digits after dot
2292             {
2293 16         92 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
2294             $shifty = $LIB->_num($y->{_e}); # no more digits after dot
2295 16         75 # 123 => 1230, $y->{_m} is already 25
2296             $x->{_m} = $LIB->_lsft($x->{_m}, $y->{_e}, 10);
2297             }
2298             # $ym is now mantissa of $y based on exponent 0
2299 405         607  
2300 405 100       856 my $shiftx = 0; # correct _e of $x by this
2301             if ($x->{_es} eq '-') # has digits after dot
2302             {
2303 27         97 # 123.4 % 20 => 1234 % 200
2304 27         98 $shiftx = $LIB->_num($x->{_e}); # no more digits after dot
2305             $ym = $LIB->_lsft($ym, $x->{_e}, 10); # 123 => 1230
2306             }
2307 405 100 100     1410 # 123e1 % 20 => 1230 % 20
2308 117         385 if ($x->{_es} eq '+' && !$LIB->_is_zero($x->{_e})) {
2309             $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # es => '+' here
2310             }
2311 405         1060  
2312 405         796 $x->{_e} = $LIB->_new($shiftx);
2313 405 100 100     1508 $x->{_es} = '+';
2314 405 100       870 $x->{_es} = '-' if $shiftx != 0 || $shifty != 0;
2315             $x->{_e} = $LIB->_add($x->{_e}, $LIB->_new($shifty)) if $shifty != 0;
2316              
2317             # now mantissas are equalized, exponent of $x is adjusted, so calc result
2318 405         1339  
2319             $x->{_m} = $LIB->_mod($x->{_m}, $ym);
2320 405 100       1039  
2321 405         1022 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # fix sign for -0
2322             $x = $x->bnorm();
2323              
2324 405 100 66     1101 # if one of them negative => correct in place
2325 51         174 if ($neg != 0 && ! $x -> is_zero()) {
2326 51         121 my $r = $y - $x;
2327 51         117 $x->{_m} = $r->{_m};
2328 51         104 $x->{_e} = $r->{_e};
2329 51 50       147 $x->{_es} = $r->{_es};
2330 51         140 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # fix sign for -0
2331             $x = $x->bnorm();
2332             }
2333 405         1719  
2334 405 0 0     1180 $x = $x->round($r[0], $r[1], $r[2], $y);
      33        
2335             return $downgrade -> new($x -> bdstr(), @r)
2336 405         3376 if defined($downgrade) && ($x->is_int() || $x->is_inf() || $x->is_nan());
2337             return $x;
2338             }
2339              
2340             sub bmodpow {
2341             # takes a very large number to a very large exponent in a given very
2342             # large modulus, quickly, thanks to binary exponentiation. Supports
2343 20 50 33 20 1 427 # negative exponents.
2344             my ($class, $num, $exp, $mod, @r)
2345             = ref($_[0]) && ref($_[0]) eq ref($_[1]) && ref($_[1]) eq ref($_[2])
2346             ? (ref($_[0]), @_)
2347             : objectify(3, @_);
2348 20 50       85  
2349             return $num if $num->modify('bmodpow');
2350 20 50 33     67  
      33        
2351             return $num -> bnan(@r)
2352             if $mod->is_nan() || $exp->is_nan() || $mod->is_nan();
2353              
2354 20 50 33     100 # check modulus for valid values
2355             return $num->bnan(@r) if $mod->{sign} ne '+' || $mod->is_zero();
2356              
2357 20 50       89 # check exponent for valid values
2358             if ($exp->{sign} =~ /\w/) {
2359 0         0 # i.e., if it's NaN, +inf, or -inf...
2360             return $num->bnan(@r);
2361             }
2362 20 50       76  
2363             $num = $num->bmodinv($mod, @r) if $exp->{sign} eq '-';
2364              
2365 20 50       77 # check num for valid values (also NaN if there was no inverse but $exp < 0)
2366             return $num->bnan(@r) if $num->{sign} !~ /^[+-]$/;
2367              
2368             # $mod is positive, sign on $exp is ignored, result also positive
2369              
2370 20         65 # XXX TODO: speed it up when all three numbers are integers
2371             $num = $num->bpow($exp)->bmod($mod);
2372 20 0 0     69  
      33        
2373             return $downgrade -> new($num -> bdstr(), @r) if defined($downgrade)
2374 20         66 && ($num->is_int() || $num->is_inf() || $num->is_nan());
2375             return $num -> round(@r);
2376             }
2377              
2378             sub bpow {
2379             # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
2380             # compute power of two numbers, second arg is used as integer
2381             # modifies first argument
2382              
2383 1042     1042 1 12950 # set up parameters
2384             my ($class, $x, $y, $a, $p, $r) = (ref($_[0]), @_);
2385 1042 100 100     4645 # objectify is costly, so avoid it
2386 110         416 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) {
2387             ($class, $x, $y, $a, $p, $r) = objectify(2, @_);
2388             }
2389 1042 50       3234  
2390             return $x if $x -> modify('bpow');
2391              
2392 1042 100 100     2678 # $x and/or $y is a NaN
2393             return $x -> bnan() if $x -> is_nan() || $y -> is_nan();
2394              
2395 926 100       2511 # $x and/or $y is a +/-Inf
    100          
    100          
    100          
2396 60 100       196 if ($x -> is_inf("-")) {
2397 32 100       144 return $x -> bzero() if $y -> is_negative();
2398 28 100       104 return $x -> bnan() if $y -> is_zero();
2399 20         94 return $x if $y -> is_odd();
2400             return $x -> bneg();
2401 60 100       204 } elsif ($x -> is_inf("+")) {
2402 32 100       105 return $x -> bzero() if $y -> is_negative();
2403 28         369 return $x -> bnan() if $y -> is_zero();
2404             return $x;
2405 44 100       223 } elsif ($y -> is_inf("-")) {
2406 40 100 100     165 return $x -> bnan() if $x -> is_one("-");
2407 28 100       122 return $x -> binf("+") if $x > -1 && $x < 1;
2408 24         153 return $x -> bone() if $x -> is_one("+");
2409             return $x -> bzero();
2410 44 100       269 } elsif ($y -> is_inf("+")) {
2411 40 100 100     242 return $x -> bnan() if $x -> is_one("-");
2412 28 100       130 return $x -> bzero() if $x > -1 && $x < 1;
2413 24         121 return $x -> bone() if $x -> is_one("+");
2414             return $x -> binf("+");
2415             }
2416 718 100       2583  
2417 48 100       135 if ($x -> is_zero()) {
2418 44 100       155 return $x -> bone() if $y -> is_zero();
2419 24         289 return $x -> binf() if $y -> is_negative();
2420             return $x;
2421             }
2422              
2423             # We don't support complex numbers, so upgrade or return NaN.
2424 670 100 100     2450  
2425 80 50       196 if ($x -> is_negative() && !$y -> is_int()) {
2426 80         227 return $upgrade -> bpow($x, $y, $a, $p, $r) if defined $upgrade;
2427             return $x -> bnan();
2428             }
2429 590 100 100     1649  
2430 119         1169 if ($x -> is_one("+") || $y -> is_one()) {
2431             return $x;
2432             }
2433 471 100       1206  
2434 24 100       80 if ($x -> is_one("-")) {
2435 12         48 return $x if $y -> is_odd();
2436             return $x -> bneg();
2437             }
2438 447 100       1298  
2439             return $x -> _pow($y, $a, $p, $r) if !$y -> is_int();
2440 333         1050  
2441             my $y1 = $y -> as_int()->{value}; # make MBI part
2442 333         935  
2443 333 100       1171 my $new_sign = '+';
    100          
2444             $new_sign = $LIB -> _is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
2445              
2446 333         1291 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
2447 333         1184 $x->{_m} = $LIB -> _pow($x->{_m}, $y1);
2448             $x->{_e} = $LIB -> _mul($x->{_e}, $y1);
2449 333         720  
2450 333         1049 $x->{sign} = $new_sign;
2451             $x = $x -> bnorm();
2452              
2453             # x ** (-y) = 1 / (x ** y)
2454 333 100       1229  
2455             if ($y->{sign} eq '-') {
2456 105         347 # modify $x in place!
2457 105         386 my $z = $x -> copy();
2458             $x = $x -> bone();
2459 105         408 # round in one go (might ignore y's A!)
2460             return scalar $x -> bdiv($z, $a, $p, $r);
2461             }
2462 228         930  
2463             $x = $x -> round($a, $p, $r, $y);
2464 228 50 66     833  
      66        
2465             return $downgrade -> new($x)
2466 226         2417 if defined($downgrade) && ($x->is_int() || $x->is_inf() || $x->is_nan());
2467             return $x;
2468             }
2469              
2470             sub blog {
2471             # Return the logarithm of the operand. If a second operand is defined, that
2472             # value is used as the base, otherwise the base is assumed to be Euler's
2473             # constant.
2474 256     256 1 1742  
2475             my ($class, $x, $base, @r);
2476              
2477             # Only objectify the base if it is defined, since an undefined base, as in
2478             # $x->blog() or $x->blog(undef) signals that the base is Euler's number =
2479             # 2.718281828...
2480 256 100 66     1200  
2481             if (!ref($_[0]) && $_[0] =~ /^[A-Za-z]|::/) {
2482 18 100       86 # E.g., Math::BigFloat->blog(256, 2)
2483             ($class, $x, $base, @r) =
2484             defined $_[2] ? objectify(2, @_) : objectify(1, @_);
2485             } else {
2486 238 100       1161 # E.g., $x->blog(2) or the deprecated Math::BigFloat::blog(256, 2)
2487             ($class, $x, $base, @r) =
2488             defined $_[1] ? objectify(2, @_) : objectify(1, @_);
2489             }
2490 256 50       1210  
2491             return $x if $x->modify('blog');
2492              
2493             # Handle all exception cases and all trivial cases. I have used Wolfram
2494             # Alpha (http://www.wolframalpha.com) as the reference for these cases.
2495 256 100       797  
2496             return $x -> bnan(@r) if $x -> is_nan();
2497 252 100       815  
2498 42 50 33     529 if (defined $base) {
2499             $base = $class -> new($base)
2500 42 100 66     139 unless defined(blessed($base)) && $base -> isa(__PACKAGE__);
    100 66        
    100          
2501 8         29 if ($base -> is_nan() || $base -> is_one()) {
2502             return $x -> bnan(@r);
2503 4 50 33     27 } elsif ($base -> is_inf() || $base -> is_zero()) {
2504 4         36 return $x -> bnan(@r) if $x -> is_inf() || $x -> is_zero();
2505             return $x -> bzero(@r);
2506 4 50       15 } elsif ($base -> is_negative()) { # -inf < base < 0
2507 4 50       32 return $x -> bzero(@r) if $x -> is_one(); # x = 1
2508             return $x -> bone('+', @r) if $x == $base; # x = base
2509 4 50       19 # we can't handle these cases, so upgrade, if we can
2510 4         19 return $upgrade -> blog($x, $base, @r) if defined $upgrade;
2511             return $x -> bnan(@r);
2512 26 100       98 }
2513             return $x -> bone(@r) if $x == $base; # 0 < base && 0 < x < inf
2514             }
2515 225 100       731  
    100          
    100          
    100          
2516 8 50 33     48 if ($x -> is_inf()) { # x = +/-inf
2517 8         32 my $sign = defined($base) && $base < 1 ? '-' : '+';
2518             return $x -> binf($sign, @r);
2519 16 50       59 } elsif ($x -> is_neg()) { # -inf < x < 0
2520 16         59 return $upgrade -> blog($x, $base, @r) if defined $upgrade;
2521             return $x -> bnan(@r);
2522 16         116 } elsif ($x -> is_one()) { # x = 1
2523             return $x -> bzero(@r);
2524 8 50 33     51 } elsif ($x -> is_zero()) { # x = 0
2525 8         33 my $sign = defined($base) && $base < 1 ? '+' : '-';
2526             return $x -> binf($sign, @r);
2527             }
2528              
2529 177         566 # we need to limit the accuracy to protect against overflow
2530 177         385 my $fallback = 0;
2531 177         804 my ($scale, @params);
2532             ($x, @params) = $x->_find_round_parameters(@r);
2533              
2534 177 100       726 # no rounding at all, so must use fallback
2535             if (scalar @params == 0) {
2536 91         389 # simulate old behaviour
2537 91         208 $params[0] = $class->div_scale(); # and round to it as accuracy
2538 91         168 $params[1] = undef; # P = undef
2539 91         163 $scale = $params[0]+4; # at least four more for proper round
2540 91         163 $params[2] = $r[2]; # round mode by caller or undef
2541             $fallback = 1; # to clear a/p afterwards
2542             } else {
2543             # the 4 below is empirical, and there might be cases where it is not
2544 86   33     687 # enough...
2545             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2546             }
2547              
2548             # when user set globals, they would interfere with our calculation, so
2549 43     43   508 # disable them and later re-enable them
  43         146  
  43         25793  
2550 177         510 no strict 'refs';
2551 177         448 my $abr = "$class\::accuracy";
2552 177         400 my $ab = $$abr;
2553 177         429 $$abr = undef;
2554 177         375 my $pbr = "$class\::precision";
2555 177         390 my $pb = $$pbr;
2556             $$pbr = undef;
2557             # we also need to disable any set A or P on $x (_find_round_parameters took
2558 177         448 # them already into account), since these would interfere, too
2559 177         372 $x->{_a} = undef;
2560             $x->{_p} = undef;
2561 177         319  
2562             my $done = 0;
2563              
2564             # If both $x and $base are integers, try to calculate an integer result
2565             # first. This is very fast, and if the exact result was found, we are done.
2566 177 50 66     609  
      66        
2567 11         38 if (defined($base) && $base -> is_int() && $x -> is_int()) {
2568 11         51 my $x_lib = $LIB -> _new($x -> bdstr());
2569 11         56 my $b_lib = $LIB -> _new($base -> bdstr());
2570 11 100       42 ($x_lib, my $exact) = $LIB -> _log_int($x_lib, $b_lib);
2571 10         27 if ($exact) {
2572 10         34 $x->{_m} = $x_lib;
2573 10         43 $x->{_e} = $LIB -> _zero();
2574 10         47 $x = $x -> bnorm();
2575             $done = 1;
2576             }
2577             }
2578              
2579             # If the integer result was not accurate, compute the natural logarithm
2580             # log($x) (using reduction by 10 and possibly also by 2), and if a
2581             # different base was requested, convert the result with log($x)/log($base).
2582 177 100       524  
2583 167         717 unless ($done) {
2584 167 100       615 $x = $x -> _log_10($scale);
2585             if (defined $base) {
2586 1         13 # log_b(x) = ln(x) / ln(b), so compute ln(b)
2587 1         5 my $base_log_e = $base -> copy() -> _log_10($scale);
2588             $x = $x -> bdiv($base_log_e, $scale);
2589             }
2590             }
2591              
2592             # shortcut to not run through _find_round_parameters again
2593 177 50       572  
2594 177         626 if (defined $params[0]) {
2595             $x = $x -> bround($params[0], $params[2]); # then round accordingly
2596 0         0 } else {
2597             $x = $x -> bfround($params[1], $params[2]); # then round accordingly
2598 177 100       886 }
2599             if ($fallback) {
2600 91         275 # clear a/p after round, since user did not request it
2601 91         212 $x->{_a} = undef;
2602             $x->{_p} = undef;
2603             }
2604 177         695 # restore globals
2605 177         465 $$abr = $ab;
2606             $$pbr = $pb;
2607 177 100 100     641  
2608             return $downgrade -> new($x -> bdstr(), @r)
2609 176         2537 if defined($downgrade) && $x -> is_int();
2610             return $x;
2611             }
2612              
2613             sub bexp {
2614 20 100   20 1 126 # Calculate e ** X (Euler's number to the power of X)
2615             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
2616 20 50       103  
2617             return $x if $x->modify('bexp');
2618 20 50       94  
2619 20 50       69 return $x->bnan(@r) if $x -> is_nan();
2620 20 50       50 return $x->binf(@r) if $x->{sign} eq '+inf';
2621             return $x->bzero(@r) if $x->{sign} eq '-inf';
2622              
2623 20         41 # we need to limit the accuracy to protect against overflow
2624 20         36 my $fallback = 0;
2625 20         74 my ($scale, @params);
2626             ($x, @params) = $x->_find_round_parameters(@r);
2627              
2628 20 50       1003 # error in _find_round_parameters?
2629             return $x->bnan(@r) if $x->{sign} eq 'NaN';
2630              
2631 20 100       71 # no rounding at all, so must use fallback
2632             if (scalar @params == 0) {
2633 11         60 # simulate old behaviour
2634 11         22 $params[0] = $class->div_scale(); # and round to it as accuracy
2635 11         23 $params[1] = undef; # P = undef
2636 11         16 $scale = $params[0]+4; # at least four more for proper round
2637 11         18 $params[2] = $r[2]; # round mode by caller or undef
2638             $fallback = 1; # to clear a/p afterwards
2639             } else {
2640             # the 4 below is empirical, and there might be cases where it's not
2641 9   33     36 # enough ...
2642             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2643             }
2644 20 50       97  
2645             return $x->bone(@params) if $x->is_zero();
2646 20 50       63  
2647 0         0 if (!$x->isa('Math::BigFloat')) {
2648 0         0 $x = Math::BigFloat->new($x);
2649             $class = ref($x);
2650             }
2651              
2652             # when user set globals, they would interfere with our calculation, so
2653 43     43   409 # disable them and later re-enable them
  43         157  
  43         42718  
2654 20         78 no strict 'refs';
2655 20         51 my $abr = "$class\::accuracy";
2656 20         48 my $ab = $$abr;
2657 20         44 $$abr = undef;
2658 20         44 my $pbr = "$class\::precision";
2659 20         50 my $pb = $$pbr;
2660             $$pbr = undef;
2661             # we also need to disable any set A or P on $x (_find_round_parameters took
2662 20         41 # them already into account), since these would interfere, too
2663 20         41 $x->{_a} = undef;
2664             $x->{_p} = undef;
2665              
2666             # Disabling upgrading and downgrading is no longer necessary to avoid an
2667             # infinite recursion, but it avoids unnecessary upgrading and downgrading in
2668             # the intermediate computations.
2669              
2670             # Temporarily disable downgrading
2671 20         78  
2672 20         67 my $dng = Math::BigFloat -> downgrade();
2673             Math::BigFloat -> downgrade(undef);
2674 20         70  
2675             my $x_org = $x->copy();
2676              
2677             # We use the following Taylor series:
2678              
2679             # x x^2 x^3 x^4
2680             # e = 1 + --- + --- + --- + --- ...
2681             # 1! 2! 3! 4!
2682              
2683             # The difference for each term is X and N, which would result in:
2684             # 2 copy, 2 mul, 2 add, 1 inc, 1 div operations per term
2685              
2686             # But it is faster to compute exp(1) and then raising it to the
2687             # given power, esp. if $x is really big and an integer because:
2688              
2689             # * The numerator is always 1, making the computation faster
2690             # * the series converges faster in the case of x == 1
2691             # * We can also easily check when we have reached our limit: when the
2692             # term to be added is smaller than "1E$scale", we can stop - f.i.
2693             # scale == 5, and we have 1/40320, then we stop since 1/40320 < 1E-5.
2694             # * we can compute the *exact* result by simulating bigrat math:
2695              
2696             # 1 1 gcd(3, 4) = 1 1*24 + 1*6 5
2697             # - + - = ---------- = --
2698             # 6 24 6*24 24
2699              
2700             # We do not compute the gcd() here, but simple do:
2701             # 1 1 1*24 + 1*6 30
2702             # - + - = --------- = --
2703             # 6 24 6*24 144
2704              
2705             # In general:
2706             # a c a*d + c*b and note that c is always 1 and d = (b*f)
2707             # - + - = ---------
2708             # b d b*d
2709              
2710             # This leads to: which can be reduced by b to:
2711             # a 1 a*b*f + b a*f + 1
2712             # - + - = --------- = -------
2713             # b b*f b*b*f b*f
2714              
2715             # The first terms in the series are:
2716              
2717             # 1 1 1 1 1 1 1 1 13700
2718             # -- + -- + -- + -- + -- + --- + --- + ---- = -----
2719             # 1 1 2 6 24 120 720 5040 5040
2720              
2721             # Note that we cannot simply reduce 13700/5040 to 685/252, but must keep
2722             # the numerator and the denominator!
2723 20 100       88  
2724             if ($scale <= 75) {
2725 17         94 # set $x directly from a cached string form
2726             $x->{_m} = $LIB->_new("2718281828459045235360287471352662497757" .
2727 17         63 "2470936999595749669676277240766303535476");
2728 17         33 $x->{sign} = '+';
2729 17         98 $x->{_es} = '-';
2730             $x->{_e} = $LIB->_new(79);
2731             } else {
2732             # compute A and B so that e = A / B.
2733              
2734             # After some terms we end up with this, so we use it as a starting
2735 3         13 # point:
2736             my $A = $LIB->_new("9093339520860578540197197" .
2737 3         12 "0164779391644753259799242");
2738 3         8 my $F = $LIB->_new(42);
2739             my $step = 42;
2740              
2741             # Compute how many steps we need to take to get $A and $B sufficiently
2742 3         12 # big
2743             my $steps = _len_to_steps($scale - 4);
2744 3         26 # print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
2745             while ($step++ <= $steps) {
2746 79         152 # calculate $a * $f + 1
2747 79         182 $A = $LIB->_mul($A, $F);
2748             $A = $LIB->_inc($A);
2749 79         155 # increment f
2750             $F = $LIB->_inc($F);
2751             }
2752              
2753             # Compute $B as factorial of $steps (this is faster than doing it
2754 3         16 # manually)
2755             my $B = $LIB->_fac($LIB->_new($steps));
2756              
2757             # print "A ", $LIB->_str($A), "\nB ", $LIB->_str($B), "\n";
2758              
2759 3         10 # compute A/B with $scale digits in the result (truncate, not round)
2760 3         14 $A = $LIB->_lsft($A, $LIB->_new($scale), 10);
2761             $A = $LIB->_div($A, $B);
2762 3         12  
2763 3         10 $x->{_m} = $A;
2764 3         7 $x->{sign} = '+';
2765 3         11 $x->{_es} = '-';
2766             $x->{_e} = $LIB->_new($scale);
2767             }
2768              
2769             # $x contains now an estimate of e, with some surplus digits, so we can
2770 20 100       72 # round
2771             if (!$x_org->is_one()) {
2772 10         24 # Reduce size of fractional part, followup with integer power of two.
2773 10   66     66 my $lshift = 0;
2774 21         71 while ($lshift < 30 && $x_org->bacmp(2 << $lshift) > 0) {
2775             $lshift++;
2776             }
2777 10 100       41 # Raise $x to the wanted power and round it.
2778 5         23 if ($lshift == 0) {
2779             $x = $x->bpow($x_org, @params);
2780 5         27 } else {
2781 5         29 my($mul, $rescale) = (1 << $lshift, $scale+1+$lshift);
2782             $x = $x -> bpow(scalar $x_org->bdiv($mul, $rescale), $rescale)
2783             -> bpow($mul, @params);
2784             }
2785             } else {
2786 10         27 # else just round the already computed result
2787 10         23 $x->{_a} = undef;
2788             $x->{_p} = undef;
2789 10 50       45 # shortcut to not run through _find_round_parameters again
2790 10         40 if (defined $params[0]) {
2791             $x = $x->bround($params[0], $params[2]); # then round accordingly
2792 0         0 } else {
2793             $x = $x->bfround($params[1], $params[2]); # then round accordingly
2794             }
2795             }
2796 20 100       210  
2797             if ($fallback) {
2798 11         31 # clear a/p after round, since user did not request it
2799 11         25 $x->{_a} = undef;
2800             $x->{_p} = undef;
2801             }
2802              
2803 20         72 # Restore globals
2804 20         52 $$abr = $ab;
2805             $$pbr = $pb;
2806              
2807             # Restore downgrading.
2808 20         88  
2809             Math::BigFloat -> downgrade($dng);
2810 20 50 33     219  
2811             return $downgrade -> new($x -> bdstr(), @r)
2812 20         193 if defined($downgrade) && $x -> is_int();
2813             $x;
2814             }
2815              
2816             sub bnok {
2817             # Calculate n over k (binomial coefficient or "choose" function) as integer.
2818 60 50 33 60 1 1010 # set up parameters
2819             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
2820             ? (ref($_[0]), @_)
2821             : objectify(2, @_);
2822 60 50       169  
2823             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
2824 60 50       202  
2825             return $x if $x->modify('bnok');
2826 60 100 100     170  
2827 48 50 66     174 return $x->bnan() if $x->is_nan() || $y->is_nan();
      33        
      33        
2828             return $x->bnan() if (($x->is_finite() && !$x->is_int()) ||
2829             ($y->is_finite() && !$y->is_int()));
2830 48         185  
2831 48         158 my $xint = Math::BigInt -> new($x -> bsstr());
2832 48         195 my $yint = Math::BigInt -> new($y -> bsstr());
2833             $xint = $xint -> bnok($yint);
2834 48 50       129  
2835             return $xint if defined $downgrade;
2836 48         138  
2837             my $xflt = Math::BigFloat -> new($xint);
2838 48         175  
2839 48         94 $x->{_m} = $xflt->{_m};
2840 48         79 $x->{_e} = $xflt->{_e};
2841 48         95 $x->{_es} = $xflt->{_es};
2842             $x->{sign} = $xflt->{sign};
2843 48         652  
2844             return $x;
2845             }
2846              
2847             sub bsin {
2848 40 50   40 1 548 # Calculate a sinus of x.
2849             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
2850              
2851             # taylor: x^3 x^5 x^7 x^9
2852             # sin = x - --- + --- - --- + --- ...
2853             # 3! 5! 7! 9!
2854 40 50       166  
2855             return $x if $x->modify('bsin');
2856 40 100       101  
2857 32 100 100     102 return $x -> bzero(@r) if $x->is_zero();
2858             return $x -> bnan(@r) if $x->is_nan() || $x->is_inf();
2859              
2860 20         52 # we need to limit the accuracy to protect against overflow
2861 20         39 my $fallback = 0;
2862 20         58 my ($scale, @params);
2863             ($x, @params) = $x->_find_round_parameters(@r);
2864              
2865 20 50       73 # error in _find_round_parameters?
2866             return $x->bnan(@r) if $x->is_nan();
2867              
2868 20 50       78 # no rounding at all, so must use fallback
2869             if (scalar @params == 0) {
2870 0         0 # simulate old behaviour
2871 0         0 $params[0] = $class->div_scale(); # and round to it as accuracy
2872 0         0 $params[1] = undef; # disable P
2873 0         0 $scale = $params[0]+4; # at least four more for proper round
2874 0         0 $params[2] = $r[2]; # round mode by caller or undef
2875             $fallback = 1; # to clear a/p afterwards
2876             } else {
2877             # the 4 below is empirical, and there might be cases where it is not
2878 20   33     62 # enough...
2879             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2880             }
2881              
2882             # when user set globals, they would interfere with our calculation, so
2883 43     43   400 # disable them and later re-enable them
  43         131  
  43         24523  
2884 20         51 no strict 'refs';
2885 20         63 my $abr = "$class\::accuracy";
2886 20         42 my $ab = $$abr;
2887 20         41 $$abr = undef;
2888 20         42 my $pbr = "$class\::precision";
2889 20         52 my $pb = $$pbr;
2890             $$pbr = undef;
2891             # we also need to disable any set A or P on $x (_find_round_parameters took
2892 20         42 # them already into account), since these would interfere, too
2893 20         38 $x->{_a} = undef;
2894             $x->{_p} = undef;
2895              
2896             # Disabling upgrading and downgrading is no longer necessary to avoid an
2897             # infinite recursion, but it avoids unnecessary upgrading and downgrading in
2898             # the intermediate computations.
2899 20         52  
2900 20         38 local $Math::BigInt::upgrade = undef;
2901             local $Math::BigFloat::downgrade = undef;
2902 20         64  
2903 20         72 my $over = $x * $x; # X ^ 2
2904 20         65 my $x2 = $over->copy(); # X ^ 2; difference between terms
2905 20         47 $over = $over->bmul($x); # X ^ 3 as starting value
2906 20         57 my $sign = 1; # start with -=
2907 20         85 my $below = $class->new(6);
2908 20         76 my $factorial = $class->new(4);
2909 20         41 $x->{_a} = undef;
2910             $x->{_p} = undef;
2911 20         73  
2912 20         66 my $limit = $class->new("1E-". ($scale-1));
2913             while (1) {
2914             # we calculate the next term, and add it to the last
2915             # when the next term is below our limit, it won't affect the outcome
2916 196         471 # anymore, so we stop:
2917 196 100       507 my $next = $over->copy()->bdiv($below, $scale);
2918             last if $next->bacmp($limit) <= 0;
2919 176 100       367  
2920 80         196 if ($sign == 0) {
2921             $x = $x->badd($next);
2922 96         248 } else {
2923             $x = $x->bsub($next);
2924 176         270 }
2925             $sign = 1-$sign; # alternate
2926 176         397 # calculate things for the next term
2927 176         385 $over = $over->bmul($x2); # $x*$x
2928 176         414 $below = $below->bmul($factorial); # n*(n+1)
2929 176         392 $factorial = $factorial->binc();
2930 176         368 $below = $below -> bmul($factorial); # n*(n+1)
2931             $factorial = $factorial->binc();
2932             }
2933              
2934 20 50       67 # shortcut to not run through _find_round_parameters again
2935 20         90 if (defined $params[0]) {
2936             $x = $x->bround($params[0], $params[2]); # then round accordingly
2937 0         0 } else {
2938             $x = $x->bfround($params[1], $params[2]); # then round accordingly
2939 20 50       77 }
2940             if ($fallback) {
2941 0         0 # clear a/p after round, since user did not request it
2942 0         0 $x->{_a} = undef;
2943             $x->{_p} = undef;
2944             }
2945 20         67 # restore globals
2946 20         74 $$abr = $ab;
2947             $$pbr = $pb;
2948 20 50 33     126  
2949             return $downgrade -> new($x -> bdstr(), @r)
2950 20         432 if defined($downgrade) && $x -> is_int();
2951             $x;
2952             }
2953              
2954             sub bcos {
2955 36 50   36 1 603 # Calculate a cosinus of x.
2956             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
2957              
2958             # Taylor: x^2 x^4 x^6 x^8
2959             # cos = 1 - --- + --- - --- + --- ...
2960             # 2! 4! 6! 8!
2961              
2962 36         60 # we need to limit the accuracy to protect against overflow
2963 36         74 my $fallback = 0;
2964 36         125 my ($scale, @params);
2965             ($x, @params) = $x->_find_round_parameters(@r);
2966              
2967 36 100 66     222 # constant object or error in _find_round_parameters?
2968 32 100       96 return $x if $x->modify('bcos') || $x->is_nan();
2969 24 100       93 return $x->bnan() if $x->is_inf();
2970             return $x->bone(@r) if $x->is_zero();
2971              
2972 16 50       82 # no rounding at all, so must use fallback
2973             if (scalar @params == 0) {
2974 0         0 # simulate old behaviour
2975 0         0 $params[0] = $class->div_scale(); # and round to it as accuracy
2976 0         0 $params[1] = undef; # disable P
2977 0         0 $scale = $params[0]+4; # at least four more for proper round
2978 0         0 $params[2] = $r[2]; # round mode by caller or undef
2979             $fallback = 1; # to clear a/p afterwards
2980             } else {
2981             # the 4 below is empirical, and there might be cases where it is not
2982 16   33     64 # enough...
2983             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2984             }
2985              
2986             # when user set globals, they would interfere with our calculation, so
2987 43     43   369 # disable them and later re-enable them
  43         129  
  43         37641  
2988 16         41 no strict 'refs';
2989 16         52 my $abr = "$class\::accuracy";
2990 16         40 my $ab = $$abr;
2991 16         35 $$abr = undef;
2992 16         34 my $pbr = "$class\::precision";
2993 16         43 my $pb = $$pbr;
2994             $$pbr = undef;
2995             # we also need to disable any set A or P on $x (_find_round_parameters took
2996 16         40 # them already into account), since these would interfere, too
2997 16         24 $x->{_a} = undef;
2998             $x->{_p} = undef;
2999 16         57  
3000 16         50 my $over = $x * $x; # X ^ 2
3001 16         31 my $x2 = $over->copy(); # X ^ 2; difference between terms
3002 16         68 my $sign = 1; # start with -=
3003 16         71 my $below = $class->new(2);
3004 16         114 my $factorial = $class->new(3);
3005 16         36 $x = $x->bone();
3006 16         35 $x->{_a} = undef;
3007             $x->{_p} = undef;
3008 16         59  
3009             my $limit = $class->new("1E-". ($scale-1));
3010 16         53 #my $steps = 0;
3011             while (3 < 5) {
3012             # we calculate the next term, and add it to the last
3013             # when the next term is below our limit, it won't affect the outcome
3014 156         388 # anymore, so we stop:
3015 156 100       396 my $next = $over->copy()->bdiv($below, $scale);
3016             last if $next->bacmp($limit) <= 0;
3017 140 100       271  
3018 68         197 if ($sign == 0) {
3019             $x = $x->badd($next);
3020 72         188 } else {
3021             $x = $x->bsub($next);
3022 140         234 }
3023             $sign = 1-$sign; # alternate
3024 140         328 # calculate things for the next term
3025 140         311 $over = $over->bmul($x2); # $x*$x
3026 140         349 $below = $below->bmul($factorial); # n*(n+1)
3027 140         313 $factorial = $factorial -> binc();
3028 140         316 $below = $below->bmul($factorial); # n*(n+1)
3029             $factorial = $factorial -> binc();
3030             }
3031              
3032 16 50       49 # shortcut to not run through _find_round_parameters again
3033 16         50 if (defined $params[0]) {
3034             $x = $x->bround($params[0], $params[2]); # then round accordingly
3035 0         0 } else {
3036             $x = $x->bfround($params[1], $params[2]); # then round accordingly
3037 16 50       62 }
3038             if ($fallback) {
3039 0         0 # clear a/p after round, since user did not request it
3040 0         0 $x->{_a} = undef;
3041             $x->{_p} = undef;
3042             }
3043 16         65 # restore globals
3044 16         45 $$abr = $ab;
3045             $$pbr = $pb;
3046 16 50 33     58  
3047             return $downgrade -> new($x -> bdstr(), @r)
3048 16         326 if defined($downgrade) && $x -> is_int();
3049             $x;
3050             }
3051              
3052             sub batan {
3053 175 50   175 1 1616 # Calculate a arcus tangens of x.
3054             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
3055              
3056             # taylor: x^3 x^5 x^7 x^9
3057             # atan = x - --- + --- - --- + --- ...
3058             # 3 5 7 9
3059 175 50       675  
3060             return $x if $x->modify('batan');
3061 175 100       536  
3062             return $x -> bnan(@r) if $x->is_nan();
3063              
3064             # We need to limit the accuracy to protect against overflow.
3065 171         364  
3066 171         326 my $fallback = 0;
3067 171         534 my ($scale, @params);
3068             ($x, @params) = $x->_find_round_parameters(@r);
3069              
3070             # Error in _find_round_parameters?
3071 171 50       640  
3072             return $x -> bnan(@r) if $x->is_nan();
3073 171 100       617  
3074             if ($x->{sign} =~ /^[+-]inf\z/) {
3075             # +inf result is PI/2
3076             # -inf result is -PI/2
3077 16         74 # calculate PI/2
3078             my $pi = $class->bpi(@r);
3079 16         88 # modify $x in place
3080 16         51 $x->{_m} = $pi->{_m};
3081 16         43 $x->{_e} = $pi->{_e};
3082             $x->{_es} = $pi->{_es};
3083 16         73 # -y => -PI/2, +y => PI/2
3084 16         68 $x->{sign} = substr($x->{sign}, 0, 1); # "+inf" => "+"
3085 16         186 $x -> {_m} = $LIB->_div($x->{_m}, $LIB->_new(2));
3086             return $x;
3087             }
3088 155 100       445  
3089             return $x->bzero(@r) if $x->is_zero();
3090              
3091 129 50       576 # no rounding at all, so must use fallback
3092             if (scalar @params == 0) {
3093 0         0 # simulate old behaviour
3094 0         0 $params[0] = $class->div_scale(); # and round to it as accuracy
3095 0         0 $params[1] = undef; # disable P
3096 0         0 $scale = $params[0]+4; # at least four more for proper round
3097 0         0 $params[2] = $r[2]; # round mode by caller or undef
3098             $fallback = 1; # to clear a/p afterwards
3099             } else {
3100             # the 4 below is empirical, and there might be cases where it is not
3101 129   33     530 # enough...
3102             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3103             }
3104              
3105             # 1 or -1 => PI/4
3106 129 100 100     477 # inlined is_one() && is_one('-')
3107 27         155 if ($LIB->_is_one($x->{_m}) && $LIB->_is_zero($x->{_e})) {
3108             my $pi = $class->bpi($scale - 3);
3109 27         86 # modify $x in place
3110 27         73 $x->{_m} = $pi->{_m};
3111 27         73 $x->{_e} = $pi->{_e};
3112             $x->{_es} = $pi->{_es};
3113 27         110 # leave the sign of $x alone (+1 => +PI/4, -1 => -PI/4)
3114 27         212 $x->{_m} = $LIB->_div($x->{_m}, $LIB->_new(4));
3115             return $x;
3116             }
3117              
3118             # This series is only valid if -1 < x < 1, so for other x we need to
3119 102         218 # calculate PI/2 - atan(1/x):
3120 102 100       339 my $pi = undef;
3121             if ($x->bacmp($x->copy()->bone) >= 0) {
3122 40         249 # calculate PI/2
3123 40         190 $pi = $class->bpi($scale - 3);
3124             $pi->{_m} = $LIB->_div($pi->{_m}, $LIB->_new(2));
3125 40         186 # calculate 1/$x:
3126             my $x_copy = $x->copy();
3127 40         157 # modify $x in place
3128 40         207 $x = $x->bone();
3129             $x = $x->bdiv($x_copy, $scale);
3130             }
3131 102         445  
3132 102         476 my $fmul = 1;
3133 174         323 foreach (0 .. int($scale / 20)) {
3134 174         489 $fmul *= 2;
3135             $x = $x->bdiv($x->copy()->bmul($x)->binc()->bsqrt($scale + 4)->binc(),
3136             $scale + 4);
3137             }
3138              
3139             # When user set globals, they would interfere with our calculation, so
3140 43     43   408 # disable them and later re-enable them.
  43         141  
  43         51635  
3141 102         386 no strict 'refs';
3142 102         366 my $abr = "$class\::accuracy";
3143 102         290 my $ab = $$abr;
3144 102         251 $$abr = undef;
3145 102         259 my $pbr = "$class\::precision";
3146 102         219 my $pb = $$pbr;
3147             $$pbr = undef;
3148             # We also need to disable any set A or P on $x (_find_round_parameters
3149 102         268 # took them already into account), since these would interfere, too
3150 102         206 $x->{_a} = undef;
3151             $x->{_p} = undef;
3152              
3153             # Disabling upgrading and downgrading is no longer necessary to avoid an
3154             # infinite recursion, but it avoids unnecessary upgrading and downgrading in
3155             # the intermediate computations.
3156 102         235  
3157 102         196 local $Math::BigInt::upgrade = undef;
3158             local $Math::BigFloat::downgrade = undef;
3159 102         372  
3160 102         412 my $over = $x * $x; # X ^ 2
3161 102         412 my $x2 = $over->copy(); # X ^ 2; difference between terms
3162 102         265 $over = $over->bmul($x); # X ^ 3 as starting value
3163 102         354 my $sign = 1; # start with -=
3164 102         520 my $below = $class->new(3);
3165 102         492 my $two = $class->new(2);
3166 102         216 $x->{_a} = undef;
3167             $x->{_p} = undef;
3168 102         411  
3169             my $limit = $class->new("1E-". ($scale-1));
3170 102         353 #my $steps = 0;
3171             while (1) {
3172             # We calculate the next term, and add it to the last. When the next
3173             # term is below our limit, it won't affect the outcome anymore, so we
3174 990         2491 # stop:
3175 990 100       2777 my $next = $over->copy()->bdiv($below, $scale);
3176             last if $next->bacmp($limit) <= 0;
3177 888 100       1897  
3178 416         1177 if ($sign == 0) {
3179             $x = $x->badd($next);
3180 472         1277 } else {
3181             $x = $x->bsub($next);
3182 888         1459 }
3183             $sign = 1-$sign; # alternatex
3184 888         2136 # calculate things for the next term
3185 888         2267 $over = $over->bmul($x2); # $x*$x
3186             $below = $below->badd($two); # n += 2
3187 102         464 }
3188             $x = $x->bmul($fmul);
3189 102 100       495  
3190 40         305 if (defined $pi) {
3191             my $x_copy = $x->copy();
3192 40         202 # modify $x in place
3193 40         102 $x->{_m} = $pi->{_m};
3194 40         93 $x->{_e} = $pi->{_e};
3195             $x->{_es} = $pi->{_es};
3196 40         112 # PI/2 - $x
3197             $x = $x->bsub($x_copy);
3198             }
3199              
3200 102 50       398 # Shortcut to not run through _find_round_parameters again.
3201 102         363 if (defined $params[0]) {
3202             $x = $x->bround($params[0], $params[2]); # then round accordingly
3203 0         0 } else {
3204             $x = $x->bfround($params[1], $params[2]); # then round accordingly
3205 102 50       461 }
3206             if ($fallback) {
3207 0         0 # Clear a/p after round, since user did not request it.
3208 0         0 $x->{_a} = undef;
3209             $x->{_p} = undef;
3210             }
3211              
3212 102         379 # restore globals
3213 102         277 $$abr = $ab;
3214             $$pbr = $pb;
3215 102 0 0     311  
      33        
3216             return $downgrade -> new($x -> bdstr(), @r)
3217 102         1989 if defined($downgrade) && ($x -> is_int() || $x -> is_inf());
3218             $x;
3219             }
3220              
3221             sub batan2 {
3222             # $y -> batan2($x) returns the arcus tangens of $y / $x.
3223              
3224 213 50 33 213 1 2986 # Set up parameters.
3225             my ($class, $y, $x, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3226             ? (ref($_[0]), @_)
3227             : objectify(2, @_);
3228              
3229 213 50       770 # Quick exit if $y is read-only.
3230             return $y if $y -> modify('batan2');
3231              
3232 213 100 100     964 # Handle all NaN cases.
3233             return $y -> bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
3234              
3235 201         375 # We need to limit the accuracy to protect against overflow.
3236 201         327 my $fallback = 0;
3237 201         616 my ($scale, @params);
3238             ($y, @params) = $y -> _find_round_parameters(@r);
3239              
3240 201 50       682 # Error in _find_round_parameters?
3241             return $y if $y->is_nan();
3242              
3243 201 100       604 # No rounding at all, so must use fallback.
3244             if (scalar @params == 0) {
3245 45         149 # Simulate old behaviour
3246 45         86 $params[0] = $class -> div_scale(); # and round to it as accuracy
3247 45         89 $params[1] = undef; # disable P
3248 45         85 $scale = $params[0] + 4; # at least four more for proper round
3249 45         74 $params[2] = $r[2]; # round mode by caller or undef
3250             $fallback = 1; # to clear a/p afterwards
3251             } else {
3252             # The 4 below is empirical, and there might be cases where it is not
3253 156   33     475 # enough ...
3254             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3255             }
3256 201 100       554  
    100          
    100          
    100          
3257 20 100       95 if ($x -> is_inf("+")) { # x = inf
    100          
3258 4         26 if ($y -> is_inf("+")) { # y = inf
3259             $y = $y -> bpi($scale) -> bmul("0.25"); # pi/4
3260 4         41 } elsif ($y -> is_inf("-")) { # y = -inf
3261             $y = $y -> bpi($scale) -> bmul("-0.25"); # -pi/4
3262 12         67 } else { # -inf < y < inf
3263             return $y -> bzero(@r); # 0
3264             }
3265 20 100       64 } elsif ($x -> is_inf("-")) { # x = -inf
    100          
    100          
3266 4         53 if ($y -> is_inf("+")) { # y = inf
3267             $y = $y -> bpi($scale) -> bmul("0.75"); # 3/4 pi
3268 4         33 } elsif ($y -> is_inf("-")) { # y = -inf
3269             $y = $y -> bpi($scale) -> bmul("-0.75"); # -3/4 pi
3270 8         33 } elsif ($y >= 0) { # y >= 0
3271             $y = $y -> bpi($scale); # pi
3272 4         29 } else { # y < 0
3273             $y = $y -> bpi($scale) -> bneg(); # -pi
3274             }
3275 87 100       294 } elsif ($x > 0) { # 0 < x < inf
    100          
3276 4         40 if ($y -> is_inf("+")) { # y = inf
3277             $y = $y -> bpi($scale) -> bmul("0.5"); # pi/2
3278 4         29 } elsif ($y -> is_inf("-")) { # y = -inf
3279             $y = $y -> bpi($scale) -> bmul("-0.5"); # -pi/2
3280 79         422 } else { # -inf < y < inf
3281             $y = $y -> bdiv($x, $scale) -> batan($scale); # atan(y/x)
3282             }
3283 20         80 } elsif ($x < 0) { # -inf < x < 0
3284 20 100       62 my $pi = $class -> bpi($scale);
3285 12         51 if ($y >= 0) { # y >= 0
3286             $y = $y -> bdiv($x, $scale) -> batan() # atan(y/x) + pi
3287             -> badd($pi);
3288 8         35 } else { # y < 0
3289             $y = $y -> bdiv($x, $scale) -> batan() # atan(y/x) - pi
3290             -> bsub($pi);
3291             }
3292 54 100       167 } else { # x = 0
    100          
3293 25         158 if ($y > 0) { # y > 0
3294             $y = $y -> bpi($scale) -> bmul("0.5"); # pi/2
3295 22         119 } elsif ($y < 0) { # y < 0
3296             $y = $y -> bpi($scale) -> bmul("-0.5"); # -pi/2
3297 7         80 } else { # y = 0
3298             return $y -> bzero(@r); # 0
3299             }
3300             }
3301 182         676  
3302             $y = $y -> round(@r);
3303 182 100       523  
3304 42         102 if ($fallback) {
3305 42         80 $y->{_a} = undef;
3306             $y->{_p} = undef;
3307             }
3308 182         2252  
3309             return $y;
3310             }
3311              
3312             sub bsqrt {
3313 442 100   442 1 2828 # calculate square root
3314             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
3315 442 50       1593  
3316             return $x if $x->modify('bsqrt');
3317              
3318             # Handle trivial cases.
3319 442 100       1250  
3320 434 100       1360 return $x -> bnan(@r) if $x->is_nan();
3321 430 100 100     999 return $x -> binf("+", @r) if $x->{sign} eq '+inf';
3322             return $x -> round(@r) if $x->is_zero() || $x->is_one();
3323              
3324             # We don't support complex numbers.
3325 420 100       1409  
3326 20 50       69 if ($x -> is_neg()) {
3327 20         63 return $upgrade -> bsqrt($x, @r) if defined($upgrade);
3328             return $x -> bnan(@r);
3329             }
3330              
3331 400         852 # we need to limit the accuracy to protect against overflow
3332 400         660 my $fallback = 0;
3333 400         1133 my (@params, $scale);
3334             ($x, @params) = $x->_find_round_parameters(@r);
3335              
3336 400 50       1163 # error in _find_round_parameters?
3337             return $x -> bnan(@r) if $x->is_nan();
3338              
3339 400 100       1171 # no rounding at all, so must use fallback
3340             if (scalar @params == 0) {
3341 123         407 # simulate old behaviour
3342 123         277 $params[0] = $class->div_scale(); # and round to it as accuracy
3343 123         240 $scale = $params[0]+4; # at least four more for proper round
3344 123         221 $params[2] = $r[2]; # round mode by caller or undef
3345             $fallback = 1; # to clear a/p afterwards
3346             } else {
3347             # the 4 below is empirical, and there might be cases where it is not
3348 277   100     749 # enough...
3349             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3350             }
3351              
3352             # when user set globals, they would interfere with our calculation, so
3353 43     43   427 # disable them and later re-enable them
  43         127  
  43         42100  
3354 400         937 no strict 'refs';
3355 400         890 my $abr = "$class\::accuracy";
3356 400         855 my $ab = $$abr;
3357 400         759 $$abr = undef;
3358 400         808 my $pbr = "$class\::precision";
3359 400         761 my $pb = $$pbr;
3360             $$pbr = undef;
3361             # we also need to disable any set A or P on $x (_find_round_parameters took
3362 400         747 # them already into account), since these would interfere, too
3363 400         741 $x->{_a} = undef;
3364             $x->{_p} = undef;
3365              
3366             # Disabling upgrading and downgrading is no longer necessary to avoid an
3367             # infinite recursion, but it avoids unnecessary upgrading and downgrading in
3368             # the intermediate computations.
3369 400         687  
3370 400         587 local $Math::BigInt::upgrade = undef;
3371             local $Math::BigFloat::downgrade = undef;
3372 400         1227  
3373 400 100       1360 my $i = $LIB->_copy($x->{_m});
3374 400         2027 $i = $LIB->_lsft($i, $x->{_e}, 10) unless $LIB->_is_zero($x->{_e});
3375 400         925 my $xas = Math::BigInt->bzero();
3376             $xas->{value} = $i;
3377 400         1535  
3378             my $gs = $xas->copy()->bsqrt(); # some guess
3379 400 100 100     1985  
3380             if (($x->{_es} ne '-') # guess can't be accurate if there are
3381             # digits after the dot
3382             && ($xas->bacmp($gs * $gs) == 0)) # guess hit the nail on the head?
3383             {
3384 50         190 # exact result, copy result over to keep $x
3385 50         197 $x->{_m} = $gs->{value};
3386 50         110 $x->{_e} = $LIB->_zero();
3387 50         132 $x->{_es} = '+';
3388             $x = $x->bnorm();
3389 50 50       153 # shortcut to not run through _find_round_parameters again
3390 50         154 if (defined $params[0]) {
3391             $x = $x->bround($params[0], $params[2]); # then round accordingly
3392 0         0 } else {
3393             $x = $x->bfround($params[1], $params[2]); # then round accordingly
3394 50 100       133 }
3395             if ($fallback) {
3396 33         62 # clear a/p after round, since user did not request it
3397 33         61 $x->{_a} = undef;
3398             $x->{_p} = undef;
3399             }
3400 50         96 # re-enable A and P, upgrade is taken care of by "local"
  50         157  
3401 50         76 ${"$class\::accuracy"} = $ab;
  50         114  
3402 50         541 ${"$class\::precision"} = $pb;
3403             return $x;
3404             }
3405              
3406             # sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the
3407             # accuracy of the result by multiplying the input by 100 and then divide the
3408             # integer result of sqrt(input) by 10. Rounding afterwards returns the real
3409             # result.
3410              
3411 350         1251 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
3412             my $y1 = $LIB->_copy($x->{_m});
3413 350         1053  
3414             my $length = $LIB->_len($y1);
3415              
3416 350         1033 # Now calculate how many digits the result of sqrt(y1) would have
3417             my $digits = int($length / 2);
3418              
3419 350         579 # But we need at least $scale digits, so calculate how many are missing
3420             my $shift = $scale - $digits;
3421              
3422             # This happens if the input had enough digits
3423 350 100       824 # (we take care of integer guesses above)
3424             $shift = 0 if $shift < 0;
3425              
3426 350         585 # Multiply in steps of 100, by shifting left two times the "missing" digits
3427             my $s2 = $shift * 2;
3428              
3429             # We now make sure that $y1 has the same odd or even number of digits than
3430             # $x had. So when _e of $x is odd, we must shift $y1 by one digit left,
3431             # because we always must multiply by steps of 100 (sqrt(100) is 10) and not
3432             # steps of 10. The length of $x does not count, since an even or odd number
3433             # of digits before the dot is not changed by adding an even number of digits
3434 350 100       1083 # after the dot (the result is still odd or even digits long).
3435             $s2++ if $LIB->_is_odd($x->{_e});
3436 350         1054  
3437             $y1 = $LIB->_lsft($y1, $LIB->_new($s2), 10);
3438              
3439 350         1363 # now take the square root and truncate to integer
3440             $y1 = $LIB->_sqrt($y1);
3441              
3442             # By "shifting" $y1 right (by creating a negative _e) we calculate the final
3443             # result, which is than later rounded to the desired scale.
3444              
3445             # calculate how many zeros $x had after the '.' (or before it, depending
3446 350         1267 # on sign of $dat, the result should have half as many:
3447 350 100       1283 my $dat = $LIB->_num($x->{_e});
3448 350         620 $dat = -$dat if $x->{_es} eq '-';
3449             $dat += $length;
3450 350 100       798  
3451             if ($dat > 0) {
3452             # no zeros after the dot (e.g. 1.23, 0.49 etc)
3453             # preserve half as many digits before the dot than the input had
3454 320         840 # (but round this "up")
3455             $dat = int(($dat+1)/2);
3456 30         107 } else {
3457             $dat = int(($dat)/2);
3458 350         922 }
3459 350 50       840 $dat -= $LIB->_len($y1);
3460 350         817 if ($dat < 0) {
3461 350         951 $dat = abs($dat);
3462 350         840 $x->{_e} = $LIB->_new($dat);
3463             $x->{_es} = '-';
3464 0         0 } else {
3465 0         0 $x->{_e} = $LIB->_new($dat);
3466             $x->{_es} = '+';
3467 350         747 }
3468 350         1071 $x->{_m} = $y1;
3469             $x = $x->bnorm();
3470              
3471 350 100       1096 # shortcut to not run through _find_round_parameters again
3472 328         1072 if (defined $params[0]) {
3473             $x = $x->bround($params[0], $params[2]); # then round accordingly
3474 22         88 } else {
3475             $x = $x->bfround($params[1], $params[2]); # then round accordingly
3476 350 100       1116 }
3477             if ($fallback) {
3478 90         201 # clear a/p after round, since user did not request it
3479 90         184 $x->{_a} = undef;
3480             $x->{_p} = undef;
3481             }
3482 350         1224 # restore globals
3483 350         893 $$abr = $ab;
3484             $$pbr = $pb;
3485 350 0 0     863  
      33        
3486             return $downgrade -> new($x -> bdstr(), @r)
3487 350         3152 if defined($downgrade) && ($x -> is_int() || $x -> is_inf());
3488             $x;
3489             }
3490              
3491             sub broot {
3492             # calculate $y'th root of $x
3493              
3494 186 100 100 186 1 3084 # set up parameters
3495             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3496             ? (ref($_[0]), @_)
3497             : objectify(2, @_);
3498 186 50       648  
3499             return $x if $x->modify('broot');
3500              
3501             # Handle trivial cases.
3502 186 100 100     525  
3503             return $x -> bnan(@r) if $x->is_nan() || $y->is_nan();
3504 150 100       432  
3505             if ($x -> is_neg()) {
3506 44 100 66     147 # -27 ** (1/3) = -3
      100        
3507             return $x -> broot($y -> copy() -> bneg(), @r) -> bneg()
3508 40 50       114 if $x -> is_int() && $y -> is_int() && $y -> is_neg();
3509 40         111 return $upgrade -> broot($x, $y, @r) if defined $upgrade;
3510             return $x -> bnan(@r);
3511             }
3512              
3513             # NaN handling: $x ** 1/0, x or y NaN, or y inf/-inf or y == 0
3514 106 100 66     561 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
      100        
3515             $y->{sign} !~ /^\+$/;
3516 82 100 100     185  
      100        
      100        
3517             return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
3518              
3519 66         175 # we need to limit the accuracy to protect against overflow
3520 66         114 my $fallback = 0;
3521 66         213 my (@params, $scale);
3522             ($x, @params) = $x->_find_round_parameters(@r);
3523 66 50       174  
3524             return $x if $x->is_nan(); # error in _find_round_parameters?
3525              
3526 66 50       179 # no rounding at all, so must use fallback
3527             if (scalar @params == 0) {
3528 66         202 # simulate old behaviour
3529 66         125 $params[0] = $class->div_scale(); # and round to it as accuracy
3530 66         107 $scale = $params[0]+4; # at least four more for proper round
3531 66         121 $params[2] = $r[2]; # round mode by caller or undef
3532             $fallback = 1; # to clear a/p afterwards
3533             } else {
3534             # the 4 below is empirical, and there might be cases where it is not
3535 0   0     0 # enough...
3536             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3537             }
3538              
3539             # when user set globals, they would interfere with our calculation, so
3540 43     43   394 # disable them and later re-enable them
  43         142  
  43         609879  
3541 66         140 no strict 'refs';
3542 66         143 my $abr = "$class\::accuracy";
3543 66         147 my $ab = $$abr;
3544 66         118 $$abr = undef;
3545 66         132 my $pbr = "$class\::precision";
3546 66         130 my $pb = $$pbr;
3547             $$pbr = undef;
3548             # we also need to disable any set A or P on $x (_find_round_parameters took
3549 66         112 # them already into account), since these would interfere, too
3550 66         117 $x->{_a} = undef;
3551             $x->{_p} = undef;
3552              
3553             # Disabling upgrading and downgrading is no longer necessary to avoid an
3554             # infinite recursion, but it avoids unnecessary upgrading and downgrading in
3555             # the intermediate computations.
3556 66         129  
3557 66         105 local $Math::BigInt::upgrade = undef;
3558             local $Math::BigFloat::downgrade = undef;
3559              
3560 66         96 # remember sign and make $x positive, since -4 ** (1/2) => -2
3561 66 50       176 my $sign = 0;
3562 66         116 $sign = 1 if $x->{sign} eq '-';
3563             $x->{sign} = '+';
3564 66         93  
3565 66 50       168 my $is_two = 0;
3566             if ($y->isa('Math::BigFloat')) {
3567 66   66     345 $is_two = $y->{sign} eq '+' && $LIB->_is_two($y->{_m})
3568             && $LIB->_is_zero($y->{_e});
3569 0         0 } else {
3570             $is_two = $y == 2;
3571             }
3572              
3573 66 100       205 # normal square root if $y == 2:
    50          
3574 46         137 if ($is_two) {
3575             $x = $x->bsqrt($scale+4);
3576             } elsif ($y->is_one('-')) {
3577 0         0 # $x ** -1 => 1/$x
3578             my $u = $class->bone()->bdiv($x, $scale);
3579 0         0 # copy private parts over
3580 0         0 $x->{_m} = $u->{_m};
3581 0         0 $x->{_e} = $u->{_e};
3582             $x->{_es} = $u->{_es};
3583             } else {
3584             # calculate the broot() as integer result first, and if it fits, return
3585             # it rightaway (but only if $x and $y are integer):
3586 20         47  
3587 20 50 33     53 my $done = 0; # not yet
3588 20         80 if ($y->is_int() && $x->is_int()) {
3589 20 50       76 my $i = $LIB->_copy($x->{_m});
3590 20         121 $i = $LIB->_lsft($i, $x->{_e}, 10) unless $LIB->_is_zero($x->{_e});
3591 20         64 my $int = Math::BigInt->bzero();
3592 20         73 $int->{value} = $i;
3593             $int = $int->broot($y->as_number());
3594 20 100       124 # if ($exact)
3595             if ($int->copy()->bpow($y) == $x) {
3596 14         44 # found result, return it
3597 14         58 $x->{_m} = $int->{value};
3598 14         46 $x->{_e} = $LIB->_zero();
3599 14         71 $x->{_es} = '+';
3600 14         53 $x = $x->bnorm();
3601             $done = 1;
3602             }
3603 20 100       81 }
3604 6         67 if ($done == 0) {
3605 6         42 my $u = $class->bone()->bdiv($y, $scale+4);
3606 6         20 $u->{_a} = undef;
3607 6         21 $u->{_p} = undef;
3608             $x = $x->bpow($u, $scale+4); # el cheapo
3609             }
3610 66 50       188 }
3611             $x = $x->bneg() if $sign == 1;
3612              
3613 66 50       146 # shortcut to not run through _find_round_parameters again
3614 66         180 if (defined $params[0]) {
3615             $x = $x->bround($params[0], $params[2]); # then round accordingly
3616 0         0 } else {
3617             $x = $x->bfround($params[1], $params[2]); # then round accordingly
3618 66 50       224 }
3619             if ($fallback) {
3620 66         128 # clear a/p after round, since user did not request it
3621 66         119 $x->{_a} = undef;
3622             $x->{_p} = undef;
3623             }
3624 66         169 # restore globals
3625 66         145 $$abr = $ab;
3626             $$pbr = $pb;
3627 66 0 0     210  
      33        
3628             return $downgrade -> new($x -> bdstr(), @r)
3629 66         964 if defined($downgrade) && ($x -> is_int() || $x -> is_inf());
3630             $x;
3631             }
3632              
3633             sub bfac {
3634             # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
3635             # compute factorial number, modifies first argument
3636              
3637 80 50   80 1 882 # set up parameters
3638             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
3639              
3640 80 50       284 # inf => inf
3641             return $x if $x->modify('bfac');
3642 80 100 100     224  
3643 72 100       237 return $x -> bnan(@r) if $x->is_nan() || $x->is_inf("-");
3644 68 100 100     230 return $x -> binf("+", @r) if $x->is_inf("+");
3645             return $x -> bone(@r) if $x->is_zero() || $x->is_one();
3646 60 100 66     183  
3647 4 50       17 if ($x -> is_neg() || !$x -> is_int()) {
3648 4         22 return $upgrade -> bfac($x, @r) if defined($upgrade);
3649             return $x -> bnan(@r);
3650             }
3651 56 100       173  
3652 8         39 if (! $LIB->_is_zero($x->{_e})) {
3653 8         39 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # change 12e1 to 120e0
3654 8         26 $x->{_e} = $LIB->_zero(); # normalize
3655             $x->{_es} = '+';
3656 56         196 }
3657             $x->{_m} = $LIB->_fac($x->{_m}); # calculate factorial
3658 56         162  
3659             $x = $x->bnorm()->round(@r); # norm again and round result
3660 56 0 0     142  
      33        
3661             return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade)
3662 56         585 && ($x -> is_int() || $x -> is_inf());
3663             $x;
3664             }
3665              
3666             sub bdfac {
3667             # compute double factorial
3668              
3669 72 50   72 1 892 # set up parameters
3670             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
3671 72 50       250  
3672             return $x if $x->modify('bdfac');
3673 72 100 100     198  
3674 64 100       226 return $x -> bnan(@r) if $x->is_nan() || $x->is_inf("-");
3675             return $x -> binf("+", @r) if $x->is_inf("+");
3676 60 100 66     232  
3677 4 50       17 if ($x <= -2 || !$x -> is_int()) {
3678 4         21 return $upgrade -> bdfac($x, @r) if defined($upgrade);
3679             return $x -> bnan(@r);
3680             }
3681 56 100       166  
3682             return $x->bone() if $x <= 1;
3683 44 50       314  
3684             croak("bdfac() requires a newer version of the $LIB library.")
3685             unless $LIB->can('_dfac');
3686 44 100       134  
3687 4         45 if (! $LIB->_is_zero($x->{_e})) {
3688 4         41 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # change 12e1 to 120e0
3689 4         21 $x->{_e} = $LIB->_zero(); # normalize
3690             $x->{_es} = '+';
3691 44         159 }
3692             $x->{_m} = $LIB->_dfac($x->{_m}); # calculate factorial
3693 44         134  
3694             $x = $x->bnorm()->round(@r); # norm again and round result
3695 44 50 33     120  
3696             return $downgrade -> new($x -> bdstr(), @r)
3697 44         443 if defined($downgrade) && $x -> is_int();
3698             return $x;
3699             }
3700              
3701             sub btfac {
3702             # compute triple factorial
3703              
3704 76 50   76 1 822 # set up parameters
3705             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
3706 76 50       251  
3707             return $x if $x->modify('btfac');
3708 76 100 100     220  
3709 68 100       241 return $x -> bnan(@r) if $x->is_nan() || $x->is_inf("-");
3710             return $x -> binf("+", @r) if $x->is_inf("+");
3711 64 100 66     219  
3712 4 50       19 if ($x <= -3 || !$x -> is_int()) {
3713 4         18 return $upgrade -> btfac($x, @r) if defined($upgrade);
3714             return $x -> bnan(@r);
3715             }
3716 60         224  
3717 60 50       244 my $k = $class -> new("3");
3718             return $x->bnan(@r) if $x <= -$k;
3719 60         275  
3720 60 100       144 my $one = $class -> bone();
3721             return $x->bone(@r) if $x <= $one;
3722 44         132  
3723 44         135 my $f = $x -> copy();
3724 60         196 while ($f -> bsub($k) > $one) {
3725             $x = $x -> bmul($f);
3726             }
3727 44         138  
3728             $x = $x->round(@r);
3729 44 50 33     116  
3730             return $downgrade -> new($x -> bdstr(), @r)
3731 44         664 if defined($downgrade) && $x -> is_int();
3732             return $x;
3733             }
3734              
3735 364 50 33 364 1 5324 sub bmfac {
3736             my ($class, $x, $k, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3737             ? (ref($_[0]), @_)
3738             : objectify(2, @_);
3739 364 50       1185  
3740             return $x if $x->modify('bmfac');
3741 364 100 100     987  
      100        
3742 308 100       797 return $x -> bnan(@r) if $x->is_nan() || $x->is_inf("-") || !$k->is_pos();
3743             return $x -> binf("+", @r) if $x->is_inf("+");
3744 288 100 66     956  
      100        
      100        
3745             if ($x <= -$k || !$x -> is_int() ||
3746             ($k -> is_finite() && !$k -> is_int()))
3747 24 50       85 {
3748 24         78 return $upgrade -> bmfac($x, $k, @r) if defined($upgrade);
3749             return $x -> bnan(@r);
3750             }
3751 264         1280  
3752 264 100       633 my $one = $class -> bone();
3753             return $x->bone(@r) if $x <= $one;
3754 184         476  
3755 184         491 my $f = $x -> copy();
3756 284         915 while ($f -> bsub($k) > $one) {
3757             $x = $x -> bmul($f);
3758             }
3759 184         531  
3760             $x = $x->round(@r);
3761 184 50 33     490  
3762             return $downgrade -> new($x -> bdstr(), @r)
3763 184         3093 if defined($downgrade) && $x -> is_int();
3764             return $x;
3765             }
3766              
3767             sub blsft {
3768             # shift left by $y in base $b, i.e., multiply by $b ** $y
3769              
3770 33 50 66 33 1 541 # set up parameters
3771             my ($class, $x, $y, $b, @r)
3772             = ref($_[0]) && ref($_[0]) eq ref($_[1]) && ref($_[1]) eq ref($_[2])
3773             ? (ref($_[0]), @_)
3774             : objectify(2, @_);
3775 33 50       156  
3776             return $x if $x -> modify('blsft');
3777 33 100 66     94  
3778             return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan();
3779 29 50       91  
3780 29 50 33     140 $b = 2 if !defined $b;
3781             $b = $class -> new($b)
3782 29 50       146 unless defined(blessed($b)) && $b -> isa(__PACKAGE__);
3783             return $x -> bnan(@r) if $b -> is_nan();
3784              
3785             # There needs to be more checking for special cases here. Fixme!
3786              
3787 29 50       119 # shift by a negative amount?
3788             return $x -> brsft($y -> copy() -> babs(), $b) if $y -> {sign} =~ /^-/;
3789 29         87  
3790             $x = $x -> bmul($b -> bpow($y), $r[0], $r[1], $r[2], $y);
3791 29 0 0     101  
      33        
3792             return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade)
3793 29         320 && ($x -> is_int() || $x -> is_inf() || $x -> is_nan());
3794             return $x;
3795             }
3796              
3797             sub brsft {
3798             # shift right by $y in base $b, i.e., divide by $b ** $y
3799              
3800 48 50 66 48 1 677 # set up parameters
3801             my ($class, $x, $y, $b, @r)
3802             = ref($_[0]) && ref($_[0]) eq ref($_[1]) && ref($_[1]) eq ref($_[2])
3803             ? (ref($_[0]), @_)
3804             : objectify(2, @_);
3805 48 50       196  
3806             return $x if $x -> modify('brsft');
3807 48 100 66     174  
3808             return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan();
3809              
3810             # There needs to be more checking for special cases here. Fixme!
3811 44 100       131  
3812 44 50 66     262 $b = 2 if !defined $b;
3813             $b = $class -> new($b)
3814 44 50       207 unless defined(blessed($b)) && $b -> isa(__PACKAGE__);
3815             return $x -> bnan(@r) if $b -> is_nan();
3816              
3817 44 50       157 # shift by a negative amount?
3818             return $x -> blsft($y -> copy() -> babs(), $b) if $y -> {sign} =~ /^-/;
3819              
3820 44         167 # call bdiv()
3821             $x = $x -> bdiv($b -> bpow($y), $r[0], $r[1], $r[2], $y);
3822 44 0 0     223  
      33        
3823             return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade)
3824 44         613 && ($x -> is_int() || $x -> is_inf() || $x -> is_nan());
3825             return $x;
3826             }
3827              
3828             ###############################################################################
3829             # Bitwise methods
3830             ###############################################################################
3831              
3832             # Bitwise left shift.
3833              
3834 8 50 33 8 1 72 sub bblsft {
3835             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3836             ? (ref($_[0]), @_)
3837             : objectify(2, @_);
3838 8         47  
3839             my $xint = Math::BigInt -> bblsft($x, $y, @r);
3840              
3841             # disable downgrading
3842 8         89  
3843 8         58 my $dng = $class -> downgrade();
3844             $class -> downgrade(undef);
3845              
3846             # convert to Math::BigFloat without downgrading
3847 8         50  
3848             my $xflt = $class -> new($xint);
3849              
3850             # reset downgrading
3851 8         69  
3852             $class -> downgrade($dng);
3853 8         30  
3854 8         23 $x -> {sign} = $xflt -> {sign};
3855 8         17 $x -> {_m} = $xflt -> {_m};
3856 8         17 $x -> {_es} = $xflt -> {_es};
3857             $x -> {_e} = $xflt -> {_e};
3858              
3859             # now we might downgrade
3860 8 50       26  
3861 8         51 return $downgrade -> new($x) if defined($downgrade);
3862             $x -> round(@r);
3863             }
3864              
3865             # Bitwise left shift.
3866              
3867 8 50 33 8 1 78 sub bbrsft {
3868             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3869             ? (ref($_[0]), @_)
3870             : objectify(2, @_);
3871 8         47  
3872             my $xint = Math::BigInt -> bbrsft($x, $y, @r);
3873              
3874             # disable downgrading
3875 8         36  
3876 8         34 my $dng = $class -> downgrade();
3877             $class -> downgrade(undef);
3878              
3879             # convert to Math::BigFloat without downgrading
3880 8         32  
3881             my $xflt = $class -> new($xint);
3882              
3883             # reset downgrading
3884 8         65  
3885             $class -> downgrade($dng);
3886 8         24  
3887 8         23 $x -> {sign} = $xflt -> {sign};
3888 8         30 $x -> {_m} = $xflt -> {_m};
3889 8         21 $x -> {_es} = $xflt -> {_es};
3890             $x -> {_e} = $xflt -> {_e};
3891              
3892             # now we might downgrade
3893 8 50       31  
3894 8         29 return $downgrade -> new($x) if defined($downgrade);
3895             $x -> round(@r);
3896             }
3897              
3898 1 50 33 1 1 7 sub band {
3899             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3900             ? (ref($_[0]), @_)
3901             : objectify(2, @_);
3902 1 50       7  
3903             return if $x -> modify('band');
3904 1 50 33     6  
3905             return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan();
3906 1         7  
3907 1         6 my $xint = $x -> as_int(); # to Math::BigInt
3908             my $yint = $y -> as_int(); # to Math::BigInt
3909 1         10  
3910             $xint = $xint -> band($yint);
3911 1 50       5  
3912             return $xint -> round(@r) if defined $downgrade;
3913 1         4  
3914 1         3 my $xflt = $class -> new($xint); # back to Math::BigFloat
3915 1         3 $x -> {sign} = $xflt -> {sign};
3916 1         4 $x -> {_m} = $xflt -> {_m};
3917 1         4 $x -> {_es} = $xflt -> {_es};
3918             $x -> {_e} = $xflt -> {_e};
3919 1         3  
3920             return $x -> round(@r);
3921             }
3922              
3923 1 50 33 1 1 21 sub bior {
3924             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3925             ? (ref($_[0]), @_)
3926             : objectify(2, @_);
3927 1 50       7  
3928             return if $x -> modify('bior');
3929 1 50 33     18  
3930             return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan();
3931 1         7  
3932 1         8 my $xint = $x -> as_int(); # to Math::BigInt
3933             my $yint = $y -> as_int(); # to Math::BigInt
3934 1         9  
3935             $xint = $xint -> bior($yint);
3936 1 50       5  
3937             return $xint -> round(@r) if defined $downgrade;
3938 1         15  
3939 1         24 my $xflt = $class -> new($xint); # back to Math::BigFloat
3940 1         5 $x -> {sign} = $xflt -> {sign};
3941 1         4 $x -> {_m} = $xflt -> {_m};
3942 1         2 $x -> {_es} = $xflt -> {_es};
3943             $x -> {_e} = $xflt -> {_e};
3944 1         5  
3945             return $x -> round(@r);
3946             }
3947              
3948 1 50 33 1 1 11 sub bxor {
3949             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3950             ? (ref($_[0]), @_)
3951             : objectify(2, @_);
3952 1 50       7  
3953             return if $x -> modify('bxor');
3954 1 50 33     5  
3955             return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan();
3956 1         8  
3957 1         4 my $xint = $x -> as_int(); # to Math::BigInt
3958             my $yint = $y -> as_int(); # to Math::BigInt
3959 1         8  
3960             $xint = $xint -> bxor($yint);
3961 1 50       6  
3962             return $xint -> round(@r) if defined $downgrade;
3963 1         4  
3964 1         3 my $xflt = $class -> new($xint); # back to Math::BigFloat
3965 1         5 $x -> {sign} = $xflt -> {sign};
3966 1         4 $x -> {_m} = $xflt -> {_m};
3967 1         18 $x -> {_es} = $xflt -> {_es};
3968             $x -> {_e} = $xflt -> {_e};
3969 1         8  
3970             return $x -> round(@r);
3971             }
3972              
3973 4 50   4 1 20 sub bnot {
3974             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
3975 4 50       14  
3976             return if $x -> modify('bnot');
3977 4 50       10  
3978             return $x -> bnan(@r) if $x -> is_nan();
3979 4         16  
3980 4         13 my $xint = $x -> as_int(); # to Math::BigInt
3981             $xint = $xint -> bnot();
3982 4 50       10  
3983             return $xint -> round(@r) if defined $downgrade;
3984 4         11  
3985 4         12 my $xflt = $class -> new($xint); # back to Math::BigFloat
3986 4         9 $x -> {sign} = $xflt -> {sign};
3987 4         8 $x -> {_m} = $xflt -> {_m};
3988 4         6 $x -> {_es} = $xflt -> {_es};
3989             $x -> {_e} = $xflt -> {_e};
3990 4         14  
3991             return $x -> round(@r);
3992             }
3993              
3994             ###############################################################################
3995             # Rounding methods
3996             ###############################################################################
3997              
3998             sub bround {
3999             # accuracy: preserve $N digits, and overwrite the rest with 0's
4000 41467 50   41467 1 136369  
4001             my ($class, $x, @a) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4002 41467 100 100     109102  
4003 3         813 if (($a[0] || 0) < 0) {
4004             croak('bround() needs positive accuracy');
4005             }
4006 41464 50       113751  
4007             return $x if $x->modify('bround');
4008 41464         109535  
4009 41464 100       90488 my ($scale, $mode) = $x->_scale_a(@a);
4010 3889 50 66     7759 if (!defined $scale) { # no-op
      66        
4011             return $downgrade -> new($x) if defined($downgrade)
4012 3885         10920 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4013             return $x;
4014             }
4015              
4016             # Scale is now either $x->{_a}, $accuracy, or the input argument. Test
4017             # whether $x already has lower accuracy, do nothing in this case but do
4018             # round if the accuracy is the same, since a math operation might want to
4019             # round a number with A=5 to 5 digits afterwards again
4020 37575 100 100     120932  
4021 4 0 0     17 if (defined $x->{_a} && $x->{_a} < $scale) {
      33        
4022             return $downgrade -> new($x) if defined($downgrade)
4023 4         16 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4024             return $x;
4025             }
4026              
4027             # scale < 0 makes no sense
4028             # scale == 0 => keep all digits
4029             # never round a +-inf, NaN
4030 37571 100 66     163091  
4031 12 0 0     41 if ($scale <= 0 || $x->{sign} !~ /^[+-]$/) {
      33        
4032             return $downgrade -> new($x) if defined($downgrade)
4033 12         126 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4034             return $x;
4035             }
4036              
4037             # 1: never round a 0
4038 37559 100 100     85286 # 2: if we should keep more digits than the mantissa has, do nothing
4039 12805 100 100     45428 if ($x->is_zero() || $LIB->_len($x->{_m}) <= $scale) {
4040 12805 50 33     24881 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
      66        
4041             return $downgrade -> new($x) if defined($downgrade)
4042 12805         38685 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4043             return $x;
4044             }
4045              
4046 24754         96776 # pass sign to bround for '+inf' and '-inf' rounding modes
4047             my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
4048 24754         71523  
4049 24754         53240 $m = $m->bround($scale, $mode); # round mantissa
4050 24754         60138 $x->{_m} = $m->{value}; # get our mantissa back
4051 24754         38432 $x->{_a} = $scale; # remember rounding
4052             $x->{_p} = undef; # and clear P
4053              
4054 24754         61852 # bnorm() downgrades if necessary, so no need to check whether to downgrade.
4055             $x->bnorm(); # del trailing zeros gen. by bround()
4056             }
4057              
4058             sub bfround {
4059             # precision: round to the $Nth digit left (+$n) or right (-$n) from the '.'
4060             # $n == 0 means round to integer
4061             # expects and returns normalized numbers!
4062 724 50   724 1 9912  
4063             my ($class, $x, @p) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4064 724 50       2303  
4065             return $x if $x->modify('bfround'); # no-op
4066 724         2190  
4067 724 100       1641 my ($scale, $mode) = $x->_scale_p(@p);
4068 4 50 66     15 if (!defined $scale) {
      33        
4069             return $downgrade -> new($x) if defined($downgrade)
4070 0         0 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4071             return $x;
4072             }
4073              
4074             # never round a 0, +-inf, NaN
4075 720 100       1721  
4076 20 50 33     177 if ($x->is_zero()) {
4077 20 0 0     76 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
      33        
4078             return $downgrade -> new($x) if defined($downgrade)
4079 20         103 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4080             return $x;
4081             }
4082 700 100       3014  
4083 12 0 0     40 if ($x->{sign} !~ /^[+-]$/) {
      33        
4084             return $downgrade -> new($x) if defined($downgrade)
4085 12         128 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4086             return $x;
4087             }
4088              
4089 688 50 100     1811 # don't round if x already has lower precision
      66        
4090 0 0 0     0 if (defined $x->{_p} && $x->{_p} < 0 && $scale < $x->{_p}) {
      0        
4091             return $downgrade -> new($x) if defined($downgrade)
4092 0         0 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4093             return $x;
4094             }
4095 688         1378  
4096 688         1203 $x->{_p} = $scale; # remember round in any case
4097 688 100       1442 $x->{_a} = undef; # and clear A
4098             if ($scale < 0) {
4099             # round right from the '.'
4100 532 100       1193  
4101 28 0 0     85 if ($x->{_es} eq '+') { # e >= 0 => nothing to round
      33        
4102             return $downgrade -> new($x) if defined($downgrade)
4103 28         91 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4104             return $x;
4105             }
4106 504         850  
4107 504         1472 $scale = -$scale; # positive for simplicity
4108             my $len = $LIB->_len($x->{_m}); # length of mantissa
4109              
4110             # the following poses a restriction on _e, but if _e is bigger than a
4111 504         1551 # scalar, you got other problems (memory etc) anyway
4112 504         993 my $dad = -(0+ ($x->{_es}.$LIB->_num($x->{_e}))); # digits after dot
4113 504 100       1050 my $zad = 0; # zeros after dot
4114             $zad = $dad - $len if (-$dad < -$len); # for 0.00..00xxx style
4115              
4116             # print "scale $scale dad $dad zad $zad len $len\n";
4117             # number bsstr len zad dad
4118             # 0.123 123e-3 3 0 3
4119             # 0.0123 123e-4 3 1 4
4120             # 0.001 1e-3 1 2 3
4121             # 1.23 123e-2 3 0 2
4122             # 1.2345 12345e-4 5 0 4
4123              
4124             # do not round after/right of the $dad
4125 504 100       998  
4126 47 0 0     122 if ($scale > $dad) { # 0.123, scale >= 3 => exit
      33        
4127             return $downgrade -> new($x) if defined($downgrade)
4128 47         443 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4129             return $x;
4130             }
4131              
4132             # round to zero if rounding inside the $zad, but not for last zero like:
4133             # 0.0065, scale -2, round last '0' with following '65' (scale == zad
4134 457 100       1013 # case)
4135 40 0 0     128 if ($scale < $zad) {
      33        
4136             return $downgrade -> new($x) if defined($downgrade)
4137 40         148 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4138             return $x->bzero();
4139             }
4140 417 100       856  
4141 41         92 if ($scale == $zad) { # for 0.006, scale -3 and trunc
4142             $scale = -$len;
4143             } else {
4144 376 100       698 # adjust round-point to be inside mantissa
4145 78         150 if ($zad != 0) {
4146             $scale = $scale-$zad;
4147 298         433 } else {
4148 298 50       616 my $dbd = $len - $dad;
4149 298         534 $dbd = 0 if $dbd < 0; # digits before dot
4150             $scale = $dbd+$scale;
4151             }
4152             }
4153             } else {
4154             # round left from the '.'
4155              
4156             # 123 => 100 means length(123) = 3 - $scale (2) => 1
4157 156         496  
4158             my $dbt = $LIB->_len($x->{_m});
4159 156         548 # digits before dot
4160             my $dbd = $dbt + ($x->{_es} . $LIB->_num($x->{_e}));
4161 156 100       457 # should be the same, so treat it as this
4162             $scale = 1 if $scale == 0;
4163 156 100 100     576 # shortcut if already integer
4164 10 0 0     39 if ($scale == 1 && $dbt <= $dbd) {
      33        
4165             return $downgrade -> new($x) if defined($downgrade)
4166 10         36 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4167             return $x;
4168             }
4169 146         214 # maximum digits before dot
4170             ++$dbd;
4171 146 100       398  
    100          
4172             if ($scale > $dbd) {
4173 27 50       80 # not enough digits before dot, so round to zero
4174 27         96 return $downgrade -> new($x) if defined($downgrade);
4175             return $x->bzero;
4176             } elsif ($scale == $dbd) {
4177 72         140 # maximum
4178             $scale = -$dbt;
4179 47         87 } else {
4180             $scale = $dbd - $scale;
4181             }
4182             }
4183              
4184 536         2002 # pass sign to bround for rounding modes '+inf' and '-inf'
4185 536         1701 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
4186 536         1320 $m = $m->bround($scale, $mode);
4187             $x->{_m} = $m->{value}; # get our mantissa back
4188              
4189 536         1391 # bnorm() downgrades if necessary, so no need to check whether to downgrade.
4190             $x->bnorm();
4191             }
4192              
4193             sub bfloor {
4194 84 50   84 1 872 # round towards minus infinity
4195             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4196 84 50       322  
4197             return $x if $x->modify('bfloor');
4198 84 100       278  
4199             return $x -> bnan(@r) if $x -> is_nan();
4200 79 100       344  
4201             if ($x->{sign} =~ /^[+-]$/) {
4202 70 100       208 # if $x has digits after dot, remove them
4203 47         229 if ($x->{_es} eq '-') {
4204 47         173 $x->{_m} = $LIB->_rsft($x->{_m}, $x->{_e}, 10);
4205 47         152 $x->{_e} = $LIB->_zero();
4206             $x->{_es} = '+';
4207 47 100       183 # increment if negative
4208             $x->{_m} = $LIB->_inc($x->{_m}) if $x->{sign} eq '-';
4209 70         269 }
4210             $x = $x->round(@r);
4211 79 100       256 }
4212 77         506 return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade);
4213             return $x;
4214             }
4215              
4216             sub bceil {
4217 43 50   43 1 570 # round towards plus infinity
4218             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4219 43 50       163  
4220             return $x if $x->modify('bceil');
4221 43 100       145  
4222             return $x -> bnan(@r) if $x -> is_nan();
4223              
4224 38 100       192 # if $x has digits after dot, remove them
4225 29 100       98 if ($x->{sign} =~ /^[+-]$/) {
4226 17         170 if ($x->{_es} eq '-') {
4227 17         83 $x->{_m} = $LIB->_rsft($x->{_m}, $x->{_e}, 10);
4228 17         58 $x->{_e} = $LIB->_zero();
4229 17 100       69 $x->{_es} = '+';
4230 9         84 if ($x->{sign} eq '+') {
4231             $x->{_m} = $LIB->_inc($x->{_m}); # increment if positive
4232 8 100       31 } else {
4233             $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # avoid -0
4234             }
4235 29         112 }
4236             $x = $x->round(@r);
4237             }
4238 38 100       152  
4239 36         323 return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade);
4240             return $x;
4241             }
4242              
4243             sub bint {
4244 179 50   179 1 1021 # round towards zero
4245             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4246 179 50       582  
4247             return $x if $x->modify('bint');
4248 179 100       518  
4249             return $x -> bnan(@r) if $x -> is_nan();
4250 174 100       771  
4251             if ($x->{sign} =~ /^[+-]$/) {
4252 165 100       427 # if $x has digits after the decimal point
4253 13         103 if ($x->{_es} eq '-') {
4254 13         74 $x->{_m} = $LIB->_rsft($x->{_m}, $x->{_e}, 10); # remove frac part
4255 13         35 $x->{_e} = $LIB->_zero(); # truncate/normalize
4256 13 100       48 $x->{_es} = '+'; # abs e
4257             $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # avoid -0
4258 165         455 }
4259             $x = $x->round(@r);
4260             }
4261 174 100       500  
4262 172         781 return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade);
4263             return $x;
4264             }
4265              
4266             ###############################################################################
4267             # Other mathematical methods
4268             ###############################################################################
4269              
4270             sub bgcd {
4271             # (BINT or num_str, BINT or num_str) return BINT
4272             # does not modify arguments, but returns new object
4273              
4274 133 100 100 133 1 2409 # Class::method(...) -> Class->method(...)
      66        
4275             unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) ||
4276             $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i))
4277             {
4278             #carp "Using ", (caller(0))[3], "() as a function is deprecated;",
4279 1         5 # " use is as a method instead";
4280             unshift @_, __PACKAGE__;
4281             }
4282 133         469  
4283             my ($class, @args) = objectify(0, @_);
4284 133         265  
4285 133 50 33     538 my $x = shift @args;
4286             $x = defined(blessed($x)) && $x -> isa(__PACKAGE__) ? $x -> copy()
4287 133 100       377 : $class -> new($x);
4288             return $class->bnan() unless $x -> is_int();
4289 105         275  
4290 116         209 while (@args) {
4291 116 50 33     420 my $y = shift @args;
4292             $y = $class->new($y)
4293 116 100       264 unless defined(blessed($y)) && $y -> isa(__PACKAGE__);
4294             return $class->bnan() unless $y -> is_int();
4295              
4296 104         282 # greatest common divisor
4297 246         561 while (! $y->is_zero()) {
4298             ($x, $y) = ($y->copy(), $x->copy()->bmod($y));
4299             }
4300 104 100       237  
4301             last if $x -> is_one();
4302 93         294 }
4303             $x = $x -> babs();
4304 93 50 33     238  
4305             return $downgrade -> new($x)
4306 93         934 if defined $downgrade && $x->is_int();
4307             return $x;
4308             }
4309              
4310             sub blcm {
4311             # (BFLOAT or num_str, BFLOAT or num_str) return BFLOAT
4312             # does not modify arguments, but returns new object
4313             # Least Common Multiple
4314              
4315 35 100 100 35 1 1202 # Class::method(...) -> Class->method(...)
      66        
4316             unless (@_ && (defined(blessed($_[0])) && $_[0] -> isa(__PACKAGE__) ||
4317             $_[0] =~ /^[a-z]\w*(?:::[a-z]\w*)*$/i))
4318             {
4319             #carp "Using ", (caller(0))[3], "() as a function is deprecated;",
4320 1         4 # " use is as a method instead";
4321             unshift @_, __PACKAGE__;
4322             }
4323 35         140  
4324             my ($class, @args) = objectify(0, @_);
4325 35         139  
4326 35 50 33     147 my $x = shift @args;
4327             $x = defined(blessed($x)) && $x -> isa(__PACKAGE__) ? $x -> copy()
4328 35 100       173 : $class -> new($x);
4329             return $class->bnan() if $x->{sign} !~ /^[+-]$/; # x NaN?
4330 27         123  
4331 30         60 while (@args) {
4332 30 50 33     130 my $y = shift @args;
4333             $y = $class -> new($y)
4334 30 100       87 unless defined(blessed($y)) && $y -> isa(__PACKAGE__);
4335 26         81 return $x->bnan() unless $y -> is_int();
4336 26         177 my $gcd = $x -> bgcd($y);
4337             $x = $x -> bdiv($gcd) -> bmul($y);
4338             }
4339 23         75  
4340             $x = $x -> babs();
4341 23 50 33     82  
4342             return $downgrade -> new($x)
4343 23         293 if defined $downgrade && $x->is_int();
4344             return $x;
4345             }
4346              
4347             ###############################################################################
4348             # Object property methods
4349             ###############################################################################
4350              
4351 24 50   24 1 390 sub length {
4352             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4353 24 50       77  
4354             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4355 24 100       110  
4356             return 1 if $LIB->_is_zero($x->{_m});
4357 20         97  
4358 20 50       116 my $len = $LIB->_len($x->{_m});
4359 20 50       62 $len += $LIB->_num($x->{_e}) if $x->{_es} eq '+';
4360 0         0 if (wantarray()) {
4361 0 0       0 my $t = 0;
4362 0         0 $t = $LIB->_num($x->{_e}) if $x->{_es} eq '-';
4363             return ($len, $t);
4364 20         159 }
4365             $len;
4366             }
4367              
4368             sub mantissa {
4369 36 50   36 1 514 # return a copy of the mantissa
4370             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4371              
4372             # The following line causes a lot of noise in the test suits for
4373             # the Math-BigRat and bignum distributions. Fixme!
4374             #carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4375 36 100       108  
4376             return $x -> bnan(@r) if $x -> is_nan();
4377 32 100       143  
4378 8         26 if ($x->{sign} !~ /^[+-]$/) {
4379 8         36 my $s = $x->{sign};
4380 8         35 $s =~ s/^\+//;
4381             return Math::BigInt->new($s, undef, undef); # -inf, +inf => +inf
4382 24         106 }
4383 24 100       143 my $m = Math::BigInt->new($LIB->_str($x->{_m}), undef, undef);
4384 24         133 $m = $m->bneg() if $x->{sign} eq '-';
4385             $m;
4386             }
4387              
4388             sub exponent {
4389 36 50   36 1 490 # return a copy of the exponent
4390             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4391              
4392             # The following line causes a lot of noise in the test suits for
4393             # the Math-BigRat and bignum distributions. Fixme!
4394             #carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4395 36 100       117  
4396             return $x -> bnan(@r) if $x -> is_nan();
4397 32 100       137  
4398 8         51 if ($x->{sign} !~ /^[+-]$/) {
4399 8         45 my $s = $x->{sign};
4400 8         37 $s =~ s/^[+-]//;
4401             return Math::BigInt->new($s, undef, undef); # -inf, +inf => +inf
4402 24         121 }
4403             Math::BigInt->new($x->{_es} . $LIB->_str($x->{_e}), undef, undef);
4404             }
4405              
4406             sub parts {
4407 32 50   32 1 465 # return a copy of both the exponent and the mantissa
4408             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4409 32 50       87  
4410             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4411 32 100       142  
4412 12         45 if ($x->{sign} !~ /^[+-]$/) {
4413 12         48 my $s = $x->{sign};
4414 12         33 $s =~ s/^\+//;
4415 12         37 my $se = $s;
4416             $se =~ s/^-//;
4417 12         74 # +inf => inf and -inf, +inf => inf
4418             return ($class->new($s), $class->new($se));
4419 20         80 }
4420 20         84 my $m = Math::BigInt->bzero();
4421 20 100       100 $m->{value} = $LIB->_copy($x->{_m});
4422 20         85 $m = $m->bneg() if $x->{sign} eq '-';
4423             ($m, Math::BigInt->new($x->{_es} . $LIB->_num($x->{_e})));
4424             }
4425              
4426             # Parts used for scientific notation with significand/mantissa and exponent as
4427             # integers. E.g., "12345.6789" is returned as "123456789" (mantissa) and "-4"
4428             # (exponent).
4429              
4430 0 0   0 1 0 sub sparts {
4431             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4432 0 0       0  
4433             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4434              
4435             # Not-a-number.
4436 0 0       0  
4437 0         0 if ($x -> is_nan()) {
4438 0 0       0 my $mant = $class -> bnan(); # mantissa
4439 0         0 return $mant unless wantarray; # scalar context
4440 0         0 my $expo = $class -> bnan(); # exponent
4441             return ($mant, $expo); # list context
4442             }
4443              
4444             # Infinity.
4445 0 0       0  
4446 0         0 if ($x -> is_inf()) {
4447 0 0       0 my $mant = $class -> binf($x->{sign}); # mantissa
4448 0         0 return $mant unless wantarray; # scalar context
4449 0         0 my $expo = $class -> binf('+'); # exponent
4450             return ($mant, $expo); # list context
4451             }
4452              
4453             # Finite number.
4454 0         0  
4455 0         0 my $mant = $class -> new($x);
4456 0         0 $mant->{_es} = '+';
4457 0 0       0 $mant->{_e} = $LIB->_zero();
4458 0 0       0 $mant = $downgrade -> new($mant) if defined $downgrade;
4459             return $mant unless wantarray;
4460 0         0  
4461 0 0       0 my $expo = $class -> new($x -> {_es} . $LIB->_str($x -> {_e}));
4462 0         0 $expo = $downgrade -> new($expo) if defined $downgrade;
4463             return ($mant, $expo);
4464             }
4465              
4466             # Parts used for normalized notation with significand/mantissa as either 0 or a
4467             # number in the semi-open interval [1,10). E.g., "12345.6789" is returned as
4468             # "1.23456789" and "4".
4469              
4470 0 0   0 1 0 sub nparts {
4471             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4472 0 0       0  
4473             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4474              
4475             # Not-a-number and Infinity.
4476 0 0 0     0  
4477             return $x -> sparts() if $x -> is_nan() || $x -> is_inf();
4478              
4479             # Finite number.
4480 0         0  
4481             my ($mant, $expo) = $x -> sparts();
4482 0 0       0  
4483 0         0 if ($mant -> bcmp(0)) {
4484 0         0 my ($ndigtot, $ndigfrac) = $mant -> length();
4485             my $expo10adj = $ndigtot - $ndigfrac - 1;
4486 0 0       0  
4487 0         0 if ($expo10adj > 0) { # if mantissa is not an integer
4488 0 0       0 $mant = $mant -> brsft($expo10adj, 10);
4489 0         0 return $mant unless wantarray;
4490 0         0 $expo = $expo -> badd($expo10adj);
4491             return ($mant, $expo);
4492             }
4493             }
4494 0 0       0  
4495 0         0 return $mant unless wantarray;
4496             return ($mant, $expo);
4497             }
4498              
4499             # Parts used for engineering notation with significand/mantissa as either 0 or a
4500             # number in the semi-open interval [1,1000) and the exponent is a multiple of 3.
4501             # E.g., "12345.6789" is returned as "12.3456789" and "3".
4502              
4503 0 0   0 1 0 sub eparts {
4504             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4505 0 0       0  
4506             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4507              
4508             # Not-a-number and Infinity.
4509 0 0 0     0  
4510             return $x -> sparts() if $x -> is_nan() || $x -> is_inf();
4511              
4512             # Finite number.
4513 0         0  
4514             my ($mant, $expo) = $x -> nparts();
4515 0         0  
4516 0         0 my $c = $expo -> copy() -> bmod(3);
4517 0 0       0 $mant = $mant -> blsft($c, 10);
4518             return $mant unless wantarray;
4519 0         0  
4520 0         0 $expo = $expo -> bsub($c);
4521             return ($mant, $expo);
4522             }
4523              
4524             # Parts used for decimal notation, e.g., "12345.6789" is returned as "12345"
4525             # (integer part) and "0.6789" (fraction part).
4526              
4527 0 0   0 1 0 sub dparts {
4528             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4529 0 0       0  
4530             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4531              
4532             # Not-a-number.
4533 0 0       0  
4534 0         0 if ($x -> is_nan()) {
4535 0 0       0 my $int = $class -> bnan();
4536 0         0 return $int unless wantarray;
4537 0         0 my $frc = $class -> bzero(); # or NaN?
4538             return ($int, $frc);
4539             }
4540              
4541             # Infinity.
4542 0 0       0  
4543 0         0 if ($x -> is_inf()) {
4544 0 0       0 my $int = $class -> binf($x->{sign});
4545 0         0 return $int unless wantarray;
4546 0         0 my $frc = $class -> bzero();
4547             return ($int, $frc);
4548             }
4549              
4550             # Finite number.
4551 0         0  
4552 0         0 my $int = $x -> copy();
4553             my $frc;
4554              
4555             # If the input is an integer.
4556 0 0       0  
4557 0         0 if ($int->{_es} eq '+') {
4558             $frc = $class -> bzero();
4559             }
4560              
4561             # If the input has a fraction part
4562              
4563 0         0 else {
4564 0         0 $int->{_m} = $LIB -> _rsft($int->{_m}, $int->{_e}, 10);
4565 0         0 $int->{_e} = $LIB -> _zero();
4566 0 0       0 $int->{_es} = '+';
4567 0 0       0 $int->{sign} = '+' if $LIB->_is_zero($int->{_m}); # avoid -0
4568 0         0 return $int unless wantarray;
4569 0         0 $frc = $x -> copy() -> bsub($int);
4570             return ($int, $frc);
4571             }
4572 0 0       0  
4573 0 0       0 $int = $downgrade -> new($int) if defined $downgrade;
4574 0         0 return $int unless wantarray;
4575             return $int, $frc;
4576             }
4577              
4578             # Fractional parts with the numerator and denominator as integers. E.g.,
4579             # "123.4375" is returned as "1975" and "16".
4580              
4581 0 0   0 1 0 sub fparts {
4582             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4583 0 0       0  
4584             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4585              
4586             # NaN => NaN/NaN
4587 0 0       0  
4588 0 0       0 if ($x -> is_nan()) {
4589 0         0 return $class -> bnan() unless wantarray;
4590             return $class -> bnan(), $class -> bnan();
4591             }
4592              
4593             # ±Inf => ±Inf/1
4594 0 0       0  
4595 0         0 if ($x -> is_inf()) {
4596 0 0       0 my $numer = $class -> binf($x->{sign});
4597 0         0 return $numer unless wantarray;
4598 0         0 my $denom = $class -> bone();
4599             return $numer, $denom;
4600             }
4601              
4602             # Finite number.
4603              
4604             # If we get here, we know that the output is an integer.
4605 0 0       0  
4606             $class = $downgrade if defined $downgrade;
4607 0         0  
4608 0         0 my @flt_parts = ($x->{sign}, $x->{_m}, $x->{_es}, $x->{_e});
4609 0         0 my @rat_parts = $class -> _flt_lib_parts_to_rat_lib_parts(@flt_parts);
4610 0         0 my $num = $class -> new($LIB -> _str($rat_parts[1]));
4611 0 0       0 my $den = $class -> new($LIB -> _str($rat_parts[2]));
4612 0 0       0 $num = $num -> bneg() if $rat_parts[0] eq "-";
4613 0         0 return $num unless wantarray;
4614             return $num, $den;
4615             }
4616              
4617             # Given "123.4375", returns "1975", since "123.4375" is "1975/16".
4618              
4619 0 0   0 1 0 sub numerator {
4620             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4621 0 0       0  
4622             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4623 0 0       0  
4624 0 0       0 return $class -> bnan() if $x -> is_nan();
4625 0 0       0 return $class -> binf($x -> sign()) if $x -> is_inf();
4626             return $class -> bzero() if $x -> is_zero();
4627              
4628             # If we get here, we know that the output is an integer.
4629 0 0       0  
4630             $class = $downgrade if defined $downgrade;
4631 0 0       0  
    0          
4632 0         0 if ($x -> {_es} eq '-') { # exponent < 0
4633 0         0 my $numer_lib = $LIB -> _copy($x -> {_m});
4634 0         0 my $denom_lib = $LIB -> _1ex($x -> {_e});
4635 0         0 my $gcd_lib = $LIB -> _gcd($LIB -> _copy($numer_lib), $denom_lib);
4636 0         0 $numer_lib = $LIB -> _div($numer_lib, $gcd_lib);
4637             return $class -> new($x -> {sign} . $LIB -> _str($numer_lib));
4638             }
4639              
4640 0         0 elsif (! $LIB -> _is_zero($x -> {_e})) { # exponent > 0
4641 0         0 my $numer_lib = $LIB -> _copy($x -> {_m});
4642 0         0 $numer_lib = $LIB -> _lsft($numer_lib, $x -> {_e}, 10);
4643             return $class -> new($x -> {sign} . $LIB -> _str($numer_lib));
4644             }
4645              
4646 0         0 else { # exponent = 0
4647             return $class -> new($x -> {sign} . $LIB -> _str($x -> {_m}));
4648             }
4649             }
4650              
4651             # Given "123.4375", returns "16", since "123.4375" is "1975/16".
4652              
4653 0 0   0 1 0 sub denominator {
4654             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4655 0 0       0  
4656             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4657 0 0       0  
4658             return $class -> bnan() if $x -> is_nan();
4659              
4660             # If we get here, we know that the output is an integer.
4661 0 0       0  
4662             $class = $downgrade if defined $downgrade;
4663 0 0       0  
4664 0         0 if ($x -> {_es} eq '-') { # exponent < 0
4665 0         0 my $numer_lib = $LIB -> _copy($x -> {_m});
4666 0         0 my $denom_lib = $LIB -> _1ex($x -> {_e});
4667 0         0 my $gcd_lib = $LIB -> _gcd($LIB -> _copy($numer_lib), $denom_lib);
4668 0         0 $denom_lib = $LIB -> _div($denom_lib, $gcd_lib);
4669             return $class -> new($LIB -> _str($denom_lib));
4670             }
4671              
4672 0         0 else { # exponent >= 0
4673             return $class -> bone();
4674             }
4675             }
4676              
4677             ###############################################################################
4678             # String conversion methods
4679             ###############################################################################
4680              
4681             sub bstr {
4682             # (ref to BFLOAT or num_str) return num_str
4683             # Convert number from internal format to (non-scientific) string format.
4684 8480 100   8480 1 1837771 # internal format is always normalized (no leading zeros, "-0" => "+0")
4685             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4686 8480 50       23666  
4687             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4688              
4689             # Inf and NaN
4690 8480 100 100     32891  
4691 2477 100       22071 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
4692 549         5385 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
4693             return 'inf'; # +inf
4694             }
4695              
4696             # Finite number
4697 6003         11732  
4698 6003         9063 my $es = '0';
4699 6003         8758 my $len = 1;
4700 6003         8490 my $cad = 0;
4701             my $dot = '.';
4702              
4703 6003   100     27318 # $x is zero?
4704 6003 100       13485 my $not_zero = !($x->{sign} eq '+' && $LIB->_is_zero($x->{_m}));
4705 4725         14657 if ($not_zero) {
4706 4725         9272 $es = $LIB->_str($x->{_m});
4707 4725         12599 $len = CORE::length($es);
4708 4725 100       12078 my $e = $LIB->_num($x->{_e});
4709 4725 100       12940 $e = -$e if $x->{_es} eq '-';
    100          
4710 1586         3157 if ($e < 0) {
4711             $dot = '';
4712 1586 100       3810 # if _e is bigger than a scalar, the following will blow your memory
4713 576         1258 if ($e <= -$len) {
4714 576         1667 my $r = abs($e) - $len;
4715 576         1307 $es = '0.'. ('0' x $r) . $es;
4716             $cad = -($len+$r);
4717 1010         2538 } else {
4718 1010         2845 substr($es, $e, 0) = '.';
4719 1010 50       3047 $cad = $LIB->_num($x->{_e});
4720             $cad = -$cad if $x->{_es} eq '-';
4721             }
4722             } elsif ($e > 0) {
4723 678         1891 # expand with zeros
4724 678         1240 $es .= '0' x $e;
4725 678         1225 $len += $e;
4726             $cad = 0;
4727             }
4728             } # if not zero
4729 6003 100       14797  
4730             $es = '-'.$es if $x->{sign} eq '-';
4731 6003 100 100     27722 # if set accuracy or precision, pad with zeros on the right side
    100 100        
4732             if ((defined $x->{_a}) && ($not_zero)) {
4733 694         1289 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
4734 694 50       1591 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
4735 694 100       1587 $zeros = $x->{_a} - $len if $cad != $len;
4736             $es .= $dot.'0' x $zeros if $zeros > 0;
4737             } elsif ((($x->{_p} || 0) < 0)) {
4738 502         899 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
4739 502 100       1148 my $zeros = -$x->{_p} + $cad;
4740             $es .= $dot.'0' x $zeros if $zeros > 0;
4741 6003         80277 }
4742             $es;
4743             }
4744              
4745             # Decimal notation, e.g., "12345.6789" (no exponent).
4746              
4747 48 50   48 1 178 sub bdstr {
4748             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4749 48 50       106  
4750             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4751              
4752             # Inf and NaN
4753 48 100 100     159  
4754 9 100       31 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
4755 7         32 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
4756             return 'inf'; # +inf
4757             }
4758              
4759             # Upgrade?
4760 39 50 33     100  
4761             return $upgrade -> bdstr($x, @r)
4762             if defined($upgrade) && !$x -> isa(__PACKAGE__);
4763              
4764             # Finite number
4765 39         117  
4766 39         108 my $mant = $LIB->_str($x->{_m});
4767 39         100 my $esgn = $x->{_es};
4768             my $eabs = $LIB -> _num($x->{_e});
4769 39         64  
4770             my $uintmax = ~0;
4771 39         67  
4772 39 50       82 my $str = $mant;
4773             if ($esgn eq '+') {
4774 39 50       89  
4775             croak("The absolute value of the exponent is too large")
4776             if $eabs > $uintmax;
4777 39         87  
4778             $str .= "0" x $eabs;
4779              
4780 0         0 } else {
4781 0         0 my $mlen = CORE::length($mant);
4782             my $c = $mlen - $eabs;
4783 0         0  
4784 0 0       0 my $intmax = ($uintmax - 1) / 2;
4785             croak("The absolute value of the exponent is too large")
4786             if (1 - $c) > $intmax;
4787 0 0       0  
4788 0         0 $str = "0" x (1 - $c) . $str if $c <= 0;
4789             substr($str, -$eabs, 0) = '.';
4790             }
4791 39 100       218  
4792             return $x->{sign} eq '-' ? '-' . $str : $str;
4793             }
4794              
4795             # Scientific notation with significand/mantissa and exponent as integers, e.g.,
4796             # "12345.6789" is written as "123456789e-4".
4797              
4798 175 100   175 1 12912 sub bsstr {
4799             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4800 175 50       415  
4801             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4802              
4803             # Inf and NaN
4804 175 100 100     654  
4805 28 100       231 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
4806 12         107 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
4807             return 'inf'; # +inf
4808             }
4809              
4810             # Upgrade?
4811 147 50 33     410  
4812             return $upgrade -> bsstr($x, @r)
4813             if defined($upgrade) && !$x -> isa(__PACKAGE__);
4814              
4815             # Finite number
4816              
4817 147 100       559 ($x->{sign} eq '-' ? '-' : '') . $LIB->_str($x->{_m})
4818             . 'e' . $x->{_es} . $LIB->_str($x->{_e});
4819             }
4820              
4821             # Normalized notation, e.g., "12345.6789" is written as "1.23456789e+4".
4822              
4823 483 50   483 1 1326 sub bnstr {
4824             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4825 483 50       991  
4826             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4827              
4828             # Inf and NaN
4829 483 50 66     1353  
4830 0 0       0 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
4831 0         0 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
4832             return 'inf'; # +inf
4833             }
4834              
4835             # Upgrade?
4836 483 50 33     1088  
4837             return $upgrade -> bnstr($x, @r)
4838             if defined($upgrade) && !$x -> isa(__PACKAGE__);
4839              
4840             # Finite number
4841 483 100       1122  
4842             my $str = $x->{sign} eq '-' ? '-' : '';
4843              
4844             # Get the mantissa and the length of the mantissa.
4845 483         1466  
4846 483         918 my $mant = $LIB->_str($x->{_m});
4847             my $mantlen = CORE::length($mant);
4848 483 100       1021  
4849             if ($mantlen == 1) {
4850              
4851             # Not decimal point when the mantissa has length one, i.e., return the
4852             # number 2 as the string "2", not "2.".
4853 65         288  
4854             $str .= $mant . 'e' . $x->{_es} . $LIB->_str($x->{_e});
4855              
4856             } else {
4857              
4858             # Compute new exponent where the original exponent is adjusted by the
4859             # length of the mantissa minus one (because the decimal point is after
4860             # one digit).
4861              
4862 418         1406 my ($eabs, $esgn) = $LIB -> _sadd($LIB -> _copy($x->{_e}), $x->{_es},
4863 418         1390 $LIB -> _new($mantlen - 1), "+");
4864 418         1429 substr $mant, 1, 0, ".";
4865             $str .= $mant . 'e' . $esgn . $LIB->_str($eabs);
4866              
4867             }
4868 483         2829  
4869             return $str;
4870             }
4871              
4872             # Engineering notation, e.g., "12345.6789" is written as "12.3456789e+3".
4873              
4874 0 0   0 1 0 sub bestr {
4875             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4876 0 0       0  
4877             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4878              
4879             # Inf and NaN
4880 0 0 0     0  
4881 0 0       0 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
4882 0         0 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
4883             return 'inf'; # +inf
4884             }
4885              
4886             # Upgrade?
4887 0 0 0     0  
4888             return $upgrade -> bestr($x, @r)
4889             if defined($upgrade) && !$x -> isa(__PACKAGE__);
4890              
4891             # Finite number
4892 0 0       0  
4893             my $str = $x->{sign} eq '-' ? '-' : '';
4894              
4895             # Get the mantissa, the length of the mantissa, and adjust the exponent by
4896             # the length of the mantissa minus 1 (because the dot is after one digit).
4897 0         0  
4898 0         0 my $mant = $LIB->_str($x->{_m});
4899             my $mantlen = CORE::length($mant);
4900 0         0 my ($eabs, $esgn) = $LIB -> _sadd($LIB -> _copy($x->{_e}), $x->{_es},
4901             $LIB -> _new($mantlen - 1), "+");
4902 0         0  
4903 0         0 my $dotpos = 1;
4904 0 0       0 my $mod = $LIB -> _mod($LIB -> _copy($eabs), $LIB -> _new("3"));
4905 0 0       0 unless ($LIB -> _is_zero($mod)) {
4906 0         0 if ($esgn eq '+') {
4907 0         0 $eabs = $LIB -> _sub($eabs, $mod);
4908             $dotpos += $LIB -> _num($mod);
4909 0         0 } else {
4910 0         0 my $delta = $LIB -> _sub($LIB -> _new("3"), $mod);
4911 0         0 $eabs = $LIB -> _add($eabs, $delta);
4912             $dotpos += $LIB -> _num($delta);
4913             }
4914             }
4915 0 0       0  
    0          
4916 0         0 if ($dotpos < $mantlen) {
4917             substr $mant, $dotpos, 0, ".";
4918 0         0 } elsif ($dotpos > $mantlen) {
4919             $mant .= "0" x ($dotpos - $mantlen);
4920             }
4921 0         0  
4922             $str .= $mant . 'e' . $esgn . $LIB->_str($eabs);
4923 0         0  
4924             return $str;
4925             }
4926              
4927             # Fractional notation, e.g., "123.4375" is written as "1975/16".
4928              
4929 0 0   0 1 0 sub bfstr {
4930             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_);
4931 0 0       0  
4932             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4933              
4934             # Inf and NaN
4935 0 0 0     0  
4936 0 0       0 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
4937 0         0 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
4938             return 'inf'; # +inf
4939             }
4940              
4941             # Upgrade?
4942 0 0 0     0  
4943             return $upgrade -> bfstr($x, @r)
4944             if defined($upgrade) && !$x -> isa(__PACKAGE__);
4945              
4946             # Finite number
4947 0 0       0  
4948             my $str = $x->{sign} eq '-' ? '-' : '';
4949 0 0       0  
4950 0         0 if ($x->{_es} eq '+') {
4951             $str .= $LIB -> _str($x->{_m}) . ("0" x $LIB -> _num($x->{_e}));
4952 0         0 } else {
4953 0         0 my @flt_parts = ($x->{sign}, $x->{_m}, $x->{_es}, $x->{_e});
4954 0         0 my @rat_parts = $class -> _flt_lib_parts_to_rat_lib_parts(@flt_parts);
4955 0 0       0 $str = $LIB -> _str($rat_parts[1]) . "/" . $LIB -> _str($rat_parts[2]);
4956             $str = "-" . $str if $rat_parts[0] eq "-";
4957             }
4958 0         0  
4959             return $str;
4960             }
4961              
4962             sub to_hex {
4963 36 50   36 1 525 # return number as hexadecimal string (only for integers defined)
4964             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_);
4965 36 50       90  
4966             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4967              
4968             # Inf and NaN
4969 36 100 100     177  
4970 12 100       102 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
4971 4         41 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
4972             return 'inf'; # +inf
4973             }
4974              
4975             # Upgrade?
4976 24 50 33     64  
4977             return $upgrade -> to_hex($x, @r)
4978             if defined($upgrade) && !$x -> isa(__PACKAGE__);
4979              
4980             # Finite number
4981 24 100       62  
4982             return '0' if $x->is_zero();
4983 16 50       55  
4984             return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex?
4985 16         56  
4986 16 50       52 my $z = $LIB->_copy($x->{_m});
4987 0         0 if (! $LIB->_is_zero($x->{_e})) { # > 0
4988             $z = $LIB->_lsft($z, $x->{_e}, 10);
4989 16         100 }
4990 16 100       189 my $str = $LIB->_to_hex($z);
4991             return $x->{sign} eq '-' ? "-$str" : $str;
4992             }
4993              
4994             sub to_oct {
4995 40 50   40 1 587 # return number as octal digit string (only for integers defined)
4996             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_);
4997 40 50       130  
4998             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4999              
5000             # Inf and NaN
5001 40 100 100     166  
5002 12 100       118 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
5003 4         39 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
5004             return 'inf'; # +inf
5005             }
5006              
5007             # Upgrade?
5008 28 50 33     378  
5009             return $upgrade -> to_hex($x, @r)
5010             if defined($upgrade) && !$x -> isa(__PACKAGE__);
5011              
5012             # Finite number
5013 28 100       75  
5014             return '0' if $x->is_zero();
5015 20 50       91  
5016             return $nan if $x->{_es} ne '+'; # how to do 1e-1 in octal?
5017 20         66  
5018 20 50       61 my $z = $LIB->_copy($x->{_m});
5019 0         0 if (! $LIB->_is_zero($x->{_e})) { # > 0
5020             $z = $LIB->_lsft($z, $x->{_e}, 10);
5021 20         102 }
5022 20 100       261 my $str = $LIB->_to_oct($z);
5023             return $x->{sign} eq '-' ? "-$str" : $str;
5024             }
5025              
5026             sub to_bin {
5027 40 50   40 1 569 # return number as binary digit string (only for integers defined)
5028             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), $_[0]) : objectify(1, @_);
5029 40 50       111  
5030             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
5031              
5032             # Inf and NaN
5033 40 100 100     160  
5034 12 100       121 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
5035 4         49 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
5036             return 'inf'; # +inf
5037             }
5038              
5039             # Upgrade?
5040 28 50 33     77  
5041             return $upgrade -> to_hex($x, @r)
5042             if defined($upgrade) && !$x -> isa(__PACKAGE__);
5043              
5044             # Finite number
5045 28 100       94  
5046             return '0' if $x->is_zero();
5047 20 50       65  
5048             return $nan if $x->{_es} ne '+'; # how to do 1e-1 in binary?
5049 20         63  
5050 20 50       57 my $z = $LIB->_copy($x->{_m});
5051 0         0 if (! $LIB->_is_zero($x->{_e})) { # > 0
5052             $z = $LIB->_lsft($z, $x->{_e}, 10);
5053 20         101 }
5054 20 100       229 my $str = $LIB->_to_bin($z);
5055             return $x->{sign} eq '-' ? "-$str" : $str;
5056             }
5057              
5058 0 0   0 1 0 sub to_ieee754 {
5059             my ($class, $x, $format, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
5060 0 0       0  
5061             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
5062 0         0  
5063             my $enc; # significand encoding (applies only to decimal)
5064 0         0 my $k; # storage width in bits
5065             my $b; # base
5066 0 0       0  
    0          
    0          
    0          
    0          
    0          
    0          
    0          
5067 0         0 if ($format =~ /^binary(\d+)\z/) {
5068 0         0 $k = $1;
5069             $b = 2;
5070 0         0 } elsif ($format =~ /^decimal(\d+)(dpd|bcd)?\z/) {
5071 0         0 $k = $1;
5072 0   0     0 $b = 10;
5073             $enc = $2 || 'dpd'; # default is dencely-packed decimals (DPD)
5074 0         0 } elsif ($format eq 'half') {
5075 0         0 $k = 16;
5076             $b = 2;
5077 0         0 } elsif ($format eq 'single') {
5078 0         0 $k = 32;
5079             $b = 2;
5080 0         0 } elsif ($format eq 'double') {
5081 0         0 $k = 64;
5082             $b = 2;
5083 0         0 } elsif ($format eq 'quadruple') {
5084 0         0 $k = 128;
5085             $b = 2;
5086 0         0 } elsif ($format eq 'octuple') {
5087 0         0 $k = 256;
5088             $b = 2;
5089 0         0 } elsif ($format eq 'sexdecuple') {
5090 0         0 $k = 512;
5091             $b = 2;
5092             }
5093 0 0       0  
5094             if ($b == 2) {
5095              
5096             # Get the parameters for this format.
5097 0         0  
5098             my $p; # precision (in bits)
5099 0         0 my $t; # number of bits in significand
5100             my $w; # number of bits in exponent
5101 0 0       0  
    0          
    0          
5102 0         0 if ($k == 16) { # binary16 (half-precision)
5103 0         0 $p = 11;
5104 0         0 $t = 10;
5105             $w = 5;
5106 0         0 } elsif ($k == 32) { # binary32 (single-precision)
5107 0         0 $p = 24;
5108 0         0 $t = 23;
5109             $w = 8;
5110 0         0 } elsif ($k == 64) { # binary64 (double-precision)
5111 0         0 $p = 53;
5112 0         0 $t = 52;
5113             $w = 11;
5114 0 0 0     0 } else { # binaryN (quadruple-precition and above)
5115 0         0 if ($k < 128 || $k != 32 * sprintf('%.0f', $k / 32)) {
5116             croak "Number of bits must be 16, 32, 64, or >= 128 and",
5117             " a multiple of 32";
5118 0         0 }
5119 0         0 $p = $k - sprintf('%.0f', 4 * log($k) / log(2)) + 13;
5120 0         0 $t = $p - 1;
5121             $w = $k - $t - 1;
5122             }
5123              
5124             # The maximum exponent, minimum exponent, and exponent bias.
5125 0         0  
5126 0         0 my $emax = $class -> new(2) -> bpow($w - 1) -> bdec();
5127 0         0 my $emin = 1 - $emax;
5128             my $bias = $emax;
5129              
5130             # Get numerical sign, exponent, and mantissa/significand for bit
5131             # string.
5132 0         0  
5133 0         0 my $sign = 0;
5134             my $expo;
5135             my $mant;
5136 0 0       0  
    0          
    0          
5137 0         0 if ($x -> is_nan()) { # nan
5138 0         0 $sign = 1;
5139 0         0 $expo = $emax -> copy() -> binc();
5140             $mant = $class -> new(2) -> bpow($t - 1);
5141 0 0       0 } elsif ($x -> is_inf()) { # inf
5142 0         0 $sign = 1 if $x -> is_neg();
5143 0         0 $expo = $emax -> copy() -> binc();
5144             $mant = $class -> bzero();
5145 0         0 } elsif ($x -> is_zero()) { # zero
5146 0         0 $expo = $emin -> copy() -> bdec();
5147             $mant = $class -> bzero();
5148             } else { # normal and subnormal
5149 0 0       0  
5150             $sign = 1 if $x -> is_neg();
5151              
5152             # Now we need to compute the mantissa and exponent in base $b.
5153 0         0  
5154 0         0 my $binv = $class -> new("0.5");
5155 0         0 my $b = $class -> new(2);
5156             my $one = $class -> bone();
5157              
5158             # We start off by initializing the exponent to zero and the
5159             # mantissa to the input value. Then we increase the mantissa and
5160             # decrease the exponent, or vice versa, until the mantissa is in
5161             # the desired range or we hit one of the limits for the exponent.
5162 0         0  
5163             $mant = $x -> copy() -> babs();
5164              
5165             # We need to find the base 2 exponent. First make an estimate of
5166             # the base 2 exponent, before adjusting it below. We could skip
5167             # this estimation and go straight to the while-loops below, but the
5168             # loops are slow, especially when the final exponent is far from
5169             # zero and even more so if the number of digits is large. This
5170             # initial estimation speeds up the computation dramatically.
5171             #
5172             # log2($m * 10**$e) = log10($m + 10**$e) * log(10)/log(2)
5173             # = (log10($m) + $e) * log(10)/log(2)
5174             # = (log($m)/log(10) + $e) * log(10)/log(2)
5175 0         0  
5176 0         0 my ($m, $e) = $x -> nparts();
5177 0         0 my $ms = $m -> numify();
5178             my $es = $e -> numify();
5179 0         0  
5180 0         0 my $expo_est = (log(abs($ms))/log(10) + $es) * log(10)/log(2);
5181             $expo_est = int($expo_est);
5182              
5183             # Limit the exponent.
5184 0 0       0  
    0          
5185 0         0 if ($expo_est > $emax) {
5186             $expo_est = $emax;
5187 0         0 } elsif ($expo_est < $emin) {
5188             $expo_est = $emin;
5189             }
5190              
5191             # Don't multiply by a number raised to a negative exponent. This
5192             # will cause a division, whose result is truncated to some fixed
5193             # number of digits. Instead, multiply by the inverse number raised
5194             # to a positive exponent.
5195 0         0  
5196 0 0       0 $expo = $class -> new($expo_est);
    0          
5197 0         0 if ($expo_est > 0) {
5198             $mant = $mant -> bmul($binv -> copy() -> bpow($expo));
5199 0         0 } elsif ($expo_est < 0) {
5200 0         0 my $expo_abs = $expo -> copy() -> bneg();
5201             $mant = $mant -> bmul($b -> copy() -> bpow($expo_abs));
5202             }
5203              
5204             # Final adjustment of the estimate above.
5205 0   0     0  
5206 0         0 while ($mant >= $b && $expo <= $emax) {
5207 0         0 $mant = $mant -> bmul($binv);
5208             $expo = $expo -> binc();
5209             }
5210 0   0     0  
5211 0         0 while ($mant < $one && $expo >= $emin) {
5212 0         0 $mant = $mant -> bmul($b);
5213             $expo = $expo -> bdec();
5214             }
5215              
5216             # This is when the magnitude is larger than what can be represented
5217             # in this format. Encode as infinity.
5218 0 0       0  
    0          
5219 0         0 if ($expo > $emax) {
5220 0         0 $mant = $class -> bzero();
5221             $expo = $emax -> copy() -> binc();
5222             }
5223              
5224             # This is when the magnitude is so small that the number is encoded
5225             # as a subnormal number.
5226             #
5227             # If the magnitude is smaller than that of the smallest subnormal
5228             # number, and rounded downwards, it is encoded as zero. This works
5229             # transparently and does not need to be treated as a special case.
5230             #
5231             # If the number is between the largest subnormal number and the
5232             # smallest normal number, and the value is rounded upwards, the
5233             # value must be encoded as a normal number. This must be treated as
5234             # a special case.
5235              
5236             elsif ($expo < $emin) {
5237              
5238             # Scale up the mantissa (significand), and round to integer.
5239 0         0  
5240 0         0 my $const = $class -> new($b) -> bpow($t - 1);
5241 0         0 $mant = $mant -> bmul($const);
5242             $mant = $mant -> bfround(0);
5243              
5244             # If the mantissa overflowed, encode as the smallest normal
5245             # number.
5246 0 0       0  
5247 0         0 if ($mant == $const -> bmul($b)) {
5248 0         0 $mant = $mant -> bzero();
5249             $expo = $expo -> binc();
5250             }
5251             }
5252              
5253             # This is when the magnitude is within the range of what can be
5254             # encoded as a normal number.
5255              
5256             else {
5257              
5258             # Remove implicit leading bit, scale up the mantissa
5259             # (significand) to an integer, and round.
5260 0         0  
5261 0         0 $mant = $mant -> bdec();
5262 0         0 my $const = $class -> new($b) -> bpow($t);
5263             $mant = $mant -> bmul($const) -> bfround(0);
5264              
5265             # If the mantissa overflowed, encode as the next larger value.
5266             # This works correctly also when the next larger value is
5267             # infinity.
5268 0 0       0  
5269 0         0 if ($mant == $const) {
5270 0         0 $mant = $mant -> bzero();
5271             $expo = $expo -> binc();
5272             }
5273             }
5274             }
5275 0         0  
5276             $expo = $expo -> badd($bias); # add bias
5277 0         0  
5278             my $signbit = "$sign";
5279 0         0  
5280 0         0 my $mantbits = $mant -> to_bin();
5281             $mantbits = ("0" x ($t - CORE::length($mantbits))) . $mantbits;
5282 0         0  
5283 0         0 my $expobits = $expo -> to_bin();
5284             $expobits = ("0" x ($w - CORE::length($expobits))) . $expobits;
5285 0         0  
5286 0         0 my $bin = $signbit . $expobits . $mantbits;
5287             return pack "B*", $bin;
5288             }
5289 0         0  
5290             croak("The format '$format' is not yet supported.");
5291             }
5292              
5293             sub as_hex {
5294             # return number as hexadecimal string (only for integers defined)
5295 36 50   36 1 476  
5296             my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
5297 36 50       87  
5298             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
5299 36 100       143  
5300 24 100       69 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
5301             return '0x0' if $x->is_zero();
5302 16 50       47  
5303             return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex?
5304 16         52  
5305 16 50       53 my $z = $LIB->_copy($x->{_m});
5306 0         0 if (! $LIB->_is_zero($x->{_e})) { # > 0
5307             $z = $LIB->_lsft($z, $x->{_e}, 10);
5308 16         71 }
5309 16 100       196 my $str = $LIB->_as_hex($z);
5310             return $x->{sign} eq '-' ? "-$str" : $str;
5311             }
5312              
5313             sub as_oct {
5314             # return number as octal digit string (only for integers defined)
5315 40 50   40 1 549  
5316             my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
5317 40 50       104  
5318             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
5319 40 100       158  
5320 28 100       75 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
5321             return '00' if $x->is_zero();
5322 20 50       70  
5323             return $nan if $x->{_es} ne '+'; # how to do 1e-1 in octal?
5324 20         60  
5325 20 50       63 my $z = $LIB->_copy($x->{_m});
5326 0         0 if (! $LIB->_is_zero($x->{_e})) { # > 0
5327             $z = $LIB->_lsft($z, $x->{_e}, 10);
5328 20         82 }
5329 20 100       221 my $str = $LIB->_as_oct($z);
5330             return $x->{sign} eq '-' ? "-$str" : $str;
5331             }
5332              
5333             sub as_bin {
5334             # return number as binary digit string (only for integers defined)
5335 40 50   40 1 615  
5336             my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
5337 40 50       103  
5338             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
5339 40 100       194  
5340 28 100       76 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
5341             return '0b0' if $x->is_zero();
5342 20 50       58  
5343             return $nan if $x->{_es} ne '+'; # how to do 1e-1 in binary?
5344 20         63  
5345 20 50       60 my $z = $LIB->_copy($x->{_m});
5346 0         0 if (! $LIB->_is_zero($x->{_e})) { # > 0
5347             $z = $LIB->_lsft($z, $x->{_e}, 10);
5348 20         68 }
5349 20 100       236 my $str = $LIB->_as_bin($z);
5350             return $x->{sign} eq '-' ? "-$str" : $str;
5351             }
5352              
5353             sub numify {
5354             # Make a Perl scalar number from a Math::BigFloat object.
5355 483 50   483 1 1697  
5356             my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
5357 483 50       1132  
5358             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
5359 483 50       1235  
5360 0         0 if ($x -> is_nan()) {
5361 0         0 require Math::Complex;
5362 0         0 my $inf = $Math::Complex::Inf;
5363             return $inf - $inf;
5364             }
5365 483 50       1181  
5366 0         0 if ($x -> is_inf()) {
5367 0         0 require Math::Complex;
5368 0 0       0 my $inf = $Math::Complex::Inf;
5369             return $x -> is_negative() ? -$inf : $inf;
5370             }
5371              
5372             # Create a string and let Perl's atoi()/atof() handle the rest.
5373 483         1328  
5374             return 0 + $x -> bnstr();
5375             }
5376              
5377             ###############################################################################
5378             # Private methods and functions.
5379             ###############################################################################
5380              
5381 55     55   1990 sub import {
5382 55         123 my $class = shift;
5383             $IMPORT++; # remember we did import()
5384 55         610  
5385 55         112 my @import = ('objectify');
5386             my @a; # unrecognized arguments
5387 55         296  
5388 17         44 while (@_) {
5389             my $param = shift;
5390              
5391             # Enable overloading of constants.
5392 17 50       74  
5393             if ($param eq ':constant') {
5394             overload::constant
5395              
5396 0     0   0 integer => sub {
5397             $class -> new(shift);
5398             },
5399              
5400 0     0   0 float => sub {
5401             $class -> new(shift);
5402             },
5403              
5404             binary => sub {
5405             # E.g., a literal 0377 shall result in an object whose value
5406 0 0   0   0 # is decimal 255, but new("0377") returns decimal 377.
5407 0         0 return $class -> from_oct($_[0]) if $_[0] =~ /^0_*[0-7]/;
5408 0         0 $class -> new(shift);
5409 0         0 };
5410             next;
5411             }
5412              
5413             # Upgrading.
5414 17 100       53  
5415 2         38 if ($param eq 'upgrade') {
5416 2         8 $class -> upgrade(shift);
5417             next;
5418             }
5419              
5420             # Downgrading.
5421 15 100       40  
5422 1         9 if ($param eq 'downgrade') {
5423 1         3 $class -> downgrade(shift);
5424             next;
5425             }
5426              
5427             # Accuracy.
5428 14 50       38  
5429 0         0 if ($param eq 'accuracy') {
5430 0         0 $class -> accuracy(shift);
5431             next;
5432             }
5433              
5434             # Precision.
5435 14 50       44  
5436 0         0 if ($param eq 'precision') {
5437 0         0 $class -> precision(shift);
5438             next;
5439             }
5440              
5441             # Rounding mode.
5442 14 50       37  
5443 0         0 if ($param eq 'round_mode') {
5444 0         0 $class -> round_mode(shift);
5445             next;
5446             }
5447              
5448             # Backend library.
5449 14 100       173  
5450 11         37 if ($param =~ /^(lib|try|only)\z/) {
5451 11 50       37 push @import, $param;
5452 11         45 push @import, shift() if @_;
5453             next;
5454             }
5455 3 50       10  
5456             if ($param eq 'with') {
5457             # alternative class for our private parts()
5458             # XXX: no longer supported
5459             # $LIB = shift() || 'Calc';
5460 3         5 # carp "'with' is no longer supported, use 'lib', 'try', or 'only'";
5461 3         12 shift;
5462             next;
5463             }
5464              
5465             # Unrecognized parameter.
5466 0         0  
5467             push @a, $param;
5468             }
5469 55         1715  
5470             Math::BigInt -> import(@import);
5471              
5472 55         351 # find out which one was actually loaded
5473             $LIB = Math::BigInt -> config('lib');
5474 55         63406  
5475             $class->export_to_level(1, $class, @a); # export wanted functions
5476             }
5477              
5478             sub _len_to_steps {
5479             # Given D (digits in decimal), compute N so that N! (N factorial) is
5480 3     3   5 # at least D digits long. D should be at least 50.
5481             my $d = shift;
5482              
5483 3         7 # two constants for the Ramanujan estimate of ln(N!)
5484 3         5 my $lg2 = log(2 * 3.14159265) / 2;
5485             my $lg10 = log(10);
5486              
5487 3         7 # D = 50 => N => 42, so L = 40 and R = 50
5488 3         6 my $l = 40;
5489             my $r = $d;
5490              
5491             # Otherwise this does not work under -Mbignum and we do not yet have "no
5492 3 50       10 # bignum;" :(
5493 3 50       8 $l = $l->numify if ref($l);
5494 3 50       9 $r = $r->numify if ref($r);
5495 3 50       8 $lg2 = $lg2->numify if ref($lg2);
5496             $lg10 = $lg10->numify if ref($lg10);
5497              
5498             # binary search for the right value (could this be written as the reverse of
5499 3         12 # lg(n!)?)
5500 15         29 while ($r - $l > 1) {
5501 15         46 my $n = int(($r - $l) / 2) + $l;
5502             my $ramanujan
5503             = int(($n * log($n) - $n + log($n * (1 + 4*$n*(1+2*$n))) / 6 + $lg2)
5504 15 100       42 / $lg10);
5505             $ramanujan > $d ? $r = $n : $l = $n;
5506 3         8 }
5507             $l;
5508             }
5509              
5510             sub _log {
5511             # internal log function to calculate ln() based on Taylor series.
5512 131     131   402 # Modifies $x in place.
5513 131         308 my ($x, $scale) = @_;
5514             my $class = ref $x;
5515              
5516 131 100       377 # in case of $x == 1, result is 0
5517             return $x -> bzero() if $x -> is_one();
5518              
5519             # XXX TODO: rewrite this in a similar manner to bexp()
5520              
5521             # http://www.efunda.com/math/taylor_series/logarithmic.cfm?search_string=log
5522              
5523             # u = x-1, v = x+1
5524             # _ _
5525             # Taylor: | u 1 u^3 1 u^5 |
5526             # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 0
5527             # |_ v 3 v^3 5 v^5 _|
5528              
5529             # This takes much more steps to calculate the result and is thus not used
5530             # u = x-1
5531             # _ _
5532             # Taylor: | u 1 u^2 1 u^3 |
5533             # ln (x) = 2 | --- + - * --- + - * --- + ... | x > 1/2
5534             # |_ x 2 x^2 3 x^3 _|
5535              
5536 108         469 # scale used in intermediate computations
5537             my $scaleup = $scale + 4;
5538 108         306  
5539             my ($v, $u, $numer, $denom, $factor, $f);
5540 108         309  
5541 108         767 $v = $x -> copy();
5542 108         752 $v = $v -> binc(); # v = x+1
5543 108         529 $x = $x -> bdec();
5544             $u = $x -> copy(); # u = x-1; x = x-1
5545 108         612  
5546             $x = $x -> bdiv($v, $scaleup); # first term: u/v
5547 108         508  
5548 108         511 $numer = $u -> copy(); # numerator
5549             $denom = $v -> copy(); # denominator
5550 108         636  
5551 108         671 $u = $u -> bmul($u); # u^2
5552             $v = $v -> bmul($v); # v^2
5553 108         549  
5554 108         616 $numer = $numer -> bmul($u); # u^3
5555             $denom = $denom -> bmul($v); # v^3
5556 108         622  
5557 108         612 $factor = $class -> new(3);
5558             $f = $class -> new(2);
5559 108         379  
5560 3114         7767 while (1) {
5561             my $next = $numer -> copy() -> bround($scaleup)
5562             -> bdiv($denom -> copy() -> bmul($factor) -> bround($scaleup), $scaleup);
5563 3114         12259  
5564 3114         5103 $next->{_a} = undef;
5565 3114         6951 $next->{_p} = undef;
5566 3114         7531 my $x_prev = $x -> copy();
5567             $x = $x -> badd($next);
5568 3114 100       7984  
5569             last if $x -> bacmp($x_prev) == 0;
5570              
5571 3006         7934 # calculate things for the next term
5572 3006         6846 $numer = $numer -> bmul($u);
5573 3006         8215 $denom = $denom -> bmul($v);
5574             $factor = $factor -> badd($f);
5575             }
5576 108         636  
5577 108         691 $x = $x -> bmul($f); # $x *= 2
5578             $x = $x -> bround($scale);
5579             }
5580              
5581             sub _log_10 {
5582             # Internal log function based on reducing input to the range of 0.1 .. 9.99
5583 168     168   443 # and then "correcting" the result to the proper one. Modifies $x in place.
5584 168         366 my ($x, $scale) = @_;
5585             my $class = ref $x;
5586              
5587             # Taking blog() from numbers greater than 10 takes a *very long* time, so we
5588             # break the computation down into parts based on the observation that:
5589             # blog(X*Y) = blog(X) + blog(Y)
5590             # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
5591             # $x is the faster it gets. Since 2*$x takes about 10 times as
5592             # long, we make it faster by about a factor of 100 by dividing $x by 10.
5593              
5594             # The same observation is valid for numbers smaller than 0.1, e.g. computing
5595             # log(1) is fastest, and the further away we get from 1, the longer it
5596             # takes. So we also 'break' this down by multiplying $x with 10 and subtract
5597             # the log(10) afterwards to get the correct result.
5598              
5599             # To get $x even closer to 1, we also divide by 2 and then use log(2) to
5600             # correct for this. For instance if $x is 2.4, we use the formula:
5601             # blog(2.4 * 2) == blog(1.2) + blog(2)
5602             # and thus calculate only blog(1.2) and blog(2), which is faster in total
5603             # than calculating blog(2.4).
5604              
5605             # In addition, the values for blog(2) and blog(10) are cached.
5606              
5607             # Calculate the number of digits before the dot, i.e., 1 + floor(log10(x)):
5608             # x = 123 => dbd = 3
5609             # x = 1.23 => dbd = 1
5610             # x = 0.0123 => dbd = -1
5611             # x = 0.000123 => dbd = -3
5612             # etc.
5613 168         654  
5614 168 100       824 my $dbd = $LIB->_num($x->{_e});
5615 168         580 $dbd = -$dbd if $x->{_es} eq '-';
5616             $dbd += $LIB->_len($x->{_m});
5617              
5618             # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
5619             # infinite recursion
5620 168         441  
5621             my $calc = 1; # do some calculation?
5622              
5623             # No upgrading or downgrading in the intermediate computations.
5624 168         380  
5625 168         309 local $Math::BigInt::upgrade = undef;
5626             local $Math::BigFloat::downgrade = undef;
5627              
5628             # disable the shortcut for 10, since we need log(10) and this would recurse
5629 168 100 66     898 # infinitely deep
      100        
5630             if ($x->{_es} eq '+' && # $x == 10
5631             ($LIB->_is_one($x->{_e}) &&
5632             $LIB->_is_one($x->{_m})))
5633 7         23 {
5634             $dbd = 0; # disable shortcut
5635 7 50       41 # we can use the cached value in these cases
5636 7         39 if ($scale <= $LOG_10_A) {
5637 7         30 $x = $x->bzero();
5638 7         25 $x = $x->badd($LOG_10); # modify $x in place
5639             $calc = 0; # no need to calc, but round
5640             }
5641             # if we can't use the shortcut, we continue normally
5642             } else {
5643 161 100 100     619 # disable the shortcut for 2, since we maybe have it cached
5644             if (($LIB->_is_zero($x->{_e}) && # $x == 2
5645             $LIB->_is_two($x->{_m})))
5646 32         96 {
5647             $dbd = 0; # disable shortcut
5648 32 50       156 # we can use the cached value in these cases
5649 32         111 if ($scale <= $LOG_2_A) {
5650 32         222 $x = $x->bzero();
5651 32         112 $x = $x->badd($LOG_2); # modify $x in place
5652             $calc = 0; # no need to calc, but round
5653             }
5654             # if we can't use the shortcut, we continue normally
5655             }
5656             }
5657              
5658 168 100 100     1272 # if $x = 0.1, we know the result must be 0-log(10)
      100        
      100        
5659             if ($calc != 0 &&
5660             ($x->{_es} eq '-' && # $x == 0.1
5661             ($LIB->_is_one($x->{_e}) &&
5662             $LIB->_is_one($x->{_m}))))
5663 2         5 {
5664             $dbd = 0; # disable shortcut
5665 2 50       10 # we can use the cached value in these cases
5666 2         8 if ($scale <= $LOG_10_A) {
5667 2         19 $x = $x->bzero();
5668 2         4 $x = $x->bsub($LOG_10);
5669             $calc = 0; # no need to calc, but round
5670             }
5671             }
5672 168 100       577  
5673             return $x if $calc == 0; # already have the result
5674              
5675 127         317 # default: these correction factors are undef and thus not used
5676             my $l_10; # value of ln(10) to A of $scale
5677             my $l_2; # value of ln(2) to A of $scale
5678 127         424  
5679             my $two = $class->new(2);
5680              
5681             # $x == 2 => 1, $x == 13 => 2, $x == 0.1 => 0, $x == 0.01 => -1
5682 127 100 100     1367 # so don't do this shortcut for 1 or 0
5683             if (($dbd > 1) || ($dbd < 0)) {
5684             # convert our cached value to an object if not already (avoid doing this
5685 48 100       235 # at import() time, since not everybody needs this)
5686             $LOG_10 = $class->new($LOG_10, undef, undef) unless ref $LOG_10;
5687              
5688             # got more than one digit before the dot, or more than one zero after
5689             # the dot, so do:
5690             # log(123) == log(1.23) + log(10) * 2
5691             # log(0.0123) == log(1.23) - log(10) * 2
5692 48 100       225  
5693             if ($scale <= $LOG_10_A) {
5694 47         148 # use cached value
5695             $l_10 = $LOG_10->copy(); # copy for mul
5696             } else {
5697             # else: slower, compute and cache result
5698              
5699             # shorten the time to calculate log(10) based on the following:
5700             # log(1.25 * 8) = log(1.25) + log(8)
5701             # = log(1.25) + log(2) + log(2) + log(2)
5702              
5703 1 50       5 # first get $l_2 (and possible compute and cache log(2))
5704 1 50       5 $LOG_2 = $class->new($LOG_2, undef, undef) unless ref $LOG_2;
5705             if ($scale <= $LOG_2_A) {
5706 1         4 # use cached value
5707             $l_2 = $LOG_2->copy(); # copy() for the mul below
5708             } else {
5709 0         0 # else: slower, compute and cache result
5710 0         0 $l_2 = $two->copy();
5711 0         0 $l_2 = $l_2->_log($scale); # scale+4, actually
5712             $LOG_2 = $l_2->copy(); # cache the result for later
5713 0         0 # the copy() is for mul below
5714             $LOG_2_A = $scale;
5715             }
5716              
5717 1         5 # now calculate log(1.25):
5718 1         6 $l_10 = $class->new('1.25');
5719             $l_10 = $l_10->_log($scale); # scale+4, actually
5720              
5721 1         6 # log(1.25) + log(2) + log(2) + log(2):
5722 1         6 $l_10 = $l_10->badd($l_2);
5723 1         9 $l_10 = $l_10->badd($l_2);
5724 1         7 $l_10 = $l_10->badd($l_2);
5725             $LOG_10 = $l_10->copy(); # cache the result for later
5726 1         7 # the copy() is for mul below
5727             $LOG_10_A = $scale;
5728 48 100       228 }
5729 48         159 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
5730 48         287 $l_10 = $l_10->bmul($class->new($dbd)); # log(10) * (digits_before_dot-1)
5731 48 100       186 my $dbd_sign = '+';
5732 7         20 if ($dbd < 0) {
5733 7         19 $dbd = -$dbd;
5734             $dbd_sign = '-';
5735             }
5736 48         225 ($x->{_e}, $x->{_es}) =
5737             $LIB -> _ssub($x->{_e}, $x->{_es}, $LIB->_new($dbd), $dbd_sign);
5738             }
5739              
5740             # Now: 0.1 <= $x < 10 (and possible correction in l_10)
5741              
5742             ### Since $x in the range 0.5 .. 1.5 is MUCH faster, we do a repeated div
5743             ### or mul by 2 (maximum times 3, since x < 10 and x > 0.1)
5744 127 100       564  
5745             $HALF = $class->new($HALF) unless ref($HALF);
5746 127         321  
5747 127         444 my $twos = 0; # default: none (0 times)
5748 40         103 while ($x->bacmp($HALF) <= 0) { # X <= 0.5
5749 40         121 $twos--;
5750             $x = $x->bmul($two);
5751 127         419 }
5752 52         143 while ($x->bacmp($two) >= 0) { # X >= 2
5753 52         279 $twos++;
5754             $x = $x->bdiv($two, $scale+4); # keep all digits
5755 127         688 }
5756             $x = $x->bround($scale+4);
5757             # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
5758 127 100       651 # So calculate correction factor based on ln(2):
5759 72 100       462 if ($twos != 0) {
5760 72 100       318 $LOG_2 = $class->new($LOG_2, undef, undef) unless ref $LOG_2;
5761             if ($scale <= $LOG_2_A) {
5762 69         222 # use cached value
5763             $l_2 = $LOG_2->copy(); # copy() for the mul below
5764             } else {
5765 3         12 # else: slower, compute and cache result
5766 3         31 $l_2 = $two->copy();
5767 3         12 $l_2 = $l_2->_log($scale); # scale+4, actually
5768             $LOG_2 = $l_2->copy(); # cache the result for later
5769 3         25 # the copy() is for mul below
5770             $LOG_2_A = $scale;
5771 72         460 }
5772             $l_2 = $l_2->bmul($twos); # * -2 => subtract, * 2 => add
5773 55         158 } else {
5774             undef $l_2;
5775             }
5776 127         828  
5777 127 100       677 $x = $x->_log($scale); # need to do the "normal" way
5778 127 100       731 $x = $x->badd($l_10) if defined $l_10; # correct it by ln(10)
5779             $x = $x->badd($l_2) if defined $l_2; # and maybe by ln(2)
5780              
5781 127         1311 # all done, $x contains now the result
5782             $x;
5783             }
5784              
5785             sub _pow {
5786 114     114   438 # Calculate a power where $y is a non-integer, like 2 ** 0.3
5787 114         287 my ($x, $y, @r) = @_;
5788             my $class = ref($x);
5789              
5790 114 100       401 # if $y == 0.5, it is sqrt($x)
5791 114 100       429 $HALF = $class->new($HALF) unless ref($HALF);
5792             return $x->bsqrt(@r, $y) if $y->bcmp($HALF) == 0;
5793              
5794             # Using:
5795             # a ** x == e ** (x * ln a)
5796              
5797             # u = y * ln x
5798             # _ _
5799             # Taylor: | u u^2 u^3 |
5800             # x ** y = 1 + | --- + --- + ----- + ... |
5801             # |_ 1 1*2 1*2*3 _|
5802              
5803 69         183 # we need to limit the accuracy to protect against overflow
5804 69         177 my $fallback = 0;
5805 69         303 my ($scale, @params);
5806             ($x, @params) = $x->_find_round_parameters(@r);
5807 69 50       299  
5808             return $x if $x->is_nan(); # error in _find_round_parameters?
5809              
5810 69 100       319 # no rounding at all, so must use fallback
5811             if (scalar @params == 0) {
5812 52         195 # simulate old behaviour
5813 52         192 $params[0] = $class->div_scale(); # and round to it as accuracy
5814 52         139 $params[1] = undef; # disable P
5815 52         139 $scale = $params[0]+4; # at least four more for proper round
5816 52         119 $params[2] = $r[2]; # round mode by caller or undef
5817             $fallback = 1; # to clear a/p afterwards
5818             } else {
5819             # the 4 below is empirical, and there might be cases where it is not
5820 17   33     71 # enough...
5821             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
5822             }
5823              
5824             # when user set globals, they would interfere with our calculation, so
5825 43     43   545 # disable them and later re-enable them
  43         120  
  43         23052  
5826 69         190 no strict 'refs';
5827 69         162 my $abr = "$class\::accuracy";
5828 69         164 my $ab = $$abr;
5829 69         158 $$abr = undef;
5830 69         160 my $pbr = "$class\::precision";
5831 69         148 my $pb = $$pbr;
5832             $$pbr = undef;
5833             # we also need to disable any set A or P on $x (_find_round_parameters took
5834 69         170 # them already into account), since these would interfere, too
5835 69         167 $x->{_a} = undef;
5836             $x->{_p} = undef;
5837              
5838             # Disabling upgrading and downgrading is no longer necessary to avoid an
5839             # infinite recursion, but it avoids unnecessary upgrading and downgrading in
5840             # the intermediate computations.
5841 69         152  
5842 69         221 local $Math::BigInt::upgrade = undef;
5843             local $Math::BigFloat::downgrade = undef;
5844 69         202  
5845             my ($limit, $v, $u, $below, $factor, $next, $over);
5846 69         192  
5847 69         440 $u = $x->copy()->blog(undef, $scale)->bmul($y);
5848 69 100       367 my $do_invert = ($u->{sign} eq '-');
5849 69         422 $u = $u->bneg() if $do_invert;
5850 69         384 $v = $class->bone(); # 1
5851 69         463 $factor = $class->new(2); # 2
5852             $x = $x->bone(); # first term: 1
5853 69         347  
5854 69         376 $below = $v->copy();
5855             $over = $u->copy();
5856 69         661  
5857 69         288 $limit = $class->new("1E-". ($scale-1));
5858             while (3 < 5) {
5859             # we calculate the next term, and add it to the last
5860             # when the next term is below our limit, it won't affect the outcome
5861 3277         8423 # anymore, so we stop:
5862 3277 100       14291 $next = $over->copy()->bdiv($below, $scale);
5863 3208         8364 last if $next->bacmp($limit) <= 0;
5864             $x = $x->badd($next);
5865 3208         10453 # calculate things for the next term
5866 3208         8604 $over *= $u;
5867 3208         8153 $below *= $factor;
5868             $factor = $factor->binc();
5869 3208 50       13150  
5870             last if $x->{sign} !~ /^[-+]$/;
5871             }
5872 69 100       470  
5873 32         130 if ($do_invert) {
5874 32         255 my $x_copy = $x->copy();
5875             $x = $x->bone->bdiv($x_copy, $scale);
5876             }
5877              
5878 69 50       471 # shortcut to not run through _find_round_parameters again
5879 69         333 if (defined $params[0]) {
5880             $x = $x->bround($params[0], $params[2]); # then round accordingly
5881 0         0 } else {
5882             $x = $x->bfround($params[1], $params[2]); # then round accordingly
5883 69 100       474 }
5884             if ($fallback) {
5885 52         183 # clear a/p after round, since user did not request it
5886 52         123 $x->{_a} = undef;
5887             $x->{_p} = undef;
5888             }
5889 69         290 # restore globals
5890 69         203 $$abr = $ab;
5891 69         2303 $$pbr = $pb;
5892             $x;
5893             }
5894              
5895             # These functions are only provided for backwards compabibility so that old
5896             # version of Math::BigRat etc. don't complain about missing them.
5897              
5898 0     0     sub _e_add {
5899 0           my ($x, $y, $xs, $ys) = @_;
5900             return $LIB -> _sadd($x, $xs, $y, $ys);
5901             }
5902              
5903 0     0     sub _e_sub {
5904 0           my ($x, $y, $xs, $ys) = @_;
5905             return $LIB -> _ssub($x, $xs, $y, $ys);
5906             }
5907              
5908             1;
5909              
5910             __END__