Вычисление логарифма подсчетом числа сдвигов

Сравнение с версией программы floor в файле fast_floor.cpp

В программе, находящейся в файле fast_floor.cpp (приложение 3) обнаружено две ошибки:

1. При E – bias >= m – 2 (или d >= w) сдвигаемое влево значение выходит за пределы m – 1 битов. В частности, при d = w старшая единица попадает в позицию знака в представлении целого числа, в результате чего результат floor становится отрицательным.

В реализации следовало бы запретить использовать функцию floor при d > 0, т.е. при E – bias – p + 1 > 0 из-за экспоненциально растущей погрешности.

2. Ненормализованное сверхмалое число обрабатывается так же, как нормализованное (с добавлением 1 в разряд p). Тем не менее это не меняет значения функции floor из-за малости числа. Однако при вычислении других функций результат был бы ошибочным.

Программа в файле fast_floor.cpp с точностью до обозначений и за исключением незначительных деталей подобна конечной программе (5) без первого оператора if (E = 0) { if (S=0) y = 0 else y = -1 }, который можно было бы опустить. Операторы int d = E – bias – p + 1; nat z = 2^(p-1) or T, присутствующие в (5), разнесены в fast_floor.cpp по разным ветвям, за счет чего программа в fast_floor.cpp, возможно, работает чуть быстрее.

4.Алгоритмы вычисления целочисленного
двоичного логарифма

Рассмотрим алгоритмы вычисления целой части двоичного логарифма, т.е. функции y = ilog2(x) = floor(log₂(x)), где floor ¾ функция вычисления целой части вещественного аргумента. Функцию ilog2 можно представить следующей спецификацией:

log2(real x: real y) pre x >0 post y = ln(x) / ln(2);

ilog2(nat a: nat r) pre a >= 1 post r = floor(log2(a));

Известно, что целочисленный двоичный логарифм ¾ это номер старшей единицы в двоичном разложении числа a. Простейший способ вычисления этого номера заключается в определении числа сдвигов (на один бит) вправо числа a до полного его обнуления. Результат сдвига вправо числа a на один бит есть div(a, 2). Используя приведенные соображения нетрудно построить следующий алгоритм:

ilog2(nat a: nat r) pre a >= 1

{ if (a = 1) r = 0 else r = ilog2(div(a, 2)) +1

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

Однако далеко не просто доказать корректность данного алгоритма.

Сначала построим программу с хвостовой рекурсией, используя следующее обобщение задачи:

ilg(nat a, m: nat r) pre a >= 1 post r = m + floor(log2(a));

Здесь параметр m используется в качестве накопителя числа сдвигов.

Программа представлена ниже в виде двух определений предикатов:

ilog2(nat a: nat r) pre a >= 1

{ ilg(a, 0: r)

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

ilg(nat a, m: nat r) pre a >= 1

{ if (a = 1) r = m else ilg(div(a, 2), m +1: r)

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

measure e(a) = a - 1

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

ilg: r <- m;

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

ilog2(nat a: nat r)

{ ilg(a, 0: r) }

ilg(nat a, r: nat r)

{ if (a = 1) r = r else ilg(div(a, 2), r +1: r) }

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

ilog2(nat a: nat r)

{ for (r = 0; a!= 1; r = r +1) a = div(a, 2); }

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

ilog2: THEORYBEGINIMPORTING ilg, loga, r: VAR natP_ilog2(a): bool = a >= 1Q_ilog2(a, r): bool = r = floor(log(2, a))satisfy_spec: LEMMA P_ilog2(a) IMPLIES EXISTS r: Q_ilog2(a, r)

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

P_ilog2(a) & Q_ilog2(a, r)├ LS(ilg(a, 0: r)).

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

FB1: LEMMA

P_ilog2(a) & Q_ilog2(a, r) IMPLIES P_ilg(a) & Q_ilg(a, 0, r).

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

ilg: THEORYBEGINa, m, r: VAR natP_ilg(a): bool = a >= 1Q_ilg(a, m, r): bool = r = m + floor(log(2, a))satisfy_spec: LEMMA P_ilg(a) IMPLIES EXISTS r: Q_ilg(a, m, r)e(a): nat = a - 1

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

FR1: Induct(a) & P_ilg(a) & Q_ilg(a, m, r)├

LS(if (a = 1) r = m else ilg(div(a, 2), m +1: r)).

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

FC1: LEMMA

P_ilg(a) & Q_ilg(a, m, r) & a = 1 IMPLIES r = m.

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

FC2: Induct(a) & P_ilg(a) & Q_ilg(a, m, r)& NOT a = 1 ├

LS(ilg(div(a, 2), m +1: r)).

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

FB3: LEMMA P_ilg(a) & Q_ilg(a, m, r) & NOT a = 1 IMPLIES

e(div(a, 2)) < e(a) & P_ilg(div(a, 2)) & Q_ilg(div(a, 2), m +1, r).

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