2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




На страницу Пред.  1, 2
 
 Re: Оценить скорость сходимости формулы Герона
Сообщение27.08.2020, 23:09 
Eule_A, правильно помните, проверил. (Это одна из нескольких книг по которым мы учились. По вычислительным вопросам в курсе начал анализа часто отсылали к этой книге.) Но там вопрос об эффективности выбора начального приближения не обсуждается. Имеет ли такой подход практическое применение?

 
 
 
 Re: Оценить скорость сходимости формулы Герона
Сообщение27.08.2020, 23:36 
Аватара пользователя
GAA
Сейчас была тоже возможность заглянуть в книгу. Там глава про пределы заканчивается приложением как раз про обсуждаемую последовательность. И выбор именно
GAA в сообщении #1481036 писал(а):
В качестве начального приближения можно взять $x_1 = 2^k \left(\frac {2^l} {3} M + 17/24\right).$
там обосновывается (или Вы не об этом?..). А уж применяется ли это сейчас - вот тут ничего не могу сказать. Я от вычислительных вопросов довольно далеко...

 
 
 
 Re: Оценить скорость сходимости формулы Герона
Сообщение27.08.2020, 23:41 
Я это и по памяти помню, и сейчас проверил. Я интересовался практическим нахождением начального приближения.

-- Thu 27.08.2020 22:53:26 --

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

 
 
 
 Re: Оценить скорость сходимости формулы Герона
Сообщение28.08.2020, 10:48 
Всем спасибо!

 
 
 
 Re: Оценить скорость сходимости формулы Герона
Сообщение28.08.2020, 11:47 
Аватара пользователя
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$

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

 
 
 
 Re: Оценить скорость сходимости формулы Герона
Сообщение28.08.2020, 22:35 
В выражении для начального приближения $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. Не так просто вопрос решается. :)

 
 
 [ Сообщений: 21 ]  На страницу Пред.  1, 2


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group