Один способ вычисления логарифма по основанию 2

Вычисление логарифмов довольно распространённая операция в цифровой обработке сигналов. Чаще пожалуй приходится считать только свёртки (умножение с накоплением) и амплитуды с фазами. Как правило для вычисления логарифмов на FPGA применяется алгоритм CORDIC в гиперболическом варианте, требующий только таблицы и простых арифметических операций. Однако это не всегда бывает удобно, особенно если проект большой, кристалл маленький и начинаются танцы с оптимизацией. Именно с такой ситуацией и пришлось мне однажды столкнуться. Оба порта блока RAM (Cyclone IV) уже плотненько были в работе, не оставляя свободных окон. Использовать ещё один блок под гиперболический CORDIC не хотелось. Зато был умножитель, для которого во временной диаграмме получалось приличное свободное окно. Денёк подумав, я сочинил следующий алгоритм, в котором не используется таблиц, но есть умножение, точнее возведение в квадрат. И поскольку схемотехнически возведение в квадрат проще общего случая умножения, возможно этот алгоритм представляет интерес для специализированных микросхем, хотя для FPGA разницы конечно нет. Подробнее под катом.

Объяснять что к чему, тут проще для действительных чисел. С них и начнём. К целочисленной реализации перейдём потом.
Пусть есть число X. Требуется найти число Y такое, чтобы $$.
Положим так же что X лежит в интервале от 1 до 2. Это не слишком ограничивает общность, поскольку X всегда можно перевести в этот интервал умножением или делением на степень двойки. Для Y это будет означать добавление или вычитание целого числа, что делается легко. Итак X лежит в интервале от 1 до 2. Тогда Y будет лежать в интервале от 0 до 1. Запишем Y как бесконечную двоичную дробь:

$$

Коэффициенты $$ в этой записи есть ни что иное, как просто биты двоичного представления числа Y. Причём поскольку Y меньше 1, очевидно что $$=0.

Возведём наше первое уравнение в квадрат: $$ и как и ранее, запишем двоичное представление числа 2Y. Очевидно, что

$$

Т.е. биты $$ остались теми же, просто сдвинулись степени двойки. Бита $$ в представлении нет, потому что он равен нулю.
Возможны два случая:
1) $$, Y > 1, $$
2) $$, Y < 1, $$
В первом случае в качестве нового значения X примем $$, во втором — $$.
В итоге задача свелась к прежней. Новое X опять лежит в интервале от 1 до 2, новое Y от 0 до 1. Но мы узнали один бит результата. Делая такие же шаги в дальнейшем, мы можем получить сколько угодно битов Y.
Посмотрим как это работает в программе на С:

#include <stdio.h>
#include <math.h>
int main()
{
	double w=1.4;
	double s=0.0;
	double a=0.5;
	double u=w;
	for(int i=0; i<16; i++)
	{
		u=u*u;
		if(u>2)
		{
			u=u/2;
			s+=a;
		}
		a*=0.5;
	}
	w=log2(w);
	double err=100*abs(2*(s-w)/(s+w));
	printf("res=%f, log=%f, err=%f%c\n",s,w,err,'%');
	return 0;
}

Мы посчитали логарифм с точностью до 16 бит и сравнили с тем, что даёт математическая библиотека. Программа вывела:
res=0.485413, log=0.485427, err=0.002931%
Результат совпал с библиотекой с точностью 0.003%, что показывает работоспособность нашего алгоритма.
Перейдём к целочисленной реализации.
Пусть N-разрядные двоичные беззнаковые числа, отображают интервал [0, 1]. Для удобства будем считать единицей число $$, а не $$, и соответственно двойкой число $$. Напишем программу по образу и подобию предыдущей, но работающую с целыми числами:

#include <stdio.h>
#include <math.h>
#define DIG 	18            //разрядность чисел
#define N_BITS  16            //число бит которое считаем
unsigned ONE=1<<(DIG-1);      //единица
unsigned TWO=ONE<<1;          //двойка
unsigned SCALE=1<<(N_BITS+1);   //масштаб логарифмов
unsigned myLog(unsigned w)
{
	unsigned s=0;
	unsigned long long u=w;
	for(int i=0; i<N_BITS+1; i++)
	{
		s<<=1;
		u=(u*u)>>(DIG-1);  //берём только старшие разряды !
		if(u&TWO)   //тут достаточно одного бита для сравнения
		{
			u>>=1;
			s+=1;
		}
		printf("%X\n", (int)u);
	}
	return s;
}
int main()
{
	double w=1.2345678;
	unsigned iw=(unsigned)(ONE*w);
	double dlog=log2(w);
	unsigned ilog=myLog(iw);
	unsigned test=(unsigned)(SCALE*dlog);
	int err=abs((int)(ilog-test));
	printf("val=0x%X, res=0x%X, log=0x%X, err=%d\n",iw,ilog,test,err);
	return 0;
}

Поиграв в программе с разными разрядностью чисел (DIG), точностью вычислений (N_BITS), и аргументами логарифма (w), видим, что всё вычисляется правильно. В частности с теми параметрами, которые заданны в этом исходнике, программа выдаёт:
val=0x27819, res=0x9BA5, log=0x9BA6, err=1

Теперь всё готово для того, чтобы реализовать на верилоге железку, делающую точно то же самое, что и функция myLog на С. Переменные s и u в нашей функции можно распечатать в цикле и сравнивать с тем, что выдает симулятор верилога. Соответствие этих переменных железной реализации весьма прозрачно и понятно. u это рабочий регистр, принимающий во время итераций новые значения X. s это сдвиговый регистр, в котором накапливается результат. Интерфейс нашего модуля будет выглядеть вот так:

module logarithm(
		input                 clk,      //клок
		input                 wr,       //сигнал записи данных
		input[17:0]           din,      //входная шина данных
		output[nbits-1:0]     dout,     //выходная шина данных
		output                rdy       //сигнал готовности		
		);
	parameter nbits=16;       //точность вычисления   

Входная шина принята 18-разрядной, соответственно разрядности умножителей в Cyclone IV. Числа на наш модуль должны поступать нормализованные. Т.е. со старшим битом, равным единице. В моём проекте это выполнялось автоматически. Но в случае чего реализовать нормализатор, думаю никому труда не составит. Точность вычислений задаётся параметром nbits, по умолчанию равным 16. Модуль считает по одному биту за такт, и за 16 тактов вычислит логарифм с точностью до 16 битов. Если надо быстрее с той же точностью или точнее с той же скоростью, надеюсь никому не составит особого труда разделить модуль на несколько устройств и конвейеризовать.

Запустим тест вот таким скриптом:

#!/bin/sh
rm -f *.vvp
rm -f *.vcd
iverilog -o testbench.vvp logarithm.v testbench.v
vvp testbench.vvp
gtkwave testbench.vcd testbench.gtkw 

Запустив тест, видим окончательный вывод симулятора — value=27819, result=9ba5. Верилог выдал то же самое что С. Временная диаграмма тут довольно тривиальна и особого интереса не представляет. Поэтому её не привожу.

Убеждаемся, что они совпадают бит в бит. Итого, реализация на верилоге повторяет с точностью до бита модель на С. Это и есть тот результат, которого следует достигать, реализуя алгоритмы в железе.

На этом пожалуй всё. Надеюсь кому-то этот мой опыт окажется полезен.