CORDIC

CORDIC (Метод CORDIC от англ. COordinate Rotation DIgital Computer — цифровой вычислитель поворота системы координат; метод «цифра за цифрой», алгоритм Волдера) — итерационный метод сведения прямых вычислений сложных функций к выполнению простых операций сложения и сдвига.

Идея метода заключается в сведении вычисления значений сложных (например, гиперболических) функций к набору простых шагов — сложению и сдвигу.

Такой подход особенно полезен при вычислении функций на устройствах с ограниченными вычислительными возможностями, такими как микроконтроллеры или программируемые логические матрицы (FPGA). Кроме того, поскольку шаги однотипны, то при аппаратной реализации алгоритм поддаётся развёртыванию в конвейер либо свертыванию в цикл.

Метод впервые был описан в применении к вычислению тригонометрических функций и операций преобразования координат Джеком Волдером в 1956 году, затем в 1959 году. В 1956 году Акушский и Юдицкий выдвинули аналогичную идею для вычисления показательных и логарифмических функций. Первоначально же близкая к этому идея была предложена Генри Бриггсом в 1624 году при составлении им таблиц логарифмов.

Метод был использован при реализации некоторых видов инструкций с плавающей точкой в математических сопроцессорах Intel, включая 8087^{[1][2][3][4][5]}, 80287^[5][6], 80387^[5][6] и вплоть до 80486^[1], а также в сопроцессорах Motorola 68881^[1][2] и 68882. Это позволило уменьшить число логических элементов (и сложность) сопроцессора.

Общая идея метода сводится к следующему. Последовательным умножением аргумента на заранее выбранные константы, приблизить аргумент с заданной точностью для одних функций к единице, для других функций — к нулю. Однако, для того, чтобы само значение функции при этом оставалось неизменным, необходимо одновременно совершать эквивалентные действия над выбранными константами. При выборе особым образом значений констант удаётся существенно упростить вычисления значений функции.

Например, Бриггс умножал значение аргумента функции десятичного логарифма ${\displaystyle \log X}$ на константы вида: ${\displaystyle 1+10^{-i}}$ либо ${\displaystyle 1-10^{-i}}$.

При этом выбор первого множителя (${\displaystyle 1+10^{-i}}$) осуществлялся в том случае, если текущее значение величины ${\displaystyle X}$ оказывалось меньше единицы, а второго — в том случае, если текущее ${\displaystyle X}$ было больше единицы. Выбирая последовательно значения показателя степени от 1 до ${\displaystyle N}$, где ${\displaystyle 10^{-N}}$ являлось максимальной допустимой ошибкой вычислений, Бриггс сводил значение ${\displaystyle X}$ с требуемой точностью к единице, а тем самым значение функции ${\displaystyle \log X}$ к нулю.

Однако, для равносильности указанных преобразований необходимо одновременно с умножением текущего ${\displaystyle X}$ на ${\displaystyle 1+10^{-i}}$ делить указанное значение на ту же самую величину. Но, как известно, логарифм частного равен разности логарифмов числителя и знаменателя. Следовательно, одновременно с умножением текущего ${\displaystyle X}$ на ${\displaystyle 1+10^{-i}}$ необходимо заранее вычисленную функцию логарифма значения ${\displaystyle 1+10^{-i}}$ отнимать от произведения аргумента ${\displaystyle X}$ на множитель ${\displaystyle 1+10^{-i}}$ .

Значения всех констант вида ${\displaystyle 1+10^{-i}}$ и ${\displaystyle 1-10^{-i}}$ могут быть вычислены заранее, поскольку их относительно немного. Например, при допустимой ошибке ${\displaystyle 10^{-6}}$ их всего двенадцать.

Остаётся пояснить, что умножение на константы вида ${\displaystyle 1+10^{-i}}$ и ${\displaystyle 1-10^{-i}}$ сводится к операциям сложения, вычитания и переноса запятой (сдвига).

Следовательно, процедура вычисления функции десятичного логарифма по Бриггсу сводится к операциям сложения, вычитания и десятичного сдвига.

Обобщение идеи метода Бриггса на комплексные числа было осуществлено в середине-конце пятидесятых годов Джеком Волдером и почти одновременно с ним Акушским и Юдицким. Это позволило вычислять тригонометрические функции.

