view release on metacpan or  search on metacpan

cpan/Math-BigInt/lib/Math/BigInt/Calc.pm  view on Meta::CPAN

        $step = [$step];
        while ($c->_acmp($step, $n) <= 0) {
            if ($cx->[0] == 0) {
                $zero_elements ++;
                shift @$cx;
            }
            $c->_mul($cx, $step);
            $c->_inc($step);
        }
    } else {
        # Yes, so we can speed it up slightly

        #    print "# left over steps $n\n";

        my $base_4 = int(sqrt(sqrt($BASE))) - 2;
        #print STDERR "base_4: $base_4\n";
        my $n4 = $n - 4;
        while ($step < $n4 && $step < $base_4) {
            if ($cx->[0] == 0) {
                $zero_elements ++;
                shift @$cx;
            }
            my $b = $step * ($step + 1);
            $step += 2;
            $b *= $step * ($step + 1);
            $step += 2;
            $c->_mul($cx, [$b]);
        }
        my $base_2 = int(sqrt($BASE)) - 1;
        my $n2 = $n - 2;
        #print STDERR "base_2: $base_2\n";
        while ($step < $n2 && $step < $base_2) {
            if ($cx->[0] == 0) {
                $zero_elements ++;
                shift @$cx;
            }
            my $b = $step * ($step + 1);
            $step += 2;
            $c->_mul($cx, [$b]);
        }
        # do what's left over
        while ($step <= $n) {
            $c->_mul($cx, [$step]);
            $step++;
            if ($cx->[0] == 0) {
                $zero_elements ++;
                shift @$cx;
            }
        }
    }
    # multiply in the zeros again
    unshift @$cx, (0) x $zero_elements;
    $cx;                        # return result
}

sub _log_int {
    # calculate integer log of $x to base $base
    # ref to array, ref to array - return ref to array
    my ($c, $x, $base) = @_;

    # X == 0 => NaN
    return if @$x == 1 && $x->[0] == 0;

    # BASE 0 or 1 => NaN
    return if @$base == 1 && $base->[0] < 2;

    # X == 1 => 0 (is exact)
    if (@$x == 1 && $x->[0] == 1) {
        @$x = 0;
        return $x, 1;
    }

    my $cmp = $c->_acmp($x, $base);

    # X == BASE => 1 (is exact)
    if ($cmp == 0) {
        @$x = 1;
        return $x, 1;
    }

    # 1 < X < BASE => 0 (is truncated)
    if ($cmp < 0) {
        @$x = 0;
        return $x, 0;
    }

    my $x_org = $c->_copy($x);  # preserve x

    # Compute a guess for the result based on:
    # $guess = int ( length_in_base_10(X) / ( log(base) / log(10) ) )
    my $len = $c->_len($x_org);
    my $log = log($base->[-1]) / log(10);

    # for each additional element in $base, we add $BASE_LEN to the result,
    # based on the observation that log($BASE, 10) is BASE_LEN and
    # log(x*y) == log(x) + log(y):
    $log += (@$base - 1) * $BASE_LEN;

    # calculate now a guess based on the values obtained above:
    my $res = $c->_new(int($len / $log));

    @$x = @$res;
    my $trial = $c->_pow($c->_copy($base), $x);
    my $acmp = $c->_acmp($trial, $x_org);

    # Did we get the exact result?

    return $x, 1 if $acmp == 0;

    # Too small?

    while ($acmp < 0) {
        $c->_mul($trial, $base);
        $c->_inc($x);
        $acmp = $c->_acmp($trial, $x_org);
    }

    # Too big?

    while ($acmp > 0) {
        $c->_div($trial, $base);
        $c->_dec($x);
        $acmp = $c->_acmp($trial, $x_org);
    }

    return $x, 1 if $acmp == 0;         # result is exact
    return $x, 0;                       # result is too small
}

sub _ilog2 {
    # calculate int(log2($x))

    # There is virtually nothing to gain from computing this any differently
    # than _log_int(), but it is important that we don't use the method
    # inherited from the parent, because that method is very slow for backend
    # libraries whose internal representation uses base 10.

    my ($c, $x) = @_;
    ($x, my $is_exact) = $c -> _log_int($x, $c -> _two());
    return wantarray ? ($x, $is_exact) : $x;
}

