view release on metacpan or  search on metacpan

lib/Math/BigInt/Lib.pm  view on Meta::CPAN

        while ($class -> _acmp($r, $y) >= 0) {
            $r = $class -> _sub($r, $y);
        }
        return $r;
    }
}

##############################################################################
# shifts

sub _rsft {
    my ($class, $x, $n, $b) = @_;
    $b = $class -> _new($b) unless ref $b;
    return scalar $class -> _div($x, $class -> _pow($class -> _copy($b), $n));
}

sub _lsft {
    my ($class, $x, $n, $b) = @_;
    $b = $class -> _new($b) unless ref $b;
    return $class -> _mul($x, $class -> _pow($class -> _copy($b), $n));
}

sub _pow {
    # power of $x to $y
    my ($class, $x, $y) = @_;

    if ($class -> _is_zero($y)) {
        return $class -> _one();        # y == 0 => x => 1
    }

    if (($class -> _is_one($x)) ||      #    x == 1
        ($class -> _is_one($y)))        # or y == 1
    {
        return $x;
    }

    if ($class -> _is_zero($x)) {
        return $class -> _zero();       # 0 ** y => 0 (if not y <= 0)
    }

    my $pow2 = $class -> _one();

    my $y_bin = $class -> _as_bin($y);
    $y_bin =~ s/^0b//;
    my $len = length($y_bin);

    while (--$len > 0) {
        $pow2 = $class -> _mul($pow2, $x) if substr($y_bin, $len, 1) eq '1';
        $x = $class -> _mul($x, $x);
    }

    $x = $class -> _mul($x, $pow2);
    return $x;
}

sub _nok {
    # Return binomial coefficient (n over k).
    my ($class, $n, $k) = @_;

    # If k > n/2, or, equivalently, 2*k > n, compute nok(n, k) as
    # nok(n, n-k), to minimize the number if iterations in the loop.

    {
        my $twok = $class -> _mul($class -> _two(), $class -> _copy($k));
        if ($class -> _acmp($twok, $n) > 0) {
            $k = $class -> _sub($class -> _copy($n), $k);
        }
    }

    # Example:
    #
    # / 7 \       7!       1*2*3*4 * 5*6*7   5 * 6 * 7
    # |   | = --------- =  --------------- = --------- = ((5 * 6) / 2 * 7) / 3
    # \ 3 /   (7-3)! 3!    1*2*3*4 * 1*2*3   1 * 2 * 3
    #
    # Equivalently, _nok(11, 5) is computed as
    #
    # (((((((7 * 8) / 2) * 9) / 3) * 10) / 4) * 11) / 5

    if ($class -> _is_zero($k)) {
        return $class -> _one();
    }

    # Make a copy of the original n, in case the subclass modifies n in-place.

    my $n_orig = $class -> _copy($n);

    # n = 5, f = 6, d = 2 (cf. example above)

    $n = $class -> _sub($n, $k);
    $n = $class -> _inc($n);

    my $f = $class -> _copy($n);
    $f = $class -> _inc($f);

    my $d = $class -> _two();

    # while f <= n (the original n, that is) ...

    while ($class -> _acmp($f, $n_orig) <= 0) {
        $n = $class -> _mul($n, $f);
        $n = $class -> _div($n, $d);
        $f = $class -> _inc($f);
        $d = $class -> _inc($d);
    }

    return $n;
}

#sub _fac {
#    # factorial
#    my ($class, $x) = @_;
#
#    my $two = $class -> _two();
#
#    if ($class -> _acmp($x, $two) < 0) {
#        return $class -> _one();
#    }
#
#    my $i = $class -> _copy($x);
#    while ($class -> _acmp($i, $two) > 0) {

