("long read":))

Here's an attempt to solve this Collatz challenge in purely vectorized fashion. As it (slowly) progressed, it soon became apparent the emerging result is somewhat impractical, merely a curiosity; now it's finished (i.e. I can't find how to improve it) and ready to amaze/amuse/frighten the public.

The "PDL" was mentioned in this (or 2 related) threads, but it served just as binary container, "set/at" having same role as "substr + pack/unpack", so it wasn't really "PDL", was it. This node sticks to original task as stated:

calculate the sequence length for all starting numbers up to 1000000 (1e6), and output the starting number and sequence length for the longest 20 sequences

The naive straightforward implementation is simple, it should read without much comment (as "prose":)):

use strict;
use warnings;
use feature 'say';
use PDL;

use Time::HiRes 'time';
my $t = time;

use constant MAX => 1e6;
use constant TOP => MAX < 20 ? MAX : 20;

my $seqs = 1 + sequence( longlong, MAX );
my $lengths = ones( short, MAX );

while ( any my $good_mask = $seqs-> inplace
                                 -> setvaltobad( 1 )
                                 -> isgood ) {
    my $odd_mask = $seqs & 1;
    $lengths-> where( $odd_mask ) ++;
    $lengths-> where( $good_mask ) ++;
    ( $seqs-> where( $odd_mask ) *= 3 ) ++; 
    $seqs >>= 1; 
}

my $top_i = $lengths-> qsorti
                    -> slice([ MAX - 1, MAX - TOP ]);

say $lengths-> index( $top_i )
            -> longlong
            -> cat( $top_i + 1 )
            -> transpose;
say time - $t;

__END__

[
 [   525 837799]
 [   509 626331]
...
 [   445 938143]
 [   445 906175]
 [   445 922525]
 [   445 922526]
]

7.98023009300232
[download]

Above, completed sequences are marked as "BAD", which is just agreed-upon value, treated specially. (1) There are built-in methods to check for good/bad, it saves us comparisons while creating masks. (2) More important: BADs kind of taint whatever they interact with (cf. NaN in FP math), which is quite useful.

The timing is sloooooow, yet decent among non-caching solutions around here. Another issue stems from QuickSort being unstable (I recently ranted, in this thread, about neat ordering and "correct 20" extraction).

To fix the latter, there's qsortvec (and qsortveci companion) method to sort 2D array (as opposed to "qsort" for 1D, used above), i.e. 1st on 1st axis, then on 2nd. But here's dilemma: (1) build full-height (2 x 1e6) array, qsortvec, extract top-20. Possible, but, for speed, I'd prefer (2): qsort lengths (as above), extract "many enough" but close to 20, build small "2 x N" array, qsortvec, extract correct (and correctly arranged) top-20.

For that, find value at MAX - TOP (445), look left, find how much to extract (22). More fuss: 1st column is to descend, 2nd to ascend -- thus temporarily negate (flip bits) one of them. So, huge and unpleasantly looking new "tail" after main loop, in script below, is there to fix top-20 extraction. But in fact it adds almost nothing to consumed time.

To improve speed, there are couple of tricks. (1) Where has neat ability to work in list assignment -- same mask for several client piddles. (2) Marking 1's as "BAD" on each iteration is redundant. Piddle can be told to treat any value as bad, automatically. Sadly, setting this "any" to "1" won't work here, because $seqs & 1 would then result in something like [BAD 0 BAD 0 BAD 0 ...], regardless of having BADs in $seqs already. Let's mark "2" as bad, so stopper value in sequence would now be 2 instead of 1:

[ BAD BAD 3 4 5 ... ]       # initial $seqs
[ 1 2 2 2 2 ... ]           # initial $lengths
[download]

Other than that, it's the same non-caching approach, with original clarity and simplicity somewhat spoiled by fixes/optimizations:

use strict;
use warnings;
use feature 'say';
use PDL;

use Time::HiRes 'time';
my $t = time;

use constant MAX => 1e6;
use constant TOP => MAX < 20 ? MAX : 20;

my $seqs = 1 + sequence( longlong, MAX );
$seqs-> setbadat( 0 );
$seqs-> badvalue( 2 );

my $lengths = ones( short, MAX );
$lengths <<= 1;
$lengths-> set( 0, 1 );

while ( any my $good_mask = $seqs-> isgood ) {

    my ( $seqs_odd, $lengths_odd_masked )
        = where( $seqs, $lengths, $seqs & 1 );

    $lengths_odd_masked ++;
    $lengths-> where( $good_mask ) ++;
    ( $seqs_odd *= 3 ) ++;
    $seqs >>= 1; 
}

