File Coverage

blib/lib/ConstantCalculus/CircleConstant.pm

Criterion	Covered	Total	%
statement	37	226	16.3
branch	4	56	7.1
condition	3	3	100.0
subroutine	9	28	32.1
pod	15	21	71.4
total	68	334	20.3

package ConstantCalculus::CircleConstant;

# Set the VERSION.
our $VERSION = "0.03";

# Load the Perl pragmas.
use 5.008009;
use strict;
use warnings;

# Load the Perl pragma Exporter.
use vars qw(@ISA @EXPORT @EXPORT_OK);
use Exporter 'import';

# Base class of this module.
our @ISA = qw(Exporter);

# Exporting the implemented subroutine.
our @EXPORT = qw(
    chudnovsky_algorithm
    borwein25_algorithm
    pi_borwein25
    tau_borwein25
    pi_chudnovsky
    tau_chudnovsky
    bbp_algorithm
    bbp_digits
    bbp_digit
    S
    $PI_DEC
    $PI_HEX
);

# Exporting subroutines on demand basis.
@EXPORT_OK = qw(
    factorial
    sqrtroot
    modexp
    estimate_terms
    truncate_places
);

# Disable a warning.
no warnings 'recursion';

# Load the Perl module.
use bignum;

# Set the precision for bignum.
bignum -> precision(-14);

# Initialise pre-defined values of Pi.
our $PI_DEC = undef;
our $PI_HEX = undef;

# Initialise the array with the hexadecimal numbers.
our @HX = ("0", "1", "2", "3", "4", "5", "6", "7",
           "8", "9", "A", "B", "C", "D", "E", "F");

# ---------------------------------------------------------------------------- #
# Subroutine modexp()                                                          #
#                                                                              #
# Description:                                                                 #
# The subroutine returns the result of a modular exponentiation. A modular     #
# exponentiation is an exponentiation applied over a modulus. The result of    #
# the modular exponentiation is the remainder when an integer b (the base)     #
# is exponentiated by e (the exponent) and divided by a positive integer m     #
# (modulus).                                                                   #
#                                                                              #
# @params   $b  Value of the base                                              #
#           $e  Value of the exponent                                          #
#           $m  Value of the modulus                                           #
# @returns  $y  Result of the modular exponentiation                           #
# ---------------------------------------------------------------------------- #
sub modexp {
    # Assign the subroutine arguments to the local variables.
    my($b, $e, $m) = @_;
    # Initialise the return variable.
    my $y = 1;
    # Perfom the modular exponentiation.
    do {
        ($y *= $b) %= $m if ($e % 2);
        ($b *= $b) %= $m;
    } while ($e = ($e/2) - (($e/2) % 1));
    # Return the result of the modular exponentiation.
    return $y;
};

# ---------------------------------------------------------------------------- #
# Subroutine sqrtroot()                                                        #
#                                                                              #
# Description:                                                                 #
# The subroutine calculates the square root of given number.                   #
#                                                                              #
# @params   $x  Number for which the square root is sought                     #
# @returns  $y  Square root of the given number $x                             #
# ---------------------------------------------------------------------------- #
sub sqrtroot {
    # Assign the subroutine argument to the local variable.
    my $x = $_[0];
    # Set the start values for the iteration.
    my $y = 1;
    my $y0 = 0;
    # Initialise the loop variables.
    my $i = 0;
    my $run = 1;
    # Run an infinite loop until the exit conditions is reached.
    while ($run == 1) {
        # Calculate the approximation.
        $y -= ($y * $y - $x) / (2 * $y);
        # Check the exit condition.
        $run = ($y0 == $y ? 0 : 1);
        # Save the approximation.
        $y0 = $y;
        # Increment the running variable.
        $i++;
    };
    # Return the square root.
    return $y;
};