lib/Math/BigInt/Lib.pm  view on Meta::CPAN

    my $DEBUG = 0;

    # Split y into mantissa and exponent in base 10, so that
    #
    #   y = xm * 10^xe, where 0 < xm < 1 and xe is an integer

    my $y_str  = $class -> _str($y);
    my $ym = "." . $y_str;
    my $ye = length($y_str);

    # From this compute the approximate base 10 logarithm of y
    #
    #   log_10(y) = log_10(ym) + log_10(ye^10)
    #             = log(ym)/log(10) + ye

    my $log10y = log($ym) / log(10) + $ye;

    # And from this compute the approximate base 10 logarithm of x, where
    # x = y^(1/n)
    #
    #   log_10(x) = log_10(y)/n

    my $log10x = $log10y / $class -> _num($n);

    # From this compute xm and xe, the mantissa and exponent (in base 10) of x,
    # where 1 < xm <= 10 and xe is an integer.

    my $xe = int $log10x;
    my $xm = 10 ** ($log10x - $xe);

    # Scale the mantissa and exponent to increase the integer part of ym, which
    # gives us better accuracy.

    if ($DEBUG) {
        print "\n";
        print "y_str  = $y_str\n";
        print "ym     = $ym\n";
        print "ye     = $ye\n";
        print "log10y = $log10y\n";
        print "log10x = $log10x\n";
        print "xm     = $xm\n";
        print "xe     = $xe\n";
    }

    my $d = $xe < 15 ? $xe : 15;
    $xm *= 10 ** $d;
    $xe -= $d;

    if ($DEBUG) {
        print "\n";
        print "xm     = $xm\n";
        print "xe     = $xe\n";
    }

    # If the mantissa is not an integer, round up to nearest integer, and then
    # convert the number to a string. It is important to always round up due to
    # how Newton's method behaves in this case. If the initial guess is too
    # small, the next guess will be too large, after which every succeeding
    # guess converges the correct value from above. Now, if the initial guess
    # is too small and n is large, the next guess will be much too large and
    # require a large number of iterations to get close to the solution.
    # Because of this, we are likely to find the solution faster if we make
    # sure the initial guess is not too small.

    my $xm_int = int($xm);
    my $x_str = sprintf '%.0f', $xm > $xm_int ? $xm_int + 1 : $xm_int;
    $x_str .= "0" x $xe;

    my $x = $class -> _new($x_str);

    if ($DEBUG) {
        print "xm     = $xm\n";
        print "xe     = $xe\n";
        print "\n";
        print "x_str  = $x_str (initial guess)\n";
        print "\n";
    }

    # Use Newton's method for computing n'th root of y.
    #
    # x(i+1) = x(i) - f(x(i)) / f'(x(i))
    #        = x(i) - (x(i)^n - y) / (n * x(i)^(n-1))   # use if x(i)^n > y
    #        = x(i) + (y - x(i)^n) / (n * x(i)^(n-1))   # use if x(i)^n < y

    # Determine if x, our guess, is too small, correct, or too large. Rather
    # than computing x(i)^n and x(i)^(n-1) directly, compute x(i)^(n-1) and
    # then the same value multiplied by x.

    my $nm1     = $class -> _dec($class -> _copy($n));           # n-1
    my $xpownm1 = $class -> _pow($class -> _copy($x), $nm1);     # x(i)^(n-1)
    my $xpown   = $class -> _mul($class -> _copy($xpownm1), $x); # x(i)^n
    my $acmp    = $class -> _acmp($xpown, $y);                   # x(i)^n <=> y

    if ($DEBUG) {
        print "\n";
        print "x      = ", $class -> _str($x), "\n";
        print "x^n    = ", $class -> _str($xpown), "\n";
        print "y      = ", $class -> _str($y), "\n";
        print "acmp   = $acmp\n";
    }

    # If x is too small, do one iteration of Newton's method. Since the
    # function f(x) = x^n - y is concave and monotonically increasing, the next
    # guess for x will either be correct or too large.

    if ($acmp < 0) {

        # x(i+1) = x(i) + (y - x(i)^n) / (n * x(i)^(n-1))

        my $numer = $class -> _sub($class -> _copy($y), $xpown);    # y - x(i)^n
        my $denom = $class -> _mul($class -> _copy($n), $xpownm1);  # n * x(i)^(n-1)
        my $delta = $class -> _div($numer, $denom);

        if ($DEBUG) {
            print "\n";
            print "numer  = ", $class -> _str($numer), "\n";
            print "denom  = ", $class -> _str($denom), "\n";
            print "delta  = ", $class -> _str($delta), "\n";
        }

        unless ($class -> _is_zero($delta)) {