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

Для вычисления квадратного корня в этом примере использован метод Ньютона. Рассмотрены машинно-зависимый и машинно-независимый варианты. Перед применением алгоритма Ньютона область определения исходного числа сужается до [0.5,2] (во втором варианте до [1,16]). Второй вариант машинно-независим, но работает дольше.

Вообще, это не самый быстрый вариант, но один из самых быстрых. Основная проблема заключается в том, что нет универсального машинного представления чисел с плавающей точкой, поэтому разделение числа на мантиссу и двоичную экспоненту как составляющих его компьютерного представления нельзя записать единым образом. Имено поэтому подключено описание математической библиотеки, из которой, впрочем, используются только frexp() и ldexp(). Конкретная реализация этих функций очень проста, но машинно-зависима. При разделении числа на мантиссу и экспоненту, мантисса оказывается в пределах [0.5,1). Если при этом экспонента нечетная, мантисса умножается на 2, тем самым область определения становится [0.5,2). К мантиссе применяется алгоритм Ньютона с начальным приближением 1. Окончательный результат получается с помощью применения ldexp() с половинным значением экспоненты.

Для машинно-независимого варианта выделение экспоненты заменяется серией последовательных делений исходного значения на 16, пока аргумент не станет определяться на интервале [1,16]. Если исходный аргумент был меньше 1, то перед серией делений обращаем его. Алгоритм Ньютона применяется с начальным приближением 2. Окончательный результат получается серией последовательных умножений на 4 и дальнейшим обращением, если аргумент обращался.

Сам алгоритм Ньютона для вычисления a=Sqroot(x) представляет быстро сходящуюся (при хорошем начальном приближении) серию итераций: a_i+1=0.5*(a_i+x/a_i), где i -- номер итерации.

 

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

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

/* Вычисление квадратного корня почти без использования мат. библиотеки :).
   Оптимизировано для чисел с плавающей точкой одинарной точности.
   Фукнции 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);
}