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