Алгоритм вычисления логарифма сдвигами по степеням двойки

Для ускорения вычисления целочисленного двоичного логарифма вместо сдвига по одному биту можно использовать более длинные сдвиги. Допустим, требуется вычислить ilog2(a) для a < 2^(2^p) и a >= 2^(2^(p-1)), где p > 1. Тогда в соответствии со свойствами логарифма

log2(a) = log2(a / 2^(2^(p-1))) + log2(2^(2^(p-1))) =

= log2(a / 2^(2^(p-1))) + 2^(p-1).

Поскольку a / 2^(2^(p-1)) < 2^(2^(p-1)), то приведенное тождество определяет схему вычисления двоичного логарифма. Можно показать, что при вычислении ilog2 вещественное деление можно заменить целочисленным. Это дает возможность представить алгоритм в виде следующего определения предиката:

ilogB(nat a, p: nat r) pre a >= 1 & a < 2^(2^p)

{ if (p = 0) r = 0

elseif (a >= 2^(2^(p-1))) r = ilogB(div(a, 2^(2^(p-1))), p - 1) + 2^(p-1)

else ilogB(a, p - 1: r)

} post r = floor(log2(a));

Рекурсия в данном определении не является хвостовой. Введем дополнительный параметр ¾ накопитель m. Предикатная программа с хвостовой рекурсией представлена следующими двумя определениями:

ilogB(nat a, p: nat r) pre a >= 1 & a < 2^(2^p)

{ lgB(a, p, 0: r) } post r = floor(log2(a));

lgB(nat a, p, m: nat r) pre a >= 1 & a < 2^(2^p)

{ if (p = 0) r = m

elseif (a >= 2^(2^(p-1))) lgB(div(a, 2^(2^(p-1))), p - 1, m + 2^(p-1): r)

else lgB(a, p - 1, m: r)

} post r = floor(log2(a)) + m;

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

lgB: r <- m;

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

ilogB(nat a, p: nat r)

{ lgB(a, p, 0: r) }

lgB(nat a, p, r: nat r)

{ if p = 0 then r = r

elseif (a >= 2^(2^(p-1))) lgB(div(a, 2^(2^(p-1))), p - 1, r + 2^(p-1): r)

else lgB(a, p - 1, r: r)

}

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

ilogB(nat a, p: nat r)

{ r = 0;

for (; p!= 0; p = p - 1)

if (a >= 2^(2^(p-1))) {a = div(a, 2^(2^(p-1))); r = r + 2^(p-1)}

}

Для p = 5, что соответствует 32-битным натуральным числам, проведем развертку цикла. Получим программу:

r = 0;

if (a >= 2^16) {a = div(a, 2^16); r = r + 16}

if (a >= 2^8) {a = div(a, 2^8); r = r + 8}

if (a >= 2^4) {a = div(a,2^4); r = r + 4}

if (a >= 2^2) {a = div(a, 2^2); r = r + 2}

if (a >= 2^1) {a = div(a, 2^1); r = r + 1}

Наконец, используя операции сдвига, получим окончательную программу:

r = 0;

if (a >= 1<<16) {a = a >>16; r = r + 16}

if (a >= 1<< 8) {a = a >>8); r = r + 8}

if (a >= 1<< 4) {a = a >>4; r = r + 4}

if (a >= 1<< 2) {a = a >>2; r = r + 2}

if (a >= 1<< 1) { r = r + 1}

Вернемся к предикатной программе и рассмотрим ее верификацию относительно спецификаций.

ilogB(nat a, p: nat r) pre a >= 1 & a < 2^(2^p)

{ lgB(a, p, 0: r) } post r = floor(log(2, a));

lgB(nat a, p, m: nat r) pre a >= 1 & a < 2^(2^p)

{ if p = 0 then r = m

elseif a >= 2^(2^(p-1)) then lgB(div(a, 2^(2^(p-1))), p - 1, m + 2^(p-1): r)

else lgB(a, p - 1, m: r)

} post r = floor(log(2, a)) + m;

