Алгоритмы с квадратным корнем с плавающей точкой

Мне интересно, кто-нибудь может дать мне несколько примеров алгоритмов квадратного корня с плавающей точкой, которые могут использовать аппаратный делитель.

Дополнительные детали:
У меня есть модуль с плавающей запятой, который я разрабатываю и который оснащен аппаратным 32-разрядным умножителем, сумматором и делителем IEEE-754 с плавающей запятой. Я уже реализовал квадратный корень, используя метод Ньютона-Рафсона, используя только умножение и сложение / вычитание, но теперь я хочу сравнить пропускную способность квадратного корня, если у меня есть аппаратный делитель.

Одним конкретным входом, который трудно точно вычислить, является квадратный корень из 0x7F7FFFFF (3.4028234663852886E38).

1

Решение

Решение, предоставленное @tmyklebu, безусловно, соответствует вашим требованиям.

r = input value
s(0) = initial estimate of sqrt(r).  Example: r with its exponent halved.
s(n) = sqrt(r)

s <- (s + r / s) / 2

Имеет квадратичную сходимость, выполняет запрошенное деление. N = 3 или 4 должны сделать это для 32-разрядного числа с плавающей запятой.

[Изменить N = 2 для 32-разрядного числа с плавающей запятой, N = 3 (возможно, 4) для двойного]
[Редактировать за запрос OP] [Редактировать Добавленные комментарии для запроса OP]
// Initial estimate
static double S0(double R) {
double OneOverRoot2 = 0.70710678118654752440084436210485;
double Root2        = 1.4142135623730950488016887242097;
int Expo;
// Break R into mantissa and exponent parts.
double Mantissa = frexp(R, &Expo);
int j;
printf("S0 %le %d %le\n", Mantissa, Expo, frexp(sqrt(R), &j));
// If exponent is odd ...
if (Expo & 1) {
// Pretend the mantissa [0.5 ... 1.0) is multiplied by 2 as Expo is odd,
//   so it now has the value [1.0 ... 2.0)
// Estimate the sqrt(mantissa) as [1.0 ... sqrt(2))
// IOW: linearly map (0.5 ... 1.0) to (1.0 ... sqrt(2))
Mantissa = (Root2 - 1.0)/(1.0 - 0.5)*(Mantissa - 0.5) + 1.0;
}
else {
// The mantissa is in range [0.5 ... 1.0)
// Estimate the sqrt(mantissa) as [1/sqrt(2) ... 1.0)
// IOW: linearly map (0.5 ... 1.0) to (1/sqrt(2) ... 1.0)
Mantissa = (1.0 - OneOverRoot2)/(1.0 - 0.5)*(Mantissa - 0.5) + OneOverRoot2;
}
// Form initial estimate by using the above mantissa estimate and exponent/2
return ldexp(Mantissa, Expo/2);
}

// S = (S + R/S)/2 method
double Sqrt(double R) {
double S = S0(R);
int i = 5;  // May be reduced to 3 or 4 for double and 2 for float
do {
printf("S  %u %le %le\n", 5-i, S, (S-sqrt(R))/sqrt(R));
S = (S + R/S)/2;
} while (--i);
return S;
}

void STest(double x) {
printf("T  %le %le %le\n", x, Sqrt(x), sqrt(x));
}

int main(void) {
STest(612000000000.0);
return 0;
}

Сходится после 3 итераций за двойное.

S0 5,566108e-01 40 7,460635e-01
S 0 7,762279e + 05 -7,767318e-03
S 1 7,823281e + 05 3.040175e-05
S 2 7,823043e + 05 4,621193e-10
S 3 7,823043e + 05 0,000000e + 00
S 4 7,823043e + 05 0,000000e + 00
Т 6.120000e + 11 7.823043e + 05 7.823043e + 05

1

Другие решения

Других решений пока нет …