Beefy Boxes and Bandwidth Generously Provided by pair Networks
XP is just a number
 
PerlMonks  

Re^4: porting C code to Perl (Inline::C)

by marioroy (Priest)
on Nov 06, 2017 at 02:36 UTC ( #1202793=note: print w/replies, xml ) Need Help??


in reply to Re^3: porting C code to Perl (Inline::C)
in thread porting C code to Perl

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]); }

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";

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

Regards, Mario

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://1202793]
help
Chatterbox?
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others romping around the Monastery: (2)
As of 2018-07-19 21:31 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?
    It has been suggested to rename Perl 6 in order to boost its marketing potential. Which name would you prefer?















    Results (420 votes). Check out past polls.

    Notices?