CORDIC может быть использован для расчета ряда различных функций. Это объяснение показывает, как использовать CORDIC в режиме вращения для расчета синуса и косинуса угла. Предполагается, что желаемый угол задается в радианах и результаты представлены в формате с фиксированной запятой. Чтобы определить синус или косинус угла ${\displaystyle \beta }$, должны быть найдены координаты точки у или х на единичной окружности в соответствии с желаемым углом. Используя CORDIC, мы начинаем с вектора ${\displaystyle v_{0}}$:

${\displaystyle v_{0}={\begin{pmatrix}1\\0\end{pmatrix}}}$

В первой итерации этот вектор будет вращаться на 45° против часовой стрелки, чтобы получить вектор ${\displaystyle v_{1}}$. Последовательные итерации будут вращать вектор в одном или другом направлении с уменьшающимся шагом, пока желаемый угол не будет достигнут. Величина i-го шага — arctg(1/(2ⁱ⁻¹)), для i = 1, 2, 3, ….

Более формально, на каждой итерации вычисляется вращение, которое осуществляется путём умножения вектора ${\displaystyle v_{i-1}}$ на матрицу поворота ${\displaystyle R_{i}}$:

${\displaystyle v_{i}=R_{i}v_{i-1}\ }$

Матрицы вращения R определяются по формуле:

${\displaystyle R_{i}=\left({\begin{array}{rr}\cos \gamma _{i}&-\sin \gamma _{i}\\\sin \gamma _{i}&\cos \gamma _{i}\end{array}}\right)}$

Используя следующие два тригонометрические тождества

${\displaystyle {\begin{aligned}\cos \alpha &=&{1 \over {\sqrt {1+\operatorname {tg} ^{2}\alpha }}}\\\sin \alpha &=&{{\operatorname {tg} \alpha } \over {\sqrt {1+\operatorname {tg} ^{2}\alpha }}}\end{aligned}}}$

получаем матрицу поворота:

${\displaystyle R_{i}={1 \over {\sqrt {1+\operatorname {tg} ^{2}\gamma _{i}}}}{\begin{pmatrix}1&-\operatorname {tg} \gamma _{i}\\\operatorname {tg} \gamma _{i}&1\end{pmatrix}}}$

Выражение для поворачиваемого вектора ${\displaystyle v_{i}=R_{i}v_{i-1}}$:

${\displaystyle v_{i}={1 \over {\sqrt {1+\operatorname {tg} ^{2}\gamma _{i}}}}{\begin{pmatrix}x_{i-1}&-&y_{i-1}\operatorname {tg} \gamma _{i}\\x_{i-1}\operatorname {tg} \gamma _{i}&+&y_{i-1}\end{pmatrix}}}$

где ${\displaystyle x_{i-1}}$ и ${\displaystyle y_{i-1}}$ компоненты ${\displaystyle v_{i-1}}$. Ограничивая значения углов ${\displaystyle \gamma _{i}}$ так что ${\displaystyle \operatorname {tg} \gamma _{i}}$ принимает значения ${\displaystyle \pm 2^{-i}}$ умножение на тангенс можно заменить делением на степень двойки, которое эффективно реализовано в цифровом аппаратном обеспечении компьютера с помощью битового сдвига. Получаем выражение:

${\displaystyle v_{i}=K_{i}{\begin{pmatrix}x_{i-1}&-&\sigma _{i}2^{-i}y_{i-1}\\\sigma _{i}2^{-i}x_{i-1}&+&y_{i-1}\end{pmatrix}}}$

где

${\displaystyle K_{i}={1 \over {\sqrt {1+2^{-2i}}}}}$

и ${\displaystyle \sigma _{i}}$ может иметь значения −1 или 1 и используется для определения направления вращения: если угол ${\displaystyle \gamma _{i}}$ является положительным, то ${\displaystyle \sigma _{i}}$ равняется 1, в противном случае он равен −1.

Мы можем игнорировать ${\displaystyle K_{i}}$ в итерационном процессе, а затем применить его потом для получения коэффициента масштабирования:

${\displaystyle K(n)=\prod _{i=0}^{n-1}K_{i}=\prod _{i=0}^{n-1}1/{\sqrt {1+2^{-2i}}}}$

который рассчитывается заранее и хранится в таблице, или как одна константа, если число итераций фиксировано. Эта поправка также может быть сделана заранее, путём масштабирования ${\displaystyle v_{0}}$.

