view release on metacpan or  search on metacpan

lib/Math/BigFloat.pm  view on Meta::CPAN

051870721134999999837297804995105973173281609631859502445945534690830264252
230825334468503526193118817101000313783875288658753320838142061717766914730
359825349042875546873115956286388235378759375195778185778053217122680661300
192787661119590921642019893809525720106548586327886593615338182796823030195
EOF

        # Should we round up?

        my $round_up;

        # From the string above, we need to extract the number of digits we
        # want plus extra characters for the newlines.

        my $nchrs = $n + int($n / 75);

        # Extract the digits we want.

        my $digits = substr($all_digits, 0, $nchrs);

        # Find out whether we should round up or down. Rounding is easy, since
        # pi is trancendental. With directed rounding, it doesn't matter what
        # the following digits are. With rounding to nearest, we only have to
        # look at one extra digit.

        if ($rmode eq 'trunc') {
            $round_up = 0;
        } else {
            my $next_digit = substr($all_digits, $nchrs, 1);
            $round_up = $next_digit lt '5' ? 0 : 1;
        }

        # Remove the newlines.

        $digits =~ tr/0-9//cd;

        # Now do the rounding. We could easily make the regex substitution
        # handle all cases, but we avoid using the regex engine when it is
        # simple to avoid it.

        if ($round_up) {
            my $last_digit = substr($digits, -1, 1);
            if ($last_digit lt '9') {
                substr($digits, -1, 1) = ++$last_digit;
            } else {
                $digits =~ s{([0-8])(9+)$}
                            { ($1 + 1) . ("0" x CORE::length($2)) }e;
            }
        }

        # Convert to an object.

        $pi = bless {
                     sign => '+',
                     _m   => $LIB -> _new($digits),
                     _es  => CORE::length($digits) > 1 ? '-' : '+',
                     _e   => $LIB -> _new($n - 1),
                    }, $class;

    } else {

        # For large accuracy, the arctan formulas become very inefficient with
        # Math::BigFloat, so use Brent-Salamin (aka AGM or Gauss-Legendre).

        # Use a few more digits in the intermediate computations.
        $n += 8;

        $HALF = $class -> new($HALF) unless ref($HALF);
        my ($an, $bn, $tn, $pn)
          = ($class -> bone, $HALF -> copy() -> bsqrt($n),
             $HALF -> copy() -> bmul($HALF), $class -> bone);
        while ($pn < $n) {
            my $prev_an = $an -> copy();
            $an -> badd($bn) -> bmul($HALF, $n);
            $bn -> bmul($prev_an) -> bsqrt($n);
            $prev_an -> bsub($an);
            $tn -> bsub($pn * $prev_an * $prev_an);
            $pn -> badd($pn);
        }
        $an -> badd($bn);
        $an -> bmul($an, $n) -> bdiv(4 * $tn, $n);

        $an -> round(@r);
        $pi = $an;
    }

    if (defined $r[0]) {
        $pi -> accuracy($r[0]);
    } elsif (defined $r[1]) {
        $pi -> precision($r[1]);
    }

    $pi -> _dng() if ($pi -> is_int() ||
                      $pi -> is_inf() ||
                      $pi -> is_nan());

    %$self = %$pi;
    bless $self, ref($pi);
    return $self;
}

sub copy {
    my ($x, $class);
    if (ref($_[0])) {           # $y = $x -> copy()
        $x = shift;
        $class = ref($x);
    } else {                    # $y = Math::BigInt -> copy($y)
        $class = shift;
        $x = shift;
    }

    carp "Rounding is not supported for ", (caller(0))[3], "()" if @_;

    my $copy = bless {}, $class;

    $copy->{sign} = $x->{sign};
    $copy->{_es}  = $x->{_es};
    $copy->{_m}   = $LIB->_copy($x->{_m});
    $copy->{_e}   = $LIB->_copy($x->{_e});

    $copy->{accuracy}  = $x->{accuracy} if exists $x->{accuracy};
    $copy->{precision} = $x->{precision} if exists $x->{precision};

lib/Math/BigFloat.pm  view on Meta::CPAN

    $v = $x -> copy();
    $v = $v -> binc();                  # v = x+1
    $x = $x -> bdec();
    $u = $x -> copy();                  # u = x-1; x = x-1

    $x = $x -> bdiv($v, $scaleup);        # first term: u/v

    $numer = $u -> copy();              # numerator
    $denom = $v -> copy();              # denominator

    $u = $u -> bmul($u);                # u^2
    $v = $v -> bmul($v);                # v^2

    $numer = $numer -> bmul($u);        # u^3
    $denom = $denom -> bmul($v);        # v^3

    $factor = $class -> new(3);
    $f = $class -> new(2);

    while (1) {
        my $next = $numer -> copy() -> bround($scaleup)
          -> bdiv($denom -> copy() -> bmul($factor) -> bround($scaleup), $scaleup);

        $next->{accuracy} = undef;
        $next->{precision} = undef;
        my $x_prev = $x -> copy();
        $x = $x -> badd($next);

        last if $x -> bacmp($x_prev) == 0;

        # calculate things for the next term
        $numer  = $numer -> bmul($u);
        $denom  = $denom -> bmul($v);
        $factor = $factor -> badd($f);
    }

    $x = $x -> bmul($f);             # $x *= 2
    $x = $x -> bround($scale);
}

