Оценить скорость сходимости формулы Герона

[GAA]

Eule_A, правильно помните, проверил. (Это одна из нескольких книг по которым мы учились. По вычислительным вопросам в курсе начал анализа часто отсылали к этой книге.) Но там вопрос об эффективности выбора начального приближения не обсуждается. Имеет ли такой подход практическое применение?

[Eule_A]

GAA
Сейчас была тоже возможность заглянуть в книгу. Там глава про пределы заканчивается приложением как раз про обсуждаемую последовательность. И выбор именно

GAA в сообщении #1481036 писал(а):

В качестве начального приближения можно взять $x_1 = 2^k \left(\frac {2^l} {3} M + 17/24\right).$

там обосновывается (или Вы не об этом?..). А уж применяется ли это сейчас - вот тут ничего не могу сказать. Я от вычислительных вопросов довольно далеко...

[GAA]

Я это и по памяти помню, и сейчас проверил. Я интересовался практическим нахождением начального приближения.

-- Thu 27.08.2020 22:53:26 --

Другими словами, стоит ли утомлять рассмотрением начального приближения студентов, если это не применяется?

[VoprosT]

Всем спасибо!

[TOTAL]

VoprosT в сообщении #1480977 писал(а):

Необходимо оценить величину $\left\lvert x_n - \sqrt{a}\right\rvert$ в зависимости от $n$.

$x_0=(a+1)/2, \; a>1$
$\left\lvert x_{n+1} - \sqrt{a}\right\rvert < |x_{n+1}-x_n|^2/2$

Метод с квадратичной сходимостью, поэтому трудиться над начальным приближением нет смысла

[GAA]

В выражении для начального приближения $x_1 = 2^k \left(\frac {2^l} {3} M + 17/24\right)$ величина $M$ — это мантисса, $2k+1$ — несмещённый порядок. Кажется, что можно быстро построить начальное приближение $x_1$. Аппаратно массовыми платформами поддерживаются binary32 и binary64. Если предположить, что уже реализованы арифметические операции над 128 битными, то можно пробовать реализовать извлечение корня. Под рукой у меня компилятора с поддержкой binary128 сейчас нет, поэтому я попробовал набросать черновик извлечения корня для 32 битных, но не использовать специфические инструкции fxtract, fscale. Остальные инструкции предположим макросы. Представление в памяти всех binary подобно.
Функция с начальным приближением из книги Ильина и Позняка

Используется синтаксис Delphi

 function Hero(a: Single): Single;
const
 c1div3  : Double = 1/3;
 c17div24: Double = 17/24;
 cOneHalf: Double = 0.5;
var
 x1:     Single;
 count : Word;
begin
 asm // Определение x1 - Result
// Получение k, l и M
  mov eax, dword ptr a
  mov ecx, eax
  and ecx, 007FFFFFh    // 11111111111111111111111 Обнуляем порядок
  or  ecx, 3F800000h    // 0011111110000000000000000000000b Задаём порядок = 127
  mov dword ptr x1, ecx // ecx = M, Занесли в x1 мантиссу
  shr eax, 23
  sub ax,  127          // Вычитаем из порядка смещение
  shr ax,  1            // делим на 2: al = k, CF = l
// Вычисление 1/3*2^l*M+17/24
  fld x1                // Мантиссу загружаем в вершину стека
  jae @@                // Eсли l=0, то удваивать мантиссу не надо
  fadd st, st           // Если l=1, то удваиваем мантиссу st(0) = 2*M
@@:
  fmul c1div3           // st(0) = 1/3*2^l*M
  fadd c17div24         // st(0) = 1/3*2^l*M+17/24
  fstp Result           // Сохраняем 1/3*2^l*M+17/24 в Result
// Умножение (1/3*2^l*M+17/24) на 2^k
  mov ecx, dword ptr Result
  mov edx, ecx
  shr edx, 23
  add dx,  ax
  shl edx, 23
  and ecx, 007FFFFFh    // Обнуляем порядок
  or  ecx, edx
  mov dword ptr Result, ecx // Сохраняем x1 в Result для последующего уточнения
 end;