measure e(p) = p

Представим теорию для определения предиката ilogB.

ilogB: THEORYBEGINIMPORTING lgb, ilog2a: VAR posnatr, p: VAR natP_ilogB(a, p): bool = a < 2^(2^p)Q_ilogB(a, r): bool = r = ilog2(a))satisfy_spec: LEMMA P_ilogB(a, p) IMPLIES EXISTS r: Q_ilogB(a, r)

Целью является доказательство истинности правила:

P_ilogB(a, p) & Q_ilogB(a, r) ├ LS(lgB(a, p, 0: r)).

В соответствии с правилом FB1 генерируется лемма:

FB1: LEMMA

P_ilogB(a, p) & Q_ilogB(a, r) IMPLIES P_lgB(a, p) & Q_lgB(a, p, 0, r).

Эта лемма завершает теорию ilogB. Представим начальную часть теории lgB.

lgB: THEORYBEGINp, m, r: VAR nata: VAR posnatP_lgB(a, p): bool = a < 2^(2^p)Q_lgB(a, m, r): bool = r = m + ilog2(a)satisfy_spec: LEMMA P_lgB(a, p) IMPLIES EXISTS r: Q_lgB(a, m, r)e(p): nat = p

Целью является доказательство истинности правила

FR1: Induct(e) & P_lgB(a, p) & Q_lgB(a, m, r)├ LS(< тело lgB >).

В соответствии с правилом FC1 для первой альтернативы условного оператора генерируется лемма:

FC1: LEMMA

P_lgB(a, p) & Q_lgB(a, m, r) & p = 0 IMPLIES r = m.

В соответствии с правилом FC2 для второй альтернативы условного оператора генерируется цель:

FC2: Induct(e) & P_lgB(a, p) & Q_lgB(a, m, r) & NOT p = 0 ├

LS(<внутренний условный оператор>).

В соответствии с правилом FC1 для первой альтернативы условного оператора генерируется цель:

FC1: Induct(e) & P_lgB(a, p) & Q_lgB(a, m, r) & NOT p = 0 & a >= 2^(2^(p-1))

├ LS(lgB(div(a, 2^(2^(p-1))), p - 1, m + 2^(p-1): r)).

Для первого рекурсивного вызова lgB по правилу FB3 генерируется лемма:

FB3: LEMMA P_lgB(a, p) & Q_lgB(a, m, r) & NOT p = 0 & a >= 2^(2^(p-1))

IMPLIES e(p - 1) < e(p) & P_lgB(div(a, 2^(2^(p-1))), p - 1) &

Q_lgB(div(a, 2^(2^(p-1))), m + 2^(p-1), r).

В соответствии с правилом FC2 для второй альтернативы условного оператора генерируется цель:

FC2: Induct(e) & P_lgB(a, p) & Q_lgB(a, m, r) &

NOT p = 0 & NOT a >= 2^(2^(p-1)) ├ LS(lgB(a, p - 1, m: r)).

Для второго рекурсивного вызова lgB по правилу FB3 генерируется лемма:

FB3_1: LEMMA P_lgB(a, p) & Q_lgB(a, m, r) & NOT p = 0 & NOT a >= 2^(2^(p-1))

IMPLIES e(p - 1) < e(p) & P_lgB(a, p - 1) & Q_lgB(a, m, r).

Автоматическая верификация приведенных условий самих корректности не потребовала много усилий. Однако для их доказательства потребовалось построить теории log, log2, ilog2 и доказать 32 леммы в них. Вот эти теории:

log: THEORY

BEGIN

p: VAR {y: nat | y > 1}

IMPORTING lnexp@ln_exp

x, y: VAR posreal

n: VAR posnat

i: VAR int

log(p, x): real = ln(x) / ln(p)

log_1: LEMMA log(p, 1) = 0

