kinh 3 мая в 14:20

Быстрое приближённое умножение и деление чисел с плавающей точкой

Сложный

27 мин

6.6K

C++*Алгоритмы*Математика*Программирование микроконтроллеров*Разработка игр*

Туториал

Содержание

Введение
Умножение положительных чисел
Умножение на отрицательное число
Антипереполнение (Underflow)
Переполнение (Overflow)
Умножение на ноль и денормализованные числа
Умножение на бесконечность, неопределённость и нечисла
Умножение без операторов сравнения и переходов
Деление
Приближённое вычисление обратного значения
Антипереполнение при делении
Переполнение при делении
Специальные значения делимого и делителя
Деление без операторов сравнения и переходов
Итоговая программа
Повышение точности умножения и деления с помощью таблиц поправок

Введение

Особенности формата хранения чисел с плавающей точкой позволяют быстро находить приближённое значение логарифма, и, за счёт этого, выполнять умножение и деление. Результат при этом будет неточным, однако может быть применимым там, где особая точность не требуется.

В статье будет использоваться формат хранения чисел float языка C++, который предполагается одинаковым с форматом binary32 стандарта IEEE 754 - 2008 (IEEE Standard for Floating-Point Arithmetic).
Для краткости, в данной статье то, что хранится в поле порядка числа с плавающей точкой, будет называться «порядком» (Exponent), а то, что хранится в поле мантиссы – «мантиссой» (Mantissa). Следует заметить, что эти значения – не то, что в классической математике понимается под мантиссой и порядком. Но так как в этой статье эти термины ни для чего другого использоваться не будут, то неоднозначностей не должно возникнуть.

Умножение

Умножение положительных чисел

В формате float число хранится в 32 битах, при этом младшие 23 бита хранят мантиссу M, к ним примыкают 8 битов, которые хранят порядок E, и самым старшим битом является бит знака S. Если бит знака равен 1, то число – отрицательное. Пока мы будем рассматривать лишь положительные числа, для которых бит знака = 0. Значения порядка, которые полностью состоят из нулевых битов, или полностью из единичных, являются специальными, и мы их тоже рассмотрим позднее.

	Знак	Порядок	Мантисса
Обозначение	S	E	M
Размер в битах	1	8	23

Мантиссу удобно для наших целей рассматривать как число с фиксированной точкой, причём эта точка расположена сразу слева от мантиссы. Таким образом, если самый старший бит мантиссы хранит 1, то мантисса M = 0.5
Представим некое положительное вещественное число r > 0 в следующем виде:

$r = m 2^{P}$

При этом, пусть $P = \left \lfloor \log_2 r \right \rfloor$
Нижними полуквадратными скобками $\lfloor x \rfloor$ обозначается округление числа x вниз до ближайшего целого. Для любого вещественного x справедливо тождество:
$x = \left \lfloor x \right \rfloor + \{x\}$ , где дробная часть числа: $0 ≤ \{ x \} < 1$

Рассчитаем m:

$m = r 2^{-P}=2^{\log_2 r - P}=2^{\log_2 r - \left \lfloor \log_2 r \right \rfloor }=2^{ \left \{ \log_2 r \right \} }$ $1 \leqslant m < 2$

Если мы храним число r в формате float, то справедливы следующие соотношения:
Мантисса M = m - 1
Порядок E = P + BIAS, где BIAS – специальная константа, которая зависит от формата. Для формата float: BIAS = 0x7F = 127
Выразим r через мантиссу и порядок:

$r=(1+M) 2^{E-BIAS}$

При этом:

$E = \left \lfloor \log_2 r \right \rfloor + BIAS$

Рассмотрим два вещественных числа: x > 0 и y > 0. Выразим их через мантиссы и порядки:

$x=(1+M_x) 2^{E_x-BIAS}$ $y=(1+M_y) 2^{E_y-BIAS}$

Как известно из алгебры, для произведения положительных чисел справедливо:

$x y=2^{\log_2 (x y)}=2^{\log_2 x +\log_2 y}$

Подставим сюда мантиссы и порядки:

$x y =2^{\log_2 \left ((1+M_x) 2^{E_x - BIAS} \right ) +\log_2 \left ((1+M_y) 2^{E_y-BIAS} \right )} = 2^{\log_2 (1+M_x) +\log_2 (1+M_y) + E_x + E_y - 2 BIAS}$

Так как диапазон изменения мантиссы в нашем формате $0 \leqslant M < 1$ , то найдём значения функции $\log_2 (1+M)$ на обоих концах этого диапазона:

$\log_2 (1+M) \vert_{M = 0} = 0$ $\log_2 (1+M) \vert_{M = 1} = 1$

То есть мы видим, что мантисса совпадает со значениями $\log_2 (1+M)$ на обоих концах своего диапазона. Следовательно, мантисса является линейной интерполяцией этой функции на данном диапазоне:

$\log_2 (1+M) \approx {M}$

Оценим точность этого представления. Для начала найдём значение M, при котором будет наибольшее отклонение. Для этого решим уравнение (log_2(1+M) - M)’=0 . Решив его, получаем, что экстремум достигается при $M= (1/{\ln 2}) – 1 \approx 0.443$ . В этой точке абсолютное отклонение приблизительно равно 0.086, что не так уж и много.
Подставим нашу интерполяцию в формулу для произведения:

$x y \approx 2^{M_x + M_y + E_x + E_y - 2 BIAS}$

Обозначим r = x y . Порядок этого произведения будет равен:

$E = \left \lfloor \log_2 r \right \rfloor + BIAS \approx \left \lfloor M_x + M_y + E_x + E_y - 2 BIAS \right \rfloor + BIAS = \left \lfloor M_x + M_y \right \rfloor + E_x + E_y - BIAS$

Так как мантисса принимает значения $0 \leqslant M < 1$ , то сумма двух мантисс содержится в интервале 0 ≤ M_x + M_y < 2 . Рассмотрим два случая.
Случай 1: M_x + M_y < 1
Округление вниз суммы обнулится: $\left \lfloor M_x + M_y \right \rfloor = 0$
Тогда порядок произведения: $E = E_x + E_y – {BIAS}$
Случай 2: 1 ≤ M_x + M_y < 2
Округление вниз суммы: $\left \lfloor M_x + M_y \right \rfloor = 1$
Тогда порядок произведения: E = E_x + E_y – BIAS + 1
И здесь нужно учесть довольно удачное расположение полей мантиссы и порядка в формате float. Если мы сложим два таких числа, побитово интерпретируя их как целые числа, то если при сложении мантисс результат будет больше 1, то автоматически возникнет перенос единицы в поле порядка результата. То есть, нам, по сути дела, ничего специально делать не нужно: оба случая суммы мантисс обрабатываются одинаковым выражением.
Теперь найдём мантиссу произведения. Как мы уже установили:

$M=m-1=2^{ \left \{ \log_2 r \right \} }-1$

Подставим сюда нашу аппроксимацию:

$r = x y \approx 2^{M_x + M_y + E_x + E_y - 2 BIAS}$

Тогда:

$M≈2^{ \left \{ M_x + M_y + E_x + E_y - 2 BIAS \right \} }-1$

Так как E_x, E_y и BIAS представляют собой целые числа, то их можно выбросить из взятия дробной части. Остаётся:

$M≈2^{ \left \{ M_x + M_y \right \} }-1$

Для дробной части справедливо неравенство: $0 ≤ \left \{ M_x + M_y \right \} < 1$ . Следовательно, мы можем применить нашу аппроксимацию:

$\left \{ M_x + M_y \right \} \approx \log_2 (1+ \left \{ M_x + M_y \right \})$

Следовательно, мантисса произведения:

$M \approx \left \{ M_x + M_y \right \}$

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

typedef signed __int32 i32;
typedef unsigned __int32 u32;

/*
	Число одинарной точности
	Размер: 32 бита
	Мантисса: 23 бита
	Порядок: 8 битов
*/

