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