Модулярная арифметика - Вклад участника [ru]

Метод умножения Шёнхаге — Штрассена (New)

2014-02-17T07:27:57Z

Myachikov:

'''Алгоритм Шёнхаге - Штрассена''' - асимптотически быстрый алгоритм умножения больших чисел. Он был разработан Арнольдом Шёнхаге и Фолькером Штрассеном в 1971. Сложность алгоритма в нотации "О - большого" составляет O(''N'' log ''N'' log log ''N''). Алгоритм
рекурсивно использует быстрое преобразование Фурье над кольцом из 22''n'' + 1 элементов, специальный тип дискретного преобразования Фурье - теоретико-числовое преобразование.

Алгоритм Шёнхаге - Штрассена являлся наиболее быстрым асимптотически известным методом умножения с 1971 до 2007, когда Фюрером был предложен более быстрый метод. Тем не менее, алгоритм Фюрера в настоящее время имеет преимущество только для астрономически больших чисел и не используется на практике.

На практике алгоритм Шёнхаге - Штрассена начинает превосходить в скорости старые методы, такие как алгоритм Карацубы и алгоритм Тоома-Кука для чисел между 2215 и 2217 (от 10,000 до 40,000 десятичных разрядов).

Приложения алгоритма Шёнхаге - Штрассена включают в себя как чисто математические задачи, такие как поиск простых чисел Мерсенна и вычисление цифр числа ''π'', так и практические, такие как вычисление подстановки Кронекера, в которой умножение многочленов с целыми коэффициентами можно эффективно заменить умножением больших чисел; это используется на практике by GMP-ECM для факторизации целых чисел методом эллиптических кривых Ленстры.

== Теория ==

=== Свертка ===
Допустим, что мы перемножаем числа 123 и 456 «в столбик», однако без выполнения переноса. Результат будет выглядеть так:

{|width=300
| || || 1 || 2 || 3
|-
| || ×|| 4 || 5 || 6
|-
|colspan=5|<hr />
|-
| || || 6 || 12 || 18
|-
| || 5 || 10 || 15 ||
|-
| 4 || 8 || 12 || ||
|-
|colspan=5|<hr />
|-
| 4 || 13 || 28 || 27 || 18
|}

Эта последовательность (4, 13, 28, 27, 18) называется ''ациклической'' или ''линейной свёрткой'' от последовательностей (1,2,3) и (4,5,6). Зная ациклическую свёртку двух последовательностей, рассчитать произведение несложно: достаточно выполнить перенос (например, в самом правом столбце, мы оставляем 8 и добавляем 1 к столбцу, содержащему 27). В нашем примере это приводит к результату 56088.

Есть и другие типы свёрток, которые могут быть полезны. Допустим, что входящие последовательности содержат ''n'' элементов (в примере — 3). Тогда результирующая линейная свёртка содержит ''n'' + ''n'' − 1 элементов; если мы возьмём самый правый столбец ''n'' элементов и добавим самый левый столбец ''n'' − 1 ', в результате мы получим циклическую свёртку:

{|width=300
| || 28 || 27 || 18
|-
| + || || 4 || 13
|-
|colspan=5|<hr />
|-
| || 28 || 31 || 31
|}

Если мы произведём перенос при циклическом свёртывании, результат будет тот же, что и при произведении чисел по модулю B''n'' − 1 (в данном примере это 103 − 1 = 999). Выполним перенос по (28, 31, 31) и получим 3141, при этом 3141 ≡ 56088 (mod 999).

Наоборот, если мы возьмём самый правый столбец ''n'' элементов и ''вычтем'' самый левый столбец ''n''−1 элементов, то в результате мы получим обратную свёртку:

{|width=300
| || 28 || 27 || 18
|-
| − || || 4 || 13
|-
|colspan=5|<hr />
|-
| || 28 || 23 || 5
|}

Если мы произведём перенос при обратном свёртывании, то результат будет тот же, что и при произведении чисел по модулю B''n'' + 1. В данном примере, 103 + 1 = 1001, выполним перенос по (28, 23, 5) и получим 3035, при этом 3035 ≡ 56088 (mod 1001). Обратная свёртка может содержать отрицательные числа, которые могут быть убраны во время переноса, используя ту же технику, что и при длинных вычитаниях.

=== Теорема о свёртке ===

Как и другие методы, основанные на быстром преобразовании Фурье, алгоритм Шёнхаге - Штрассена в корне зависит от теоремы о свёртке, которая обеспечивает эффективный способ посчитать циклическую свёртку двух последовательностей. Её идея состоит в следующем:

: Циклическая свёртка двух векторов может быть найдена через [[дискретное преобразование Фурье]] (ДПФ) каждого из них, путём произведения результирующих векторов элемент за элементом, с последующим обратным преобразованием Фурье (ОДПФ).

Или через формулы:

: CyclicConvolution(X, Y) = IDFT(DFT(X) · DFT(Y)), где:
:: CyclicConvolution — ''циклическая свертка'',
:: DFT — ''дискретное преобразование Фурье'',
:: IDFT — ''обратное дискретное преобразование Фурье''.

Если мы посчитаем ДПФ и ОДПФ, используя быстрое преобразование Фурье, и вызовем наш алгоритм перемножения рекурсивно, чтобы перемножить входы(?) преобразованных векторов DFT(X) и DFT(Y), то в результате мы получим эффективный алгоритм для расчёта циклической свёртки.

В этом алгоритме, гораздо эффективней считать ''обратную циклическую'' свёртку; как оказывается, немного модифицированная версия теоремы о свёртке может позволить и это. Предположим, что вектора X и Y имеют длину ''n'', и ''a'' — примитивный корень порядка 2''n'' (это означает, что ''a''2''n'' = 1 и все меньшие степени ''a'' не равны 1). Таким образом мы можем определить третий вектор ''A'', называемый ''вектор веса'', обладающий следующими свойствами:
[[Файл:DIT-FFT-butterfly.png|мини|Операция «[[бабочка (БПФ)|бабочка]]».]]

