TL; DR Да.
В зависимости от возможностей и недостатков аппаратного обеспечения, цепочки инструментов и математической библиотеки вашей платформы, это может не обязательно предусматривать самый быстрый или наименее болезненный способ вычисления квадратного корня двойной точности. Ниже я показываю простой подход, основанный на связанной итерации Арнольда Шенхаге для квадратного корня и обратного квадратного корня:
Начиная с аппроксимации обратного квадратного корня r ок ~ = 1 / √a, мы вычисляем s 0 = a * r ок и r 0 = r приблизительно / 2, затем итерация:
s i + 1 = s i + r i * (a - s i * s i )
r i + 1 = r i + r i * (1 - r i * 2 * s i + 1 )
где s i - приближения к √a, а r i - приближения к 1 / (2√a). Эта итерация представляет собой итерацию Ньютона-Рафсона, искусно переставленную и имеющую квадратичную сходимость, что означает, что каждый шаг примерно удваивает количество правильных битов. Начиная с одинарной точности r ок , для достижения точности двойной точности требуется всего два шага.
Если мы теперь воспользуемся объединенной операцией умножения-сложения (FMA), поддерживаемой обычными современными процессорами и обычно доступной через функцию fma(), для каждого полшага требуется только две FMA. В качестве дополнительного бонуса нам не нужна специальная логика округления, поскольку FMA вычисляет a*b+c, используя полный продукт a*b, без применения усечения или округления. Результирующий код, приведенный здесь в версии ISO C99, короткий и приятный:
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <fenv.h>
#include <math.h>
double my_sqrt (double a)
{
double b, r, v, w;
float bb, rr, ss;
int e, t, f;
if ((a <= 0) || isinf (a) || isnan (a)) {
if (a < 0) {
r = 0.0 / 0.0;
} else {
r = a + a;
}
} else {
/* compute exponent adjustments */
b = frexp (a, &e);
t = e - 2*512;
f = t / 2;
t = t - 2 * f;
f = f + 512;
/* map argument into the primary approximation interval [0.25,1) */
b = ldexp (b, t);
bb = (float)b;
/* compute reciprocal square root */
ss = 1.0f / bb;
rr = sqrtf (ss);
r = (double)rr;
/* Use A. Schoenhage's coupled iteration for the square root */
v = 0.5 * r;
w = b * r;
w = fma (fma (w, -w, b), v, w);
v = fma (fma (r, -w, 1), v, v);
w = fma (fma (w, -w, b), v, w);
/* map back from primary approximation interval by jamming exponent */
r = ldexp (w, f);
}
return r;
}
/* Professor George Marsaglia's 64-bit KISS PRNG */
static uint64_t xx = 1234567890987654321ULL;
static uint64_t cc = 123456123456123456ULL;
static uint64_t yy = 362436362436362436ULL;
static uint64_t zz = 1066149217761810ULL;
static uint64_t tt;
#define MWC64 (tt = (xx << 58) + cc, cc = (xx >> 6), xx += tt, cc += (xx < tt), xx)
#define XSH64 (yy ^= (yy << 13), yy ^= (yy >> 17), yy ^= (yy << 43))
#define CNG64 (zz = 6906969069ULL * zz + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)
int main (void)
{
volatile union {
double f;
unsigned long long int i;
} arg, res, ref;
unsigned long long int count = 0ULL;
do {
arg.i = KISS64;
ref.f = sqrt (arg.f);
res.f = my_sqrt (arg.f);
if (res.i != ref.i) {
printf ("\n!!!! arg=% 23.16e %016llx res=% 23.16e %016llx ref=% 23.16e %016llx\n",
arg.f, arg.i, res.f, res.i, ref.f, ref.i);
}
count++;
if ((count & 0xffffff) == 0) printf ("\rtests = %llu", count);
} while (1);
return EXIT_SUCCESS;
}
Для исчерпывающего тестирования этого кода на двух последовательных бинардах потребуется около недели или около того небольшого кластера машин, здесь я включил быстрый тест "дыма" с использованием случайных операндов.
На оборудовании, которое не поддерживает операцию FMA, fma() будет основано на эмуляции. Это медленно, и несколько таких эмуляций оказались ошибочными. Итерация Шёнхаге будет прекрасно работать без FMA, но в этом случае необходимо будет добавить дополнительную логику округления. В тех случаях, когда поддерживаются усеченные (с округлением до нуля) умножения с плавающей запятой, наиболее простым решением является использование округления Таккермана . В противном случае, вероятно, потребуется повторно интерпретировать аргумент двойной точности и предварительный результат в 64-разрядное целое число и выполнить округление с помощью целочисленных операций.