Сказка о потерянном времени

Если честно, то не совсем и сказка, а суровая жизнь. Но время ведь потеряно совершенно настоящее, хоть и с пользой. А началось всё совершенно случайно. На одном сайте один умный товарищ написал пост ^[1] о гипотезе Эйлера ^[2]. Суть достаточно проста. Гипотеза Эйлера утверждает, что для любого натурального числа n>2 никакую n-ю степень натурального числа нельзя представить в виде суммы (n-1) n-х степеней других натуральных чисел. То есть, уравнения:

не имеют решения в натуральных числах.

Ну собственно так оно и было до 1966 года…

Пока Л. Ландер (L. J. Lander), Т. Паркин (T. R. Parkin) и Дж. Селфридж ( J. L. Selfridge) не нашли первый контрпример для n = 5. И сделали они это на суперкомпьютере того времени — CDC 6600 ^[3], разработанного под командованием не безызвестного Сеймура Крэя ^[4] (Seymour Roger Cray) и имел этот суперкомпьютер производительность аж 3 MFLOPS. Их научная работа выглядела так:

То есть простым перебором на суперкомпьютере они нашли числа степени 5, которые опровергали гипотезу Эйлера: 27⁵ + 84⁵ + 110⁵ + 133⁵ = 144⁵.

И всё бы ничего, но другой умный товарищ спросил ^[5]: "интересно, а вот кто–нибудь из программистов может набросать код для супер–современного i5 по поиску еще таких совпадений?...".

Как вам уже понятно, предложение сработало как красная тряпка. Первое решение оказалось достаточно красивым и написанным с умом. Суть в том, что вначале считаем пятые степени для 1-N чисел, заносим их в таблицу и рекурсивно начинаем перебирать с самого низу четыре слагаемых пятых степеней, попутно ища в таблице сумму получившихся значений. Если нашли — вот оно и решение (индекс в таблице будет числом, степень которого мы нашли).

Вот этот первый вариант пользователя 2Alias ^[6]:

	#include <iostream>
	#include <algorithm>
	#include <stdlib.h>
	#include <vector>
	#include <unordered_map>
	 
	using namespace std;
	typedef long long LL;
	const int N = 250;
	 
	LL p5(int x)
	{
	    int t = x * x;
	    return t * (LL) (t * x);
	}
	 
	vector<LL> p;
	std::unordered_map<LL, int> all;
	int res[5];
	 
	void rec(int pr, LL sum, int n)
	{
	    if (n == 4)
	    {
	        if (all.find(sum) != all.end())
	        {
	            cout << "Okn";
	            res[4] = all[sum];
	            for (int i = 0; i < 5; ++i)
                cout << res[i] << " ";
	            cout << "n";
	            exit(0);
	        }
	        return;
	    }
	    for (int i = pr + 1; i < N; ++i)
	    {
	        res[n] = i;
	        rec(i, sum + p[i], n + 1);
	    }
	}
	 
	int main()
	{
	    for (int i = 0; i < N; ++i)
	    {
	        p.push_back(p5(i));
	        all[p.back()] = i;
	    }
	    rec(1, 0, 0);
	    return 0;
	}

И как оно обычно бывает, я подумал, а можно ли быстрее? Заодно у людей возник вопрос, а что будет если проверить в этом деле C#. Я в лоб переписал решение на C# и программа показала примерно такой же результат по времени. Уже интересно! Но оптимизировать всё же будем C++. Ведь потом всё легко перенести в C#.

Первое что пришло в голову — убрать рекурсию. Хорошо, просто введём 4 переменные и будем их перебирать с инкрементом старших при переполнении младших.

	// N - до какого числа ищем 1 - N в степени 5
	// powers - массив с посчитанными заранее степенями
	// all = unordered_map< key=степень числа, value=само число> 
	uint32 ind0 = 0x02; // ищем начиная с 2
	uint32 ind1 = 0x02;
	uint32 ind2 = 0x02;
	uint32 ind3 = 0x02;
	uint64 sum = 0;

	while (true)
	{
		sum = powers[ind0] + powers[ind1] + powers[ind2] + powers[ind3];
		if (all.find(sum) != all.end())
		{
			// нашли совпадение - ура!!
			foundVal = all[sum];
			...
		}

		// следующий
		++ind0;
		if (ind0 < N)
		{
			continue;
		}
		else
		{
			ind0 = 0x02;
			++ind1;
		}
		if (ind1 >= N)
		{
			ind1 = 0x02;
			++ind2;
		}
		if (ind2 >= N)
		{
			ind2 = 0x02;
			++ind3;
		}
		if (ind3 >= N)
		{
			break;
		}
	}

