Если ваше оборудование имеет достаточно быстрый целочисленный множитель (скажем, 4-5 циклов), использование итеративного подхода для вычисления rece = 1 / делитель, вероятно, будет самым быстрым.Затем вы бы вычислили частное как дивиденд * получатель.Это очень полезно для необходимых вычислений с фиксированной запятой, если аппаратное обеспечение предлагает либо целочисленное умножение с результатом двойной ширины (т. Е. Полный продукт), либо инструкцию «mulhi», которая доставляет старшие 32 бита произведения двух 32-битных целых чисел.
Вам нужно будет масштабировать вычисления с фиксированной запятой на лету, чтобы сохранить промежуточные результаты, близкие к 32 битам, чтобы получить результат с точностью до 24 бит, которые требуются для конечного результата после округления.
Самый быстрый подход, вероятно, заключается в создании 9-разрядного начального приближения "приблизительно" к 1 / делителю с последующей кубически сходящейся итерацией для обратной величины:
e = 1 - divisor * approx
e = e * e + e
recip = e * approx + approx
Проще всего предварительно вычислитьначальное приближение и сохранить его в массиве из 256 байтов, проиндексированных битами 23:16 делителя (т.е. 8 старших значащих битов мантиссы).Поскольку все значения аппроксимации являются числами в диапазоне 0x100 ... 0x1FF (что соответствует 0,5–0,998046875), нет необходимости сохранять старший значащий бит каждого значения, поскольку он будет равен «1».Просто добавьте 0x100 к полученному элементу таблицы, чтобы получить окончательное значение начального приближения.
Если вы не можете позволить себе 256 байтов хранения таблицы, альтернативным способом создания начального приближения с точностью до 9 бит будет полином.степени 3, которая приближается к 1 / (1 + f), где f - дробная часть делителя мантиссы, m.Поскольку в IEEE-754 известно, что m находится в [1.0,2.0), f находится в [0.0,1.0).
Корректное округление до ближайшего или даже (если требуется) может быть реализовано с помощьюУмножение обратного множителя предварительного отношения делителем для установления остатка и выбор окончательного отношения так, чтобы остаток был минимальным.
Следующий код демонстрирует принципы аппроксимации, обсужденные выше, используя более простой случай обратнойс надлежащим округлением в соответствии с IEEE-754 режимом ближнего или четного и с соответствующей обработкой особых случаев (ноль, денормали, бесконечность, NaNs).Предполагается, что 32-разрядное значение IEEE-754 с плавающей запятой одинарной точности может быть передано по битам из 32-разрядного целого без знака и обратно.Затем все операции выполняются над 32-разрядными целыми числами.
unsigned int frcp_rn_core (unsigned int z)
{
unsigned int x, y;
int sign;
int expo;
sign = z & 0x80000000;
expo = (z >> 23) & 0xff;
x = expo - 1;
if (x > 0xfd) {
/* handle special cases */
x = z << 1;
/* zero or small denormal returns infinity of like sign */
if (x <= 0x00400000) {
return sign | 0x7f800000;
}
/* infinity returns zero of like sign */
else if (x == 0xff000000) {
return sign;
}
/* convert SNaNs to QNaNs */
else if (x > 0xff000000) {
return z | 0x00400000;
}
/* large denormal, normalize it */
else {
y = x < 0x00800000;
z = x << y;
expo = expo - y;
}
}
y = z & 0x007fffff;
#if USE_TABLE
z = 256 + rcp_tab[y >> 15]; /* approx */
#else
x = y << 3;
z = 0xe39ad7a0;
z = mul_32_hi (x, z) + 0x0154bde4;
z = mul_32_hi (x, z) + 0xfff87521;
z = mul_32_hi (x, z) + 0x00001ffa;
z = z >> 4;
#endif /* USE_TABLE */
y = y | 0x00800000;
/* cubically convergent approximation to reciprocal */
x = (unsigned int)-(int)(y * z); /* x = 1 - arg * approx */
x = mul_32_hi (x, x) + x; /* x = x * x + x */
z = z << 15;
z = mul_32_hi (x, z) + z; /* approx = x * approx + approx */
/* compute result exponent */
expo = 252 - expo;
if (expo >= 0) {
/* result is a normal */
x = y * z + y;
z = (expo << 23) + z;
z = z | sign;
/* round result */
if ((int)x <= (int)(y >> 1)) {
z++;
}
return z;
}
/* result is a denormal */
expo = -expo;
z = z >> expo;
x = y * z + y;
z = z | sign;
/* round result */
if ((int)x <= (int)(y >> 1)) {
z++;
}
return z;
}
Функция mul_32_high () является заполнителем для машинно-зависимой операции, которая возвращает старшие 32 бита умножения со знаком на два 32-разрядных бита.целые числа.Полупереносимая реализация вместо машинной версии:
/* 32-bit int, 64-bit long long int */
int mul_32_hi (int a, int b)
{
return (int)(unsigned int)(((unsigned long long)(((long long)a)*b)) >> 32);
}
. Взаимную таблицу с 256 записями, используемую в табличном варианте, можно построить следующим образом:
static unsigned char rcp_tab[256];
for (int i = 0; i < 256; i++) {
rcp_tab[i] = (unsigned char)(((1./(1.+((i+.5)/256.)))*512.+.5)-256.);
}