The test with Inline::C was a collateral experiment. cpan -fi Inline::C still failed with strawberry perl 5.14 I suppose while building prerequisite Win32::Mutex or it's preprerequisite Win32::IPC

Instead strawberrry Perl 5.24 64bit installed all modules smoothly. This is good.

So i retried the Inline::C test but again my ignorance won: as I understand from Inline::C docs I must have subroutine in the C code to be exported into Perl. I was not able to tranform the C code into a sub. More: the first error (or well the first that seems an error to me), complains with int len = floor(10 * N/3) + 1 pointing to the f of foor I assumed that the math.h was not imported at that time.

thanks again, but probably I have no much fu to proceed further with C, or I must take it seriously and study from basics..

L*

Hi Discipulus,

Strawberry Perl with PDL editions include Inline::C. Moreover, GMP libraries are included as well. The following is the Chudnovsky algorithm (found here) for computing Pi. It runs fine with Strawberry Perl. On Unix platforms, change the inc and libs paths to point to the GMP location.

use strict;
use warnings;

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Inline::C
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

use Inline 'C' => Config => ccflagsex => '-O2', inc => '-I/usr/local/i
+nclude', libs => '-L/usr/local/lib -lgmp -lm', clean_after_build => 0
+;

use Inline 'C' => <<'END_C';

/* http://beej.us/blog/data/pi-chudnovsky-gmp/
 *
 * Copyright (c) 2012 Brian "Beej Jorgensen" Hall <beej@beej.us>
 *
 * Permission is hereby granted, free of charge, to any person obtaini
+ng
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject 
+to
 * the following conditions:
 *
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEME
+NT.
 * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR AN
+Y
 * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT
+,
 * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */

#include <stdlib.h>
#include <string.h>
#include <gmp.h>

// how many decimal digits the algorithm generates per iteration:
#define DIGITS_PER_ITERATION 14.1816474627254776555

char *chudnovsky (unsigned long digits)
{
    mpf_t con, A, B, F, sum;
    mpz_t a, b, c, d, e;
    char *output;
    mp_exp_t exp;
    double bits_per_digit;

    unsigned long int k, threek;
    unsigned long iterations = (digits/DIGITS_PER_ITERATION)+1;
    unsigned long precision_bits;

    // roughly compute how many bits of precision we need for
    // this many digit:
    bits_per_digit = 3.32192809488736234789; // log2(10)
    precision_bits = (digits * bits_per_digit) + 1;

    mpf_set_default_prec(precision_bits);

    // allocate GMP variables
    mpf_inits(con, A, B, F, sum, NULL);
    mpz_inits(a, b, c, d, e, NULL);

    mpf_set_ui(sum, 0); // sum already zero at this point, so just FYI

    // first the constant sqrt part
    mpf_sqrt_ui(con, 10005);
    mpf_mul_ui(con, con, 426880);

    // now the fun bit
    for (k = 0; k < iterations; k++) {
        threek = 3*k;

        mpz_fac_ui(a, 6*k);  // (6k)!

        mpz_set_ui(b, 545140134); // 13591409 + 545140134k
        mpz_mul_ui(b, b, k);
        mpz_add_ui(b, b, 13591409);

        mpz_fac_ui(c, threek);  // (3k)!

        mpz_fac_ui(d, k);  // (k!)^3
        mpz_pow_ui(d, d, 3);

        mpz_ui_pow_ui(e, 640320, threek); // -640320^(3k)

        if ((threek&1) == 1) { mpz_neg(e, e); }

        // numerator (in A)
        mpz_mul(a, a, b);
        mpf_set_z(A, a);

        // denominator (in B)
        mpz_mul(c, c, d);
        mpz_mul(c, c, e);
        mpf_set_z(B, c);

        // result
        mpf_div(F, A, B);

        // add on to sum
        mpf_add(sum, sum, F);
    }

    // final calculations (solve for pi)
    mpf_ui_div(sum, 1, sum); // invert result
    mpf_mul(sum, sum, con); // multiply by constant sqrt part

    // get result base-10 in a string:
    output = mpf_get_str(NULL, &exp, 10, digits, sum); // calls malloc
+()

    // free GMP variables
    mpf_clears(con, A, B, F, sum, NULL);
    mpz_clears(a, b, c, d, e, NULL);

    return output;
}

END_C

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Perl
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

my $pi = chudnovsky(shift || 15);

print "3.", substr($pi, 1), "\n";
[download]

Regards, Mario

The Chudnovsky algorithm allows one to run parallel. C calls perl_add which adds to a shared Math::BigFloat object.

my $sum = MCE::Shared->share({ module => 'Math::BigFloat' }, 0);

sub perl_add { $sum->badd($_[0]); }
[download]

Parallel Example

use strict;
use warnings;

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Inline::C
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

use Inline 'C' => Config => ccflagsex => '-O2', inc => '-I/usr/local/i
+nclude', libs => '-L/usr/local/lib -lgmp -lm', clean_after_build => 0
+;

use Inline 'C' => <<'END_C';

/* http://beej.us/blog/data/pi-chudnovsky-gmp/
 *
 * Copyright (c) 2012 Brian "Beej Jorgensen" Hall <beej@beej.us>
 *
 * Permission is hereby granted, free of charge, to any person obtaini
+ng
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject 
+to
 * the following conditions:
 *
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEME
+NT.
 * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR AN
+Y
 * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT
+,
 * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */

