Нахождение квадратного корня через двоичное разложение

Еще одна оптимизация алгоритма

Следующее улучшение алгоритма

Вычисление 2* k + 1 реализуется быстрее, чем (k + 1)^2, В свою очередь, вычисление 2* k + 1 можно ускорить, используя простой факт, что его значение увеличивается на 2 на каждом шаге. Рассмотрим спецификацию:

sq2(nat x, k, n, d: nat m)

pre k^2 <= x & n = k^2 & d = 2 * k + 1

post m^2 <= x < (m + 1)^2;

Дополнительный параметр d обеспечивает возможность ускоренного вычисления 2 * k + 1. Соответствующая программа приведена ниже:

isqrt(nat x: nat m)

{ sq2(x, 0, 0, 1: m) }

post m^2 <= x < (m + 1)^2;

sq2(nat x, k, n, d: nat m)

pre k^2 <= x & n = k^2 & d = 2 * k + 1

{ nat p = n + d;

if (x < p) m = k else sq2(x, k + 1, p, d + 2: m)

}

post m^2 <= x < (m + 1)^2;

Приведем снова алгоритм sq2.

sq2(nat x, k, n, d: nat m)

pre k^2 <= x & n = k^2 & d = 2 * k + 1

{ nat p = n + d;

if (x < p) m = k else sq2(x, k + 1, p, d + 2: m)

}

Здесь p используется только в операции сравнения x < p, которая эквивалентна x – p < 0. Мы можем пересчитывать на каждом шаге значение
x – n. В этом случае переменные p и n будут не нужны. Определим спецификацию:

sq3(nat x, y, k, d: nat m)

pre k^2 <= x & d = 2 * k + 1 & y = x - k^2

post m^2 <= x < (m + 1)^2;

Программа получается соответствующими модификациями программы sq2.

isqrt(nat x: nat m)

{ sq3(x, x, 0, 1: m) }

post m^2 <= x < (m + 1)^2;

sq3(nat x, y, k, d: nat m)

pre k^2 <= x & d = 2 * k + 1 & y = x - k^2

{ if (y < d) m = k else sq3(x, y – d, k + 1, d + 2: m) }

post m^2 <= x < (m + 1)^2;

Построим императивную программу для данной предикатной программы. На первом этапе трансформации приведенной предикатной программы в эффективную императивную реализуется следующие склеивания переменных.

sq3: m <- k; x <- y;

В результате склеивания получим:

sq3(nat x, x, m, d: nat m)

{ if (x < d) m = m else sq3(x, x – d, m + 1, d + 2: m) }

На втором этапе реализуется замена хвостовой рекурсии циклами. Результатом проведения двух этапов трансформации является следующая программа:

isqrt(nat x: nat m)

{ sq3(x, x, 0, 1: m) }

sq3(nat x, x, m, d: nat m)

{

for (;;) if (x < d) breakelse { x = x – d; m = m + 1; d = d + 2 }

}

На третьем этапе проведем подстановку телa определения sq3 на место вызова:

isqrt(nat x: nat m)

{ m = 0; nat d = 1;

for (;;) if (x < d) breakelse { x = x – d; m = m + 1; d = d + 2 }

}

Наконец, реализуем вынесение оператора x = x – d перед условным оператором:

isqrt(nat x: nat m)

{ m = 0; nat d = 1;

for (;;) { x = x – d; if (x < 0) break; m = m + 1; d = d + 2 }

}

Последовательная оптимизация простого переборного алгоритма sq0 привела к алгоритму sq3. Все рассмотренные версии этого алгоритма имеют серьезный недостаток: число шагов алгоритма равно значению квадратного корня. Существуют алгоритмы, в которых число шагов имеет логарифмическую оценку. Один из них базируется на двоичном разложении результата, и на каждом шаге находит очередную двоичную цифру квадратного корня. Предположим, что аргумент n ограничен значением 2^2*p. Представим исходную задачу спецификацией:

isqrt(nat n, p: nat s)

pre p > 0 & n < 2^(2*p)

post s^2 <= n < (s + 1)^2;

Результат s будет иметь не более p двоичных цифр. Значение старшей цифры равно 0, если n < 2^(2*(p-1)) и 1 в противном случае. Допустим, мы нашли для результата s очередное приближение q в виде старших двоичных цифр с номерами от p до k. Тогда должно выполняться условие q^2 <= n < (q + 2^k)^2. Представим это спецификацией:

sq4(nat n, p, k, q: nat s)

pre p > 0 & n < 2^(2*p) & k <= p & q^2 <= n < (q + 2^k)^2

post s^2 <= n < (s + 1)^2;

