Определение квадратного корня с плавающей точкой

Question

Как определить квадратный корень числа с плавающей запятой? Является ли метод Ньютона-Рафсона хорошим способом? У меня тоже нет аппаратного квадратного корня. У меня также нет аппаратного деления (но я реализовал деление с плавающей запятой).

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

Кроме того, какой должна быть первоначальная догадка, чтобы уменьшить общее количество итераций ???

Большое спасибо!

Answer 1

Когда вы используете Ньютона-Рафсона для вычисления квадратного корня, вы на самом деле хотите использовать итерацию, чтобы найти обратный квадратный корень (после чего вы можете просто умножить на вход - с некоторой тщательностью для округления - для полученияквадратный корень).

Точнее: мы используем функцию f(x) = x^-2 - n.Понятно, что если f(x) = 0, то x = 1/sqrt(n).Это приводит к итерации Ньютона:

x_(i+1) = x_i - f(x_i)/f'(x_i)
        = x_i - (x_i^-2 - n)/(-2x_i^-3)
        = x_i + (x_i - nx_i^3)/2
        = x_i*(3/2 - 1/2 nx_i^2)

Обратите внимание, что (в отличие от итерации для квадратного корня), эта итерация для обратного квадратного корня не включает делений, поэтому в целом она намного эффективнее.

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

---

Редактировать: спрашивающий, похоже, все еще сбит с толку, поэтому давайте поработаем с примером: sqrt(612).612 - это 1.1953125 x 2^9 (или b1.0011001 x 2^9, если вы предпочитаете двоичный файл).Вытяните четную часть показателя степени (9), чтобы записать входные данные как f * 2^(2m), где m является целым числом, а f находится в диапазоне [1,4).Тогда у нас будет:

sqrt(n) = sqrt(f * 2^2m) = sqrt(f)*2^m

, применяя это сокращение к нашему примеру, получим f = 1.1953125 * 2 = 2.390625 (b10.011001) и m = 4.Теперь выполните итерацию Ньютона-Рафсона, чтобы найти x = 1/sqrt(f), используя начальное предположение 0,5 (как я отметил в комментарии, это предположение сходится для всех f, но вы можете значительно улучшить результаты, используя линейное приближение в качестве исходного предположения):

x_0 = 0.5
x_1 = x_0*(3/2 - 1/2 * 2.390625 * x_0^2)
    = 0.6005859...
x_2 = x_1*(3/2 - 1/2 * 2.390625 * x_1^2)
    = 0.6419342...
x_3 = 0.6467077...
x_4 = 0.6467616...

Так что даже при (относительно плохом) начальном предположении мы получаем быструю сходимость к истинному значению 1/sqrt(f) = 0.6467616600226026.

Теперь мы просто собираем конечный результат:

sqrt(f) = x_n * f = 1.5461646...
sqrt(n) = sqrt(f) * 2^m = 24.738633...

И проверьте: sqrt (612) = 24.738633 ...

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

Если вы хотите исправить округление без явной проверки остатка, вам нужно вычислить sqrt (f) с точностью до 2p + 2 бит (где pточность источника и типа назначения).Тем не менее, вы также можете применить стратегию вычисления sqrt (f) к чуть более p битам, возвести в квадрат это значение и при необходимости отрегулировать конечный бит на единицу (что часто дешевле).

sqrt isприятно то, что это унарная функция, которая делает исчерпывающее тестирование на одинарную точность осуществимым на стандартном оборудовании.

Функцию OS X soft-float sqrtf можно найти на opensource.apple.com, которая используеталгоритм, описанный выше (я написал, как это происходит).Он лицензирован в соответствии с APSL, который может или не может соответствовать вашим потребностям.

Answer 2

Вероятно (все еще) самая быстрая реализация для поиска обратного квадратного корня и 10 строк кода, которые я обожаю больше всего.

Это основано на приближении Ньютона, но с несколькими причудами. Вокруг этого есть даже замечательная история .

Answer 3

Легче всего реализовать (вы даже можете реализовать это в калькуляторе):

def sqrt(x, TOL=0.000001):
    y=1.0
    while( abs(x/y -y) > TOL ):
        y= (y+x/y)/2.0
    return y

Это точно равно Ньютону Рафсону:

y (новый) = y - f (y) / f '(y)

f (y) = y ^ 2-x и f '(y) = 2y

Подставляя эти значения:

y (новый) = y - (y ^ 2-x) / 2y = (y ^ 2 + x) / 2y = (y + x / y) / 2

Если деление дорого, вы должны учитывать: http://en.wikipedia.org/wiki/Shifting_nth-root_algorithm.

Алгоритмы сдвига:

Предположим, у вас есть два числа a и b, такие, что младшая значащая цифра (равная 1) больше, чем b, а b имеет только один бит, равный (например, a = 1000 и b = 10). Пусть s (b) = log_2 (b) (это просто местоположение бита со значением 1 в b).

Предположим, мы уже знаем значение ^ 2. Теперь (a + b) ^ 2 = a ^ 2 + 2ab + b ^ 2. a ^ 2 уже известно, 2ab: сдвиг a на s (b) +1, b ^ 2: сдвиг b на s (b).

Алгоритм:

Initialize a such that a has only one bit equal to one and a^2<= n < (2*a)^2. 
Let q=s(a).    
b=a
sqra = a*a

For i = q-1 to -10 (or whatever significance you want):
    b=b/2
    sqrab = sqra + 2ab + b^2
    if sqrab > n:
        continue
    sqra = sqrab
    a=a+b

n=612
a=10000 (16)

sqra = 256

Iteration 1:
    b=01000 (8) 
    sqrab = (a+b)^2 = 24^2 = 576
    sqrab < n => a=a+b = 24

Iteration 2:
    b = 4
    sqrab = (a+b)^2 = 28^2 = 784
    sqrab > n => a=a

Iteration 3:
    b = 2
    sqrab = (a+b)^2 = 26^2 = 676
    sqrab > n => a=a

Iteration 4:
    b = 1
    sqrab = (a+b)^2 = 25^2 = 625
    sqrab > n => a=a

Iteration 5:
    b = 0.5
    sqrab = (a+b)^2 = 24.5^2 = 600.25
    sqrab < n => a=a+b = 24.5

Iteration 6:
    b = 0.25
    sqrab = (a+b)^2 = 24.75^2 = 612.5625
    sqrab < n => a=a


Iteration 7:
    b = 0.125
    sqrab = (a+b)^2 = 24.625^2 = 606.390625
    sqrab < n => a=a+b = 24.625

and so on.

Answer 4

Хорошее приближение к квадратному корню в диапазоне [1,4):

def sqrt(x):
  y = x*-0.000267
  y = x*(0.004686+y)
  y = x*(-0.034810+y)
  y = x*(0.144780+y)
  y = x*(-0.387893+y)
  y = x*(0.958108+y)
  return y+0.315413

Нормализуйте число с плавающей запятой, чтобы мантисса находилась в диапазоне [1,4), используйте для него приведенный выше алгоритм, а затем разделите показатель степени на 2. Никаких делений с плавающей запятой нигде.

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