// Итерации до прекращения изменения
  repeat
   x1    := Result;
   Result:= cOneHalf*(x1+a/x1);
  until Result = x1;
end;

Код особо не вылизывался. Возможно, можно ускорить, но что-то не верится в значительное ускорение. Или я где-то крепко ошибся.

В качестве начального приближения можно взять «Быстрый обратный квадратный корень», см. сообщение post1480711.html#p1480711 photon и далее. И затем применить те же итерации для уточнения корня.

Используется синтаксис Delphi

function FastSqrt(a: Single):Single;
const
  cOneHalf: Double = 0.5;
var
 x: Single;
 count : Word;
begin
 asm
  mov eax, dword ptr a
  shr eax, 1
  neg eax
  add eax, 5F3759DFh
  mov dword ptr x, eax
 end;
 Result:= a*x;
 repeat
  x := Result;
  Result:= cOneHalf*(x+a/x);
 until Result = x;
end;

Код значительно короче, а начальное приближение для последующих уточнений методом Ньютона для опробованных мною нескольких значений оказывалось лучше (требовалось на одну итерацию метода Ньютона меньше). Сравнив два кода, понял, что подсчитывать старательно времена нет особого смысла. Начальное приближение из книги Ильина и Позняка имеет простое и ясное обоснование, но, увы, для вычислений мало пригодно.
Выборочное сравнение времён выполнения

Используется синтаксис Delphi

program fsqrt;
uses sysutils;

const
 RepN = 100000000;
var
 Time1, Time2: TDateTime;

{$I sqrt.inc}

function TOTAL(a: Single):Single;
const
  cOneHalf: Double = 0.5;
var
 x: Single;
begin
  Result:= (a+1)/2;
  repeat
   x := Result;
   Result:= cOneHalf*(x+a/x);
  until Result = x;
end;

function TOTALRep: TDateTime;
var
 Rep: longint;
begin
  Time1:= Time;
  for Rep:= 1 to RepN
   do begin
       TOTAL(Rep);
       TOTAL(Rep);
       TOTAL(Rep);
       TOTAL(Rep);
       TOTAL(Rep);
       TOTAL(Rep);
       TOTAL(Rep);
       TOTAL(Rep);
       TOTAL(Rep);
       TOTAL(Rep);
      end;
  Time2:= Time;
  Result:= Time2-Time1;

end;

function HeroRep: TDateTime;
var
 Rep: longint;
begin
  Time1:= Time;
  for Rep:= 1 to RepN
   do begin
       Hero(Rep);
       Hero(Rep);
       Hero(Rep);
       Hero(Rep);
       Hero(Rep);
       Hero(Rep);
       Hero(Rep);
       Hero(Rep);
       Hero(Rep);
       Hero(Rep);
      end;
  Time2:= Time;
  Result:= Time2-Time1;
end;

function FastSqrtRep: TDateTime;
var
 Rep: longint;
begin
  Time1:= Time;
  for Rep:= 1 to RepN
   do begin
       FastSqrt(Rep);
       FastSqrt(Rep);
       FastSqrt(Rep);
       FastSqrt(Rep);
       FastSqrt(Rep);
       FastSqrt(Rep);
       FastSqrt(Rep);
       FastSqrt(Rep);
       FastSqrt(Rep);
       FastSqrt(Rep);
      end;
  Time2:= Time;
  Result:= Time2-Time1;
end;
var
 d1, d2, d3: TDateTime;
begin
  d1:= HeroRep; d2:= FastSqrtRep; d3:= TOTALRep;
  writeln(d1, ' ', d2, ' ', d1/d2, d3/d2);
end.
 

Функция с начальным приближением из Ильина и Позняка приблизительно в 1.1 дольше выполняется по сравнению с начальным приближением Fast…, а функция с начальным приближением TOTAL приблизительно в 5.6 раз дольше Fast….

Погуглив, нашёл кучу работ, содержащих, в том числе, и вычисление корня для binary128. Не так просто вопрос решается. :)