Тогда следует повторить исходный вопрос
1. Квадратный корень вычисляется в большинстве библиотек при помощи метода Ньютона для обратного корня, чтобы избежать деления. (Выше об этом много писали). Если начальное приближение очень точное или достаточно не очень точного результата, то нужно выполнить 2-3 итерации. Поэтому цикл «разворачивается»[см. ниже пример TI]. (В некоторых библиотеках используется очень точное начальное приближение. В этом случае итераций может быть одна или две и вместо использования метода Ньютона для обратного корня используется метод Ньютона для самого корня.)
В основном отличаются подходы в получении начального приближения. Эти подходы делятся на те, что явно сводят значение аргумента к некоторому диапазону (редукция аргумента) и те, что не сводят. После редукции аргумента используют различные варианты поиска по таблице (например Ken Turkowski “Computing the Inverse Square Root”, Apple Technical Report No. 95, 1994) или интерполяции полиномами, примером может служить реализация cephes (из архива single.tgz — для digits32, т.е. для одинарной точности)
/* sqrtf.c
*
* Square root
*
*
*
* SYNOPSIS:
*
* float x, y, sqrtf();
*
* y = sqrtf( x );
*
*
*
* DESCRIPTION:
*
* Returns the square root of x.
*
* Range reduction involves isolating the power of two of the
* argument and using a polynomial approximation to obtain
* a rough value for the square root. Then Heron's iteration
* is used three times to converge to an accurate value.
*
*
*
* ACCURACY:
*
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,1.e38 100000 8.7e-8 2.9e-8
*
*
* ERROR MESSAGES:
*
* message condition value returned
* sqrtf domain x < 0 0.0
*
*/
_
/*
Cephes Math Library Release 2.2: June, 1992
Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* Single precision square root
* test interval: [sqrt(2)/2, sqrt(2)]
* trials: 30000
* peak relative error: 8.8e-8
* rms relative error: 3.3e-8
*
* test interval: [0.01, 100.0]
* trials: 50000
* peak relative error: 8.7e-8
* rms relative error: 3.3e-8
*
* Copyright (C) 1989 by Stephen L. Moshier. All rights reserved.
*/
#include "mconf.h"
#ifdef ANSIC
float frexpf( float, int * );
float ldexpf( float, int );
float sqrtf( float xx )
#else
float frexpf(), ldexpf();
float sqrtf(xx)
float xx;
#endif
{
float f, x, y;
int e;
f = xx;
if( f <= 0.0 )
{
if( f < 0.0 )
mtherr( "sqrtf", DOMAIN );
return( 0.0 );
}
x = frexpf( f, &e ); /* f = x * 2**e, 0.5 <= x < 1.0 */
/* If power of 2 is odd, double x and decrement the power of 2. */
if( e & 1 )
{
x = x + x;
e -= 1;
}
e >>= 1; /* The power of 2 of the square root. */
if( x > 1.41421356237 )
{
/* x is between sqrt(2) and 2. */
x = x - 2.0;
y =
((((( -9.8843065718E-4 * x
+ 7.9479950957E-4) * x
- 3.5890535377E-3) * x
+ 1.1028809744E-2) * x
- 4.4195203560E-2) * x
+ 3.5355338194E-1) * x
+ 1.41421356237E0;
goto sqdon;
}
if( x > 0.707106781187 )
{
/* x is between sqrt(2)/2 and sqrt(2). */
x = x - 1.0;
y =
((((( 1.35199291026E-2 * x
- 2.26657767832E-2) * x
+ 2.78720776889E-2) * x
- 3.89582788321E-2) * x
+ 6.24811144548E-2) * x
- 1.25001503933E-1) * x * x
+ 0.5 * x
+ 1.0;
goto sqdon;
}
/* x is between 0.5 and sqrt(2)/2. */
x = x - 0.5;
y =
((((( -3.9495006054E-1 * x
+ 5.1743034569E-1) * x
- 4.3214437330E-1) * x
+ 3.5310730460E-1) * x
- 3.5354581892E-1) * x
+ 7.0710676017E-1) * x
+ 7.07106781187E-1;
sqdon:
y = ldexpf( y, e ); /* y = y * 2**e */
return( y);
}
Некоторые процессоры поддерживают только вычисление начального приближения для обратного квадратного корня. См., например, TMS320C28x Floating Point Unit and Instruction Set Reference Guide, 2008 2007 Revised 2008
pdf, см. с. 50. Обратный квадратный корень находится приближенно при помощи итераций Ньютона
; Y = sqrt(X)
; Ye = Estimate(1/sqrt(X))
; Ye = Ye*(1.5 - Ye*Ye*X*0.5)
; Ye = Ye*(1.5 - Ye*Ye*X*0.5)
; Y = X*Ye
_sqrt: ; R0H = X on entry
EISQRTF32 R1H, R0H ; R1H = Ye = Estimate(1/sqrt(X))
MPYF32 R2H, R0H, #0.5 ; R2H = X*0.5
MPYF32 R3H, R1H, R1H ; R3H = Ye*Ye
NOP
MPYF32 R3H, R3H, R2H ; R3H = Ye*Ye*X*0.5
NOP
SUBF32 R3H, #1.5, R3H ; R3H = 1.5 - Ye*Ye*X*0.5
NOP
MPYF32 R1H, R1H, R3H ; R2H = Ye = Ye*(1.5 - Ye*Ye*X*0.5)
NOP
MPYF32 R3H, R1H, R2H ; R3H = Ye*X*0.5
NOP
MPYF32 R3H, R1H, R3H ; R3H = Ye*Ye*X*0.5
NOP
SUBF32 R3H, #1.5, R3H ; R3H = 1.5 - Ye*Ye*X*0.5
CMPF32 R0H, #0.0 ; Check if X == 0
MPYF32 R1H, R1H, R3H ; R2H = Ye = Ye*(1.5 - Ye*Ye*X*0.5)
NOP
MOV32 R1H, R0H, EQ ; If X is zero, change the Ye estimate to 0
MPYF32 R0H, R0H, R1H ; R0H = Y = X*Ye = sqrt(X)
LRETR
Для тестирования функций / аппаратных реализаций иногда используется сведение к возведению в степень, например Камкин А. C., Чупилко М. M. Тестирование модулей арифметики с плавающей точкой микропроцессоров на соответствие стандарту IEEE 754, 2008
Цитата:
4.1.2. Эталонная реализация извлечения квадратного корня
Извлечение квадратного корня реализовано следующим образом. Знак результата совпадает со знаком операнда (отрицательный знак операнда возможен только в случае, когда его значение равно –0), порядок равен половине порядка операнда6. Мантисса результата вычисляется с помощью старшего к младшему биту мантиссы производятся следующие действия:
* текущий бит мантиссы результата устанавливается в единицу;
* полученная мантисса возводится в квадрат и сравнивается с исходным числом;
* если исходное число меньше, установленная единица сбрасывается в нуль.
2. Многие функции вычисляются при помощи многочленов Чебышева, а для менее точных вычислений используются Паде аппроксимации. Вычисление экспоненты digits64 самыми популярными методами я попробовал рассмотреть в теме
«Быстрое вычисление функции exp(x)». В частности, в связи с таблицами см. указанную в той теме статью Tang P.T.P. Table-Driven Implementation of the Exponential Function in IEEE Floating-Point Arithmetic.
для итерационной формулы Ньютона в виде: 
) нам достаточно массива длиной 
с элементами типа binary32 (Single): 131072 байт или 128 кбайт. Это большая таблица. Зачем вместо этой таблицы делать две?
или не менее 29 верных битов (или даже 30 если допустить ошибку в обе стороны, а не только с недостачей).