#! perl -slw
use strict;
no warnings 'recursion';

$|++;

sub qsortpd {
    my( $ref, $lo, $hi ) = @_;

    my $pivot = unpack( 'd', substr $$ref, 8*int( ( $lo + $hi ) / 2 ), 8 );
    my( $lb, $rb ) = ( $lo, $hi );
    while( 1 ) {
        $lb++ while unpack( 'd', substr $$ref, 8*$lb, 8 ) < $pivot;
        $rb-- while unpack( 'd', substr $$ref, 8*$rb, 8 ) > $pivot;
        $lb++ while $lb != $rb
            and unpack( 'd', substr $$ref, 8*$lb, 8 )
            ==  unpack( 'd', substr $$ref, 8*$rb, 8 );
        last if $lb == $rb;
        my $temp = substr( $$ref, 8*$lb, 8 );
        substr( $$ref, 8*$lb, 8 ) = substr( $$ref, 8*$rb, 8 );
        substr( $$ref, 8*$rb, 8 ) = $temp;
    }
    qsortpd( $ref, $lo, $lb ) if $lo < --$lb;
    qsortpd( $ref, $rb, $hi ) if $hi > ++$rb;
}

our $N ||= 1e6;

my $packedDoubles = '';

open RAM, '>', \$packedDoubles;
seek RAM, $N*8, 0;
print RAM chr(0);
seek RAM, 0, 0;
printf RAM "%s", pack 'd1000', map{ rand 32767 } 1 .. 1000 for 1 .. $N / 1000;
close RAM;

print "Starting sort of $N packed doubles: ", time;
qsortpd( \$packedDoubles, 0, $N-1 );
print "Finished sort of $N packed doubles: ", time;

<>;
print for unpack 'd*', $packedDoubles;

##</code><code>##

#! perl -slw
use strict;
use Devel::Peek;
use Inline C => Config => BUILD_NOISY => 1;
use Inline C => <<'END_C',  NAME => 'QsortPD', CLEAN_AFTER_BUILD => 0;

#include <stdlib.h>
#include <string.h>
#include <stdio.h>

int compare( const void *arg1, const void *arg2 ) {
    double diff = *(double*)arg1 - *(double*)arg2;
    return diff < 0 ? -1 : diff > 0 ? 1 : 0;
}

int qsortPackedDoubles( SV *data, SV *n ) {
    qsort( (void *)SvPVX( data ), SvUV( n ), sizeof( double ), compare );
}

END_C
our $N ||= 1e6; ## Warning N less that 1000 will give funny results

my $packedDoubles = '';

## Preallocate the scalar and random fill with packed doubles
open RAM, '>', \$packedDoubles;
seek RAM, $N*8, 0;
print RAM chr(0);
seek RAM, 0, 0;
printf RAM "%s", pack 'd1000', map{ rand 32767 } 1 .. 1000 for 1 .. $N / 1000;
close RAM;

print "Starting sort of $N packed doubles: ", time;
qsortPackedDoubles( $packedDoubles, $N );
print "Finished sort of $N packed doubles: ", time;

<>;
print for unpack 'd*', $packedDoubles;