: ''A'' = (''a''''j''), 0 ≤ ''j'' < ''n''
: ''A''−1 = (''a''−j''), 0 ≤ ''j'' < ''n''

Теперь мы можем записать:

: NegacyclicConvolution(''X'', ''Y'') = ''A''−1 · IDFT(DFT(''A'' · ''X'') · DFT(''A'' · ''Y'')), где
:: NegacyclicConvolution — ''Обратная циклическая свертка'',
:: DFT — ''дискретное преобразование Фурье'',
:: IDFT — ''обратное дискретное преобразование Фурье''.

Другими словами, это то же самое за исключением того, что входящие векторы умножены на ''A'', а результат умножен на ''A''−1.

=== Выбор кольца ===

Дискретное преобразование Фурье — абстрактная операция, которая может быть выполнена в любом алгебраическом [[Кольцо (математика)|кольце]]; обычно оно берётся из поля комплексных чисел, но фактически использовать комплексную арифметику с достаточной точностью, чтобы обеспечить точные результаты, медленно и неэффективно. Вместо этого мы можем использовать теоретико-числовое преобразование, которое производит преобразование в поле целых чисел по модулю N для некоторого целого N.

Так же как есть примитивные корни единицы любого порядка на комплексной плоскости, при любом заданном ''n'' мы можем выбрать подходящее N такое, что ''b'' — примитивный корень единицы порядка ''n'' в поле целых чисел по модулю N (другими словами, ''b''''n'' ≡ 1 (mod N), и все меньшие степени ''b'' не равны 1 mod N).

Алгоритм тратит большую часть времени на рекурсивное выполнение произведения меньших чисел; в простом варианте алгоритма это происходит в ряде мест:

# Внутри алгоритма быстрого преобразования Фурье, примитивный корень единицы ''b'' неоднократно возводится в степень и умножается на другие числа.
# При возведении в степень примитивного корня единицы ''a'' для получения вектора веса A с последующим умножением векторов A или A−1 на другие вектора.
# При выполнении последовательного перемножения преобразованных векторов.

Ключевой момент — выбрать N, модуль, равный 2''n'' + 1 для некоторого целого ''n''. У этого способа есть ряд преимуществ в ряде стандартных систем, в которых большие целые числа представлены в двоичном виде:

* Любое число может быть быстро уменьшено по модулю 2''n'' + 1 используя только сдвиг и сложение.
* Любые примитивные корни единицы в этом кольце могут быть записаны в форме 2''k''; соответственно мы можем умножать или делить любое число на корень из единицы используя сдвиг.
* Поэлементное рекурсивное перемножение преобразованный векторов может быть выполнено, используя обратную свёртку, которая работает быстрее, чем ациклическая свёртка, и в которой уже есть уменьшение результата
по модулю 2''n'' + 1.

=== Пример ===
Следующий пример иллюстрирует работу алгоритма.
Положим задано два числа: X и Y, каждое длиной 8 бит. Требуется найти значение RES = (X*Y) mod 2N+1
Возьмём для определённости X = Y = 12910 = 1000 00012.
Выберем N = 8 + 8 = 16.
* Этап 1. Разбиение чисел на K = 4 слов длины M = 4 бита:
**Дополним числа нулями до длины N = 16 и разобьём на слова. Получим две последовательности (в данном примере идентичные): {0001, 1000, 0000, 0000}. Или в десятичном виде: {3, 8, 0, 0} (слова расположены в порядке от младших бит к старшим). Все операции далее выполняются по модулю 2n+1, где n выбирается из условия n >= 2*M + log2(K) и является степенью двойки, то есть в нашем случае n >= 2*4 + 2 = 10, n = 16.

* Этап 2. Вычисление свёртки двух последовательностей:
** Этап 2.1 Умножение последовательностей на весовые коэффициенты. 2K-ый корень из единицы в кольце из 2n+1 элементов равен 22n/2K = 24 = 16. Результат умножения на вектор весовых коэффициентов в десятичном виде: {1, 128, 0, 0}.
**Этап 2.2 Прямое преобразование Фурье над кольцом из 2n элементов. K-ый корень из единицы в кольце вычетов по модулю 2n+1 равен 22n/2K = 24 = 256. Формула преобразования: <math>b_i = \sum^{K-1}_{j=0}{a_j\omega^{ij}}</math>, <math> \omega </math> - К-тый корень из единицы в кольце вычетов по модулю 2n+1. Результат: {129, 32769, 65410, 32770}.
**Этап 2.3 Покомпонентное произведение последовательностей в кольце из 2n элементов. Результат: {16641, 49153, 16129, 49155}.
**Этап 2.4 Обратное преобразование Фурье над кольцом из 2n элементов. Формула обратного преобразования: <math>b_i = \frac {1}{K} \sum^{K-1}_{j=0}{a_j\omega^{ij}}</math>, <math> \frac {1}{K} </math> - обратный к <math>K</math> по модулю 2n+1. Результат: {1, 256, 16384, 0}.
**Этап 2.5 Умножение последовательностей на коэффициенты, обратные к используемым в 2.1. Результат: {1, 16, 64, 0}.
* Этап 3. Выполнение операции переноса. Результат: {1, 0, 1, 4}. В двоичном виде: {0001, 0000, 0001, 0100}
* Этап 5. "Склейка" полученной последовательности, окончательный результат: 10100000010002 = 16641 = 129 * 129.
* Этап 5. Вычисление остатка по модулю 2N+1. В данном случае 16641 = 16641 mod 216+1

=== Описание ===
* [http://gmplib.org/manual/FFT-Multiplication.html#FFT-Multiplication Документация GMP LIB]
* [http://hal.archives-ouvertes.fr/docs/00/12/64/62/PDF/fft.pdf Статья A GMP-BASED IMPLEMENTATION OF SCHONHAGE-STRASSEN’S LARGE INTEGER MULTIPLICATION ALGORITHM]
* [http://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm Описание в английской Wikipedia]

=== Исходники ===
[http://www.ginac.de/CLN/cln.git/?p=cln.git;a=blob;f=src/base/digitseq/cl_DS_mul_fftm.h;h=a984802ce039e8179a26fda2681a2e133dba6fbd;hb=refs/heads/master CLN - Class Library for Numbers]

[https://github.com/myachik/ss_mult/blob/master/ss_mult.cpp Вариант реализации алгоритма]

Метод умножения Шёнхаге — Штрассена (New)

2014-02-06T08:57:56Z

Myachikov:

'''Алгоритм Шёнхаге - Штрассена''' - асимптотически быстрый алгоритм умножения больших чисел. Он был разработан Арнольдом Шёнхаге и Фолькером Штрассеном в 1971. Сложность алгоритма в нотации "О - большого" составляет O(''N'' log ''N'' log log ''N''). Алгоритм
рекурсивно использует быстрое преобразование Фурье над кольцом из 22''n'' + 1 элементов, специальный тип дискретного преобразования Фурье - теоретико-числовое преобразование.

Алгоритм Шёнхаге - Штрассена являлся наиболее быстрым асимптотически известным методом умножения с 1971 до 2007, когда Фюрером был предложен более быстрый метод. Тем не менее, алгоритм Фюрера в настоящее время имеет преимущество только для астрономически больших чисел и не используется на практике.

На практике алгоритм Шёнхаге - Штрассена начинает превосходить в скорости старые методы, такие как алгоритм Карацубы и алгоритм Тоома-Кука для чисел между 2215 и 2217 (от 10,000 до 40,000 десятичных разрядов).

Приложения алгоритма Шёнхаге - Штрассена включают в себя как чисто математические задачи, такие как поиск простых чисел Мерсенна и вычисление цифр числа ''π'', так и практические, такие как вычисление подстановки Кронекера, в которой умножение многочленов с целыми коэффициентами можно эффективно заменить умножением больших чисел; это используется на практике by GMP-ECM для факторизации целых чисел методом эллиптических кривых Ленстры.

== Теория ==

=== Свертка ===
Допустим, что мы перемножаем числа 123 и 456 «в столбик», однако без выполнения переноса. Результат будет выглядеть так:

{|width=300
| || || 1 || 2 || 3
|-
| || ×|| 4 || 5 || 6
|-
|colspan=5|<hr />
|-
| || || 6 || 12 || 18
|-
| || 5 || 10 || 15 ||
|-
| 4 || 8 || 12 || ||
|-
|colspan=5|<hr />
|-
| 4 || 13 || 28 || 27 || 18
|}

Эта последовательность (4, 13, 28, 27, 18) называется ''ациклической'' или ''линейной свёрткой'' от последовательностей (1,2,3) и (4,5,6). Зная ациклическую свёртку двух последовательностей, рассчитать произведение несложно: достаточно выполнить перенос (например, в самом правом столбце, мы оставляем 8 и добавляем 1 к столбцу, содержащему 27). В нашем примере это приводит к результату 56088.

Есть и другие типы свёрток, которые могут быть полезны. Допустим, что входящие последовательности содержат ''n'' элементов (в примере — 3). Тогда результирующая линейная свёртка содержит ''n'' + ''n'' − 1 элементов; если мы возьмём самый правый столбец ''n'' элементов и добавим самый левый столбец ''n'' − 1 ', в результате мы получим циклическую свёртку:

{|width=300
| || 28 || 27 || 18
|-
| + || || 4 || 13
|-
|colspan=5|<hr />
|-
| || 28 || 31 || 31
|}

Если мы произведём перенос при циклическом свёртывании, результат будет тот же, что и при произведении чисел по модулю B''n'' − 1 (в данном примере это 103 − 1 = 999). Выполним перенос по (28, 31, 31) и получим 3141, при этом 3141 ≡ 56088 (mod 999).

Наоборот, если мы возьмём самый правый столбец ''n'' элементов и ''вычтем'' самый левый столбец ''n''−1 элементов, то в результате мы получим обратную свёртку:

{|width=300
| || 28 || 27 || 18
|-
| − || || 4 || 13
|-
|colspan=5|<hr />
|-
| || 28 || 23 || 5
|}

Если мы произведём перенос при обратном свёртывании, то результат будет тот же, что и при произведении чисел по модулю B''n'' + 1. В данном примере, 103 + 1 = 1001, выполним перенос по (28, 23, 5) и получим 3035, при этом 3035 ≡ 56088 (mod 1001). Обратная свёртка может содержать отрицательные числа, которые могут быть убраны во время переноса, используя ту же технику, что и при длинных вычитаниях.

=== Теорема о свёртке ===

Как и другие методы, основанные на быстром преобразовании Фурье, алгоритм Шёнхаге - Штрассена в корне зависит от теоремы о свёртке, которая обеспечивает эффективный способ посчитать циклическую свёртку двух последовательностей. Её идея состоит в следующем:

: Циклическая свёртка двух векторов может быть найдена через [[дискретное преобразование Фурье]] (ДПФ) каждого из них, путём произведения результирующих векторов элемент за элементом, с последующим обратным преобразованием Фурье (ОДПФ).

Или через формулы:

: CyclicConvolution(X, Y) = IDFT(DFT(X) · DFT(Y)), где:
:: CyclicConvolution — ''циклическая свертка'',
:: DFT — ''дискретное преобразование Фурье'',
:: IDFT — ''обратное дискретное преобразование Фурье''.

Если мы посчитаем ДПФ и ОДПФ, используя быстрое преобразование Фурье, и вызовем наш алгоритм перемножения рекурсивно, чтобы перемножить входы(?) преобразованных векторов DFT(X) и DFT(Y), то в результате мы получим эффективный алгоритм для расчёта циклической свёртки.

В этом алгоритме, гораздо эффективней считать ''обратную циклическую'' свёртку; как оказывается, немного модифицированная версия теоремы о свёртке может позволить и это. Предположим, что вектора X и Y имеют длину ''n'', и ''a'' — примитивный корень порядка 2''n'' (это означает, что ''a''2''n'' = 1 и все меньшие степени ''a'' не равны 1). Таким образом мы можем определить третий вектор ''A'', называемый ''вектор веса'', обладающий следующими свойствами:
[[Файл:DIT-FFT-butterfly.png|мини|Операция «[[бабочка (БПФ)|бабочка]]».]]

: ''A'' = (''a''''j''), 0 ≤ ''j'' < ''n''
: ''A''−1 = (''a''−j''), 0 ≤ ''j'' < ''n''

Теперь мы можем записать:

: NegacyclicConvolution(''X'', ''Y'') = ''A''−1 · IDFT(DFT(''A'' · ''X'') · DFT(''A'' · ''Y'')), где
:: NegacyclicConvolution — ''Обратная циклическая свертка'',
:: DFT — ''дискретное преобразование Фурье'',
:: IDFT — ''обратное дискретное преобразование Фурье''.

Другими словами, это то же самое за исключением того, что входящие векторы умножены на ''A'', а результат умножен на ''A''−1.

=== Выбор кольца ===

Дискретное преобразование Фурье — абстрактная операция, которая может быть выполнена в любом алгебраическом [[Кольцо (математика)|кольце]]; обычно оно берётся из поля комплексных чисел, но фактически использовать комплексную арифметику с достаточной точностью, чтобы обеспечить точные результаты, медленно и неэффективно. Вместо этого мы можем использовать теоретико-числовое преобразование, которое производит преобразование в поле целых чисел по модулю N для некоторого целого N.

Так же как есть примитивные корни единицы любого порядка на комплексной плоскости, при любом заданном ''n'' мы можем выбрать подходящее N такое, что ''b'' — примитивный корень единицы порядка ''n'' в поле целых чисел по модулю N (другими словами, ''b''''n'' ≡ 1 (mod N), и все меньшие степени ''b'' не равны 1 mod N).

Алгоритм тратит большую часть времени на рекурсивное выполнение произведения меньших чисел; в простом варианте алгоритма это происходит в ряде мест:

# Внутри алгоритма быстрого преобразования Фурье, примитивный корень единицы ''b'' неоднократно возводится в степень и умножается на другие числа.
# При возведении в степень примитивного корня единицы ''a'' для получения вектора веса A с последующим умножением векторов A или A−1 на другие вектора.
# При выполнении последовательного перемножения преобразованных векторов.

Ключевой момент — выбрать N, модуль, равный 2''n'' + 1 для некоторого целого ''n''. У этого способа есть ряд преимуществ в ряде стандартных систем, в которых большие целые числа представлены в двоичном виде:

* Любое число может быть быстро уменьшено по модулю 2''n'' + 1 используя только сдвиг и сложение.
* Любые примитивные корни единицы в этом кольце могут быть записаны в форме 2''k''; соответственно мы можем умножать или делить любое число на корень из единицы используя сдвиг.
* Поэлементное рекурсивное перемножение преобразованный векторов может быть выполнено, используя обратную свёртку, которая работает быстрее, чем ациклическая свёртка, и в которой уже есть уменьшение результата
по модулю 2''n'' + 1.

=== Пример ===
Следующий пример иллюстрирует работу алгоритма.
Положим задано два числа: X и Y, каждое длиной 8 бит. Требуется найти значение RES = (X*Y) mod 2N+1
Возьмём для определённости X = Y = 12910 = 1000 00012.
Выберем N = 8 + 8 = 16.
* Этап 1. Разбиение чисел на K = 4 слов длины M = 4 бита:
**Дополним числа нулями до длины N = 16 и разобьём на слова. Получим две последовательности (в данном примере идентичные): {0001, 1000, 0000, 0000}. Или в десятичном виде: {3, 8, 0, 0} (слова расположены в порядке от младших бит к старшим). Все операции далее выполняются по модулю 2n+1, где n выбирается из условия n >= 2*M + log2(K) и является степенью двойки, то есть в нашем случае n >= 2*4 + 2 = 10, n = 16.

* Этап 2. Вычисление свёртки двух последовательностей:
** Этап 2.1 Умножение последовательностей на весовые коэффициенты. 2K-ый корень из единицы в кольце из 2n+1 элементов равен 22n/2K = 24 = 16. Результат умножения на вектор весовых коэффициентов в десятичном виде: {1, 128, 0, 0}.
**Этап 2.2 Прямое преобразование Фурье над кольцом из 2n элементов. K-ый корень из единицы в кольце вычетов по модулю 2n+1 равен 22n/2K = 24 = 256. Формула преобразования: <math>b_i = \sum^{K-1}_{j=0}{a_j\omega^{ij}}</math>, <math> \omega </math> - К-тый корень в кольце вычетов по модулю 2n+1. Результат: {129, 32769, 65410, 32770}.
**Этап 2.3 Покомпонентное произведение последовательностей в кольце из 2n элементов. Результат: {16641, 49153, 16129, 49155}.
**Этап 2.4 Обратное преобразование Фурье над кольцом из 2n элементов. Формула обратного преобразования: <math>b_i = \frac {1}{K} \sum^{K-1}_{j=0}{a_j\omega^{ij}}</math>, <math> \frac {1}{K} </math> - обратный к <math>K</math> по модулю 2n+1. Результат: {1, 256, 16384, 0}.
**Этап 2.5 Умножение последовательностей на коэффициенты, обратные к используемым в 2.1. Результат: {1, 16, 64, 0}.
* Этап 3. Выполнение операции переноса. Результат: {1, 0, 1, 4}. В двоичном виде: {0001, 0000, 0001, 0100}
* Этап 5. "Склейка" полученной последовательности, окончательный результат: 10100000010002 = 16641 = 129 * 129.
* Этап 5. Вычисление остатка по модулю 2N+1. В данном случае 16641 = 16641 mod 216+1

=== Описание ===
* [http://gmplib.org/manual/FFT-Multiplication.html#FFT-Multiplication Документация GMP LIB]
* [http://hal.archives-ouvertes.fr/docs/00/12/64/62/PDF/fft.pdf Статья A GMP-BASED IMPLEMENTATION OF SCHONHAGE-STRASSEN’S LARGE INTEGER MULTIPLICATION ALGORITHM]
* [http://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm Описание в английской Wikipedia]

=== Исходники ===
[http://www.ginac.de/CLN/cln.git/?p=cln.git;a=blob;f=src/base/digitseq/cl_DS_mul_fftm.h;h=a984802ce039e8179a26fda2681a2e133dba6fbd;hb=refs/heads/master CLN - Class Library for Numbers]

[https://github.com/myachik/ss_mult/blob/master/ss_mult.cpp Вариант реализации алгоритма]

Метод умножения Шёнхаге — Штрассена (New)

2014-01-13T13:12:48Z

Myachikov:

'''Алгоритм Шёнхаге - Штрассена''' - асимптотически быстрый алгоритм умножения больших чисел. Он был разработан Арнольдом Шёнхаге и Фолькером Штрассеном в 1971. Сложность алгоритма в нотации "О - большого" составляет O(''N'' log ''N'' log log ''N''). Алгоритм
рекурсивно использует быстрое преобразование Фурье над кольцом из 22''n'' + 1 элементов, специальный тип дискретного преобразования Фурье - теоретико-числовое преобразование.

Алгоритм Шёнхаге - Штрассена являлся наиболее быстрым асимптотически известным методом умножения с 1971 до 2007, когда Фюрером был предложен более быстрый метод. Тем не менее, алгоритм Фюрера в настоящее время имеет преимущество только для астрономически больших чисел и не используется на практике.

На практике алгоритм Шёнхаге - Штрассена начинает превосходить в скорости старые методы, такие как алгоритм Карацубы и алгоритм Тоома-Кука для чисел между 2215 и 2217 (от 10,000 до 40,000 десятичных разрядов).

Приложения алгоритма Шёнхаге - Штрассена включают в себя как чисто математические задачи, такие как поиск простых чисел Мерсенна и вычисление цифр числа ''π'', так и практические, такие как вычисление подстановки Кронекера, в которой умножение многочленов с целыми коэффициентами можно эффективно заменить умножением больших чисел; это используется на практике by GMP-ECM для факторизации целых чисел методом эллиптических кривых Ленстры.

== Теория ==

=== Свертка ===
Допустим, что мы перемножаем числа 123 и 456 «в столбик», однако без выполнения переноса. Результат будет выглядеть так:

{|width=300
| || || 1 || 2 || 3
|-
| || ×|| 4 || 5 || 6
|-
|colspan=5|<hr />
|-
| || || 6 || 12 || 18
|-
| || 5 || 10 || 15 ||
|-
| 4 || 8 || 12 || ||
|-
|colspan=5|<hr />
|-
| 4 || 13 || 28 || 27 || 18
|}

Эта последовательность (4, 13, 28, 27, 18) называется ''ациклической'' или ''линейной свёрткой'' от последовательностей (1,2,3) и (4,5,6). Зная ациклическую свёртку двух последовательностей, рассчитать произведение несложно: достаточно выполнить перенос (например, в самом правом столбце, мы оставляем 8 и добавляем 1 к столбцу, содержащему 27). В нашем примере это приводит к результату 56088.

Есть и другие типы свёрток, которые могут быть полезны. Допустим, что входящие последовательности содержат ''n'' элементов (в примере — 3). Тогда результирующая линейная свёртка содержит ''n'' + ''n'' − 1 элементов; если мы возьмём самый правый столбец ''n'' элементов и добавим самый левый столбец ''n'' − 1 ', в результате мы получим циклическую свёртку:

{|width=300
| || 28 || 27 || 18
|-
| + || || 4 || 13
|-
|colspan=5|<hr />
|-
| || 28 || 31 || 31
|}

Если мы произведём перенос при циклическом свёртывании, результат будет тот же, что и при произведении чисел по модулю B''n'' − 1 (в данном примере это 103 − 1 = 999). Выполним перенос по (28, 31, 31) и получим 3141, при этом 3141 ≡ 56088 (mod 999).

Наоборот, если мы возьмём самый правый столбец ''n'' элементов и ''вычтем'' самый левый столбец ''n''−1 элементов, то в результате мы получим обратную свёртку:

{|width=300
| || 28 || 27 || 18
|-
| − || || 4 || 13
|-
|colspan=5|<hr />
|-
| || 28 || 23 || 5
|}

Если мы произведём перенос при обратном свёртывании, то результат будет тот же, что и при произведении чисел по модулю B''n'' + 1. В данном примере, 103 + 1 = 1001, выполним перенос по (28, 23, 5) и получим 3035, при этом 3035 ≡ 56088 (mod 1001). Обратная свёртка может содержать отрицательные числа, которые могут быть убраны во время переноса, используя ту же технику, что и при длинных вычитаниях.

=== Теорема о свёртке ===

Как и другие методы, основанные на быстром преобразовании Фурье, алгоритм Шёнхаге - Штрассена в корне зависит от теоремы о свёртке, которая обеспечивает эффективный способ посчитать циклическую свёртку двух последовательностей. Её идея состоит в следующем:

: Циклическая свёртка двух векторов может быть найдена через [[дискретное преобразование Фурье]] (ДПФ) каждого из них, путём произведения результирующих векторов элемент за элементом, с последующим обратным преобразованием Фурье (ОДПФ).

Или через формулы:

: CyclicConvolution(X, Y) = IDFT(DFT(X) · DFT(Y)), где:
:: CyclicConvolution — ''циклическая свертка'',
:: DFT — ''дискретное преобразование Фурье'',
:: IDFT — ''обратное дискретное преобразование Фурье''.

Если мы посчитаем ДПФ и ОДПФ, используя быстрое преобразование Фурье, и вызовем наш алгоритм перемножения рекурсивно, чтобы перемножить входы(?) преобразованных векторов DFT(X) и DFT(Y), то в результате мы получим эффективный алгоритм для расчёта циклической свёртки.

В этом алгоритме, гораздо эффективней считать ''обратную циклическую'' свёртку; как оказывается, немного модифицированная версия теоремы о свёртке может позволить и это. Предположим, что вектора X и Y имеют длину ''n'', и ''a'' — примитивный корень порядка 2''n'' (это означает, что ''a''2''n'' = 1 и все меньшие степени ''a'' не равны 1). Таким образом мы можем определить третий вектор ''A'', называемый ''вектор веса'', обладающий следующими свойствами:
[[Файл:DIT-FFT-butterfly.png|мини|Операция «[[бабочка (БПФ)|бабочка]]».]]

: ''A'' = (''a''''j''), 0 ≤ ''j'' < ''n''
: ''A''−1 = (''a''−j''), 0 ≤ ''j'' < ''n''

Теперь мы можем записать:

: NegacyclicConvolution(''X'', ''Y'') = ''A''−1 · IDFT(DFT(''A'' · ''X'') · DFT(''A'' · ''Y'')), где
:: NegacyclicConvolution — ''Обратная циклическая свертка'',
:: DFT — ''дискретное преобразование Фурье'',
:: IDFT — ''обратное дискретное преобразование Фурье''.

Другими словами, это то же самое за исключением того, что входящие векторы умножены на ''A'', а результат умножен на ''A''−1.

=== Выбор кольца ===

Дискретное преобразование Фурье — абстрактная операция, которая может быть выполнена в любом алгебраическом [[Кольцо (математика)|кольце]]; обычно оно берётся из поля комплексных чисел, но фактически использовать комплексную арифметику с достаточной точностью, чтобы обеспечить точные результаты, медленно и неэффективно. Вместо этого мы можем использовать теоретико-числовое преобразование, которое производит преобразование в поле целых чисел по модулю N для некоторого целого N.

Так же как есть примитивные корни единицы любого порядка на комплексной плоскости, при любом заданном ''n'' мы можем выбрать подходящее N такое, что ''b'' — примитивный корень единицы порядка ''n'' в поле целых чисел по модулю N (другими словами, ''b''''n'' ≡ 1 (mod N), и все меньшие степени ''b'' не равны 1 mod N).

Алгоритм тратит большую часть времени на рекурсивное выполнение произведения меньших чисел; в простом варианте алгоритма это происходит в ряде мест:

# Внутри алгоритма быстрого преобразования Фурье, примитивный корень единицы ''b'' неоднократно возводится в степень и умножается на другие числа.
# При возведении в степень примитивного корня единицы ''a'' для получения вектора веса A с последующим умножением векторов A или A−1 на другие вектора.
# При выполнении последовательного перемножения преобразованных векторов.

Ключевой момент — выбрать N, модуль, равный 2''n'' + 1 для некоторого целого ''n''. У этого способа есть ряд преимуществ в ряде стандартных систем, в которых большие целые числа представлены в двоичном виде:

* Любое число может быть быстро уменьшено по модулю 2''n'' + 1 используя только сдвиг и сложение.
* Любые примитивные корни единицы в этом кольце могут быть записаны в форме 2''k''; соответственно мы можем умножать или делить любое число на корень из единицы используя сдвиг.
* Поэлементное рекурсивное перемножение преобразованный векторов может быть выполнено, используя обратную свёртку, которая работает быстрее, чем ациклическая свёртка, и в которой уже есть уменьшение результата
по модулю 2''n'' + 1.

=== Пример ===
Следующий пример иллюстрирует работу алгоритма.
Положим задано два числа: X и Y, каждое длиной 8 бит. Требуется найти значение RES = (X*Y) mod 2N+1
Возьмём для определённости X = Y = 12910 = 1000 00012.
Выберем N = 8 + 8 = 16.
* Этап 1. Разбиение чисел на K = 4 слов длины M = 4 бита:
**Дополним числа нулями до длины N = 16 и разобьём на слова. Получим две последовательности (в данном примере идентичные): {0001, 1000, 0000, 0000}. Или в десятичном виде: {3, 8, 0, 0} (слова расположены в порядке от младших бит к старшим). Все операции далее выполняются по модулю 2n+1, где n выбирается из условия n >= 2*M + log2(K) и является степенью двойки, то есть в нашем случае n >= 2*4 + 2 = 10, n = 16.

* Этап 2. Вычисление свёртки двух последовательностей:
** Этап 2.1 Умножение последовательностей на весовые коэффициенты. 2K-ый корень из единицы в кольце из 2n+1 элементов равен 22n/2K = 24 = 16. Результат умножения на вектор весовых коэффициентов в десятичном виде: {1, 128, 0, 0}.
**Этап 2.2 Прямое преобразование Фурье над кольцом из 2n элементов. Результат: {129, 32769, 65410, 32770}.
**Этап 2.3 Покомпонентное произведение последовательностей в кольце из 2n элементов. Результат: {16641, 49153, 16129, 49155}.
**Этап 2.4 Обратное преобразование Фурье над кольцом из 2n элементов. Результат: {1, 256, 16384, 0}.
**Этап 2.5 Умножение последовательностей на коэффициенты, обратные к используемым в 2.1. Результат: {1, 16, 64, 0}.
* Этап 3. Выполнение операции переноса. Результат: {1, 0, 1, 4}. В двоичном виде: {0001, 0000, 0001, 0100}
* Этап 5. "Склейка" полученной последовательности, окончательный результат: 10100000010002 = 16641 = 129 * 129.
* Этап 5. Вычисление остатка по модулю 2N+1. В данном случае 16641 = 16641 mod 216+1

=== Описание ===
* [http://gmplib.org/manual/FFT-Multiplication.html#FFT-Multiplication Документация GMP LIB]
* [http://hal.archives-ouvertes.fr/docs/00/12/64/62/PDF/fft.pdf Статья A GMP-BASED IMPLEMENTATION OF SCHONHAGE-STRASSEN’S LARGE INTEGER MULTIPLICATION ALGORITHM]
* [http://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm Описание в английской Wikipedia]

=== Исходники ===
[http://www.ginac.de/CLN/cln.git/?p=cln.git;a=blob;f=src/base/digitseq/cl_DS_mul_fftm.h;h=a984802ce039e8179a26fda2681a2e133dba6fbd;hb=refs/heads/master CLN - Class Library for Numbers]

[https://github.com/myachik/ss_mult/blob/master/ss_mult.cpp Вариант реализации алгоритма]

Метод умножения Шёнхаге — Штрассена (New)

2014-01-13T11:31:08Z

Myachikov:

'''Алгоритм Шёнхаге - Штрассена''' - асимптотически быстрый алгоритм умножения больших чисел. Он был разработан Арнольдом Шёнхаге и Волкером Штрассеном в 1971. Сложность алгоритма в нотации "О - большого" составляет O(''N'' log ''N'' log log ''N''). Алгоритм
рекурсивно использует быстрое преобразование Фурье над кольцом из 22''n'' + 1 элементов, специальный тип дискретного преобразования Фурье - теоретико-числовое преобразование.

Алгоритм Шёнхаге - Штрассена являлся наиболее быстрым асимптотически известным методом умножения с 1971 до 2007, когда Фюрером был предложен более быстрый метод. Тем не менее, алгоритм Фюрера в настоящее время имеет преимущество только для астрономически больших чисел и не используется на практике.

На практике алгоритм Шёнхаге - Штрассена начинает превосходить в скорости старые методы, такие как алгоритм Карацубы и алгоритм Тоома-Кука для чисел между 2215 и 2217 (от 10,000 до 40,000 десятичных разрядов).

Приложения алгоритма Шёнхаге - Штрассена включают в себя как чисто математические задачи, такие как поиск простых чисел Мерсенна и вычисление цифр числа ''π'', так и практические, такие как вычисление подстановки Кронекера, в которой умножение многочленов с целыми коэффициентами можно эффективно заменить умножением больших чисел; это используется на практике by GMP-ECM для факторизации целых чисел методом эллиптических кривых Ленстры.

== Теория ==

=== Свертка ===
Допустим, что мы перемножаем числа 123 и 456 «в столбик», однако без выполнения переноса. Результат будет выглядеть так:

{|width=300
| || || 1 || 2 || 3
|-
| || ×|| 4 || 5 || 6
|-
|colspan=5|<hr />
|-
| || || 6 || 12 || 18
|-
| || 5 || 10 || 15 ||
|-
| 4 || 8 || 12 || ||
|-
|colspan=5|<hr />
|-
| 4 || 13 || 28 || 27 || 18
|}

Эта последовательность (4, 13, 28, 27, 18) называется ''ациклической'' или ''линейной свёрткой'' от последовательностей (1,2,3) и (4,5,6). Зная ациклическую свёртку двух последовательностей, рассчитать произведение несложно: достаточно выполнить перенос (например, в самом правом столбце, мы оставляем 8 и добавляем 1 к столбцу, содержащему 27). В нашем примере это приводит к результату 56088.

Есть и другие типы свёрток, которые могут быть полезны. Допустим, что входящие последовательности содержат ''n'' элементов (в примере — 3). Тогда результирующая линейная свёртка содержит ''n'' + ''n'' − 1 элементов; если мы возьмём самый правый столбец ''n'' элементов и добавим самый левый столбец ''n'' − 1 ', в результате мы получим циклическую свёртку:

{|width=300
| || 28 || 27 || 18
|-
| + || || 4 || 13
|-
|colspan=5|<hr />
|-
| || 28 || 31 || 31
|}

Если мы произведём перенос при циклическом свёртывании, результат будет тот же, что и при произведении чисел по модулю B''n'' − 1 (в данном примере это 103 − 1 = 999). Выполним перенос по (28, 31, 31) и получим 3141, при этом 3141 ≡ 56088 (mod 999).

Наоборот, если мы возьмём самый правый столбец ''n'' элементов и ''вычтем'' самый левый столбец ''n''−1 элементов, то в результате мы получим обратную свёртку:

{|width=300
| || 28 || 27 || 18
|-
| − || || 4 || 13
|-
|colspan=5|<hr />
|-
| || 28 || 23 || 5
|}

Если мы произведём перенос при обратном свёртывании, то результат будет тот же, что и при произведении чисел по модулю B''n'' + 1. В данном примере, 103 + 1 = 1001, выполним перенос по (28, 23, 5) и получим 3035, при этом 3035 ≡ 56088 (mod 1001). Обратная свёртка может содержать отрицательные числа, которые могут быть убраны во время переноса, используя ту же технику, что и при длинных вычитаниях.

=== Теорема о свёртке ===

Как и другие методы, основанные на быстром преобразовании Фурье, алгоритм Шёнхаге - Штрассена в корне зависит от теоремы о свёртке, которая обеспечивает эффективный способ посчитать циклическую свёртку двух последовательностей. Её идея состоит в следующем:

: Циклическая свёртка двух векторов может быть найдена через [[дискретное преобразование Фурье]] (ДПФ) каждого из них, путём произведения результирующих векторов элемент за элементом, с последующим обратным преобразованием Фурье (ОДПФ).

Или через формулы:

: CyclicConvolution(X, Y) = IDFT(DFT(X) · DFT(Y)), где:
:: CyclicConvolution — ''циклическая свертка'',
:: DFT — ''дискретное преобразование Фурье'',
:: IDFT — ''обратное дискретное преобразование Фурье''.

Если мы посчитаем ДПФ и ОДПФ, используя быстрое преобразование Фурье, и вызовем наш алгоритм перемножения рекурсивно, чтобы перемножить входы(?) преобразованных векторов DFT(X) и DFT(Y), то в результате мы получим эффективный алгоритм для расчёта циклической свёртки.

В этом алгоритме, гораздо эффективней считать ''обратную циклическую'' свёртку; как оказывается, немного модифицированная версия теоремы о свёртке может позволить и это. Предположим, что вектора X и Y имеют длину ''n'', и ''a'' — примитивный корень порядка 2''n'' (это означает, что ''a''2''n'' = 1 и все меньшие степени ''a'' не равны 1). Таким образом мы можем определить третий вектор ''A'', называемый ''вектор веса'', обладающий следующими свойствами:
[[Файл:DIT-FFT-butterfly.png|мини|Операция «[[бабочка (БПФ)|бабочка]]».]]

: ''A'' = (''a''''j''), 0 ≤ ''j'' < ''n''
: ''A''−1 = (''a''−j''), 0 ≤ ''j'' < ''n''

Теперь мы можем записать:

: NegacyclicConvolution(''X'', ''Y'') = ''A''−1 · IDFT(DFT(''A'' · ''X'') · DFT(''A'' · ''Y'')), где
:: NegacyclicConvolution — ''Обратная циклическая свертка'',
:: DFT — ''дискретное преобразование Фурье'',
:: IDFT — ''обратное дискретное преобразование Фурье''.

Другими словами, это то же самое за исключением того, что входящие векторы умножены на ''A'', а результат умножен на ''A''−1.

=== Выбор кольца ===

Дискретное преобразование Фурье — абстрактная операция, которая может быть выполнена в любом алгебраическом [[Кольцо (математика)|кольце]]; обычно оно берётся из поля комплексных чисел, но фактически использовать комплексную арифметику с достаточной точностью, чтобы обеспечить точные результаты, медленно и неэффективно. Вместо этого мы можем использовать теоретико-числовое преобразование, которое производит преобразование в поле целых чисел по модулю N для некоторого целого N.

Так же как есть примитивные корни единицы любого порядка на комплексной плоскости, при любом заданном ''n'' мы можем выбрать подходящее N такое, что ''b'' — примитивный корень единицы порядка ''n'' в поле целых чисел по модулю N (другими словами, ''b''''n'' ≡ 1 (mod N), и все меньшие степени ''b'' не равны 1 mod N).

Алгоритм тратит большую часть времени на рекурсивное выполнение произведения меньших чисел; в простом варианте алгоритма это происходит в ряде мест:

# Внутри алгоритма быстрого преобразования Фурье, примитивный корень единицы ''b'' неоднократно возводится в степень и умножается на другие числа.
# При возведении в степень примитивного корня единицы ''a'' для получения вектора веса A с последующим умножением векторов A или A−1 на другие вектора.
# При выполнении последовательного перемножения преобразованных векторов.

Ключевой момент — выбрать N, модуль, равный 2''n'' + 1 для некоторого целого ''n''. У этого способа есть ряд преимуществ в ряде стандартных систем, в которых большие целые числа представлены в двоичном виде:

* Любое число может быть быстро уменьшено по модулю 2''n'' + 1 используя только сдвиг и сложение.
* Любые примитивные корни единицы в этом кольце могут быть записаны в форме 2''k''; соответственно мы можем умножать или делить любое число на корень из единицы используя сдвиг.
* Поэлементное рекурсивное перемножение преобразованный векторов может быть выполнено, используя обратную свёртку, которая работает быстрее, чем ациклическая свёртка, и в которой уже есть уменьшение результата
по модулю 2''n'' + 1.

=== Пример ===
Следующий пример иллюстрирует работу алгоритма.
Положим задано два числа: X и Y, каждое длиной 8 бит. Требуется найти значение RES = (X*Y) mod 2N+1
Возьмём для определённости X = Y = 12910 = 1000 00012.
Выберем N = 8 + 8 = 16.
* Этап 1. Разбиение чисел на K = 4 слов длины M = 4 бита:
**Дополним числа нулями до длины N = 16 и разобьём на слова. Получим две последовательности (в данном примере идентичные): {0001, 1000, 0000, 0000}. Или в десятичном виде: {3, 8, 0, 0} (слова расположены в порядке от младших бит к старшим). Все операции далее выполняются по модулю 2n+1, где n выбирается из условия n >= 2*M + log2(K) и является степенью двойки, то есть в нашем случае n >= 2*4 + 2 = 10, n = 16.

* Этап 2. Вычисление свёртки двух последовательностей:
** Этап 2.1 Умножение последовательностей на весовые коэффициенты. 2K-ый корень из единицы в кольце из 2n+1 элементов равен 22n/2K = 24 = 16. Результат умножения на вектор весовых коэффициентов в десятичном виде: {1, 128, 0, 0}.
**Этап 2.2 Прямое преобразование Фурье над кольцом из 2n элементов. Результат: {129, 32769, 65410, 32770}.
**Этап 2.3 Покомпонентное произведение последовательностей в кольце из 2n элементов. Результат: {16641, 49153, 16129, 49155}.
**Этап 2.4 Обратное преобразование Фурье над кольцом из 2n элементов. Результат: {1, 256, 16384, 0}.
**Этап 2.5 Умножение последовательностей на коэффициенты, обратные к используемым в 2.1. Результат: {1, 16, 64, 0}.
* Этап 3. Выполнение операции переноса. Результат: {1, 0, 1, 4}. В двоичном виде: {0001, 0000, 0001, 0100}
* Этап 5. "Склейка" полученной последовательности, окончательный результат: 10100000010002 = 16641 = 129 * 129.
* Этап 5. Вычисление остатка по модулю 2N+1. В данном случае 16641 = 16641 mod 216+1

=== Описание ===
* [http://gmplib.org/manual/FFT-Multiplication.html#FFT-Multiplication Документация GMP LIB]
* [http://hal.archives-ouvertes.fr/docs/00/12/64/62/PDF/fft.pdf Статья A GMP-BASED IMPLEMENTATION OF SCHONHAGE-STRASSEN’S LARGE INTEGER MULTIPLICATION ALGORITHM]
* [http://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm Описание в английской Wikipedia]

=== Исходники ===
[http://www.ginac.de/CLN/cln.git/?p=cln.git;a=blob;f=src/base/digitseq/cl_DS_mul_fftm.h;h=a984802ce039e8179a26fda2681a2e133dba6fbd;hb=refs/heads/master CLN - Class Library for Numbers]

[https://github.com/myachik/ss_mult/blob/master/ss_mult.cpp Вариант реализации алгоритма]

Метод умножения Шёнхаге — Штрассена (New)

2014-01-13T11:30:29Z

Myachikov:

'''Алгоритм Шёнхаге - Штрассена''' - асимптотически быстрый алгоритм умножения больших чисел. Он был разработан Арнольдом Шёнхаге и Волкером Штрассеном в 1971. Сложность алгоритма в нотации "О - большого" составляет O(''N'' log ''N'' log log ''N''). Алгоритм
рекурсивно использует быстрое преобразование Фурье над кольцом из 22''n'' + 1 элементов, специальный тип дискретного преобразования Фурье - теоретико-числовое преобразование.

Алгоритм Шёнхаге - Штрассена являлся наиболее быстрым асимптотически известным методом умножения с 1971 до 2007, когда Фюрером был предложен более быстрый метод. Тем не менее, алгоритм Фюрера в настоящее время имеет преимущество только для астрономически больших чисел и не используется на практике.

На практике алгоритм Шёнхаге - Штрассена начинает превосходить в скорости старые методы, такие как алгоритм Карацубы и алгоритм Тоома-Кука для чисел между 2215 и 2217 (от 10,000 до 40,000 десятичных разрядов).

Приложения алгоритма Шёнхаге - Штрассена включают в себя как чисто математические задачи, такие как поиск простых чисел Мерсенна и вычисление цифр числа ''π'', так и практические, такие как вычисление подстановки Кронекера, в которой умножение многочленов с целыми коэффициентами можно эффективно заменить умножением больших чисел; это используется на практике by GMP-ECM для факторизации целых чисел методом эллиптических кривых Ленстры.

== Теория ==

=== Свертка ===
Допустим, что мы перемножаем числа 123 и 456 «в столбик», однако без выполнения переноса. Результат будет выглядеть так:

{|width=300
| || || 1 || 2 || 3
|-
| || ×|| 4 || 5 || 6
|-
|colspan=5|<hr />
|-
| || || 6 || 12 || 18
|-
| || 5 || 10 || 15 ||
|-
| 4 || 8 || 12 || ||
|-
|colspan=5|<hr />
|-
| 4 || 13 || 28 || 27 || 18
|}

Эта последовательность (4, 13, 28, 27, 18) называется ''ациклической'' или ''линейной свёрткой'' от последовательностей (1,2,3) и (4,5,6). Зная ациклическую свёртку двух последовательностей, рассчитать произведение несложно: достаточно выполнить перенос (например, в самом правом столбце, мы оставляем 8 и добавляем 1 к столбцу, содержащему 27). В нашем примере это приводит к результату 56088.

Есть и другие типы свёрток, которые могут быть полезны. Допустим, что входящие последовательности содержат ''n'' элементов (в примере — 3). Тогда результирующая линейная свёртка содержит ''n'' + ''n'' − 1 элементов; если мы возьмём самый правый столбец ''n'' элементов и добавим самый левый столбец ''n'' − 1 ', в результате мы получим циклическую свёртку:

{|width=300
| || 28 || 27 || 18
|-
| + || || 4 || 13
|-
|colspan=5|<hr />
|-
| || 28 || 31 || 31
|}

Если мы произведём перенос при циклическом свёртывании, результат будет тот же, что и при произведении чисел по модулю B''n'' − 1 (в данном примере это 103 − 1 = 999). Выполним перенос по (28, 31, 31) и получим 3141, при этом 3141 ≡ 56088 (mod 999).

Наоборот, если мы возьмём самый правый столбец ''n'' элементов и ''вычтем'' самый левый столбец ''n''−1 элементов, то в результате мы получим обратную свёртку:

{|width=300
| || 28 || 27 || 18
|-
| − || || 4 || 13
|-
|colspan=5|<hr />
|-
| || 28 || 23 || 5
|}

Если мы произведём перенос при обратном свёртывании, то результат будет тот же, что и при произведении чисел по модулю B''n'' + 1. В данном примере, 103 + 1 = 1001, выполним перенос по (28, 23, 5) и получим 3035, при этом 3035 ≡ 56088 (mod 1001). Обратная свёртка может содержать отрицательные числа, которые могут быть убраны во время переноса, используя ту же технику, что и при длинных вычитаниях.

=== Теорема о свёртке ===

Как и другие методы, основанные на быстром преобразовании Фурье, алгоритм Шёнхаге - Штрассена в корне зависит от теоремы о свёртке, которая обеспечивает эффективный способ посчитать циклическую свёртку двух последовательностей. Её идея состоит в следующем:

: Циклическая свёртка двух векторов может быть найдена через [[дискретное преобразование Фурье]] (ДПФ) каждого из них, путём произведения результирующих векторов элемент за элементом, с последующим обратным преобразованием Фурье (ОДПФ).

Или через формулы:

: CyclicConvolution(X, Y) = IDFT(DFT(X) · DFT(Y)), где:
:: CyclicConvolution — ''циклическая свертка'',
:: DFT — ''дискретное преобразование Фурье'',
:: IDFT — ''обратное дискретное преобразование Фурье''.

Если мы посчитаем ДПФ и ОДПФ, используя быстрое преобразование Фурье, и вызовем наш алгоритм перемножения рекурсивно, чтобы перемножить входы(?) преобразованных векторов DFT(X) и DFT(Y), то в результате мы получим эффективный алгоритм для расчёта циклической свёртки.

В этом алгоритме, гораздо эффективней считать ''обратную циклическую'' свёртку; как оказывается, немного модифицированная версия теоремы о свёртке может позволить и это. Предположим, что вектора X и Y имеют длину ''n'', и ''a'' — примитивный корень порядка 2''n'' (это означает, что ''a''2''n'' = 1 и все меньшие степени ''a'' не равны 1). Таким образом мы можем определить третий вектор ''A'', называемый ''вектор веса'', обладающий следующими свойствами:
[[Файл:DIT-FFT-butterfly.png|мини|Операция «[[бабочка (БПФ)|бабочка]]».]]

: ''A'' = (''a''''j''), 0 ≤ ''j'' < ''n''
: ''A''−1 = (''a''−j''), 0 ≤ ''j'' < ''n''

Теперь мы можем записать:

: NegacyclicConvolution(''X'', ''Y'') = ''A''−1 · IDFT(DFT(''A'' · ''X'') · DFT(''A'' · ''Y'')), где
:: NegacyclicConvolution — ''Обратная циклическая свертка'',
:: DFT — ''дискретное преобразование Фурье'',
:: IDFT — ''обратное дискретное преобразование Фурье''.

Другими словами, это то же самое за исключением того, что входящие векторы умножены на ''A'', а результат умножен на ''A''−1.

=== Выбор кольца ===

Дискретное преобразование Фурье — абстрактная операция, которая может быть выполнена в любом алгебраическом [[Кольцо (математика)|кольце]]; обычно оно берётся из поля комплексных чисел, но фактически использовать комплексную арифметику с достаточной точностью, чтобы обеспечить точные результаты, медленно и неэффективно. Вместо этого мы можем использовать теоретико-числовое преобразование, которое производит преобразование в поле целых чисел по модулю N для некоторого целого N.

Так же как есть примитивные корни единицы любого порядка на комплексной плоскости, при любом заданном ''n'' мы можем выбрать подходящее N такое, что ''b'' — примитивный корень единицы порядка ''n'' в поле целых чисел по модулю N (другими словами, ''b''''n'' ≡ 1 (mod N), и все меньшие степени ''b'' не равны 1 mod N).

Алгоритм тратит большую часть времени на рекурсивное выполнение произведения меньших чисел; в простом варианте алгоритма это происходит в ряде мест:

# Внутри алгоритма быстрого преобразования Фурье, примитивный корень единицы ''b'' неоднократно возводится в степень и умножается на другие числа.
# При возведении в степень примитивного корня единицы ''a'' для получения вектора веса A с последующим умножением векторов A или A−1 на другие вектора.
# При выполнении последовательного перемножения преобразованных векторов.

Ключевой момент — выбрать N, модуль, равный 2''n'' + 1 для некоторого целого ''n''. У этого способа есть ряд преимуществ в ряде стандартных систем, в которых большие целые числа представлены в двоичном виде:

* Любое число может быть быстро уменьшено по модулю 2''n'' + 1 используя только сдвиг и сложение.
* Любые примитивные корни единицы в этом кольце могут быть записаны в форме 2''k''; соответственно мы можем умножать или делить любое число на корень из единицы используя сдвиг.
* Поэлементное рекурсивное перемножение преобразованный векторов может быть выполнено, используя обратную свёртку, которая работает быстрее, чем ациклическая свёртка, и в которой уже есть уменьшение результата
по модулю 2''n'' + 1.

=== Пример ===
Следующий пример иллюстрирует работу алгоритма.
Положим задано два числа: X и Y, каждое длиной 8 бит. Требуется найти значение RES = (X*Y) mod 2N+1
Возьмём для определённости X = Y = 12910 = 1000 00012.
Выберем N = 8 + 8 = 16.
* Этап 1. Разбиение чисел на K = 4 слов длины M = 4 бита:
**Дополним числа нулями до длины N = 16 и разобьём на слова. Получим две последовательности (в данном примере идентичные): {0001, 1000, 0000, 0000}. Или в десятичном виде: {3, 8, 0, 0} (слова расположены в порядке от младших бит к старшим). Все операции далее выполняются по модулю 2n+1, где n выбирается из условия n >= 2*M + log2(K) и является степенью двойки, то есть в нашем случае n >= 2*4 + 2 = 10, n = 16.

* Этап 2. Вычисление свёртки двух последовательностей:
** Этап 2.1 Умножение последовательностей на весовые коэффициенты. 2K-ый корень из единицы в кольце из 2n+1 элементов равен 22n/2K = 24 = 16. Результат умножения на вектор весовых коэффициентов в десятичном виде: {1, 128, 0, 0}.
**Этап 2.2 Прямое преобразование Фурье над кольцом из 2n элементов. Результат: {129, 32769, 65410, 32770}.
**Этап 2.3 Покомпонентное произведение последовательностей в кольце из 2n элементов. Результат: {16641, 49153, 16129, 49155}.
**Этап 2.4 Обратное преобразование Фурье над кольцом из 2n элементов. Результат: {1, 256, 16384, 0}.
**Этап 2.5 Умножение последовательностей на коэффициенты, обратные к используемым в 2.1. Результат: {1, 16, 64, 0}.
* Этап 3. Выполнение операции переноса. Результат: {1, 0, 1, 4}. В двоичном виде: {0001, 0000, 0001, 0100}
* Этап 5. "Склейка" полученной последовательности, получение результата: 10100000010002 = 16641 = 129 * 129.
* Этап 5. Вычисление остатка по модулю 2N+1. В данном случае 16641 = 16641 mod 216+1

=== Описание ===
* [http://gmplib.org/manual/FFT-Multiplication.html#FFT-Multiplication Документация GMP LIB]
* [http://hal.archives-ouvertes.fr/docs/00/12/64/62/PDF/fft.pdf Статья A GMP-BASED IMPLEMENTATION OF SCHONHAGE-STRASSEN’S LARGE INTEGER MULTIPLICATION ALGORITHM]
* [http://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm Описание в английской Wikipedia]

=== Исходники ===
[http://www.ginac.de/CLN/cln.git/?p=cln.git;a=blob;f=src/base/digitseq/cl_DS_mul_fftm.h;h=a984802ce039e8179a26fda2681a2e133dba6fbd;hb=refs/heads/master CLN - Class Library for Numbers]

[https://github.com/myachik/ss_mult/blob/master/ss_mult.cpp Вариант реализации алгоритма]

Метод умножения Шёнхаге — Штрассена (New)

2014-01-13T09:25:23Z

Myachikov: Новая страница: «'''Алгоритм Шёнхаге - Штрассена''' - асимптотически быстрый алгоритм умножения больших чис…»

'''Алгоритм Шёнхаге - Штрассена''' - асимптотически быстрый алгоритм умножения больших чисел. Он был разработан Арнольдом Шёнхаге и Волкером Штрассеном в 1971. Сложность алгоритма в нотации "О - большого" составляет O(''N'' log ''N'' log log ''N''). Алгоритм
рекурсивно использует быстрое преобразование Фурье над кольцом из 22''n'' + 1 элементов, специальный тип дискретного преобразования Фурье - теоретико-числовое преобразование.

Алгоритм Шёнхаге - Штрассена являлся наиболее быстрым асимптотически известным методом умножения с 1971 до 2007, когда Фюрером был предложен более быстрый метод. Тем не менее, алгоритм Фюрера в настоящее время имеет преимущество только для астрономически больших чисел и не используется на практике.

На практике алгоритм Шёнхаге - Штрассена начинает превосходить в скорости старые методы, такие как алгоритм Карацубы и алгоритм Тоома-Кука для чисел между 2215 и 2217 (от 10,000 до 40,000 десятичных разрядов).

== Теория ==

=== Свертка ===
Допустим, что мы перемножаем числа 123 и 456 «в столбик», однако без выполнения переноса. Результат будет выглядеть так:

{|width=300
| || || 1 || 2 || 3
|-
| || ×|| 4 || 5 || 6
|-
|colspan=5|<hr />
|-
| || || 6 || 12 || 18
|-
| || 5 || 10 || 15 ||
|-
| 4 || 8 || 12 || ||
|-
|colspan=5|<hr />
|-
| 4 || 13 || 28 || 27 || 18
|}

Эта последовательность (4, 13, 28, 27, 18) называется ''ациклической'' или ''линейной свёрткой'' от последовательностей (1,2,3) и (4,5,6). Зная ациклическую свёртку двух последовательностей, рассчитать произведение несложно: достаточно выполнить перенос (например, в самом правом столбце, мы оставляем 8 и добавляем 1 к столбцу, содержащему 27). В нашем примере это приводит к результату 56088.

Есть и другие типы свёрток, которые могут быть полезны. Допустим, что входящие последовательности содержат ''n'' элементов (в примере — 3). Тогда результирующая линейная свёртка содержит ''n'' + ''n'' − 1 элементов; если мы возьмём самый правый столбец ''n'' элементов и добавим самый левый стоблец ''n'' − 1 ', в результате мы получим циклическую свёртку:

{|width=300
| || 28 || 27 || 18
|-
| + || || 4 || 13
|-
|colspan=5|<hr />
|-
| || 28 || 31 || 31
|}

Если мы произведём перенос при циклическом свёртывании, результат будет тот же, что и при произведении чисел по модулю B''n'' − 1 (в данном примере это 103 − 1 = 999). Выполним перенос по (28, 31, 31) и получим 3141, при этом 3141 ≡ 56088 (mod 999).

Наоборот, если мы возьмём самый правый столбец ''n'' элементов и ''вычтем'' самый левый столбец ''n''−1 элементов, то в результате мы получим обратную свёртку:

{|width=300
| || 28 || 27 || 18
|-
| − || || 4 || 13
|-
|colspan=5|<hr />
|-
| || 28 || 23 || 5
|}

Если мы произведём перенос при обратном свёртывании, то результат будет тот же, что и при произведении чисел по модулю B''n'' + 1. В данном примере, 103 + 1 = 1001, выполним перенос по (28, 23, 5) и получим 3035, при этом 3035 ≡ 56088 (mod 1001). Обратная свёртка может содержать отрицательные числа, которые могут быть убраны во время переноса, используя ту же технику, что и при длинных вычитаниях.

=== Теорема о свёртке ===

Как и другие методы, основанные на быстром преобразовании Фурье, алгоритм Фюрера в корне зависит от теоремы о свёртке, которая обеспечивает эффективный способ посчитать циклическую свёртку двух последовательностей. Её идея состоит в следующем:

: Циклическая свёртка двух векторов может быть найдена через [[дискретное преобразование Фурье]] (ДПФ) каждого из них, путём произведения результирующих векторов элемент за элементом, с последующим обратным преобразованием Фурье (ОДПФ).

Или через формулы:

: CyclicConvolution(X, Y) = IDFT(DFT(X) · DFT(Y)), где:
:: CyclicConvolution — ''циклическая свертка'',
:: DFT — ''дискретное преобразование Фурье'',
:: IDFT — ''обратное дискретное преобразование Фурье''.

Если мы посчитаем ДПФ и ОДПФ, используя быстрое преобразование Фурье, и вызовем наш алгоритм перемножения рекурсивно, чтобы перемножить входы(?) преобразованных векторов DFT(X) и DFT(Y), то в результате мы получим эффективный алгоритм для расчёта циклической свёртки.

В этом алгоритме, гораздо эффективней считать ''обратную циклическую'' свёртку; как оказывается, немного модифицированная версия теоремы о свёртке может позволить и это. Предположим, что вектора X и Y имеют длину ''n'', и ''a'' — примитивный корень порядка 2''n'' (это означает, что ''a''2''n'' = 1 и все меньшие степени ''a'' не равны 1). Таким образом мы можем определить третий вектор ''A'', называемый ''вектор веса'', обладающий следующими свойствами:
[[Файл:DIT-FFT-butterfly.png|мини|Операция «[[бабочка (БПФ)|бабочка]]».]]

: ''A'' = (''a''''j''), 0 ≤ ''j'' < ''n''
: ''A''−1 = (''a''−j''), 0 ≤ ''j'' < ''n''

Теперь мы можем записать:

: NegacyclicConvolution(''X'', ''Y'') = ''A''−1 · IDFT(DFT(''A'' · ''X'') · DFT(''A'' · ''Y'')), где
:: NegacyclicConvolution — ''Обратная циклическая сертка'',
:: DFT — ''дискретное преобразование Фурье'',
:: IDFT — ''обратное дискретное преобразование Фурье''.

Другими словами, это то же самое за исключением того, что входящие векторы умножены на ''A'', а результат умножен на ''A''−1.

=== Выбор кольца ===

Дискретное преобразование Фурье — абстрактная операция, которая может быть выполнена в любом алгебраическом [[Кольцо (математика)|кольце]]; обычно оно берётся из поля комплексных чисел, но фактически использовать комплексную арифметику с достаточной точностью, чтобы обеспечить точные результаты, медленно и неэффективно. Вместо этого мы можем использовать теоретико-числовое преобразование, которое производит преобразование в поле целых чисел по модулю N для некоторого целого N.

Так же как есть примитивные корни единицы любого порядка на комплексной плоскости, при любом заданном ''n'' мы можем выбрать подходящее N такое, что ''b'' — примитивный корень единицы порядка ''n'' в поле целых чисел по модулю N (другими словами, ''b''''n'' ≡ 1 (mod N), и все меньшие степени ''b'' не равны 1 mod N).

Алгоритм тратит большую часть времени на рекурсивное выполнение произведения меньших чисел; в простом варианте алгоритма это происходит в ряде мест:

# Внутри алгоритма быстрого преобразования Фурье, примитивный корень единицы ''b'' неоднократно возводится в степень и умножается на другие числа.
# При возведении в степень примитивного корня единицы ''a'' для получения вектора веса A с последующим умножением векторов A или A−1 на другие вектора.
# При выполнении последовательного перемножения преобразованных векторов.

Ключевой момент — выбрать N, модуль, равный 2''n'' + 1 для некоторого целого ''n''. У этого способа есть ряд преимуществ в ряде стандартных систем, в которых большие целые числа представлены в двоичном виде:

* Любое число может быть быстро уменьшено по модулю 2''n'' + 1 используя только сдвиг и сложение.
* Любые примитивные корни единицы в этом кольце могут быть записаны в форме 2''k''; соответственно мы можем умножать или делить любое число на корень из единицы используя сдвиг.
* Поэлементное рекурсивное перемножение преобразованный векторов может быть выполнено, используя обратную свёртку, которая работает быстрее, чем ациклическая свёртка, и в которой уже есть уменьшение результата
по модулю 2''n'' + 1.

=== Пример ===
Положим задано два числа: X и Y, каждое длиной 512 бит. Требуется найти значение RES = (X*Y) mod 2N+1
* Этап 1. Основные параметры алгоритма:
** len1 = 512
** len2 = 512
** N = 1024
** 2k ~ sqrt(N) = 32, k = 5
** n >= 2*N/2k+k = (2048/32)+5 = 69. Поскольку n должно делиться на 2k, то можно выбрать n = 128.
* Этап 2 (возникает на умножении DFT(X)*DFT(Y) mod 2128+1):
** len1 = 129
** len2 = 129
** N = 128
** 2k ~ (sqrt(N) = 11.3). Выберем k = 4.
** n >= 2*N/2k+k = (256/16)+4 = 20. Поскольку n должно делиться на 2k, то можно выбрать n = 32.
* Этап 3 (возникает на умножении DFT(X)*DFT(Y) mod 232+1):
** len1 = 33
** len2 = 33
это умножение можно делать аппаратно, не углубляюсь в рекурсию.

=== Описание ===
* [http://gmplib.org/manual/FFT-Multiplication.html#FFT-Multiplication Документация GMP LIB]
* [http://hal.archives-ouvertes.fr/docs/00/12/64/62/PDF/fft.pdf Статья A GMP-BASED IMPLEMENTATION OF SCHONHAGE-STRASSEN’S LARGE INTEGER MULTIPLICATION ALGORITHM]
* [http://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm Описание в английской Wikipedia]

=== Исходники ===
[http://www.ginac.de/CLN/cln.git/?p=cln.git;a=blob;f=src/base/digitseq/cl_DS_mul_fftm.h;h=a984802ce039e8179a26fda2681a2e133dba6fbd;hb=refs/heads/master CLN - Class Library for Numbers]

Схема коррекции ошибок Rajendra S. Katti

2013-07-31T12:53:54Z

Myachikov:

== Введение ==
Метод, описанный автором, основан на использовании избыточной системы остаточных классов. Недостатки представления чисел в такой системе с целью коррекции ошибок при вычислениях покажем на примере.

Пусть избыточная СОК определяется следующими модулями <math> m_1 = 3, m_2 = 4, m_3 = 5, m_4 = 7</math>, а динамический диапазон определяется <math>m_1</math> и <math>m_2</math> и соответственно равен <math>\{0,...,11\}</math>, а модули <math>m_3, m_4</math> являются дополнительными и необходимы для проверки на ошибки.

Натуральное число из динамического диапазона может быть восстановлено по любым трём из четырёх остатков по представленным модулям. Например, число <math>x=5=(2,1,0,2)</math> может быть представлено как <math>(2,1,0)</math> в системе модулей <math>\{m_1, m_2, m_3\}</math> или <math>(2,0,2), (1,0,2)</math> в системах <math>\{m_1, m_3, m_4\}</math> и <math>\{m_2, m_3, m_4\}</math> соответственно. Поэтому, если одна цифра отбрасывается в представлении <math>(2,1,0,2)</math>, правильный результат восстановления будет получен тогда и только тогда, когда отбрасываемая цифра ошибочна. Итак, обнаружение ошибок и коррекция выполняется следующим образом:
# Если число лежит вне динамического диапазона - произошла ошибка.
# Остатки отбрасываются по одному и если в результате одного из отбрасываний полученное число лежит в динамическом диапазоне, то отброшенный остаток является ошибочным.
# Корректное значение остатка получается расширением полученного на шаге 2 множества модулей, остатки по которым не содержат ошибок.

Такой подход требует много времени и вычислительных ресурсов. Ниже будет предложена новая схема обнаружения и коррекции ошибок на основе СОК.

== Модули, имеющие общие делители ==
Каждый набор остатков необязательно соответствует числу из динамического диапазона, если модули не являются попарно взаимно простыми. Легко показать, что такое соответствие есть тогда и только тогда, когда равенство
:<math>|r_i|_k = |r_j|_k</math>,
выполняется для любых <math>i, j</math>, где <math>k = gcd(r_i, r_j)</math>. Таким образом, ошибочный разряд может быть найден выполнением только простой операции вычисления остатка.

Приведём пример, демонстрирующий, как можно проверить набор остатков на наличие ошибок. Проверим набор <math>(2, 6, 30, 42)</math> в системе модулей <math>\{4, 15, 36, 48\}</math>. Наибольшие общие делители всевозможных пар модулей соответственно равны: <math>gcd(4, 15) = 1, gcd(4, 36) = 4, gcd(4, 48) = 4, gcd(15, 36) = 3, gcd(15, 48) = 3, gcd(36, 48) = 12</math>.
Используя описанное выше соотношение для проверки, получаем: <math>|2|_4=|30|_4, |2|_4=|42|_4, |6|_3 = |30|_3, |6|_3=|42|_3</math> и <math> |30|_{12}=|42|_{12}</math>. То есть ошибок нет и число может быть восстановлено с использованием китайской теоремы об остатках (КТО) для модулей, имеющих нетривиальные общие делители.

'''КТО для системы модулей, попарно не взаимно простых'''. Пусть есть представление <math>(r_1, ..., r_n)</math> натурального числа <math>x</math>, определяемое системой модулей <math>(m_1, ..., m_n)</math>. Тогда
:<math>|x|_M = \left | \sum^{n}_{i=1} {{{\alpha}_i}r_i} \right |</math>,
где <math> {\alpha}_i </math>- число, такое что <math>|{\alpha}_i|_M = 0, |{\alpha}_i|_{\frac{M}{{\mu}_i}} = 1 </math>,
где <math> {\mu}_i </math> - числа, такие что <math> M = \prod^n_{i=1}{{\mu}_i} </math>, <math> {\mu}_i </math> делит <math> m_i </math> и <math> M = lcm(m_1, ..., m_n) </math>. Если для некоторого <math>i</math> такое <math>{\alpha}_i</math> не существует, то соответствующее слагаемое в сумме опускается.

Операции над числами остаются корректными и в системе не взаимно простых попарно модулей.

== Алгоритм ==
Если в множестве модулей <math>(m_1, ..., m_n)</math> есть имеющие нетривиальный общий делитель, то каждый набор остатков <math>(r_1, ..., r_n)</math> необязательно соответствует некоторому числу. Каждый набор, для которого такое соответствие имеет место, будем называть '''кодовым словом'''.

'''Расстоянием''' между двумя представлениями чисел в СОК(кодовыми словами) будем считать количество разрядов, в которых кодовые слова различаются между собой.

'''Кодовым расстоянием''' назовём минимальное расстояние среди всех пар кодовых слов.

Чтобы обнаружить и исправить <math>d-1</math> или меньше ошибок, необходимо и достаточно, чтобы кодовое расстояние было не меньше <math>d</math>. Не более <math>t</math> ошибок может быть обнаружено и исправлено тогда и только тогда, когда кодовое расстояние больше, чем <math>(2t+1)</math>. Если кодовое расстояние больше или равно, чем <math>(t+d+1)</math>, то любая комбинация из <math>t</math> ошибок может быть исправлена и до <math>d</math>ошибок могут быть обнаружены.
Для множества из <math>n</math> модулей в соответствующей системе могут быть представлены числа <math>(0, ..., M-1)</math>, где <math>M = lcm(m_1, ..., m_n)</math>.

'''Циклическим числом''' <math>c_i</math> для каждого модуля <math>m_i</math> назовём <math>c_i = M/m_i</math>.
=== Модули с попарно взаимно простыми циклическими числами ===
Построим теперь систему модулей с попарно взаимно простыми циклическими числами, которая поможет получит эффективный алгоритм обнаружения и коррекции ошибок.

'''Построение 1'''.
#Выберем <math>n</math> попарно взаимно простых циклических чисел, не равных единице <math>(c_1, ..., c_n)</math>.
#Получим систему модулей <math>(m_1, ..., m_n)</math> следующим образом:
:<math>m_i = \frac{1}{c_i}\prod^n_{k=1}{c_k}</math>.

Полученные модули обладают двумя свойствами:
#Наименьшие общие кратны всевозможных пар модулей равны между собой.
#Наибольшие общие делители всевозможных пар модулей различны.

'''Теорема 1'''. Для построенной системы модулей кодовое расстояние равно <math>n-1</math>, где <math>n</math> - количество модулей.

==== Алгоритм обнаружения ошибки ====
#Если <math>r_i \geqslant m_i</math>, то произошла ошибка в <math>i</math>-том разряде.
#Если не выполняется хотя бы одно из следующей системы равенств:
:<math>|r_i|_{k_{ij}}=|r_j|_{k_{ij}}, k_{ij} = gcd(m_i, m_j)</math>,
то произошла ошибка.

'''Теорема 2'''. Для системы из четырёх модулей <math>(m_1, m_2, m_3, m_4)</math> ошибка в одном разряде влечёт за собой невыполнением как минимум двух, а как максимум трёх равенств из системы выше. Причём, если ошибка в <math>r_i</math>, то равенства, которые не выполняются, содержат <math>r_i</math>.

==== Алгоритм локализации ошибки ====
Теорема 2 даёт алгоритм локализации ошибки. Ошибочным является разряд, который содержится во всех равенствах, которые не выполняются.
==== Алгоритм коррекции ошибки ====
Предложим, что ошибка в <math>r_i</math> и <math>\bar{r_i}</math> - правильное значение. Тогда <math>\bar{r_i}</math> можно найти, решая следующую систему уравнений:
:<math>|\bar{r_i}|_{k_a} = |r_j|_{k_a}, k_a = gcd(m_i, m_j)</math>,
:<math>|\bar{r_i}|_{k_b} = |r_k|_{k_b}, k_b = gcd(m_i, m_k)</math>.

Решение заключается в составлении множеств <math>A</math> и <math>B</math>:
:<math>A = \left \{ f = |r_j|_{k_a} + k_a k_x : f < m_i, k_x \geq 0 \right \}</math>,
:<math>B = \left \{ g = |r_k|_{k_b} + k_b k_y : g < m_i, k_y \geq 0 \right \}</math>.
и поиске пересечения <math>A \cap B</math>. Можно доказать, что такое пересечение всегда непусто и состоит ровно из одного элемента.

==== Примеры ====
'''Пример обнаружения одной ошибки.'''
Воспользуемся Построением 1 для получения системы модулей. Используя попарно взаимно простые циклические числа <math>(c_1, c_2, c_3, c_4)=(2, 3, 5, 7)</math>, получим систему модулей <math>(m_1, m_2, m_3, m_4)=(105, 70, 42, 30)</math>. Пусть теперь в остаточном представлении <math>(r_1, r_2, r_3, r_4)=(30, 60, 30, 0)</math> содержится одна ошибка. Проверяя совместность представления, найдём, что следующие равенства не выполняются:
:<math>|r_1|_{35} = |r_2|_{35}, 35 = gcd(105, 70)</math>,
:<math>|r_2|_{14} = |r_3k|_{14}, 14 = gcd(70, 42)</math>.

Так как <math>r_2</math> является общим в этих равенствах, ошибочным является остаток <math>60</math>. Для того чтобы получить правильное значение, составим множества <math>A</math> и <math>B</math> как было описано выше:
:<math>A = \left \{ f = |r_1|_{35} + 35 k_x : f < 70, k_x \geq 0 \right \}</math>,
:<math>B = \left \{ g = |r_3|_{14} + 14 k_y : g < 70, k_y \geq 0 \right \}</math>.

Получим:
:<math>A = \left \{ 30, 65 \right \}</math>,
:<math>B = \left \{ 2,16,30,44,58\right \}</math>.

Таким образом скорректированное значение <math>\{ \bar {r_2} \}= A \cap B = \{ 30 \}</math>.

Необязательно вычислять множество <math>B</math> полностью. Так как <math>|r_1|_{35} > 14</math>, то <math>k_y \geq 2</math>. Вычисления элементов множеств следует прекратить, как только их пересечение больше не пусто.
'''Пример обнаружения двух и коррекции одной ошибки.'''
Для обнаружения двух и коррекции одной ошибки требуется по крайней мере 5 модулей и кодовое расстояние равное 4. Если мы выберем циклические числа <math>(2, 3, 5, 7, 11)</math> мы получим набор модулей <math>(m_1, m_2, m_3, m_4, m_5) = (1155, 770, 462, 330, 210)</math>. Чтобы проверить, есть ли ошибки в остаточном представлении <math>(r_1, r_2, r_3, r_4, r_5)</math>, нужно выяснить, какие из следующих равенств не выполняются:
:<math>|r_1|_{385} = |r_2|_{385}, 385 = gcd(1155, 770)</math>,
:<math>|r_1|_{231} = |r_3|_{231}, 231 = gcd(1155, 462)</math>,
:<math>|r_1|_{105} = |r_5|_{105}, 105 = gcd(1155, 210)</math>,
:<math>|r_2|_{154} = |r_3|_{154}, 154 = gcd(462, 770)</math>,
:<math>|r_2|_{110} = |r_4|_{110}, 110 = gcd(330, 770)</math>,
:<math>|r_2|_{70} = |r_5|_{70}, 70 = gcd(210, 770)</math>,
:<math>|r_3|_{66} = |r_4|_{66}, 66 = gcd(462, 330)</math>,
:<math>|r_3|_{42} = |r_5|_{42}, 42 = gcd(462, 210)</math>,
:<math>|r_4|_{30} = |r_5|_{30}, 30 = gcd(330, 210)</math>.
Следует отметить, что если три или четыре из равенств выше не выполняются, то обнаружена одна ошибка. Если их пять и больше, то ошибок две. Аналогично предыдущему примеру, можно построить три множества <math>A, B, C</math>. Исправленное значение ошибочного разряда получим, найдя пересечение этих множеств:
<math>\{ \bar {r_i} \}= A \cap B \cap C</math>.

=== Модули с произвольными циклическими числами ===
Если мы уберём ограничение на взаимную простоту циклических числе, то обнаружение и коррекция ошибки может быть всё ещё осуществлена сокращением динамического диапазона. Следующие теорема и следствие связывают кодовое расстояние и динамический диапазон.

'''Теорема 3'''. Пусть есть множество модулей <math>(m_1, m_2, ..., m_n)</math> и <math>(l_1, l_2, ..., l_a)</math> - множество наименьших общих кратных всевозможных различных пар модулей, где <math>a = C^2_n</math>. Тогда, если <math>l = \min{(l_1, l_2, ..., l_a)}</math>, то кодовое расстояние множества всех представлений чисел <math>\{0, ..., l-1\}</math> в системе остаточных классов равно <math>n-1</math>.

'''Следствие'''. Пусть есть множество модулей <math>(m_1, m_2, ..., m_n)</math> и <math>(l_1, l_2, ..., l_a)</math> - множество наименьших общих кратных всевозможных кортежей из <math>i</math> различных модулей, где <math>a = C^i_n</math>. Тогда, если <math>l = \min{(l_1, l_2, ..., l_a)}</math>, то кодовое расстояние множества всех представлений чисел <math>\{0, ..., l-1\}</math> в системе остаточных классов равно <math>n - 1 + i</math>.

Далее будем предполагать, что ошибочный разряд только один. Из теоремы 3 следует, что для коррекции одной ошибки требуется множество из 4 модулей <math>(m_1, m_2, m_3, m_4)</math>. Пусть <math>l_1, l_2, l_3</math> - соответственно наименьшее из НОК всех пар, наименьшее из НОК всех троек и НОК всех модулей. Из теоремы 3 и следствия из неё следует, что для множества чисел <math>\{0, ..., l_1-1\}</math> кодовое расстояние соответствующего множества кодовых слов равно <math>3</math>, а для множеств <math>\{0, ..., l_2-1\}</math> и <math>\{0, ..., l_3-1\}</math> - не меньше <math>2</math> и не меньше <math>1</math> соответственно. Разделим множество всех кодовых слов на три подмножества:
:<math>A = \{x : 0 \leq x \leq l_1-1\}</math>,
:<math>B = \{x : l_1 \leq x \leq l_2-1\}</math>,
:<math>C = \{x : l_2 \leq x \leq l_3-1\}</math>.

Для коррекции одной ошибки необходимо кодовое расстояние, равное трём, то есть все слова из <math>A</math> и <math>B</math> выходят за границы динамического диапазона. Пусть необходимо обнаружить/исправить одну ошибку в остаточном представлении.

==== Алгоритм обнаружения ошибки ====
#Если <math>r_i \geqslant m_i</math>, то произошла ошибка в <math>i</math>-том разряде.
#Если остаточное представление не является совместным и слово не лежит в динамическом диапазоне, то произошла ошибка. Для проверки нужно воспользоваться КТО для перевода числа из остаточного представления и проверить, выполняется ли <math>x \leq l_1-1</math>.
#Если остаточное представление не является совместным, то произошла ошибка.
#Если остаточное представление является совместным и лежит в динамическом диапазоне, то ошибок нет.
==== Алгоритм коррекции ошибки ====
Начнём описание алгоритма со следующей леммы.

'''Лемма 1'''. Существует и единственно кодовое слово в множестве <math>A</math>, такое что кодовое расстояние между этим словом и остаточным представлением <math>(r_1, r_2, r_3, r_4)</math> равно <math>1</math>.

Опишем теперь коррекцию ошибки для каждого варианта:
#Кодовое слово за пределами динамического диапазона и совместно.
#Кодовое слово не совместно.

'''Вариант 1'''. Пусть остаточное представление <math>(r_1, r_2, r_3, r_4)</math> соответствует целому <math>R</math>. Так как только один разряд ошибочный, то <math>R</math> принадлежит множеству <math>C</math>. Коррекция выполняется следующим образом:
#Для всех кортежей из трёх остатков вычислить по КТО соответствующее целое.
#Если полученное целое принадлежит динамическому диапазону, то ошибка в исключённом разряде.
#Предположим без потери общности, что ошибочным является <math>r_3</math>. Найдём правильное значение <math>\bar {r_3}</math> из следующих уравнений:
:<math>|\bar{r_3}|_{k_a} = |r_1|_{k_a}, k_a = gcd(m_3, m_1)</math>,
:<math>|\bar{r_3}|_{k_b} = |r_2|_{k_b}, k_b = gcd(m_3, m_2)</math>,
:<math>|\bar{r_3}|_{k_c} = |r_4|_{k_c}, k_c = gcd(m_3, m_4)</math>.

Если полученных решений несколько, то согласно Лемме 1 только одно из них принадлежит динамическому диапазону.

'''Вариант 2'''. Пусть кодовое слово <math>(r_1, r_2, r_3, r_4)</math> не является совместным, то есть не выполняется равенство:
:<math>|r_i|_{k_a} = |r_j|_{k_a}, k_a = gcd(m_i, m_j)</math>.

Значит ошибочным является либо <math>r_i</math>, либо <math>r_j</math>. Пусть ошибка в <math>r_i</math>, найдём правильное значение, как в шаге 3 выше. Если для полученных решений, соответствующие целые не лежат в динамическом диапазоне, то ошибка в <math>r_j</math>, и её можно исправить тем же способом, описанным в шаге 3.

==== Примеры ====
Определим множество модулей <math>(m_1, m_2, m_3, m_4) = (8,6,4,2)</math>. Все кодовые слова и их целые эквиваленты приведены ниже:
:<math>(0,0,0,0) = 0</math>
:<math>(1,1,1,1) = 1</math>
:<math>(2,2,2,0) = 2</math>
:<math>(3,3,3,1) = 3</math>
:<math>(4,4,0,0) = 4</math>
:<math>(5,5,1,3) = 5</math>
:<math>(6,0,2,0) = 6</math>
:<math>(7,1,3,1) = 7</math>
:<math>(0,2,0,0) = 8</math>
:<math>(1,3,1,1) = 9</math>
:<math>(2,4,2,0) = 10</math>
:<math>(3,5,3,1) = 11</math>
:<math>(4,0,0,0) = 12</math>
:<math>(5,1,1,1) = 13</math>
:<math>(6,2,2,0) = 14</math>
:<math>(7,3,3,1) = 15</math>
:<math>(0,4,0,0) = 16</math>
:<math>(1,5,1,1) = 17</math>
:<math>(2,0,2,0) = 18</math>
:<math>(3,1,3,1) = 19</math>
:<math>(4,2,0,0) = 20</math>
:<math>(5,3,1,1) = 21</math>
:<math>(6,4,2,0) = 22</math>
:<math>(7,5,3,1) = 23</math>

Множества <math>A, B, C</math>, определённые выше, равны <math>A = \{0, ..., 3\}, B = \{4, ..., 7\}, C = \{8, ..., 23\}</math>.

Рассмотрим теперь несколько примеров обнаружения и коррекции.

'''Пример 1 '''. Остаточное представление <math>(r_1, r_2, r_3, r_4) = (2,0,2,0)</math> является совместным и не принадлежит динамическому диапазону. Четыре кортежа из трёх остатков и соответствующие им целые:
:<math>(2,0,2) = 18</math>
:<math>(0,2,0) = 6</math>
:<math>(2,0,0) = 18</math>
:<math>(2,2,0) = 2</math>

Единственное представление, лежащее в <math>\{0, ..., 3\}</math> - это <math>(2,2,0)</math>, следовательно <math>r_2 = 0</math> является ошибочным. Эта ошибка может быть исправлена, правильное значение можно найти из следующих уравнений:
:<math>|\bar {r_2}|_2 = |r_1|_2 = 0</math>,
:<math>|\bar {r_2}|_2 = |r_3|_2 = 0</math>,
:<math>|\bar {r_2}|_2 = |r_4|_2 = 0</math>,
:<math>\bar {r_2} \leq 5</math>.

Им удовлетворяют три решения: <math>0, 2, 4</math>. Правильным является <math> \bar {r_2} = 2</math>, так оставшиеся два соответствуют представлениям, лежащим вне динамического диапазона.

'''Пример 2 '''. Остаточное представление <math>(r_1, r_2, r_3, r_4) = (0,0,0,5)</math>. Так как <math>r_4 \geq 2</math>, то <math>r_4</math> - ошибочный. Найдём правильное значение из следующих уравнений:
:<math>|\bar {r_4}|_2 = |r_1|_2 = 0</math>,
:<math>|\bar {r_4}|_2 = |r_2|_2 = 0</math>,
:<math>|\bar {r_4}|_2 = |r_3|_2 = 0</math>,
:<math>\bar {r_4} \leq 1</math>.

Получаем <math>\bar {r_4} = 1</math>.

'''Пример 3 '''. Остаточное представление <math>(r_1, r_2, r_3, r_4) = (1,3,3,1)</math> является несовместным, так как не выполняется следующее равенство:
:<math>|r_1|_4 = |r_3|_4</math>.

Пусть <math>r_3</math> - ошибочный. Тогда скорректированное значение <math>r_3</math> может быть найдено из следующих уравнений:
:<math>|\bar {r_3}|_4 = |r_1|_4 = 1</math>,
:<math>|\bar {r_3}|_2 = |r_2|_2 = 1</math>,
:<math>|\bar {r_3}|_2 = |r_4|_2 = 1</math>,
:<math>\bar {r_3} \leq 3</math>.

Решая, находим <math> \bar {r_3} = 1 </math>. Скорректированное остаточное представление <math>(1,3,1,1) = 9</math> лежит вне динамического диапазона. Значит, ошибочным является <math>r_1</math>. Проводя аналогичную процедуру, находим <math>r_1 = 3</math>. Полученное остаточное представление <math>(3,3,3,1) = 3</math> лежит в динамическом диапазоне, коррекция завершена. Заметим, что мы пытались обнаружить ошибку в остаточном представлении <math>(1,3,3,1)</math>, которое имеет кодовое расстояние равное <math>1</math> с совместными представлениями <math>(1,3,1,1)</math> и <math>(3,3,3,1)</math>. По Лемме 1 только одно из них лежит в динамическом диапазоне.
== Литература ==
Rajendra S. Katti, A New Residue Arithmetic Error Correction Scheme, IEEE Trans. Computers, 1996, том 45, стр. 13-19

Схема коррекции ошибок Rajendra S. Katti

2013-07-31T12:34:08Z

Myachikov:

Схема коррекции ошибок Rajendra S. Katti

2013-07-31T10:50:13Z

Myachikov:

== Введение ==
Метод, описанный автором, основан на использовании избыточной системы остаточных классов. Недостатки представления чисел в такой системе с целью коррекции ошибок при вычислениях покажем на примере.

Пусть избыточная СОК определяется следующими модулями <math> m_1 = 3, m_2 = 4, m_3 = 5, m_4 = 7</math>, а динамический диапазон определяется <math>m_1</math> и <math>m_2</math> и соответственно равен <math>\{0,...,11\}</math>, а модули <math>m_3, m_4</math> являются дополнительными и необходимы для проверки на ошибки.

Натуральное число из динамического диапазона может быть восстановлено по любым трём из четырёх остатков по представленным модулям. Например, число <math>x=5=(2,1,0,2)</math> может быть представлено как <math>(2,1,0)</math> в системе модулей <math>\{m_1, m_2, m_3\}</math> или <math>(2,0,2), (1,0,2)</math> в системах <math>\{m_1, m_3, m_4\}</math> и <math>\{m_2, m_3, m_4\}</math> соответственно. Поэтому, если одна цифра отбрасывается в представлении <math>(2,1,0,2)</math>, правильный результат восстановления будет получен тогда и только тогда, когда отбрасываемая цифра ошибочна. Итак, обнаружение ошибок и коррекция выполняется следующим образом:
# Если число лежит вне динамического диапазона - произошла ошибка.
# Остатки отбрасываются по одному и если в результате одного из отбрасываний полученное число лежит в динамическом диапазоне, то отброшенный остаток является ошибочным.
# Корректное значение остатка получается расширением полученного на шаге 2 множества модулей, остатки по которым не содержат ошибок.

Такой подход требует много времени и вычислительных ресурсов. Ниже будет предложена новая схема обнаружения и коррекции ошибок на основе СОК.

== Модули, имеющие общие делители ==
Каждый набор остатков необязательно соответствует числу из динамического диапазона, если модули не являются попарно взаимно простыми. Легко показать, что такое соответствие есть тогда и только тогда, когда
:<math>|r_i|_k = |r_j|_k</math>,
выполняется для любых <math>i, j</math>, где <math>k = gcd(r_i, r_j)</math>. Таким образом, ошибочный разряд может быть найден выполнением только простой операции вычисления остатка.

Приведём пример, демонстрирующий, как можно проверить набор остатков на наличие ошибок. Проверим набор <math>(2, 6, 30, 42)</math> в системе модулей <math>\{4, 15, 36, 48\}</math>. Наибольшие общие делители всевозможных пар модулей соответственно равны: <math>gcd(4, 15) = 1, gcd(4, 36) = 4, gcd(4, 48) = 4, gcd(15, 36) = 3, gcd(15, 48) = 3, gcd(36, 48) = 12</math>.
Используя описанное выше соотношение для проверки, получаем: <math>|2|_4=|30|_4, |2|_4=|42|_4, |6|_3 = |30|_3, |6|_3=|42|_3</math> и <math> |30|_{12}=|42|_{12}</math>. То есть ошибок нет и число может быть восстановлено с использованием китайской теоремы об остатках (КТО) для модулей, имеющих нетривиальные общие делители.

'''КТО для системы модулей, попарно не взаимно простых'''. Пусть есть представление <math>(r_1, ..., r_n)</math> натурального числа <math>x</math>, определяемое системой модулей <math>(m_1, ..., m_n)</math>. Тогда
:<math>|x|_M = \left | \sum^{n}_{i=1} {{{\alpha}_i}r_i} \right |</math>,
где <math> {\alpha}_i </math>- число, такое что <math>|{\alpha}_i|_M = 0, |{\alpha}_i|_{\frac{M}{{\mu}_i}} = 1 </math>,
где <math> {\mu}_i </math> - числа, такие что <math> M = \prod^n_{i=1}{{\mu}_i} </math>, <math> {\mu}_i </math> делит <math> m_i </math> и <math> M = lcm(m_1, ..., m_n) </math>. Если для некоторого <math>i</math> такое <math>{\alpha}_i</math> не существует, то соответствующее слагаемое в сумме опускается.

Операции над числами остаются корректными и в системе не взаимно простых попарно модулей.

== Алгоритм ==
Если в множестве модулей <math>(m_1, ..., m_n)</math> есть имеющие нетривиальный общий делитель, то каждый набор остатков <math>(r_1, ..., r_n)</math> необязательно соответствует некоторому числу. Каждый набор, для которого такое соответствие имеет место, будем называть '''кодовым словом'''.

'''Расстоянием''' между двумя представлениями чисел в СОК(кодовыми словами) будем считать количество разрядов, в которых кодовые слова различаются между собой.

'''Кодовым расстоянием''' назовём минимальное расстояние среди всех пар кодовых слов.

Чтобы обнаружить и исправить <math>d-1</math> или меньше ошибок, необходимо и достаточно, чтобы кодовое расстояние было не меньше <math>d</math>. Не более <math>t</math> ошибок может быть обнаружено и исправлено тогда и только тогда, когда кодовое расстояние больше, чем <math>(2t+1)</math>. Если кодовое расстояние больше или равно, чем <math>(t+d+1)</math>, то любая комбинация из <math>t</math> ошибок может быть исправлена и до <math>d</math>ошибок могут быть обнаружены.
Для множества из <math>n</math> модулей в соответствующей системе могут быть представлены числа <math>(0, ..., M-1)</math>, где <math>M = lcm(m_1, ..., m_n)</math>.

'''Циклическим числом''' <math>c_i</math> для каждого модуля <math>m_i</math> назовём <math>c_i = M/m_i</math>.
=== Модули с попарно взаимно простыми циклическими числами ===
Построим теперь систему модулей с попарно взаимно простыми циклическими числами, которая поможет получит эффективный алгоритм обнаружения и коррекции ошибок.

'''Построение 1'''.
#Выберем <math>n</math> попарно взаимно простых циклических чисел, не равных единице <math>(c_1, ..., c_n)</math>.
#Получим систему модулей <math>(m_1, ..., m_n)</math> следующим образом:
:<math>m_i = \frac{1}{c_i}\prod^n_{k=1}{c_k}</math>.

Полученные модули обладают двумя свойствами:
#Наименьшие общие кратны всевозможных пар модулей равны между собой.
#Наибольшие общие делители всевозможных пар модулей различны.

'''Теорема 1'''. Для построенной системы модулей кодовое расстояние равно <math>n-1</math>, где <math>n</math> - количество модулей.

==== Алгоритм обнаружения ошибки ====
#Если <math>r_i \geqslant m_i</math>, то произошла ошибка в <math>i</math>-том разряде.
#Если не выполняется хотя бы одно из следующей системы равенств:
:<math>|r_i|_{k_{ij}}=|r_j|_{k_{ij}}, k_{ij} = gcd(m_i, m_j)</math>,
то произошла ошибка.

'''Теорема 2'''. Для системы из четырёх модулей <math>(m_1, m_2, m_3, m_4)</math> ошибка в одном разряде влечёт за собой невыполнением как минимум двух, а как максимум трёх равенств из системы выше. Причём, если ошибка в <math>r_i</math>, то равенства, которые не выполняются, содержат <math>r_i</math>.

==== Алгоритм локализации ошибки ====
Теорема 2 даёт алгоритм локализации ошибки. Ошибочным является разряд, который содержится во всех равенствах, которые не выполняются.
==== Алгоритм коррекции ошибки ====
Предложим, что ошибка в <math>r_i</math> и <math>\bar{r_i}</math> - правильное значение. Тогда <math>\bar{r_i}</math> можно найти, решая следующую систему уравнений:
:<math>|\bar{r_i}|_{k_a} = |r_j|_{k_a}, k_a = gcd(m_i, m_j)</math>,
:<math>|\bar{r_i}|_{k_b} = |r_k|_{k_b}, k_b = gcd(m_i, m_k)</math>.

Решение заключается в составлении множеств <math>A</math> и <math>B</math>:
:<math>A = \left \{ f = |r_j|_{k_a} + k_a k_x : f < m_i, k_x \geq 0 \right \}</math>,
:<math>B = \left \{ g = |r_k|_{k_b} + k_b k_y : g < m_i, k_y \geq 0 \right \}</math>.
и поиске пересечения <math>A \cap B</math>. Можно доказать, что такое пересечение всегда непусто и состоит ровно из одного элемента.

==== Примеры ====
'''Пример обнаружения одной ошибки.'''
Воспользуемся Построением 1 для получения системы модулей. Используя попарно взаимно простые циклические числа <math>(c_1, c_2, c_3, c_4)=(2, 3, 5, 7)</math>, получим систему модулей <math>(m_1, m_2, m_3, m_4)=(105, 70, 42, 30)</math>. Пусть теперь в остаточном представлении <math>(r_1, r_2, r_3, r_4)=(30, 60, 30, 0)</math> содержится одна ошибка. Проверяя совместность представления, найдём, что следующие равенства не выполняются:
:<math>|r_1|_{35} = |r_2|_{35}, 35 = gcd(105, 70)</math>,
:<math>|r_2|_{14} = |r_3k|_{14}, 14 = gcd(70, 42)</math>.

Так как <math>r_2</math> является общим в этих равенствах, ошибочным является остаток <math>60</math>. Для того чтобы получить правильное значение, составим множества <math>A</math> и <math>B</math> как было описано выше:
:<math>A = \left \{ f = |r_1|_{35} + 35 k_x : f < 70, k_x \geq 0 \right \}</math>,
:<math>B = \left \{ g = |r_3|_{14} + 14 k_y : g < 70, k_y \geq 0 \right \}</math>.

Получим:
:<math>A = \left \{ 30, 65 \right \}</math>,
:<math>B = \left \{ 2,16,30,44,58\right \}</math>.

Таким образом скорректированное значение <math>\{ \bar {r_2} \}= A \cap B = \{ 30 \}</math>.

Необязательно вычислять множество <math>B</math> полностью. Так как <math>|r_1|_{35} > 14</math>, то <math>k_y \geq 2</math>. Вычисления элементов множеств следует прекратить, как только их пересечение больше не пусто.
'''Пример обнаружения двух и коррекции одной ошибки.'''
Для обнаружения двух и коррекции одной ошибки требуется по крайней мере 5 модулей и кодовое расстояние равное 4. Если мы выберем циклические числа <math>(2, 3, 5, 7, 11)</math> мы получим набор модулей <math>(m_1, m_2, m_3, m_4, m_5) = (1155, 770, 462, 330, 210)</math>. Чтобы проверить, есть ли ошибки в остаточном представлении <math>(r_1, r_2, r_3, r_4, r_5)</math>, нужно выяснить, какие из следующих равенств не выполняются:
:<math>|r_1|_{385} = |r_2|_{385}, 385 = gcd(1155, 770)</math>,
:<math>|r_1|_{231} = |r_3|_{231}, 231 = gcd(1155, 462)</math>,
:<math>|r_1|_{105} = |r_5|_{105}, 105 = gcd(1155, 210)</math>,
:<math>|r_2|_{154} = |r_3|_{154}, 154 = gcd(462, 770)</math>,
:<math>|r_2|_{110} = |r_4|_{110}, 110 = gcd(330, 770)</math>,
:<math>|r_2|_{70} = |r_5|_{70}, 70 = gcd(210, 770)</math>,
:<math>|r_3|_{66} = |r_4|_{66}, 66 = gcd(462, 330)</math>,
:<math>|r_3|_{42} = |r_5|_{42}, 42 = gcd(462, 210)</math>,
:<math>|r_4|_{30} = |r_5|_{30}, 30 = gcd(330, 210)</math>.
Следует отметить, что если три или четыре из равенств выше не выполняются, то обнаружена одна ошибка. Если их пять и больше, то ошибок две. Аналогично предыдущему примеру, можно построить три множества <math>A, B, C</math>. Исправленное значение ошибочного разряда получим, найдя пересечение этих множеств:
<math>\{ \bar {r_i} \}= A \cap B \cap C</math>.

=== Модули с произвольными циклическими числами ===
Если мы уберём ограничение на взаимную простоту циклических числе, то обнаружение и коррекция ошибки может быть всё ещё осуществлена сокращением динамического диапазона. Следующие теорема и следствие связывают кодовое расстояние и динамический диапазон.

'''Теорема 3'''. Пусть есть множество модулей <math>(m_1, m_2, ..., m_n)</math> и <math>(l_1, l_2, ..., l_a)</math> - множество наименьших общих кратных всевозможных различных пар модулей, где <math>a = C^2_n</math>. Тогда, если <math>l = \min{(l_1, l_2, ..., l_a)}</math>, то кодовое расстояние множества всех представлений чисел <math>\{0, ..., l-1\}</math> в системе остаточных классов равно <math>n-1</math>.

'''Следствие'''. Пусть есть множество модулей <math>(m_1, m_2, ..., m_n)</math> и <math>(l_1, l_2, ..., l_a)</math> - множество наименьших общих кратных всевозможных кортежей из <math>i</math> различных модулей, где <math>a = C^i_n</math>. Тогда, если <math>l = \min{(l_1, l_2, ..., l_a)}</math>, то кодовое расстояние множества всех представлений чисел <math>\{0, ..., l-1\}</math> в системе остаточных классов равно <math>n - 1 + i</math>.

Далее будем предполагать, что ошибочный разряд только один. Из теоремы 3 следует, что для коррекции одной ошибки требуется множество из 4 модулей <math>(m_1, m_2, m_3, m_4)</math>. Пусть <math>l_1, l_2, l_3</math> - соответственно наименьшее из НОК всех пар, наименьшее из НОК всех троек и НОК всех модулей. Из теоремы 3 и следствия из неё следует, что для множества чисел <math>\{0, ..., l_1-1\}</math> кодовое расстояние соответствующего множества кодовых слов равно <math>3</math>, а для множеств <math>\{0, ..., l_2-1\}</math> и <math>\{0, ..., l_3-1\}</math> - не меньше <math>2</math> и не меньше <math>1</math> соответственно. Разделим множество всех кодовых слов на три подмножества:
:<math>A = \{x : 0 \leq x \leq l_1-1\}</math>,
:<math>B = \{x : l_1 \leq x \leq l_2-1\}</math>,
:<math>C = \{x : l_2 \leq x \leq l_3-1\}</math>.

Для коррекции одной ошибки необходимо кодовое расстояние, равное трём, то есть все слова из <math>A</math> и <math>B</math> выходят за границы динамического диапазона. Пусть необходимо обнаружить/исправить одну ошибку в остаточном представлении.

==== Алгоритм обнаружения ошибки ====
#Если <math>r_i \geqslant m_i</math>, то произошла ошибка в <math>i</math>-том разряде.
#Если остаточное представление не является совместным и слово не лежит в динамическом диапазоне, то произошла ошибка. Для проверки нужно воспользоваться КТО для перевода числа из остаточного представления и проверить, выполняется ли <math>x \leq l_1-1</math>.
#Если остаточное представление не является совместным, то произошла ошибка.
#Если остаточное представление является совместным и лежит в динамическом диапазоне, то ошибок нет.
==== Алгоритм коррекции ошибки ====
Начнём описание алгоритма со следующей леммы.

'''Лемма 1'''. Существует и единственно кодовое слово в множестве <math>A</math>

Схема коррекции ошибок Rajendra S. Katti

2013-07-22T13:12:04Z

Myachikov:

Схема коррекции ошибок Rajendra S. Katti

2013-07-22T12:55:22Z

Myachikov:

Схема коррекции ошибок Rajendra S. Katti

2013-07-22T09:41:05Z

Myachikov: Новая страница: «== Введение == Метод, описанный автором, основан на использовании избыточной системы оста…»

Фильтр с конечной импульсной характеристикой

2013-07-15T14:02:23Z

Myachikov:

'''Фильтр с конечной импульсной характеристикой''' ('''Нерекурсивный фильтр''', '''КИХ-фильтр''') или FIR-фильтр (FIR сокр. от finite impulse response — конечная импульсная характеристика) — один из видов линейных цифровых фильтров, характерной особенностью которого является ограниченность по времени его импульсной характеристики (с какого-то момента времени она становится точно равной нулю). Такой фильтр называют ещё нерекурсивным из-за отсутствия обратной связи. Знаменатель передаточной функции такого фильтра — некая константа.

== Динамические характеристики ==
Разностное уравнение, описывающее связь между входным и выходным сигналами фильтра:
<math>y\left (n\right)=b_0 x\left(n\right)+ b_1 x\left(n-1 \right)+...+b_P x\left(n-P \right)</math>
где <math>P</math> — порядок фильтра, <math>x(n)</math> — входной сигнал, <math>y(n)</math> — выходной сигнал, а <math>b_{i}</math> — коэффициенты фильтра. Иными словами, значение любого отсчета выходного сигнала определяется суммой масштабированных значений <math>P</math> предыдущих отсчетов. Можно сказать иначе: значение выхода фильтра в любой момент времени есть значение отклика на мгновенное значение входа и сумма всех постепенно затухающих откликов <math>P</math> предыдущих отсчетов сигнала, которые всё ещё оказывают влияние на выход (после <math>P</math>-отсчетов импульсная переходная функция становится равной нулю, как уже было сказано, поэтому все члены после <math>P</math>-го тоже станут равными нулю). Запишем предыдущее уравнение в более ёмком виде:
: <math>y \left( n \right) = \sum_{i=0}^{P} b_i x \left( n-i \right)</math>
Для того, чтобы найти ядро фильтра положим
: <math>x(n) = \delta(n)</math>
где <math>\delta(n)</math> — дельта-функция. Тогда импульсная характеристика КИХ-фильтра может быть записана как:
: <math>h\left (n\right)=\sum_{i=0}^{P}b_i \delta\left(n-i\right)</math>
Z-преобразование импульсной характеристики даёт нам передаточную функцию КИХ-фильтра:
: <math>H\left(z\right)=\sum_{i=0}^{P}b_i z^{-i}</math>

== Свойства ==
КИХ-фильтр обладает рядом полезных свойств, из-за которых он иногда более предпочтителен в использовании, чем БИХ-фильтр. Вот некоторые из них:
* КИХ-фильтры устойчивы.
* КИХ-фильтры при реализации не требуют наличия обратной связи.
* Фаза КИХ-фильтров может быть сделана линейной

== Прямая форма КИХ фильтра ==
КИХ фильтры могут быть реализованы с использованием трех элементов: умножитель, сумматор и блок задержки. Вариант, показанный на рисунке, есть прямая реализация КИХ-фильтров типа 1.
[[Файл:FIR Filter.png|center|300px|thumb|Реализация прямой формы КИХ фильтра]]

== Пример программы ==
Ниже приведен пример программы КИХ-фильтра, написанный на C :
<pre>

#define N 100 // порядок фильтра
float h[N] = { #include “f1.h” }; //вставка файла с известными коэффициентами фильтра
float x[N];
float y[N];

short my_FIR(short sample_data)
{
float result = 0;

for ( int i = N - 2 ; i >= 0 ; i-- )
{
x[i + 1] = x[i];
y[i + 1] = y[i];
}

x[0] = (float)sample_data;

for (int k = 0; k < N; k++)
{
result = result + x[k]*h[k];
}
y[0] = result;

return ((short)result);
}

</pre>

== Построение КИХ-фильтров ==
Построение КИХ-фильтров заключается в выборе коэффициентов фильтра таким образом, чтобы полученная система имела требуемые характеристики. Чаще всего фильтр строят по заданной амплитудно-частотной характеристике. Существуют разные методы построения КИХ-фильтров:
# Построение методом окна
# Метод частотной выборки
# Метод наименьших квадратов
# Метод Паркса-Макклиллана (также известный как Equiripple, Optimal, или Minimax метод). Обычно используется алгоритм Ремеза, чтобы найти оптимальный набор equiripple-коэффициентов. Задаётся желаемая амплитудно-частотная характеристика, функция весов ошибки отклика и порядок фильтра ''N''. Затем алгоритм находит множество из <math>\scriptstyle (N \,+\, 1)</math> коэффициентов, которые минимизируют максимум отклонения от идеального отклика. Интуитивно понятно, что этот алгоритм находит фильтр, который имеет АЧХ близкую к желаемой с той точностью, которая может быть достигнута с использованием <math>\scriptstyle (N \,+\, 1)</math> коэффициентов. Этот метод легко реализуется поскольку по крайней мере одна статья<ref>Rabiner, Lawrence R., and Gold, Bernard, 1975: Theory and Application of Digital Signal Processing (Englewood Cliffs, New Jersey: Prentice-Hall, Inc.) ISBN 0-13-914101-4</ref> включает в себя программу, которая принимает требуемую АЧХ и ''N'', и возвращает оптимальные коэффициенты фильтра.
# Фильтры Чебышева могут быть с тем же успехом построены с использованием алгоритмов БПФ. Алгоритм итеративен по своей природе. Вычисляется ДПФ от начальной характеристики фильтра с помощью БПФ (в качестве начального приближения можно взять h[n]=delta[n]). В частотной области проводится коррекция АЧХ в соответствии со спецификациями и вычисляется ОДПФ. Во временной области оставляют только N коэффициентов(обнулив оставшиеся). Далее повторяют БПФ, коррекцию и т.д.

Программные пакеты такие как MATLAB, GNU Octave, Scilab, и SciPy реализуют описанные выше методы.

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

=== Построение методом окна ===

В данном методе сначала строится идеальный БИХ-фильтр, затем к нему применяется Функция окна – во временной области, умножением бесконечной импульсной характеристики на функцию окна. Затем выполняется свёртка полученного результата с откликом на функцию окна. Если идеальная АЧХ достаточно проста, например, как прямоугольный отклик, результат свёртки можно сравнительно легко определить. На самом деле обычно сначала определяют требуемый результат и решают обратную задачу определения подходящей функции окна. Для этого метода особенно хорошо подходят окна Кайзера.

== Пример КИХ-фильтра ==
Примером простого КИХ-фильтра может служить скользящее среднее, используемое для обработки сигналов и изображений, системах автоматического управления и для других прикладных целей.
Разностное уравнение, которое характеризует фильтр скользящее среднее, является уравнением КИХ-фильтра. Пусть <math>x \left( n \right)</math> — входной сигнал фильтра, <math>y \left( n \right)</math> — выходной сигнал. Тогда разностное уравнение будет иметь вид:
: <math>y \left( n \right) = \sum_{i=0}^{P} b_i x \left( n-i \right)</math>

Отличительной особенностью скользящего среднего является равенство единице суммы коэффициентов <math>\ b_i </math>:
: <math> \sum_{i=0}^{P} b_i = 1</math>,

Последнее равенство отличает скользящее среднее от любого другого КИХ-фильтра. В частности для простого скользящего среднего:
: <math>b_{i}=\frac{1}{P+1} </math> для <math>i=0,1,\dots,P</math>

Для того, чтобы найти импульсную переходную функцию скользящего среднего делается следующее допущение:
: <math>\ x(n) = \delta(n)</math>

где <math>\ \delta(n)</math> — дельта-функция. Тогда импульсная характеристика такого фильтра может быть записана как:
: <math>h\left (n\right)=\sum_{i=0}^{P}b_i \delta\left(n-i\right)</math>

Z-преобразование импульсной характеристики даёт передаточную функцию:
: <math>H\left(z\right)=\sum_{i=0}^{P}b_i z^{-i}</math>

==Литература==

<references/>

== Ссылки ==
* [http://www.dsplib.ru/content/filters/fir/fir.html Расчет КИХ фильтра с линейной фазочастотной характеристикой методом частотной выборки]
* [http://www.dspguru.com/dsp/faqs/fir FIR FAQ]

Фильтр с конечной импульсной характеристикой

2013-07-15T13:59:09Z

Myachikov:

'''Фильтр с конечной импульсной характеристикой''' ('''Нерекурсивный фильтр''', '''КИХ-фильтр''') или FIR-фильтр (FIR сокр. от finite impulse response — конечная импульсная характеристика) — один из видов линейных цифровых фильтров, характерной особенностью которого является ограниченность по времени его импульсной характеристики (с какого-то момента времени она становится точно равной нулю). Такой фильтр называют ещё нерекурсивным из-за отсутствия обратной связи. Знаменатель передаточной функции такого фильтра — некая константа.

== Динамические характеристики ==
Разностное уравнение, описывающее связь между входным и выходным сигналами фильтра:
<math>y\left (n\right)=b_0 x\left(n\right)+ b_1 x\left(n-1 \right)+...+b_P x\left(n-P \right)</math>
где <math>P</math> — порядок фильтра, <math>x(n)</math> — входной сигнал, <math>y(n)</math> — выходной сигнал, а <math>b_{i}</math> — коэффициенты фильтра. Иными словами, значение любого отсчета выходного сигнала определяется суммой масштабированных значений <math>P</math> предыдущих отсчетов. Можно сказать иначе: значение выхода фильтра в любой момент времени есть значение отклика на мгновенное значение входа и сумма всех постепенно затухающих откликов <math>P</math> предыдущих отсчетов сигнала, которые всё ещё оказывают влияние на выход (после <math>P</math>-отсчетов импульсная переходная функция становится равной нулю, как уже было сказано, поэтому все члены после <math>P</math>-го тоже станут равными нулю). Запишем предыдущее уравнение в более ёмком виде:
: <math>y \left( n \right) = \sum_{i=0}^{P} b_i x \left( n-i \right)</math>
Для того, чтобы найти ядро фильтра положим
: <math>x(n) = \delta(n)</math>
где <math>\delta(n)</math> — дельта-функция. Тогда импульсная характеристика КИХ-фильтра может быть записана как:
: <math>h\left (n\right)=\sum_{i=0}^{P}b_i \delta\left(n-i\right)</math>
Z-преобразование импульсной характеристики даёт нам передаточную функцию КИХ-фильтра:
: <math>H\left(z\right)=\sum_{i=0}^{P}b_i z^{-i}</math>

== Свойства ==
КИХ-фильтр обладает рядом полезных свойств, из-за которых он иногда более предпочтителен в использовании, чем БИХ-фильтр. Вот некоторые из них:
* КИХ-фильтры устойчивы.
* КИХ-фильтры при реализации не требуют наличия обратной связи.
* Фаза КИХ-фильтров может быть сделана линейной

== Прямая форма КИХ фильтра ==
КИХ фильтры могут быть реализованы с использованием трех элементов: умножитель, сумматор и блок задержки. Вариант, показанный на рисунке, есть прямая реализация КИХ-фильтров типа 1.
[[Файл:FIR Filter.png|center|300px|thumb|Реализация прямой формы КИХ фильтра]]

== Пример программы ==
Ниже приведен пример программы КИХ-фильтра, написанный на C :
<pre>

#define N 100 // порядок фильтра
float h[N] = { #include “f1.h” }; //вставка файла с известными коэффициентами фильтра
float x[N];
float y[N];

short my_FIR(short sample_data)
{
float result = 0;

for ( int i = N - 2 ; i >= 0 ; i-- )
{
x[i + 1] = x[i];
y[i + 1] = y[i];
}

x[0] = (float)sample_data;

for (int k = 0; k < N; k++)
{
result = result + x[k]*h[k];
}
y[0] = result;

return ((short)result);
}

</pre>

== Построение КИХ-фильтров ==
Построение КИХ-фильтров заключается в выборе коэффициентов фильтра таким образом, чтобы полученная система имела требуемые характеристики. Чаще всего фильтр строят по заданной амплитудно-частотной характеристике. Существуют разные методы построения КИХ-фильтров:
# Построение методом окна
# Метод частотной выборки
# Метод наименьших квадратов
# Метод Паркса-Макклиллана (также известный как Equiripple, Optimal, или Minimax метод). Алгоритм Ремеза обычно используется, чтобы найти оптимальный набор equiripple-коэффициентов. Задаётся желаемая амплитудно-частотная характеристика, функция весов ошибки отклика и порядок фильтра ''N''. Затем алгоритм находит множество из <math>\scriptstyle (N \,+\, 1)</math> коэффициентов, которые минимизируют максимум отклонения от идеального отклика. Интуитивно понятно, что этот алгоритм находит фильтр, который имеет АЧХ близкую к желаемой с той точностью, которая может быть достигнута с использованием <math>\scriptstyle (N \,+\, 1)</math> коэффициентов. Этот метод легко реализуется поскольку по крайней мере одна статья<ref>Rabiner, Lawrence R., and Gold, Bernard, 1975: Theory and Application of Digital Signal Processing (Englewood Cliffs, New Jersey: Prentice-Hall, Inc.) ISBN 0-13-914101-4</ref> включает в себя программу, которая принимает требуемую АЧХ и ''N'', и возвращает оптимальные коэффициенты фильтра.
# Фильтры Чебышева могут быть с тем же успехом построены с использованием алгоритмов БПФ. Алгоритм итеративен по своей природе. Вычисляется ДПФ от начальной характеристики фильтра с помощью БПФ (в качестве начального приближения можно взять h[n]=delta[n]). В частотной области проводится коррекция АЧХ в соответствии со спецификациями и вычисляется ОДПФ. Во временной области оставляют только N коэффициентов(обнулив оставшиеся). Далее повторяют БПФ, коррекцию и т.д.

Программные пакеты такие как MATLAB, GNU Octave, Scilab, и SciPy реализуют описанные выше методы.

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

=== Построение методом окна ===

В данном методе сначала строится идеальный БИХ-фильтр, затем к нему применяется Функция окна – во временной области, умножением бесконечной импульсной характеристики на функцию окна. Затем выполняется свёртка полученного результата с откликом на функцию окна. Если идеальная АЧХ достаточно проста, например, как прямоугольный отклик, результат свёртки можно сравнительно легко определить. На самом деле обычно сначала определяют требуемый результат и решают обратную задачу определения подходящей функции окна. Для этого метода особенно хорошо подходят окна Кайзера.

== Пример КИХ-фильтра ==
Примером простого КИХ-фильтра может служить скользящее среднее, используемое для обработки сигналов и изображений, системах автоматического управления и для других прикладных целей.
Разностное уравнение, которое характеризует фильтр скользящее среднее, является уравнением КИХ-фильтра. Пусть <math>x \left( n \right)</math> — входной сигнал фильтра, <math>y \left( n \right)</math> — выходной сигнал. Тогда разностное уравнение будет иметь вид:
: <math>y \left( n \right) = \sum_{i=0}^{P} b_i x \left( n-i \right)</math>

Отличительной особенностью скользящего среднего является равенство единице суммы коэффициентов <math>\ b_i </math>:
: <math> \sum_{i=0}^{P} b_i = 1</math>,

Последнее равенство отличает скользящее среднее от любого другого КИХ-фильтра. В частности для простого скользящего среднего:
: <math>b_{i}=\frac{1}{P+1} </math> для <math>i=0,1,\dots,P</math>

Для того, чтобы найти импульсную переходную функцию скользящего среднего делается следующее допущение:
: <math>\ x(n) = \delta(n)</math>

где <math>\ \delta(n)</math> — [[дельта-функция]]. Тогда импульсная характеристика такого фильтра может быть записана как:
: <math>h\left (n\right)=\sum_{i=0}^{P}b_i \delta\left(n-i\right)</math>

Z-преобразование импульсной характеристики даёт передаточную функцию:
: <math>H\left(z\right)=\sum_{i=0}^{P}b_i z^{-i}</math>

==Литература==

<references/>

== Ссылки ==
* [http://www.dsplib.ru/content/filters/fir/fir.html Расчет КИХ фильтра с линейной фазочастотной характеристикой методом частотной выборки]
* [http://www.dspguru.com/dsp/faqs/fir FIR FAQ]

Фильтр с конечной импульсной характеристикой

2013-07-15T13:48:42Z

Myachikov:

'''Фильтр с конечной импульсной характеристикой''' ('''Нерекурсивный фильтр''', '''КИХ-фильтр''') или FIR-фильтр (FIR сокр. от finite impulse response — конечная импульсная характеристика) — один из видов линейных цифровых фильтров, характерной особенностью которого является ограниченность по времени его импульсной характеристики (с какого-то момента времени она становится точно равной нулю). Такой фильтр называют ещё нерекурсивным из-за отсутствия обратной связи. Знаменатель передаточной функции такого фильтра — некая константа.

== Динамические характеристики ==
Разностное уравнение, описывающее связь между входным и выходным сигналами фильтра:
<math>y\left (n\right)=b_0 x\left(n\right)+ b_1 x\left(n-1 \right)+...+b_P x\left(n-P \right)</math>
где <math>P</math> — порядок фильтра, <math>x(n)</math> — входной сигнал, <math>y(n)</math> — выходной сигнал, а <math>b_{i}</math> — коэффициенты фильтра. Иными словами, значение любого отсчета выходного сигнала определяется суммой масштабированных значений <math>P</math> предыдущих отсчетов. Можно сказать иначе: значение выхода фильтра в любой момент времени есть значение отклика на мгновенное значение входа и сумма всех постепенно затухающих откликов <math>P</math> предыдущих отсчетов сигнала, которые всё ещё оказывают влияние на выход (после <math>P</math>-отсчетов импульсная переходная функция становится равной нулю, как уже было сказано, поэтому все члены после <math>P</math>-го тоже станут равными нулю). Запишем предыдущее уравнение в более ёмком виде:
: <math>y \left( n \right) = \sum_{i=0}^{P} b_i x \left( n-i \right)</math>
Для того, чтобы найти ядро фильтра положим
: <math>x(n) = \delta(n)</math>
где <math>\delta(n)</math> — дельта-функция. Тогда импульсная характеристика КИХ-фильтра может быть записана как:
: <math>h\left (n\right)=\sum_{i=0}^{P}b_i \delta\left(n-i\right)</math>
Z-преобразование импульсной характеристики даёт нам передаточную функцию КИХ-фильтра:
: <math>H\left(z\right)=\sum_{i=0}^{P}b_i z^{-i}</math>

== Свойства ==
КИХ-фильтр обладает рядом полезных свойств, из-за которых он иногда более предпочтителен в использовании, чем БИХ-фильтр. Вот некоторые из них:
* КИХ-фильтры устойчивы.
* КИХ-фильтры при реализации не требуют наличия обратной связи.
* Фаза КИХ-фильтров может быть сделана линейной

== Прямая форма КИХ фильтра ==
КИХ фильтры могут быть реализованы с использованием трех элементов: умножитель, сумматор и блок задержки. Вариант, показанный на рисунке, есть прямая реализация КИХ-фильтров типа 1.
[[Файл:FIR Filter.png|center|300px|thumb|Реализация прямой формы КИХ фильтра]]

== Пример программы ==
Ниже приведен пример программы КИХ-фильтра, написанный на C :
<pre>

#define N 100 // порядок фильтра
float h[N] = { #include “f1.h” }; //вставка файла с известными коэффициентами фильтра
float x[N];
float y[N];

short my_FIR(short sample_data)
{
float result = 0;

for ( int i = N - 2 ; i >= 0 ; i-- )
{
x[i + 1] = x[i];
y[i + 1] = y[i];
}

x[0] = (float)sample_data;

for (int k = 0; k < N; k++)
{
result = result + x[k]*h[k];
}
y[0] = result;

return ((short)result);
}

</pre>

== Построение КИХ-фильтров ==
Построение КИХ-фильтров заключается в выборе коэффициентов фильтра таким образом, чтобы полученная система имела требуемые характеристики. Чаще всего фильтр строят по заданной амплитудно-частотной характеристике. Существуют разные методы построения КИХ-фильтров:
# Построение методом окна
# Метод частотной выборки
# Метод наименьших квадратов
# Метод Паркса-Макклиллана (также известный как Equiripple, Optimal, или Minimax метод). Алгоритм Ремеза обычно используется, чтобы найти оптимальный набор equiripple-коэффициентов. Задаётся желаемая амплитудно-частотная характеристика, функция весов ошибки отклика и порядок фильтра ''N''. Затем алгоритм находит множество из <math>\scriptstyle (N \,+\, 1)</math> коэффициентов, которые минимизируют максимум отклонения от идеального отклика. Интуитивно понятно, что этот алгоритм находит фильтр, который имеет АЧХ близкую к желаемой с той точностью, которая может быть достигнута с использованием <math>\scriptstyle (N \,+\, 1)</math> коэффициентов. Этот метод легко реализуется поскольку по крайней мере одна статья<ref>Rabiner, Lawrence R., and Gold, Bernard, 1975: Theory and Application of Digital Signal Processing (Englewood Cliffs, New Jersey: Prentice-Hall, Inc.) ISBN 0-13-914101-4</ref> включает в себя программу, которая принимает требуемую АЧХ и ''N'', и возвращает оптимальные коэффициенты фильтра.
# Фильтры Чебышева могут быть с тем же успехом построены с использованием алгоритмов БПФ. Алгоритм итеративен по своей природе. Вычисляется ДПФ от начальной характеристики фильтра с помощью БПФ (в качестве начального приближения можно взять h[n]=delta[n]). В частотной области проводится коррекция АЧХ в соответствии со спецификациями и вычисляется ОДПФ. Во временной области оставляют только N коэффициентов(обнулив оставшиеся). Далее повторяют БПФ, коррекцию и т.д.

Программные пакеты такие как MATLAB, GNU Octave, Scilab, и SciPy реализуют описанные выше методы.

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

=== Построение методом окна ===

В данном методе сначала строится идеальный БИХ-фильтр, затем к нему применяется Функция окна – во временной области, умножением бесконечной импульсной характеристики на функцию окна. Затем выполняется свёртка полученного результата с откликом на функцию окна. Если идеальная АЧХ достаточно проста, например, как прямоугольный отклик, результат свёртки можно сравнительно легко определить. На самом деле обычно сначала определяют требуемый результат и решают обратную задачу определения подходящей функции окна. Для этого метода особенно хорошо подходят окна Кайзера.

==Литература==

<references/>

== Ссылки ==
* [http://www.dsplib.ru/content/filters/fir/fir.html Расчет КИХ фильтра с линейной фазочастотной характеристикой методом частотной выборки]
* [http://www.dspguru.com/dsp/faqs/fir FIR FAQ]

Фильтр с конечной импульсной характеристикой

2013-07-15T13:47:01Z

Myachikov:

'''Фильтр с конечной импульсной характеристикой''' ('''Нерекурсивный фильтр''', '''КИХ-фильтр''') или FIR-фильтр (FIR сокр. от finite impulse response — конечная импульсная характеристика) — один из видов линейных цифровых фильтров, характерной особенностью которого является ограниченность по времени его импульсной характеристики (с какого-то момента времени она становится точно равной нулю). Такой фильтр называют ещё нерекурсивным из-за отсутствия обратной связи. Знаменатель передаточной функции такого фильтра — некая константа.

== Динамические характеристики ==
Разностное уравнение, описывающее связь между входным и выходным сигналами фильтра:
<math>y\left (n\right)=b_0 x\left(n\right)+ b_1 x\left(n-1 \right)+...+b_P x\left(n-P \right)</math>
где <math>P</math> — порядок фильтра, <math>x(n)</math> — входной сигнал, <math>y(n)</math> — выходной сигнал, а <math>b_{i}</math> — коэффициенты фильтра. Иными словами, значение любого отсчета выходного сигнала определяется суммой масштабированных значений <math>P</math> предыдущих отсчетов. Можно сказать иначе: значение выхода фильтра в любой момент времени есть значение отклика на мгновенное значение входа и сумма всех постепенно затухающих откликов <math>P</math> предыдущих отсчетов сигнала, которые всё ещё оказывают влияние на выход (после <math>P</math>-отсчетов импульсная переходная функция становится равной нулю, как уже было сказано, поэтому все члены после <math>P</math>-го тоже станут равными нулю). Запишем предыдущее уравнение в более ёмком виде:
: <math>y \left( n \right) = \sum_{i=0}^{P} b_i x \left( n-i \right)</math>
Для того, чтобы найти ядро фильтра положим
: <math>x(n) = \delta(n)</math>
где <math>\delta(n)</math> — дельта-функция. Тогда импульсная характеристика КИХ-фильтра может быть записана как:
: <math>h\left (n\right)=\sum_{i=0}^{P}b_i \delta\left(n-i\right)</math>
Z-преобразование импульсной характеристики даёт нам передаточную функцию КИХ-фильтра:
: <math>H\left(z\right)=\sum_{i=0}^{P}b_i z^{-i}</math>

== Свойства ==
КИХ-фильтр обладает рядом полезных свойств, из-за которых он иногда более предпочтителен в использовании, чем БИХ-фильтр. Вот некоторые из них:
* КИХ-фильтры устойчивы.
* КИХ-фильтры при реализации не требуют наличия обратной связи.
* Фаза КИХ-фильтров может быть сделана линейной

== Прямая форма КИХ фильтра ==
КИХ фильтры могут быть реализованы с использованием трех элементов: умножитель, сумматор и блок задержки. Вариант, показанный на рисунке, есть прямая реализация КИХ-фильтров типа 1.
[[Файл:FIR Filter.png|center|300px|thumb|Реализация прямой формы КИХ фильтра]]

== Пример программы ==
Ниже приведен пример программы КИХ-фильтра, написанный на C :
<pre>

#define N 100 // порядок фильтра
float h[N] = { #include “f1.h” }; //вставка файла с известными коэффициентами фильтра
float x[N];
float y[N];

short my_FIR(short sample_data)
{
float result = 0;

for ( int i = N - 2 ; i >= 0 ; i-- )
{
x[i + 1] = x[i];
y[i + 1] = y[i];
}

x[0] = (float)sample_data;

for (int k = 0; k < N; k++)
{
result = result + x[k]*h[k];
}
y[0] = result;

return ((short)result);
}

</pre>

== Построение КИХ-фильтров ==
Построение КИХ-фильтров заключается в выборе коэффициентов фильтра таким образом, чтобы полученная система имела требуемые характеристики. Чаще всего фильтр строят по заданной амплитудно-частотной характеристике. Существуют разные методы построения КИХ-фильтров:
# Window design method
# Метод частотной выборки
# Метод наименьших квадратов
# Метод Паркса-Макклиллана (также известный как Equiripple, Optimal, или Minimax метод). Алгоритм Ремеза обычно используется, чтобы найти оптимальный набор equiripple-коэффициентов. Задаётся желаемая амплитудно-частотная характеристика, функция весов ошибки отклика и порядок фильтра ''N''. Затем алгоритм находит множество из <math>\scriptstyle (N \,+\, 1)</math> коэффициентов, которые минимизируют максимум отклонения от идеального отклика. Интуитивно понятно, что этот алгоритм находит фильтр, который имеет АЧХ близкую к желаемой с той точностью, которая может быть достигнута с использованием <math>\scriptstyle (N \,+\, 1)</math> коэффициентов. Этот метод легко реализуется поскольку по крайней мере одна статья<ref>Rabiner, Lawrence R., and Gold, Bernard, 1975: Theory and Application of Digital Signal Processing (Englewood Cliffs, New Jersey: Prentice-Hall, Inc.) ISBN 0-13-914101-4</ref> включает в себя программу, которая принимает требуемую АЧХ и ''N'', и возвращает оптимальные коэффициенты фильтра.
# Фильтры Чебышева могут быть с тем же успехом построены с использованием алгоритмов БПФ. Алгоритм итеративен по своей природе. Вычисляется ДПФ от начальной характеристики фильтра с помощью БПФ (в качестве начального приближения можно взять h[n]=delta[n]). В частотной области проводится коррекция АЧХ в соответствии со спецификациями и вычисляется ОДПФ. Во временной области оставляют только N коэффициентов(обнулив оставшиеся). Далее повторяют БПФ, коррекцию и т.д.

Программные пакеты такие как MATLAB, GNU Octave, Scilab, и SciPy реализуют описанные выше методы.

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

=== Window design method ===

В методе Window Design Method сначала строится идеальный БИХ-фильтр, затем к нему применяется Функция окна – во временной области, умножением бесконечной импульсной характеристики на функцию окна. Затем выполняется свёртка полученного результата с откликом на функцию окна. Если идеальная АЧХ достаточно проста, например, как прямоугольный отклик, результат свёртки можно сравнительно легко определить. На самом деле обычно сначала определяют требуемый результат и решают обратную задачу определения подходящей функции окна. Для этого метода особенно хорошо подходят окна Кайзера.

==References==

<references/>

== Ссылки ==
* [http://www.dsplib.ru/content/filters/fir/fir.html Расчет КИХ фильтра с линейной фазочастотной характеристикой методом частотной выборки]
* [http://www.dspguru.com/dsp/faqs/fir FIR FAQ]

Фильтр с конечной импульсной характеристикой

2013-07-15T12:44:26Z

Myachikov:

Фильтр с конечной импульсной характеристикой

2013-07-15T12:34:43Z

Myachikov: Новая страница: «'''Фильтр с конечной импульсной характеристикой''' ('''Нерекурсивный фильтр''', '''КИХ-фильтр'…»

Файл:FIR Filter.png

2013-07-15T12:31:41Z

Myachikov: Myachikov загружена новая версия «Файл:FIR Filter.png»: Версия с более низким разрешением

Реализация прямой формы КИХ фильтра

Файл:FIR Filter.png

2013-07-15T12:30:02Z

Myachikov: Реализация прямой формы КИХ фильтра

Реализация прямой формы КИХ фильтра