/*
 * Run parallel using Perl and MCE::Shared by Mario Roy - Nov 02, 2017
+.
 */

#include <stdlib.h>
#include <gmp.h>

// how many decimal digits the algorithm generates per iteration:
#define DIGITS_PER_ITERATION 14.1816474627254776555

unsigned long c_init (unsigned long digits)
{
    unsigned long iterations = (digits/DIGITS_PER_ITERATION) + 1;
    unsigned long precision_bits;
    double bits_per_digit;

    // roughly compute how many bits of precision we need for
    // this many digit:
    bits_per_digit = 3.32192809488736234789; // log2(10)
    precision_bits = (digits * bits_per_digit) + 1;

    mpf_set_default_prec(precision_bits);

    return iterations;
}

void c_loop (unsigned long digits, unsigned long beg, unsigned long en
+d)
{
    unsigned long int k, threek;

    mpz_t a, b, c, d, e;
    mpf_t A, B, F, sum;

    // allocate GMP variables
    mpz_inits(a, b, c, d, e, NULL);
    mpf_inits(A, B, F, sum, NULL);

    // now the fun bit
    for (k = beg; k <= end; k++) {
        threek = 3*k;

        mpz_fac_ui(a, 6*k);  // (6k)!

        mpz_set_ui(b, 545140134); // 13591409 + 545140134k
        mpz_mul_ui(b, b, k);
        mpz_add_ui(b, b, 13591409);

        mpz_fac_ui(c, threek);  // (3k)!

        mpz_fac_ui(d, k);  // (k!)^3
        mpz_pow_ui(d, d, 3);

        mpz_ui_pow_ui(e, 640320, threek); // -640320^(3k)

        if ((threek&1) == 1) { mpz_neg(e, e); }

        // numerator (in A)
        mpz_mul(a, a, b);
        mpf_set_z(A, a);

        // denominator (in B)
        mpz_mul(c, c, d);
        mpz_mul(c, c, e);
        mpf_set_z(B, c);

        // result
        mpf_div(F, A, B);

        // add on to sum
        mpf_add(sum, sum, F);
    }

    // free GMP variables
    mpz_clears(a, b, c, d, e, NULL);
    mpf_clears(A, B, F, NULL);

    char *output = (char *) malloc(digits + 30);
    gmp_sprintf(output, "%0.*Fe", digits, sum);
    mpf_clears(sum, NULL);

    char *args[] = { output, NULL };
    call_argv("perl_add", G_DISCARD, args);

    free(output);
}

char* c_final (unsigned long digits, char *str)
{
    mpf_t sum;
    mpf_inits(sum, NULL);
    mpf_set_str(sum, str, 10);

    char *output;
    mp_exp_t exp;

    // first the constant sqrt part
    mpf_t con;
    mpf_inits(con, NULL);
    mpf_sqrt_ui(con, 10005);
    mpf_mul_ui(con, con, 426880);

    // final calculations (solve for pi)
    mpf_ui_div(sum, 1, sum); // invert result
    mpf_mul(sum, sum, con); // multiply by constant sqrt part

    // get result base-10 in a string:
    output = mpf_get_str(NULL, &exp, 10, digits, sum); // calls malloc
+()

    mpf_clears(con, sum, NULL);

    return output;
}

END_C

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Perl
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

use MCE::Hobo;
use MCE::Shared;

use Math::BigFloat lib => 'GMP';
use Signal::Safety;

my $digits     = shift || 15;
my $iterations = c_init($digits);
my $max_hobos  = MCE::Util::get_ncpu();
my $chunk_size = int($iterations / $max_hobos / 16 + 1);

my $sum = MCE::Shared->share({ module => 'Math::BigFloat' }, 0);

my $seq = MCE::Shared->sequence(
    { chunk_size => $chunk_size, bounds_only => 1 },
    0, $iterations - 1
);

sub perl_add {
    $sum->badd($_[0]);
}

{
    local $Signal::Safety = 0;

    my $task = sub {
        while ( my ($beg, $end) = $seq->next() ) {
            c_loop($digits, $beg, $end);
        }
    };

    MCE::Hobo->create($task) for 1 .. $max_hobos;
    MCE::Hobo->wait_all();
}

my $pi = c_final($digits, $sum->bstr());

print "3.", substr($pi, 1), "\n";
[download]

I'm using Signal::Safety so that pressing Ctrl-C responds immediately while inside C.

Regards, Mario

This seems to be much slower than AGM. 10k digits takes 0.13s but 100k digits takes 33.8s. Even the Machin-type formula are faster.

The paper Exploration in π Calculation Using Various Methods Implemented in Python has some insight, and is useful to read in general as well. We certainly shouldn't be doing all those factorials in the loop. Those are pretty easy to pull out. But the killer is the huge-precision divide done for every loop iteration. Binary splitting is what's needed, though it is a non-trivial change. You can see in their Table 2 just how significant the difference is. If I recall correctly, fast implementations use fine control of the precision on every loop iteration, where in contrast AGM has to use full precision for each iteration.

I got caught up with computing Pi. Chudnovsky-Pi lives here on Github. It runs well on Cygwin and Unix, of course ;-). Strawberry Perl 64-bit is not built with true 64-bit integer support, unfortunately. Anyway, the Chudnovsky-Pi demonstration runs with Linux Bash Shell on Windows 10. The m4 binary is needed and optionally yasm, needed for MPIR. Hence, sudo apt-get install m4 yasm.

Cheers, Mario