# ---------------------------------------------------------------------------- #
# Subroutine factorial()                                                       #
#                                                                              #
# Description:                                                                 #
# The subroutine calculates the factorial of given number.                     #
#                                                                              #
# @params   $n     Number for which the factorial is sought                    #
# @returns  $fact  Factorial of the given number $n                            #
# ---------------------------------------------------------------------------- #
sub factorial {
    # Assign the subroutine argument to the local variable.
    my $n = $_[0];
    # Calculate the factorial of a number by recursion.
    my $fact = ($n == 0 ? 1 : factorial($n - 1) * $n);
    # Return the value of the factorial.
    return $fact;
};

# ---------------------------------------------------------------------------- #
# Subroutine truncate_places()                                                 #
# ---------------------------------------------------------------------------- #
sub truncate_places {
    # Assign the subroutine arguments to the local variables.
    my $decimal = $_[0];
    my $places = $_[1];
    # Calculate the multiplication factor.
    my $factor = 10**$places;
    # Truncate the places using the multiplication factor.
    $decimal = int($decimal*$factor) / ($factor);
    # Truncate the decimal using the places plus two chars.
    $decimal = substr($decimal, 0, $places + 2);
    # Return the truncated decimal number.
    return $decimal;
};

# ---------------------------------------------------------------------------- #
# Subroutine estimate_terms()                                                  #
# ---------------------------------------------------------------------------- #
sub estimate_terms {
    # Assign the subroutine arguments to the local variables.
    my $places = (defined $_[0] ? $_[0] : 0);
    my $new_places = (defined $_[1] ? $_[1] : 0);
    # Declare the return variable.
    my $number_terms = undef;
    # Determine the estimated number of terms.
    if ($places <= 1) {
        $number_terms = $places;
    } elsif ($places <= $new_places) {
        $number_terms = 1;
    } else {
        $number_terms = int($places / $new_places) + 1;
    };
    # Return the estimated number of terms.
    return $number_terms;
};

#------------------------------------------------------------------------------#
# Subroutine is_unsigned_int()                                                 #
#------------------------------------------------------------------------------#
sub is_unsigned_int {
    # Assign the subroutine argument to the local variable.
    my $arg = (defined $_[0] ? $_[0] : '');
    # Set the Perl conform regex pattern.
    my $re = qr/^(([1-9][0-9]*)|0)$/;
    # Check the argument with the regex pattern.
    my $is_unsigned_int = (($arg =~ $re) ? 1 : 0);
    # Return 0 (false) or 1 (true).
    return $is_unsigned_int;
};

# ---------------------------------------------------------------------------- #
# Subroutine check_places()                                                    #
# ---------------------------------------------------------------------------- #
sub check_places {
    # Assign the subroutine argument to the local variable.
    my $places = (defined $_[0] ? $_[0] : 0);
    # Set the string for the farewell message.
    my $msg = "The number of places is not valid. Terminating processing. Bye!";
    # Check if the given places are a valid number.
    if (is_unsigned_int($places) != 1) {
        # Print a message into the terminal window.
        print $msg . "\n";
        # Exit the script.
        exit;
    };
    # Return 1 (true).
    return 1;
};

# ---------------------------------------------------------------------------- #
# Subroutine S()                                                               #
# ---------------------------------------------------------------------------- #
sub S {
    # Assign the subroutine arguments to local variables.
    my ($j, $n, $p) = @_;
    # Set the precision for bignum.
    bignum -> precision(-$p);
    # Initialise the local variables.
    my $sum = 0;
    my $l = 0;
    my $r = 0;
    my $d = 0;
    my $k = 0;
    my $r0 = 0;
    my $run = 1;
    # Calculate the left sum.
    while ($k <= $n) {
        # Calculate the value of the denominator.
        $d = 8*$k+$j;
        # Calculate the left partial sum.
        $l = ($l + modexp(16, $n-$k, $d) / $d) % 1;
        # Increment the loop counter $k.
        $k += 1;
    };
    # Reset the loop counter.
    $k = $n + 1;
    # Calculate the right sum.
    while ($run == 1) {
        # Calculate the value of the denominator.
        $d = 8*$k+$j;
        # Calculate the right partial sum.
        $r = $r0 + (16**($n-$k)) / $d;
        # Check the exit condition of the loop.
        $run = ($r0 == $r ? 0 : 1);
        # Save the old partial sum.
        $r0 = $r;
        # Increment the loop counter $k.
        $k += 1;
    };
    # Calculate the total sum.
    $sum = $l + $r;
    # Return the total sum.
    return $sum;
};

