суббота, 19 февраля 2011 г.

CORDIC

Имею задачу как можно быстрее вычислить десятичный логарифм с заданной точностью на микроконтроллерной платформе. Есть куча методов: табличный, таблично-интерполяционный, метод Бриггса (этот дядя в 1620г посчитал логарифм до 16 знака!!!). Обнаружил такой замечательный алгоритм, как CORDIC, для вычисления иррациональных функций типа логарифмов, арктангенсов и прочих синусов да косинусов, делений и умножений. http://www.ee.siue.edu/~gengel/pdf/cordic.pdf. Вкратце суть вычисления логарифма при помощи CORDIC:

log10(x*B1*...*Bn) = log10(x) + log10(B1) + ... + log10(Bn)

x*B1*...*Bn -> 1 => log10(x) -> -(log10(B1) + ... + log10(Bn))


Таким образом итерационно подбирая, коэффициенты
Bi, вычисляем искомое значение. Есть набор коэффициентов, при котором метод сходится, а именно:

Bn = (1 - 2^(-i)) если x*B1*...*Bn-1 больше 1

Bn = (1 + 2^(-i)) если x*B1*...*Bn-1 меньше 1


Точность метода оперделяется количеством итерраций. 10 итераций дают ошибку максимум в 3-м знаке. Значений
log10(Bi) хранятся в таблице. Такой выбор Bi позволяет расчитать логарифм для значений аргумента 0.4194 до 3.4627, что не может не печалить. Так как надо считать для целых значений аргумента от 1 до 4095. Есть несколько вариантов расширения границ применимости:
  • Найти правильный набор коэффициентов
  • Использовать хитрость y = (x-1)/(x+1), log10(x) = log10(1+y) - log10(1-y)
По первому пункту я ничего не придумал толкового, второй пункт не устраивает производительностью: два раза CORDIC + еще медленное деление. В голову спришла ледующая мысль:

log10(x) = log10(A) + log10(x/A)


Коэффициент A можно подобрать таким, чтобы x/A лежал в диапазоне применимости CORDIC, а чтобы еще быстро делить, можно взять в качестве А степени двойки. Ну еще добавляется таблица для значений log10(A).

5 комментариев:

  1. где то я уже видел похожую формулу

    ОтветитьУдалить
  2. А вы учли необходимость использования двойных итераций для вычисления логарифма?
    В.Д.Байков
    http://baykov.de/CORDIC1972.htm
    http://baykov.de/Cordic1975.htm
    Эту книжку можно полностью бесплатно в сети скачать.

    ОтветитьУдалить
  3. эммм, не совсем понял что за двойные итерации, если поясните буду очень признателен, просто сейчас пока отпала надобность считать логарифм, поэтому как-то CORDIC далее этого и следующего поста не продвинулся

    ОтветитьУдалить
  4. Начал я заниматься методами CORDIC, он же алгоритм Волдера или метод "цифра за цифрой"
    в 1963 году:
    http://proza.ru/2010/07/13/764

    Составил и отладил программы для машины "Проминь М". И тут стал замечать, что в ряде точек для функций lnX, expX, arcsinX, arccosX, shX, chX
    пробиваются очень большие ошибки.
    Начал разбираться и понял, что для них не выполняется правило "самоисправляемости итераций". Т.е., когда мы в середине итераций попадаем почти точно в результат, но вынуждены не останавливтаься на этом, а делать очередной шаг, хоть в одну, хоть в другую сторону, то этот шаг должен быть суммой последующих итераций "исправлен", т.е. компенсирован.
    Но если мы меняем модуль итерационного шага - (2**(-i)) -
    ВСЯКИЙ раз, то как показывает теория, этого сделать не удается. Поэтому в году так 1967-1968, когда и писал на эту тему диплом в ЛЭТИ, я и предложил изменять для перечисленных
    выше функций значение показателя степени (i)
    не каждую итерацию, а Ч Е Р Е З итерацию,
    т.е. дважды выполнять шаги с одним и тем же значением i, т.е. делать двойные итерации.
    Именно это и гарантирует сходимость метода во во все диапазоне.
    В 1972 году я защитил на тему CORDIC кандидатскую. Американцы предложили этот же способ двойных итераций только в 1992 году, было опуьликовано в IEEE Trans. on Computers,
    в 2006 году француз Muller описал этот способ в своей книге.

    ОтветитьУдалить
  5. Спасибо, очень интересные пояснения. Как раз назрела работа с вычислением иррациональных функций, так что вскоре опять буду смотреть в сторону CORDIC и возможно Ваши пояснения очень пригодятся

    ОтветитьУдалить