Beefy Boxes and Bandwidth Generously Provided by pair Networks
Problems? Is your data what you think it is?
 
PerlMonks  

porting C code to Perl

by Discipulus (Canon)
on Oct 22, 2017 at 19:27 UTC ( [id://1201848]=perlquestion: print w/replies, xml ) Need Help??

Discipulus has asked for the wisdom of the Perl Monks concerning the following question:

Hello nuns and monks,

if you know me or not, I'm completely unaware of other programming languages; I just know a little Perl but I found myself in the rare situation where I need to translate a little code from C to Perl.

Infact I rapidly (*) understood that is a task to be done by hand, and I've been told many times that these two languages share a lot in their syntax.

My attempt is below and does not produces the output I expected (**).

I looked a bit to some description of the C syntax to try to understand if there were some difference between the C operator and the correspective Perl's one. For example for ++ autoincrement or arithmentics ones. I found nothing relevant: many operators seems to act the same. Doubts remain about the C array syntax ( int A[len] ?? that I read as the elelment len of the array A is an int but..).

Here below my attempt: can someone be so kind to point me where I lost in the translation? After a plain translation, when I possibly end with some working Perl code I'll arrange it into a more perlish version.

use strict; use warnings; # #include <math.h> # #include <stdio.h> # #define N 100 my $n=100; # int len = floor(10 * N/3) + 1; my $len = 1 + int (10 * $n / 3); # int A[len]; my @a; $#a=$len-1; #? -1 ???? # for(int i = 0; i < len; ++i) {A[i] = 2;} for (my $i = 0; $i < $len; $i++){ $a[$i]= 2; } # int nines = 0; my $nines = 0; # int predigit = 0; my $predigit = 0; # for(int j = 1; j < N + 1; ++j) { for (my $j = 1; $j < $n + 1; ++$j){ # int q = 0; my $q = 0; # for(int i = len; i > 0; --i) { for(my $i = $len; $i > 0; $i--){ # int x = 10 * A[i-1] + q*i; my $x = 10 * $a[$i-1] + $q * $i; # A[i-1] = x % (2*i - 1); $a[$i-1] = $x % (2 * $i - 1); # q = x / (2*i - 1); } $q = $x / (2 * $i - 1); } # A[0] = q%10; $a[0]=$q%10; # q = q/10; $q=$q/10; # if (9 == q) { ++nines;} if (9 == $q){ ++$nines; } # else if (10 == q) { elsif(10 == $q){ # printf("%d", predigit + 1); printf("%d", $predigit + 1); # for (int k = 0; k < nines; ++k) { printf("%d", 0); } for (my $k = 0; $k < $nines; $k++){ printf("%d", 0);} # predigit, nines = 0; $predigit = $nines = 0; # } } # else { else{ # printf("%d", predigit); printf("%d", $predigit); # predigit = q; $predigit = $q; # if (0 != nines) { if(0 != $nines){ # for (int k = 0; k < nines; ++k) { for (my $k = 0; $k < $nines; $k++) { # printf("%d", 9); printf("%d", 9); # } } # nines = 0; $nines = 0; } } } # printf("%d", predigit); printf("%d", $predigit);

(*) Not so rapid: I found a thread here at PM with a link to a C to perl translator, but I missed the <ironic> tags and my hands called rapidly gcc -P -E file.c > file.pl using the compiler I have shipped within strawberry perl. Was not useful and the next line gcc -Larry -Wall file.c > file.pl revealed me it was an humoristic faq or well a iaq..

(**) Well i wanted to verify that the C code printed what I expected and I tried blindly using gcc to compile it (??) and using Inline::C but I had no success to not even install it on my strawberry perl.

If you need my C source you can, obviously perl -lne 'print if s/\s?#\s//' my_post.pl

Thanks

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.

Replies are listed 'Best First'.
Re: porting C code to Perl (updated!)
by haukex (Archbishop) on Oct 22, 2017 at 19:57 UTC

    All of your C variables are ints, so you need to change the two divisions in your program like so:

    $q = int( $x / (2 * $i - 1) ); ... $q = int( $q / 10 );

    And then the two programs produce the same output (Update: "03141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067"*, in case anyone was wondering). Alternatively, you can stick a use integer; at the top of the Perl file (of course this only applies when all the variables in the program are integers).

    Also, although this doesn't seem to affect your program, note that there are differences between C's modulo % and Perl's % (use integer; also affects this).

    * Update 2: The C compiler (gcc -std=c99 -Wall -Wextra -Wpedantic 1201848.c) complained about the line predigit, nines = 0;: "warning: left-hand operand of comma expression has no effect", so initially I adjusted it to your interpretation, predigit=0; nines=0;, but actually that causes a slight variation in the output (first one doesn't change predigit, second one does set it to zero):

    0314159265358979323846264338327954288419716939937510582097494459230781 +6406286208998628734825342117067 0314159265358979323846264338327950288419716939937510582097494459230781 +6406286208998628034825342117067 here ^ + and here ^

    In either case, making the equivalent change to the Perl program produces the same output. Since these are the digits of π, obviously the second one is the correct one. BTW, I guess you took the C code from here?

    Update 3: Minor ninja edits to the above. I'm done for now :-)

      thanks a lot haukex!

      your insights are exactly what I needed. I spotted the type casting of int but I have not thought it constrained all the operations too.. what a complicated world they have with C..

      Now I ported to a more perlish versions, but speaking frankly, I'm quite struggling with it's efficiency. Infact the inner for loop for a precision of 1000 cause each of it's line to be called something like 3 million of times if I read correctly the Devel::NYTProf output.

      UPDATE some output from NYTProf Performance Profile

      Line Statements Time on line Code 9 1002 49.4ms for $i(reverse 1..$l){ 10 3347682 441ms $x=10*$a[$i-1]+$q*$i; 11 3347682 474ms $a[$i-1]=$x%(2*$i-1); 12 3347682 882ms $q=int($x/(2*$i-1))

      This inner for loop if I read correctly it's not influenced by the outer $j and only set a final (?) $q and the array element $a[$i-1] so i'm trying to cut the loop from there and prebuild an array of q linkable with $j indexes and prefilling the @a accordingly. No big luck until now.

      In reality I'm working blindly not having well understood what the code do. I'll try a reverse eng. based on print statements if i do not desist sooner. Another ugly think I've done is the last line where I merely susbtitute the leading 03 with 3, any attempt to do this within the loop failed.

      Anyway thank you again!

      use strict; use warnings; my $n=1002; my $len = 1 + int (10 * $n / 3); my $pi; my @a= (2) x $len; my $nines = 0; my $predigit = 0; foreach my $j(1..$n){ my $q = 0; foreach my $i( reverse 1..$len){ my $x = 10 * $a[$i-1] + $q * $i; $a[$i-1] = $x % (2 * $i - 1); $q = int($x / (2 * $i - 1)); } $a[0]=$q%10; $q=int($q/10); if (9 == $q){ ++$nines; } elsif(10 == $q){ $pi.=$predigit + 1; for (my $k = 0; $k < $nines; $k++){$pi.=0} $predigit = $nines = 0; } else{ $pi.= int $predigit; $predigit = $q; if(0 != $nines){ for (my $k = 0; $k < $nines; $k++) { $pi.=9; } $nines = 0; } } } print $pi=~s/^03/3./r;

      PS and yes these are the Pi digits and the C source is what you spotted!

      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.
        I spotted the type casting of int but I have not thought it constrained all the operations too

        In C, when you say int q;, that variable is fixed to, depending on how your compiler defines it, a 32 (or 16 or 64) bit signed integer value, no more, no less - simplifying a bit, no matter what you try to assign to it, it will be cast back to an int.

        Now I ported to a more perlish versions, but speaking frankly, I'm quite struggling with it's efficiency.

        Since functionally the two programs are practically identical, that'll simply be the difference between a complied and interpreted language. I don't have much time for optimization right now, but you might try to compare the performance of the program with and without use integer; - I'm not sure if it'll make a significant difference but it's worth a try.

        In reality I'm working blindly not have well understood wht the code do.

        A bit of googling on the algorithm brings me to this page, while it's a bit difficult to read visually, it seems to give an example-based explanation of the "spigot algorithm" that this appears to be an implementation of.

        Update: Implemented the divisor suggestion by soonix.

        Hi Discipulus,

        To have Perl match the C output (regarding the number of characters), a missing append is needed at the end of the Perl script. The C code is found here and at stackoverflow. Regarding C, a fix is necessary described here (under Update 2) and here.

        $pi .= $predigit; # <-- missing line print $pi =~ s/^03/3./r;

        I validated using the following code, based on yours and replies by fellow monks. Notice the Perlish for my $k (0..$nines - 1) { ... } near the end of the script. In your demonstration, calling int on $predigit isn't necessary because $q contains an integer value; e.g. $q = int( ... ).

        # The spigot algorithm for calculating the digits of Pi. # http://www.cut-the-knot.org/Curriculum/Algorithms/SpigotForPi.shtml use strict; use warnings; my $n = 100; my $len = int(10 * $n / 3) + 1; my $pi; my @a = (2) x $len; my $nines = 0; my $predigit = 0; for my $j (1..$n) { my $q = 0; for my $i (reverse 0..$len - 1) { my $x = 10 * $a[$i] + $q * ($i + 1); my $divisor = 2 * $i + 1; $a[$i] = $x % $divisor; $q = int($x / $divisor); } $a[0] = $q % 10; $q = int($q / 10); if (9 == $q) { ++$nines; } elsif (10 == $q) { $pi .= $predigit + 1; for my $k (0..$nines - 1) { $pi .= 0; } $predigit = $nines = 0; } else { $pi .= $predigit; $predigit = $q; if (0 != $nines) { for my $k (0..$nines - 1) { $pi .= 9; } $nines = 0; } } } $pi .= $predigit; print $pi, "\n";

        Pi computed to 10,000 digits is found here.

        Update: The above runs 1.4 times faster with small changes. Basically, use integer. Thanks, haukex.

        2a3 > use integer; 5c6 < my $len = int(10 * $n / 3) + 1; --- > my $len = 10 * $n / 3 + 1; 20c21 < $q = int($x / $divisor); --- > $q = $x / $divisor; 24c25 < $q = int($q / 10); --- > $q = $q / 10;

        Regards, Mario

        foreach my $i( reverse 1..$len){ my $x = 10 * $a[$i-1] + $q * $i; $a[$i-1] = $x % (2 * $i - 1); $q = int($x / (2 * $i - 1)); }
        not sure if it's better or worse, but I'd change that to
        foreach my $i( reverse 0 .. $len-1){ my $x = 10 * $a[$i] + $q * ($i+1); my $divisor = 2 * $i + 1; $a[$i] = $x % $divisor; $q = int($x / $divisor); }
        not for efficiency, but for (perhaps) easier understanding
        • $i one less, so the array index is a simple variable
        • the common expression of the last two lines as a variable of its own
        but I have not thought it constrained all the operations too.. what a complicated world they have with C

        In Perl a variable is a concept - you can put whatever you want into it and Perl will try to do something suitable. In C a variable is a contract - you (and the compiler) must honor it and that includes forced conversions if necessary (and allowed). In a 'do what I mean' sense C is more complicated, in a 'what states do I need to consider' sense it is less complicated.

        If I write a library function in C then I do not have to worry if my input variables may contain an integer, a float, an array, a hash, a reference, an object or whatever. Each variable can contain exactly one of these (unless I specifically do not want to care about it). This may be a boon or a burden. It depends.


        (Actually you do not have to honor the contract. You can throw a tantrum, stomp with your feet and break it. But you should.)

Re: porting C code to Perl
by syphilis (Archbishop) on Oct 23, 2017 at 11:15 UTC
    .... using Inline::C but I had no success to not even install it on my strawberry perl

    On my own Windows builds of perl, I often find that Inline::C's t/27inline_maker.t test script fails some tests.
    I haven't bothered investigating and, since it's not a failure that interferes with anything that I want to do, I just force install it (cpan -fi Inline::C).
    Specifically, I get a test report that shows:
    Test Summary Report ------------------- t/27inline_maker.t (Wstat: 1024 Tests: 8 Failed: 4) Failed tests: 2-3, 6-7 Non-zero exit status: 4
    However, all tests usually pass when I install Inline::C on Strawberry Perl. What's the problem that you're seeing there ?

    Cheers,
    Rob
      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.

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

        Regards, Mario

Re: porting C code to Perl
by Mr. Muskrat (Canon) on Oct 23, 2017 at 16:53 UTC

    haukex and syphilis have already provided some good answers but I thought that I would point out some things that I didn't see mentioned.

    Perl's int is not the same as C's floor(). What you want is POSIX::floor. It has the identical behavior as C's floor().

    Perl has a pre-decrement just like C.

    # for(int i = len; i > 0; --i) { # C # for(my $i = $len; $i > 0; $i--){ # incorrect Perl for(my $i = $len; $i > 0; --$i){ # correct Perl

    YMMV as I haven't tried it out with these changes implemented.

      Hello Mr. Muskrat and thanks for your reply,

      > Perl's int is not the same as C's floor()..

      I already admitted my superb ignorance but i read here that in C floor turns 1.6 1.2 2.8 2.3 into 1 1 2 2 or simply 'This function returns the largest integral value not greater than x'.

      This is identical to the int function where any decimal are stripped out: say int $_ for qw(1.6 1.2 2.8 2.3)

      > Perl has a pre-decrement just like C..

      I have meditated a bit over this before trying to port the code and, sincerely I didnt read the C documents about post and pre increment, but as I have always understood the third part of the Perl's C style for loop is a statement on it's own: ie ++$var; or $var++; are identical.

      Infact for($i=0;$i<4;$i++){say $i} and for($i=0;$i<4;++$i){say $i} produce the very same, expectable results. Switching between them in the above code does not change the result.

      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.
        the third part of the Perl's C style for loop is a statement on it's own: ie ++$var; or $var++; are identical.

        Yes, that's correct*, when they're written standalone as in for(;;$var++). ++$var vs. $var++ only makes a difference when its return value is used: ++$var will first increment the variable and then return the new value, while $var++ will return the old value and then increment the variable.

        $ perl -le 'for ( my $x=0; $x<3; print $x++ ) { }' 0 1 2 $ perl -le 'for ( my $x=0; $x<3; print ++$x ) { }' 1 2 3

        * Update: In fact, note what Perl does here (edited for brevity):

        $ perl -MO=Deparse -e 'for ( my $x=0; $x<3; print $x++ ) { }' for (my $x = 0; $x < 3; print $x++) { (); } $ perl -MO=Deparse -e 'for ( my $x=0; $x<3; print ++$x ) { }' for (my $x = 0; $x < 3; print ++$x) { (); } $ perl -MO=Deparse -e 'for ( my $x=0; $x<3; ++$x ) { }' for (my $x = 0; $x < 3; ++$x) { (); } $ perl -MO=Deparse -e 'for ( my $x=0; $x<3; $x++ ) { }' for (my $x = 0; $x < 3; ++$x) { (); }
        For this particular use case, it doesn't matter, but with negative numbers, int and floor give different results (this is mentioned explicitly in the first link)
Re: porting C code to Perl
by danaj (Friar) on Oct 25, 2017 at 02:02 UTC

    I should make a blog post about this, but here are some thoughts.

    The C code you started with is really poor quality, seemingly because it got involved with code golfing. Global variables, sigh -- that was ugly in the 1980s and still is. This makes it tough to follow and translate.

    The rest of this is really all about performance. If you're happy using it, then great, the rest of this is just blathering on about details. But ever since as a kid I bought Petr Beckmann's book I've found this interesting. It was also one of my many digressions -- I needed a few million digits of Pi for a puzzle and was disappointed that the standard Perl modules were horribly slow, so went down the rat hole.

    The spigot code you're using is one of the slower methods for generating N digits of Pi. There is a modified version of the simple C spigot from Winter, Flammenkamp, et al.) that runs about 8 times faster (code in ntheory, Luschny's Java version, many other variants on the web). I added a little bit to make sure it rounds correctly at reasonable sizes (the table maker's dilemma still applies, but not for quite some time). I actually only use it for <= 100 digits or if MPFR and GMP are not available, as those methods are far faster.

    You can see a timing comparison of some C methods I measured in 2014 on RosettaCode. Spigot 1 in that table is Winter/Flammenkamp version (as mentioned, about 8x faster than the one you're using). Here is a (hastily assembled) chart of 13 algorithm/implementations I made today, all but one used from Perl. The time (in seconds) on the y-axis uses a log scale which lets us see how AGM and Chudnovsky have a favorably slope compared to the others (about 2.4x time per doubling of digits vs. 4x).

    Pi timing graph

    • y-cruncher is what you want if you really need a few billion digits for some project. It's pretty cool if you're interested in this stuff.
    • Pari's Chudnovsky method (they call it Ramanujan + binary splitting) is definitely faster than AGM. This is basically the same method as Xue's Chudnovsky program shown on the GMP site, though implemented differently. That's the one method that was called using GP rather than plugged into Perl. Putting a GMP version of this in my module has been on my to-do list for a couple years.
    • My AGM implementation in GMP (used in Math::Prime::Util aka ntheory when possible) is just a fraction slower than MPFR (I use the normal GMP mpf interface, they use MPFR of course). It's short code, quite fast, and asymptotically looks nice.
    • The Machin formulas implemented in GMP are inferior to AGM or Chudnovsky at all sizes and get worse. Using updated formulas helps only slightly. They still easily beat the spigot methods and all but one of the Perl methods. Math::BigFloat used to use these (but in Perl not C+GMP).
    • The spigot algorithms don't compete well for this task (producing N correctly rounded digits of Pi). There is no point in using GMP for these particular ones. Implementing it in Perl directly has a very large performance hit, as all these operations are really simple so C is going to shine.
    • Math::BigFloat changed from Machin to AGM in version 1.999705 (2015). rt://90033 has some performance numbers. With the GMP back end it isn't bad. The overhead at small sizes isn't ideal. With the Pari or Calc back ends it's terribly slow (assuming you want 1000+ digits).
    • I suspect if you used Math::GMPf or Math::AnyNum and wrote the simple AGM code it would be only a small penalty from the GMP code, so easily beat all the other Perl code shown. At that point you could just use the pi() function from Math::AnyNum which calls MPFR, or Pi() from ntheory which uses its best available choice from MPFR, GMP, C, or pure Perl.

      Update: Running ./configure is all one needs. GMP builds fine on the Mac and sets MPN_PATH properly. The custom configure-gpm.sh script isn't necessary.

      Update: Ditto for MPFR. Running ./configure --with-gmp=/usr/local is all one needs. Please forget the custom configure-mpfr.sh script.

      Update: It requires sudo to install the libraries under /usr/local.

      Hello danaj and like Discipulus am astonished too.

      This post describes how I build the GMP and MPFR libraries on my Mac (Haswell architecture), beneficial for ntheory. I have found that setting MPN_PATH is helpful. The -march=core2 -mtune=intel compiler options enable extra performance as well, particularly -mtune=intel.

      On my CentOS 7.3 VM, I just run sudo yum install gmp-devel mpfr-devel before installing Math::Prime::Util.

      GMP

      Below, "configure-gmp.sh" to configure the GMP library. Remove the architecture(s) not supported by your processor (MPN_64). Otherwise, make check will fail.

      PREFIX="/usr/local" OPT_FLAGS="-O3 -march=core2 -mtune=intel -fomit-frame-pointer" OS_ARCH="`uname -m`" OS_REV="`uname -r | cut -d. -f1`" TRIPLE="${OS_ARCH}-apple-darwin${OS_REV}" export CXXFLAGS="${OPT_FLAGS} -pipe -Wall -fno-common -DPIC" export CFLAGS="${OPT_FLAGS} -pipe -Wall -fno-common -DPIC" export LDFLAGS="-Wl,-no_pie" export CXX=g++ export CC=gcc # MPN_64="x86_64/coreibwl x86_64/coreihwl x86_64/coreisbr x86_64/corei +nhm x86_64/fastsse x86_64/core2 x86_64 generic" MPN_64="x86_64/coreihwl x86_64/coreisbr x86_64/coreinhm x86_64/fastsse + x86_64/core2 x86_64 generic" # Replace two options with --disable-static --enable-shared for shared + libraries MPN_PATH="${MPN_64}" \ ./configure \ --build="${TRIPLE}" --prefix="${PREFIX}" --disable-shared --enable- +static \ --enable-cxx --enable-assembly --enable-fft

      Build instructions.

      cd gmp-6.1.2 sh /path/to/configure-gmp.sh make -j 4 make check # <-- important, do not skip sudo make install-strip

      MPFR

      Below, "configure-mpfr.sh" to configure the MPFR library.

      PREFIX="/usr/local" OPT_FLAGS="-O3 -march=core2 -mtune=intel -fomit-frame-pointer" OS_ARCH="`uname -m`" OS_REV="`uname -r | cut -d. -f1`" TRIPLE="${OS_ARCH}-apple-darwin${OS_REV}" export CXXFLAGS="${OPT_FLAGS} -pipe -Wall -fno-common -DPIC" export CFLAGS="${OPT_FLAGS} -pipe -Wall -fno-common -DPIC" export LDFLAGS="-Wl,-no_pie" export CXX=g++ export CC=gcc # Replace two options with --disable-static --enable-shared for shared + libraries ./configure \ --build="${TRIPLE}" --prefix="${PREFIX}" --disable-shared --enable- +static \ --disable-dependency-tracking --disable-assert --enable-thread-safe + \ --with-gmp="${PREFIX}"

      Build instructions. Download mpfr-3.1.2 and allpatches (important).

      cd mpfr-3.1.2 patch -N -Z -p1 < /path/to/allpatches.txt sh /path/to/configure-mpfr.sh make -j 4 make check # <-- important, do not skip sudo make install-strip

      Math::Prime::Util

      Now that the two libraries are installed, (re)install Math::Prime::Util using your favorite method. It will pick dependencies such as Math::Prime::Util::GMP and Math::BigInt::GMP. Without GMP/MPFR, it takes 14 seconds to compute Pi to 50,000 digits. With GMP/MPFR, it takes 1.5 seconds to compute Pi to 1 million digits.

      $ time /opt/perl-5.26.1/bin/perl -Mntheory=:all -E 'say Pi(1e6)' | wc +-c 1000002 real 0m1.428s user 0m1.408s sys 0m0.017s

      No matter the Unix OS, ensure the GMP and MPFR development libraries are installed before installing Math::Prime::Util. Most OS vendors have packages already made.

      Regards, Mario

        GMP is installed on many systems now, partly because GCC started requiring it. MPFR is much less common. GMP is a big help for lots of things, but MPFR much less so, and I used it via Math::MPFR so you'd need that installed also. It looks like it is used in prime_count_lower and prime_count_upper for large enough values, and the floating point routines harmreal, ExponentialIntegral, LogarithmicIntegral, RiemannZeta, RiemannR, and Pi. For Pi it's a very small savings vs.the standard GMP back end. But the floating point things like RiemannZeta are just hard to do in Perl or GMP's rudimentary mpf layer (for me at least). Using Math::MPFR (thanks syphilis!) is thousands of times faster than Math::BigFloat even using the fast GMP backend. The MPFR implementors are very sharp and also better practicing numerical analysts than me.

        For AWS EC2 machines I usually just install the packages with apt-get or yum. I also do git and clone/build/install the latest versions of my modules but that's because I want the absolute latest. cpan is fine and as long as you have GMP you should just be able to cpan Math::Prime::Util Math::BigInt::GMP. But make sure it did find GMP (try cpan Math::Prime::Util::GMP to check).

        14 seconds for Pi to 50k digits... Ouch. Moar fastar! At least it's faster than the default Math::BigFloat (which doesn't have the luxury of dropping into custom C code).

        As Discipulus says I should try use integer as well to test the Perl spigot. I usually don't like that pragma because it's signed which is never what I want, but that's not relevent here. Now I'm pondering a digression into how to implement use uinteger. At first glance it doesn't look terrible hard to do in the core, but using a whole new hint bit is unlikely to fly. Maybe a test just to see if it would work...

        Update: Added before/after results for Pi(1e9) to emphasize the performance increase.

        I did an experiment on an i7-Broadwell machine running CentOS 7.3. The system gmp-devel libraries were built to support the least common denominator. Thus, leaving out optimizations for newer CPU instructions.

        # default package build MPN_PATH=" x86_64/k8 x86_64 generic" # manually on the broadwell machine, includes extra optimizations MPN_PATH=" x86_64/coreibwl x86_64/coreihwl x86_64/coreisbr x86_64/core +inhm x86_64/core2 x86_64 generic"

        I have everything point to the /opt/perl-5.26.1 folder to not impact the system in any way. Then, re-installed GMP modules including Math::Prime::Util.

        cd gmp-6.1.2 ./configure --prefix=/opt/perl-5.26.1 make -j 4 make check make install cd Math-GMP-2.15 /opt/perl-5.26.1/bin/perl Makefile.PL INC="-I/opt/perl-5.26.1/include" + LIBS="-L/opt/perl-5.26.1/lib -lgmp" make make test make install cd Math-BigInt-GMP-1.6004 /opt/perl-5.26.1/bin/perl Makefile.PL INC="-I/opt/perl-5.26.1/include" + LIBS="-L/opt/perl-5.26.1/lib -lgmp" make make test make install cd Math-Prime-Util-GMP-0.46 /opt/perl-5.26.1/bin/perl Makefile.PL INC="-I/opt/perl-5.26.1/include" + LIBS="-L/opt/perl-5.26.1/lib -lgmp" make make test make install cd Math-Prime-Util-master /opt/perl-5.26.1/bin/perl Makefile.PL INC="-I/opt/perl-5.26.1/include" + LIBS="-L/opt/perl-5.26.1/lib -lgmp" make make test make install

        Beyond 20% gain is possible by simply building GMP manually. That will set MPN_PATH to utilize the full capabilities of the processor.

        time /opt/perl-5.26.1/bin/perl -Mntheory=:all -E 'say Pi(1e6)' | wc -c 1000002 0m 1.524s using system gmp libraries 0m 1.168s using gmp libraries from /opt/perl-5.26.1/lib time /opt/perl-5.26.1/bin/perl -Mntheory=:all -E 'say Pi(1e9)' | wc -c 1000000002 91m 7.581s using system gmp libraries 67m 40.951s using gmp libraries from /opt/perl-5.26.1/lib

        Regards, Mario

      Hello danaj and.. I'm astonished!!

      You provided a lot of high quality informations, each one worth to be investigated.

      I must admit I started this as a game, following my eclectic inclination even if math fu is stalled ~at primary school (you already know my math medieval level), my C understanding is zero and my Perl is.. still improving.

      > If you're happy using it, then great..

      Obviously no! I'm fascinated by this matter and I love Perl so my next step will be to (re-)implement the Chudnovsky algorithm in Perl just for fun and self improvement.

      I'd really like to see a blog post by you about different approach to Pi calculation. I hope it will be accessible (understandable) also for not mathematicians. In this spirit dunno if you noticed the page linked by haukex in this thread. I must admit I'm still trying to decipher it.

      For completness using integer speed up the code quite a bit. See also the interesting ambrus's post about Computing pi to multiple precision

      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.
        "...astonished!!"

        She He always performs like this. See also at github. Great.

        Update: I adressed the honorable monk wrong. Sorry and thanks to choroba. But it isn't so easy: Dana. See also he

        «The Crux of the Biscuit is the Apostrophe»

        perl -MCrypt::CBC -E 'say Crypt::CBC->new(-key=>'kgb',-cipher=>"Blowfish")->decrypt_hex($ENV{KARL});'Help

Re: porting C code to Perl
by Monk::Thomas (Friar) on Oct 23, 2017 at 16:48 UTC

    Hello

    I commented some parts of the code to make it clearer what the C code is trying to do.

    #define N 100 int len = floor(10 * N/3) + 1; // define an array, in C the size must be specified during compile // time (just use @a and let the loop populate the array with <len> // elements) int A[len]; for(int i = 0; i < len; ++i) { A[i] = 2; } int nines = 0; int predigit = 0; for(int j = 1; j < N + 1; ++j) { // iterate over the reversed array, transform each value // and calculate q -> 'map {...} reverse @a' ? int q = 0; for(int i = len; i > 0; --i) { int x = 10 * A[i-1] + q*i; A[i-1] = x % (2*i - 1); q = x / (2*i - 1); } A[0] = q%10; q = q/10; if (9 == q) { ++nines; } else if (10 == q) { printf("%d", predigit + 1); // print '0'x$nines; for (int k = 0; k < nines; ++k) { printf("%d", 0); } predigit, nines = 0; } else { printf("%d", predigit); predigit = q; if (0 != nines) { // print '9'x$nines; for (int k = 0; k < nines; ++k) { printf("%d", 9); } nines = 0; } } } printf("%d", predigit);

      Hi Monk::Thomas,

      The predigit value may set to 0 magically with some C compilers, I'm not sure. It's not the case on the Mac.

      predigit, nines = 0;

      I changed the line to this.

      predigit = 0, nines = 0;

      Basing on Discipulus's example, I too am interested and will try a plain Perl version. Then, will try parallelization via a number sequence for j. Thank you for the fun, fellow monks.

      Regards, Mario

        The predigit value may set to 0 magically with some C compilers

        Since this value is defined as a global variable (not inside a function) and is an integral data type (not part of an object/struct) it is safe to assume 'predigit' is initialized as zero by the compiler. I guess the author of the C code wasn't sure and tried to initialize both variables (just in case) ... and got it wrong. However it doesn't matter, because the '= 0' doesn't actually add anything.

        test #1 - define at global scope

        #include <stdio.h> int value1, value2 = 9; void main(void) { printf("%d %d\n", value1, value2); }
        prints '0 9', not 9 9

        test #2 - define in function

        #include <stdio.h> void func(void) { int value1, value2 = 9; printf("%d %d\n", value1, value2); } void main(void) { func(); }
        may print '996873600 9' or '-1944996480 9' or '-51174016 9'

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: perlquestion [id://1201848]
Approved by haukex
Front-paged by Corion
help
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others having an uproarious good time at the Monastery: (6)
As of 2024-03-28 09:55 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found