Math-BigNum

 view release on metacpan or  search on metacpan

MANIFEST  view on Meta::CPAN

Build.PL
Changes
examples/agm_pi.pl
examples/arithmetic_coding.pl
examples/bernoulli_numbers_recursive.pl
examples/binary_arithmetic_coding.pl
examples/binradix_arithmetic_coding.pl
examples/closed_form_fib.pl
examples/coin_change.pl
examples/computing_pi.pl
examples/faulhaber_s_formula.pl
examples/fibonacci_validation.pl
examples/inverse_of_factorial.pl
examples/is_power.pl
examples/miller_rabin_primality_test.pl
examples/pi_machin.pl
examples/power_pairs.pl
examples/prime_count_approx.pl
examples/rsa_algorithm.pl
examples/tac-compressor.pl
examples/zeta_2n.pl

examples/faulhaber_s_formula.pl  view on Meta::CPAN

#!/usr/bin/perl

# The formula for calculating the sum of consecutive
# numbers raised to a given power, such as:
#    1^p + 2^p + 3^p + ... + n^p
# where p is a positive integer.

# See also: https://en.wikipedia.org/wiki/Faulhaber%27s_formula

use 5.010;
use strict;
use warnings;

use lib qw(../lib);
use Math::BigNum qw(:constant binomial);

# This function returns the nth Bernoulli number (Akiyama–Tanigawa algorithm)
# See: https://en.wikipedia.org/wiki/Bernoulli_number

examples/faulhaber_s_formula.pl  view on Meta::CPAN

        $A[$m] = 1 / ($m + 1);

        for (my $j = $m ; $j > 0 ; $j--) {
            $A[$j - 1] = $j * ($A[$j - 1] - $A[$j]);
        }
    }

    return $A[0];                    # which is Bn
}

# The Faulhaber's formula
# See: https://en.wikipedia.org/wiki/Faulhaber%27s_formula
sub faulhaber_s_formula {
    my ($p, $n) = @_;

    my $sum = 0;
    for my $j (0 .. $p) {
        $sum += binomial($p + 1, $j) * bernoulli_number($j) * ($n + 1)**($p + 1 - $j);
    }

    $sum / ($p + 1);
}

# Alternate expression using Bernoulli polynomials
# See: https://en.wikipedia.org/wiki/Faulhaber%27s_formula#Alternate_expressions
sub bernoulli_polynomials {
    my ($n, $x) = @_;

    my $sum = 0;
    for my $k (0 .. $n) {
        $sum += binomial($n, $k) * bernoulli_number($n - $k) * $x**$k;
    }

    $sum;
}

sub faulhaber_s_formula_2 {
    my ($p, $n) = @_;
    1 + (bernoulli_polynomials($p + 1, $n + 1) - bernoulli_polynomials($p + 1, 1)) / ($p + 1);
}

# Test for 1^4 + 2^4 + 3^4 + ... + 10^4
foreach my $i (0 .. 10) {
    say "$i: ", faulhaber_s_formula(4, $i)->as_rat;
    say "$i: ", faulhaber_s_formula_2(4, $i)->as_rat;
}

examples/zeta_2n.pl  view on Meta::CPAN

#!/usr/bin/perl

# Author: Daniel "Trizen" Șuteu
# License: GPLv3
# Date: 06 September 2015
# Website: https://github.com/trizen

# Calculate zeta(2n) using a closed-form formula.
# See: https://en.wikipedia.org/wiki/Riemann_zeta_function

use 5.010;
use strict;
use warnings;

use lib qw(../lib);
use Math::BigNum qw(:constant);
use constant PI => Math::BigNum->pi;

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

    bless \$r, __PACKAGE__;
}

=head2 harmfrac

    $n->harmfrac                   # => BigNum | Nan

Returns the nth-Harmonic number C<H_n>. The harmonic numbers are the sum of
reciprocals of the first C<n> natural numbers: C<1 + 1/2 + 1/3 + ... + 1/n>.

For values greater than 7000, binary splitting (Fredrik Johansson's elegant formulation) is used.

=cut

sub harmfrac {
    my ($n) = @_;

    my $ui = CORE::int(Math::GMPq::Rmpq_get_d($$n));

    $ui || return zero();
    $ui < 0 and return nan();