union Float
{
	i32 i;
	u32 u;
	float f;
};

#define BITS_PER_ALL		32
#define BITS_PER_MANTISSA	23

#define BIAS		0x7F
#define EXP_BIAS	(BIAS << BITS_PER_MANTISSA)
#define EXP_ONE		(1 << BITS_PER_MANTISSA)
#define EXP_MASK	(0xFF << BITS_PER_MANTISSA)

#define SIGN_MASK	(1 << (BITS_PER_ALL - 1))

// Приближённое умножение
float apprMul(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;
	result.u = x.u + y.u - EXP_BIAS;
	return result.f;
}

Умножение на отрицательное число

Теперь рассмотрим сомножители, которые могут быть отрицательными. Бит знака результата при этом вычисляется через операцию XOR битов знака сомножителей. Соответственно, если мы выполним XOR сомножителей целиком, как если бы они были целыми числами, то старший бит результата этой операции будет равен знаку результата произведения.
Теперь заметим, что если мы просто просуммируем два бита, то младший бит результата тоже вычисляется как результат операции XOR слагаемых.
Справа от поля знака находится поле порядка. Мы сначала к полю порядка прибавляем какое-то число, затем ещё одно число отнимаем. Если после этих двух операций не возникло ни переноса, ни заёма, то бит знака результата произведения можно вычислить как сумму битов знаков сомножителей. То есть, в этом случае алгоритм приближённого умножения, который мы уже вывели для положительных сомножителей, можно применить и для отрицательных, без каких-либо изменений.

Антипереполнение (Underflow)

Как мы только что упоминали, при операциях с полем порядка мы можем выйти за допустимый диапазон. Начнём с антипереполнения. Эта ситуация может возникнуть, когда при вычитании из суммы порядков константы BIAS результат становится меньше, либо равен нулю. Равенство нулю тоже относится к антипереполнению, так как нулевое значение поля порядка зарезервировано для специальных значений.
В случае антипереполнения мы будем обнулять результат. Денормализованные числа мы использовать не будем, так как у нас всё равно умножение неточное. Обычно, большие по модулю отрицательные степени, которые и вызывают антипереполнение, возникают в результате ошибок округления, когда вычисляется разность двух очень близких чисел. Так что обнуление результата в таких случаях вполне оправдано: скорее всего, результатом и должен быть ноль.
Признаком антипереполнения является следующая ситуация: после суммирования порядков переполнения не было (то есть, бит знака не повреждён), однако после вычитания константы BIAS произошёл заём из поля знака, либо поле порядка обнулилось.
Бит знака результата формируется в результате сложения. Для того, чтобы знак временно отключить, мы альтернативно найдём его через XOR сомножителей, и проксорим его с суммой сомножителей.
Чтобы вместо двух проверок порядка результата умножения сделать только одну, мы вычтем из порядка единицу. Это изменит заём лишь в двух случаях: когда поле порядка было равно 0 – тогда заём появится, что нам и нужно, и когда сумма порядков была такой большой, что после вычитания BIAS имела вид в двоичном представлении: 100…00, причём старший единичный бит попадал в поле знака. Однако эта ситуация невозможна, так как делает ложной нашу первую проверку: перенос в этом случае появился до вычитания, так что в любом случае, мы ничего не теряем.
Для проверки самого старшего бита, мы будем проверять число на отрицательность, как если бы оно было знаковым целым. У отрицательных чисел самый старший бит установлен в 1.
Наш алгоритм приобретает следующий вид:

// Приближённое умножение
float apprMul(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;
	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 sum = (x.u + y.u) ^ sign;
	result.u = sum - EXP_BIAS;
	if(sum >= 0 && (i32)(result.u - EXP_ONE) < 0)
	{
		return 0;
	}
	result.u ^= sign;
	return result.f;
}

Переполнение (Overflow)

Переполнение возникает, если до вычитания константы BIAS был перенос в поле знака, и вычитание это не исправило. Либо если после вычитания поле порядка полностью состоит из единичных битов – это зарезервировано для специальных значений.
При переполнении результат умножения имеет очень большую степень, и может возникнуть искушение в качестве результата использовать бесконечность. Однако так делать опасно, потому что это – не настоящая бесконечность. Если мы сначала умножили число на 2, в результате чего оно вышло за диапазон, после чего поделили на 4 (иначе говоря, умножили на 0.25), то результат должен вернуться в представимый диапазон. Однако бесконечность после деления на 4 так бесконечностью и останется. Поэтому безопаснее избегать переполнения, чем пытаться его использовать.
Для упрощения алгоритма, при переполнении мы будем считать результатом неопределённость.
Неопределённость (Indefinite) – это специальное значение, которое в полях порядка и знака имеет только единичные биты, в поле мантиссы имеет установленным в 1 самый старший бит, а все остальные биты – нулевые. Формально, неопределённость относится к «тихим» нечислам (QNAN) – то есть, не вызывает исключений.
Аналогично случаю антипереполнения, чтобы свести две проверки поля порядка к одной, мы предварительно прибавим к порядку результата 1. Это вызовет перенос в поле знака, если поле порядка состояло только из единиц. Второй случай, когда бит переноса сбросится в 0, возможен лишь тогда, когда до вычитания BIAS сумма порядков была равна 126 = 0x7E, однако это невозможно, так как блокируется первой проверкой.
Теперь умножение выглядит так:

// Неопределённость
#define INDEFINITE	(SIGN_MASK | EXP_MASK | (1 << (BITS_PER_MANTISSA - 1)))

// Приближённое умножение
float apprMul(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;
	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 sum = (x.u + y.u) ^ sign;
	result.u = sum - EXP_BIAS;
	if(sum < 0)
	{
		if((i32)(result.u + EXP_ONE) < 0)
		{
			// Переполнение (Overflow)
			result.u = INDEFINITE;
			return result.f;
		}
	}
	else
	{
		if((i32)(result.u - EXP_ONE) < 0)
		{
			// Антипереполнение (Underflow)
			return 0;
		}
	}

	// Восстановление знака

	// Если мы оказались здесь, то поле знака не повреждено порядком,
	// и в текущеий момент равно 0.
	// Следовательно, вместо:
	// result.u ^= sign;
	// мы также можем записать:
	// result.u |= sign;
	// или
	// result.u += sign;

	result.u ^= sign;
	return result.f;
}

Умножение на ноль и денормализованные числа

Ноль и денормализованные числа имеют нулевое поле порядка. Причём ноль имеет нулевую мантиссу, а денормализованные числа – ненулевую. Для упрощения алгоритма, мы будем считать денормализованные числа равными нулю. Соответственно, если хотя бы один из сомножителей имел нулевой порядок, то мы будем считать результат равным нулю.
Проверить это можно следующими способами.
Первый способ явно использует операторы сравнения. Сначала мы маской выделяем поле порядка, затем проверяем его на равенство нулю:

if(((x.u & EXP_MASK ) == 0) | ((y.u & EXP_MASK) == 0))
	return 0;

Однако если процессор использует вместо инструкций присвоения по условию – условные переходы, то такой код может быть медленным. Поэтому рассмотрим другой способ проверки. Если мы сложим переменную с константой, полностью состоящей из единичных битов: 11…11, то перенос в старший разряд произойдёт тогда и только тогда, когда переменная не была равна нулю. Соответственно, мы сложим каждый из порядков с 11…11, после чего проверим, возник ли перенос при этом сложении. Если хотя бы у одного из порядков его не было, то порядок был нулевым, и результат умножения равен нулю.

if((i32)(((x.u & EXP_MASK ) + EXP_MASK) & ((y.u & EXP_MASK) + EXP_MASK)) >= 0)
	return 0;

Проверку порядков добавляем в начале функции умножения:

// Приближённое умножение
float apprMul(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;
	
	if(((x.u & EXP_MASK ) == 0) | ((y.u & EXP_MASK ) == 0))
		return 0;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 sum = (x.u + y.u) ^ sign;
	result.u = sum - EXP_BIAS;
	if(sum < 0)
	{
		if((i32)(result.u + EXP_ONE) < 0)
		{
			// Переполнение (Overflow)
			result.u = INDEFINITE;
			return result.f;
		}
	}
	else
	{
		if((i32)(result.u - EXP_ONE) < 0)
		{
			// Антипереполнение (Underflow)
			return 0;
		}
	}

	// Восстановление знака

	result.u ^= sign;
	return result.f;
}

Умножение на бесконечность, неопределённость и нечисла

Бесконечность, неопределённость и нечисла имеют поле порядка, полностью состоящее из единичных битов: 11…11. Для упрощения алгоритма будем считать результатом умножения на такие значения неопределённостью.
Проверить порядок сомножителей на равенство 11…11 можно, например, напрямую:

if(((x.u & EXP_MASK ) == EXP_MASK) | ((y.u & EXP_MASK) == EXP_MASK))
{
	result.u = INDEFINITE;
	return result.f;
}

Либо можно использовать следующий факт. Если сложить переменную с 1, то перенос возникнет тогда и только тогда, когда переменная полностью состояла из единичных битов: 11…11. Соответственно, если хотя бы у одного из сомножителей при сложении порядка с 1 возник перенос – результатом умножения будет неопределённость.

if((i32)(((x.u & EXP_MASK ) + EXP_ONE) | ((y.u & EXP_MASK) + EXP_ONE)) < 0)
{
	result.u = INDEFINITE;
	return result.f;
}

И теперь мы можем записать итоговую версию программы, которая работает со всеми возможными сомножителями:

// Приближённое умножение
float apprMul(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;

	u32 xe = x.u & EXP_MASK;
	u32 ye = y.u & EXP_MASK;

	if((i32)((xe + EXP_ONE) | (ye + EXP_ONE)) < 0)
	{
		result.u = INDEFINITE;
		return result.f;
	}

	if((i32)((xe + EXP_MASK) & (ye + EXP_MASK)) >= 0)
		return 0;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 sum = (x.u + y.u) ^ sign;
	result.u = sum - EXP_BIAS;
	if(sum < 0)
	{
		if((i32)(result.u + EXP_ONE) < 0)
		{
			// Переполнение (Overflow)
			result.u = INDEFINITE;
			return result.f;
		}
	}
	else
	{
		if((i32)(result.u - EXP_ONE) < 0)
		{
			// Антипереполнение (Underflow)
			return 0;
		}
	}

	// Восстановление знака

	result.u ^= sign;
	return result.f;
}

Умножение без операторов сравнения и переходов

Проведём серьёзную оптимизацию в несколько этапов. Для начала вместо неравенств запишем явную проверку поля знака:

// Приближённое умножение
float apprMul(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;

	u32 xe = x.u & EXP_MASK;
	u32 ye = y.u & EXP_MASK;

	if(((xe + EXP_ONE) | (ye + EXP_ONE)) & SIGN_MASK)
	{
		result.u = INDEFINITE;
		return result.f;
	}

	if((((xe + EXP_MASK) & (ye + EXP_MASK)) & SIGN_MASK) == 0)
		return 0;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 sum = (x.u + y.u) ^ sign;
	result.u = sum - EXP_BIAS;
	if(sum & SIGN_MASK)
	{
		if((result.u + EXP_ONE) & SIGN_MASK)
		{
			// Переполнение (Overflow)
			result.u = INDEFINITE;
			return result.f;
		}
	}
	else
	{
		if((result.u - EXP_ONE) & SIGN_MASK)
		{
			// Антипереполнение (Underflow)
			return 0;
		}
	}

	// Восстановление знака

	result.u ^= sign;
	return result.f;
}

Теперь избавимся от return в середине. Для этого проверку сомножителей придется перенести в конец функции:

// Приближённое умножение
float apprMul8(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 sum = (x.u + y.u) ^ sign;
	u32 t = sum - EXP_BIAS;
	result.u = t ^ sign;
	if(sum & SIGN_MASK)
	{
		if((t + EXP_ONE) & SIGN_MASK)
		{
			// Переполнение (Overflow)
			result.u = INDEFINITE;
		}
	}
	else
	{
		if((t - EXP_ONE) & SIGN_MASK)
		{
			// Антипереполнение (Underflow)
			result.u = 0;
		}
	}

	u32 xe = x.u & EXP_MASK;
	u32 ye = y.u & EXP_MASK;

	if((((xe + EXP_MASK) & (ye + EXP_MASK)) & SIGN_MASK) == 0)
		result.u = 0;

	if(((xe + EXP_ONE) | (ye + EXP_ONE)) & SIGN_MASK)
		result.u = INDEFINITE;

	return result.f;
}

Легко заметить, что все наши условия проверяют лишь поле знака. Поэтому, чтобы не загромождать код, напишем следующую вспомогательную функцию:

// Если expr < 0, то возращает ifTrue, иначе ifFalse
u32 ifSign(i32 expr, u32 ifTrue, u32 ifFalse)
{
	// Так как expr - знаковое целое, то сдвиг происходит
	// с распространением знака
	return ((ifTrue ^ ifFalse) & (expr >> (BITS_PER_ALL - 1))) ^ ifFalse;
}

Эта функция использует тот факт, что выражение в булевых переменных ((a ^ b) & c) ^ b возвращает a, если c = 1; и возвращает b, если c = 0. Вместо c мы используем маску, которая целиком состоит либо из нулей, либо из единиц.
Теперь можно записать функцию приближённого умножения без ветвлений и операторов сравнения, которая обрабатывает все возможные сомножители, в том числе – денормализованные числа, бесконечности, неопределённости и нечисла.

float apprMul(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 sum = (x.u + y.u) ^ sign;
	u32 t = sum - EXP_BIAS;
	result.u = t ^ sign;
	result.u = ifSign(sum,
		ifSign(t + EXP_ONE, INDEFINITE, result.u),
		ifSign(t - EXP_ONE, 0, result.u));

	u32 xe = x.u & EXP_MASK;
	u32 ye = y.u & EXP_MASK;

	result.u = ifSign((xe + EXP_MASK) & (ye + EXP_MASK), result.u, 0);
	result.u = ifSign((xe + EXP_ONE) | (ye + EXP_ONE), INDEFINITE, result.u);

	return result.f;
}

Строго говоря, вызов функции вида ifSign(…, …, 0) можно упростить, так как при этом происходит XOR с нулём. Но с этим оптимизирующий компилятор и сам может справиться.

Деление

Так как деление – операция, обратная к умножению, то из нашей формулы для приближённого умножения автоматически следует формула для приближённого деления.
Первый вариант программы запишем пока без учёта специальных значений.

float apprDiv(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;
	result.u = x.u - y.u + EXP_BIAS;
	return result.f;
}

Если при делении не возникает переполнения и других особых случаев, то эта функция работает и с отрицательными числами, как это было с аналогичной функцией умножения.

Приближённое вычисление обратного значения

Чтобы найти обратное значение переменной, нужно единицу разделить на эту переменную. В формате float у единицы нулевая мантисса и нулевое поле знака, а в поле порядка содержится BIAS. Соответственно, функция приближённого обратного значения имеет следующий вид:

float apprOneDiv(float yf)
{
	Float y, result;
	y.f = yf;
	result.u = EXP_BIAS * 2 - y.u;
	return result.f;
}

Эта функция работает как с положительными, так и с отрицательными числами, но предполагает, что параметр «хороший»: то есть, не равен нулю либо другим специальным значениям, и что результат не слишком велик, и не слишком мал.

Антипереполнение при делении

