Казка про втрачений час

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


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

Ну власне так воно і було до 1966 року…

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



Тобто простим перебором на суперкомп'ютері вони знайшли числа ступеня 5, які спростовували гіпотезу Ейлера: 275 + 845 + 1105 + 1335 = 1445.

І все б нічого, але інший розумний товариш запитав: "цікаво, а ось хто–небудь з програмістів може накидати код для супер–сучасного i5 з пошуку ще таких збігів?...".

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

Ось цей перший варіант користувача 2Alias:

#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 << "Ok\n";
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 року мав максимальне число 1445 (61 917 364 224), а другий (2004 рік) вже 853595 (4 531 548 087 264 753 520 490 799) — числа перестають влазити в 64 біта…

Йдемо найпростішим шляхом: беремо boost::multiprecision::uint128 — от його нам вистачить надовго! Це з-за того, що я користуюся MS CL, а він просто не підтримує uint128, як всі інші компілятори. До речі, за час пошуку рішення проблеми uint128 і компіляторів я ще дізнався про шикарний сайт — Compiler Explorer. Прямо в онлайні можна скомпілювати код усіма популярними компіляторами різних версій і подивитися на що вони транслюють його(асемблер), причому з різними прапорами компіляції. MS CL теж є, але на бета сайті. І крім С++ Rust, D і Go. Власне, за кодом і стало зрозуміло, що MS CL зовсім не розуміє 128 складові цілі, всі компілятори транслюють перемножування двох 64 бітних 128 бітну структуру за три множення, а MS CL — за чотири. Але повернемося до коду.

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

Ну значить самі розумієте — стало прикро. Стільки вже зроблено і все даремно. Значить будемо винаходити велосипед! Тобто простий тип даних uint128. По суті, нам же потрібно тільки присвоювання, порівняння, множення і додавання. Не так і складно, але процес на початку пішов не туди, так як спочатку я взявся за множення і дійшло це до використання асемблера. Тут не чим пишатися, найкраще себе показали интринсики. Чому процес пішов не туди? А нам множення і не важливо. Адже множення бере участь тільки на етапі прорахунку таблиці ступенів і в основному циклі не бере участь. На всяк випадок в исходниках залишився файл з множенням на асемблері — раптом знадобиться.

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 &h 2);
h += _mulx_u64(d.h, s.l &h 2);
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 переборів в секунду для 64bit ( при 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 індексу, які беруть значення з таблиці. Таблиця у нас постійно зростаючих значень і сума всіх чотирьох ступенів у нас постійно зростаюча (до перемикання старших індексів). І різниця ступенів стає все більше і більше.
25 32, а 35 це вже 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 не потрібні. Але, безумовно, це дуже цікавий досвід.

Хабр любить исходники і вони є у мене. У коммитах можете навіть подивитися історію «боротьби». Там лежать обидві версії і 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 (останній варіант).

Ось така казка і мені урок. А може і ще комусь буде цікаво. Адже багато яких значень ще не знайшли


* к/ф " Казка про втрачений час, 1964р.
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.