Beefy Boxes and Bandwidth Generously Provided by pair Networks
Do you know where your variables are?
 
PerlMonks  

Longest Collatz using MCE::Flow, Inline::C, and GCC compiler intrinsics

by marioroy (Vicar)
on Apr 12, 2020 at 02:52 UTC ( #11115392=perlmeditation: print w/replies, xml ) Need Help??

Dearest Monks,

My mind has been amused with Collatz conjecture. See 1nickt's post about obtaining the top 20 sequences. Below is code for obtaining the longest progression. Here I try a GCC compiler intrinsic to further increase performance. That went well and so updated my prior post adding collatz_count_c2 there.

use strict; use warnings; use feature 'say'; use MCE::Flow; use Inline C => Config => CCFLAGSEX => '-O2 -fomit-frame-pointer', clean_after_build => 0; use Inline C => <<'END_OF_C_CODE'; #include <stdlib.h> #include <stdint.h> #if defined(_WIN32) #define strtoull _strtoui64 #endif void collatz_longest_c1( SV* _beg_n, SV* _end_n ) { uint64_t beg_n, end_n, i, n, steps; uint64_t number = 0, highest = 0; Inline_Stack_Vars; #ifdef __LP64__ beg_n = SvUV( _beg_n ); end_n = SvUV( _end_n ); #else beg_n = strtoull( SvPV_nolen( _beg_n ), NULL, 10 ); end_n = strtoull( SvPV_nolen( _end_n ), NULL, 10 ); #endif for ( i = end_n; i >= beg_n; i-- ) { n = i, steps = 0; /* count using the T(x) notation */ do { n % 2 ? ( steps += 2, n = (3 * n + 1) >> 1 ) : ( steps += 1, n = n >> 1 ); } while ( n != 1 ); if ( steps >= highest ) { number = i, highest = steps; } } Inline_Stack_Reset; Inline_Stack_Push( sv_2mortal( newSVuv(number ) ) ); Inline_Stack_Push( sv_2mortal( newSVuv(highest) ) ); Inline_Stack_Done; } void collatz_longest_c2( SV* _beg_n, SV* _end_n ) { uint64_t beg_n, end_n, i, n, steps; uint64_t number = 0, highest = 0; Inline_Stack_Vars; #ifdef __LP64__ beg_n = SvUV( _beg_n ); end_n = SvUV( _end_n ); #else beg_n = strtoull( SvPV_nolen( _beg_n ), NULL, 10 ); end_n = strtoull( SvPV_nolen( _end_n ), NULL, 10 ); #endif /* based on GCC compiler intrinsics demonstration by Alex Shirley +*/ /* https://stackoverflow.com/questions/38114205/c-collatz-conjectu +re-optimization */ /* https://www.go4expert.com/articles/builtin-gcc-functions-builti +nclz-t29238 */ for ( i = beg_n; i <= end_n; i++ ) { n = i, steps = 0; if ( n % 2 == 0 ) { steps += __builtin_ctz(n); /* account for all evens */ n >>= __builtin_ctz(n); /* always returns an odd */ } /* when we enter we're always working on an odd number */ do { n = 3 * n + 1; steps += __builtin_ctz(n) + 1; /* account for odd and even + */ n >>= __builtin_ctz(n); /* always returns an odd */ } while ( n != 1 ); if ( steps > highest ) { number = i, highest = steps; } } Inline_Stack_Reset; Inline_Stack_Push( sv_2mortal( newSVuv(number ) ) ); Inline_Stack_Push( sv_2mortal( newSVuv(highest) ) ); Inline_Stack_Done; } END_OF_C_CODE sub collatz_longest { my ( $beg_n, $end_n ) = @_; my ( $number, $highest ) = ( 0, 0 ); my ( $i, $n, $steps ); for ( $i = $end_n; $i >= $beg_n; $i-- ) { $n = $i, $steps = 0; # count using the T(x) notation $n % 2 ? ( $steps += 2, $n = (3 * $n + 1) >> 1 ) : ( $steps += 1, $n = $n >> 1 ) while $n != 1; $number = $i, $highest = $steps if ( $steps >= $highest ); } return ( $number, $highest ); } no warnings 'once'; #*collatz = \&collatz_longest; # choose collatz here #*collatz = \&collatz_longest_c1; # using T(x) notation *collatz = \&collatz_longest_c2; # using compiler intrinsics my $m = shift || '1e7'; my ( @sizes, $chunk_size ); $chunk_size = int( $m / MCE::Util::get_ncpu() / 80 + 1 ); $chunk_size += 1 if $chunk_size % 2; mce_flow_s { max_workers => MCE::Util::get_ncpu(), chunk_size => $chunk_size, gather => \@sizes, bounds_only => 1, }, sub { MCE->gather([ collatz( @{ $_[1] } ) ]); }, 1, $m; MCE::Flow->finish; # Output the longest progression for the initial starting number. # https://en.wikipedia.org/wiki/Collatz_conjecture my $highest = ( sort { $b->[1] <=> $a->[1] } @sizes )[ 0 ]->[ 1 ]; say "Longest Collatz (index and value)"; say "@$_" for ( sort { $a->[0] <=> $b->[0] } grep { $_->[1] == $highest } @sizes )[ 0..0 ];

Output

The times include launching Perl, loading modules, spawning workers, reaping workers, and output (~ 0.100 seconds).

This outputs the longest progression for the number of steps to reach 1.
These numbers are the lowest ones with the indicated step count.

1e7 : 8400511 685 1 core collatz_longest 1m16.034s collatz_longest_c1 0m01.868s collatz_longest_c2 0m00.778s 2 cores collatz_longest 0m37.912s collatz_longest_c1 0m00.965s collatz_longest_c2 0m00.422s 4 cores collatz_longest 0m19.799s collatz_longest_c1 0m00.516s collatz_longest_c2 0m00.239s 8 cores collatz_longest 0m10.042s collatz_longest_c1 0m00.285s collatz_longest_c2 0m00.147s 16 cores collatz_longest 0m05.196s collatz_longest_c1 0m00.178s collatz_longest_c2 0m00.109s 32 cores collatz_longest 0m02.717s collatz_longest_c1 0m00.137s collatz_longest_c2 0m00.105s collatz_longest_c1 (Inline C), collatz_longest (Perl) 32 cores 1e8 : 63728127 949 Inline C 0m00.738s Perl 0m30.554s 1e9 : 670617279 986 Inline C 0m07.198s Perl 5m51.938s 1e10 : 9780657630 1132 Inline C 1m17.059s 1e11 : 75128138247 1228 Inline C 13m51.122s collatz_longest_c2 (Inline C) 32 cores 1e8 : 63728127 949 Inline C 0m00.340s 1e9 : 670617279 986 Inline C 0m03.023s 1e10 : 9780657630 1132 Inline C 0m33.152s 1e11 : 75128138247 1228 Inline C 6m10.355s

I can now sleep knowing that MCE can handle this. Just be sure to use sequence generation in MCE (i.e. mce_flow_s) with the bounds_only option.

Regards, Mario

Replies are listed 'Best First'.
Re: Longest Collatz using MCE::Flow, Inline::C, and GCC compiler intrinsics (Caching)
by marioroy (Vicar) on Apr 12, 2020 at 05:47 UTC

    Hi again,

    I tried caching using an array and File::Map. All run using one core and are reasonably fast. The main goal here are attempts made at reducing memory consumption.

    Update 1: Added PDL example.
    Update 2: Improved performance.
    Update 3: Added scalar example.

    use strict; use warnings; use feature 'say'; use File::Map qw/map_anonymous unmap/; use PDL; # based on the caching demonstration by iM71 # https://stackoverflow.com/a/55361008 # https://stackoverflow.com/questions/38114205/c-collatz-conjecture-op +timization sub collatz_longest_pdl { my ( $size ) = @_; my ( $length, $number ) = ( 0, 0 ); my ( $n, $steps, $tmp, $n2 ); my $cache = zeros( short(), $size + 1 ); for my $current ( 2 .. $size ) { $n = $current, $steps = 0; # count using the T(x) notation $n % 2 ? ( $steps += 2, $n = (3 * $n + 1) >> 1 ) : ( $steps += 1, $n = $n >> 1 ) while $n != 1 && $n >= $current; $tmp = $steps + at( $cache, $n ); set( $cache, $current, $tmp ); $length = $tmp, $number = $current if $tmp > $length; } return ( $number, $length ); } sub collatz_longest_array { my ( $size ) = @_; my ( $length, $number ) = ( 0, 0 ); my ( $n, $steps, $tmp ); my @cache; for my $current ( 2 .. $size ) { $n = $current, $steps = 0; # count using the T(x) notation $n % 2 ? ( $steps += 2, $n = (3 * $n + 1) >> 1 ) : ( $steps += 1, $n = $n >> 1 ) while $n != 1 && $n >= $current; $tmp = $steps + ( $cache[ $n ] // 0 ); $cache[ $current ] = $tmp; $length = $tmp, $number = $current if $tmp > $length; } return ( $number, $length ); } sub collatz_longest_filemap { my ( $size ) = @_; my ( $length, $number ) = ( 0, 0 ); my ( $n, $steps, $tmp ); map_anonymous my $cache, $size * 2 + 2, 'shared'; # fill cache with zero's substr($cache, 0, $size * 2 + 2, pack('s', 0) x ( $size + 1 )); for my $current ( 2 .. $size ) { $n = $current, $steps = 0; # count using the T(x) notation $n % 2 ? ( $steps += 2, $n = (3 * $n + 1) >> 1 ) : ( $steps += 1, $n = $n >> 1 ) while $n != 1 && $n >= $current; $tmp = $steps + unpack('s', substr($cache, $n * 2, 2)); substr($cache, $current * 2, 2, pack('s', $tmp)); $length = $tmp, $number = $current if $tmp > $length; } unmap $cache; return ( $number, $length ); } sub collatz_longest_scalar { my ( $size ) = @_; my ( $length, $number ) = ( 0, 0 ); my ( $n, $steps, $tmp ); # fill cache with zero's my $cache = pack('s', 0) x ( $size + 1 ); for my $current ( 2 .. $size ) { $n = $current, $steps = 0; # count using the T(x) notation $n % 2 ? ( $steps += 2, $n = (3 * $n + 1) >> 1 ) : ( $steps += 1, $n = $n >> 1 ) while $n != 1 && $n >= $current; $tmp = $steps + unpack('s', substr($cache, $n * 2, 2)); substr($cache, $current * 2, 2, pack('s', $tmp)); $length = $tmp, $number = $current if $tmp > $length; } return ( $number, $length ); } #*collatz = \&collatz_longest_pdl; # choose collatz here #*collatz = \&collatz_longest_array; #*collatz = \&collatz_longest_filemap; *collatz = \&collatz_longest_scalar; my ( $number, $length ) = collatz( shift || '1e7' ); say "Longest Collatz (index and value)"; say "$number $length";

    Output

    1e7 : 8400511 685 1 core collatz_longest 1m16.034s MCE Perl code collatz_longest_pdl 0m13.691s Consumes 39 MiB collatz_longest_filemap 0m06.615s Consumes 59 MiB collatz_longest_scalar 0m06.148s Consumes 39 MiB collatz_longest_array 0m04.986s Consumes 330 MiB collatz_longest_c1 0m01.868s Inline C code collatz_longest_c2 0m00.778s Inline C code

    I've been wanting to try File::Map and pleasantly surprised.

    Regards, Mario

      Greetings,

      I am pleased to provide the final work. Computing collatz_longest for 1e9 requires ~ 3.8 GiB of available memory. Try 5e8 for lesser memory consumption.

      This is a parallel demonstration using MCE::Flow and File::Map for caching. The 2nd example uses Inline::C for counting the number of steps. Unlike the other examples, workers working on chunks 2 and higher require results from prior chunks. This is accounted for. Obviously the worker processing chunk 1 doesn't need prior results. Once the workers being processing (i.e. after the mapped cache creation), it doesn't take long before full CPU utilization kicks in.

      Q. Why does this work, especially when subsequent chunks need results from prior chunks?

      A. The magic lies with using a smaller chunk size set to 2500. The initial ramp up is a one time occurrence. One cool thing about MCE is that input IO is sequential. This applies to number sequences as well. A worker obtaining chunk 1 begins processing immediately (there's no wasting time). The worker obtaining the next chunk begins processing but may need to pause a little here and there. Eventually, the chunks as far as timing goes (starting) are spread out where workers need to pause left often. This can be seen by looking at the CPU utilization. The Power of Randomness kicks in at some point with CPU utilization near 100% until completion.

      Credit to PerlMonks choroba, Laurent_R, 1nickt, rjt, and vr. See this thread. Worthy mention goes to Leon Timmermans, author of File::Map. Wow!

      Credit for the caching technique used here is based on the caching demonstration by iM71, a response in that thread. My first attempt at parallelization failed. I tried again by maximizing on MCE's strengths described above. It's surreal :)

      Credit for reducing the number of loop iterations was from watching Notation and compressed dynamics, one minute into the video (i.e. the T(x) notation).

      Below, the minimum and maximum argument (size) is 1e6 and 1e9 respectively. The two scripts will set to limit quietly if exceeded.

      Cache miss update:

      Cache miss is less than 1%. Therefore, it is faster to compute for $n than waiting for the result.

      Final update:

      #!/usr/bin/env perl use strict; use warnings; use feature qw/say/; use File::Map qw/map_anonymous unmap/; use Time::HiRes qw/usleep/; use MCE::Flow; my $size = shift || 1e6; $size = 1e6 if $size < 1e6; # minimum $size = 1e9 if $size > 1e9; # maximum map_anonymous my $cache, $size * 2 + 2, 'shared'; # fill cache with zeroes substr($cache, 0, $size * 2 + 2, pack('s', 0) x ( $size + 1 )); # local to workers and the manager process my ( $length, $number ) = ( 0, 0 ); sub collatz_longest { my ( $chunk_id, $seq_beg, $seq_end ) = @_; my ( $n, $steps, $tmp ); for my $input ( $seq_beg .. $seq_end ) { $n = $input, $steps = 0; $n % 2 ? ( $steps += 2, $n = (3 * $n + 1) >> 1 ) : ( $steps += 1, $n = $n >> 1 ) while $n != 1 && $n >= $input; $tmp = unpack('s', substr($cache, $n * 2, 2)); # another worker with a lesser chunk_id is not yet # completed processing $n, so pause a little # if ( $tmp == 0 && $chunk_id > 1 ) { # do { # usleep 100; # $tmp = unpack('s', substr($cache, $n * 2, 2)); # } while ( $tmp == 0 ); # } # do this instead (faster): compute $n if cache miss $tmp = _collatz($n) if $tmp == 0 && $chunk_id > 1; substr($cache, $input * 2, 2, pack('s', $steps += $tmp)); $length = $steps, $number = $input if $steps > $length; } } sub _collatz { my ( $input ) = @_; my ( $n, $steps ) = ( $input, 0 ); $n % 2 ? ( $steps += 2, $n = (3 * $n + 1) >> 1 ) : ( $steps += 1, $n = $n >> 1 ) while $n != 1 && $n >= $input; my $tmp = unpack('s', substr($cache, $n * 2, 2)); $tmp = _collatz($n) if $tmp == 0; substr($cache, $input * 2, 2, pack('s', $steps += $tmp)); return $steps } my $chunk_size; $chunk_size = int( $size / MCE::Util::get_ncpu() / 80 + 1 ); $chunk_size += 1 if $chunk_size % 2; MCE::Flow->init( max_workers => MCE::Util::get_ncpu(), chunk_size => $chunk_size, # specify 2500 if pausing above bounds_only => 1, gather => sub { $length = $_[0], $number = $_[1] if $_[0] > $length; }, user_end => sub { MCE->gather($length, $number); }, ); mce_flow_s sub { my ( $mce, $seq_ref, $chunk_id ) = @_; collatz_longest($chunk_id, @{ $seq_ref }); }, 2, $size; MCE::Flow->finish; unmap $cache; say "Longest Collatz (index and value)"; say "$number $length";

      Count steps via Inline C:

      #!/usr/bin/env perl use strict; use warnings; use feature qw/say/; use File::Map qw/map_anonymous unmap/; use Time::HiRes qw/usleep/; use MCE::Flow; use Inline C => Config => CCFLAGSEX => '-O2 -fomit-frame-pointer'; use Inline C => <<'END_OF_C_CODE'; #include <stdint.h> void num_steps_c( SV* _n, SV* _s ) { uint64_t n, input; int steps = 0; n = input = SvUV(_n); while ( n != 1 && n >= input ) { n % 2 ? ( steps += 2, n = (3 * n + 1) >> 1 ) : ( steps += 1, n = n >> 1 ); } sv_setuv(_n, n); sv_setiv(_s, steps); return; } END_OF_C_CODE my $size = shift || 1e6; $size = 1e6 if $size < 1e6; # minimum $size = 1e9 if $size > 1e9; # maximum map_anonymous my $cache, $size * 2 + 2, 'shared'; # fill cache with zeroes substr($cache, 0, $size * 2 + 2, pack('s', 0) x ( $size + 1 )); # local to workers and the manager process my ( $length, $number ) = ( 0, 0 ); sub collatz_longest { my ( $chunk_id, $seq_beg, $seq_end ) = @_; my ( $n, $steps, $tmp ); for my $input ( $seq_beg .. $seq_end ) { num_steps_c($n = $input, $steps); $tmp = unpack('s', substr($cache, $n * 2, 2)); # another worker with a lesser chunk_id is not yet # completed processing $n, so pause a little # if ( $tmp == 0 && $chunk_id > 1 ) { # do { # usleep 100; # $tmp = unpack('s', substr($cache, $n * 2, 2)); # } while ( $tmp == 0 ); # } # do this instead (faster): compute $n if cache miss $tmp = _collatz($n) if $tmp == 0 && $chunk_id > 1; substr($cache, $input * 2, 2, pack('s', $steps += $tmp)); $length = $steps, $number = $input if $steps > $length; } } sub _collatz { my ( $input ) = @_; num_steps_c( my $n = $input, my $steps ); my $tmp = unpack('s', substr($cache, $n * 2, 2)); $tmp = _collatz($n) if $tmp == 0; substr($cache, $input * 2, 2, pack('s', $steps += $tmp)); return $steps } my $chunk_size; $chunk_size = int( $size / MCE::Util::get_ncpu() / 80 + 1 ); $chunk_size += 1 if $chunk_size % 2; MCE::Flow->init( max_workers => MCE::Util::get_ncpu(), chunk_size => $chunk_size, # specify 2500 if pausing above bounds_only => 1, gather => sub { $length = $_[0], $number = $_[1] if $_[0] > $length; }, user_end => sub { MCE->gather($length, $number); }, ); mce_flow_s sub { my ( $mce, $seq_ref, $chunk_id ) = @_; collatz_longest($chunk_id, @{ $seq_ref }); }, 2, $size; MCE::Flow->finish; unmap $cache; say "Longest Collatz (index and value)"; say "$number $length";

      Results: Unix time (i.e. time perl script.pl 1e9).

      1e9, 32 cores collatz_longest_final 0m26.371s collatz_longest_inline_c 0m15.634s Longest Collatz (index and value) 670617279 986

      Regards, Mario

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: perlmeditation [id://11115392]
Approved by marto
Front-paged by marto
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others about the Monastery: (9)
As of 2020-08-14 09:02 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?
    Which rocket would you take to Mars?










    Results (75 votes). Check out past polls.

    Notices?