Антипереполнение возникает, когда результат деления слишком мал, и в поле порядка результата оказывается отрицательное число: соответственно, из поля знака происходит заём. Либо если в поле порядка оказывается ноль.
Признаком антипереполнения является следующая ситуация: после вычитания порядков произошёл заём (то есть, бит знака повреждён), и прибавление константы BIAS эту ситуацию не исправило, либо поле порядка обнулилось.
Точно так же, как при умножении, отнимем от порядка 1, чтобы вызвать заём, если порядок стал равен 0.
Функция деления принимает следующий вид:

float apprDiv(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;
	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 dif = (x.u - y.u) ^ sign;
	result.u = dif + EXP_BIAS;
	if(dif < 0 && (i32)(result.u - EXP_ONE) < 0)
	{
		// Антипереполнение (Underflow)
		return 0;
	}
	result.u ^= sign;
	return result.f;
}

Переполнение при делении

Переполнение возникает, если результат слишком велик. При вычитании порядок увеличиться не может, так как в поле порядка – неотрицательные значения. Соответственно, переполнение происходит, когда мы после вычитания к порядку прибавляем константу BIAS.
Следовательно, признаком переполнения является то, что после вычитания поле знака не повреждено, однако после сложения с BIAS возникает перенос в поле знака.
Аналогично умножению, прибавим к порядку 1, чтобы спровоцировать перенос, если поле порядка полностью состояло из единичных битов.

float apprDiv(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;
	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 dif = (x.u - y.u) ^ sign;
	result.u = dif + EXP_BIAS;
	if(dif < 0)
	{
		if((i32)(result.u - EXP_ONE) < 0)
		{
			// Антипереполнение (Underflow)
			return 0;
		}
	}
	else
	{
		if((i32)(result.u + EXP_ONE) < 0)
		{
			// Переполнение (Overflow)
			result.u = INDEFINITE;
			return result.f;
		}
	}
	result.u ^= sign;
	return result.f;
}

Специальные значения делимого и делителя

Если делимое равно нулю, либо денормализовано, то результатом будем считать ноль.
Деление на ноль приводит к неоднозначности, какую бесконечность использовать: положительную или отрицательную. Если заранее известно, что результат умножения – положительный, то можно было бы использовать положительную бесконечность. Но так как мы, для упрощения алгоритма, бесконечности вообще не используем, то при делении на ноль будем возвращать неопределённость.
Точно так же, как мы поступили при умножении – результатом операций с бесконечностью при делении мы будем считать неопределённость.
Операции с неопределённостью и нечислами тоже возвращают неопределённость. Особенностью деления является то, что сразу два значения порядка делителя приводят к неопределённому результату: когда порядок делителя имеет вид 00…00 или 11…11. Было бы удобно объединить эти две проверки. Если прибавить к порядку делителя 1, после чего обнулить перенос, то эти два значения порядка станут соседними, и их можно проверить одной операцией.
Функция приближённого деления, которая обрабатывает все возможные значения делимого и делителя, включая специальные, имеет вид:

float apprDiv(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;

	u32 xe = x.u & EXP_MASK;

	if((xe == EXP_MASK) | ((i32)((y.u + EXP_ONE) & EXP_MASK) <= EXP_ONE))
	{
		result.u = INDEFINITE;
		return result.f;
	}

	if(xe == 0)
		return 0;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 dif = (x.u - y.u) ^ sign;
	result.u = dif + EXP_BIAS;
	if(dif < 0)
	{
		if((i32)(result.u - EXP_ONE) < 0)
		{
			// Антипереполнение (Underflow)
			return 0;
		}
	}
	else
	{
		if((i32)(result.u + EXP_ONE) < 0)
		{
			// Переполнение (Overflow)
			result.u = INDEFINITE;
			return result.f;
		}
	}
	result.u ^= sign;
	return result.f;
}

Деление без операторов сравнения и переходов

Почти все оптимизации, которые нам понадобятся, мы уже делали в функции умножения. Новой является лишь оптимизация проверки порядка делителя:

(i32)((y.u + EXP_ONE) & EXP_MASK) <= EXP_ONE

Чтобы избавиться в ней от оператора сравнения, вычтем из обеих частей неравенства единицу из порядка, и ещё 1, чтобы получить строгое неравенство.

(i32)(((y.u + EXP_ONE) & EXP_MASK) - EXP_ONE - 1) < 0

И теперь нам нужно проверить здесь знаковый бит. Строго говоря, в компьютерной арифметике вычитание, как и сложение, не всегда сохраняет знак неравенства – из-за того, что целые числа расположены не на прямой, а на замкнутом кольце, и следом за наименьшим отрицательным числом идёт наибольшее положительное. Но у нас знак неравенства должен сохраниться, так как до вычитания старший бит числа – нулевой.

Пропустив промежуточные выкладки, которые аналогичны случаю умножения, запишем функцию приближённого деления без сравнений и переходов:

float apprDiv(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 dif = (x.u - y.u) ^ sign;
	u32 t = dif + EXP_BIAS;
	result.u = t ^ sign;
	result.u = ifSign(dif,
		ifSign(t - EXP_ONE, 0, result.u),
		ifSign(t + EXP_ONE, INDEFINITE, result.u));

	u32 xe = x.u & EXP_MASK;

	result.u = ifSign(xe + EXP_MASK, result.u, 0);
	result.u = ifSign((xe + EXP_ONE) | (((y.u + EXP_ONE) & EXP_MASK) - EXP_ONE - 1),
		INDEFINITE, result.u);

	return result.f;
}

Итоговая программа

Теперь соберём умножение, обратное значение и деление в одной программе, и проверим, как оно всё работает. В последней колонке результатов будем выводить приближённое значение.

#include <stdio.h>
#include <stdlib.h>

typedef signed __int32 i32;
typedef unsigned __int32 u32;

/*
	Число одинарной точности
	Размер: 32 бита
	Мантисса: 23 бита
	Порядок: 8 битов
*/

union Float
{
	i32 i;
	u32 u;
	float f;
};

#define BITS_PER_ALL		32
#define BITS_PER_MANTISSA	23

#define BIAS		0x7F
#define EXP_BIAS	(BIAS << BITS_PER_MANTISSA)
#define EXP_ONE		(1 << BITS_PER_MANTISSA)
#define EXP_MASK	(0xFF << BITS_PER_MANTISSA)

#define SIGN_MASK	(1 << (BITS_PER_ALL - 1))

#define INDEFINITE	(SIGN_MASK | EXP_MASK | (1 << (BITS_PER_MANTISSA - 1)))

// Если expr < 0, то возращает ifTrue, иначе ifFalse
inline u32 ifSign(i32 expr, u32 ifTrue, u32 ifFalse)
{
	// Так как expr - знаковое целое, то сдвиг происходит
	// с распространением знака
	return ((ifTrue ^ ifFalse) & (expr >> (BITS_PER_ALL - 1))) ^ ifFalse;
}

// Приближённое умножение
inline float apprMul(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 sum = (x.u + y.u) ^ sign;
	u32 t = sum - EXP_BIAS;
	result.u = t ^ sign;
	result.u = ifSign(sum,
		ifSign(t + EXP_ONE, INDEFINITE, result.u),
		ifSign(t - EXP_ONE, 0, result.u));

	u32 xe = x.u & EXP_MASK;
	u32 ye = y.u & EXP_MASK;

	result.u = ifSign((xe + EXP_MASK) & (ye + EXP_MASK), result.u, 0);
	result.u = ifSign((xe + EXP_ONE) | (ye + EXP_ONE), INDEFINITE, result.u);

	return result.f;
}

// Приближённое обратное значение
inline float apprOneDiv(float yf)
{
	Float y, result;
	y.f = yf;
	result.u = (EXP_BIAS * 2) - y.u;
	return result.f;
}