Очевидно определение для isqrt.

isqrt(nat n, p: nat s)

pre p > 0 & n < 2^(2*p)

{ sq4(n, p, p, 0, n: s)}

post s^2 <= n < (s + 1)^2;

Спецификацию sq4 применим для вычисления цифры с номером k-1. Получим определение:

sq4(nat n, p, q, k: nat s)

pre p >= 0 & n < 2^(2*p) & k <= p & q^2 <= n < (q + 2^k)^2

{ if (k = 0) s = q

elseif (n < (q + 2^(k-1))^2) sq4(n, p, k - 1, q: s)

else sq4(n, p, k - 1, q + 2^(k-1): s)

}

post s^2 <= n < (s + 1)^2;

measure k;

В соответствии с данным алгоритмом в зависимости от условия следующим приближением для s является q или q + 2^(k-1). Чтобы упростить вычисление условия n < (q + 2^(k-1))^2, реализуется серия оптимизаций, аналогичная переходу от sq1 к sq3. Имеет место

(q + 2^(k-1))^2 = q^2 + 2^((k-1)*2) + q * 2^k.

Тогда условие n < (q + 2^(k-1))^2 эквивалентно:

n - q^2 < 2^((k-1)*2) + q * 2^k.

Обозначим y = n - q^2. Преобразуем определение sq4 к следующему виду:

sq4(nat n, p, k, q, y: nat s)

pre p >= 0 & n < 2^(2*p) & k <= p & q^2 <= n < (q + 2^k)^2 & y = n - q^2

{ if (k = 0) s = q

else { nat t = 2^((k-1)*2) + q * 2^k;

if (y < t) sq4(n, p, k - 1, q, y: s)

else sq4(n, p, k - 1, q or 2^(k-1), y – t: s)

}

}

post s^2 <= n < (s + 1)^2;

measure k;

Кроме того, мы заменили q + 2^(k-1) эквивалентным q or 2^(k-1), вычисляемым быстрее, где or ¾ побитовая операция над натуральными.

Построим императивную программу для данной предикатной программы. На первом этапе трансформации приведенной предикатной программы в эффективную императивную реализуется следующие склеивания переменных:

sq4: n <- y; s <- q; k <- p.

В результате склеивания получим:

isqrt(nat n, p: nat s)

{ sq4(n, p, p, 0, n: s)}

sq4(nat n, k, k, s, n: nat s)

{ if (k = 0) s = s

else { nat t = 2^((k-1)*2) + s * 2^k;

if (n < t) sq4(n, k, k - 1, s, n: s)

else sq4(n, k, k - 1, s or 2^(k-1), n – t: s)

}

}

На втором этапе реализуется замена хвостовой рекурсии циклами. Результатом проведения двух этапов трансформации является следующая программа:

isqrt(nat n, p: nat s)

{ sq4(n, p, p, 0, n: s) }

sq4(nat n, k, k, s, n: nat s)

{ for (;;) {

if (k = 0) break;

nat t = 2^((k-1)*2) + s * 2^k;

if (n < t) k = k - 1

else { n = n – t; s = s or 2^(k-1); k = k – 1 }

}

}

На третьем этапе проведем подстановку телa определения sq4 на место вызова и упрощения. Получим программу:

nat n, p;

nat s = 0;

for (nat k = p; k = 0; k = k - 1) {

nat t = 2^((k-1)*2) + s * 2^k;

if (n >=t) { n = n – t; s = s or 2^(k-1); }

}

Обозначим:

sq(k) = { t = 2^((k-1)*2) + s * 2^k; if (n >=t) { n = n – t; s = s or 2^(k-1); } }

Проведем развертку цикла. Получим

nat n, p;

nat s = 0;

sq(p); sq(p – 1); …; sq(2); sq(1)

Проведем подстановку тела для sq(p) и sq(1). Получим итоговую программу.

sq(k) = { t = 2^((k-1)*2) + s * 2^k; if (n >=t) { n = n – t; s = s or 2^(k-1); } }

nat n, p;

nat s = 0;

if (n >= 2^((p-1)*2)) { n = n – 2^((p-1)*2); s = 2^(p-1); }

sq(p – 1); …; sq(2);

if (n >= s * 2) s = s or 1

По сравнению с программой, представленной на сайте

http://www.finesse.demon.co.uk/steven/sqrt.html,

параметр k смещен на единицу. Однако поскольку тело sq(k) подставляется открыто на место вызовов sq(p – 1), …, sq(2), то обе программы идентичны.

На примере данной задачи видно, что проблема автоматизированной трансформации предикатной программы в эффективную императивную программу далеко не проста.

Дата добавления: 2015-06-27