#------------------------------------------------------------------------------#
# Subroutine bbp_digit()                                                       #
#------------------------------------------------------------------------------#
sub bbp_digit {
    # Assign the subroutine arguments to the local variables.
    my ($n, $p) = @_;
    # Set the precision for bignum.
    bignum -> precision(-$p);
    # Decrement the value with the digit of interest by 1.
    $n -= 1;
    # Calculate the fraction with the hex numbers from the nth digit upwards.
    my $dec = (4*S(1, $n, $p) - 2*S(4, $n, $p) -
               S(5, $n, $p) - S(6, $n, $p)) % 1;
    # Get the decimal representation of the digit.
    $dec = ($dec * 16) - ($dec * 16) % 1;
    # Get the hexadecimal representation of the digit.
    my $hex = $HX[$dec];
    # Return the hex number.
    return $hex;
};

#------------------------------------------------------------------------------#
# Subroutine bbp_digits()                                                      #
#------------------------------------------------------------------------------#
sub bbp_digits {
    # Assign the subroutine arguments to the local variables.
    my ($n, $digits, $p) = @_;
    # Set the precision for bignum.
    bignum -> precision(-$p);
    # Initialise the local variables.
    my $newdec = undef;
    my $hex_digit = '';
    my $hex_digits = '';
    # Decrement the value with the digit of interest by 1.
    $n -= 1;
    # Calculate the fraction with the hex numbers from the nth digit upwards.
    my $dec = (4*S(1, $n, $p) - 2*S(4, $n, $p) -
               S(5, $n, $p) - S(6, $n, $p)) % 1;
    # Calculate the hexadecimal number.
    for (my $j = 0; $j < length($dec); $j++) {
        $newdec = $dec * 16;
        $dec = $newdec - ($newdec % 1);
        $hex_digit = $HX[$dec];
        $hex_digits = $hex_digits . $hex_digit;
        $dec = $newdec - $dec;
    };
    # Return the hexadecimal number.
    return substr($hex_digits, 0, $digits);
};

#------------------------------------------------------------------------------#
# Subroutine bbp_algorith()                                                    #
#------------------------------------------------------------------------------#
sub bbp_algorithm {
    # Assign the subroutine arguments to the local variables.
    my ($places, $p) = @_;
    # Set the precision for bignum.
    bignum -> precision(-$p);
    # Initialise the local variable $pi_hex.
    my $pi_hex = "3.";
    # Declare the local variable $nth_digit.
    my $digits = undef;
    # Determine the hexadecimal places in a loop.
    for (my $i = 1; $i <= $places; $i++) {
        $digits = bbp_digit($i, $p);
        $pi_hex = $pi_hex . $digits;
    };
    # Return Pi in hexadecimal representation.
    return $pi_hex;
};