// Приближённое деление
inline float apprDiv(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 dif = (x.u - y.u) ^ sign;
	u32 t = dif + EXP_BIAS;
	result.u = t ^ sign;
	result.u = ifSign(dif,
		ifSign(t - EXP_ONE, 0, result.u),
		ifSign(t + EXP_ONE, INDEFINITE, result.u));

	u32 xe = x.u & EXP_MASK;

	result.u = ifSign(xe + EXP_MASK, result.u, 0);
	result.u = ifSign((xe + EXP_ONE) | (((y.u + EXP_ONE) & EXP_MASK) - EXP_ONE - 1),
		INDEFINITE, result.u);

	return result.f;
}

// Случайное число с плавающей точкой
float frand()
{
	Float result;
	result.u = (rand() << 17) + (rand() << 2);
	return result.f;
}

int main()
{
	FILE *fp = fopen("result.txt", "wt");

	float a, b;

	fputs("Умножение\n\n", fp);
	for(int i = 0; i < 20; i++)
	{
		a = frand();
		b = frand();
		float r = a * b;
		fprintf(fp, "%14g * %14g = %14g %14g\n", a, b, r, apprMul(a, b));
	}

	fputs("\n\nОбратное значение\n\n", fp);
	for(int i = 0; i < 20; i++)
	{
		b = frand();
		float r = 1 / b;
		fprintf(fp, "1 / %14g = %14g %14g\n", b, r, apprOneDiv(b));
	}

	fputs("\n\nДеление\n\n", fp);
	for(int i = 0; i < 20; i++)
	{
		a = frand();
		b = frand();
		float r = a / b;
		fprintf(fp, "%14g / %14g = %14g %14g\n", a, b, r, apprDiv(a, b));
	}

	fclose(fp);

	return 0;
}

В результате выполнения программы получили следующие результаты:

Умножение

   7.63402e-39 *    3.69062e-09 =              0              0
  -7.87436e-26 *    6.11476e+15 =   -4.81498e-10   -4.37925e-10
  -3.55394e+11 *    4.19839e-12 =       -1.49208       -1.44696
  -1.69148e-06 *    4.40465e+08 =       -745.039       -724.228
   7.45525e-25 *    3.16295e-16 =    2.35806e-40              0
  -1.18709e+37 *    1.33413e-20 =   -1.58373e+17   -1.56388e+17
   1.47492e-37 *   -4.66864e-34 =             -0              0
  -2.92024e-23 *    2.54598e-13 =   -7.43489e-36   -7.36049e-36
   1.82783e+31 *    3.82635e-30 =        69.9391        64.9598
       -277719 *   -7.08136e-36 =    1.96663e-30    1.95008e-30
  -5.51827e+19 *   -1.15262e+32 =            inf      -nan(ind)
  -7.04335e-33 *    1.02665e+33 =       -7.23104        -6.8986
  -4.19779e+17 *       -77663.9 =    3.26016e+22    3.10062e+22
  -2.31285e+38 *   -9.18401e-22 =    2.12412e+17    2.08048e+17
       666.378 *   -1.60581e+14 =   -1.07008e+17   -1.03944e+17
   5.30878e+19 *   -1.25119e-11 =   -6.64231e+08   -6.22006e+08
   1.42409e-37 *        20549.9 =    2.92649e-33    2.72499e-33
  -2.09028e-25 *   -6.52706e+13 =    1.36434e-11    1.35762e-11
   3.41468e+36 *      -0.223571 =   -7.63423e+35   -7.13149e+35
   2.40725e+32 *      -0.285579 =   -6.87459e+31   -6.59542e+31


Обратное значение

1 /   -1.13102e-23 =    -8.8416e+22   -9.75343e+22
1 /     -0.0103282 =       -96.8219       -107.391
1 /   -8.62841e-28 =   -1.15896e+27   -1.19576e+27
1 /    1.41866e-15 =    7.04891e+14    7.89668e+14
1 /    2.08521e+27 =    4.79568e-28    5.31359e-28
1 /   -2.65638e-36 =   -3.76452e+35   -4.10244e+35
1 /        -3.5733 =      -0.279853      -0.303338
1 /    6.93897e-13 =    1.44114e+12     1.6208e+12
1 /    3.92333e+37 =    2.54885e-38      2.716e-38
1 /   -3.79536e-08 =   -2.63479e+07   -2.89657e+07
1 /    3.39831e-16 =    2.94264e+15    3.30911e+15
1 /    4.15922e-17 =     2.4043e+16    2.70483e+16
1 /   -3.26845e-09 =   -3.05956e+08   -3.34273e+08
1 /    3.44574e+27 =    2.90213e-28    3.24789e-28
1 /   -4.93606e+33 =   -2.02591e-34   -2.11602e-34
1 /    2.60562e-15 =    3.83786e+14    4.31549e+14
1 /   -2.33178e+11 =   -4.28857e-12   -4.74176e-12
1 /    1.83613e-09 =    5.44624e+08    5.52156e+08
1 /   -4.80367e+10 =   -2.08174e-11   -2.33114e-11
1 /    5.85664e-15 =    1.70746e+14    1.90207e+14


Деление

   -0.00349138 /        75334.4 =   -4.63451e-08   -4.88185e-08
  -7.61755e-11 /   -9.86817e-05 =    7.71931e-07    8.06753e-07
   5.39424e-20 /   -2.36491e+13 =   -2.28094e-33   -2.53579e-33
  -2.57816e-36 /    3.02153e+27 =             -0              0
  -1.07626e+31 /   -2.69156e-31 =            inf      -nan(ind)
      3.16e+34 /   -1.85275e+15 =   -1.70557e+19   -1.73022e+19
   1.23307e+32 /   -1.68849e+36 =   -7.30279e-05   -7.62692e-05
    -0.0456871 /    1.01003e-08 =   -4.52333e+06   -4.64034e+06
  -1.25782e+14 /   -5.48995e-39 =            inf      -nan(ind)
   1.03812e-34 /   -9.34649e-28 =    -1.1107e-07   -1.14501e-07
  -3.98391e-34 /   -4.60153e+14 =              0              0
  -1.20124e-14 /   -3.30396e+14 =    3.63575e-29    3.82891e-29
   5.47548e-22 /    1.55613e+10 =    3.51864e-32    3.65164e-32
        772803 /    2.53124e-17 =    3.05306e+22    3.11686e+22
  -2.32844e+13 /   -1.75253e-23 =    1.32861e+36    1.32882e+36
       5.57153 /    9.74577e-39 =            inf      -nan(ind)
  -9.63453e+08 /    1.73026e-06 =   -5.56827e+14   -5.57395e+14
  -9.73926e-31 /    2.12961e-21 =   -4.57325e-10   -4.60422e-10
      -4.54288 /   -3.04582e-32 =    1.49152e+32    1.54162e+32
    1.2652e+29 /   -1.08604e-08 =   -1.16496e+37   -1.21145e+37

Повышение точности умножения и деления с помощью таблиц поправок

Я добавил эту главу уже после публикации статьи.

Так как для умножения и деления мы используем приближённую формулу $M≈\log_2 (1+M)$ , то мы можем повысить точность этого приближения, если разобьём диапазон значений M на интервалы, и для каждого интервала рассчитаем свою константную поправку. В качестве номера интервала будем использовать старшие биты мантиссы.
Всего нам понадобятся две таблицы: одна для приближённого преобразования M в $\log_2 (1+M)$ , вторая – для обратного преобразования. До приближённого умножения мы суммируем мантиссы множителей с соответствующими поправками, чтобы они лучше приближались к логарифмам. Затем суммируем множители – то есть, находим сумму логарифмов, затем результат суммируем с поправкой, которая от логарифма переходит назад к мантиссе.
Для нулевого интервала всегда будем использовать поправку равную нулю, чтобы не испортить умножение чисел, кратных степени двойки. На остальных интервалах будем следить за тем, чтобы поправка не вызвала переполнение, способное повредить поле порядка.
Для поиска лучшей, в среднем, поправки на заданном интервале, будем минимизировать интеграл квадрата отклонения.
Начнём с преобразования M в $\log_2 (1+M)$ . На интервале 0 < M < 1 логарифм лежит выше соответствующих значений мантиссы. Соответственно, нам понадобится прибавить к мантиссе какое-то число, которое назовём D. Обозначим концы интервала, на котором ищем поправку через a, b. Решим следующее уравнение:

$(\int\limits_a^b (x+D-\log_2 (1+x) )^2 dx)'_D = 0$

В результате получаем, что поправка на интервале [a, b]:

$D=\frac{(1+b)\ln(1+b)-(1+a)\ln(1+a)}{(b-a)\ln 2}-\frac{a+b}{2}-\frac{1}{\ln 2}$

Теперь найдём поправку для обратного преобразования. После сложения логарифмов у нас сейчас в поле мантиссы – логарифм, и нам нужно восстановить его аргумент. Иначе говоря, для $x=\log_2 (1+y)$ нам нужно вычислить y = 2^x-1 . Так как на интервале 0 < x < 1 функция лежит ниже, чем y=x, то поправку мы будем вычитать, чтобы раньше времени не иметь дела с отрицательными числами. Уже после расчётов знак поправки можно будет изменить.
Обозначим поправку обратного преобразования через F. Тогда нам нужно будет решить уравнение:

$(\int\limits_a^b (x-F-2^x+1)^2 dx)’_F=0$

Решив его, получаем, что поправка для обратного преобразования:

$F=\frac{a+b}{2}+1-\frac{2^b-2^a}{(b-a) \ln 2}$

Теперь перейдём к программной реализации умножения и деления с поправками. Добавилось две функции корректировки: одна будет корректировать параметры до умножения и деления, вторая – после этого. В самих функциях умножения и деления изменения незначительны, и связаны лишь с этими новыми функциями корректировки.
Самая сложная новая часть программы – это инициализация таблиц поправок. Для более точных вычислений я применил числа двойной точности, поэтому пришлось добавить константы и маски для работы с ними: они начинаются с префикса “DOUBLE_”. Числа двойной точности создали сложности с операцией сдвига, так как double занимает 64 бита, а в C++ оператор >> может не поддерживать сдвиг более чем на 31 бит. Поэтому пришлось написать свои функции сдвига. Я их не оптимизировал – и так сойдёт. Всё равно инициализация таблиц происходит только раз.
Чтобы можно было сравнить, как изменилась точность умножения – я оставил и функции умножения и деления без поправок, лишь заменив префиксы их названия на old.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef signed __int32 i32;
typedef unsigned __int32 u32;
typedef unsigned __int64 u64;

/*
	Число одинарной точности
	Размер: 32 бита
	Мантисса: 23 бита
	Порядок: 8 битов
*/

union Float
{
	i32 i;
	u32 u;
	float f;
};

#define BITS_PER_ALL		32
#define BITS_PER_MANTISSA	23
#define MANTISSA_MASK		((1 << BITS_PER_MANTISSA) - 1)

#define BIAS		0x7F
#define EXP_BIAS	(BIAS << BITS_PER_MANTISSA)
#define EXP_ONE		(1 << BITS_PER_MANTISSA)
#define EXP_MASK	(0xFF << BITS_PER_MANTISSA)

#define SIGN_MASK	(1 << (BITS_PER_ALL - 1))

#define INDEFINITE	(SIGN_MASK | EXP_MASK | (1 << (BITS_PER_MANTISSA - 1)))

// Количество старших битов мантиссы, которые будут использоваться
// как индекс в таблице поправок
#define BITS_PER_CORRECTION_INDEX	5
// Количество элементов в таблице поправок
#define CORRECTIONS_COUNT			(1 << BITS_PER_CORRECTION_INDEX)
#define CORR_SHIFT					(BITS_PER_MANTISSA - BITS_PER_CORRECTION_INDEX)
#define CORR_MASK					((CORRECTIONS_COUNT - 1) << CORR_SHIFT)

// Если expr < 0, то возращает ifTrue, иначе ifFalse
inline u32 ifSign(i32 expr, u32 ifTrue, u32 ifFalse)
{
	// Так как expr - знаковое целое, то сдвиг происходит
	// с распространением знака
	return ((ifTrue ^ ifFalse) & (expr >> (BITS_PER_ALL - 1))) ^ ifFalse;
}

///// Старые функции умножения и деления без поправок /////

// Приближённое умножение
inline float oldMul(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 sum = (x.u + y.u) ^ sign;
	u32 t = sum - EXP_BIAS;
	result.u = t ^ sign;
	result.u = ifSign(sum,
		ifSign(t + EXP_ONE, INDEFINITE, result.u),
		ifSign(t - EXP_ONE, 0, result.u));

	u32 xe = x.u & EXP_MASK;
	u32 ye = y.u & EXP_MASK;

	result.u = ifSign((xe + EXP_MASK) & (ye + EXP_MASK), result.u, 0);
	result.u = ifSign((xe + EXP_ONE) | (ye + EXP_ONE), INDEFINITE, result.u);

	return result.f;
}

// Приближённое обратное значение
inline float oldOneDiv(float yf)
{
	Float y, result;
	y.f = yf;
	result.u = (EXP_BIAS * 2) - y.u;
	return result.f;
}

// Приближённое деление
inline float oldDiv(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 dif = (x.u - y.u) ^ sign;
	u32 t = dif + EXP_BIAS;
	result.u = t ^ sign;
	result.u = ifSign(dif,
		ifSign(t - EXP_ONE, 0, result.u),
		ifSign(t + EXP_ONE, INDEFINITE, result.u));

	u32 xe = x.u & EXP_MASK;

	result.u = ifSign(xe + EXP_MASK, result.u, 0);
	result.u = ifSign((xe + EXP_ONE) | (((y.u + EXP_ONE) & EXP_MASK) - EXP_ONE),
		INDEFINITE, result.u);

	return result.f;
}

///// Новые функции умножения и деления с поправками /////

// Таблицы поправок
u32 tableCorrBefore[CORRECTIONS_COUNT];
u32 tableCorrAfter[CORRECTIONS_COUNT];

inline u32 corrBefore(u32 x)
{
	return x + tableCorrBefore[(x & CORR_MASK) >> CORR_SHIFT];
}

inline u32 corrAfter(u32 x)
{
	return x + tableCorrAfter[(x & CORR_MASK) >> CORR_SHIFT];
}

// Приближённое умножение
inline float apprMul(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 sum = (corrBefore(x.u) + corrBefore(y.u)) ^ sign;
	u32 t = corrAfter(sum) - EXP_BIAS;
	result.u = t ^ sign;
	result.u = ifSign(sum,
		ifSign(t + EXP_ONE, INDEFINITE, result.u),
		ifSign(t - EXP_ONE, 0, result.u));

	u32 xe = x.u & EXP_MASK;
	u32 ye = y.u & EXP_MASK;

	result.u = ifSign((xe + EXP_MASK) & (ye + EXP_MASK), result.u, 0);
	result.u = ifSign((xe + EXP_ONE) | (ye + EXP_ONE), INDEFINITE, result.u);

	return result.f;
}

// Приближённое обратное значение
inline float apprOneDiv(float yf)
{
	Float y, result;
	y.f = yf;
	result.u = corrAfter(EXP_BIAS * 2 - corrBefore(y.u));
	return result.f;
}

// Приближённое деление
inline float apprDiv(float xf, float yf)
{
	Float x, y, result;
	x.f = xf;
	y.f = yf;

	u32 sign = (x.u ^ y.u) & SIGN_MASK;
	i32 dif = (corrBefore(x.u) - corrBefore(y.u)) ^ sign;
	u32 t = corrAfter(dif) + EXP_BIAS;
	result.u = t ^ sign;
	result.u = ifSign(dif,
		ifSign(t - EXP_ONE, 0, result.u),
		ifSign(t + EXP_ONE, INDEFINITE, result.u));

	u32 xe = x.u & EXP_MASK;

	result.u = ifSign(xe + EXP_MASK, result.u, 0);
	result.u = ifSign((xe + EXP_ONE) | (((y.u + EXP_ONE) & EXP_MASK) - EXP_ONE - 1),
		INDEFINITE, result.u);

	return result.f;
}

