Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
In-place array interleaving
#include <iostream>
#include <vector>
#include <algorithm>
#include <iterator>
#include <assert.h>
typedef std::chrono::high_resolution_clock Clock;
template<typename T>
void print( const std::vector<T>& V ) {
std::copy(V.begin(), V.end(), std::ostream_iterator<T>(std::cout,","));
std::cout << "\n";
}
//! Return the highest power of 3 less or equal to the array lenght (2n)
size_t highest_leader( size_t len )
{
size_t i = 1;
while ( i * 3 < len ) i *= 3;
return i;
}
//! Permutation (sigma) function returning the next index in the cycle. Permutation indices start at 1.
//inline size_t next_index( size_t j, size_t len )
//{
// return ( j * 2 ) % ( len + 1 );
//}
// Optimization suggested by @johanofverstedt. 1291ms to 703ms with -03
inline size_t next_index( size_t j, size_t len )
{
j *= 2;
size_t m = len + 1;
while(j >= m)
j -= m;
return j;
}
//! Apply the cycle leader algorithm for the in-shuffle permutation
template<typename T>
void cycle_rotate( std::vector<T>& A, size_t leader, size_t section_offset, size_t section_len )
{
size_t i = next_index(leader, section_len);
while( i != leader ) {
std::swap( A[section_offset + i - 1], A[section_offset + leader - 1] );
i = next_index(i, section_len);
}
}
template<typename T>
void interleave_inplace( std::vector<T>& A )
{
if (A.size() <= 2)
return;
// Array size should be even
assert( A.size() % 2 == 0 );
auto it = A.begin();
size_t m=0, n=1;
while( m < n ) {
size_t section_len = std::distance( it, A.end() );
size_t h = highest_leader(section_len);
m = ( h > 1 ) ? h/2 : 1;
n = section_len/2;
//right cyclic shift
std::rotate( it + m, it + n, it + m + n );
for( size_t leader = 1; leader < 2*m; leader *= 3 ) {
cycle_rotate<T>( A, leader, std::distance( A.begin(), it ), 2*m );
}
//move iterator to remaining un-interleaved section
std::advance(it, 2*m);
}
}
template <typename T>
void interleave( std::vector<T>& A )
{
std::vector<T> B(A);
size_t half = A.size() / 2;
for( size_t i =0; i < B.size(); i += 2 ) {
A[i] = B[half + i/2];
A[i+1] = B[i/2];
}
}
int main(int argc, const char * argv[])
{
int length = 500000;
int nb_applications = 200;
std::vector<int> A(length);
for(int i=0; i<A.size(); i++) {
A[i] = ( i < A.size()/2 ) ? -i-1 : int(i - A.size()/2)+1;
}
std::vector<int> B(A);
{
std::cout << "Trivial copy algo: ";
auto start = Clock::now();
for( int i=0; i < nb_applications; ++i ) {
interleave( A );
}
auto end = Clock::now();
auto nsdur = std::chrono::duration_cast<std::chrono::nanoseconds>( end - start );
auto msdur = std::chrono::duration_cast<std::chrono::milliseconds>( end - start );
std::cout << nsdur.count() << "ns, " << msdur.count() << "ms" << std::endl;
}
{
std::cout << "Inplace algo: ";
auto start = Clock::now();
for( int i=0; i < nb_applications; ++i ) {
interleave_inplace( B );
}
auto end = Clock::now();
auto nsdur = std::chrono::duration_cast<std::chrono::nanoseconds>( end - start );
auto msdur = std::chrono::duration_cast<std::chrono::milliseconds>( end - start );
std::cout << nsdur.count() << "ns, " << msdur.count() << "ms" << std::endl;
}
assert( A == B );
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.