sub _ilog10 {
    # calculate int(log10($x))

    my ($c, $x) = @_;

    # X == 0 => NaN
    return if @$x == 1 && $x->[0] == 0;

    # X == 1 => 0 (is exact)
    if (@$x == 1 && $x->[0] == 1) {
        @$x = 0;
        return wantarray ? ($x, 1) : $x;
    }

    my $x_orig = $c -> _copy($x);
    my $nm1 = $c -> _len($x) - 1;

    my $xtmp = $c -> _new($nm1);
    @$x = @$xtmp;

    return $x unless wantarray;

    # See if the original $x is an exact power of 10, in which case all but the
    # most significan chunks are 0, and the most significant chunk is a power
    # of 10.

    my $is_pow10 = 1;
    for my $i (0 .. $#$x_orig - 1) {
        last unless $is_pow10 = $x_orig->[$i] == 0;
    }
    $is_pow10 &&= $x_orig->[-1] == 10**int(0.5 + log($x_orig->[-1]) / log(10));

    return wantarray ? ($x, 1) : $x if $is_pow10;
    return wantarray ? ($x, 0) : $x;
}

sub _clog2 {
    # calculate ceil(log2($x))

    my ($c, $x) = @_;

    # X == 0 => NaN

    return if @$x == 1 && $x->[0] == 0;

    # X == 1 => 0 (is exact)

    if (@$x == 1 && $x->[0] == 1) {
        @$x = 0;
        return wantarray ? ($x, 1) : $x;
    }

    my $base = $c -> _two();
    my $acmp = $c -> _acmp($x, $base);

    # X == BASE => 1 (is exact)

    if ($acmp == 0) {
        @$x = 1;
        return wantarray ? ($x, 1) : $x;
    }

    # 1 < X < BASE => 0 (is truncated)

    if ($acmp < 0) {
        @$x = 0;
        return wantarray ? ($x, 0) : $x;
    }

    # Compute a guess for the result based on:
    # $guess = int( length_in_base_10(X) / (log(base) / log(10)) )

    my $len    = $c -> _len($x);
    my $log    = log(2) / log(10);
    my $guess  = $c -> _new(int($len / $log));
    my $x_orig = $c -> _copy($x);
    @$x = @$guess;

    my $trial = $c -> _pow($c -> _copy($base), $x);
    $acmp = $c -> _acmp($trial, $x_orig);

    # Too big?

    while ($acmp > 0) {
        $c -> _div($trial, $base);
        $c -> _dec($x);
        $acmp = $c -> _acmp($trial, $x_orig);
    }

    # Too small?

    while ($acmp < 0) {
        $c -> _mul($trial, $base);
        $c -> _inc($x);
        $acmp = $c -> _acmp($trial, $x_orig);
    }

    return wantarray ? ($x, 1) : $x if $acmp == 0;      # result is exact
    return wantarray ? ($x, 0) : $x;                    # result is too small
}

sub _clog10 {
    # calculate ceil(log2($x))
    my ($c, $x) = @_;

    # X == 0 => NaN
    return if @$x == 1 && $x->[0] == 0;

    # X == 1 => 0 (is exact)
    if (@$x == 1 && $x->[0] == 1) {
        @$x = 0;
        return wantarray ? ($x, 1) : $x;
    }

    # Get the number of base 10 digits. $n is the desired output, except when
    # $x is an exact power of 10, in which case $n is 1 too big.

    my $n = $c -> _len($x);

    # See if $x is an exact power of 10, in which case all but the most
    # significan chunks are 0, and the most significant chunk is a power of 10.

    my $is_pow10 = 1;
    for my $i (0 .. $#$x - 1) {
        last unless $is_pow10 = $x->[$i] == 0;
    }
    $is_pow10 &&= $x->[-1] == 10**int(0.5 + log($x->[-1]) / log(10));

    $n-- if $is_pow10;

    my $xtmp = $c ->_new($n);
    @$x = @$xtmp;

    return wantarray ? ($x, 1) : $x if $is_pow10;       # result is exact
    return wantarray ? ($x, 0) : $x;                    # result is too small
}

# for debugging:
use constant DEBUG => 0;
my $steps = 0;
sub steps { $steps };