И тут же результат стал хуже. Ведь нам будут встречаться множество одинаковых сумм, большую часть которых рекурсивный алгоритм обходит. Допустим у нас числа от 1 до 3, переберём их:

111
112
113
121 < — уже было!
122
123
131 < — уже было!
132 < — уже было!
133
…

Значит придётся смотреть комбинаторику, но мне это не сильно помогло, так как достаточного объёма знаний по этой теме я не имею, пришлось выписывать примерно такой же пример на бумажку и подумать. И решение нашлось: при инкременте старших разрядов, не надо сбрасывать младшие на самый минимум, а присваивать значение старшего всем младшим — так мы отсечём лишнее.

Код увеличения индексов стал таким:

	...
	// следующий
	++ind0;
	if (ind0 < N)
	{
		continue;
	}
	else
	{
		ind0 = ++ind1;
	}
	if (ind1 >= N)
	{
		ind0 = ind1 = ++ind2;
	}
	if (ind2 >= N)
	{
		ind0 = ind1 = ind2 = ++ind3;
	}
	if (ind3 >= N)
	{
		break;
	}

Ура! И уже сразу чуть быстрее. Но что нам говорит профайлер? Большую часть времени мы сидим в unordered_map.find…

Начинаю вспоминать алгоритмы поиска и разнообразные знания(вплоть до демосцены). А что если перед проверкой в unordered_map как-то быстро отсекать часть ненужного?

Так появился ещё один массив, уже битовый (bitset). Так как числа нам в него не занести (они слишком большие), придётся быстро делать хэш из степени, приводить его к количеству бит в массиве и отмечать там. Всё это надо делать при построении таблицы степеней. В процессе написания оказалось, что std::bitset немного медленнее простого массива и минимальной логики что я набросал. Ну да ладно, это ерунда. А в целом ускорение оказалось существенным, около двух раз.

Много экспериментируя с размером bitset'a и сложностью хэша стало понятно, что по большому счёту влияет только память, причём для разных N по-разному и большая степень фильтрации обращений к unordered_map.find лучше только до определённого предела.

Выглядеть это стало примерно так:

	...
	// тут теперь мы быстро фильтруем сумму по хэшу и битовой таблице
	if (findBit(sum))
	{
		// и только потом проверяем в map, а проверять надо - ведь у нас коллизии из-за хэша
		if (all.find(sum) != all.end())
		{
			// нашли!
		}
	}
	// следующий
	...

Дальше возникла проблема номер два. Первый пример из далёкого 1966 года имел максимальное число 144⁵ (61 917 364 224), а второй (2004 год) уже 85359⁵ (4 531 548 087 264 753 520 490 799) — числа перестают влезать в 64 бита…

Идём самым простым путём: берём boost::multiprecision::uint128 — вот его нам хватит надолго! Это из-за того, что я пользуюсь MS CL, а он просто не поддерживает uint128, как все остальные компиляторы. Кстати, за время поиска решения проблемы uint128 и компиляторов я ещё узнал про шикарный сайт — Compiler Explorer ^[7]. Прямо в онлайне можно скомпилировать код всеми популярными компиляторами разных версий и посмотреть во что они транслируют его(ассемблер), причём с разными флагами компиляции. MS CL ^[8] тоже есть, но на бета сайте. И помимо С++ ^[9] есть Rust ^[10], D ^[11] и Go ^[12]. Собственно, по коду и стало понятно, что MS CL совсем не понимает 128 составные целые, все компиляторы транслируют перемножение двух 64 битных в 128 битную структуру за три умножения, а MS CL — за четыре. Но вернёмся к коду.

С boost::multiprecision::uint128 производительность упала в 25 раз. И это как-то совсем неправильно, ведь в теории должно быть не более 3х раз. Забавно, что на столько же упала производительность C# версии с типом decimal (он хоть и не совсем целочисленный, но его мантисса 96бит). А предварительная фильтрация обращений к Dictionary (своеобразный аналог unordered_map из STL) — работает хорошо, ускорение очень заметно.

