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