# ---------------------------------------------------------------------------- #
# Subroutine borwein25_algorithm()                                             #
# ---------------------------------------------------------------------------- #
sub borwein25_algorithm {
    # Assign the hash with the subroutine arguments to the local hash.
    my (%params) = @_;
    # Initialise the local variables.
    my $type0 = "pi";
    my $type1 = "tau";
    # Make a reference from the hash.
    my $args = \%params;
    # Initialise the locale variables.
    my $places = $args->{Places};
    my $n = $args->{Terms};
    my $p = $args->{Precision};
    my $trim = $args->{Trim};
    my $type = $args->{Type};
    # Check if the the local variables are defined.
    $places = (defined $places ? $places : 0);
    $n = (defined $n ? $n : 0);
    $p = (defined $p ? $p : $places);
    $trim = (defined $trim ? $trim : $places);
    $type = (defined $type ? $type : $type0);
    # Set the precision.
    bignum -> precision(-$p);
    # Initialise the return variable.
    my $circle_constant = 0;
    # Set the local constants.
    my $sqrt61 = sqrtroot(61);
    my $A = (212175710912 * $sqrt61) + 1657145277365;
    my $B = (13773980892672 * $sqrt61) + 107578229802750;
    my $C = (5280 * (236674 + 30303 * $sqrt61))**3;
    my $factor = sqrtroot($C) / 12;
    # Set the initial values.
    my $a_k = 0;
    my $a_sum = 0;
    my $numerator = 0;
    my $denominator = 0;
    # Set the sign to 1.
    my $sign = 1;
    # Run index is of type int.
    for (my $k = 0; $k <= $n; $k++) {
        # Calculate numerator and denominator.
        $numerator = factorial(6*$k)*($A +$k*$B);
        $denominator = (factorial($k)**3)*factorial(3*$k)*($C**($k));
        # Calculate the quotient.
        $a_k = $sign*($numerator / $denominator);
        # Add quotient to total sum.
        $a_sum += $a_k;
        # Change the sign.
        $sign *= -1;
    };
    # Calculate the value of the circle_constant.
    $circle_constant = $factor / $a_sum;
    # Determine the value for Pi or for Tau.
    if ($type eq $type0) {
        $circle_constant = 1.0 * $circle_constant;
    } elsif ($type eq $type1) {
        $circle_constant = 2.0 * $circle_constant;
    };
    # Cut the decimal number to the needed places.
    $circle_constant = truncate_places($circle_constant, $trim);
    # Return the value of the circle constant.
    return "$circle_constant";
};

# ---------------------------------------------------------------------------- #
# Subroutine chudnovsky_algorithm()                                            #
# ---------------------------------------------------------------------------- #
sub chudnovsky_algorithm {
    # Assign the hash with the subroutine arguments to the local hash.
    my (%params) = @_;
    # Initialise the local variables.
    my $type0 = "pi";
    my $type1 = "tau";
    # Make a reference from the hash.
    my $args = \%params;
    # Initialise the locale variables.
    my $places = $args->{Places};
    my $n = $args->{Terms};
    my $p = $args->{Precision};
    my $trim = $args->{Trim};
    my $type = $args->{Type};
    # Check if the the local variables are defined.
    $places = (defined $places ? $places : 0);
    $n = (defined $n ? $n : 0);
    $p = (defined $p ? $p : $places);
    $trim = (defined $trim ? $trim : $places);
    $type = (defined $type ? $type : $type0);
    # Set the precision.
    bignum -> precision(-$p);
    # Initialise the return variable.
    my $circle_constant = 0;
    # Set the local integer constants.
    my $c0 = 545140134;
    my $c1 = 13591409;
    my $c2 = 640320;
    my $c3 = 4270934400;
    my $c4 = 10005;
    # Set the initial values.
    my $a_k = 0;
    my $a_sum = 0;
    my $numerator = 0;
    my $denominator = 0;
    # Set the sign to 1.
    my $sign = 1;
    # Run index is of type int.
    for (my $k = 0; $k <= $n; $k++) {
        # Calculate numerator and denominator.
        $numerator = factorial(6*$k)*($c0*$k+$c1);
        $denominator = factorial(3*$k)*(factorial($k)**3)*($c2**(3*$k));
        # Calculate the quotient.
        $a_k = $sign*($numerator / $denominator);
        # Add quotient to total sum.
        $a_sum += $a_k;
        # Change the sign.
        $sign *= -1;
    };
    # Calculate the value of the circle_constant.
    $circle_constant = $c3 / (sqrt($c4)*$a_sum);
    # Determine Pi or Tau.
    if ($type eq $type0) {
        $circle_constant = 1.0 * $circle_constant;
    } elsif ($type eq $type1) {
        $circle_constant = 2.0 * $circle_constant;
    };
    # Cut the decimal number to the needed places.
    $circle_constant = truncate_places($circle_constant, $trim);
    # Return the value of the circle constant.
    return "$circle_constant";
};