sub _log_10 {
    # Internal log function based on reducing input to the range of 0.1 .. 9.99
    # and then "correcting" the result to the proper one. Modifies $x in place.
    my ($x, $scale) = @_;
    my $class = ref $x;

    # Taking blog() from numbers greater than 10 takes a *very long* time, so
    # we break the computation down into parts based on the observation that:
    #  blog(X*Y) = blog(X) + blog(Y)
    # We set Y here to multiples of 10 so that $x becomes below 1 - the smaller
    # $x is the faster it gets. Since 2*$x takes about 10 times as long, we
    # make it faster by about a factor of 100 by dividing $x by 10.

    # The same observation is valid for numbers smaller than 0.1, e.g.
    # computing log(1) is fastest, and the further away we get from 1, the
    # longer it takes. So we also 'break' this down by multiplying $x with 10
    # and subtract the log(10) afterwards to get the correct result.

    # To get $x even closer to 1, we also divide by 2 and then use log(2) to
    # correct for this. For instance if $x is 2.4, we use the formula:
    #  blog(2.4 * 2) == blog(1.2) + blog(2)
    # and thus calculate only blog(1.2) and blog(2), which is faster in total
    # than calculating blog(2.4).

    # In addition, the values for blog(2) and blog(10) are cached.

    # Calculate the number of digits before the dot, i.e., 1 + floor(log10(x)):
    #   x = 123      => dbd =  3
    #   x = 1.23     => dbd =  1
    #   x = 0.0123   => dbd = -1
    #   x = 0.000123 => dbd = -3
    #   etc.

    my $dbd = $LIB->_num($x->{_e});
    $dbd = -$dbd if $x->{_es} eq '-';
    $dbd += $LIB->_len($x->{_m});

    # more than one digit (e.g. at least 10), but *not* exactly 10 to avoid
    # infinite recursion

    my $calc = 1;               # do some calculation?

    # No upgrading or downgrading in the intermediate computations.

    my $upg = $class -> upgrade();
    my $dng = $class -> downgrade();
    $class -> upgrade(undef);
    $class -> downgrade(undef);

    # disable the shortcut for 10, since we need log(10) and this would recurse
    # infinitely deep
    if ($x->{_es} eq '+' &&                     # $x == 10
        ($LIB->_is_one($x->{_e}) &&
         $LIB->_is_one($x->{_m})))
    {
        $dbd = 0;               # disable shortcut
        # we can use the cached value in these cases
        if ($scale <= $LOG_10_A) {
            $x = $x -> bzero();
            $x = $x -> badd($LOG_10); # modify $x in place
            $calc = 0;                      # no need to calc, but round
        }
        # if we can't use the shortcut, we continue normally
    } else {
        # disable the shortcut for 2, since we maybe have it cached
        if (($LIB->_is_zero($x->{_e}) &&        # $x == 2
             $LIB->_is_two($x->{_m})))
        {
            $dbd = 0;           # disable shortcut
            # we can use the cached value in these cases
            if ($scale <= $LOG_2_A) {
                $x = $x -> bzero();
                $x = $x -> badd($LOG_2); # modify $x in place
                $calc = 0;                     # no need to calc, but round
            }
            # if we can't use the shortcut, we continue normally
        }
    }

    # if $x = 0.1, we know the result must be 0-log(10)