Вычисление квадратного корня по алгоритму Ньютона
Страница 2. Программы вычисления квадратного корня

Программы вычисления квадратного корня 

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

/* Вычисление квадратного корня почти без использования мат. библиотеки :).
 Оптимизировано для чисел с плавающей точкой одинарной точности.
 Фукнции frexp() и ldexp() из мат. библиотеки быстры и машинно-зависимы. 
 Второй вариант менее быстр, зато не использует мат. библиотеку.

 Разложение квадратного корня аргумента в мантиссу и экспонент (быстрый алгоритм):
 float Sqroot(float x);

 Вычисление квадратного корня без использования внешних фукнций (медленнее, но машинно-независим):
 float Sqroot1(float x);

 В случае некорректного промежутка (x<0.) фукнции возвращают 0 и не генерят ошибку. 
*/

#include <math.h> /* только frexp() и ldexp()*/

/* для одинарной точности нужны 4 итерации */
#define ITNUM 4

/* вариант с использованием внешних фукнций разложения/объединения на [0.5,1] */
float Sqroot(float x) {
 int expo,i;
 float a,b;
 if(x<=0.F) return(0.F);
 /* разложение x в мантиссу на промежутке [0.5,1) и экспонент. 
 Машинно-зависимые операции представлены вызовами функций. */
 x=frexp(x,&expo);
 /* нечетный экспонент: умножаем мантиссу на 2 и уменьшаем экспонент, делая его четным.
 Теперь мантисса в промежутке [0.5,2.) */
 if(expo&1) {x*=2.F;expo--;}
 /* начальное приближение */
 a=1.F;
 for(i=ITNUM;i>0;i--) {
 b=x/a; a+=b; a*=0.5F;
 }
 /* делим экспонент на 2 и объединяем результат.
 Фукнция ldexp() противоположна frexp. */
 a=ldexp(a,expo/2);
 return(a);
}

/* Вариант без использования библиотек. Промежуток уменьшен до [1,16]. 
 Используется 16 повторяющихся делений. */
float Sqroot1(float x) {
 int sp=0,i,inv=0;
 float a,b;
 if(x<=0.F) return(0.F);
 /* аргумент меньше 1 : инвертируем его */
 if(x<1.F) {x=1.F/x;inv=1;}
 /* последовательно делим на 16 пока аргумент не станет <16 */
 while(x>16.F) {sp++;x/=16.F;}
 /* начальное приближение */
 a=2.F;
 /* Алгоритм Ньютона */
 for(i=ITNUM;i>0;i--) {
 b=x/a; a+=b; a*=0.5F;
 }
 while(sp>0) {sp--;a*=4.F;}
 /* инвертируем результат для инвертированнго аргумента */
 if(inv) a=1.F/a;
 return(a);
}