Статья любезно предоставлена
Copyright © Nikitine Valeri F. 2000

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

Для вычисления квадратного корня здесь использован известный метод Ньютона. Рассмотрены машинно-зависимый и машинно-независимый варианты. Перед применением алгоритма Ньютона, область определения исходного числа сужается до [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) представляет быстро сходящуюся (при хорошем начальном приближении) серию итераций: ai+1=0.5*(ai+x/ai), где i -- номер итерации.

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

/* Square root almost without mathematic library. Optimized for floating
   point single precision. The functions frexp() and ldexp() from the
   mathematic library are very fast and machine-dependent. 
   The second variant is less fast but does not use the math library 
   at all.
   Copyright (c) Nikitin V.F. 2000

   Square root decomposing the argument into the mantisse and the 
   exponent (faster program):
   float Sqroot(float x);

   Square root without usage of any external functions (slower but 
   machine-independent):
   float Sqroot1(float x);

   In case of invalid domain (x<0.) functions return zero (do not 
   generate an error as library equivalents). 
*/

#include <math.h>   /* frexp() and ldexp() only */

/* 4 iterations needed for the single precision */
#define ITNUM 4

/* variant using external f.p. decomposition/combining to [0.5,1] */
float Sqroot(float x) {
  int expo,i;
  float a,b;
  if(x<=0.F) return(0.F);
  /* decompose x into mantisse ranged [0.5,1) and exponent. Machine-
     dependent operation is presented here as a function call. */
  x=frexp(x,&expo);
  /* odd exponent: multiply mantisse by 2 and decrease exponent 
     making it even. Now the mantisse is ranged [0.5,2.) */
  if(expo&1) {x*=2.F;expo--;}
  /* initial approximation */
  a=1.F;
  /* process ITNUM Newtonian iterations with the mantisse as 
     an argument */
  for(i=ITNUM;i>0;i--) {
    b=x/a; a+=b; a*=0.5F;
  }
  /* divide the exponent by 2 and combine the result.
     Function ldexp() is opposite to frexp. */
  a=ldexp(a,expo/2);
  return(a);
}

/* variant without math lib usage. The domain is shrunk to [1,16]. 
   Repeated divisions by 16 are used. */
float Sqroot1(float x) {
  int sp=0,i,inv=0;
  float a,b;
  if(x<=0.F) return(0.F);
  /* argument less than 1 : invert it */
  if(x<1.F) {x=1.F/x;inv=1;}
  /* process series of division by 16 until argument is <16 */
  while(x>16.F) {sp++;x/=16.F;}
  /* initial approximation */
  a=2.F;
  /* Newtonian algorithm */
  for(i=ITNUM;i>0;i--) {
    b=x/a; a+=b; a*=0.5F;
  }
  /* multiply result by 4 : as much times as divisions by 16 took place */
  while(sp>0) {sp--;a*=4.F;}
  /* invert result for inverted argument */
  if(inv) a=1.F/a;
  return(a);
}