view release on metacpan or  search on metacpan

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

#
## Ratmod
#

sub ratmod ($$) {
    my ($n, $m) = @_;

    $n = _star2obj($n) // goto &nan;
    $m = _star2mpz($m) // goto &nan;

    my $r = Math::GMPz::Rmpz_init();

    if (ref($n) eq 'Math::GMPz') {
        Math::GMPz::Rmpz_mod($r, $n, $m);
    }
    else {
        $r = _modular_rational($n, $m) // goto &nan;
        Math::GMPz::Rmpz_mod($r, $r, $m);
    }

    bless \$r;
}

#
## Powmod
#

sub powmod ($$$) {
    my ($n, $k, $m) = @_;

    $n = _star2obj($n) // goto &nan;
    $k = _star2mpz($k) // goto &nan;
    $m = _star2mpz($m) // goto &nan;

    Math::GMPz::Rmpz_sgn($m) || goto &nan;

    if (ref($n) ne 'Math::GMPz') {
        if (__is_int__($n)) {
            $n = _any2mpz($n) // goto &nan;
        }
        else {
            $n = _modular_rational($n, $m) // goto &nan;
        }
    }

    my $r = Math::GMPz::Rmpz_init();

    if (Math::GMPz::Rmpz_sgn($k) < 0) {
        Math::GMPz::Rmpz_invert($r, $n, $m) or goto &nan;
    }

    Math::GMPz::Rmpz_fits_ulong_p($k)
      ? Math::GMPz::Rmpz_powm_ui($r, $n, Math::GMPz::Rmpz_get_ui($k), $m)
      : Math::GMPz::Rmpz_powm($r, $n, $k, $m);

    bless \$r;
}

#
## Geometric summation formula
## https://en.wikipedia.org/wiki/Geometric_series
#

sub geometric_sum ($$) {
    my ($n, $r) = @_;

    $n = ref($n) eq __PACKAGE__ ? $$n : _star2obj($n);
    $r = ref($r) eq __PACKAGE__ ? $$r : _star2obj($r);

    bless \__div__(__sub__(__pow__($r, __add__($n, 1)), 1), __sub__($r, 1));
}

#
## Faulhaber's summation formula
## https://en.wikipedia.org/wiki/Faulhaber%27s_formula
#

sub faulhaber_sum ($$) {
    my ($n, $p) = @_;

    $n = _star2mpz($n) // goto &nan;
    $p = _star2ui($p)  // goto &nan;

    if ($p == 0) {
        return bless \$n;
    }

    if ($p == 1 or $p == 3) {
        my $r = Math::GMPz::Rmpz_init();
        Math::GMPz::Rmpz_add_ui($r, $n, 1);
        Math::GMPz::Rmpz_mul($r, $r, $n);
        Math::GMPz::Rmpz_div_2exp($r, $r, 1);
        Math::GMPz::Rmpz_mul($r, $r, $r) if ($p == 3);
        return bless \$r;
    }

    state $z = Math::GMPz::Rmpz_init_nobless();

    if ($p == 2) {    # n*(n+1)*(2*n+1)/6
        my $r = Math::GMPz::Rmpz_init();
        Math::GMPz::Rmpz_add_ui($z, $n, 1);
        Math::GMPz::Rmpz_mul($r, $z, $n);
        Math::GMPz::Rmpz_mul_2exp($z, $z, 1);
        Math::GMPz::Rmpz_sub_ui($z, $z, 1);
        Math::GMPz::Rmpz_mul($r, $r, $z);
        Math::GMPz::Rmpz_divexact_ui($r, $r, 6);
        return bless \$r;
    }

    # When p >= n, sum the powers directly.
    if (Math::GMPz::Rmpz_cmp_ui($n, $p) <= 0) {
        my $r = Math::GMPz::Rmpz_init_set_ui(0);
        foreach my $k (1 .. Math::GMPz::Rmpz_get_ui($n)) {
            Math::GMPz::Rmpz_ui_pow_ui($z, $k, $p);
            Math::GMPz::Rmpz_add($r, $r, $z);
        }
        return bless \$r;
    }

    my @B = _bernoulli_numbers($p);

    my $q = Math::GMPq::Rmpq_init();
    my $u = Math::GMPz::Rmpz_init_set_ui(1);

    my $sum = Math::GMPq::Rmpq_init();
    Math::GMPq::Rmpq_set_ui($sum, 0, 1);

    # Sum_{k=1..n} k^p = 1/(p+1) * Sum_{j=0..p} binomial(p+1, j) * n^(p-j+1) * bernoulli(j)
    #                  = 1/(p+1) * Sum_{j=0..p} binomial(p+1, p-j) * n^(j+1) * bernoulli(p-j)

    foreach my $j (0 .. $p - 2) {

        Math::GMPz::Rmpz_mul($u, $u, $n);

        # Skip when bernoulli(p-j) == 0