///// Инициализация таблиц поправок /////

/*
	Число двойной точности
	Размер: 64 бита
	Мантисса: 52 бита
	Порядок: 11 битов
*/
union Double
{
	u64 u;
	double f;
};

#define DOUBLE_BITS_PER_MANTISSA	52
#define DOUBLE_BIAS					0x3FF
#define DOUBLE_EXP_ONE				0x10000000000000ull
#define DOUBLE_MANTISSA_MASK		0x0FFFFFFFFFFFFFull

u64 lshr(u64 u, i32 shift)
{
	while(shift > 16)
	{
		u >>= 16;
		shift -= 16;
	}
	return u >> shift;
}

u64 lshl(u64 u, i32 shift)
{
	while(shift > 16)
	{
		u <<= 16;
		shift -= 16;
	}
	return u << shift;
}

double evalCorrectionBefore(double a, double b)
{
	return ((1 + b) * log(1 + b) - (1 + a) * log(1 + a))
		/ ((b - a) * log((double)2)) - (a + b) / 2 - 1 / log((double)2);
}

double evalCorrectionAfter(double a, double b)
{
	return (a + b) / 2 + 1 - (exp(b * log((double)2)) - exp(a * log((double)2)))
		/ ((b - a) * log((double)2));
}

// Выравнивание поправки
u32 alignCorrection(double corr)
{
	Double c;
	u32 result;
	c.f = corr;

	i32 shift = (u32)(DOUBLE_BIAS - lshr(c.u, DOUBLE_BITS_PER_MANTISSA)
		+ DOUBLE_BITS_PER_MANTISSA - BITS_PER_MANTISSA) - 1;

	/*
		Округление к чётному:
		(среди Y хотя бы один ненулевой бит)
		XXX|0		-> XXX
		XXX|0YYY	-> XXX
		XX0|1		-> XX0
		XX0|1YYY	-> XX0 + 1
		XX1|1		-> XX1 + 1
		XX1|1YYY	-> XX1 + 1
	*/

	c.u = (c.u & DOUBLE_MANTISSA_MASK) | DOUBLE_EXP_ONE;
	result = (u32)lshr(c.u, shift);
	u32 isEjectedOne = result & 1;
	result >>= 1;

	// Весь следующий if(...) ... можно заменить выражением
	// result += (((c.u & (lshl(1, shift) - 1)) != 0) | result) & isEjectedOne;

	if(isEjectedOne)
	{
		if(result & 1)
			result++;
		else if(c.u & (lshl(1, shift) - 1))
			result++;
	}
	return result;
}

i32 minof(u32 x, u32 y)
{
	return x < y? x: y;
}

void initCorrectionTables()
{
	u32 guardIndexMask = 2;
	u32 indexShift = BITS_PER_MANTISSA;
	u32 shiftedExp = (BIAS - BITS_PER_CORRECTION_INDEX) << BITS_PER_MANTISSA;

	tableCorrBefore[0] = 0;
	tableCorrAfter[0] = 0;

	for(u32 index = 1; index < CORRECTIONS_COUNT; index++)
	{
		if(index & guardIndexMask)
		{
			// Размер index в битах увеличился
			guardIndexMask <<= 1;
			indexShift--;
			shiftedExp += EXP_ONE;
		}

		// Вычисляем начало и конец интервала, к которому
		// будем находить поправку

		Float a, b;

		/*
			При инициализации поправок нам нужно учесть,
			что при сложении логарифмов мы работаем с мантиссой так,
			как если бы она была ненормализованным числом
			фиксированной точности: то есть, имела явную лидирующую единицу.

			Поэтому перед тем, как вычислять поправку для интервала –
			мы концы интервала нормализуем, сдвигая их мантиссы так,
			чтобы лидирующая единица попала в поле порядка,
			затем её срезаем операцией (... & MANTISSA_MASK).
			Вычислив поправку, в конце итерации мы восстанавливаем
			лидирующую единицу с помощью функции alignCorrection
		*/
		a.u = ((index << indexShift) & MANTISSA_MASK) + shiftedExp;
		b.u = (((((index + 1) << CORR_SHIFT) - 1)
				<< (indexShift - CORR_SHIFT)) & MANTISSA_MASK) + shiftedExp;

		// Вычисляем максимально допустимое значение поправки,
		// которое не повредит порядок.
		// Эти числа являются числами фиксированной точности
		// с явной лидирующей единицей.
		u32 maxAdd = EXP_ONE - ((index + 1) << CORR_SHIFT);
		u32 maxSub = index << CORR_SHIFT;

		tableCorrBefore[index] = minof(maxAdd, 
			alignCorrection(evalCorrectionBefore(a.f, b.f)));
		tableCorrAfter[index] = -minof(maxSub, 
			alignCorrection(evalCorrectionAfter(a.f, b.f)));
	}
}

///// Тестирование /////

// Случайное число с плавающей точкой
float frand()
{
	Float result;
	result.u = (rand() << 17) + (rand() << 2);
	return result.f;
}

// Возвращает отклонение параметров друг от друга в процентах
float er(float a, float b)
{
	float r;
	a = fabs(a);
	b = fabs(b);
	if(a < b) r = a / b; else r = b / a;
	return (1 - r) * 100;
}

int main()
{
	float a, b;

	FILE *fp = fopen("result.txt", "wt");

	initCorrectionTables();

	fputs("Умножение\n\n"
"                                                            new            old       new error       old error\n"
		, fp);
	for(int i = 0; i < 20; i++)
	{
		a = frand();
		b = frand();
		float r = a * b;
		float r2 = apprMul(a, b);
		float r3 = oldMul(a, b);
		float e2 = er(r, r2);
		float e3 = er(r, r3);
		fprintf(fp, "%14g * %14g = %14g %14g %14g %14.5g%% %14.5g%% %s\n",
				a, b, r, r2, r3, e2, e3, e2 > e3? "++": "--");
	}

	fputs("\n\nОбратное значение\n\n"
"                                               new            old       new error       old error\n"
		, fp);
	for(int i = 0; i < 20; i++)
	{
		b = frand();
		float r = 1 / b;
		float r2 = apprOneDiv(b);
		float r3 = oldOneDiv(b);
		float e2 = er(r, r2);
		float e3 = er(r, r3);
		fprintf(fp, "1 / %14g = %14g %14g %14g %14.5g%% %14.5g%% %s\n",
				b, r, r2, r3, e2, e3, e2 > e3? "++": "--");
	}

	fputs("\n\nДеление\n\n"
"                                                            new            old       new error       old error\n"
		, fp);
	for(int i = 0; i < 20; i++)
	{
		a = frand();
		b = frand();
		float r = a / b;
		float r2 = apprDiv(a, b);
		float r3 = oldDiv(a, b);
		float e2 = er(r, r2);
		float e3 = er(r, r3);
		fprintf(fp, "%14g / %14g = %14g %14g %14g %14.5g%% %14.5g%% %s\n",
				a, b, r, r2, r3, e2, e3, e2 > e3? "++": "--");
	}

	fclose(fp);

	return 0;
}

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