Единственная задача которая осталась, это определить по часовой стрелке или против часовой стрелки должно выполняться вращение на каждой итерации (выбор значения ${\displaystyle \sigma }$). Это делается путём отслеживания величины поворота на каждой итерации и вычитания из желаемого угла, а затем проверки, если ${\displaystyle \beta _{n+1}}$ является положительным, то мы должны вращаться по часовой стрелке или если он отрицательный, мы должны вращаться против часовой стрелки, чтобы приблизится к углу ${\displaystyle \beta }$.

${\displaystyle \beta _{i}=\beta _{i-1}-\sigma _{i}\gamma _{i}.\quad \gamma _{i}=\operatorname {arctg} 2^{-i},}$

Значения ${\displaystyle \gamma _{n}}$ также должны быть посчитаны заранее. Но для малых углов ${\displaystyle \operatorname {arctg} (\gamma _{n})=\gamma _{n}}$ в представлении с фиксированной точкой, что позволяет снизить размер таблицы.

Как видно на иллюстрации выше, синус угла ${\displaystyle \beta }$ является y координатой окончательного вектора ${\displaystyle v_{n}}$, а x координата — значением косинуса.

В работах ^[7] и ^[8] было предложено использовать при вычислении функций arcsinX, arccosX, lnX, expX, а также для вычисления гиперболических функций метод "двойных итераций" (double iterations). Состоит он в том, что в отличие от классического метода CORDIC, где величина шага итерации меняется ВСЯКИЙ раз, т.е. на каждой итерации, при методе двойных итераций величина шага итерации повторяется дважды и меняется лишь через одну итерацию. Отсюда появилось обозначение для показателя степени при двойных итерациях: i= 1,1,2,2,3,3... Тогда как при обычных итерациях: i=1,2,3... Метод двойных итераций гарантирует сходимость метода во всем допустимом диапазоне изменения аргументов.

Обобщение вопросов сходимости алгоритмов метода CORDIC на произвольную позиционную систему счисления с основанием R, приведённое в работе^[9], показало, что для функций sin, cos, arctg достаточно выполнения (R-1)-итераций при каждом значении i (i от 0 или 1 до n, где n - разрядность чисел), т.е. для каждого разряда результата. Для функций ln, exp, sh, ch, arth следует выполнять R итераций для каждого значения i. Для функций arcsin и arccos следует выполнять 2 ⋅ (R-1) итераций на каждый разряд, т.е. при каждом значении i.

Для функций arsh, arch количество итераций составит 2R для каждого i, т.е. для каждого разряда результата.

• Jack E. Volder, The CORDIC Trigonometric Computing Technique, IRE Transactions on Electronic Computers, September 1959
• Daggett, D. H., Decimal-Binary conversions in CORDIC, IRE Transactions on Electronic Computers, Vol. EC-8 #5, pp 335–339, IRE, September 1959.
• J. E. Meggitt, Pseudo Division and Pseudo Multiplication Processes, IBM Journal, April 1962.
• Байков В. Д. Вопросы исследования вычисления элементарных функций по методу «цифра за цифрой», автореферат диссертации на соискание учёной степени кандидата технических наук, Ленинград, ЛЭТИ, 1972
• Schmid, Hermann, Decimal computation. New York, Wiley, 1974
• Байков В. Д., Смолов В. Б. Аппаратурная реализация элементарных функций в ЦВМ, Ленинград, изд-во ЛГУ, 1975, 96 стр.
• Байков В. Д., Селютин С. А., Вычисление элементарных функций в ЭКВМ, Москва, Радио и связь, 1982, 64 стр.
• Байков В. Д., Смолов В. Б. Специализированные процессоры: итерационные алгоритмы и структуры, Москва, «Радио и связь», 1985, 288 стр.
• Andraka, Ray, A survey of CORDIC algorithms for FPGA based computers
• Henry Briggs, Arithmetica Logarithmica. London, 1624, folio
• M. Zechmeister Solving Kepler’s equation with CORDIC double iterations Institut für Astrophysik, Georg-August-Universität, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany, Preprint 10 August 2020, p.1-10.
• Tomás Lang & Elisardo Antelo CORDIC-Based Computation of ArcCos Journal of VLSI signal processing systems for signal, image and video technology volume 25, pages19–38 (2000)
• Vladimir Baykov Double iterations in CORDIC CORDIC Bibliography
• Digit by digit methods, Jacques Laporte, Paris 2006 (la filiation entre CORDIC et les anciennes méthodes notamment celle de Briggs)
• CORDIC Bibliography Site, Shaoyun Wang, July 2011
• CORDIC Bibliography Site