sub _sqrt {
    # square-root of $x in-place

    my ($c, $x) = @_;

    if (@$x == 1) {
        # fits into one Perl scalar, so result can be computed directly
        $x->[0] = int(sqrt($x->[0]));
        return $x;
    }

    # Create an initial guess for the square root.

    my $s;
    if (@$x % 2) {
        $s = [ (0) x ((@$x - 1) / 2), int(sqrt($x->[-1])) ];
    } else {
        $s = [ (0) x ((@$x - 2) / 2), int(sqrt($x->[-2] + $x->[-1] * $BASE)) ];
    }

    # Newton's method for the square root of y:
    #
    #                      x(n) * x(n) - y
    #     x(n+1) = x(n) - -----------------

cpan/Math-BigInt/lib/Math/BigInt/Calc.pm  view on Meta::CPAN

    }
    $x;
}

sub _from_bin {
    # convert a hex number to decimal (string, return ref to array)
    my ($c, $bs) = @_;

    # instead of converting X (8) bit at a time, it is faster to "convert" the
    # number to hex, and then call _from_hex.

    my $hs = $bs;
    $hs =~ s/^[+-]?0b//;                                # remove sign and 0b
    my $l = length($hs);                                # bits
    $hs = '0' x (8 - ($l % 8)) . $hs if ($l % 8) != 0;  # padd left side w/ 0
    my $h = '0x' . unpack('H*', pack ('B*', $hs));      # repack as hex

    $c->_from_hex($h);
}

##############################################################################
# special modulus functions

sub _modinv {

    # modular multiplicative inverse
    my ($c, $x, $y) = @_;

    # modulo zero
    if ($c->_is_zero($y)) {
        return;
    }

    # modulo one
    if ($c->_is_one($y)) {
        return $c->_zero(), '+';
    }

    my $u = $c->_zero();
    my $v = $c->_one();
    my $a = $c->_copy($y);
    my $b = $c->_copy($x);

    # Euclid's Algorithm for bgcd(), only that we calc bgcd() ($a) and the result
    # ($u) at the same time. See comments in BigInt for why this works.
    my $q;
    my $sign = 1;
    {
        ($a, $q, $b) = ($b, $c->_div($a, $b));          # step 1
        last if $c->_is_zero($b);

        my $t = $c->_add(                               # step 2:
                         $c->_mul($c->_copy($v), $q),   #  t =   v * q
                         $u);                           #      + u
        $u = $v;                                        #  u = v
        $v = $t;                                        #  v = t
        $sign = -$sign;
        redo;
    }

    # if the gcd is not 1, then return NaN
    return unless $c->_is_one($a);

    ($v, $sign == 1 ? '+' : '-');
}

sub _modpow {
    # modulus of power ($x ** $y) % $z
    my ($c, $num, $exp, $mod) = @_;

    # a^b (mod 1) = 0 for all a and b
    if ($c->_is_one($mod)) {
        @$num = 0;
        return $num;
    }

    # 0^a (mod m) = 0 if m != 0, a != 0
    # 0^0 (mod m) = 1 if m != 0
    if ($c->_is_zero($num)) {
        if ($c->_is_zero($exp)) {
            @$num = 1;
        } else {
            @$num = 0;
        }
        return $num;
    }

    # We could do the following, but it doesn't actually save any time. The
    # _copy() is needed in case $num and $mod are the same object.
    #$num = $c->_mod($c->_copy($num), $mod);

    my $acc = $c->_copy($num);
    my $t = $c->_one();

    my $expbin = $c->_to_bin($exp);
    my $len = length($expbin);
    while ($len--) {
        if (substr($expbin, $len, 1) eq '1') {  # if odd
            $t = $c->_mul($t, $acc);
            $t = $c->_mod($t, $mod);
        }
        $acc = $c->_mul($acc, $acc);
        $acc = $c->_mod($acc, $mod);
    }
    @$num = @$t;
    $num;
}

sub _gcd {
    # Greatest common divisor.

    my ($c, $x, $y) = @_;

    # gcd(0, 0) = 0
    # gcd(0, a) = a, if a != 0

    if (@$x == 1 && $x->[0] == 0) {
        if (@$y == 1 && $y->[0] == 0) {
            @$x = 0;
        } else {
            @$x = @$y;