Re: porting C code to Perl

Replies are listed 'Best First'.
Re^2: porting C code to Perl (Inline::C) by Discipulus (Canon) on Oct 23, 2017 at 15:41 UTC
Thanks syphilis for your suggestions, 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* There are no rules, there are no thumbs.. Reinvent the wheel, then learn The Wheel; may be one day you reinvent one of THE WHEELS.	[reply] [d/l] [select]
Re^3: porting C code to Perl (Inline::C) by marioroy (Prior) on Oct 30, 2017 at 04:43 UTC
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 = 3k; mpz_fac_ui(a, 6k); // (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	[reply] [d/l] [select]
Re^4: porting C code to Perl (Inline::C) by marioroy (Prior) on Nov 06, 2017 at 02:36 UTC
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 = 3k; mpz_fac_ui(a, 6k); // (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	[reply] [d/l] [select]
Re^4: porting C code to Perl (Inline::C) by danaj (Friar) on Nov 07, 2017 at 00:48 UTC
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.	[reply]
Re^5: porting C code to Perl (Inline::C) by marioroy (Prior) on Nov 07, 2017 at 09:15 UTC
Re^4: porting C code to Perl (Inline::C) by marioroy (Prior) on Jan 23, 2018 at 09:53 UTC
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	[reply]


Think about Loose Coupling
	PerlMonks