// #define VALARRAY_USE

#include <ivec.h>       // MMX классы и интринсики
#if defined VALARRAY_USE
#include <valarray>     // Оставил для эстетов. Раскоменьте первую строчку
#endif
#include <vector>
#include <list>
#include <iostream>
#include <cstdlib>
#include <ctime>
#include <algorithm>

/* "Среднее арифметическое без переполнения" в лоб.
   Переполнение нейтрализуется заменой сложения на вычитание, а антипереполнение -
   обменом параметров в формуле. */
unsigned avg(unsigned x, unsigned y)
{
 if (x > y) return y + (x - y) / 2;
  else return x + (y - x) / 2;
}

/* Переполнение нейтрализуется раскытием скобок.
   Компенсация возвросшей погрешности деления для случая нечётности обоих операндов. */
unsigned drgl(unsigned x, unsigned y)
{
 return (x>>1) + (y>>1) + (x&y&1);
}

/* Устранение переполнения уходом от арифметических операций к битовым.
   Для однобитных a и b справедливо a+b = c1c2, где старший бит c1 суть логическое умножение,
   а младший c2 - сложение по модулю 2. Многобитные x и y обрабатываются аналогично, а все
   их биты - одновременно, но переносы всех c1 всё равно требуется учесть арифметическим сложением.
   Деление заменено сдвигом и оптимизировано путём применения только ко всем c2, ибо незачем все c1
   сдвигать в следующий разряд, чтобы потом возвращать на место. */
unsigned mean(unsigned x, unsigned y)
{
 return (x&y) + ((x^y) >> 1);
}

/* Устранение переполнения путём увеличения диапазона промежуточного результата.
   Тут всё банально. */
unsigned port(unsigned x, unsigned y)
{
 return (static_cast<unsigned long long>(x) + y) / 2;
}

/* Устранение переполнения использованием особенностей ia32.
   CF регистра FLAGS играет роль дополнительного бита, и если его рассматривать в этом
   качестве, переполнения никогда не возникает. */
unsigned cool(unsigned x, unsigned y)
{
 __asm mov eax, [x];
 __asm add eax, [y];    // в CF самый старший (дополнительный) бит результата
 __asm rcr eax, 1;      // сдвинуть 33-битный результат (деление на 2)
}

/* Простой класс для поддержки MMX функций. Просто "деактивирует" MMX в деструкторе. */
class Guard
{
public:
 ~Guard() { empty(); }
};

/* Реализация avg() на MMX классах и интринсиках.
   Вычисления идут сразу по обоим формулам. Правильный из двух результатов выбирается в конце. */
unsigned ivec(unsigned x, unsigned y)
{
 Guard    guard;                // В конце блока будет EMMS
 Iu32vec2 m1(x, y),             // Первая формула - младший DWORD в MMX packed qword,
          m2(y, x);             // вторая         - старший DWORD.

 m1 += ((m2-m1) >> 1);

 return _MM_2UDW((x<y), m1);     // Выбор старшего или младшего DWORDа
}

/* Реализация avg() на MMX инструкциях без бранчей.
   Вычисления идут сразу по обоим формулам. Правильный из двух результатов выбирается в конце. */
unsigned mmx(unsigned x, unsigned y)
{
 __asm
 {
  mov       eax, [x]
  mov       edx, [y]
  movd      mm0, eax    // 0 x                 в mm0
  movd      mm1, edx    // 0 y                 в mm1
  cmp       eax, edx    // определение правильного DWORDа
  punpckldq mm1, mm0    // y x                 в mm1
  punpckldq mm0, mm1    // x y                 в mm0
  psubd     mm0, mm1    // x-y y-x             в mm0
  psrld     mm0, 1      // (x-y)/2 (y-x)/2     в mm0
  paddd     mm1, mm0    // y+(x-y)/2 x+(y-x)/2 в mm1
  movd      ecx, mm1    // младший DWORD
  punpckhdq mm1, mm1    // старший DWORD в младший
  movd      eax, mm1    // исходно старший DWORD
  cmovnb    eax, ecx    // выбор правильного DWORDа
  emms                  // Вуаля, бранчей нет
 }
}

#if defined VALARRAY_USE
/* Для сравнения реализация avg() на std::valarray.
   Вычисления идут сразу по обоим формулам. Правильный из двух результатов выбирается в конце.
   Реализация в лоб. Почти один в один MMXовая. */
unsigned vlr1(unsigned x, unsigned y)
{
 std::valarray<unsigned> v1(2),
                         v2(2);

 v1[0] = v2[1] = x;
 v1[1] = v2[0] = y;

 v1 += (v2-v1) / 2u;

 return v1[x>y];
}

/* Слабая попытка улучшить лобовую std::valarray-ную реализацию avg(). */
unsigned vlr2(unsigned x, unsigned y)
{
 std::valarray<unsigned> v1(2);

 v1[0] = x; v1[1] = y;

 std::valarray<unsigned> v2(v1.cshift(1));

 /* Не более одной операции за выражение, минимизация промежуточных результатов,
    все операции по месту вместо создания временных экземпляров */
 v2 -= v1;
 v2>>= 1;
 v1 += v2;

 return v1[x>y];
}
#endif

/* Далее следует тестирование вышеуказанных функций */