Умножение

                                                            new            old       new error       old error
   7.63402e-39 *    3.69062e-09 =              0              0              0      -nan(ind)%      -nan(ind)% --
  -7.87436e-26 *    6.11476e+15 =   -4.81498e-10   -4.81885e-10   -4.37925e-10       0.080359%         9.0496% --
  -3.55394e+11 *    4.19839e-12 =       -1.49208       -1.48906       -1.44696        0.20283%         3.0241% --
  -1.69148e-06 *    4.40465e+08 =       -745.039       -745.837       -724.228        0.10704%         2.7932% --
   7.45525e-25 *    3.16295e-16 =    2.35806e-40              0              0            100%            100% --
  -1.18709e+37 *    1.33413e-20 =   -1.58373e+17   -1.56691e+17   -1.56388e+17         1.0621%         1.2535% --
   1.47492e-37 *   -4.66864e-34 =             -0              0              0      -nan(ind)%      -nan(ind)% --
  -2.92024e-23 *    2.54598e-13 =   -7.43489e-36   -7.43439e-36   -7.36049e-36      0.0066876%         1.0006% --
   1.82783e+31 *    3.82635e-30 =        69.9391        69.6956        64.9598        0.34822%         7.1195% --
       -277719 *   -7.08136e-36 =    1.96663e-30    1.96221e-30    1.95008e-30        0.22472%        0.84174% --
  -5.51827e+19 *   -1.15262e+32 =            inf      -nan(ind)      -nan(ind)            nan%            nan% --
  -7.04335e-33 *    1.02665e+33 =       -7.23104       -7.23512        -6.8986       0.056398%         4.5974% --
  -4.19779e+17 *       -77663.9 =    3.26016e+22    3.25816e+22    3.10062e+22       0.061423%         4.8936% --
  -2.31285e+38 *   -9.18401e-22 =    2.12412e+17    2.12078e+17    2.08048e+17        0.15756%         2.0544% --
       666.378 *   -1.60581e+14 =   -1.07008e+17   -1.06978e+17   -1.03944e+17       0.028098%         2.8628% --
   5.30878e+19 *   -1.25119e-11 =   -6.64231e+08   -6.63697e+08   -6.22006e+08       0.080395%         6.3571% --
   1.42409e-37 *        20549.9 =    2.92649e-33    2.92712e-33    2.72499e-33       0.021404%         6.8854% --
  -2.09028e-25 *   -6.52706e+13 =    1.36434e-11    1.35676e-11    1.35762e-11        0.55561%        0.49247% ++
   3.41468e+36 *      -0.223571 =   -7.63423e+35   -7.62871e+35   -7.13149e+35       0.072306%         6.5853% --
   2.40725e+32 *      -0.285579 =   -6.87459e+31   -6.87582e+31   -6.59542e+31        0.01778%         4.0609% --


Обратное значение

                                               new            old       new error       old error
1 /   -1.13102e-23 =    -8.8416e+22   -8.82411e+22   -9.75343e+22         0.1978%         9.3488% --
1 /     -0.0103282 =       -96.8219       -96.8324       -107.391       0.010777%         9.8417% --
1 /   -8.62841e-28 =   -1.15896e+27   -1.15431e+27   -1.19576e+27        0.40098%         3.0774% --
1 /    1.41866e-15 =    7.04891e+14    7.05326e+14    7.89668e+14       0.061578%         10.736% --
1 /    2.08521e+27 =    4.79568e-28    4.79851e-28    5.31359e-28       0.059068%         9.7469% --
1 /   -2.65638e-36 =   -3.76452e+35   -3.77057e+35   -4.10244e+35        0.16052%          8.237% --
1 /        -3.5733 =      -0.279853      -0.279856      -0.303338     0.00077486%         7.7419% --
1 /    6.93897e-13 =    1.44114e+12    1.44059e+12     1.6208e+12        0.03826%         11.085% --
1 /    3.92333e+37 =    2.54885e-38    2.56079e-38      2.716e-38        0.46604%         6.1541% --
1 /   -3.79536e-08 =   -2.63479e+07   -2.63521e+07   -2.89657e+07       0.015867%         9.0373% --
1 /    3.39831e-16 =    2.94264e+15    2.94004e+15    3.30911e+15       0.088561%         11.074% --
1 /    4.15922e-17 =     2.4043e+16    2.40285e+16    2.70483e+16       0.060326%         11.111% --
1 /   -3.26845e-09 =   -3.05956e+08   -3.05699e+08   -3.34273e+08       0.083935%         8.4714% --
1 /    3.44574e+27 =    2.90213e-28    2.90243e-28    3.24789e-28        0.01027%         10.646% --
1 /   -4.93606e+33 =   -2.02591e-34   -2.01841e-34   -2.11602e-34        0.37006%         4.2585% --
1 /    2.60562e-15 =    3.83786e+14    3.83669e+14    4.31549e+14       0.030547%         11.068% --
1 /   -2.33178e+11 =   -4.28857e-12   -4.29431e-12   -4.74176e-12         0.1338%         9.5575% --
1 /    1.83613e-09 =    5.44624e+08    5.52156e+08    5.52156e+08         1.3641%         1.3641% --
1 /   -4.80367e+10 =   -2.08174e-11   -2.08222e-11   -2.33114e-11        0.02265%         10.699% --
1 /    5.85664e-15 =    1.70746e+14     1.7092e+14    1.90207e+14        0.10151%         10.231% --


Деление

                                                            new            old       new error       old error
   -0.00349138 /        75334.4 =   -4.63451e-08   -4.63703e-08   -4.88185e-08       0.054222%         5.0664% --
  -7.61755e-11 /   -9.86817e-05 =    7.71931e-07    7.71448e-07    8.06753e-07       0.062531%         4.3163% --
   5.39424e-20 /   -2.36491e+13 =   -2.28094e-33   -2.27469e-33   -2.53579e-33        0.27416%          10.05% --
  -2.57816e-36 /    3.02153e+27 =             -0              0              0      -nan(ind)%      -nan(ind)% --
  -1.07626e+31 /   -2.69156e-31 =            inf      -nan(ind)      -nan(ind)            nan%            nan% --
      3.16e+34 /   -1.85275e+15 =   -1.70557e+19   -1.70629e+19   -1.73022e+19       0.042206%         1.4247% --
   1.23307e+32 /   -1.68849e+36 =   -7.30279e-05   -7.30117e-05   -7.62692e-05       0.022089%         4.2499% --
    -0.0456871 /    1.01003e-08 =   -4.52333e+06   -4.52238e+06   -4.64034e+06       0.020891%         2.5217% --
  -1.25782e+14 /   -5.48995e-39 =            inf      -nan(ind)      -nan(ind)            nan%            nan% --
   1.03812e-34 /   -9.34649e-28 =    -1.1107e-07   -1.10736e-07   -1.14501e-07        0.30091%         2.9961% --
  -3.98391e-34 /   -4.60153e+14 =              0              0              0      -nan(ind)%      -nan(ind)% --
  -1.20124e-14 /   -3.30396e+14 =    3.63575e-29     3.6324e-29    3.82891e-29       0.092065%         5.0448% --
   5.47548e-22 /    1.55613e+10 =    3.51864e-32    3.51248e-32    3.65164e-32        0.17512%         3.6422% --
        772803 /    2.53124e-17 =    3.05306e+22    3.05673e+22    3.11686e+22        0.12006%          2.047% --
  -2.32844e+13 /   -1.75253e-23 =    1.32861e+36    1.32491e+36    1.32882e+36        0.27851%        0.01564% ++
       5.57153 /    9.74577e-39 =            inf      -nan(ind)      -nan(ind)            nan%            nan% --
  -9.63453e+08 /    1.73026e-06 =   -5.56827e+14   -5.57534e+14   -5.57395e+14        0.12693%          0.102% ++
  -9.73926e-31 /    2.12961e-21 =   -4.57325e-10    -4.5793e-10   -4.60422e-10        0.13206%        0.67251% --
      -4.54288 /   -3.04582e-32 =    1.49152e+32    1.49552e+32    1.54162e+32        0.26761%           3.25% --
    1.2652e+29 /   -1.08604e-08 =   -1.16496e+37   -1.16134e+37   -1.21145e+37        0.31089%         3.8375% --

Хабы: