// #define VALARRAY_USE

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

/* "Среднее арифметическое без переполнения" в лоб. Сравнение знаков операндов оптимизировано.
   Для отрицательной разности y-x результат может быть округлён в сторону,
   противоположную ожидаемой по Стандарту. Потому деление заменено сдвигом,
   чтобы получить округление к минус бесконечности, и применена "стандартаная"
   коррекция: если (результат отрицательный) && (операнды имеют разную чётность) */
signed avg(signed x, signed y)
{
 signed tmp;

 if ((x^y) < 0) tmp = (x + y) >> 1;
  else tmp = x + ((y - x) >> 1);

 return tmp + (tmp<0 && (x^y)&1);
}

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

 return tmp + (tmp<0 && (x^y)&1);
}

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

 return tmp + (tmp<0 && (x^y)&1);
}

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

/* Устранение переполнения использованием особенностей ia32.
   Переполнение, если возникает, инвертит знаковый бит и детектируется флагом OF во FLAGS. */
signed cool(signed x, signed y)
{
 __asm
 {
  xor  eax, eax
  mov  ecx, [x]
  mov  edx, ecx
  add  ecx, [y]
  seto al               // Получаем признак переполнения
  sar  ecx, 1           // Делим на 2 (и сохраняем текущий знак)
  ror  eax, 1           // Перемещаем признак переполнения в знаковый бит
  xor  edx, [y]         // x^y
  and  edx, 1           // (x^y)&1
  xor  eax, ecx         // Восстанавливаем правильный знак
  sets cl               // tmp<0 (весь ecx нулить не нужно: edx и так будет либо 0, либо 1)
  and  edx, ecx         // Готовая компенсация для отрицательного
  add  eax, edx         // Вуаля, бранчей нет
 }
}

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

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

 m3 += cmpgt(Is32vec2(0), m3) & (m1^m2) & Is32vec2(1, 1); // Компенсация при отрицательном

 // Выбор старшего при разных знаках операндов, младшего при наоборот
 return _MM_2DW(((x^y) < 0), m3);
}

/* Реализация avg() на MMX инструкциях без бранчей.
   Вычисления идут сразу по обоим формулам. Правильный из двух результатов выбирается в конце. */
signed mmx(signed x, signed y)
{
 __asm
 {
  mov       edx, [y]
  mov       eax, [x]
  movd      mm0, edx    // 0 y               в mm0
  movd      mm1, eax    // 0 x               в mm1
  pxor      mm2, mm2    // 0 0               в mm2
  psubd     mm2, mm1    // 0 -x              в mm2
  punpckldq mm0, mm0    // y y               в mm0
  punpckldq mm2, mm1    // x -x              в mm2
  movq      mm3, mm2    //                 и в mm3 на будущее
  mov       ecx, 1      // Подготовка единиц в mm2
  paddd     mm2, mm0    // y+x y-x           в mm2
  psrad     mm2, 1      // (y+x)/2 (y-x)/2   в mm2 (деление заменено сдвигом)
  paddd     mm1, mm2    // (y+x)/2 x+(y-x)/2 в mm1
  pxor      mm4, mm4    // 0 0               в mm4
  movd      mm2, ecx    // 0 1               в mm2
  pcmpgtd   mm4, mm1    // tmp (hi и lo) < 0 в обоих DWORDах mm4 (0 или -1)
  pxor      mm3, mm0    // y^x               в обоих DWORDах mm3 (знак x пофигу, важен младший бит)
  punpckldq mm2, mm2    // 1 1               в mm2
  pand      mm4, mm3    // tmp<0 & (y^x)     в обоих DWORDах mm4 (0 или не 0)
  pand      mm4, mm2    // компенсации в mm4 (0 или 1) готовы
  xor       eax, edx    // знаки операндов
  paddd     mm1, mm4    // компенсации из mm4 учтены
  movd      eax, mm1    // младший DWORD
  punpckhdq mm1, mm1    // старший DWORD в младший
  movd      edx, mm1    // исходно старший DWORD
  cmovs     eax, edx    // выбор правильного DWORDа
  emms                  // Вуаля, бранчей нет
 }
}

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

 v2[0] =-v2[0];
 v3[0] = x;

 v3 += (v2+v1) >> 1;
 v3 += (v3 >> 31) & (v1^v2) & 1; // Сравнения возвращают несовместимый std::valarray<bool>

 return v3[(x^y) < 0];
}

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

 v2[0] =-v2[0];
 v3[0] = x;

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

 v2 ^= v1;
 v2 &= 1;
 v1  = t1;
 v1>>= 31;
 v1 &= v2;
 v1 += t1;

 return v1[(x^y) < 0];
}
#endif

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

/* Тестируемая функция передаётся как шаблонный параметр, чтобы позволить
   её агрессивную оптимизацию.
   Пришедший параметр-пара раскладывается на параметры тестируемой функции */
template <signed F(signed, signed)>
signed doStep(const std::pair<signed, signed>& 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<signed, signed> > cases = (Sequence<std::pair<signed, signed> >(),
                                                  std::make_pair( 123, 456),
                                                  std::make_pair( 123, 457),
                                                  std::make_pair( 124, 456),
                                                  std::make_pair( 124, 457),

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

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

                                                  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( 456,-123),
                                                  std::make_pair( 457,-123),
                                                  std::make_pair( 456,-124),
                                                  std::make_pair( 457,-124),

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

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

                                                  std::make_pair( 2147483647,-456),
                                                  std::make_pair( 2147483647, 456),
                                                  std::make_pair(-2147483648,-456),
                                                  std::make_pair(-2147483648, 456),

                                                  std::make_pair( 123, 2147483647),
                                                  std::make_pair(-123, 2147483647),
                                                  std::make_pair( 123,-2147483648),
                                                  std::make_pair(-123,-2147483648),

                                                  std::make_pair( 2147483647, 2147483647),
                                                  std::make_pair(-2147483648, 2147483647),
                                                  std::make_pair( 2147483647,-2147483648),
                                                  std::make_pair(-2147483648,-2147483648),

                                                  std::make_pair( 2147483647, 0),
                                                  std::make_pair(-2147483648, 0),

                                                  std::make_pair( 0, 2147483647),
                                                  std::make_pair( 0,-2147483648)
                                                );

/* Размер кеша */
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<signed, signed>));
/* Количество итераций для увеличения замеряемого времени */
const size_t ITER_COUNT = 1000000000u / ARRAY_SIZE;

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

 total = 0;
 start = clock();
 for (size_t j = 0; j < ITER_COUNT; ++j)
  for(std::vector<std::pair<signed, signed> >::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<signed, signed> >            v(ARRAY_SIZE);
 std::vector<std::pair<signed, signed> >::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);              // удаляем избыток

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

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

 return 0;
}