/* Тестируемая функция передаётся как шаблонный параметр, чтобы позволить
   её агрессивную оптимизацию.
   Пришедший параметр-пара раскладывается на параметры тестируемой функции */
template <unsigned F(unsigned, unsigned)>
unsigned doStep(const std::pair<unsigned, unsigned>& n)
{
 return F(n.first, n.second);
}

/* Простой класс для инициализации контейнеров на манер initialization-list для POD-типов.
   T - тип элементов контейнера, C - тип контейнера, по дефолту std::list */
/*template <typename T, template <typename, typename = std::allocator<T> > class C = std::list>
class Sequence
{
 C<T> cont;
*/
template <typename T, template <typename, typename> class C = std::list>
class Sequence
{
 C<T, std::allocator<T> > cont;

public:
 /* Добавление элемента в initialization-list */
 Sequence& operator, (const T& t) { cont.push_back(t); return *this; }
 /* Получить initialization-list в качестве копии контейнера затребованного типа */
 template <typename U>
 operator U() const { return U(cont.begin(), cont.end()); }
};

/* Тестовые наборы. */
std::vector<std::pair<unsigned, unsigned> > cases = (Sequence<std::pair<unsigned, unsigned> >(),
                                                      std::make_pair(123, 456),
                                                      std::make_pair(123, 457),
                                                      std::make_pair(124, 456),
                                                      std::make_pair(124, 457),

                                                      std::make_pair(456, 123),
                                                      std::make_pair(457, 123),
                                                      std::make_pair(456, 124),
                                                      std::make_pair(457, 124),

                                                      std::make_pair(4294967295, 456),
                                                      std::make_pair(0,          456),

                                                      std::make_pair(123, 4294967295),
                                                      std::make_pair(123, 0),

                                                      std::make_pair(4294967295, 4294967295),
                                                      std::make_pair(0,          0),
                                                      std::make_pair(4294967295, 0),
                                                      std::make_pair(0,          4294967295)
                                                    );

/* Размер кеша */
const size_t CASHE_SIZE = 2u * 1024u * 1024u;
/* Размер используемого для тестирования кеша */
const size_t CASHE_USING= CASHE_SIZE * 3/4;
/* Размер вектора, заполняющего максимум, но не более, чем используемый кеш */
const size_t ARRAY_SIZE = (CASHE_USING / sizeof(std::pair<unsigned, unsigned>));
/* Количество итераций для увеличения замеряемого времени */
const size_t ITER_COUNT = 1000000000u / ARRAY_SIZE;

/* Функция тестировния
   name - имя тестируемой функции; v - набор тестовых пар */
template <unsigned F(unsigned, unsigned)>
void checkTime(const char* name, const std::vector<std::pair<unsigned, unsigned> >& v)
{
 unsigned total;  // для накопления результатов во избежания черезчур умных оптимизаторов
 clock_t  start;  // время старта
 double   result; // замеренное время

 total = 0;
 start = clock();
 for (size_t j = 0; j < ITER_COUNT; ++j)
  for(std::vector<std::pair<unsigned, unsigned> >::size_type i = 0; i < v.size(); ++i)
   total ^= doStep<F>(v[i]);
 result = 1000.0 * (clock() - start) / CLOCKS_PER_SEC;

 /* total тоже выводится во избежание переоптимизации. Выводимое значение следует игнорировать */
 std::cout << name << result << '\t' << total << std::endl;
}

int main()
{
 std::vector<std::pair<unsigned, unsigned> >            v(ARRAY_SIZE);
 std::vector<std::pair<unsigned, unsigned> >::size_type i;

 /* Заполняем тестовый вектор */
 for (i = 0; i+cases.size() < v.size(); i += cases.size())
  std::copy(cases.begin(), cases.end(), v.begin() + i);
 v.resize(i);

 /* Первый пробег холостой. На всякий случай, мало ли, как там с кешем-то... */
 unsigned total = 0;
 for (i = 0; i < v.size(); ++i)
  total ^= doStep<mean>(v[i]);
 std::cout << total << std::endl;

 /* Тестируем */
 checkTime<mean>("mean: ", v);
 checkTime<avg >("avg:  ", v);
 checkTime<drgl>("drgl: ", v);
 checkTime<port>("port: ", v);
 checkTime<cool>("cool: ", v);
 checkTime<ivec>("ivec: ", v);
 checkTime<mmx >("mmx:  ", v);
#if defined VALARRAY_USE
 checkTime<vlr1>("vlr1: ", v);
 checkTime<vlr2>("vlr2: ", v);
#endif

 /* Устраняем влияние предсказателя бранчей у процессора - портим ему статистику */
 std::random_shuffle(v.begin(), v.end());

 /* Cashe fresh */
 total = 0;
 for (i = 0; i < v.size(); ++i)
  total ^= doStep<mean>(v[i]);
 std::cout << total << std::endl;

 /* Перетестируем функции, потенциально содержащие бранчи */
 checkTime<avg >("avg  (no branch prediction): ", v);
 checkTime<ivec>("ivec (no branch prediction): ", v);
#if defined VALARRAY_USE
 checkTime<vlr1>("vlr1 (no branch prediction): ", v);
 checkTime<vlr2>("vlr2 (no branch prediction): ", v);
#endif

 return 0;
}