Доброго времени суток, господа и дамы! Иногда у некоторых людей возникает желание заняться откровенным непотребством в программировании — то, что не несет практической пользы напрямую, но помогает развлечься. И я — не исключение. В этой статье я хочу рассказать вам о лайфхаках, трюках (магических и не очень), алгоритмах на языке C!
Идея написать эту статью зародилась из моего поста. В нем я рассказал о том, что через последовательность Фибоначчи можно конвертировать мили в километры с небольшой погрешностью. Увидев, что многим понравилась, я задумался: почему бы не изучить еще какие-нибудь трюки, заодно практикуясь в программировании на C?
Всех, кто заинтересовался — прошу под кат.
Раз я уже начал рассказывать о Фибоначчи, немного процитирую свой пост:
Последовательность Фибоначчи может конвертировать мили в километры с небольшой погрешностью. 5 миль ≈ 8 км (5 и 8 — числа Фибоначчи). Реальность: 5 миль = 8.04672 км. Почему? 1 миля = 1.609344 километра (точное значение). Золотое сечение (φ) ≈ 1.618034. Погрешность возникает потому что отношение Fₙ₊₁ / Fₙ стремится к φ ≈ 1.618034, а точное соотношение миля/км = 1.609344. Относительная погрешность: (1.618034 — 1.609344) / 1.609344 * 100% ≈ 0.54%.
Фанфакт: если взять куб и шар одинакового объёма, то радиус шара в милях будет равен ребру куба в километрах. Потому что кубический корень из (4/3)pi как раз равен 1.61. @Enfriz
Давайте посмотрим, как реализовать такой конвертер на C.
Вот допустим есть такая функция:
uint64_t fibonacci(int num) {
if (num < 0)
return 0;
if (num == 0)
return 0;
uint64_t a = 0;
uint64_t b = 1;
if (num == 1)
return b;
for (int i = 2; i <= num; i++) {
uint64_t next = a + b;
a = b;
b = next;
}
return b;
}
Я думаю вы все знаете как работает данный алгоритм. Пояснить только стоит что мы должны передать как аргумент инкрементированное на 1 количество милей — то есть чтобы конвертировать 5 милей в километры, мы передаем 5+1=6 (следующее число) и получаем 8.
Также, можно сделать вычисление на основе золотого сечения — ибо оно также связано с последовательностью Фибоначчи.
float fib_golden_ratio(float miles) {
const double PHI = (1.0 + sqrt(5.0)) / 2.0;
if (miles < 1e-5) {
return 0.0F;
}
double n = log(miles * sqrt(5.0)) / log(PHI);
int k = (int)floor(n);
double Fk = (pow(PHI, k) - pow(-PHI, -k)) / sqrt(5.0);
double Fk1 = (pow(PHI, k + 1) - pow(-PHI, -k - 1)) / sqrt(5.0);
double Fk2 = (pow(PHI, k + 2) - pow(-PHI, -k - 2)) / sqrt(5.0);
if (Fk1 - Fk < DBL_EPSILON) {
return basic_miles2km(miles);
}
return Fk1 + ((miles - Fk) * ((float)(Fk2 - Fk1) / (Fk1 - Fk)));
}
На основе числа φ ((1 + √5) / 2) мы также можем конвертировать мили в километры. Золотое сечение φ ≈ 1.618 близко к 1.609, поэтому числа Фибоначчи имитируют конвертацию
Переменная n здесь — это индекс числа Фибоначчи, соответствующий значению miles, основанный на свойстве, что F(n) ≈ φ^n / √5. Используется логарифмическая функция, чтобы найти n, а затем округляется до ближайшего меньшего целого числа k.
Переменные Fk, Fk1, Fk2 — это три последовательных числа Фибоначчи, вычисляемые через формулы Бине для чисел Фибоначчи: F(n) = (φ^n — (1 — φ)^n) / √5.
Формула Бине — это явная формула для нахождения n-го числа Фибоначчи без рекурсивного вычисления. Она названа в честь французского математика, который её открыл. Формула Бине позволяет вычислять числа Фибоначчи за константное время (O(1)), тогда как рекурсивный способ может занять экспоненциальное время. Это делает формулу Бине очень полезной в вычислениях, связанных с последовательностями Фибоначчи.
Если разница между F(k+1) и F(k) меньше DBL_EPSILON, это означает что числа слишком близки друг к другу, и функция вместо вычисления дальнейшего значения просто вызывает функцию basic_miles2km.
Но почему мы вызываем в basic_miles2km? Мы же с вами затрагиваем нетипичные способы конвертации. На самом деле это — защитный механизм. Если Fk1 ≈ Fk (разница меньше DBL_EPSILON), знаменатель стремится к нулю → возникает неопределённость или NaN. Также формула Бине может некорректно работать при малых n (особенно n < 3).
Условие (Fk1 - Fk < DBL_EPSILON) означает: «Разница между соседними числами Фибоначчи меньше порога различимости чисел с плавающей точкой». DBL_EPSILON — это машинный эпсилон для чисел типа double (64 бита) из float.h, равна 2−52 ≈ 2,20e-16. Константа DBL_EPSILON используется для сравнения чисел с плавающей точкой. Например, два числа являются одинаковыми с точки зрения машинной арифметики, если их относительная разность по модулю меньше DBL_EPSILON.
float basic_miles2km(float miles) {
return miles * 1.609344f;
}
И в конце идет интерполяция для нахождения результативного числа Фибоначчи: return Fk1 + ((miles - Fk) * ((float)(Fk2 - Fk1) / (Fk1 - Fk)));. Если предыдущие условия не выполнены, функция выполняет линейную интерполяцию, чтобы получить более точное значение числа Фибоначчи, основываясь на значении miles. Но эта функция при больших n может накапливать ошибки.
Перейдем к другой похожей реализации:
float fib_interpolate(float miles) {
if (miles < 5.0F) {
return basic_miles2km(miles);
}
uint64_t prev_mile = 0;
uint64_t prev_km = 1;
uint64_t curr_mile = 1;
uint64_t curr_km = 2;
while (curr_mile <= miles) {
prev_mile = curr_mile;
prev_km = curr_km;
curr_mile = prev_km;
curr_km = prev_mile + prev_km;
if (curr_km < prev_km || curr_mile < prev_mile) {
break;
}
}
return prev_km + ((miles - prev_mile) * ((float)(curr_km - prev_km) / (curr_mile - prev_mile)));
}
В функции объявлены переменные, использующие 64-битные целые числа без знака для хранения значений миль и километров, где prev_mile и prev_km хранят предыдущие расстояния, а curr_mile и curr_km — текущие. Затем начинается основной цикл, который будет выполняться до тех пор, пока текущее значение миль (curr_mile) меньше или равно переданному значению miles. На каждой итерации цикла обновляются значения переменных через стандартные формулы Фибоначчи, что соответствует генерации последовательности (при этом текущее значение миль обновляется на значение предшествующего километра, а километры берут значение суммы двух предыдущих).
Чтобы предотвратить переполнение переменных, в теле цикла также содержится проверка: если новое значение километров или миль становится меньше предыдущего, цикл прерывается. Это условие служит защитой от некорректных данных и потенциальных ошибок исполнения.
После завершения цикла производится интерполяция для получения значения километров на основе предшествующих значений, где используется линейная интерполяция. Этот расчет позволяет более точно находить значение, чем простое преобразование, основанное только на числах Фибоначчи. Тем не менее, такая интерполяция может быть недостаточно точной для больших значений миль, что ослабляет точность результата в зависимости от диапазона входных параметров.
Временная сложность функции fib_interpolate — O(n), где n — это позиция в последовательности Фибоначчи.
Более подробно весь код можете просмотреть в репозитории fib_miles2km или прямо в репозитории-сборнике для этой статьи The Art Of Fun C.
❯ Бинарное возведение в степень
Раз уж мы использовали pow() в формуле Бине, давайте оптимизируем его через бинарное возведение.
Давайте попробуем импортозаместить стандартный pow нашим алгоритмом, а точнее — бинарным возведением в степень. Сложность — O(log n).
Широко известный алгоритм для возведения любого числа в целую степень с абсолютной точностью. Принцип действия прост: есть целая степень e, чтобы получить число b в этой степени нужно возвести это число во все степени 1, 2, 4, … 2n (в коде этому соответствует b *= b), каждый раз сдвигая биты e вправо (e >>= 1) пока оно не равно 0 и тогда, когда последний бит e не равен нулю ((e & 1) != 0), домножать результат v на полученное b. Этот метод позволяет возвести число в целую степень за логарифмическое время, что особенно важно при работе с большими числами.
Возьмём 2⁵. Вместо 2×2×2×2×2 мы представляем 5 в двоичном виде (101) и идем по битам справа налево. Первый бит (1) равен 2, второй бит (0) равен 2 в квадрате, и в последнем третьем биту (1) результат (2) умножаем на текущее значение (4²=16) → 32.
Технически, при каждом проходе цикла, мы проверяем младший бит степени e. Если этот бит равен 1 (e & 1 != 0), мы умножаем текущий результат v на основание b. После этого мы умножаем b на само себя (квадрат), что соответствует следующей степени двоичного представления. Уменьшая e с помощью побитового сдвига вправо, мы продолжаем цикл, пока e не достигает нуля, в результате чего получаем финальный результат.
double binary_pow(double b, unsigned long long e) {
double v = 1.0;
while(e != 0) {
if((e & 1) != 0) {
v *= b;
}
b *= b;
e >>= 1;
}
return v;
}
И при выполнении мы можем получить 10.00 ** 2.00 = 100.00, 10.50 ** 2.00 = 110.25. Но если экспонента не целое число, то данный способ не подойдет (10.50^2.50 вычисляется как 110.25). Так что важно отметить, что этот алгоритм применим только к целым числам.
Ну и вот пример измененной функции fib_golden_ratio из прошлого примера:
float fib_golden_ratio_binary(float miles) {
const double PHI = (1.0 + sqrt(5.0)) / 2.0;
if (miles < 1e-5) {
return 0.0F;
}
double n = log(miles * sqrt(5.0)) / log(PHI);
int k = (int)floor(n);
double sign_k = (k % 2 == 0) ? 1.0 : -1.0;
double sign_k1 = ((k+1) % 2 == 0) ? 1.0 : -1.0;
double sign_k2 = ((k+2) % 2 == 0) ? 1.0 : -1.0;
double phi_k = binary_pow(PHI, k);
double phi_k1 = binary_pow(PHI, k+1);
double phi_k2 = binary_pow(PHI, k+2);
double Fk = (phi_k - sign_k / phi_k) / sqrt(5.0);
double Fk1 = (phi_k1 - sign_k1 / phi_k1) / sqrt(5.0);
double Fk2 = (phi_k2 - sign_k2 / phi_k2) / sqrt(5.0);
if (Fk1 - Fk < DBL_EPSILON) {
return basic_miles2km(miles);
}
return Fk1 + ((miles - Fk) * ((float)(Fk2 - Fk1) / (Fk1 - Fk)));
}
❯ Генератор псевдослучайных чисел Xorshift
Давайте перейдем к не менее интересному алгоритму — генератору псевдослучайных чисел. Давайте попробуем реализовать свой генератор таких чисел на основе Xorshift. Xorshift — это семейство алгоритмов генерации псевдослучайных чисел, которые используют операцию сдвига и исключающего ИЛИ (XOR). Они просты в реализации, очень быстры и дают достаточно хорошие результаты для многих задач.
uint64_t xorshift64(uint64_t *state) {
uint64_t x = *state;
x ^= x << 13;
x ^= x >> 7;
x ^= x << 17;
*state = x;
return x;
}
Данный алгоритм работает так:
Сдвиг влево + XOR (x << 13): «Размазывает» биты в старшие разряды. Операция XOR смешивает оригинальные биты с изменёнными.
Сдвиг вправо + XOR (x >> 7): Работает как «обратный» шаг — компенсирует предыдущий сдвиг, добавляя нелинейность.
Сдвиг влево + XOR (x << 17): Фиксирует результат и обеспечивает полное перемешивание битов. Числа 13, 7, 17 подобраны экспериментально для оптимального качества.
Ключевое преимущество — всего 3 операции на число, что делает алгоритм невероятно быстрым — O(1) на генерацию числа.
Для генерации state давайте напишем функцию для получения таймстемпа в микросекундах:
uint64_t get_seed() {
struct timeval tv;
gettimeofday(&tv,NULL);
return tv.tv_sec*(uint64_t)1000000+tv.tv_usec;
}
Кроме того, давайте для полноты реализуем функции для генерации плавающего числа и генерацию по диапазону:
double rand_double(uint64_t *state) {
return (xorshift64(state) >> 11) * (1.0 / (UINT64_C(1) << 53));
}
Данная функция возвратит число в диапазоне [0.0, 1.0). Сначала генерируем 64-битное число и сдвигаем на 11 бит вправо. Это нужно чтобы отбросить младшие биты. В итоге остается 53 бита. После вычисляем 1.0 / (UINT64_C(1) << 53) (UINT64_C(1) << 53 = 2⁵³ = 9,007,199,254,740,992.). В итоге получаем минимальный шаг между числами double: 1.0 / 2⁵³ ≈ 1.110223e-16. И под конец перемножаем 53 бита на этот самый шаг.
Если вы задались вопросом, почему именно 53 бита — то в стандарте IEEE 754 тип double имеет 53 бита мантиссы. Структура числа двойной точности: знак — 1 бит, порядок — 11 бит, мантисса — 52+1 бит.
И генерация случайного числа в диапазоне:
uint64_t rand_range(uint64_t *state, uint64_t min, uint64_t max) {
return min + xorshift64(state) % (max - min + 1);
}
Он будет уже попроще: xorshift64(state) даёт случайное число. Через % (max — min + 1) узнаем остаток от деления на длину диапазона, и прибавляем к итоговому числу min.
И в итоге мы можем получить такие результаты работы:
xorshift64 random number: 8549869788877919663 # Случайное число
xorshift64 random num from 10 to 100: 72 # Случайное число в диапазоне
xorshift64 double random number: 0.461347 # Случайное дробное число
❯ Генератор случайных чисел xoshiro256pp
Давайте перейдем к более сложному, но интересному генератору псевдослучайных чисел. В отличии от xorshift, этот генератор уже может проходить статистические тесты.
Давайте разберем код:
typedef struct {
uint64_t s[4];
} xoshiro256pp_state;
static inline uint64_t rotl(const uint64_t x, int k) {
return (x << k) | (x >> (64 - k));
}
uint64_t xoshiro256pp_next(xoshiro256pp_state *state) {
uint64_t *s = state->s;
uint64_t result = rotl(s[0] + s[3], 23) + s[0];
uint64_t t = s[1] << 17;
s[2] ^= s[0];
s[3] ^= s[1];
s[1] ^= s[2];
s[0] ^= s[3];
s[2] ^= t;
s[3] = rotl(s[3], 45);
return result;
}
void xoshiro256pp_init(xoshiro256pp_state *state, uint64_t seed) {
uint64_t tmp = seed;
for (int i = 0; i < 4; i++) {
tmp ^= tmp >> 30;
tmp *= 0xbf58476d1ce4e5b9;
tmp ^= tmp >> 27;
tmp *= 0x94d049bb133111eb;
tmp ^= tmp >> 31;
state->s[i] = tmp;
}
}
Создаем структуру
xoshiro256pp_state, которая содержит массив s из четырех 64-битных целых чисел (всего 256 бит).Функция
rotlвыполняет циклический сдвиг влево для 64-битного числа x на k позиций. Это означает, что биты, которые выходят за пределы числового представления, возвращаются в начало.Функция
xoshiro256pp_initиспользуется для инициализации состояния генератора на основе сида (можно также использовать функцию get_seed из прошлого примера). Используется цикл, чтобы заполнить массив s четырьмя значениями.
tmp ^= tmp >> 30— применяется побитовый сдвиг вправо, чтобы смешать верхние биты в нижние.tmp *= 0xbf58476d1ce4e5b9иtmp *= 0x94d049bb133111eb— умножение на числа, также подобранные экспериментальным путем для наибольшей точности.Каждое из этих преобразований применяется несколько раз, обеспечивая высокую энтропию.
В функции
xoshiro256pp_nextпередается указатель на структуру состояния. Сначала мы создаем указатель s на массив состояния. Затем производим расчет result — это основное сгенерированное значение, которое относится к элементам состояния. Сначала складываем s[0] и s[3], затем поворачиваем результат на 23 бита влево и добавляем s[0]. В основном теле функции происходит несколько перестановок и XOR-операций. Хранится значение t (s[1] сдвинутый на 17 бит влево) и затем происходит несколько операций XOR для смешивания элементов массива, чтобы следующие состояние завесило от предыдущих элементов. Затем, мы вызываемrotlдля поворота на 45 бит и возвращаем результат.
Данный алгоритм лучше xorshift64 благодаря тому использованию 256 бит состояния вместо 64. Также он генерирует числа с лучшим равномерным распределением. Сама случайность более качественная и непредсказуемая, что делает данный алгоритм весьма хорошим (для криптографических дел, естественно, не подойдет). Все это достигается через более сложную логику генерации, что делает предсказание чисел более проблематичным. Работает он чуть медленее xorshift64, но результат работы намного лучше.
❯ Генератор случайных чисел Лемера
Очень быстрый генератор псевдослучайных чисел. Подробнее можно почитать на странице википедии.
Алгоритм Lehmer64 — это быстрый генератор псевдослучайных чисел (ГПСЧ), основанный на мультипликативном линейном конгруэнтном методе. Он использует 128-битное состояние для генерации 64-битных случайных чисел, обеспечивая высокую производительность и достаточную статистическую случайность для большинства некриптографических задач.
Мультипликативный линейный конгруэнтный метод (мультипликативный конгруэнтный метод) — это алгоритм для генерации псевдослучайных чисел, основанный на линейной рекуррентной формуле. Простыми словами, метод создаёт последовательность чисел, где каждое последующее число зависит от предыдущего, но при этом существует цикл, повторяющийся бесконечное число раз (период).
Принцип работы
Основная формула метода: Xn+1 = (a ⋅ Xn + c) mod m, где:
Xn+1 — следующее псевдослучайное число;
Xn — текущее псевдослучайное число;
a — мультипликативный множитель;
c — приращение;
m — модуль;
X0 — начальное значение (seed).
И вот его простейшая реализация на C:
static __uint128_t g_lehmer64_state;
uint64_t lehmer64(void) {
g_lehmer64_state *= 0xda942042e4dd58b5ULL;
return g_lehmer64_state >> 64;
}
Каждый раз, когда нужно новое число, мы умножаем состояние на константу 0xda942042e4dd58b5ULL (высчитанная экспериментальным путем). g_lehmer64_state — это состояние размером 128 бит, на старте в котором хранится нечетный сид. После этого берем старшие 64 бита и возвращаем их, младшие остаются.
Данный генератор прост в реализации, но по моим бенчмаркам он является самым медленным среди трех генераторов из статьи.
❯ Быстрый обратный квадратный корень из Quake III
В 2005 году id Software опубликовала под лицензией GPL-2 исходный код своей игры 1999 года Quake III Arena. В файле code/game/q_math.c есть функция для вычисления обратного квадратного корня числа, которая на первый взгляд выглядит очень любопытным алгоритмом:
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // what the fuck?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}
Как же работает этот алгоритм? Он выполняется в два этапа:
Получение грубой аппроксимации y обратного квадратного корня нужного числа number:
y = number;
i = * ( long * ) &y;
i = 0x5f3759df - ( i >> 1 );
y = * ( float * ) &i;
Уточнение аппроксимации при помощи одного шага метода Ньютона-Рафсона (NR):
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = y * ( threehalfs - ( x2 * y * y ) );
Алгоритм принимает 32-битное число с плавающей запятой (одинарной точности в формате IEEE 754) в качестве исходных данных. Точность алгоритма — менее 0,2% в меньшую сторону и никогда — в большую. Это не хватает для настоящих численных расчётов, но достаточно для трёхмерной графики.
Подробнее о работе данной функции можно прочитать в этой статье и на странице в википедии.
Для проверки можно немного изменить код:
float Q_rsqrt(float number) {
int32_t i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( int32_t* ) &y;
i = 0x5f3759df - ( i >> 1 );
y = * ( float* ) &i;
y = y * ( threehalfs - ( x2 * y * y ) );
return y;
}
И мы должны получить почти-что правильный результат: Q_rsqrt(25.00) = 0.199690. Для точности можно добавить еще одну итерацию Ньютона y = y * (threehalfs — (x2 * y * y));: Q_rsqrt(25.00) = 0.199999.
❯ Заключение
Спасибо за прочтение статьи! Я надеюсь вы узнали что-то новенькое, или может какой-нибудь трюк натолкнул вас на другой интересный алгоритм. Если нашли нюанс в самой статье — пишите в комментарии. Возможно я продолжу серию таких статей на тему трюков на разных языках программирования.
Если вам понравился изложенный материал, могу предложить вам подписаться на мой блог в телеграме. Если, конечно, вам статья понравилась и вы хотите видеть чуть больше.
Примеры работы кода вы можете увидеть в моем репозитории The Art Of Fun C.
Все представленные методы — скорее интеллектуальное развлечение, чем практические решения. Но они отлично демонстрируют, как математика и низкоуровневые трюки могут создавать неочевидные оптимизации. Хотя вполне себе можно использовать генераторы псевдослучайных чисел для симуляций или для геймдева или даже обратный корень из самого Quake III.
❯ Бенчмарки
В репозитории я также написал небольшие бенчмарки:
PRNG Performance (10,000,000 iterations):
-----------------------------------------
xorshift64: 14.18 ms (705.37M numbers/s)
lehmer64: 20.71 ms (482.89M numbers/s)
xoshiro256pp: 15.95 ms (626.77M numbers/s)
-----------------------------------------
Conversion Methods Performance (each method called 10000 times per point):
----------------------------------------------------------------------
Basic : 0.28 ms ( 0.001 us/call)
Fibonacci Interpolation : 1.56 ms ( 0.008 us/call)
Fibonacci Cache : 1.41 ms ( 0.007 us/call)
Golden Ratio : 16.33 ms ( 0.082 us/call)
Golden Ratio (Binary) : 3.56 ms ( 0.018 us/call)
----------------------------------------------------------------------
Accuracy Comparison (5 sample points):
Miles | Basic | INTERPOL | Cache | Golden | GoldenBin
------+-----------+----------+----------+----------+-----------
5 | 8.05 | 0.58% | 0.58% | 0.58% | 0.58%
30 | 48.28 | 0.53% | 0.53% | 0.53% | 0.53%
55 | 88.51 | 0.55% | 0.55% | 0.55% | 0.55%
80 | 128.75 | 0.54% | 0.54% | 0.54% | 0.54%
100 | 160.93 | 0.54% | 0.54% | 0.54% | 0.54%
---------------------------------------------------------------
Итог:
Генераторы
xorshift64 — самый быстрый, но менее качественный чем xoshiro256pp.
xoshiro256pp — менее быстрый чем xorshift64, но более качественный.
lehmer64 — самый медленный.
Конвертация (результаты для 200,000 вызовов (20 точек × 10,000 итераций))
Basic метод (простое умножение) — самый быстрый. Но здесь мы не за скоростью шли, а за интересным методом.
Фибоначчи-методы (Interpolation и Cache) — хорошая производительность.
Golden Ratio методы — значительно медленнее. Но binary-версия быстрее обычной.
У всех методов погрешность от 0.53 до 0.58 процентов.
❯ Источники
Новости, обзоры продуктов и конкурсы от команды Timeweb.Cloud — в нашем Telegram-канале ↩
Перед оплатой в разделе «Бонусы и промокоды» в панели управления активируйте промокод и получите кэшбэк на баланс.
AndreyDmitriev
Можно, наверное, добавить в копику быстрое вычисление приближённого значения степени:
VBDUnit
Или его улучшенный вариант :)