// The Funnelsort Project  -  Cache oblivious sorting algorithm implementation
// Copyright (C) 2005  Kristoffer Vinther
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA


#include <cassert>
#include <algorithm>
#include <utility>
#include <memory>
#include <iterator>
#include <stdexcept>

namespace iosort
{

namespace base_
{

enum tag_RECURSELIMIT { RECURSELIMIT = 24 };

template<class It, class Pred>
inline void insertion_sort(It begin, It end, Pred comp)
{
	for( It j=begin+1; j<end; j++ )
	{
		typename std::iterator_traits<It>::value_type e = *j;
		for( It i=j; begin<i; *i=e )
			if( comp(e,*(i-1)) )
				*i = *(i-1), --i;
			else
				break;
	}
}

template<typename T, class Pred>
inline T median(const T& a, const T& b, const T& c, Pred comp)
{
	if( comp(a,b) )
		if( comp(b,c) )
			return b;
		else if( comp(a,c) )
			return c;
		else
			return a;
	else
		if( comp(a,c) )
			return a;
		else if( comp(b,c) )
			return c;
		else
			return b;
}

template<class BidIt, class T, class Pred>
inline BidIt partition(BidIt begin, BidIt end, T pivot, Pred comp)
{
	for( ;; )
	{
		for( ; comp(*begin,pivot); ++begin );
		for( --end; comp(pivot,*end); --end );
		if( !(begin<end) )
			return begin;
		std::iter_swap(begin, end);
		++begin;
	}
}

template<class RanIt, class Pred>
inline void quicksort(RanIt First, RanIt Last, Pred comp)
{
	typename std::iterator_traits<RanIt>::difference_type Count = Last-First;
	while( RECURSELIMIT < Count )
	{
		RanIt part = partition(First, Last, median(*First,*(First+(Last-First)/2),*(Last-1),comp),comp);
		if( part-First < Last-part ) // loop on larger half
			quicksort(First, part, comp), First = part;
		else
			quicksort(part, Last, comp), Last = part;
		Count = Last-First;
	}
	if( Count > 1 )
        insertion_sort(First, Last, comp);
}

template<class Int>
inline Int logc(Int k)
{
	height_t h = 0;
	for( order_t i=k-1; i; h++, i/=2 );
	return h;
}

template<class RanIt, class Pred>
void batcher_sort(RanIt begin, RanIt end, Pred comp)
{
	typedef typename std::iterator_traits<RanIt>::difference_type ItDiff;
	ItDiff t = logc(end-begin)-1;
	ItDiff p = 1<<t;
	for( ; p!=1; p=p/2 )
	{
		RanIt i, j;
		for( i=begin; i<end-p-p; i+=p )
		{
			j = i+p;
			for( ; i<j; i++ )
				if( comp(*(i+p),*i) )
					std::iter_swap(i+p, i);
		}
		for( ; i<end-p; i++ )
			if( comp(*(i+p),*i) )
				std::iter_swap(i+p, i);
		ItDiff q = 1<<t;
		while( q != p )
		{
			ItDiff d=q-p;
			q = q/2;
			for( i=begin+p; i<end-d-p; i+=p )
			{
				j = i+p;
				for( ; i<j; i++ )
					if( comp(*(i+d),*i) )
						std::iter_swap(i+d, i);
			}
			for( ; i<end-d; i++ )
				if( comp(*(i+d),*i) )
					std::iter_swap(i+d, i);
		}
	}
	if( p )
	{
		for( RanIt i=begin; i<end-1; i+=2 )
			if( comp(*(i+1),*i) )
				std::iter_swap(i+1, i);
		ItDiff q = 1<<t;
		while( q != 1 )
		{
			ItDiff d=q-1;
			q = q/2;
			for( RanIt i=begin+1; i<end-d; i+=2 )
				if( comp(*(i+d),*i) )
					std::iter_swap(i+d, i);
		}
	}
}

template<class RanIt, class Pred>
void inplace_base_sort(RanIt begin, RanIt end, Pred comp)
{
#ifdef __ia64__
	batcher_sort(begin, end, comp);
#else // __ia64__
	quicksort(begin, end, comp);
#endif // __ia64__
}

} // namespace base_

template<class RanIt, class Splitter>
class stream_slices
{
	typedef typename std::iterator_traits<RanIt>::difference_type diff;
public:
	class iterator : public std::iterator<std::bidirectional_iterator_tag, std::pair<RanIt,RanIt>, diff>
	{
		friend class stream_slices;
	public:
		inline iterator& operator++()
		{
			begin = end;
			end = end+run_size;
			if( --big_runs == 0 )
				--run_size;
			return *this;
		}
		inline iterator& operator--()
		{
			end = begin;
			begin = end-run_size;
			if( ++big_runs == 0 )
				++run_size;
			return *this;
		}
		std::pair<RanIt,RanIt> operator*()
			{ return make_pair(begin,end); }
		bool operator==(const iterator& rhs)
			{ return begin == rhs.begin; }
	private:
		inline iterator(RanIt begin, diff run_size, size_t big_runs) : begin(begin), run_size(run_size), big_runs(big_runs)
		{
			if( big_runs )
				++run_size;
			end = begin+run_size;
		}
		RanIt begin, end;
		diff run_size;
		size_t big_runs;
	};
	stream_slices(RanIt begin, RanIt end) : b(begin), e(end), order_(Splitter::prefered_order_output(end-begin))
	{
		if( order_ < 2 )
			throw std::logic_error("Splitter::prefered_order_output returned less than two");
		run_size = (end-begin)/order_;
		size_t rem = static_cast<size_t>((end-begin) % order_);
		run_size += rem/order_;
		big_runs = rem % order_;
	}
	inline diff order() const
		{ return order_; }
	inline iterator begin()
		{ return iterator(b, run_size, big_runs); }
	inline iterator end()
		{ return iterator(e, run_size, big_runs); }
private:
	RanIt b, e;
	diff run_size;
	size_t order_, big_runs;
};

template<class Splitter, class Diff>
inline Diff run_size_(Diff d)
{
	size_t order = Splitter::prefered_order_output(d);
	Diff run_size = d/order;
	size_t rem = static_cast<size_t>(d % order);
	run_size += rem/order;
	size_t big_runs = rem % order;
	if( big_runs )
		run_size++;
	return run_size;
}

template<class Merger, class Splitter, class Diff, class MAlloc, class TAlloc>
inline std::pair<Merger*,Merger*> build_cache_(Diff run_size, MAlloc& malloc, TAlloc& talloc)
{
	Merger *begin_cache, *end_cache;
	int rec_calls = 0;
	for( Diff d=run_size; d>static_cast<Diff>(Splitter::switch_to_cache_aware()); d=run_size_<Splitter>(d), ++rec_calls );
	begin_cache = end_cache = malloc.allocate(rec_calls);
	for( Diff d=run_size; d>static_cast<Diff>(Splitter::switch_to_cache_aware()); d=run_size_<Splitter>(d), ++end_cache )
		new(end_cache) Merger(Splitter::prefered_order_output(d), talloc);
	return make_pair(begin_cache,end_cache);
}

template<class Merger, class Splitter, class It, class OutIt, class Alloc>
void merge_sort_(It begin, It end, OutIt dest, Merger* cache, Alloc alloc);

template<class Merger, class Splitter, class MPtr, class It, class Alloc>
inline void sort_streams_(MPtr cache, It begin, size_t order, size_t run_size, size_t big_runs, Alloc alloc)
{
	typedef typename std::iterator_traits<It>::value_type T;
	typename Alloc::template rebind<T>::other talloc(alloc);
	typedef typename Alloc::template rebind<T>::other::pointer TPtr;
	TPtr tmp = talloc.allocate((size_t)(run_size));
	It last = begin;
	size_t i;
	if( big_runs )
	{
		merge_sort_<Merger,Splitter>(last, last+run_size, std::raw_storage_iterator<TPtr,T>(tmp), cache, alloc);
		for( i=1, last+=run_size; i!=big_runs; i++, last+=run_size )
			merge_sort_<Merger,Splitter>(last, last+run_size, last-run_size, cache, alloc);
		run_size--;
		for( ; i!=order; i++, last+=run_size )
			merge_sort_<Merger,Splitter>(last, last+run_size, last-run_size-1, cache, alloc);
		run_size++;
	}
	else
	{
		merge_sort_<Merger,Splitter>(last, last+run_size, tmp, cache, alloc);
		for( --order, last+=run_size; order; --order, last+=run_size )
			merge_sort_<Merger,Splitter>(last, last+run_size, last-run_size, cache, alloc);
	}
	std::copy(tmp, tmp+run_size, last-run_size);
	for( TPtr p=tmp+run_size-1; p>=tmp; --p )
		talloc.destroy(p);
	talloc.deallocate(tmp,static_cast<size_t>(run_size));
}

template<class MPtr, class It>
inline void add_sorted_streams_(MPtr merger, It end, size_t order, size_t run_size, size_t big_runs)
{
	if( big_runs )
	{
		merger->add_stream(end-run_size, end);
		end -= run_size;
		--run_size, --big_runs;
		for( --order; order!=big_runs; --order, end-=run_size )
			merger->add_stream(end-run_size, end);
		++run_size;
		for( ; order; --order, end-=run_size )
			merger->add_stream(end-run_size, end);
	}
	else
		for( ; order; --order, end-=run_size )
			merger->add_stream(end-run_size, end);
}

template<class MPtr, class It, class Pred>
inline void add_streams_(MPtr merger, It begin, size_t order, size_t run_size, size_t big_runs, Pred comp)
{
	size_t i;
	for( i=0; i!=order-big_runs; ++i, begin+=run_size )
		base_::inplace_base_sort(begin, begin+run_size, comp);
	++run_size;
	for( ; i!=order; ++i, begin+=run_size )
		base_::inplace_base_sort(begin, begin+run_size, comp);

	for( ; i!=order-big_runs; --i, begin-=run_size )
		merger->add_stream(begin-run_size, begin);
	--run_size;
	for( ; i; --i, begin-=run_size )
		merger->add_stream(begin-run_size, begin);
}

template<class Merger, class Splitter, class It, class OutIt, class Alloc>
void merge_sort_(It begin, It end, OutIt dest, Merger* cache, Alloc alloc)
{
	typedef typename std::iterator_traits<It>::value_type T;
	typedef typename std::iterator_traits<It>::difference_type ItDiff;
	typedef typename Merger::stream Stream;

	order_t order = Splitter::prefered_order_output(end-begin);
	if( order < 2 )
		throw std::logic_error("Splitter::prefered_order_output returned less than two");
	ItDiff run_size = (end-begin)/order;
	size_t rem = (size_t)((end-begin) % order);
	run_size += rem/order;
	size_t big_runs = rem % order;

	cache->reset();
	if( run_size > static_cast<ItDiff>(Splitter::switch_to_cache_aware()) )
	{
		if( big_runs )
			run_size++;
		sort_streams_<Merger,Splitter>(cache+1,begin,order,run_size,big_runs,alloc);
		add_sorted_streams_(cache,end,order,run_size,big_runs);
	}
	else
		add_streams_(cache,begin,order,run_size,big_runs, typename Merger::predicate());
	(*cache)(dest);
}

template<class Merger, class Splitter, class It, class OutIt>
inline void merge_sort(It begin, It end, OutIt dest, const typename Merger::allocator& alloc)
{
	typedef typename std::iterator_traits<It>::difference_type ItDiff;
	if( end-begin > static_cast<ItDiff>(Splitter::switch_to_cache_aware()) )
	{
		order_t order = Splitter::prefered_order_output(end-begin);
		if( order < 2 )
			throw std::logic_error("Splitter::prefered_order_output returned less than two");
		ItDiff run_size = (end-begin)/order;
		size_t rem = (size_t)((end-begin) % order);
		run_size += rem/order;
		size_t big_runs = rem % order;
		if( run_size > static_cast<ItDiff>(Splitter::switch_to_cache_aware()) )
		{
			if( big_runs )
				run_size++;
			typename Merger::allocator::template rebind<Merger>::other malloc(alloc);
			std::pair<Merger*,Merger*> cache = build_cache_<Merger,Splitter>(run_size, malloc, alloc);
			sort_streams_<Merger,Splitter>(cache.first, begin, order, run_size, big_runs, alloc);
			for( Merger *pm=cache.second-1; pm>=cache.first; --pm )
				malloc.destroy(pm);
			malloc.deallocate(cache.first, cache.second-cache.first);
		}
		Merger merger(order,alloc);
		if( run_size > static_cast<ItDiff>(Splitter::switch_to_cache_aware()) )
			add_sorted_streams_(&merger, end, order, run_size, big_runs);
		else
			add_streams_(&merger, begin, order, run_size, big_runs, typename Merger::predicate());
#ifdef _DEBUG
		OutIt e = merger(dest, dest+(end-begin));
		merger.empty(dest, dest+(end-begin));
		size_t N = 0;
		for( typename Merger::stream_iterator i=merger.begin(); i!=merger.end(); ++i )
			N += (*i).end()-(*i).begin();
		assert( e == dest+(end-begin) );
#else // _DEBUG
		merger(dest);
#endif // _DEBUG
	}
	else
	{
		base_::inplace_base_sort(begin, end, typename Merger::predicate());
		std::copy(begin, end, dest);
	}
}
} // namespace iosort