Ну значит сами понимаете — стало обидно. Столько уже сделано и всё зря. Значит будем изобретать велосипед! То есть простой тип данных uint128. По сути, нам же нужно только присваивание, сравнение, умножение и сложение. Не так и сложно, но процесс в начале пошёл не туда, так как сначала я взялся за умножение и дошло это до использования ассемблера. Тут не чем гордиться, лучше всего себя показали интринсики ^[13]. Почему процесс пошёл не туда? А нам умножение-то и не важно. Ведь умножение участвует только на этапе просчёта таблицы степеней и в основном цикле не участвует. На всякий случай в исходниках остался файл с умножением на ассемблере — вдруг пригодится.

	friend uint128 operator*(const uint128& s, const uint128& d)
	{
		// intristic use
		uint64 h = 0;
		uint64 l = 0;
		uint64 h2 = 0;
		l = _mulx_u64(d.l, s.l, &h);
		h += _mulx_u64(d.l, s.h, &h2);
		h += _mulx_u64(d.h, s.l, &h2);
		return uint128( h, l);
	}

Со своим uint128 производительность тоже, конечно, просела, но всего на 30% и это отлично! Радости полно, но профайлер не забываем. А что если совсем убрать unordered_map и сделать из самописного bitset'a подобие map'a? То есть после вычисления хэша суммы мы можем уже использовать это число как индекс в ещё одной таблице(unordered_map тогда не нужен совсем).

	// вот так выглядит самописный map
	boost::container::vector<CompValue*> setMap[ SEARCHBITSETSIZE * 8 ];
	...

	// ComValue просто контейнер для степени и числа
	struct CompValue
	{
	...
		mainType fivePower;
		uint32 number;
	};

	// Поиск по битовому массиву и самописному map
	inline uint32 findBit(mainType fivePower)
	{
		uint32 bitval = (((uint32)((fivePower >> 32) ^ fivePower)));
		bitval = (((bitval >> 16) ^ bitval) & SEARCHBITSETSIZEMASK);
		uint32 b = 1 << (bitval & 0x1F);
		uint32 index = bitval >> 5;
		if((bitseta[index] & b) > 0)
		{
			for (auto itm : setMap[bitval])
			{
				if (itm->fivePower == fivePower)
				{
					return itm->number;
				}
			}
		}
		return 0;
	}

Так как проект совершенно несерьёзный и никакой полезной нагрузки не нёс, я сохранял все способы перебора, поиска и разных значений через жуткий набор дефайнов и mainType как раз один из них — это тип куда пишется степень числа, чтобы подменять его при смене только один раз в коде. Уже на этом этапе все тесты можно проводить с uint64, uint128 и boost::multiprecision::uint128 в зависимости от дефайнов — получается очень интересно.

И знаете, введение своего map'а — помогло! Но не на долго. Ведь понятно, что map не просто так придуман и используется везде, где только можно. Опыты — это, конечно, подтверждают. При определённом N (ближе к 1 000 000), когда все алгоритмы уже тормозят, голый map(без предварительного bitset'a) уже обходит самописный аналог и спасает только значительное увеличение битового массива и массива где хранятся наши значения степеней и чисел, а это огромное количество памяти. Примерный мультипликатор около N * 192, то есть для N = 1 000 000 нам нужно больше 200мб. А дальше ещё больше. К этому моменту ещё не пришло понимание, почему так падает скорость с ростом массива степеней, и я продолжил искать узкие места вместе с профайлером.

Пока происходило обдумывание, я сделал все испробованные способы переключаемыми. Ибо мало ли что.

Одна из последних оптимизаций оказалось простой, но действенной. Скорость C++ версии уже перевалила за 400 000 000 переборов в секунду для 64бит ( при N = 500 ), 300 000 000 переборов для 128 бит и всего 24 000 000 для 128 бит из boost, и влиять на скорость уже могло практически всё. Если перевести в Гб, то чтение получается около 20Гб в секунду. Ну может я где-то ошибся…

Вдруг стало понятно, что пересчитывать всю сумму на каждом проходе не надо и можно ввести промежуточную. Вместо трёх сложений будет одно. А пересчитывать промежуточную только при увеличении старших разрядов. Для больших типов это, конечно, сильнее заметно.

	mainType sum = 0U, baseSum = 0U;

	baseSum = powers[ind1] + powers[ind2] + powers[ind3];

	while(true)
	{
		sum = baseSum + powers[ind0];
		...

		// refresh without ind0
		baseSum = powers[ind1] + powers[ind2] + powers[ind3];
	}

Тут уже задача начинала надоедать, так как быстрее уже не получалось и я занялся C# версией. Всё перенёс туда. Нашёл уже готовый, написанный другим человеком UInt128 — примерно такой же простой, как и мой для C++. Ну и, конечно, скорость сильно подскочила. Разница оказалась меньше чем в два раза по сравнению с 64 битами. И это у меня ещё VS2013, то есть не roslyn (может он быстрее?).