my $sorted_i = $lengths-> qsorti;
my $sorted   = $lengths-> index( $sorted_i );
my $value    = $sorted-> at( MAX - TOP );
my $pos      = vsearch_insert_leftmost( $value, $sorted );
my $top_i    = $sorted_i-> slice([ MAX - 1 , $pos ]);

( my $result = $lengths
    -> index( $top_i )
    -> longlong
    -> bitnot
    -> cat( $top_i + 1 )
    -> transpose
    -> qsortvec
    -> slice([], [ 0, TOP - 1 ])

)-> slice([ 0 ], [])
 -> inplace
 -> bitnot;

say $result;
say time - $t;

__END__

[
 [   525 837799]
 [   509 626331]
...
 [   445 886953]
 [   445 906175]
 [   445 922524]
 [   445 922525]
]

6.0809600353241
[download]

So we have correct output and improved time. Good.

Now to something more interesting -- let's add caching/looking-up. Because we are to use $seqs as index into $lengths, and indexing starts from 0, let's prepend a dummy 0th element. To kick-start indexing, and because value "2" is occupied to mark "BAD" i.e. sequence stopper, we'll add one more seed element to lengths. Further, lengths will now all start as BAD, and switched to computed values as we go:

[ BAD BAD BAD 3 4 5 6 ... ]     # initial $seqs
[ BAD 1 2 BAD 3 BAD BAD ... ]   # initial $lengths
[download]

(by the way, BAD in $lengths is still the default, for short, -32768)

We'll also maintain $current lengths helper piddle, incremented by 1 or 2 depending on oddity mask of current sequences state. Observe, further, how where calls in list context are (over)-abused in code below. (I'm quite aware this code is no longer "a prose to read". Set MAX to 10 and dump primary piddles on each iteration to see what's going on. There are same 3. Other vars are masked views into them.)

use strict;
use warnings;
use feature 'say';
use PDL;

use Time::HiRes 'time';
my $t = time;

use constant MAX => 1e6;
use constant TOP => MAX < 20 ? MAX : 20;

my $seqs = sequence( longlong, 1 + MAX );
$seqs-> setbadat( 0 );
$seqs-> setbadat( 1 );
$seqs-> badvalue( 2 );

my $lengths = ones( short, 1 + MAX );
$lengths-> inplace-> setvaltobad( 1 );
$lengths-> set( 1, 1 );
$lengths-> set( 2, 2 );
$lengths-> set( 4, 3 );

my $current = zeroes( short, 1 + MAX );

while ( any $seqs-> isgood ) {      # sic

    my ( $seqs_odd, $current_odd_masked )
        = where( $seqs, $current, $seqs & 1 );

    $current_odd_masked ++;
    $current ++;
    ( $seqs_odd *= 3 ) ++;
    $seqs >>= 1; 
    
    my ( $seqs_cap, $lengths_cap, $current_cap )
        = where( $seqs, $lengths, $current, $seqs <= MAX );

    my $lut = $lengths-> index( $seqs_cap );
    
    # "_f" is for "finished"
    
    my ( $seqs_f, $lengths_f, $lut_f, $current_f ) 
        = where( $seqs_cap, $lengths_cap, $lut, $current_cap, 
                 $lut-> isgood );

    $lengths_f .= $lut_f + $current_f;
    $seqs_f    .= 2;                    # i.e. BAD
}

$lengths-> badflag( 0 );

my $sorted_i = $lengths-> qsorti;
my $sorted   = $lengths-> index( $sorted_i );
my $value    = $sorted-> at( MAX + 1 - TOP );
my $pos      = vsearch_insert_leftmost( $value, $sorted );
my $top_i    = $sorted_i-> slice([ MAX, $pos ]);

( my $result = $lengths
    -> index( $top_i )
    -> longlong
    -> bitnot
    -> cat( $top_i )
    -> transpose
    -> qsortvec
    -> slice([], [ 0, TOP - 1 ])

)-> slice([ 0 ], [])
 -> inplace
 -> bitnot;

say $result;
say time - $t;

__END__

[
 [   525 837799]
 [   509 626331]
...
 [   445 886953]
 [   445 906175]
 [   445 922524]
 [   445 922525]
]

2.88385105133057
[download]

And that's (~2x faster) I'm afraid is as good as it will go. As I understand, cache hits are significantly more rare than with consequential element after element array processing. For comparison, with the same hardware, Laurent_R's final/polished caching solution runs at 1.63s here, if I disable use of "magic number" 400 in there ("magic" constants to crank up performance aren't fair:)), and at 0.88s otherwise.

For 1e7 numbers, running time becomes ~65s, i.e. it gets impractical, like I said. Switching on parallel processing made no difference. Use of multiple cores can be observed for only ~first second, then "viewports" into piddles become increasingly fragmented, work can't be split.

Maybe I did something wrong, and certainly someone can improve even if a little bit, but I'm glad this mini-project is finally off my shoulders.