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   1074036 use 5.006001;
  43         301  
16 43     43   254 use strict;
  43         94  
  43         1152  
17 43     43   228 use warnings;
  43         119  
  43         1646  
18              
19 43     43   260 use Carp qw< carp croak >;
  43         87  
  43         3126  
20 43     43   326 use Scalar::Util qw< blessed >;
  43         109  
  43         2503  
21 43     43   31185 use Math::BigInt qw< >;
  43         146  
  43         77446  
22              
23             our $VERSION = '1.999840';
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   178 '+' => sub { $_[0] -> copy() -> badd($_[1]); },
39              
40 73     73   1417 '-' => sub { my $c = $_[0] -> copy();
41 73 50       384 $_[2] ? $c -> bneg() -> badd($_[1])
42             : $c -> bsub($_[1]); },
43              
44 153     153   5810 '*' => sub { $_[0] -> copy() -> bmul($_[1]); },
45              
46 12 50   12   3170 '/' => 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   1271 '**' => 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   362 '+=' => sub { $_[0] -> badd($_[1]); },
64              
65 28     28   319 '-=' => sub { $_[0] -> bsub($_[1]); },
66              
67 6424     6424   16836 '*=' => sub { $_[0] -> bmul($_[1]); },
68              
69 8     8   109 '/=' => sub { scalar $_[0] -> bdiv($_[1]); },
70              
71 8     8   161 '%=' => sub { $_[0] -> bmod($_[1]); },
72              
73 4     4   102 '**=' => sub { $_[0] -> bpow($_[1]); },
74              
75 8     8   137 '<<=' => sub { $_[0] -> bblsft($_[1]); },
76              
77 8     8   1601 '>>=' => sub { $_[0] -> bbrsft($_[1]); },
78              
79             # 'x=' => sub { },
80              
81             # '.=' => sub { },
82              
83             # overload key: num_comparison
84              
85 301 50   301   1372 '<' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> blt($_[0])
86             : $_[0] -> blt($_[1]); },
87              
88 857 50   857   3132 '<=' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> ble($_[0])
89             : $_[0] -> ble($_[1]); },
90              
91 879 50   879   3890 '>' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bgt($_[0])
92             : $_[0] -> bgt($_[1]); },
93              
94 177 100   177   739 '>=' => sub { $_[2] ? ref($_[0]) -> new($_[1]) -> bge($_[0])
95             : $_[0] -> bge($_[1]); },
96              
97 186     186   5338 '==' => sub { $_[0] -> beq($_[1]); },
98              
99 9     9   376 '!=' => 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   1512724 '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   988 '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   139 '++' => 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   555 'log' => sub { $_[0] -> copy() -> blog(); },
186              
187 0     0   0 'sqrt' => sub { $_[0] -> copy() -> bsqrt(); },
188              
189 140     140   424 'int' => sub { $_[0] -> copy() -> bint(); },
190              
191             # overload key: conversion
192              
193 280 50   280   614 'bool' => sub { $_[0] -> is_zero() ? '' : 1; },
194              
195 742     742   80172 '""' => sub { $_[0] -> bstr(); },
196              
197 8     8   44 '0+' => sub { $_[0] -> numify(); },
198              
199 0     0   0 '=' => sub { $_[0] -> copy(); },
200              
201 43     43   408 ;
  43         92  
  43         4051  
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   197 my ($class) = @_;
247 43         206 bless \$round_mode, $class;
248             }
249              
250             sub FETCH {
251 1     1   685 return $round_mode;
252             }
253              
254             sub STORE {
255 1     1   768 $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   34059 $rnd_mode = 'even';
262 43         271 tie $rnd_mode, 'Math::BigFloat';
263              
264 43         4434 *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   10 my $name = $AUTOLOAD;
274 3         58 $name =~ s/(.*):://; # split package
275 3   50     22 my $c = $1 || __PACKAGE__;
276 43     43   370 no strict 'refs';
  43         99  
  43         178144  
277 3 50       24 $c->import() if $IMPORT == 0;
278 3 50       19 if (!_method_alias($name)) {
279 3 50       10 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         28  
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   26 sub _method_alias { exists $methods{$_[0]||''}; }
318 3   50 3   19 sub _method_hand_up { exists $hand_ups{$_[0]||''}; }
319             }
320              
321             sub isa {
322 24688     24688 0 49149 my ($self, $class) = @_;
323 24688 100       59202 return if $class =~ /^Math::BigInt/; # we aren't one of these
324 23557         94695 UNIVERSAL::isa($self, $class);
325             }
326              
327             sub config {
328             # return (later set?) configuration data as hash ref
329 55   50 55 1 43313 my $class = shift || 'Math::BigFloat';
330              
331             # Getter/accessor.
332              
333 55 100 100     300 if (@_ == 1 && ref($_[0]) ne 'HASH') {
334 23         46 my $param = shift;
335 23 100       76 return $class if $param eq 'class';
336 21 100       95 return $LIB if $param eq 'with';
337 16         135 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         59 $cfg->{class} = $class;
346 31         52 $cfg->{with} = $LIB;
347 31         115 $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 6272709 my $self = shift;
362 17120         30579 my $selfref = ref $self;
363 17120   33     60695 my $class = $selfref || $self;
364              
365             # Make "require" work.
366              
367 17120 100       37386 $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       38079 return $class -> bzero() unless @_;
373              
374 17112         36829 my ($wanted, @r) = @_;
375              
376 17112 50       36411 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     60732 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       47336 $self = bless {}, $class unless $selfref;
396              
397             # Math::BigFloat or subclass
398              
399 17112 100 100     60766 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         1076 $self -> {sign} = $wanted -> {sign};
405 344         1028 $self -> {_m} = $LIB -> _copy($wanted -> {_m});
406 344         808 $self -> {_es} = $wanted -> {_es};
407 344         827 $self -> {_e} = $LIB -> _copy($wanted -> {_e});
408 344 50 66     1784 $self = $self->round(@r)
      66        
409             unless @r >= 2 && !defined($r[0]) && !defined($r[1]);
410 344         1877 return $self;
411             }
412              
413             # Shortcut for Math::BigInt and its subclasses. This should be improved.
414              
415 16768 100       41112 if (defined(blessed($wanted))) {
416 523 100       1574 if ($wanted -> isa('Math::BigInt')) {
417 522         1662 $self->{sign} = $wanted -> {sign};
418 522         1679 $self->{_m} = $LIB -> _copy($wanted -> {value});
419 522         1121 $self->{_es} = '+';
420 522         1391 $self->{_e} = $LIB -> _zero();
421 522         1521 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       85901 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       17996 return $downgrade -> new($1 . $2) if defined $downgrade;
447 7877   100     35559 $self->{sign} = $1 || '+';
448 7877         27406 $self->{_m} = $LIB -> _new($2);
449 7877         16387 $self->{_es} = '+';
450 7877         18986 $self->{_e} = $LIB -> _zero();
451 7877 100 100     36379 $self = $self->round(@r)
      100        
452             unless @r >= 2 && !defined $r[0] && !defined $r[1];
453 7877         73603 return $self;
454             }
455              
456             # Handle Infs.
457              
458 8361 100       26348 if ($wanted =~ / ^
459             \s*
460             ( [+-]? )
461             inf (?: inity )?
462             \s*
463             \z
464             /ix)
465             {
466 1727   100     6462 my $sgn = $1 || '+';
467 1727         5346 return $class -> binf($sgn, @r);
468             }
469              
470             # Handle explicit NaNs (not the ones returned due to invalid input).
471              
472 6634 100       17670 if ($wanted =~ / ^
473             \s*
474             ( [+-]? )
475             nan
476             \s*
477             \z
478             /ix)
479             {
480 475         1887 return $class -> bnan(@r);
481             }
482              
483 6159         10171 my @parts;
484              
485 6159 100 66     63044 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         24234 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts;
528              
529 5708 100 100     24749 $self = $self->round(@r)
      100        
530             unless @r >= 2 && !defined($r[0]) && !defined($r[1]);
531              
532 5708 100 100     15526 return $downgrade -> new($self -> bdstr(), @r)
533             if defined($downgrade) && $self -> is_int();
534 5706         60411 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         1410 return $class -> bnan(@r);
541             }
542              
543             sub from_dec {
544 1     1 1 2152 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     5 return $self if $selfref && $self->modify('from_dec');
551              
552 1         3 my $str = shift;
553 1         3 my @r = @_;
554              
555             # If called as a class method, initialize a new object.
556              
557 1 50       5 $self = bless {}, $class unless $selfref;
558              
559 1 50       7 if (my @parts = $class -> _dec_str_to_flt_lib_parts($str)) {
560 1         8 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts;
561              
562 1 0 33     8 $self = $self->round(@r)
      33        
563             unless @r >= 2 && !defined($r[0]) && !defined($r[1]);
564              
565 1 50 33     9 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 2227 my $self = shift;
575 1         6 my $selfref = ref $self;
576 1   33     6 my $class = $selfref || $self;
577              
578             # Don't modify constant (read-only) objects.
579              
580 1 50 33     6 return $self if $selfref && $self->modify('from_hex');
581              
582 1         2 my $str = shift;
583 1         13 my @r = @_;
584              
585             # If called as a class method, initialize a new object.
586              
587 1 50       17 $self = bless {}, $class unless $selfref;
588              
589 1 50       15 if (my @parts = $class -> _hex_str_to_flt_lib_parts($str)) {
590 1         8 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts;
591              
592 1 0 33     6 $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 2200 my $self = shift;
605 1         4 my $selfref = ref $self;
606 1   33     6 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         2 my @r = @_;
614              
615             # If called as a class method, initialize a new object.
616              
617 1 50       4 $self = bless {}, $class unless $selfref;
618              
619 1 50       10 if (my @parts = $class -> _oct_str_to_flt_lib_parts($str)) {
620 1         5 ($self->{sign}, $self->{_m}, $self->{_es}, $self->{_e}) = @parts;
621              
622 1 0 33     14 $self = $self->round(@r)
      33        
623             unless @r >= 2 && !defined($r[0]) && !defined($r[1]);
624              
625 1 50 33     18 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 2156 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     8 return $self if $selfref && $self->modify('from_bin');
641              
642 3         11 my $str = shift;
643 3         6 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       15 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     14 $self = $self->round(@r)
      33        
653             unless @r >= 2 && !defined($r[0]) && !defined($r[1]);
654              
655 3 50 33     12 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 2241 my $self = shift;
665 1         3 my $selfref = ref $self;
666 1   33     7 my $class = $selfref || $self;
667              
668             # Don't modify constant (read-only) objects.
669              
670 1 50 33     4 return $self if $selfref && $self->modify('from_ieee754');
671              
672 1         2 my $in = shift; # input string (or raw bytes)
673 1         2 my $format = shift; # format ("binary32", "decimal64" etc.)
674 1         3 my $enc; # significand encoding (applies only to decimal)
675             my $k; # storage width in bits
676 1         0 my $b; # base
677 1         3 my @r = @_; # rounding parameters, if any
678              
679 1 50       8 if ($format =~ /^binary(\d+)\z/) {
    0          
    0          
    0          
    0          
    0          
    0          
    0          
680 1         3 $k = $1;
681 1         3 $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         4 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       7 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         9 my $emax = Math::BigFloat -> new(2) -> bpow($w - 1) -> bdec();
739 1         167 my $emin = 1 - $emax;
740 1         4 my $bias = $emax;
741              
742             # Undefined input.
743              
744 1 50       6 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         3 my $len = CORE::length $in;
752 1 50       5 if (8 * $len == $k) { # bytes
    0          
    0          
753 1         7 $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       6 my $sign = substr($in, 0, 1) eq '1' ? '-' : '+';
770 1         7 my $expo = $class -> from_bin(substr($in, 1, $w));
771 1         9 my $mant = $class -> from_bin(substr($in, $w + 1));
772              
773 1         4 my $x;
774              
775 1         4 $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         4 $mant = $class -> new(2) -> bpow($t) -> badd($mant);
798 1 50       5 if ($expo < $t) {
799             # compute (1/$b)**(N) rather than ($b)**(-N)
800 1         5 $x = $class -> new("0.5"); # 1/$b
801 1         6 $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       12 $x = $x -> bneg() if $sign eq '-';
807             }
808              
809 1 50       5 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         5 $self = $x;
816             }
817              
818 1 50 33     6 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 20840 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         1522 my $self = shift;
839 603         1163 my $selfref = ref $self;
840 603   66     1554 my $class = $selfref || $self;
841              
842 603 100       1389 $self->import() if $IMPORT == 0; # make require work
843              
844             # Don't modify constant (read-only) objects.
845              
846 603 50 66     2574 return $self if $selfref && $self->modify('bzero');
847              
848             # Get the rounding parameters, if any.
849              
850 603         1237 my @r = @_;
851              
852 603 100       1410 return $downgrade -> bzero(@r) if defined $downgrade;
853              
854             # If called as a class method, initialize a new object.
855              
856 602 100       1418 $self = bless {}, $class unless $selfref;
857              
858 602         1364 $self -> {sign} = '+';
859 602         1862 $self -> {_m} = $LIB -> _zero();
860 602         1276 $self -> {_es} = '+';
861 602         1437 $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       1619 if (@r) {
870 64 50 100     353 croak "can't specify both accuracy and precision"
      66        
871             if @r >= 2 && defined($r[0]) && defined($r[1]);
872 64         152 $self->{_a} = $r[0];
873 64         157 $self->{_p} = $r[1];
874             } else {
875 538 100       1252 unless($selfref) {
876 101         322 $self->{_a} = $class -> accuracy();
877 101         342 $self->{_p} = $class -> precision();
878             }
879             }
880              
881 602         4534 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 23057 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         3916 my $self = shift;
897 1589         2894 my $selfref = ref $self;
898 1589   66     4500 my $class = $selfref || $self;
899              
900 1589 100       3461 $self->import() if $IMPORT == 0; # make require work
901              
902             # Don't modify constant (read-only) objects.
903              
904 1589 50 66     4922 return $self if $selfref && $self->modify('bone');
905              
906 1589 100       3155 return $downgrade -> bone(@_) if defined $downgrade;
907              
908             # Get the sign.
909              
910 1588         2522 my $sign = '+'; # default is to return +1
911 1588 100 100     5040 if (defined($_[0]) && $_[0] =~ /^\s*([+-])\s*$/) {
912 163         509 $sign = $1;
913 163         294 shift;
914             }
915              
916             # Get the rounding parameters, if any.
917              
918 1588         2902 my @r = @_;
919              
920             # If called as a class method, initialize a new object.
921              
922 1588 100       4167 $self = bless {}, $class unless $selfref;
923              
924 1588         3816 $self -> {sign} = $sign;
925 1588         4774 $self -> {_m} = $LIB -> _one();
926 1588         3295 $self -> {_es} = '+';
927 1588         3974 $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       3863 if (@r) {
936 29 50 100     172 croak "can't specify both accuracy and precision"
      66        
937             if @r >= 2 && defined($r[0]) && defined($r[1]);
938 29         83 $self->{_a} = $_[0];
939 29         65 $self->{_p} = $_[1];
940             } else {
941 1559 100       3390 unless($selfref) {
942 1000         2988 $self->{_a} = $class -> accuracy();
943 1000         3078 $self->{_p} = $class -> precision();
944             }
945             }
946              
947 1588         6789 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 22249 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         4236 my $self = shift;
963 2071         3360 my $selfref = ref $self;
964 2071   66     5562 my $class = $selfref || $self;
965              
966             {
967 43     43   447 no strict 'refs';
  43         108  
  43         20850  
  2071         3174  
968 2071 100       2674 if (${"${class}::_trap_inf"}) {
  2071         8554  
969 5         537 croak("Tried to create +-inf in $class->binf()");
970             }
971             }
972              
973 2066 100       4336 $self->import() if $IMPORT == 0; # make require work
974              
975             # Don't modify constant (read-only) objects.
976              
977 2066 50 66     5544 return $self if $selfref && $self->modify('binf');
978              
979 2066 100       4098 return $downgrade -> binf(@_) if $downgrade;
980              
981             # Get the sign.
982              
983 2059         3516 my $sign = '+'; # default is to return positive infinity
984 2059 100 100     10108 if (defined($_[0]) && $_[0] =~ /^\s*([+-])(inf|$)/i) {
985 1968         4560 $sign = $1;
986 1968         2899 shift;
987             }
988              
989             # Get the rounding parameters, if any.
990              
991 2059         3807 my @r = @_;
992              
993             # If called as a class method, initialize a new object.
994              
995 2059 100       5245 $self = bless {}, $class unless $selfref;
996              
997 2059         6338 $self -> {sign} = $sign . 'inf';
998 2059         6708 $self -> {_m} = $LIB -> _zero();
999 2059         4215 $self -> {_es} = '+';
1000 2059         4834 $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       5271 if (@r) {
1009 530 50 66     2557 croak "can't specify both accuracy and precision"
      33        
1010             if @r >= 2 && defined($r[0]) && defined($r[1]);
1011 530         1123 $self->{_a} = $r[0];
1012 530         1046 $self->{_p} = $r[1];
1013             } else {
1014 1529 100       3361 unless($selfref) {
1015 1234         3863 $self->{_a} = $class -> accuracy();
1016 1234         3416 $self->{_p} = $class -> precision();
1017             }
1018             }
1019              
1020 2059         22409 return $self;
1021             }
1022              
1023             sub bnan {
1024             # create/assign a 'NaN'
1025              
1026             # Class::method(...) -> Class->method(...)
1027 2115 50 66 2115 1 21989 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         4983 my $self = shift;
1036 2115         3701 my $selfref = ref $self;
1037 2115   66     5326 my $class = $selfref || $self;
1038              
1039             {
1040 43     43   400 no strict 'refs';
  43         146  
  43         414179  
  2115         3122  
1041 2115 100       2940 if (${"${class}::_trap_nan"}) {
  2115         8317  
1042 3         339 croak("Tried to create NaN in $class->bnan()");
1043             }
1044             }
1045              
1046 2112 100       4484 $self->import() if $IMPORT == 0; # make require work
1047              
1048             # Don't modify constant (read-only) objects.
1049              
1050 2112 50 66     7224 return $self if $selfref && $self->modify('bnan');
1051              
1052 2112 100       4251 return $downgrade -> bnan(@_) if defined $downgrade;
1053              
1054             # Get the rounding parameters, if any.
1055              
1056 2098         3935 my @r = @_;
1057              
1058             # If called as a class method, initialize a new object.
1059              
1060 2098 100       4891 $self = bless {}, $class unless $selfref;
1061              
1062 2098         4879 $self -> {sign} = $nan;
1063 2098         6159 $self -> {_m} = $LIB -> _zero();
1064 2098         3987 $self -> {_es} = '+';
1065 2098         4465 $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       4851 if (@r) {
1074 296 50 66     1308 croak "can't specify both accuracy and precision"
      33        
1075             if @r >= 2 && defined($r[0]) && defined($r[1]);
1076 296         689 $self->{_a} = $r[0];
1077 296         582 $self->{_p} = $r[1];
1078             } else {
1079 1802 100       3855 unless($selfref) {
1080 751         2160 $self->{_a} = $class -> accuracy();
1081 751         2253 $self->{_p} = $class -> precision();
1082             }
1083             }
1084              
1085 2098         21032 return $self;
1086             }
1087              
1088             sub bpi {
1089              
1090             # Class::method(...) -> Class->method(...)
1091 205 100 66 205 1 3805 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         7 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         459 my $self = shift;
1118 205         370 my $selfref = ref $self;
1119 205   66     578 my $class = $selfref || $self;
1120 205         463 my @r = @_; # rounding paramters
1121              
1122 205 100       415 if ($selfref) { # bpi() called as an instance method
1123 83 50       269 return $self if $self -> modify('bpi');
1124             } else { # bpi() called as a class method
1125 122         294 $self = bless {}, $class; # initialize new instance
1126             }
1127              
1128 205         608 ($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       913 my $n = defined $r[0] ? $r[0]
    100          
1134             : defined $r[1] ? 1 - $r[1]
1135             : $self -> div_scale();
1136              
1137 205 100       493 my $rmode = defined $r[2] ? $r[2] : $self -> round_mode();
1138              
1139 205         313 my $pi;
1140              
1141 205 50       430 if ($n <= 1000) {
1142              
1143             # 75 x 14 = 1050 digits
1144              
1145 205         379 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         297 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         547 my $nchrs = $n + int($n / 75);
1170              
1171             # Extract the digits we want.
1172              
1173 205         499 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       450 if ($rmode eq 'trunc') {
1181 0         0 $round_up = 0;
1182             } else {
1183 205         381 my $next_digit = substr($all_digits, $nchrs, 1);
1184 205 100       517 $round_up = $next_digit lt '5' ? 0 : 1;
1185             }
1186              
1187             # Remove the newlines.
1188              
1189 205         496 $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       498 if ($round_up) {
1196 119         594 my $last_digit = substr($digits, -1, 1);
1197 119 100       290 if ($last_digit lt '9') {
1198 108         292 substr($digits, -1, 1) = ++$last_digit;
1199             } else {
1200 11         144 $digits =~ s{([0-8])(9+)$}
  11         83  
1201             { ($1 + 1) . ("0" x CORE::length($2)) }e;
1202             }
1203             }
1204              
1205             # Convert to an object.
1206 205         737  
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       662  
    50          
1241 191         665 if (defined $r[0]) {
1242             $pi -> accuracy($r[0]);
1243 0         0 } elsif (defined $r[1]) {
1244             $pi -> precision($r[1]);
1245             }
1246 205         484  
1247 1230         2443 for my $key (qw/ sign _m _es _e _a _p /) {
1248             $self -> {$key} = $pi -> {$key};
1249             }
1250 205 50 33     611  
1251             return $downgrade -> new($self -> bdstr(), @r)
1252 205         1278 if defined($downgrade) && $self->is_int();
1253             return $self;
1254             }
1255              
1256 17365     17365 1 173675 sub copy {
1257 17365 50       34824 my ($x, $class);
1258 17365         27059 if (ref($_[0])) { # $y = $x -> copy()
1259 17365         27413 $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       35179  
1266             carp "Rounding is not supported for ", (caller(0))[3], "()" if @_;
1267 17365         36550  
1268             my $copy = bless {}, $class;
1269 17365         41834  
1270 17365         34233 $copy->{sign} = $x->{sign};
1271 17365         45429 $copy->{_es} = $x->{_es};
1272 17365         39675 $copy->{_m} = $LIB->_copy($x->{_m});
1273 17365 100       43317 $copy->{_e} = $LIB->_copy($x->{_e});
1274 17365 100       37823 $copy->{_a} = $x->{_a} if exists $x->{_a};
1275             $copy->{_p} = $x->{_p} if exists $x->{_p};
1276 17365         45885  
1277             return $copy;
1278             }
1279              
1280             sub as_int {
1281 650 50   650 1 3741 # return copy as a bigint representation of this Math::BigFloat number
1282 650 50       1612 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       1710  
1285             return $x -> copy() if $x -> isa("Math::BigInt");
1286              
1287             # disable upgrading and downgrading
1288 650         4506  
1289 650         2286 require Math::BigInt;
1290 650         1875 my $upg = Math::BigInt -> upgrade();
1291 650         1994 my $dng = Math::BigInt -> downgrade();
1292 650         1858 Math::BigInt -> upgrade(undef);
1293             Math::BigInt -> downgrade(undef);
1294 650         1038  
1295 650 100       1792 my $y;
    100          
1296 8         46 if ($x -> is_inf()) {
1297             $y = Math::BigInt -> binf($x->sign());
1298 4         43 } elsif ($x -> is_nan()) {
1299             $y = Math::BigInt -> bnan();
1300 638         2112 } else {
1301 638 100       2527 $y = $LIB->_copy($x->{_m});
    100          
1302 183         721 if ($x->{_es} eq '-') { # < 0
1303             $y = $LIB->_rsft($y, $x->{_e}, 10);
1304 14         89 } elsif (! $LIB->_is_zero($x->{_e})) { # > 0
1305             $y = $LIB->_lsft($y, $x->{_e}, 10);
1306 638         2250 }
1307             $y = Math::BigInt->new($x->{sign} . $LIB->_str($y));
1308             }
1309              
1310             # reset upgrading and downgrading
1311 650         3008  
1312 650         2129 Math::BigInt -> upgrade($upg);
1313             Math::BigInt -> downgrade($dng);
1314 650         3083  
1315             return $y;
1316             }
1317              
1318 2 50   2 1 15 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       7  
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 242334 # return true if arg (BFLOAT or num_str) is zero
1348             my (undef, $x) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
1349 110809 100 100     362217  
1350             ($x->{sign} eq '+' && $LIB->_is_zero($x->{_m})) ? 1 : 0;
1351             }
1352              
1353             sub is_one {
1354 2958 100   2958 1 10559 # 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     9188  
1357             $sign = '+' if !defined $sign || $sign ne '-';
1358              
1359             ($x->{sign} eq $sign &&
1360 2958 100 100     10996 $LIB->_is_zero($x->{_e}) &&
1361             $LIB->_is_one($x->{_m})) ? 1 : 0;
1362             }
1363              
1364             sub is_odd {
1365 104 50   104 1 883 # 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     631 ($LIB->_is_zero($x->{_e})) &&
1370             ($LIB->_is_odd($x->{_m}))) ? 1 : 0;
1371             }
1372              
1373             sub is_even {
1374 72 50   72 1 972 # 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     785 ($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 6649 # 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     17312 (($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 16233 # 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       6506  
1402             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
1403              
1404             # Handle all 'nan' cases.
1405 2913 100 100     11303  
1406             return if ($x->{sign} eq $nan) || ($y->{sign} eq $nan);
1407              
1408             # Handle all '+inf' and '-inf' cases.
1409              
1410 2855 100 100     12576 return 0 if ($x->{sign} eq '+inf' && $y->{sign} eq '+inf' ||
      100        
      100        
1411 2847 100       6082 $x->{sign} eq '-inf' && $y->{sign} eq '-inf');
1412 2791 100       5699 return +1 if $x->{sign} eq '+inf'; # x = +inf and y < +inf
1413 2727 50       5276 return -1 if $x->{sign} eq '-inf'; # x = -inf and y > -inf
1414 2727 100       5041 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     9341  
1419 2267 100 100     6626 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         4216  
1424 1961         4174 my $xz = $x->is_zero();
1425 1961 100 100     5599 my $yz = $y->is_zero();
1426 1827 100 66     4776 return 0 if $xz && $yz; # 0 <=> 0
1427 1707 100 66     6027 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         1809  
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         3431  
1438 1255         3151 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       3097  
1445             if ($mxl == $myl) {
1446              
1447             # First handle the two cases where the exponents have different signs.
1448 1134 50 66     5607  
    100 100        
1449 0         0 if ($x->{_es} eq '+' && $y->{_es} eq '-') {
1450             $cmp = +1;
1451 16         52 } 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         2965 else {
1458 1118 100       2787 $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       2444  
1465 1134 100       2520 $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         1886  
1476             my $ex;
1477             my $ey;
1478 1187 100       2413  
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       1955  
1484 1042         2725 if ($y->{_es} eq '+') {
1485 1042         2513 $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         38 else {
1493 9         49 $ex = $LIB->_copy($x->{_e});
1494 9         47 $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       427  
1503 25         111 if ($y->{_es} eq '+') {
1504 25         106 $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         357 else {
1513 111         382 $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         3393  
1521 1187         4070 $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         3587  
1526 1187 100       3167 $cmp = $LIB->_acmp($ex, $ey);
1527 1187 100       2898 $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         1740  
1534 1097         1644 my $mx = $x->{_m};
1535             my $my = $y->{_m};
1536 1097 100       2909  
    100          
1537 26         96 if ($mxl > $myl) {
1538             $my = $LIB->_lsft($LIB->_copy($my), $LIB->_new($mxl - $myl), 10);
1539 5         30 } elsif ($mxl < $myl) {
1540             $mx = $LIB->_lsft($LIB->_copy($mx), $LIB->_new($myl - $mxl), 10);
1541             }
1542 1097         2450  
1543 1097 100       2722 $cmp = $LIB->_acmp($mx, $my);
1544 1097         3937 $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 45643 # 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       20150  
1558             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
1559              
1560 8979 100 100     49781 # handle +-inf and NaN's
1561 92 100 100     653 if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/) {
1562 64 100 100     193 return if (($x->{sign} eq $nan) || ($y->{sign} eq $nan));
1563 48 100 66     123 return 0 if ($x->is_inf() && $y->is_inf());
1564 16         175 return 1 if ($x->is_inf() && !$y->is_inf());
1565             return -1;
1566             }
1567              
1568 8887         21273 # shortcut
1569 8887         16881 my $xz = $x->is_zero();
1570 8887 100 100     20484 my $yz = $y->is_zero();
1571 8883 100 66     19748 return 0 if $xz && $yz; # 0 <=> 0
1572 8843 100 66     17988 return -1 if $xz && !$yz; # 0 <=> +y
1573             return 1 if $yz && !$xz; # +x <=> 0
1574              
1575 8803         21492 # adjust so that exponents are equal
1576 8803         20099 my $lxm = $LIB->_len($x->{_m});
1577 8803         15728 my $lym = $LIB->_len($y->{_m});
1578 8803 100       19832 my ($xes, $yes) = (1, 1);
1579 8803 100       17512 $xes = -1 if $x->{_es} ne '+';
1580             $yes = -1 if $y->{_es} ne '+';
1581 8803         21118 # the numify somewhat limits our length, but makes it much faster
1582 8803         20504 my $lx = $lxm + $xes * $LIB->_num($x->{_e});
1583 8803         15268 my $ly = $lym + $yes * $LIB->_num($y->{_e});
1584 8803 100       26142 my $l = $lx - $ly;
1585             return $l <=> 0 if $l != 0;
1586              
1587             # lengths (corrected by exponent) are equal
1588 3821         5497 # so make mantissa equal-length by padding with zero (shift left)
1589 3821         5852 my $diff = $lxm - $lym;
1590 3821         5875 my $xm = $x->{_m}; # not yet copy it
1591 3821 100       9610 my $ym = $y->{_m};
    100          
1592 578         1946 if ($diff > 0) {
1593 578         2044 $ym = $LIB->_copy($y->{_m});
1594             $ym = $LIB->_lsft($ym, $LIB->_new($diff), 10);
1595 258         984 } elsif ($diff < 0) {
1596 258         1128 $xm = $LIB->_copy($x->{_m});
1597             $xm = $LIB->_lsft($xm, $LIB->_new(-$diff), 10);
1598 3821         10210 }
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 1959 # negate number or make a negated number from string
1609             my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
1610 475 50       1434  
1611             return $x if $x->modify('bneg');
1612 475 100       1198  
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     2106 $x->{sign} =~ tr/+-/-+/
1617             unless $x->{sign} eq '+' && $LIB->_is_zero($x->{_m});
1618 466 100 66     1245  
      66        
1619             return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade)
1620 463         1301 && ($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 190781 # adjust m and e so that m is smallest possible
1629             my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
1630 72362 50       145168  
1631             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
1632              
1633 72362 100       261634 # inf, nan etc
1634 6 100       42 if ($x->{sign} !~ /^[+-]$/) {
1635 4         16 return $downgrade -> new($x) if defined $downgrade;
1636             return $x;
1637             }
1638 72356         190327  
1639 72356 100       140159 my $zeros = $LIB->_zeros($x->{_m}); # correct for trailing zeros
1640 29435         71984 if ($zeros != 0) {
1641 29435         93153 my $z = $LIB->_new($zeros);
1642 29435 100       65363 $x->{_m} = $LIB->_rsft($x->{_m}, $z, 10);
1643 27727 100       74647 if ($x->{_es} eq '-') {
1644 27387         77079 if ($LIB->_acmp($x->{_e}, $z) >= 0) {
1645 27387 100       70563 $x->{_e} = $LIB->_sub($x->{_e}, $z);
1646             $x->{_es} = '+' if $LIB->_is_zero($x->{_e});
1647 340         1193 } else {
1648 340         1068 $x->{_e} = $LIB->_sub($LIB->_copy($z), $x->{_e});
1649             $x->{_es} = '+';
1650             }
1651 1708         5269 } 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       108374 # zeros). So, for something like 0Ey, set y to 0, and -0 => +0
1657 399         969 if ($LIB->_is_zero($x->{_m})) {
1658 399         746 $x->{sign} = '+';
1659 399         986 $x->{_es} = '+';
1660             $x->{_e} = $LIB->_zero();
1661             }
1662             }
1663 72356 100 100     168104  
1664             return $downgrade -> new($x)
1665 72336         242138 if defined($downgrade) && $x->is_int();
1666             return $x;
1667             }
1668              
1669             sub binc {
1670 4362 50   4362 1 12733 # increment arg by one
1671             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
1672 4362 50       11903  
1673             return $x if $x->modify('binc');
1674              
1675             # Inf and NaN
1676 4362 100       11013  
1677 4357 100       10258 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       10252  
1682 461         1530 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       9920  
1689 350         1506 if (!$LIB->_is_zero($x->{_e})) {
1690 350         1475 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # 1e2 => 100
1691 350         897 $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       8522 # now $x->{_e} == 0
    50          
1698 3871         9396 if ($x->{sign} eq '+') {
1699 3871         8660 $x->{_m} = $LIB->_inc($x->{_m});
1700             return $x->bnorm()->bround(@r);
1701 16         54 } elsif ($x->{sign} eq '-') {
1702 16 100       65 $x->{_m} = $LIB->_dec($x->{_m});
1703 16         52 $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 1303 # decrement arg by one
1714             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
1715 168 50       669  
1716             return $x if $x->modify('bdec');
1717              
1718             # Inf and NaN
1719 168 100       581  
1720 163 100       554 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       614  
1725 113         451 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       149  
1732 8         36 if (!$LIB->_is_zero($x->{_e})) {
1733 8         34 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # 1e2 => 100
1734 8         27 $x->{_e} = $LIB->_zero(); # normalize
1735             $x->{_es} = '+';
1736             }
1737              
1738 41         124 # now $x->{_e} == 0
1739 41 100 100     238 my $zero = $x->is_zero();
    50          
1740 21         83 if (($x->{sign} eq '-') || $zero) { # x <= 0
1741 21 100       101 $x->{_m} = $LIB->_inc($x->{_m});
1742 21 50       61 $x->{sign} = '-' if $zero; # 0 => 1 => -1
1743 21         80 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # -1 +1 => -0 => +0
1744             return $x->bnorm()->round(@r);
1745             }
1746 20         81 elsif ($x->{sign} eq '+') { # x > 0
1747 20         67 $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 62633 # 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       39540  
1762             return $x if $x->modify('badd');
1763              
1764 13458 100 100     71385 # inf and NaN handling
1765             if ($x->{sign} !~ /^[+-]$/ || $y->{sign} !~ /^[+-]$/) {
1766              
1767 226 100 100     1610 # $x is NaN and/or $y is NaN
    100 100        
    100          
1768 110         289 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       222 # +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         87 elsif ($y->{sign} =~ /^[+-]inf$/) {
1780             $x->{sign} = $y->{sign};
1781             }
1782 226 100       552  
1783 218         686 return $downgrade -> new($x -> bdstr(), @r) if defined $downgrade;
1784             return $x -> round(@r);
1785             }
1786 13232 100       27070  
1787             return $upgrade->badd($x, $y, @r) if defined $upgrade;
1788 13231         21938  
1789             $r[3] = $y; # no push!
1790              
1791 13231 100       26428 # for speed: no add for $x + 0
    100          
1792 24         112 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         555 # make copy, clobbering up x (modify in place!)
1799 102         340 $x->{_e} = $LIB->_copy($y->{_e});
1800 102         349 $x->{_es} = $y->{_es};
1801 102   33     425 $x->{_m} = $LIB->_copy($y->{_m});
1802 102         349 $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         22926 # take lower of the two e's and adapt m1 to it to match m2
1810 13105 50       25583 my $e = $y->{_e};
1811 13105         30535 $e = $LIB->_zero() if !defined $e; # if no BFLOAT?
1812             $e = $LIB->_copy($e); # make copy (didn't do it yet)
1813 13105         20191  
1814             my $es;
1815 13105   50     48593  
1816             ($e, $es) = $LIB -> _ssub($e, $y->{_es} || '+', $x->{_e}, $x->{_es});
1817 13105         37095  
1818             my $add = $LIB->_copy($y->{_m});
1819 13105 100       31864  
    100          
1820 7475         21115 if ($es eq '-') { # < 0
1821 7475         24444 $x->{_m} = $LIB->_lsft($x->{_m}, $e, 10);
1822             ($x->{_e}, $x->{_es}) = $LIB -> _sadd($x->{_e}, $x->{_es}, $e, $es);
1823 789         2466 } 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       32507  
1829 11511         28659 if ($x->{sign} eq $y->{sign}) {
1830             $x->{_m} = $LIB->_add($x->{_m}, $add);
1831             } else {
1832 1594         4750 ($x->{_m}, $x->{sign}) =
1833             $LIB -> _sadd($x->{_m}, $x->{sign}, $add, $y->{sign});
1834             }
1835              
1836 13105         31493 # delete trailing zeros, then round
1837             $x = $x->bnorm()->round(@r);
1838             }
1839 13231 100 66     34619  
1840             return $downgrade -> new($x -> bdstr(), @r)
1841 13224         44950 if defined($downgrade) && $x -> is_int();
1842             return $x; # rounding already done above
1843             }
1844              
1845             sub bsub {
1846 1655 100 100 1655 1 10556 # 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       5175  
1851             return $x if $x -> modify('bsub');
1852 1655 100       3323  
1853 42         181 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         3351  
1861 1613         3284 my $xsign = $x -> {sign};
1862 1613 100       3881 $y -> {sign} =~ tr/+-/-+/; # does nothing for NaN
1863             if ($xsign ne $x -> {sign}) {
1864 24 100       108 # special case of $x -> bsub($x) results in 0
1865 16         53 if ($xsign =~ /^[+-]$/) {
1866             $x = $x -> bzero(@r);
1867 8         27 } else {
1868             $x = $x -> bnan(); # NaN, -inf, +inf
1869 24 50       59 }
1870 24         78 return $downgrade -> new($x -> bdstr(), @r) if defined $downgrade;
1871             return $x -> round(@r);
1872 1589         4086 }
1873 1589         3930 $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     3531  
      66        
1877             return $downgrade -> new($x -> bdstr(), @r)
1878 1623         5904 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 91784 # 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       55328  
1890             return $x if $x->modify('bmul');
1891 19512 100 100     69326  
1892             return $x->bnan(@r) if ($x->{sign} eq $nan) || ($y->{sign} eq $nan);
1893              
1894 19455 100 100     64923 # inf handling
1895 85 100 100     245 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     431 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1900 50 100 100     313 return $x->binf(@r) if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1901 36         108 return $x->binf(@r) if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1902             return $x->binf('-', @r);
1903             }
1904 19370 100       36185  
1905             return $upgrade->bmul($x, $y, @r) if defined $upgrade;
1906              
1907 19369         55765 # aEb * cEd = (a*c)E(b+d)
1908             $x->{_m} = $LIB->_mul($x->{_m}, $y->{_m});
1909 19369         67059 ($x->{_e}, $x->{_es})
1910             = $LIB -> _sadd($x->{_e}, $x->{_es}, $y->{_e}, $y->{_es});
1911 19369         35717  
1912             $r[3] = $y; # no push!
1913              
1914 19369 100       44408 # adjust sign:
1915 19369         43164 $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1916             $x = $x->bnorm->round(@r);
1917 19369 0 33     40514  
      66        
1918             return $downgrade -> new($x -> bdstr(), @r)
1919 19365         48819 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 4160 # 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       819  
1932             return $x if $x->modify('bmuladd');
1933              
1934             return $x->bnan(@r) if (($x->{sign} eq $nan) ||
1935 247 100 100     1363 ($y->{sign} eq $nan) ||
      100        
1936             ($z->{sign} eq $nan));
1937              
1938 214 100 66     902 # inf handling
1939 18 50 33     56 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     268 # +inf * -/-inf => -inf, -inf * +/+inf => -inf
1944 12 100 100     151 return $x->binf(@r) if ($x->{sign} =~ /^\+/ && $y->{sign} =~ /^\+/);
1945 8         38 return $x->binf(@r) if ($x->{sign} =~ /^-/ && $y->{sign} =~ /^-/);
1946             return $x->binf('-', @r);
1947             }
1948              
1949 196         712 # aEb * cEd = (a*c)E(b+d)
1950             $x->{_m} = $LIB->_mul($x->{_m}, $y->{_m});
1951 196         769 ($x->{_e}, $x->{_es})
1952             = $LIB -> _sadd($x->{_e}, $x->{_es}, $y->{_e}, $y->{_es});
1953 196         401  
1954             $r[3] = $y; # no push!
1955              
1956 196 100       495 # adjust sign:
1957             $x->{sign} = $x->{sign} ne $y->{sign} ? '-' : '+';
1958              
1959 196 100       460 # z=inf handling (z=NaN handled above)
1960 1         4 if ($z->{sign} =~ /^[+-]inf$/) {
1961 1 50       8 $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         361 # take lower of the two e's and adapt m1 to it to match m2
1967 195 50       397 my $e = $z->{_e};
1968 195         480 $e = $LIB->_zero() if !defined $e; # if no BFLOAT?
1969             $e = $LIB->_copy($e); # make copy (didn't do it yet)
1970 195         299  
1971             my $es;
1972 195   50     694  
1973             ($e, $es) = $LIB -> _ssub($e, $z->{_es} || '+', $x->{_e}, $x->{_es});
1974 195         583  
1975             my $add = $LIB->_copy($z->{_m});
1976 195 100       592  
    100          
1977             if ($es eq '-') # < 0
1978 4         47 {
1979 4         41 $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         86 {
1983             $add = $LIB->_lsft($add, $e, 10);
1984             }
1985             # else: both e are the same, so just leave them
1986 195 100       459  
1987             if ($x->{sign} eq $z->{sign}) {
1988 151         387 # add
1989             $x->{_m} = $LIB->_add($x->{_m}, $add);
1990             } else {
1991 44         166 ($x->{_m}, $x->{sign}) =
1992             $LIB -> _sadd($x->{_m}, $x->{sign}, $add, $z->{sign});
1993             }
1994              
1995 195         542 # delete trailing zeros, then round
1996             $x = $x->bnorm()->round(@r);
1997 195 0 33     450  
      66        
1998             return $downgrade -> new($x -> bdstr(), @r)
1999 192         2430 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 29316 # set up parameters
2008             my ($class, $x, $y, @r) = (ref($_[0]), @_);
2009 9074 100 100     36291 # objectify is costly, so avoid it
2010 137         464 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) {
2011             ($class, $x, $y, @r) = objectify(2, @_);
2012             }
2013 9074 50       24218  
2014             return $x if $x->modify('bdiv');
2015 9074         15283  
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     22615  
2021 68 100       270 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       19571  
2030 55         134 if ($y -> is_zero()) {
2031 55 100       151 my ($quo, $rem);
2032 16         78 if ($wantarray) {
2033 16 50 33     85 $rem = $x -> copy() -> round(@r);
2034             $rem = $downgrade -> new($rem, @r)
2035             if defined($downgrade) && $rem -> is_int();
2036 55 100       119 }
2037 21         96 if ($x -> is_zero()) {
2038             $quo = $x -> bnan(@r);
2039 34         142 } else {
2040             $quo = $x -> binf($x -> {sign}, @r);
2041 52 100       482 }
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       21512  
2049 40         129 if ($x -> is_inf()) {
2050 40 100       168 my ($quo, $rem);
2051 40 100       138 $rem = $class -> bnan(@r) if $wantarray;
2052 16         55 if ($y -> is_inf()) {
2053             $quo = $x -> bnan(@r);
2054 24 100       105 } else {
2055 24         101 my $sign = $x -> bcmp(0) == $y -> bcmp(0) ? '+' : '-';
2056             $quo = $x -> binf($sign, @r);
2057 40 100       208 }
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       18220  
2067 40         87 if ($y -> is_inf()) {
2068 40 100       120 my ($quo, $rem);
2069 16 100 100     35 if ($wantarray) {
2070 12         48 if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) {
2071 12 50 33     50 $rem = $x -> copy() -> round(@r);
2072             $rem = $downgrade -> new($rem, @r)
2073 12         43 if defined($downgrade) && $rem -> is_int();
2074             $quo = $x -> bzero(@r);
2075 4         15 } else {
2076 4         21 $rem = $class -> binf($y -> {sign}, @r);
2077             $quo = $x -> bone('-', @r);
2078 16         67 }
2079             return ($quo, $rem);
2080 24 50       65 } else {
2081 24 50 33     69 if ($y -> is_inf()) {
2082 0         0 if ($x -> is_nan() || $x -> is_inf()) {
2083             return $x -> bnan(@r);
2084 24         78 } 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       16598 # x == 0?
2095 58         171 if ($x->is_zero()) {
2096 58         245 my ($quo, $rem);
2097 58 100 66     330 $quo = $x->round(@r);
2098             $quo = $downgrade -> new($quo, @r)
2099 58 100       180 if defined($downgrade) && $quo -> is_int();
2100 13         58 if ($wantarray) {
2101 13         69 $rem = $class -> bzero(@r);
2102             return $quo, $rem;
2103 45         312 }
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     20343 return $upgrade -> bdiv($x, $y, @r)
      33        
2111             if defined($upgrade) && !wantarray && !$LIB -> _is_one($y -> {_m});
2112              
2113 8813         12547 # we need to limit the accuracy to protect against overflow
2114 8813         12922 my $fallback = 0;
2115 8813         34031 my (@params, $scale);
2116             ($x, @params) = $x->_find_round_parameters($r[0], $r[1], $r[2], $y);
2117 8813 50       29215  
2118             return $x -> round(@r) if $x->is_nan(); # error in _find_round_parameters?
2119              
2120 8813 100       19421 # no rounding at all, so must use fallback
2121             if (scalar @params == 0) {
2122 506         1641 # simulate old behaviour
2123 506         963 $params[0] = $class->div_scale(); # and round to it as accuracy
2124 506         789 $scale = $params[0]+4; # at least four more for proper round
2125 506         833 $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     19741 # enough...
2130             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2131             }
2132 8813         12779  
2133 8813 100       17249 my $rem;
2134             $rem = $class -> bzero() if wantarray;
2135 8813 50       19524  
2136             $y = $class->new($y) unless $y->isa('Math::BigFloat');
2137 8813         29218  
2138 8813         20404 my $lx = $LIB -> _len($x->{_m});
2139 8813 100       19002 my $ly = $LIB -> _len($y->{_m});
2140 8813 100       16272 $scale = $lx if $lx > $scale;
2141 8813         12441 $scale = $ly if $ly > $scale;
2142 8813 100       16479 my $diff = $ly - $lx;
2143             $scale += $diff if $diff > 0; # if lx << ly, but not if ly << lx!
2144              
2145 8813   100     20774 # 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         16932 # case of $x->bsub($x); so we can catch it below:
2150 8813         18520 my $xsign = $x->{sign};
2151             $y->{sign} =~ tr/+-/-+/;
2152 8813 100       19120  
2153             if ($xsign ne $x->{sign}) {
2154 8         46 # 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         13251 # 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     19373 # calculation
2163 25         70 if (wantarray && $y_not_one) {
2164             $rem = $x->copy();
2165             }
2166 8805 100       23461  
2167             $x->{sign} = $x->{sign} ne $y->sign() ? '-' : '+';
2168              
2169 8805 100       19613 # check for / +-1 (+/- 1E0)
2170             if ($y_not_one) {
2171             # promote Math::BigInt and its subclasses (except when already a
2172 8636 50       16068 # 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         28055 # a * 10 ** b / c * 10 ** d => a/c * 10 ** (b-d)
2177 8636         30425 $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         31952 ($x->{_e}, $x->{_es})
2182             = $LIB -> _ssub($x->{_e}, $x->{_es}, $y->{_e}, $y->{_es});
2183             # correct for 10**scale
2184 8636         27317 ($x->{_e}, $x->{_es})
2185 8636         26428 = $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       18001 # shortcut to not run through _find_round_parameters again
2191 8780         14542 if (defined $params[0]) {
2192 8780         22575 $x->{_a} = undef; # clear before round
2193             $x = $x->bround($params[0], $params[2]); # then round accordingly
2194 33         73 } else {
2195 33         112 $x->{_p} = undef; # clear before round
2196             $x = $x->bfround($params[1], $params[2]); # then round accordingly
2197 8813 100       20200 }
2198             if ($fallback) {
2199 506         1400 # clear a/p after round, since user did not request it
2200 506         896 $x->{_a} = undef;
2201             $x->{_p} = undef;
2202             }
2203 8813 100       16918  
2204 49 100       135 if (wantarray) {
2205 25         112 if ($y_not_one) {
2206 25         113 $x = $x -> bfloor();
2207             $rem = $rem->bmod($y, @params); # copy already done
2208 49 50       120 }
2209             if ($fallback) {
2210 49         96 # clear a/p after round, since user did not request it
2211 49         84 $rem->{_a} = undef;
2212             $rem->{_p} = undef;
2213 49 50 33     152 }
2214             $x = $downgrade -> new($x -> bdstr(), @r)
2215 49 50 33     167 if defined($downgrade) && $x -> is_int();
2216             $rem = $downgrade -> new($rem -> bdstr(), @r)
2217 49         274 if defined($downgrade) && $rem -> is_int();
2218             return ($x, $rem);
2219             }
2220 8764 50 66     17871  
2221             $x = $downgrade -> new($x, @r)
2222 8764         30547 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 8531 # 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       2477  
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     2038  
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       1706  
2243 40         138 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       1726  
2250 48         179 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       1406  
2257 40 100 100     148 if ($y -> is_inf()) {
2258 28         118 if ($x -> is_zero() || $x -> bcmp(0) == $y -> bcmp(0)) {
2259             return $x -> round(@r);
2260 12         68 } 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     1093 # check that $y == +1 or $y == -1:
      100        
      100        
2268             ($LIB->_is_zero($y->{_e}) && $LIB->_is_one($y->{_m})));
2269 483         1438  
2270 483 100       1162 my $cmp = $x->bacmp($y); # equal or $x < $y?
2271 30         153 if ($cmp == 0) { # $x == $y => result 0
2272             return $x -> bzero(@r);
2273             }
2274              
2275 453 100       1124 # only $y of the operands negative?
2276             my $neg = $x->{sign} ne $y->{sign} ? 1 : 0;
2277 453         824  
2278 453 100 100     1233 $x->{sign} = $y->{sign}; # calc sign first
2279 48         184 if ($cmp < 0 && $neg == 0) { # $x < $y => result $x
2280             return $x -> round(@r);
2281             }
2282 405         1072  
2283             my $ym = $LIB->_copy($y->{_m});
2284              
2285             # 2e1 => 20
2286 405 100 100     1497 $ym = $LIB->_lsft($ym, $y->{_e}, 10)
2287             if $y->{_es} eq '+' && !$LIB->_is_zero($y->{_e});
2288              
2289 405         687 # if $y has digits after dot
2290 405 100       883 my $shifty = 0; # correct _e of $x by this
2291             if ($y->{_es} eq '-') # has digits after dot
2292             {
2293 16         74 # 123 % 2.5 => 1230 % 25 => 5 => 0.5
2294             $shifty = $LIB->_num($y->{_e}); # no more digits after dot
2295 16         91 # 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         624  
2300 405 100       877 my $shiftx = 0; # correct _e of $x by this
2301             if ($x->{_es} eq '-') # has digits after dot
2302             {
2303 27         117 # 123.4 % 20 => 1234 % 200
2304 27         119 $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     1433 # 123e1 % 20 => 1230 % 20
2308 117         461 if ($x->{_es} eq '+' && !$LIB->_is_zero($x->{_e})) {
2309             $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # es => '+' here
2310             }
2311 405         1129  
2312 405         846 $x->{_e} = $LIB->_new($shiftx);
2313 405 100 100     1486 $x->{_es} = '+';
2314 405 100       871 $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         1246  
2319             $x->{_m} = $LIB->_mod($x->{_m}, $ym);
2320 405 100       1016  
2321 405         1042 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # fix sign for -0
2322             $x = $x->bnorm();
2323              
2324 405 100 66     1182 # if one of them negative => correct in place
2325 51         260 if ($neg != 0 && ! $x -> is_zero()) {
2326 51         107 my $r = $y - $x;
2327 51         122 $x->{_m} = $r->{_m};
2328 51         101 $x->{_e} = $r->{_e};
2329 51 50       166 $x->{_es} = $r->{_es};
2330 51         131 $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # fix sign for -0
2331             $x = $x->bnorm();
2332             }
2333 405         1701  
2334 405 0 0     1300 $x = $x->round($r[0], $r[1], $r[2], $y);
      33        
2335             return $downgrade -> new($x -> bdstr(), @r)
2336 405         3440 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 442 # 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       87  
2349             return $num if $num->modify('bmodpow');
2350 20 50 33     79  
      33        
2351             return $num -> bnan(@r)
2352             if $mod->is_nan() || $exp->is_nan() || $mod->is_nan();
2353              
2354 20 50 33     109 # check modulus for valid values
2355             return $num->bnan(@r) if $mod->{sign} ne '+' || $mod->is_zero();
2356              
2357 20 50       98 # 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       80  
2363             $num = $num->bmodinv($mod, @r) if $exp->{sign} eq '-';
2364              
2365 20 50       80 # 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         76 # XXX TODO: speed it up when all three numbers are integers
2371             $num = $num->bpow($exp)->bmod($mod);
2372 20 0 0     80  
      33        
2373             return $downgrade -> new($num -> bdstr(), @r) if defined($downgrade)
2374 20         64 && ($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 12681 # set up parameters
2384             my ($class, $x, $y, $a, $p, $r) = (ref($_[0]), @_);
2385 1042 100 100     4627 # objectify is costly, so avoid it
2386 110         398 if ((!ref($_[0])) || (ref($_[0]) ne ref($_[1]))) {
2387             ($class, $x, $y, $a, $p, $r) = objectify(2, @_);
2388             }
2389 1042 50       3418  
2390             return $x if $x -> modify('bpow');
2391              
2392 1042 100 100     2781 # $x and/or $y is a NaN
2393             return $x -> bnan() if $x -> is_nan() || $y -> is_nan();
2394              
2395 926 100       2608 # $x and/or $y is a +/-Inf
    100          
    100          
    100          
2396 60 100       211 if ($x -> is_inf("-")) {
2397 32 100       136 return $x -> bzero() if $y -> is_negative();
2398 28 100       117 return $x -> bnan() if $y -> is_zero();
2399 20         92 return $x if $y -> is_odd();
2400             return $x -> bneg();
2401 60 100       212 } elsif ($x -> is_inf("+")) {
2402 32 100       101 return $x -> bzero() if $y -> is_negative();
2403 28         349 return $x -> bnan() if $y -> is_zero();
2404             return $x;
2405 44 100       199 } elsif ($y -> is_inf("-")) {
2406 40 100 100     172 return $x -> bnan() if $x -> is_one("-");
2407 28 100       118 return $x -> binf("+") if $x > -1 && $x < 1;
2408 24         144 return $x -> bone() if $x -> is_one("+");
2409             return $x -> bzero();
2410 44 100       262 } elsif ($y -> is_inf("+")) {
2411 40 100 100     220 return $x -> bnan() if $x -> is_one("-");
2412 28 100       127 return $x -> bzero() if $x > -1 && $x < 1;
2413 24         129 return $x -> bone() if $x -> is_one("+");
2414             return $x -> binf("+");
2415             }
2416 718 100       2708  
2417 48 100       134 if ($x -> is_zero()) {
2418 44 100       158 return $x -> bone() if $y -> is_zero();
2419 24         274 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     2455  
2425 80 50       209 if ($x -> is_negative() && !$y -> is_int()) {
2426 80         201 return $upgrade -> bpow($x, $y, $a, $p, $r) if defined $upgrade;
2427             return $x -> bnan();
2428             }
2429 590 100 100     1529  
2430 119         1313 if ($x -> is_one("+") || $y -> is_one()) {
2431             return $x;
2432             }
2433 471 100       1196  
2434 24 100       81 if ($x -> is_one("-")) {
2435 12         55 return $x if $y -> is_odd();
2436             return $x -> bneg();
2437             }
2438 447 100       1335  
2439             return $x -> _pow($y, $a, $p, $r) if !$y -> is_int();
2440 333         1074  
2441             my $y1 = $y -> as_int()->{value}; # make MBI part
2442 333         1009  
2443 333 100       1203 my $new_sign = '+';
    100          
2444             $new_sign = $LIB -> _is_odd($y1) ? '-' : '+' if $x->{sign} ne '+';
2445              
2446 333         1271 # calculate $x->{_m} ** $y and $x->{_e} * $y separately (faster)
2447 333         1091 $x->{_m} = $LIB -> _pow($x->{_m}, $y1);
2448             $x->{_e} = $LIB -> _mul($x->{_e}, $y1);
2449 333         854  
2450 333         1001 $x->{sign} = $new_sign;
2451             $x = $x -> bnorm();
2452              
2453             # x ** (-y) = 1 / (x ** y)
2454 333 100       1264  
2455             if ($y->{sign} eq '-') {
2456 105         373 # modify $x in place!
2457 105         369 my $z = $x -> copy();
2458             $x = $x -> bone();
2459 105         415 # round in one go (might ignore y's A!)
2460             return scalar $x -> bdiv($z, $a, $p, $r);
2461             }
2462 228         895  
2463             $x = $x -> round($a, $p, $r, $y);
2464 228 50 66     805  
      66        
2465             return $downgrade -> new($x)
2466 226         2456 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 1701  
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     1258  
2481             if (!ref($_[0]) && $_[0] =~ /^[A-Za-z]|::/) {
2482 18 100       103 # E.g., Math::BigFloat->blog(256, 2)
2483             ($class, $x, $base, @r) =
2484             defined $_[2] ? objectify(2, @_) : objectify(1, @_);
2485             } else {
2486 238 100       1105 # 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       1282  
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       780  
2496             return $x -> bnan(@r) if $x -> is_nan();
2497 252 100       946  
2498 42 50 33     442 if (defined $base) {
2499             $base = $class -> new($base)
2500 42 100 66     152 unless defined(blessed($base)) && $base -> isa(__PACKAGE__);
    100 66        
    100          
2501 8         49 if ($base -> is_nan() || $base -> is_one()) {
2502             return $x -> bnan(@r);
2503 4 50 33     31 } elsif ($base -> is_inf() || $base -> is_zero()) {
2504 4         40 return $x -> bnan(@r) if $x -> is_inf() || $x -> is_zero();
2505             return $x -> bzero(@r);
2506 4 50       23 } elsif ($base -> is_negative()) { # -inf < base < 0
2507 4 50       17 return $x -> bzero(@r) if $x -> is_one(); # x = 1
2508             return $x -> bone('+', @r) if $x == $base; # x = base
2509 4 50       23 # we can't handle these cases, so upgrade, if we can
2510 4         33 return $upgrade -> blog($x, $base, @r) if defined $upgrade;
2511             return $x -> bnan(@r);
2512 26 100       155 }
2513             return $x -> bone(@r) if $x == $base; # 0 < base && 0 < x < inf
2514             }
2515 225 100       671  
    100          
    100          
    100          
2516 8 50 33     49 if ($x -> is_inf()) { # x = +/-inf
2517 8         28 my $sign = defined($base) && $base < 1 ? '-' : '+';
2518             return $x -> binf($sign, @r);
2519 16 50       66 } elsif ($x -> is_neg()) { # -inf < x < 0
2520 16         62 return $upgrade -> blog($x, $base, @r) if defined $upgrade;
2521             return $x -> bnan(@r);
2522 16         127 } elsif ($x -> is_one()) { # x = 1
2523             return $x -> bzero(@r);
2524 8 50 33     72 } elsif ($x -> is_zero()) { # x = 0
2525 8         32 my $sign = defined($base) && $base < 1 ? '+' : '-';
2526             return $x -> binf($sign, @r);
2527             }
2528              
2529 177         562 # we need to limit the accuracy to protect against overflow
2530 177         364 my $fallback = 0;
2531 177         740 my ($scale, @params);
2532             ($x, @params) = $x->_find_round_parameters(@r);
2533              
2534 177 100       708 # no rounding at all, so must use fallback
2535             if (scalar @params == 0) {
2536 91         374 # simulate old behaviour
2537 91         202 $params[0] = $class->div_scale(); # and round to it as accuracy
2538 91         174 $params[1] = undef; # P = undef
2539 91         163 $scale = $params[0]+4; # at least four more for proper round
2540 91         205 $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     431 # 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   478 # disable them and later re-enable them
  43         120  
  43         26307  
2550 177         519 no strict 'refs';
2551 177         464 my $abr = "$class\::accuracy";
2552 177         402 my $ab = $$abr;
2553 177         414 $$abr = undef;
2554 177         374 my $pbr = "$class\::precision";
2555 177         422 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         393 # them already into account), since these would interfere, too
2559 177         319 $x->{_a} = undef;
2560             $x->{_p} = undef;
2561 177         331  
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     583  
      66        
2567 11         47 if (defined($base) && $base -> is_int() && $x -> is_int()) {
2568 11         51 my $x_lib = $LIB -> _new($x -> bdstr());
2569 11         72 my $b_lib = $LIB -> _new($base -> bdstr());
2570 11 100       44 ($x_lib, my $exact) = $LIB -> _log_int($x_lib, $b_lib);
2571 10         33 if ($exact) {
2572 10         36 $x->{_m} = $x_lib;
2573 10         51 $x->{_e} = $LIB -> _zero();
2574 10         29 $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       498  
2583 167         651 unless ($done) {
2584 167 100       639 $x = $x -> _log_10($scale);
2585             if (defined $base) {
2586 1         5 # log_b(x) = ln(x) / ln(b), so compute ln(b)
2587 1         4 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       527  
2594 177         585 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       1035 }
2599             if ($fallback) {
2600 91         242 # clear a/p after round, since user did not request it
2601 91         190 $x->{_a} = undef;
2602             $x->{_p} = undef;
2603             }
2604 177         680 # restore globals
2605 177         491 $$abr = $ab;
2606             $$pbr = $pb;
2607 177 100 100     695  
2608             return $downgrade -> new($x -> bdstr(), @r)
2609 176         2082 if defined($downgrade) && $x -> is_int();
2610             return $x;
2611             }
2612              
2613             sub bexp {
2614 20 100   20 1 115 # 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       96  
2617             return $x if $x->modify('bexp');
2618 20 50       78  
2619 20 50       65 return $x->bnan(@r) if $x -> is_nan();
2620 20 50       58 return $x->binf(@r) if $x->{sign} eq '+inf';
2621             return $x->bzero(@r) if $x->{sign} eq '-inf';
2622              
2623 20         34 # we need to limit the accuracy to protect against overflow
2624 20         40 my $fallback = 0;
2625 20         78 my ($scale, @params);
2626             ($x, @params) = $x->_find_round_parameters(@r);
2627              
2628 20 50       918 # error in _find_round_parameters?
2629             return $x->bnan(@r) if $x->{sign} eq 'NaN';
2630              
2631 20 100       83 # no rounding at all, so must use fallback
2632             if (scalar @params == 0) {
2633 11         68 # simulate old behaviour
2634 11         29 $params[0] = $class->div_scale(); # and round to it as accuracy
2635 11         19 $params[1] = undef; # P = undef
2636 11         19 $scale = $params[0]+4; # at least four more for proper round
2637 11         17 $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     39 # enough ...
2642             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
2643             }
2644 20 50       55  
2645             return $x->bone(@params) if $x->is_zero();
2646 20 50       79  
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   423 # disable them and later re-enable them
  43         121  
  43         43644  
2654 20         67 no strict 'refs';
2655 20         57 my $abr = "$class\::accuracy";
2656 20         49 my $ab = $$abr;
2657 20         48 $$abr = undef;
2658 20         52 my $pbr = "$class\::precision";
2659 20         44 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         44 # them already into account), since these would interfere, too
2663 20         43 $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         69  
2672 20         89 my $dng = Math::BigFloat -> downgrade();
2673             Math::BigFloat -> downgrade(undef);
2674 20         92  
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       75  
2724             if ($scale <= 75) {
2725 17         125 # set $x directly from a cached string form
2726             $x->{_m} = $LIB->_new("2718281828459045235360287471352662497757" .
2727 17         45 "2470936999595749669676277240766303535476");
2728 17         37 $x->{sign} = '+';
2729 17         52 $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         16 "0164779391644753259799242");
2738 3         7 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         17 # big
2743             my $steps = _len_to_steps($scale - 4);
2744 3         10 # print STDERR "# Doing $steps steps for ", $scale-4, " digits\n";
2745             while ($step++ <= $steps) {
2746 79         152 # calculate $a * $f + 1
2747 79         197 $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         9 # manually)
2755             my $B = $LIB->_fac($LIB->_new($steps));
2756              
2757             # print "A ", $LIB->_str($A), "\nB ", $LIB->_str($B), "\n";
2758              
2759 3         12 # compute A/B with $scale digits in the result (truncate, not round)
2760 3         10 $A = $LIB->_lsft($A, $LIB->_new($scale), 10);
2761             $A = $LIB->_div($A, $B);
2762 3         10  
2763 3         9 $x->{_m} = $A;
2764 3         7 $x->{sign} = '+';
2765 3         36 $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       69 # round
2771             if (!$x_org->is_one()) {
2772 10         31 # Reduce size of fractional part, followup with integer power of two.
2773 10   66     59 my $lshift = 0;
2774 21         73 while ($lshift < 30 && $x_org->bacmp(2 << $lshift) > 0) {
2775             $lshift++;
2776             }
2777 10 100       42 # Raise $x to the wanted power and round it.
2778 5         34 if ($lshift == 0) {
2779             $x = $x->bpow($x_org, @params);
2780 5         17 } else {
2781 5         21 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         32 # else just round the already computed result
2787 10         23 $x->{_a} = undef;
2788             $x->{_p} = undef;
2789 10 50       44 # 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       212  
2797             if ($fallback) {
2798 11         29 # clear a/p after round, since user did not request it
2799 11         23 $x->{_a} = undef;
2800             $x->{_p} = undef;
2801             }
2802              
2803 20         83 # Restore globals
2804 20         46 $$abr = $ab;
2805             $$pbr = $pb;
2806              
2807             # Restore downgrading.
2808 20         78  
2809             Math::BigFloat -> downgrade($dng);
2810 20 50 33     249  
2811             return $downgrade -> new($x -> bdstr(), @r)
2812 20         204 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 993 # 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       159  
2823             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
2824 60 50       197  
2825             return $x if $x->modify('bnok');
2826 60 100 100     166  
2827 48 50 66     197 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         165  
2831 48         135 my $xint = Math::BigInt -> new($x -> bsstr());
2832 48         171 my $yint = Math::BigInt -> new($y -> bsstr());
2833             $xint = $xint -> bnok($yint);
2834 48 50       131  
2835             return $xint if defined $downgrade;
2836 48         148  
2837             my $xflt = Math::BigFloat -> new($xint);
2838 48         131  
2839 48         107 $x->{_m} = $xflt->{_m};
2840 48         84 $x->{_e} = $xflt->{_e};
2841 48         81 $x->{_es} = $xflt->{_es};
2842             $x->{sign} = $xflt->{sign};
2843 48         654  
2844             return $x;
2845             }
2846              
2847             sub bsin {
2848 40 50   40 1 550 # 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       142  
2855             return $x if $x->modify('bsin');
2856 40 100       95  
2857 32 100 100     95 return $x -> bzero(@r) if $x->is_zero();
2858             return $x -> bnan(@r) if $x->is_nan() || $x->is_inf();
2859              
2860 20         43 # we need to limit the accuracy to protect against overflow
2861 20         42 my $fallback = 0;
2862 20         60 my ($scale, @params);
2863             ($x, @params) = $x->_find_round_parameters(@r);
2864              
2865 20 50       76 # error in _find_round_parameters?
2866             return $x->bnan(@r) if $x->is_nan();
2867              
2868 20 50       56 # 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     54 # 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   405 # disable them and later re-enable them
  43         144  
  43         24751  
2884 20         65 no strict 'refs';
2885 20         48 my $abr = "$class\::accuracy";
2886 20         43 my $ab = $$abr;
2887 20         50 $$abr = undef;
2888 20         39 my $pbr = "$class\::precision";
2889 20         39 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         41 $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         34  
2900 20         31 local $Math::BigInt::upgrade = undef;
2901             local $Math::BigFloat::downgrade = undef;
2902 20         55  
2903 20         64 my $over = $x * $x; # X ^ 2
2904 20         55 my $x2 = $over->copy(); # X ^ 2; difference between terms
2905 20         66 $over = $over->bmul($x); # X ^ 3 as starting value
2906 20         55 my $sign = 1; # start with -=
2907 20         88 my $below = $class->new(6);
2908 20         78 my $factorial = $class->new(4);
2909 20         42 $x->{_a} = undef;
2910             $x->{_p} = undef;
2911 20         91  
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         508 # anymore, so we stop:
2917 196 100       516 my $next = $over->copy()->bdiv($below, $scale);
2918             last if $next->bacmp($limit) <= 0;
2919 176 100       344  
2920 80         193 if ($sign == 0) {
2921             $x = $x->badd($next);
2922 96         224 } else {
2923             $x = $x->bsub($next);
2924 176         301 }
2925             $sign = 1-$sign; # alternate
2926 176         439 # calculate things for the next term
2927 176         377 $over = $over->bmul($x2); # $x*$x
2928 176         428 $below = $below->bmul($factorial); # n*(n+1)
2929 176         389 $factorial = $factorial->binc();
2930 176         375 $below = $below -> bmul($factorial); # n*(n+1)
2931             $factorial = $factorial->binc();
2932             }
2933              
2934 20 50       63 # shortcut to not run through _find_round_parameters again
2935 20         68 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       84 }
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         68 # restore globals
2946 20         67 $$abr = $ab;
2947             $$pbr = $pb;
2948 20 50 33     100  
2949             return $downgrade -> new($x -> bdstr(), @r)
2950 20         395 if defined($downgrade) && $x -> is_int();
2951             $x;
2952             }
2953              
2954             sub bcos {
2955 36 50   36 1 575 # 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         64 # we need to limit the accuracy to protect against overflow
2963 36         61 my $fallback = 0;
2964 36         102 my ($scale, @params);
2965             ($x, @params) = $x->_find_round_parameters(@r);
2966              
2967 36 100 66     172 # constant object or error in _find_round_parameters?
2968 32 100       92 return $x if $x->modify('bcos') || $x->is_nan();
2969 24 100       100 return $x->bnan() if $x->is_inf();
2970             return $x->bone(@r) if $x->is_zero();
2971              
2972 16 50       70 # 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     71 # 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   424 # disable them and later re-enable them
  43         110  
  43         36994  
2988 16         42 no strict 'refs';
2989 16         43 my $abr = "$class\::accuracy";
2990 16         36 my $ab = $$abr;
2991 16         47 $$abr = undef;
2992 16         38 my $pbr = "$class\::precision";
2993 16         34 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         44 # them already into account), since these would interfere, too
2997 16         35 $x->{_a} = undef;
2998             $x->{_p} = undef;
2999 16         71  
3000 16         68 my $over = $x * $x; # X ^ 2
3001 16         54 my $x2 = $over->copy(); # X ^ 2; difference between terms
3002 16         54 my $sign = 1; # start with -=
3003 16         74 my $below = $class->new(2);
3004 16         115 my $factorial = $class->new(3);
3005 16         39 $x = $x->bone();
3006 16         28 $x->{_a} = undef;
3007             $x->{_p} = undef;
3008 16         64  
3009             my $limit = $class->new("1E-". ($scale-1));
3010 16         75 #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         386 # anymore, so we stop:
3015 156 100       445 my $next = $over->copy()->bdiv($below, $scale);
3016             last if $next->bacmp($limit) <= 0;
3017 140 100       300  
3018 68         158 if ($sign == 0) {
3019             $x = $x->badd($next);
3020 72         177 } else {
3021             $x = $x->bsub($next);
3022 140         214 }
3023             $sign = 1-$sign; # alternate
3024 140         331 # calculate things for the next term
3025 140         318 $over = $over->bmul($x2); # $x*$x
3026 140         336 $below = $below->bmul($factorial); # n*(n+1)
3027 140         298 $factorial = $factorial -> binc();
3028 140         332 $below = $below->bmul($factorial); # n*(n+1)
3029             $factorial = $factorial -> binc();
3030             }
3031              
3032 16 50       71 # shortcut to not run through _find_round_parameters again
3033 16         48 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       60 }
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         60 # restore globals
3044 16         38 $$abr = $ab;
3045             $$pbr = $pb;
3046 16 50 33     75  
3047             return $downgrade -> new($x -> bdstr(), @r)
3048 16         372 if defined($downgrade) && $x -> is_int();
3049             $x;
3050             }
3051              
3052             sub batan {
3053 175 50   175 1 1641 # 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       649  
3060             return $x if $x->modify('batan');
3061 175 100       515  
3062             return $x -> bnan(@r) if $x->is_nan();
3063              
3064             # We need to limit the accuracy to protect against overflow.
3065 171         383  
3066 171         341 my $fallback = 0;
3067 171         502 my ($scale, @params);
3068             ($x, @params) = $x->_find_round_parameters(@r);
3069              
3070             # Error in _find_round_parameters?
3071 171 50       602  
3072             return $x -> bnan(@r) if $x->is_nan();
3073 171 100       675  
3074             if ($x->{sign} =~ /^[+-]inf\z/) {
3075             # +inf result is PI/2
3076             # -inf result is -PI/2
3077 16         71 # calculate PI/2
3078             my $pi = $class->bpi(@r);
3079 16         65 # modify $x in place
3080 16         52 $x->{_m} = $pi->{_m};
3081 16         44 $x->{_e} = $pi->{_e};
3082             $x->{_es} = $pi->{_es};
3083 16         64 # -y => -PI/2, +y => PI/2
3084 16         69 $x->{sign} = substr($x->{sign}, 0, 1); # "+inf" => "+"
3085 16         211 $x -> {_m} = $LIB->_div($x->{_m}, $LIB->_new(2));
3086             return $x;
3087             }
3088 155 100       506  
3089             return $x->bzero(@r) if $x->is_zero();
3090              
3091 129 50       504 # 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     439 # 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     499 # inlined is_one() && is_one('-')
3107 27         124 if ($LIB->_is_one($x->{_m}) && $LIB->_is_zero($x->{_e})) {
3108             my $pi = $class->bpi($scale - 3);
3109 27         122 # modify $x in place
3110 27         90 $x->{_m} = $pi->{_m};
3111 27         70 $x->{_e} = $pi->{_e};
3112             $x->{_es} = $pi->{_es};
3113 27         108 # leave the sign of $x alone (+1 => +PI/4, -1 => -PI/4)
3114 27         222 $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         243 # calculate PI/2 - atan(1/x):
3120 102 100       317 my $pi = undef;
3121             if ($x->bacmp($x->copy()->bone) >= 0) {
3122 40         325 # calculate PI/2
3123 40         242 $pi = $class->bpi($scale - 3);
3124             $pi->{_m} = $LIB->_div($pi->{_m}, $LIB->_new(2));
3125 40         196 # calculate 1/$x:
3126             my $x_copy = $x->copy();
3127 40         140 # modify $x in place
3128 40         206 $x = $x->bone();
3129             $x = $x->bdiv($x_copy, $scale);
3130             }
3131 102         446  
3132 102         461 my $fmul = 1;
3133 174         328 foreach (0 .. int($scale / 20)) {
3134 174         515 $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   447 # disable them and later re-enable them.
  43         135  
  43         52748  
3141 102         325 no strict 'refs';
3142 102         364 my $abr = "$class\::accuracy";
3143 102         293 my $ab = $$abr;
3144 102         273 $$abr = undef;
3145 102         216 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         241 # took them already into account), since these would interfere, too
3150 102         190 $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         211  
3157 102         224 local $Math::BigInt::upgrade = undef;
3158             local $Math::BigFloat::downgrade = undef;
3159 102         386  
3160 102         357 my $over = $x * $x; # X ^ 2
3161 102         397 my $x2 = $over->copy(); # X ^ 2; difference between terms
3162 102         298 $over = $over->bmul($x); # X ^ 3 as starting value
3163 102         359 my $sign = 1; # start with -=
3164 102         549 my $below = $class->new(3);
3165 102         467 my $two = $class->new(2);
3166 102         270 $x->{_a} = undef;
3167             $x->{_p} = undef;
3168 102         387  
3169             my $limit = $class->new("1E-". ($scale-1));
3170 102         325 #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         2594 # stop:
3175 990 100       2688 my $next = $over->copy()->bdiv($below, $scale);
3176             last if $next->bacmp($limit) <= 0;
3177 888 100       1875  
3178 416         1126 if ($sign == 0) {
3179             $x = $x->badd($next);
3180 472         1254 } else {
3181             $x = $x->bsub($next);
3182 888         1448 }
3183             $sign = 1-$sign; # alternatex
3184 888         2161 # calculate things for the next term
3185 888         2304 $over = $over->bmul($x2); # $x*$x
3186             $below = $below->badd($two); # n += 2
3187 102         491 }
3188             $x = $x->bmul($fmul);
3189 102 100       497  
3190 40         291 if (defined $pi) {
3191             my $x_copy = $x->copy();
3192 40         186 # modify $x in place
3193 40         106 $x->{_m} = $pi->{_m};
3194 40         89 $x->{_e} = $pi->{_e};
3195             $x->{_es} = $pi->{_es};
3196 40         109 # PI/2 - $x
3197             $x = $x->bsub($x_copy);
3198             }
3199              
3200 102 50       426 # Shortcut to not run through _find_round_parameters again.
3201 102         377 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       457 }
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         360 # restore globals
3213 102         256 $$abr = $ab;
3214             $$pbr = $pb;
3215 102 0 0     280  
      33        
3216             return $downgrade -> new($x -> bdstr(), @r)
3217 102         1939 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 3235 # 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       848 # Quick exit if $y is read-only.
3230             return $y if $y -> modify('batan2');
3231              
3232 213 100 100     991 # Handle all NaN cases.
3233             return $y -> bnan() if $x->{sign} eq $nan || $y->{sign} eq $nan;
3234              
3235 201         368 # We need to limit the accuracy to protect against overflow.
3236 201         350 my $fallback = 0;
3237 201         674 my ($scale, @params);
3238             ($y, @params) = $y -> _find_round_parameters(@r);
3239              
3240 201 50       707 # Error in _find_round_parameters?
3241             return $y if $y->is_nan();
3242              
3243 201 100       576 # No rounding at all, so must use fallback.
3244             if (scalar @params == 0) {
3245 45         156 # Simulate old behaviour
3246 45         114 $params[0] = $class -> div_scale(); # and round to it as accuracy
3247 45         75 $params[1] = undef; # disable P
3248 45         81 $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     411 # enough ...
3254             $scale = abs($params[0] || $params[1]) + 4; # take whatever is defined
3255             }
3256 201 100       506  
    100          
    100          
    100          
3257 20 100       76 if ($x -> is_inf("+")) { # x = inf
    100          
3258 4         41 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         108 } else { # -inf < y < inf
3263             return $y -> bzero(@r); # 0
3264             }
3265 20 100       97 } elsif ($x -> is_inf("-")) { # x = -inf
    100          
    100          
3266 4         23 if ($y -> is_inf("+")) { # y = inf
3267             $y = $y -> bpi($scale) -> bmul("0.75"); # 3/4 pi
3268 4         50 } elsif ($y -> is_inf("-")) { # y = -inf
3269             $y = $y -> bpi($scale) -> bmul("-0.75"); # -3/4 pi
3270 8         39 } elsif ($y >= 0) { # y >= 0
3271             $y = $y -> bpi($scale); # pi
3272 4         45 } else { # y < 0
3273             $y = $y -> bpi($scale) -> bneg(); # -pi
3274             }
3275 87 100       398 } 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         52 } elsif ($y -> is_inf("-")) { # y = -inf
3279             $y = $y -> bpi($scale) -> bmul("-0.5"); # -pi/2
3280 79         416 } else { # -inf < y < inf
3281             $y = $y -> bdiv($x, $scale) -> batan($scale); # atan(y/x)
3282             }
3283 20         77 } elsif ($x < 0) { # -inf < x < 0
3284 20 100       70 my $pi = $class -> bpi($scale);
3285 12         77 if ($y >= 0) { # y >= 0
3286             $y = $y -> bdiv($x, $scale) -> batan() # atan(y/x) + pi
3287             -> badd($pi);
3288 8         49 } else { # y < 0
3289             $y = $y -> bdiv($x, $scale) -> batan() # atan(y/x) - pi
3290             -> bsub($pi);
3291             }
3292 54 100       173 } else { # x = 0
    100          
3293 25         181 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         73 } else { # y = 0
3298             return $y -> bzero(@r); # 0
3299             }
3300             }
3301 182         733  
3302             $y = $y -> round(@r);
3303 182 100       553  
3304 42         105 if ($fallback) {
3305 42         84 $y->{_a} = undef;
3306             $y->{_p} = undef;
3307             }
3308 182         2294  
3309             return $y;
3310             }
3311              
3312             sub bsqrt {
3313 442 100   442 1 2794 # calculate square root
3314             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
3315 442 50       1496  
3316             return $x if $x->modify('bsqrt');
3317              
3318             # Handle trivial cases.
3319 442 100       1234  
3320 434 100       1266 return $x -> bnan(@r) if $x->is_nan();
3321 430 100 100     1074 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       1486  
3326 20 50       78 if ($x -> is_neg()) {
3327 20         58 return $upgrade -> bsqrt($x, @r) if defined($upgrade);
3328             return $x -> bnan(@r);
3329             }
3330              
3331 400         773 # we need to limit the accuracy to protect against overflow
3332 400         677 my $fallback = 0;
3333 400         1138 my (@params, $scale);
3334             ($x, @params) = $x->_find_round_parameters(@r);
3335              
3336 400 50       1182 # error in _find_round_parameters?
3337             return $x -> bnan(@r) if $x->is_nan();
3338              
3339 400 100       1204 # no rounding at all, so must use fallback
3340             if (scalar @params == 0) {
3341 123         393 # simulate old behaviour
3342 123         286 $params[0] = $class->div_scale(); # and round to it as accuracy
3343 123         223 $scale = $params[0]+4; # at least four more for proper round
3344 123         227 $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     852 # 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   441 # disable them and later re-enable them
  43         165  
  43         42356  
3354 400         959 no strict 'refs';
3355 400         913 my $abr = "$class\::accuracy";
3356 400         915 my $ab = $$abr;
3357 400         794 $$abr = undef;
3358 400         813 my $pbr = "$class\::precision";
3359 400         801 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         794 # them already into account), since these would interfere, too
3363 400         674 $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         691  
3370 400         622 local $Math::BigInt::upgrade = undef;
3371             local $Math::BigFloat::downgrade = undef;
3372 400         1178  
3373 400 100       1204 my $i = $LIB->_copy($x->{_m});
3374 400         1991 $i = $LIB->_lsft($i, $x->{_e}, 10) unless $LIB->_is_zero($x->{_e});
3375 400         984 my $xas = Math::BigInt->bzero();
3376             $xas->{value} = $i;
3377 400         1496  
3378             my $gs = $xas->copy()->bsqrt(); # some guess
3379 400 100 100     2007  
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         139 # exact result, copy result over to keep $x
3385 50         176 $x->{_m} = $gs->{value};
3386 50         110 $x->{_e} = $LIB->_zero();
3387 50         137 $x->{_es} = '+';
3388             $x = $x->bnorm();
3389 50 50       128 # shortcut to not run through _find_round_parameters again
3390 50         162 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       115 }
3395             if ($fallback) {
3396 33         71 # clear a/p after round, since user did not request it
3397 33         67 $x->{_a} = undef;
3398             $x->{_p} = undef;
3399             }
3400 50         78 # re-enable A and P, upgrade is taken care of by "local"
  50         155  
3401 50         81 ${"$class\::accuracy"} = $ab;
  50         110  
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         1281 # The following steps will transform 123.456 (in $x) into 123456 (in $y1)
3412             my $y1 = $LIB->_copy($x->{_m});
3413 350         1039  
3414             my $length = $LIB->_len($y1);
3415              
3416 350         979 # Now calculate how many digits the result of sqrt(y1) would have
3417             my $digits = int($length / 2);
3418              
3419 350         581 # 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       783 # (we take care of integer guesses above)
3424             $shift = 0 if $shift < 0;
3425              
3426 350         560 # 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       1143 # after the dot (the result is still odd or even digits long).
3435             $s2++ if $LIB->_is_odd($x->{_e});
3436 350         1076  
3437             $y1 = $LIB->_lsft($y1, $LIB->_new($s2), 10);
3438              
3439 350         1381 # 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         1296 # on sign of $dat, the result should have half as many:
3447 350 100       1236 my $dat = $LIB->_num($x->{_e});
3448 350         613 $dat = -$dat if $x->{_es} eq '-';
3449             $dat += $length;
3450 350 100       799  
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         826 # (but round this "up")
3455             $dat = int(($dat+1)/2);
3456 30         96 } else {
3457             $dat = int(($dat)/2);
3458 350         1005 }
3459 350 50       867 $dat -= $LIB->_len($y1);
3460 350         837 if ($dat < 0) {
3461 350         997 $dat = abs($dat);
3462 350         848 $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         834 }
3468 350         1147 $x->{_m} = $y1;
3469             $x = $x->bnorm();
3470              
3471 350 100       1094 # shortcut to not run through _find_round_parameters again
3472 328         1079 if (defined $params[0]) {
3473             $x = $x->bround($params[0], $params[2]); # then round accordingly
3474 22         86 } else {
3475             $x = $x->bfround($params[1], $params[2]); # then round accordingly
3476 350 100       1169 }
3477             if ($fallback) {
3478 90         214 # clear a/p after round, since user did not request it
3479 90         189 $x->{_a} = undef;
3480             $x->{_p} = undef;
3481             }
3482 350         1088 # restore globals
3483 350         810 $$abr = $ab;
3484             $$pbr = $pb;
3485 350 0 0     865  
      33        
3486             return $downgrade -> new($x -> bdstr(), @r)
3487 350         3048 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 3049 # 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       679  
3499             return $x if $x->modify('broot');
3500              
3501             # Handle trivial cases.
3502 186 100 100     503  
3503             return $x -> bnan(@r) if $x->is_nan() || $y->is_nan();
3504 150 100       465  
3505             if ($x -> is_neg()) {
3506 44 100 66     171 # -27 ** (1/3) = -3
      100        
3507             return $x -> broot($y -> copy() -> bneg(), @r) -> bneg()
3508 40 50       97 if $x -> is_int() && $y -> is_int() && $y -> is_neg();
3509 40         110 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     579 return $x->bnan() if $x->{sign} !~ /^\+/ || $y->is_zero() ||
      100        
3515             $y->{sign} !~ /^\+$/;
3516 82 100 100     199  
      100        
      100        
3517             return $x if $x->is_zero() || $x->is_one() || $x->is_inf() || $y->is_one();
3518              
3519 66         173 # we need to limit the accuracy to protect against overflow
3520 66         156 my $fallback = 0;
3521 66         237 my (@params, $scale);
3522             ($x, @params) = $x->_find_round_parameters(@r);
3523 66 50       190  
3524             return $x if $x->is_nan(); # error in _find_round_parameters?
3525              
3526 66 50       191 # no rounding at all, so must use fallback
3527             if (scalar @params == 0) {
3528 66         206 # simulate old behaviour
3529 66         118 $params[0] = $class->div_scale(); # and round to it as accuracy
3530 66         127 $scale = $params[0]+4; # at least four more for proper round
3531 66         120 $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   430 # disable them and later re-enable them
  43         143  
  43         629426  
3541 66         139 no strict 'refs';
3542 66         140 my $abr = "$class\::accuracy";
3543 66         168 my $ab = $$abr;
3544 66         123 $$abr = undef;
3545 66         134 my $pbr = "$class\::precision";
3546 66         119 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         126 # them already into account), since these would interfere, too
3550 66         115 $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         162  
3557 66         108 local $Math::BigInt::upgrade = undef;
3558             local $Math::BigFloat::downgrade = undef;
3559              
3560 66         101 # remember sign and make $x positive, since -4 ** (1/2) => -2
3561 66 50       180 my $sign = 0;
3562 66         128 $sign = 1 if $x->{sign} eq '-';
3563             $x->{sign} = '+';
3564 66         101  
3565 66 50       164 my $is_two = 0;
3566             if ($y->isa('Math::BigFloat')) {
3567 66   66     343 $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       212 # normal square root if $y == 2:
    50          
3574 46         134 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         59  
3587 20 50 33     67 my $done = 0; # not yet
3588 20         93 if ($y->is_int() && $x->is_int()) {
3589 20 50       80 my $i = $LIB->_copy($x->{_m});
3590 20         95 $i = $LIB->_lsft($i, $x->{_e}, 10) unless $LIB->_is_zero($x->{_e});
3591 20         63 my $int = Math::BigInt->bzero();
3592 20         95 $int->{value} = $i;
3593             $int = $int->broot($y->as_number());
3594 20 100       164 # if ($exact)
3595             if ($int->copy()->bpow($y) == $x) {
3596 14         57 # found result, return it
3597 14         48 $x->{_m} = $int->{value};
3598 14         45 $x->{_e} = $LIB->_zero();
3599 14         77 $x->{_es} = '+';
3600 14         84 $x = $x->bnorm();
3601             $done = 1;
3602             }
3603 20 100       104 }
3604 6         38 if ($done == 0) {
3605 6         37 my $u = $class->bone()->bdiv($y, $scale+4);
3606 6         18 $u->{_a} = undef;
3607 6         26 $u->{_p} = undef;
3608             $x = $x->bpow($u, $scale+4); # el cheapo
3609             }
3610 66 50       198 }
3611             $x = $x->bneg() if $sign == 1;
3612              
3613 66 50       146 # shortcut to not run through _find_round_parameters again
3614 66         179 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       229 }
3619             if ($fallback) {
3620 66         137 # clear a/p after round, since user did not request it
3621 66         124 $x->{_a} = undef;
3622             $x->{_p} = undef;
3623             }
3624 66         173 # restore globals
3625 66         144 $$abr = $ab;
3626             $$pbr = $pb;
3627 66 0 0     148  
      33        
3628             return $downgrade -> new($x -> bdstr(), @r)
3629 66         971 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 990 # set up parameters
3638             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
3639              
3640 80 50       290 # inf => inf
3641             return $x if $x->modify('bfac');
3642 80 100 100     230  
3643 72 100       262 return $x -> bnan(@r) if $x->is_nan() || $x->is_inf("-");
3644 68 100 100     216 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       21 if ($x -> is_neg() || !$x -> is_int()) {
3648 4         21 return $upgrade -> bfac($x, @r) if defined($upgrade);
3649             return $x -> bnan(@r);
3650             }
3651 56 100       167  
3652 8         41 if (! $LIB->_is_zero($x->{_e})) {
3653 8         52 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # change 12e1 to 120e0
3654 8         23 $x->{_e} = $LIB->_zero(); # normalize
3655             $x->{_es} = '+';
3656 56         185 }
3657             $x->{_m} = $LIB->_fac($x->{_m}); # calculate factorial
3658 56         160  
3659             $x = $x->bnorm()->round(@r); # norm again and round result
3660 56 0 0     131  
      33        
3661             return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade)
3662 56         590 && ($x -> is_int() || $x -> is_inf());
3663             $x;
3664             }
3665              
3666             sub bdfac {
3667             # compute double factorial
3668              
3669 72 50   72 1 860 # set up parameters
3670             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
3671 72 50       243  
3672             return $x if $x->modify('bdfac');
3673 72 100 100     210  
3674 64 100       204 return $x -> bnan(@r) if $x->is_nan() || $x->is_inf("-");
3675             return $x -> binf("+", @r) if $x->is_inf("+");
3676 60 100 66     242  
3677 4 50       24 if ($x <= -2 || !$x -> is_int()) {
3678 4         18 return $upgrade -> bdfac($x, @r) if defined($upgrade);
3679             return $x -> bnan(@r);
3680             }
3681 56 100       161  
3682             return $x->bone() if $x <= 1;
3683 44 50       351  
3684             croak("bdfac() requires a newer version of the $LIB library.")
3685             unless $LIB->can('_dfac');
3686 44 100       143  
3687 4         39 if (! $LIB->_is_zero($x->{_e})) {
3688 4         34 $x->{_m} = $LIB->_lsft($x->{_m}, $x->{_e}, 10); # change 12e1 to 120e0
3689 4         24 $x->{_e} = $LIB->_zero(); # normalize
3690             $x->{_es} = '+';
3691 44         170 }
3692             $x->{_m} = $LIB->_dfac($x->{_m}); # calculate factorial
3693 44         143  
3694             $x = $x->bnorm()->round(@r); # norm again and round result
3695 44 50 33     146  
3696             return $downgrade -> new($x -> bdstr(), @r)
3697 44         464 if defined($downgrade) && $x -> is_int();
3698             return $x;
3699             }
3700              
3701             sub btfac {
3702             # compute triple factorial
3703              
3704 76 50   76 1 969 # set up parameters
3705             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
3706 76 50       296  
3707             return $x if $x->modify('btfac');
3708 76 100 100     233  
3709 68 100       243 return $x -> bnan(@r) if $x->is_nan() || $x->is_inf("-");
3710             return $x -> binf("+", @r) if $x->is_inf("+");
3711 64 100 66     230  
3712 4 50       29 if ($x <= -3 || !$x -> is_int()) {
3713 4         30 return $upgrade -> btfac($x, @r) if defined($upgrade);
3714             return $x -> bnan(@r);
3715             }
3716 60         198  
3717 60 50       236 my $k = $class -> new("3");
3718             return $x->bnan(@r) if $x <= -$k;
3719 60         286  
3720 60 100       155 my $one = $class -> bone();
3721             return $x->bone(@r) if $x <= $one;
3722 44         150  
3723 44         144 my $f = $x -> copy();
3724 60         190 while ($f -> bsub($k) > $one) {
3725             $x = $x -> bmul($f);
3726             }
3727 44         127  
3728             $x = $x->round(@r);
3729 44 50 33     141  
3730             return $downgrade -> new($x -> bdstr(), @r)
3731 44         804 if defined($downgrade) && $x -> is_int();
3732             return $x;
3733             }
3734              
3735 364 50 33 364 1 5554 sub bmfac {
3736             my ($class, $x, $k, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3737             ? (ref($_[0]), @_)
3738             : objectify(2, @_);
3739 364 50       1298  
3740             return $x if $x->modify('bmfac');
3741 364 100 100     1048  
      100        
3742 308 100       808 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     972  
      100        
      100        
3745             if ($x <= -$k || !$x -> is_int() ||
3746             ($k -> is_finite() && !$k -> is_int()))
3747 24 50       93 {
3748 24         82 return $upgrade -> bmfac($x, $k, @r) if defined($upgrade);
3749             return $x -> bnan(@r);
3750             }
3751 264         1357  
3752 264 100       658 my $one = $class -> bone();
3753             return $x->bone(@r) if $x <= $one;
3754 184         531  
3755 184         514 my $f = $x -> copy();
3756 284         787 while ($f -> bsub($k) > $one) {
3757             $x = $x -> bmul($f);
3758             }
3759 184         583  
3760             $x = $x->round(@r);
3761 184 50 33     467  
3762             return $downgrade -> new($x -> bdstr(), @r)
3763 184         3237 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 540 # 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       135  
3776             return $x if $x -> modify('blsft');
3777 33 100 66     108  
3778             return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan();
3779 29 50       106  
3780 29 50 33     188 $b = 2 if !defined $b;
3781             $b = $class -> new($b)
3782 29 50       154 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       106 # shift by a negative amount?
3788             return $x -> brsft($y -> copy() -> babs(), $b) if $y -> {sign} =~ /^-/;
3789 29         90  
3790             $x = $x -> bmul($b -> bpow($y), $r[0], $r[1], $r[2], $y);
3791 29 0 0     96  
      33        
3792             return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade)
3793 29         299 && ($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 728 # 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       206  
3806             return $x if $x -> modify('brsft');
3807 48 100 66     132  
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       135  
3812 44 50 66     254 $b = 2 if !defined $b;
3813             $b = $class -> new($b)
3814 44 50       203 unless defined(blessed($b)) && $b -> isa(__PACKAGE__);
3815             return $x -> bnan(@r) if $b -> is_nan();
3816              
3817 44 50       227 # shift by a negative amount?
3818             return $x -> blsft($y -> copy() -> babs(), $b) if $y -> {sign} =~ /^-/;
3819              
3820 44         162 # call bdiv()
3821             $x = $x -> bdiv($b -> bpow($y), $r[0], $r[1], $r[2], $y);
3822 44 0 0     141  
      33        
3823             return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade)
3824 44         649 && ($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 74 sub bblsft {
3835             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3836             ? (ref($_[0]), @_)
3837             : objectify(2, @_);
3838 8         44  
3839             my $xint = Math::BigInt -> bblsft($x, $y, @r);
3840              
3841             # disable downgrading
3842 8         60  
3843 8         50 my $dng = $class -> downgrade();
3844             $class -> downgrade(undef);
3845              
3846             # convert to Math::BigFloat without downgrading
3847 8         41  
3848             my $xflt = $class -> new($xint);
3849              
3850             # reset downgrading
3851 8         53  
3852             $class -> downgrade($dng);
3853 8         24  
3854 8         28 $x -> {sign} = $xflt -> {sign};
3855 8         24 $x -> {_m} = $xflt -> {_m};
3856 8         25 $x -> {_es} = $xflt -> {_es};
3857             $x -> {_e} = $xflt -> {_e};
3858              
3859             # now we might downgrade
3860 8 50       37  
3861 8         37 return $downgrade -> new($x) if defined($downgrade);
3862             $x -> round(@r);
3863             }
3864              
3865             # Bitwise left shift.
3866              
3867 8 50 33 8 1 75 sub bbrsft {
3868             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3869             ? (ref($_[0]), @_)
3870             : objectify(2, @_);
3871 8         52  
3872             my $xint = Math::BigInt -> bbrsft($x, $y, @r);
3873              
3874             # disable downgrading
3875 8         36  
3876 8         38 my $dng = $class -> downgrade();
3877             $class -> downgrade(undef);
3878              
3879             # convert to Math::BigFloat without downgrading
3880 8         36  
3881             my $xflt = $class -> new($xint);
3882              
3883             # reset downgrading
3884 8         57  
3885             $class -> downgrade($dng);
3886 8         25  
3887 8         27 $x -> {sign} = $xflt -> {sign};
3888 8         20 $x -> {_m} = $xflt -> {_m};
3889 8         17 $x -> {_es} = $xflt -> {_es};
3890             $x -> {_e} = $xflt -> {_e};
3891              
3892             # now we might downgrade
3893 8 50       27  
3894 8         38 return $downgrade -> new($x) if defined($downgrade);
3895             $x -> round(@r);
3896             }
3897              
3898 1 50 33 1 1 21 sub band {
3899             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3900             ? (ref($_[0]), @_)
3901             : objectify(2, @_);
3902 1 50       10  
3903             return if $x -> modify('band');
3904 1 50 33     5  
3905             return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan();
3906 1         13  
3907 1         3 my $xint = $x -> as_int(); # to Math::BigInt
3908             my $yint = $y -> as_int(); # to Math::BigInt
3909 1         14  
3910             $xint = $xint -> band($yint);
3911 1 50       7  
3912             return $xint -> round(@r) if defined $downgrade;
3913 1         5  
3914 1         4 my $xflt = $class -> new($xint); # back to Math::BigFloat
3915 1         5 $x -> {sign} = $xflt -> {sign};
3916 1         3 $x -> {_m} = $xflt -> {_m};
3917 1         4 $x -> {_es} = $xflt -> {_es};
3918             $x -> {_e} = $xflt -> {_e};
3919 1         4  
3920             return $x -> round(@r);
3921             }
3922              
3923 1 50 33 1 1 9 sub bior {
3924             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3925             ? (ref($_[0]), @_)
3926             : objectify(2, @_);
3927 1 50       6  
3928             return if $x -> modify('bior');
3929 1 50 33     4  
3930             return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan();
3931 1         3  
3932 1         4 my $xint = $x -> as_int(); # to Math::BigInt
3933             my $yint = $y -> as_int(); # to Math::BigInt
3934 1         8  
3935             $xint = $xint -> bior($yint);
3936 1 50       8  
3937             return $xint -> round(@r) if defined $downgrade;
3938 1         19  
3939 1         4 my $xflt = $class -> new($xint); # back to Math::BigFloat
3940 1         3 $x -> {sign} = $xflt -> {sign};
3941 1         3 $x -> {_m} = $xflt -> {_m};
3942 1         3 $x -> {_es} = $xflt -> {_es};
3943             $x -> {_e} = $xflt -> {_e};
3944 1         4  
3945             return $x -> round(@r);
3946             }
3947              
3948 1 50 33 1 1 9 sub bxor {
3949             my ($class, $x, $y, @r) = ref($_[0]) && ref($_[0]) eq ref($_[1])
3950             ? (ref($_[0]), @_)
3951             : objectify(2, @_);
3952 1 50       14  
3953             return if $x -> modify('bxor');
3954 1 50 33     9  
3955             return $x -> bnan(@r) if $x -> is_nan() || $y -> is_nan();
3956 1         6  
3957 1         4 my $xint = $x -> as_int(); # to Math::BigInt
3958             my $yint = $y -> as_int(); # to Math::BigInt
3959 1         7  
3960             $xint = $xint -> bxor($yint);
3961 1 50       5  
3962             return $xint -> round(@r) if defined $downgrade;
3963 1         5  
3964 1         5 my $xflt = $class -> new($xint); # back to Math::BigFloat
3965 1         4 $x -> {sign} = $xflt -> {sign};
3966 1         2 $x -> {_m} = $xflt -> {_m};
3967 1         3 $x -> {_es} = $xflt -> {_es};
3968             $x -> {_e} = $xflt -> {_e};
3969 1         3  
3970             return $x -> round(@r);
3971             }
3972              
3973 4 50   4 1 22 sub bnot {
3974             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
3975 4 50       16  
3976             return if $x -> modify('bnot');
3977 4 50       13  
3978             return $x -> bnan(@r) if $x -> is_nan();
3979 4         18  
3980 4         12 my $xint = $x -> as_int(); # to Math::BigInt
3981             $xint = $xint -> bnot();
3982 4 50       19  
3983             return $xint -> round(@r) if defined $downgrade;
3984 4         14  
3985 4         13 my $xflt = $class -> new($xint); # back to Math::BigFloat
3986 4         14 $x -> {sign} = $xflt -> {sign};
3987 4         8 $x -> {_m} = $xflt -> {_m};
3988 4         9 $x -> {_es} = $xflt -> {_es};
3989             $x -> {_e} = $xflt -> {_e};
3990 4         16  
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 133378  
4001             my ($class, $x, @a) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4002 41467 100 100     109437  
4003 3         743 if (($a[0] || 0) < 0) {
4004             croak('bround() needs positive accuracy');
4005             }
4006 41464 50       112632  
4007             return $x if $x->modify('bround');
4008 41464         104408  
4009 41464 100       91566 my ($scale, $mode) = $x->_scale_a(@a);
4010 3889 50 66     7575 if (!defined $scale) { # no-op
      66        
4011             return $downgrade -> new($x) if defined($downgrade)
4012 3885         11165 && ($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     122402  
4021 4 0 0     36 if (defined $x->{_a} && $x->{_a} < $scale) {
      33        
4022             return $downgrade -> new($x) if defined($downgrade)
4023 4         15 && ($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     163085  
4031 12 0 0     59 if ($scale <= 0 || $x->{sign} !~ /^[+-]$/) {
      33        
4032             return $downgrade -> new($x) if defined($downgrade)
4033 12         141 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4034             return $x;
4035             }
4036              
4037             # 1: never round a 0
4038 37559 100 100     87194 # 2: if we should keep more digits than the mantissa has, do nothing
4039 12805 100 100     46470 if ($x->is_zero() || $LIB->_len($x->{_m}) <= $scale) {
4040 12805 50 33     24619 $x->{_a} = $scale if !defined $x->{_a} || $x->{_a} > $scale;
      66        
4041             return $downgrade -> new($x) if defined($downgrade)
4042 12805         39117 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4043             return $x;
4044             }
4045              
4046 24754         98932 # pass sign to bround for '+inf' and '-inf' rounding modes
4047             my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
4048 24754         70439  
4049 24754         55345 $m = $m->bround($scale, $mode); # round mantissa
4050 24754         59653 $x->{_m} = $m->{value}; # get our mantissa back
4051 24754         39998 $x->{_a} = $scale; # remember rounding
4052             $x->{_p} = undef; # and clear P
4053              
4054 24754         61860 # 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 10366  
4063             my ($class, $x, @p) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4064 724 50       2338  
4065             return $x if $x->modify('bfround'); # no-op
4066 724         2060  
4067 724 100       1722 my ($scale, $mode) = $x->_scale_p(@p);
4068 4 50 66     18 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       1737  
4076 20 50 33     187 if ($x->is_zero()) {
4077 20 0 0     79 $x->{_p} = $scale if !defined $x->{_p} || $x->{_p} < $scale; # -3 < -2
      33        
4078             return $downgrade -> new($x) if defined($downgrade)
4079 20         94 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4080             return $x;
4081             }
4082 700 100       2628  
4083 12 0 0     41 if ($x->{sign} !~ /^[+-]$/) {
      33        
4084             return $downgrade -> new($x) if defined($downgrade)
4085 12         127 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4086             return $x;
4087             }
4088              
4089 688 50 100     1883 # 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         1311  
4096 688         1197 $x->{_p} = $scale; # remember round in any case
4097 688 100       1520 $x->{_a} = undef; # and clear A
4098             if ($scale < 0) {
4099             # round right from the '.'
4100 532 100       1184  
4101 28 0 0     100 if ($x->{_es} eq '+') { # e >= 0 => nothing to round
      33        
4102             return $downgrade -> new($x) if defined($downgrade)
4103 28         104 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4104             return $x;
4105             }
4106 504         834  
4107 504         1528 $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         1562 # scalar, you got other problems (memory etc) anyway
4112 504         920 my $dad = -(0+ ($x->{_es}.$LIB->_num($x->{_e}))); # digits after dot
4113 504 100       1165 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       979  
4126 47 0 0     119 if ($scale > $dad) { # 0.123, scale >= 3 => exit
      33        
4127             return $downgrade -> new($x) if defined($downgrade)
4128 47         468 && ($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       1009 # case)
4135 40 0 0     132 if ($scale < $zad) {
      33        
4136             return $downgrade -> new($x) if defined($downgrade)
4137 40         147 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4138             return $x->bzero();
4139             }
4140 417 100       875  
4141 41         100 if ($scale == $zad) { # for 0.006, scale -3 and trunc
4142             $scale = -$len;
4143             } else {
4144 376 100       756 # adjust round-point to be inside mantissa
4145 78         154 if ($zad != 0) {
4146             $scale = $scale-$zad;
4147 298         498 } else {
4148 298 50       1060 my $dbd = $len - $dad;
4149 298         502 $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         555  
4158             my $dbt = $LIB->_len($x->{_m});
4159 156         537 # digits before dot
4160             my $dbd = $dbt + ($x->{_es} . $LIB->_num($x->{_e}));
4161 156 100       482 # should be the same, so treat it as this
4162             $scale = 1 if $scale == 0;
4163 156 100 100     566 # shortcut if already integer
4164 10 0 0     36 if ($scale == 1 && $dbt <= $dbd) {
      33        
4165             return $downgrade -> new($x) if defined($downgrade)
4166 10         34 && ($x->is_int() || $x->is_inf() || $x->is_nan());
4167             return $x;
4168             }
4169 146         219 # maximum digits before dot
4170             ++$dbd;
4171 146 100       413  
    100          
4172             if ($scale > $dbd) {
4173 27 50       75 # not enough digits before dot, so round to zero
4174 27         93 return $downgrade -> new($x) if defined($downgrade);
4175             return $x->bzero;
4176             } elsif ($scale == $dbd) {
4177 72         134 # maximum
4178             $scale = -$dbt;
4179 47         92 } else {
4180             $scale = $dbd - $scale;
4181             }
4182             }
4183              
4184 536         2111 # pass sign to bround for rounding modes '+inf' and '-inf'
4185 536         1724 my $m = bless { sign => $x->{sign}, value => $x->{_m} }, 'Math::BigInt';
4186 536         1403 $m = $m->bround($scale, $mode);
4187             $x->{_m} = $m->{value}; # get our mantissa back
4188              
4189 536         1483 # 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 837 # round towards minus infinity
4195             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4196 84 50       333  
4197             return $x if $x->modify('bfloor');
4198 84 100       267  
4199             return $x -> bnan(@r) if $x -> is_nan();
4200 79 100       361  
4201             if ($x->{sign} =~ /^[+-]$/) {
4202 70 100       254 # if $x has digits after dot, remove them
4203 47         243 if ($x->{_es} eq '-') {
4204 47         158 $x->{_m} = $LIB->_rsft($x->{_m}, $x->{_e}, 10);
4205 47         158 $x->{_e} = $LIB->_zero();
4206             $x->{_es} = '+';
4207 47 100       185 # increment if negative
4208             $x->{_m} = $LIB->_inc($x->{_m}) if $x->{sign} eq '-';
4209 70         237 }
4210             $x = $x->round(@r);
4211 79 100       251 }
4212 77         533 return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade);
4213             return $x;
4214             }
4215              
4216             sub bceil {
4217 43 50   43 1 568 # round towards plus infinity
4218             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4219 43 50       177  
4220             return $x if $x->modify('bceil');
4221 43 100       138  
4222             return $x -> bnan(@r) if $x -> is_nan();
4223              
4224 38 100       168 # if $x has digits after dot, remove them
4225 29 100       102 if ($x->{sign} =~ /^[+-]$/) {
4226 17         110 if ($x->{_es} eq '-') {
4227 17         77 $x->{_m} = $LIB->_rsft($x->{_m}, $x->{_e}, 10);
4228 17         50 $x->{_e} = $LIB->_zero();
4229 17 100       61 $x->{_es} = '+';
4230 9         57 if ($x->{sign} eq '+') {
4231             $x->{_m} = $LIB->_inc($x->{_m}); # increment if positive
4232 8 100       28 } else {
4233             $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # avoid -0
4234             }
4235 29         100 }
4236             $x = $x->round(@r);
4237             }
4238 38 100       168  
4239 36         315 return $downgrade -> new($x -> bdstr(), @r) if defined($downgrade);
4240             return $x;
4241             }
4242              
4243             sub bint {
4244 179 50   179 1 978 # round towards zero
4245             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4246 179 50       573  
4247             return $x if $x->modify('bint');
4248 179 100       483  
4249             return $x -> bnan(@r) if $x -> is_nan();
4250 174 100       749  
4251             if ($x->{sign} =~ /^[+-]$/) {
4252 165 100       425 # if $x has digits after the decimal point
4253 13         86 if ($x->{_es} eq '-') {
4254 13         56 $x->{_m} = $LIB->_rsft($x->{_m}, $x->{_e}, 10); # remove frac part
4255 13         50 $x->{_e} = $LIB->_zero(); # truncate/normalize
4256 13 100       51 $x->{_es} = '+'; # abs e
4257             $x->{sign} = '+' if $LIB->_is_zero($x->{_m}); # avoid -0
4258 165         456 }
4259             $x = $x->round(@r);
4260             }
4261 174 100       485  
4262 172         725 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 2331 # 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         487  
4283             my ($class, @args) = objectify(0, @_);
4284 133         264  
4285 133 50 33     496 my $x = shift @args;
4286             $x = defined(blessed($x)) && $x -> isa(__PACKAGE__) ? $x -> copy()
4287 133 100       386 : $class -> new($x);
4288             return $class->bnan() unless $x -> is_int();
4289 105         280  
4290 116         216 while (@args) {
4291 116 50 33     438 my $y = shift @args;
4292             $y = $class->new($y)
4293 116 100       282 unless defined(blessed($y)) && $y -> isa(__PACKAGE__);
4294             return $class->bnan() unless $y -> is_int();
4295              
4296 104         257 # greatest common divisor
4297 246         529 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         277 }
4303             $x = $x -> babs();
4304 93 50 33     237  
4305             return $downgrade -> new($x)
4306 93         982 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 1259 # 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         5 # " use is as a method instead";
4321             unshift @_, __PACKAGE__;
4322             }
4323 35         154  
4324             my ($class, @args) = objectify(0, @_);
4325 35         148  
4326 35 50 33     170 my $x = shift @args;
4327             $x = defined(blessed($x)) && $x -> isa(__PACKAGE__) ? $x -> copy()
4328 35 100       189 : $class -> new($x);
4329             return $class->bnan() if $x->{sign} !~ /^[+-]$/; # x NaN?
4330 27         83  
4331 30         63 while (@args) {
4332 30 50 33     135 my $y = shift @args;
4333             $y = $class -> new($y)
4334 30 100       97 unless defined(blessed($y)) && $y -> isa(__PACKAGE__);
4335 26         101 return $x->bnan() unless $y -> is_int();
4336 26         216 my $gcd = $x -> bgcd($y);
4337             $x = $x -> bdiv($gcd) -> bmul($y);
4338             }
4339 23         89  
4340             $x = $x -> babs();
4341 23 50 33     75  
4342             return $downgrade -> new($x)
4343 23         305 if defined $downgrade && $x->is_int();
4344             return $x;
4345             }
4346              
4347             ###############################################################################
4348             # Object property methods
4349             ###############################################################################
4350              
4351 24 50   24 1 345 sub length {
4352             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4353 24 50       71  
4354             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4355 24 100       97  
4356             return 1 if $LIB->_is_zero($x->{_m});
4357 20         97  
4358 20 50       117 my $len = $LIB->_len($x->{_m});
4359 20 50       71 $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         157 }
4365             $len;
4366             }
4367              
4368             sub mantissa {
4369 36 50   36 1 521 # 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       119  
4376             return $x -> bnan(@r) if $x -> is_nan();
4377 32 100       150  
4378 8         34 if ($x->{sign} !~ /^[+-]$/) {
4379 8         38 my $s = $x->{sign};
4380 8         38 $s =~ s/^\+//;
4381             return Math::BigInt->new($s, undef, undef); # -inf, +inf => +inf
4382 24         107 }
4383 24 100       121 my $m = Math::BigInt->new($LIB->_str($x->{_m}), undef, undef);
4384 24         131 $m = $m->bneg() if $x->{sign} eq '-';
4385             $m;
4386             }
4387              
4388             sub exponent {
4389 36 50   36 1 473 # 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       119  
4396             return $x -> bnan(@r) if $x -> is_nan();
4397 32 100       150  
4398 8         24 if ($x->{sign} !~ /^[+-]$/) {
4399 8         31 my $s = $x->{sign};
4400 8         31 $s =~ s/^[+-]//;
4401             return Math::BigInt->new($s, undef, undef); # -inf, +inf => +inf
4402 24         97 }
4403             Math::BigInt->new($x->{_es} . $LIB->_str($x->{_e}), undef, undef);
4404             }
4405              
4406             sub parts {
4407 32 50   32 1 449 # return a copy of both the exponent and the mantissa
4408             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4409 32 50       84  
4410             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4411 32 100       121  
4412 12         30 if ($x->{sign} !~ /^[+-]$/) {
4413 12         30 my $s = $x->{sign};
4414 12         26 $s =~ s/^\+//;
4415 12         28 my $se = $s;
4416             $se =~ s/^-//;
4417 12         38 # +inf => inf and -inf, +inf => inf
4418             return ($class->new($s), $class->new($se));
4419 20         77 }
4420 20         68 my $m = Math::BigInt->bzero();
4421 20 100       92 $m->{value} = $LIB->_copy($x->{_m});
4422 20         77 $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 1684324 # internal format is always normalized (no leading zeros, "-0" => "+0")
4685             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4686 8480 50       22737  
4687             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4688              
4689             # Inf and NaN
4690 8480 100 100     33176  
4691 2477 100       22499 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
4692 549         5723 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
4693             return 'inf'; # +inf
4694             }
4695              
4696             # Finite number
4697 6003         10644  
4698 6003         9384 my $es = '0';
4699 6003         8767 my $len = 1;
4700 6003         8985 my $cad = 0;
4701             my $dot = '.';
4702              
4703 6003   100     27368 # $x is zero?
4704 6003 100       13946 my $not_zero = !($x->{sign} eq '+' && $LIB->_is_zero($x->{_m}));
4705 4725         14142 if ($not_zero) {
4706 4725         8597 $es = $LIB->_str($x->{_m});
4707 4725         12214 $len = CORE::length($es);
4708 4725 100       12028 my $e = $LIB->_num($x->{_e});
4709 4725 100       12010 $e = -$e if $x->{_es} eq '-';
    100          
4710 1586         2662 if ($e < 0) {
4711             $dot = '';
4712 1586 100       3749 # if _e is bigger than a scalar, the following will blow your memory
4713 576         1174 if ($e <= -$len) {
4714 576         1690 my $r = abs($e) - $len;
4715 576         1167 $es = '0.'. ('0' x $r) . $es;
4716             $cad = -($len+$r);
4717 1010         2526 } else {
4718 1010         2618 substr($es, $e, 0) = '.';
4719 1010 50       3004 $cad = $LIB->_num($x->{_e});
4720             $cad = -$cad if $x->{_es} eq '-';
4721             }
4722             } elsif ($e > 0) {
4723 678         1968 # expand with zeros
4724 678         1201 $es .= '0' x $e;
4725 678         1174 $len += $e;
4726             $cad = 0;
4727             }
4728             } # if not zero
4729 6003 100       14332  
4730             $es = '-'.$es if $x->{sign} eq '-';
4731 6003 100 100     28141 # if set accuracy or precision, pad with zeros on the right side
    100 100        
4732             if ((defined $x->{_a}) && ($not_zero)) {
4733 694         1353 # 123400 => 6, 0.1234 => 4, 0.001234 => 4
4734 694 50       1702 my $zeros = $x->{_a} - $cad; # cad == 0 => 12340
4735 694 100       1540 $zeros = $x->{_a} - $len if $cad != $len;
4736             $es .= $dot.'0' x $zeros if $zeros > 0;
4737             } elsif ((($x->{_p} || 0) < 0)) {
4738 502         881 # 123400 => 6, 0.1234 => 4, 0.001234 => 6
4739 502 100       1375 my $zeros = -$x->{_p} + $cad;
4740             $es .= $dot.'0' x $zeros if $zeros > 0;
4741 6003         80908 }
4742             $es;
4743             }
4744              
4745             # Decimal notation, e.g., "12345.6789" (no exponent).
4746              
4747 48 50   48 1 149 sub bdstr {
4748             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4749 48 50       112  
4750             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4751              
4752             # Inf and NaN
4753 48 100 100     157  
4754 9 100       39 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
4755 7         37 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
4756             return 'inf'; # +inf
4757             }
4758              
4759             # Upgrade?
4760 39 50 33     137  
4761             return $upgrade -> bdstr($x, @r)
4762             if defined($upgrade) && !$x -> isa(__PACKAGE__);
4763              
4764             # Finite number
4765 39         112  
4766 39         78 my $mant = $LIB->_str($x->{_m});
4767 39         99 my $esgn = $x->{_es};
4768             my $eabs = $LIB -> _num($x->{_e});
4769 39         66  
4770             my $uintmax = ~0;
4771 39         66  
4772 39 50       112 my $str = $mant;
4773             if ($esgn eq '+') {
4774 39 50       86  
4775             croak("The absolute value of the exponent is too large")
4776             if $eabs > $uintmax;
4777 39         86  
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       222  
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 10452 sub bsstr {
4799             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4800 175 50       404  
4801             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4802              
4803             # Inf and NaN
4804 175 100 100     666  
4805 28 100       224 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
4806 12         106 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
4807             return 'inf'; # +inf
4808             }
4809              
4810             # Upgrade?
4811 147 50 33     442  
4812             return $upgrade -> bsstr($x, @r)
4813             if defined($upgrade) && !$x -> isa(__PACKAGE__);
4814              
4815             # Finite number
4816              
4817 147 100       628 ($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 1381 sub bnstr {
4824             my ($class, $x, @r) = ref($_[0]) ? (ref($_[0]), @_) : objectify(1, @_);
4825 483 50       1015  
4826             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4827              
4828             # Inf and NaN
4829 483 50 66     1349  
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     1107  
4837             return $upgrade -> bnstr($x, @r)
4838             if defined($upgrade) && !$x -> isa(__PACKAGE__);
4839              
4840             # Finite number
4841 483 100       1051  
4842             my $str = $x->{sign} eq '-' ? '-' : '';
4843              
4844             # Get the mantissa and the length of the mantissa.
4845 483         1426  
4846 483         894 my $mant = $LIB->_str($x->{_m});
4847             my $mantlen = CORE::length($mant);
4848 483 100       1039  
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         310  
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         1383 my ($eabs, $esgn) = $LIB -> _sadd($LIB -> _copy($x->{_e}), $x->{_es},
4863 418         1412 $LIB -> _new($mantlen - 1), "+");
4864 418         1371 substr $mant, 1, 0, ".";
4865             $str .= $mant . 'e' . $esgn . $LIB->_str($eabs);
4866              
4867             }
4868 483         2792  
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 502 # 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     172  
4970 12 100       161 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
4971 4         45 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
4972             return 'inf'; # +inf
4973             }
4974              
4975             # Upgrade?
4976 24 50 33     63  
4977             return $upgrade -> to_hex($x, @r)
4978             if defined($upgrade) && !$x -> isa(__PACKAGE__);
4979              
4980             # Finite number
4981 24 100       104  
4982             return '0' if $x->is_zero();
4983 16 50       66  
4984             return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex?
4985 16         59  
4986 16 50       59 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         106 }
4990 16 100       183 my $str = $LIB->_to_hex($z);
4991             return $x->{sign} eq '-' ? "-$str" : $str;
4992             }
4993              
4994             sub to_oct {
4995 40 50   40 1 548 # 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       108  
4998             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
4999              
5000             # Inf and NaN
5001 40 100 100     174  
5002 12 100       147 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
5003 4         43 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
5004             return 'inf'; # +inf
5005             }
5006              
5007             # Upgrade?
5008 28 50 33     299  
5009             return $upgrade -> to_hex($x, @r)
5010             if defined($upgrade) && !$x -> isa(__PACKAGE__);
5011              
5012             # Finite number
5013 28 100       78  
5014             return '0' if $x->is_zero();
5015 20 50       63  
5016             return $nan if $x->{_es} ne '+'; # how to do 1e-1 in octal?
5017 20         68  
5018 20 50       65 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         92 }
5022 20 100       216 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       100  
5030             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
5031              
5032             # Inf and NaN
5033 40 100 100     174  
5034 12 100       111 if ($x->{sign} ne '+' && $x->{sign} ne '-') {
5035 4         42 return $x->{sign} unless $x->{sign} eq '+inf'; # -inf, NaN
5036             return 'inf'; # +inf
5037             }
5038              
5039             # Upgrade?
5040 28 50 33     74  
5041             return $upgrade -> to_hex($x, @r)
5042             if defined($upgrade) && !$x -> isa(__PACKAGE__);
5043              
5044             # Finite number
5045 28 100       69  
5046             return '0' if $x->is_zero();
5047 20 50       69  
5048             return $nan if $x->{_es} ne '+'; # how to do 1e-1 in binary?
5049 20         57  
5050 20 50       55 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         96 }
5054 20 100       242 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 502  
5296             my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
5297 36 50       90  
5298             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
5299 36 100       167  
5300 24 100       77 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
5301             return '0x0' if $x->is_zero();
5302 16 50       52  
5303             return $nan if $x->{_es} ne '+'; # how to do 1e-1 in hex?
5304 16         69  
5305 16 50       61 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         76 }
5309 16 100       198 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 550  
5316             my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
5317 40 50       102  
5318             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
5319 40 100       171  
5320 28 100       78 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
5321             return '00' if $x->is_zero();
5322 20 50       76  
5323             return $nan if $x->{_es} ne '+'; # how to do 1e-1 in octal?
5324 20         81  
5325 20 50       68 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         81 }
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 591  
5336             my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
5337 40 50       96  
5338             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
5339 40 100       179  
5340 28 100       93 return $x->bstr() if $x->{sign} !~ /^[+-]$/; # inf, nan etc
5341             return '0b0' if $x->is_zero();
5342 20 50       74  
5343             return $nan if $x->{_es} ne '+'; # how to do 1e-1 in binary?
5344 20         63  
5345 20 50       58 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       272 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 1694  
5356             my (undef, $x, @r) = ref($_[0]) ? (undef, @_) : objectify(1, @_);
5357 483 50       1080  
5358             carp "Rounding is not supported for ", (caller(0))[3], "()" if @r;
5359 483 50       1238  
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       1177  
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         1395  
5374             return 0 + $x -> bnstr();
5375             }
5376              
5377             ###############################################################################
5378             # Private methods and functions.
5379             ###############################################################################
5380              
5381 55     55   1824 sub import {
5382 55         119 my $class = shift;
5383             $IMPORT++; # remember we did import()
5384 55         724  
5385 55         106 my @import = ('objectify');
5386             my @a; # unrecognized arguments
5387 55         263  
5388 17         39 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       51  
5415 2         15 if ($param eq 'upgrade') {
5416 2         6 $class -> upgrade(shift);
5417             next;
5418             }
5419              
5420             # Downgrading.
5421 15 100       41  
5422 1         11 if ($param eq 'downgrade') {
5423 1         4 $class -> downgrade(shift);
5424             next;
5425             }
5426              
5427             # Accuracy.
5428 14 50       39  
5429 0         0 if ($param eq 'accuracy') {
5430 0         0 $class -> accuracy(shift);
5431             next;
5432             }
5433              
5434             # Precision.
5435 14 50       49  
5436 0         0 if ($param eq 'precision') {
5437 0         0 $class -> precision(shift);
5438             next;
5439             }
5440              
5441             # Rounding mode.
5442 14 50       36  
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       131  
5450 11         33 if ($param =~ /^(lib|try|only)\z/) {
5451 11 50       44 push @import, $param;
5452 11         42 push @import, shift() if @_;
5453             next;
5454             }
5455 3 50       11  
5456             if ($param eq 'with') {
5457             # alternative class for our private parts()
5458             # XXX: no longer supported
5459             # $LIB = shift() || 'Calc';
5460 3         4 # carp "'with' is no longer supported, use 'lib', 'try', or 'only'";
5461 3         11 shift;
5462             next;
5463             }
5464              
5465             # Unrecognized parameter.
5466 0         0  
5467             push @a, $param;
5468             }
5469 55         341  
5470             Math::BigInt -> import(@import);
5471              
5472 55         365 # find out which one was actually loaded
5473             $LIB = Math::BigInt -> config('lib');
5474 55         64384  
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   7 # at least D digits long. D should be at least 50.
5481             my $d = shift;
5482              
5483 3         8 # 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         6 # 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       12 # bignum;" :(
5493 3 50       9 $l = $l->numify if ref($l);
5494 3 50       10 $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         10 # lg(n!)?)
5500 15         28 while ($r - $l > 1) {
5501 15         51 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       39 / $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   372 # Modifies $x in place.
5513 131         314 my ($x, $scale) = @_;
5514             my $class = ref $x;
5515              
5516 131 100       344 # 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         318  
5539             my ($v, $u, $numer, $denom, $factor, $f);
5540 108         319  
5541 108         652 $v = $x -> copy();
5542 108         703 $v = $v -> binc(); # v = x+1
5543 108         543 $x = $x -> bdec();
5544             $u = $x -> copy(); # u = x-1; x = x-1
5545 108         639  
5546             $x = $x -> bdiv($v, $scaleup); # first term: u/v
5547 108         542  
5548 108         582 $numer = $u -> copy(); # numerator
5549             $denom = $v -> copy(); # denominator
5550 108         513  
5551 108         605 $u = $u -> bmul($u); # u^2
5552             $v = $v -> bmul($v); # v^2
5553 108         616  
5554 108         648 $numer = $numer -> bmul($u); # u^3
5555             $denom = $denom -> bmul($v); # v^3
5556 108         636  
5557 108         669 $factor = $class -> new(3);
5558             $f = $class -> new(2);
5559 108         442  
5560 3114         7860 while (1) {
5561             my $next = $numer -> copy() -> bround($scaleup)
5562             -> bdiv($denom -> copy() -> bmul($factor) -> bround($scaleup), $scaleup);
5563 3114         12202  
5564 3114         5304 $next->{_a} = undef;
5565 3114         7131 $next->{_p} = undef;
5566 3114         7938 my $x_prev = $x -> copy();
5567             $x = $x -> badd($next);
5568 3114 100       8153  
5569             last if $x -> bacmp($x_prev) == 0;
5570              
5571 3006         8048 # calculate things for the next term
5572 3006         7099 $numer = $numer -> bmul($u);
5573 3006         7880 $denom = $denom -> bmul($v);
5574             $factor = $factor -> badd($f);
5575             }
5576 108         627  
5577 108         672 $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   459 # and then "correcting" the result to the proper one. Modifies $x in place.
5584 168         337 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         614  
5614 168 100       770 my $dbd = $LIB->_num($x->{_e});
5615 168         627 $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         408  
5621             my $calc = 1; # do some calculation?
5622              
5623             # No upgrading or downgrading in the intermediate computations.
5624 168         368  
5625 168         320 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     858 # 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       29 # we can use the cached value in these cases
5636 7         38 if ($scale <= $LOG_10_A) {
5637 7         30 $x = $x->bzero();
5638 7         29 $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     550 # 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         103 {
5647             $dbd = 0; # disable shortcut
5648 32 50       146 # we can use the cached value in these cases
5649 32         111 if ($scale <= $LOG_2_A) {
5650 32         225 $x = $x->bzero();
5651 32         146 $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     1178 # 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       7 # we can use the cached value in these cases
5666 2         8 if ($scale <= $LOG_10_A) {
5667 2         12 $x = $x->bzero();
5668 2         13 $x = $x->bsub($LOG_10);
5669             $calc = 0; # no need to calc, but round
5670             }
5671             }
5672 168 100       553  
5673             return $x if $calc == 0; # already have the result
5674              
5675 127         308 # 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         390  
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     1245 # 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       267 # 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       217  
5693             if ($scale <= $LOG_10_A) {
5694 47         150 # 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       6 # 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         4 # log(1.25) + log(2) + log(2) + log(2):
5722 1         4 $l_10 = $l_10->badd($l_2);
5723 1         4 $l_10 = $l_10->badd($l_2);
5724 1         4 $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       239 }
5729 48         173 $dbd-- if ($dbd > 1); # 20 => dbd=2, so make it dbd=1
5730 48         264 $l_10 = $l_10->bmul($class->new($dbd)); # log(10) * (digits_before_dot-1)
5731 48 100       234 my $dbd_sign = '+';
5732 7         28 if ($dbd < 0) {
5733 7         26 $dbd = -$dbd;
5734             $dbd_sign = '-';
5735             }
5736 48         256 ($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       570  
5745             $HALF = $class->new($HALF) unless ref($HALF);
5746 127         297  
5747 127         435 my $twos = 0; # default: none (0 times)
5748 40         104 while ($x->bacmp($HALF) <= 0) { # X <= 0.5
5749 40         172 $twos--;
5750             $x = $x->bmul($two);
5751 127         450 }
5752 52         169 while ($x->bacmp($two) >= 0) { # X >= 2
5753 52         234 $twos++;
5754             $x = $x->bdiv($two, $scale+4); # keep all digits
5755 127         657 }
5756             $x = $x->bround($scale+4);
5757             # $twos > 0 => did mul 2, < 0 => did div 2 (but we never did both)
5758 127 100       616 # So calculate correction factor based on ln(2):
5759 72 100       380 if ($twos != 0) {
5760 72 100       278 $LOG_2 = $class->new($LOG_2, undef, undef) unless ref $LOG_2;
5761             if ($scale <= $LOG_2_A) {
5762 69         192 # 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         40 $l_2 = $two->copy();
5767 3         13 $l_2 = $l_2->_log($scale); # scale+4, actually
5768             $LOG_2 = $l_2->copy(); # cache the result for later
5769 3         40 # the copy() is for mul below
5770             $LOG_2_A = $scale;
5771 72         343 }
5772             $l_2 = $l_2->bmul($twos); # * -2 => subtract, * 2 => add
5773 55         165 } else {
5774             undef $l_2;
5775             }
5776 127         812  
5777 127 100       673 $x = $x->_log($scale); # need to do the "normal" way
5778 127 100       660 $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         1382 # all done, $x contains now the result
5782             $x;
5783             }
5784              
5785             sub _pow {
5786 114     114   461 # Calculate a power where $y is a non-integer, like 2 ** 0.3
5787 114         263 my ($x, $y, @r) = @_;
5788             my $class = ref($x);
5789              
5790 114 100       387 # if $y == 0.5, it is sqrt($x)
5791 114 100       378 $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         199 # we need to limit the accuracy to protect against overflow
5804 69         166 my $fallback = 0;
5805 69         275 my ($scale, @params);
5806             ($x, @params) = $x->_find_round_parameters(@r);
5807 69 50       310  
5808             return $x if $x->is_nan(); # error in _find_round_parameters?
5809              
5810 69 100       288 # no rounding at all, so must use fallback
5811             if (scalar @params == 0) {
5812 52         205 # simulate old behaviour
5813 52         165 $params[0] = $class->div_scale(); # and round to it as accuracy
5814 52         119 $params[1] = undef; # disable P
5815 52         116 $scale = $params[0]+4; # at least four more for proper round
5816 52         114 $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     68 # 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   531 # disable them and later re-enable them
  43         158  
  43         23289  
5826 69         215 no strict 'refs';
5827 69         204 my $abr = "$class\::accuracy";
5828 69         189 my $ab = $$abr;
5829 69         189 $$abr = undef;
5830 69         155 my $pbr = "$class\::precision";
5831 69         165 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         193 # them already into account), since these would interfere, too
5835 69         162 $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         159  
5842 69         134 local $Math::BigInt::upgrade = undef;
5843             local $Math::BigFloat::downgrade = undef;
5844 69         178  
5845             my ($limit, $v, $u, $below, $factor, $next, $over);
5846 69         261  
5847 69         479 $u = $x->copy()->blog(undef, $scale)->bmul($y);
5848 69 100       372 my $do_invert = ($u->{sign} eq '-');
5849 69         481 $u = $u->bneg() if $do_invert;
5850 69         361 $v = $class->bone(); # 1
5851 69         501 $factor = $class->new(2); # 2
5852             $x = $x->bone(); # first term: 1
5853 69         378  
5854 69         397 $below = $v->copy();
5855             $over = $u->copy();
5856 69         609  
5857 69         279 $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         8021 # anymore, so we stop:
5862 3277 100       14539 $next = $over->copy()->bdiv($below, $scale);
5863 3208         8283 last if $next->bacmp($limit) <= 0;
5864             $x = $x->badd($next);
5865 3208         10335 # calculate things for the next term
5866 3208         7993 $over *= $u;
5867 3208         8157 $below *= $factor;
5868             $factor = $factor->binc();
5869 3208 50       12660  
5870             last if $x->{sign} !~ /^[-+]$/;
5871             }
5872 69 100       465  
5873 32         168 if ($do_invert) {
5874 32         265 my $x_copy = $x->copy();
5875             $x = $x->bone->bdiv($x_copy, $scale);
5876             }
5877              
5878 69 50       484 # shortcut to not run through _find_round_parameters again
5879 69         323 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       484 }
5884             if ($fallback) {
5885 52         162 # clear a/p after round, since user did not request it
5886 52         130 $x->{_a} = undef;
5887             $x->{_p} = undef;
5888             }
5889 69         286 # restore globals
5890 69         205 $$abr = $ab;
5891 69         2147 $$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__