А вот самописный map проигрывает по всем статьям Dictionary. Видимо проверки границ массивов дают о себе знать, ибо увеличение памяти ничего не даёт.

Дальше уже пошла совсем ерунда, даже была попытка оптимизировать сложение интринсиками, но чисто C++ версия оказалась самой быстрой. У меня почему-то не получилось заставить инлайниться ассемблерный код.

И всё же постоянно не отпускало чувство, что я что-то не вижу. Почему при росте массива всё начинает тормозить? При N = 1 000 000 производительность падает в 30 раз. Приходит в голову кэш процессора. Даже попробовал интринсик prefetch, результата — ноль. Пришла мысль кэшировать часть перебираемого массива, но при 1 000 000 значений (по 20 байт) как-то глупо выглядит. И тут начинает вырисовываться полная картина.

Так как числа у нас 4, есть 4 индекса, которые берут значения из таблицы. Таблица у нас постоянно возрастающих значений и сумма всех четырёх степеней у нас постоянно возрастающая (до переключения старших индексов). И разность степеней становится всё больше и больше.
2⁵ это 32, а 3⁵ это уже 243. А что если искать прямо в том же массиве просчитанных степеней обычным линейным поиском, сначала выставляя начальное значение на самый большой индекс и сохраняя индекс последнего найденного меньшего значения чем наша сумма (следующая будет больше) и использовать этот сохранённый индекс как начальную точку при следующем поиске, ведь значения не будут сильно меняться… Бинго!

Что в итоге?

	uint32 lastRangeIndex = 0;

	// линейный поиск в массиве степеней
	inline uint32 findInRange(mainType fivePower, uint32 startIndex)
	{
		while (startIndex < N)
		{
			lastRangeIndex = startIndex;
			if (powers[startIndex] > fivePower)
			{
				return 0;
			}
			if (powers[startIndex] == fivePower)
			{
				return startIndex;
			}
			++startIndex;
		}
		return 0;
	}

	...

	// главный цикл поиска
	baseSum = powers[ind1] + powers[ind2] + powers[ind3];
	while (true)
	{
		sum = baseSum + powers[ind0];

		foundVal = findInRange( sum, lastRangeIndex);
		if (foundVal > 0)
		{
			// нашли!
		}

		// следующий
		++ind0;
		if (ind0 < N)
		{
			continue;
		}
		else
		{		
			ind0 = ++ind1;
		}
		if (ind1 >= N)
		{
			ind0 = ind1 = ++ind2;
		}
		if (ind2 >= N)
		{
			ind0 = ind1 = ind2 = ++ind3;
		}
		if (ind3 >= N)
		{
			break;
		}
		// сброс индекса на начальное значение при изменении старших индексов
		lastRangeIndex = 0x02;
		if (ind1 > lastRangeIndex)
		{
			lastRangeIndex = ind1;
		}
		if (ind2 > lastRangeIndex)
		{
			lastRangeIndex = ind2;
		}
		if (ind3 > lastRangeIndex)
		{
			lastRangeIndex = ind3;
		}
		// refresh without ind0
		baseSum = powers[ind1] + powers[ind2] + powers[ind3];
	}

Скорость на маленьких значениях N немного уступает самописному map, но как только растёт N — скорость работы даже начинает расти! Ведь промежутки между степенями больших N растут чем дальше, тем больше и линейный поиск работает всё меньше! Сложность получается лучше O(1).

Вот вам и потеря времени. А всё почему? Не надо бросаться грудью на амбразуру, посиди — подумай. Как оказалось, самый быстрый алгоритм — это линейный поиск и никакие map/bitset не нужны. Но, безусловно, это очень интересный опыт.

Хабр любит исходники и они есть у меня ^[14]. В коммитах можете даже посмотреть историю «борьбы». Там лежат обе версии и C++, и C#, в котором этот трюк, конечно, так же отлично работает. Проекты хоть и вложены один в другой, но, конечно, никак не связаны.

Исходники ужасны, в начале находятся дефайны, где можно задать основное значение (uint64, uint128, boost::uin128/decimal(для C#), библиотеку можно переключать std/boost (boost::unordered_map оказался быстрее std::unordered_map примерно на 10%). Так же выбирается алгоритм поиска (правда сейчас вижу, что предварительный фильтр для unordered_map в версии C++ не пережил правок, но он есть в коммитах и в C# версии) unordered_map, самописный bitset и range (последний вариант).

Вот такая вот сказка и мне урок. А может и ещё кому будет интересно. Ведь много каких значений ещё не нашли ^[15]…

* к/ф Сказка о потерянном времени, 1964г.

Автор: crea7or

Источник ^[16]