view release on metacpan or  search on metacpan

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

        my $p = Math::MPFR::Rmpfr_init2($prec);
        Math::MPFR::Rmpfr_const_pi($p, $ROUND);                        # p = PI
        Math::MPFR::Rmpfr_pow_ui($p, $p, $n, $ROUND);                  # p = p^n
        Math::MPFR::Rmpfr_div($f, $f, $p, $ROUND);                     # f = f/p

        Math::GMPz::Rmpz_set_ui($z, 1);                                # z = 1
        Math::GMPz::Rmpz_mul_2exp($z, $z, $n + 1);                     # z = 2^(n+1)
        Math::GMPz::Rmpz_sub_ui($z, $z, 2);                            # z = z-2

        Math::MPFR::Rmpfr_mul_z($f, $f, $z, $ROUND);                   # f = f*z
        Math::MPFR::Rmpfr_round($f, $f);                               # f = [f]

        my $q = Math::GMPq::Rmpq_init();
        Math::MPFR::Rmpfr_get_q($q, $f);                               # q = f
        Math::GMPq::Rmpq_set_den($q, $z);                              # q = q/z
        Math::GMPq::Rmpq_canonicalize($q);                             # remove common factors

        Math::GMPq::Rmpq_neg($q, $q) if $n % 4 == 0;                   # q = -q    (iff 4|n)
        return bless \$q, __PACKAGE__;
    }

#<<<
    my @D = (
        Math::GMPz::Rmpz_init_set_ui(0),
        Math::GMPz::Rmpz_init_set_ui(1),
        map { Math::GMPz::Rmpz_init_set_ui(0) } (1 .. $n/2 - 1)
    );
#>>>

    my ($h, $w) = (1, 1);
    foreach my $i (0 .. $n - 1) {
        if ($w ^= 1) {
            Math::GMPz::Rmpz_add($D[$_], $D[$_], $D[$_ - 1]) for (1 .. $h - 1);
        }
        else {
            $w = $h++;
            Math::GMPz::Rmpz_add($D[$w], $D[$w], $D[$w + 1]) while --$w;
        }
    }

    my $den = Math::GMPz::Rmpz_init_set($ONE_Z);
    Math::GMPz::Rmpz_mul_2exp($den, $den, $n + 1);
    Math::GMPz::Rmpz_sub_ui($den, $den, 2);
    Math::GMPz::Rmpz_neg($den, $den) if $n % 4 == 0;

    my $r = Math::GMPq::Rmpq_init();
    Math::GMPq::Rmpq_set_num($r, $D[$h - 1]);
    Math::GMPq::Rmpq_set_den($r, $den);
    Math::GMPq::Rmpq_canonicalize($r);

    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();

    # Use binary splitting for large values of n. (by Fredrik Johansson)
    # http://fredrik-j.blogspot.ro/2009/02/how-not-to-compute-harmonic-numbers.html
    if ($ui > 7000) {

        my $num = Math::GMPz::Rmpz_init_set_ui(1);

        my $den = Math::GMPz::Rmpz_init();
        Math::GMPz::Rmpz_set_q($den, $$n);
        Math::GMPz::Rmpz_add_ui($den, $den, 1);

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

        # Inspired by Dana Jacobsen's code from Math::Prime::Util::{PP,GMP}.
        #   https://metacpan.org/pod/Math::Prime::Util::PP
        #   https://metacpan.org/pod/Math::Prime::Util::GMP
        my $sub;
        $sub = sub {
            my ($num, $den) = @_;
            Math::GMPz::Rmpz_sub($temp, $den, $num);

            if (Math::GMPz::Rmpz_cmp_ui($temp, 1) == 0) {
                Math::GMPz::Rmpz_set($den, $num);
                Math::GMPz::Rmpz_set_ui($num, 1);
            }
            elsif (Math::GMPz::Rmpz_cmp_ui($temp, 2) == 0) {
                Math::GMPz::Rmpz_set($den, $num);
                Math::GMPz::Rmpz_mul_2exp($num, $num, 1);
                Math::GMPz::Rmpz_add_ui($num, $num, 1);
                Math::GMPz::Rmpz_addmul($den, $den, $den);
            }
            else {
                Math::GMPz::Rmpz_add($temp, $num, $den);
                Math::GMPz::Rmpz_tdiv_q_2exp($temp, $temp, 1);
                my $q = Math::GMPz::Rmpz_init_set($temp);
                my $r = Math::GMPz::Rmpz_init_set($temp);
                $sub->($num, $q);
                $sub->($r,   $den);
                Math::GMPz::Rmpz_mul($num,  $num, $den);
                Math::GMPz::Rmpz_mul($temp, $q,   $r);
                Math::GMPz::Rmpz_add($num, $num, $temp);
                Math::GMPz::Rmpz_mul($den, $den, $q);
            }
        };

        $sub->($num, $den);

        my $q = Math::GMPq::Rmpq_init();
        Math::GMPq::Rmpq_set_num($q, $num);
        Math::GMPq::Rmpq_set_den($q, $den);