[Неочевидные алгоритмы очевидных вещей] Алгоритм 1. Корень квадратный

Серия постов [Неочевидные алгоритмы очевидных вещей] будет содержать алгоритмы действий, которые кажутся очевидными и простыми, но если задать себе вопрос «как это делается?», то ответ является далеко не очевидным. Разумеется, все эти алгоритмы можно найти в литературе. Под катом располагается алгоритм вычисления корня квадратного числа X.


Корень квадаратный числа Х


Идея


Формируется прямоугольник со сторонами a=1 и b=X. Площадь этого прямоугольника равна S = a * b = X * 1 = X. Преобразовав прямоугольник в квадрат так, что его площадь останется прежней, получим длину стороны, равную корню квадратному из площади фигуры, которая равна Х.
Каждая итерация преобразования прямоугольника в квадрат производится следующим образом:
S = a0 b0;
Создаётся новый четырёхугольник, который имеет одну сторону, равную среднему арифметическому сторон текущего прямоугольника, но площадь такую же:
S = a1 b1, где a1 = (a0+b0) / 2, а b1=S / a1
S = a2 b2, где a2 = (a1+b1) / 2, а b2=S / a2

S = an bn, где an = (an-1+bn-1) / 2, а bn=S / an
И так до тех пор пока не |an-bn| < Eps.

Код функции


#define EPS 1e-10
float my_sqrt(float x){
  float S=x, a=1,b=x;
  while(fabs(a-b)>EPS){  a=(a+b)/2; b = S / a;}
  return (a+b)/2;
}


Ads
AdBlock has stolen the banner, but banners are not teeth — they will be back

More

Comments 25

    +24
    Лучше уж сразу назвать серию постов «Алгоритмы древнего Вавилона»
    • UFO just landed and posted this here
        0
        Хорошо. Ничего критичного для жизни людей не программирует.
      +5
      while(fabs(a-b)<EPS) тут наверное обратное условие,
      преобразуем, пока разница между сторонами больше ε?
        0
        Точно! Очень извиняюсь! Спасибо! Исправил
        +5
        Ожидал увидеть очередную статью про 0x5f3759df.
        Раз уж решили искать корень перебором, почему бы обычную дихотомию не использовать? Или бином Ньютона, на худой конец.
          +7
          По-моему это называется итерационный аналитический алгоритм извлечения квадратного корня на основе формулы Герона.
            +2
            И за реализацией я бы обратился например к этой замечательной статье.
              +1
              А можете объяснить почему стандартный sqrt,

              sqrt13:
              double sqrt13(double n)
              {
                 __asm {
                   fld n
                   fsqrt
                 }
              }   
              


              и sqrt14

              double inline __declspec (naked) __fastcall sqrt14(double n)
              {
              	_asm fld qword ptr [esp+4]
              	_asm fsqrt
              	_asm ret 8
              } 
              


              работают разное время? Каждый следующий раза в два быстрее предыдущего (при этом точность одинаковая).

              Очень уж не хочется верить в то, что в наши дни нужно писать код как в примере 14, который судя по всему ещё не будет работать на x64. Как-то я больше доверял оптимизациям компилятора.
                +1
                может из-за __declspec (naked)?
                  –2
                  Ну по идее обе функции должны бы заинлайниться. Да и пролог/эпилог хороший компилятор должен сам суметь решить когда выкинуть.

                    +4
                    анализировать __asm{} блок компилятору запрещено. именно поэтому 13й вариант имеет все необходимые стандартные заголовки и не имеет права инлайниться — мало ли что ты там в asm{} блоке хочешь делать.

                    кстати, 14й тоже НЕ инлайнится, несмотря на __inline. иначе функция работала бы СОВСЕМ неправильно из-за ret 8 унутре
            +9
            А я ожидал увидеть алгоритмы для компьютеров — быстрые, с маленькой погрешностью, эффективные.
              +2
              Да? codepad.org/m6myk89P
              Автор, ты сам запускал этот код перед тем как постить?
                +1
                У автора после while фигурных скобок не хватает :)
                  –5
                  Прошу прощения=) Исправил=)
                  +3
                  И ещё вдогонку, в приведённой функции после того как поставить скобки после условия в while, всё разно есть ещё одна ошибка,
                  Она связана с неучтённой погрешностью float.

                  Вход 1e18: codepad.org/qK4Bs558 — выполняется
                  Вход 1e20: codepad.org/7Rwyo0ghбесконечный цикл
                  Вход 1e22: codepad.org/Ry9iRGbZ — выполняется
                    –2
                    У меня при компиляции на gcc все представленные входы нормально считаются. Там float какого размера используется?
                      +1
                      4 байта.
                      Проявление этой ошибки зависит от инструкций и их порядка. И будет скорее всего разное при разных настройках компилятора.
                      Можешь попробовать повычислять с числами такого же порядка, в цикле.
                +1
                Если вдруг интересно, то проект Эйлер еще жив.

                Задача чуть повеселее: euler.jakumo.org/problems/view/80.html.
                  0
                  Вот еще в копилку Формула Герона:
                  chipenable.ru/index.php/programming-avr/item/144-sqrt-root.html
                  unsigned int root(unsigned int a)
                  {
                  unsigned int x;
                  x = (a/0x3f + 0x3f)>>1;
                  x = (a/x + x)>>1;
                  x = (a/x + x)>>1;
                  return(x);
                  
                    0
                    а почему не метод Ньютона? тем более есть обобщение на случай n > 2.
                      0
                      вернее этот метод является его частным случаем.

                      Only users with full accounts can post comments. Log in, please.