log_p: LEMMA log(p, p) = 1

log_mult: LEMMA log(p, x * y) = log(p, x) + log(p, y)

log_div: LEMMA log(p, x / y) = log(p, x) - log(p, y)

log_expt: LEMMA log(p, x^i) = i * log(p, x)

log_eq_0: LEMMA log(p, x) = 0 IMPLIES x = 1

log_ge_0: LEMMA x >= 1 IMPLIES log(p, x) >= 0

log_gt_0: LEMMA x > 1 IMPLIES log(p, x) > 0

log_strict_incr: LEMMA FORALL x, y: x < y IMPLIES log(p, x) < log(p, y)

log_incr: LEMMA FORALL x, y: x <= y IMPLIES log(p, x) <= log(p, y)

z, u: VAR real

exp(p, z): real = exp(z * ln(p))

log_exp: LEMMA log(p,exp(p, z)) = z

exp_log: LEMMA exp(p,log(p, x)) = x

exp_strict_incr: LEMMA FORALL z, u: z < u IMPLIES exp(p, z) < exp(p, u)

exp_incr: LEMMA FORALL z, u: z <= u IMPLIES exp(p, z) <= exp(p, u)

END log

log2: THEORY

BEGIN

IMPORTING log

x, y: VAR posreal

i, n: VAR posnat

log2(x): real = log(2, x)

log2_1: LEMMA log2(1) = 0

log2_2: LEMMA log2(2) = 1

log2_mult: LEMMA log2(x * y) = log2(x) + log2(y)

log2_div: LEMMA log2(x / y) = log2(x) - log2(y)

log2_expt: LEMMA log2(x^i) = i * log2(x)

log2_exp2: LEMMA log2(2^n) = n

log2_div_eq: LEMMA log2(x) = log2(x / 2^n) + n

log2_eq_0: LEMMA log2(x) = 0 IMPLIES x = 1

log2_ge_0: LEMMA x >= 1 IMPLIES log2(x) >= 0

log2_gt_0: LEMMA x > 1 IMPLIES log2(x) > 0

log2_lt1: LEMMA x < 2 IMPLIES log2(x) < 1

log2_strict_incr: LEMMA FORALL x, y: x < y IMPLIES log2(x) < log2(y)

log2_incr: LEMMA FORALL x, y: x <= y IMPLIES log2(x) <= log2(y)

z, u: VAR real

k: VAR nat

ex2(z): real = exp(2, z)

log2_ex2: LEMMA log2(ex2(z)) = z

ex2_log2: LEMMA ex2(log2(x)) = x

ex2_exp2: LEMMA ex2(k) = exp2(k)

ex2_strict_incr: LEMMA FORALL z, u: z < u IMPLIES ex2(z) < ex2(u)

ex2_incr: LEMMA FORALL z, u: z <= u IMPLIES ex2(z) <= ex2(u)

END log2

ilog2: THEORY

BEGIN

IMPORTING log2, ints@div_nat, fl

a, m: VAR posnat

x: VAR posreal

y: VAR { z: real | z >= 1}

ilog2(y): nat = floor(log2(y))

ilog2_1: LEMMA ilog2(1) = 0

ilog2_floor: LEMMA ilog2(y) = ilog2(floor(y))

ilog2_ge0: LEMMA ilog2(a) >= 0

ilog2_next: LEMMA x < 1 IMPLIES ilog2(a) = ilog2(a + x)

ilog2_div: LEMMA a >= 2^m IMPLIES ilog2(div(a, 2^m)) = ilog2(a / 2^m)

ilog2_divm: LEMMA a >= 2^m IMPLIES ilog2(div(a, 2^m)) = ilog2(a) - m

ilog2_div2: LEMMA a >= 2 IMPLIES ilog2(div(a, 2)) = ilog2(a) - 1

END ilog2

Ключевой и наиболее сложной для доказательства здесь была лемма ilog2_div.

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