# ---------------------------------------------------------------------------- #
# Subroutine chudnovsky_caller()                                               #
# ---------------------------------------------------------------------------- #
sub chudnovsky_caller {
    # Assign the subroutine arguments to the local variables.
    my $places = $_[0];
    my $method = $_[1];
    # Check if the places are valid. If not exit the script.
    check_places($places);
    # Calculate the required number of terms.
    my $terms = estimate_terms($places, 14);
    # Set the precision for the calculation.
    my $precision = $places + 15;
    # Get the value of the circle constant from the calculation.
    my $circle_constant = chudnovsky_algorithm(
        Places => $places,
        Terms => $terms,
        Precision => $precision,
        Trim => $places,
        Type => $method
    );
    # Return the value of the circle constant Pi or Tau.
    return "$circle_constant";
};

# ---------------------------------------------------------------------------- #
# Subroutine borwein25_caller()                                                #
# ---------------------------------------------------------------------------- #
sub borwein25_caller {
    # Assign the subroutine arguments to the local variables.
    my $places = $_[0];
    my $method = $_[1];
    # Check if the places are valid. If not exit the script.
    check_places($places);
    # Calculate the required number of terms.
    my $terms = estimate_terms($places, 25);
    # Set the precision for the calculation.
    my $precision = $places + 26;
    # Get the value of the circle constant from the calculation.
    my $circle_constant = borwein25_algorithm(
        Places => $places,
        Terms => $terms,
        Precision => $precision,
        Trim => $places,
        Type => $method
    );
    # Return the value of the circle constant Pi or Tau.
    return "$circle_constant";
};

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #
# Call the caller subroutines                                                  #
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #
sub pi_chudnovsky { return chudnovsky_caller($_[0], "pi") };
sub tau_chudnovsky { return chudnovsky_caller($_[0], "tau") };
sub pi_borwein25 { return borwein25_caller($_[0], "pi") };
sub tau_borwein25 { return borwein25_caller($_[0], "tau") };

# ============================================================================ #
# Subroutine read_pi()                                                         #
# ============================================================================ #
sub read_pi {
    # Assign the subroutine argument to the local variable.
    my $type = $_[0];
    # Initialise the local variables.
    my $line = "";
    my $pi = "";
    # Assemble the begin and end pattern.
    my $pattern0 = "# Begin Pi ${type}";
    my $pattern1 = "# End Pi ${type}";
    # Read the data from the DATA section.
    while () {
        # Read data only between the two given pattern.
        if (/$pattern0/../$pattern1/) {
            next if /$pattern0/ || /$pattern1/;
            $line = $_;
            $line =~ s/^\s+|\s+$//g;
            $pi = $pi . $line;
        };
    };
    # Return Pi.
    return $pi;
};

# ============================================================================ #
# Subroutine read_data()                                                       #
#                                                                              #
# Description:                                                                 #
# The subroutine reads Pi in decimal and hexadecimal representation from the   #
# DATA section.                                                                #
# ============================================================================ #
sub read_data {
    # Save the position of DATA.
    my $data_start = tell DATA;
    # Read decimal Pi with 1000 places.
    my $pi_dec = uc(read_pi("dec"));
    # Reposition the filehandle right past to __DATA__
    seek DATA, $data_start, 0;
    # Read hexadecimal Pi with 1000 places.
    my $pi_hex = uc(read_pi("hex"));
    # Return PI.
    return $pi_dec, $pi_hex;
};

# Read data from DATA.
($PI_DEC, $PI_HEX) = read_data();

1;

__DATA__