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