CC BY
7Математика и математическое моделирование. 2020. № 04. С. 1 - 51
DPI: 10.24108/mathm.0420.0000223
шя
© И.В. Матюшкин, П.Д. Рубис, 2020 УДК 519.7
Четыре клеточно-автоматных алгоритма пермутаций матриц
Матюшкин И.В.1'2'*' Рубис П.Д.2
Научно-исследовательский институт молекулярной электроники, Москва, Зеленоград, Россия 2Национальный исследовательский университет «Московский институт электронной техники», Москва, Зеленоград , Россия
imatvuEhkin@niimeju
С помощью численного расчета описывается работа алгоритмов пермутаций матриц, основанных на циклических сдвигах строк и столбцов. Такой выбор алгоритмов дискретного преобразования обоснован удобством клеточно-автоматных формулировок, которые и приводятся. Получены эмпирические формулы для периода пермутаций; для последнего алгоритма формула периода носит рекуррентный характер. Для базовой и наиболее простой схемы период N (п) имеет асимптотику ехр (2п) / п для матрицы п х п с попарно различными элементами. Несмотря на усложнение схемы алгоритма, две другие модификации дают лишь полиномиальный рост степени не выше 3. Четвертая схема имеет нетривиальную зависимость периода, но не выше экспоненциальной. В ряде случаев алгоритмы порождают особые пермутации: поворот, отражение и перестановку блоков для матрицы 2к х 2к . Эти формулы тесно связаны с индивидуальными траекториями элементов, а они - с влиянием границ, что дает ветвление порядка матрицы по классу вычета по модулю 3,4 или 12. Визуализации этих траекторий приводятся в расширенном поле КА. В качестве параметра динамики КА анализируются две «метрики пере-мешанности» на пермутациях матрицы (по сравнению с начальной). Для всех схем и большинства ветвей поведение этих метрик представлено на графиках и гистограммах (условно: плотности распределения), показывающих, как часто встречаются по периоду пермутации с заданным интервалом значений метрик. Практическое значение работы состоит в оценке применения КА в областях генерации псевдослучайных чисел и криптографии.
Ключевые слова: клеточные автоматы, пермутация, случайные числа, криптография, метрика
Представлена в редакцию: 10.07.2020, исправлена 24.07.2020
* Работа выполнена при финансовой поддержке РФФИ, проект №20-07-00409
Введение
Пермутация, или перестановка элементов, матрицы — это один из способов генерации псевдослучайного числа. Простейшими пермутациями являются поворот, отражение и транспонирование матрицы. Перестановочные алгоритмы легко проектируемы на кле-точно-автоматном формализме, где уже относительно давно исследуются способы созда-
Математика ft Математическс
моделирование
Сетевое научное издание http://mathmelpub.ru ISSN 2412-591
ния псевдослучайных последовательностей, начиная с работ Стивена Вольфрама в ранних 1980-х гг. [1,2]. Клеточные автоматы (КА) обладают высокой скоростью вычислений при реализации в параллельной архитектуре [2], что обуславливает большое число исследований [3-8] с акцентом на аппаратную реализацию. Например, в [5] для эволюционирующих КА, ячейки которых связаны функцией XOR и XNOR, проверялась псевдослучайность числа в пакетном тесте DIEHARD.
Актуальность поиска эффективных способов псевдослучайной пермутации матрицы обусловлена применением этой операции в различных областях. Генерация тестовых последовательностей в логических блоках современных сверхбольших интегральных схем может требовать случайной перестановки данных [9-11], например, схемы c программируемой логикой. Помимо генерации случайных чисел перестановка матричных элементов непосредственно используется в криптографии [12-16]. Обычно вводится псевдослучайная функция, где первая часть декартова произведения относится к k-битному ключу, а вторая - к n-битному сообщению. При составлении стойких криптографических алгоритмов по сей день руководствуются утверждением Шеннона о том, что хорошее перемешивающее преобразование часто является результатом многократного применения произведения двух простых некоммутирующих операций. Такими простейшими операциями являются сдвиги строк и столбцов матрицы [17-18]. В [18] решение близко к предлагаемому нами в [ 19-20] и является дискретизированным аналогом непрерывного двухпараметриче-ского преобразования; оно применялось для кодирования серого 8-битного рисунка. Автор [17] дал более строгую формализацию и точную оценку периода пермутации для алгоритма [19], основанную на функции Ландау в теории групп.
1. КА-формализация схемы «1:2»
Рассмотрим сначала идею модифицируемого алгоритма. Квадратная матрица n х n помещена в поле КА с замыканием границ (тороидальная решётка). Алгоритм состоит из 4 периодически повторяющихся этапов, каждый из которых представляет циклический сдвиг, т.е. первый элемент ставится на место последнего, и наоборот (рис. 1):
1. Нечётные столбцы сдвигаются вверх, чётные — вниз.
2. Нечётные строки сдвигаются вправо, чётные — влево.
3. Нечётные столбцы сдвигаются вниз, чётные — вверх.
4. Нечётные строки сдвигаются влево, чётные — вправо.
Преобразование матрицы за один этап здесь и далее будем называть шагом или тактом, а за все этапы — проходом.
Рис. 1. Визуализация схемы «1:1»
Проход схемы «1:2» состоит из 6 шагов, за счет повторения этапов 2 и 4 схемы 1:1 на следующем после них шаге алгоритма.
Рис. 2. Визуализация схемы «1:2»
Приведём КА-формулировку схемы 1:2. Как и в схеме 1:1, используется квадратная решётка с замкнутыми границами (нумерация, традиционная для индексов матрицы). Окрестность клетки типа фон Неймана радиуса единица. Состояние клетки определяется также двумя компонентами (/, ^ (табл. 1).
Таблица 1. Описание компонент клетки (/, ^ схемы «1:2»
Название Флаг перемещения Регистр данных
Обозначение /
Множество значений / е {Т, 1, Ц = {0,1,2,3,5,7} Произвольное
Описание Отвечает за перемещение элементов матрицы. Стрелки в основном указывают направление движения. Хранит элементы матрицы M= |11.
Множество значений флага перемещения ячейки f дополняется двумя новыми элементами - 5 и 7. Их обоснование приведем после формулировки локальной функции перехода (ЛФП). Введем, как и ранее в схеме 1:1, понятие согласованности. Две клетки называются согласованными, если чётности их флагов совпадают. Число согласованности cons клетки определяется как число согласованных с ней соседей: cons = |{|| f = f mod2}|, где |
— локальный индекс клетки. Его значение, как и f, можно указать стрелкой. Например, « t» — сосед сверху, « ^ » — сосед справа, а f — значение флага состояния клетки снизу. Сама центральная клетка не индексируется. Знак «:=» имеет смысл присвоения значения, то есть операнд справа от него берётся до перехода, а слева — после перехода. Операция mod возвращает остаток от целочисленного деления (например, 9mod4 = 1, - lmod4 = 3 ).
Формальные правила перехода, записанные в виде последовательно выполняемых инструкций, приведены в табл. 2. Знаки л, V обозначают логические операции конъюнкции и дизъюнкции соответственно.
Таблица 2. ЛФП схемы «1:2»
№ Условие перехода Формула перехода
1 (cons = 0) Л (f = 0) v {{cons Ф 0) Л ((f = 3) v (f = 7))) s := ^
((cons = 0) л (f = 1)) v ((cons Ф 0) л (f = 0)) s :=
((cons = 0) л (f = 2)) v ((cons Ф 0) л ((f = 1) v (f = 5))) s :=
((cons = 0) л (f = 3)) v ((cons Ф 0) л (f = 2)) s := sT
2 (cons = 0) v (cons = 4) f :=( f +1) mod 4
—. ((cons = 0) v (cons = 4)) f :=( f -1) mod 4
(f = 1)Л( f_ = 1)Л( Д = 1) f := 5
(f = 3 )л( f_ = 3 )Л( Д = 3) f := 7
Без ограничения общности и в предположении попарного различия элементов положим начальное состояние слоя данных упорядоченным по столбцам, что соответствует соглашению МЛТЬЛВ, и рассмотрим эволюцию КА размера 8 х8 на первом проходе поэтапно. на всех шести этапах алгоритма (табл.3):
Таблица 3. Проход схемы «1:2» для матрицы 8 х 8
Номер итерации Слой данных Слой флагов
1 9 17 25 33 41 49 57 0 2 0 2 0 2 0 2
2 10 18 26 34 42 50 58 0 3 0 3 0 3 0 3
3 11 19 27 35 43 51 59 0 2 0 2 0 2 0 2
0 (исходная матрица) 4 12 20 28 36 44 52 60 0 3 0 3 0 3 0 3
5 13 21 29 37 45 53 61 0 2 0 2 0 2 0 2
6 14 22 30 38 46 54 62 0 3 0 3 0 3 0 3
7 15 23 31 39 47 55 63 0 2 0 2 0 2 0 2
8 16 24 32 40 48 56 64 0 3 0 3 0 3 0 3
2 16 18 32 34 48 50 64 1 1 1 1 1 1 1 1
3 9 19 25 35 41 51 57 3 0 3 0 3 0 3 0
4 10 20 26 36 42 52 58 1 1 1 1 1 1 1 1
1 5 11 21 27 37 43 53 59 3 0 3 0 3 0 3 0
6 12 22 28 38 44 54 60 1 1 1 1 1 1 1 1
7 13 23 29 39 45 55 61 3 0 3 0 3 0 3 0
8 14 24 30 40 46 56 62 1 1 1 1 1 1 1 1
1 15 17 31 33 47 49 63 3 0 3 0 3 0 3 0
2 64 2 16 18 32 34 48 50 5 5 5 5 5 5 5 5
9 19 25 35 41 51 57 3 3 0 3 0 3 0 3 0
58 4 10 20 26 36 42 52 5 5 5 5 5 5 5 5
11 21 27 37 43 53 59 5 3 0 3 0 3 0 3 0
60 6 12 22 28 38 44 54 5 5 5 5 5 5 5 5
13 23 29 39 45 55 61 7 3 0 3 0 3 0 3 0
62 8 14 24 30 40 46 56 5 5 5 5 5 5 5 5
15 17 31 33 47 49 63 1 3 0 3 0 3 0 3 0
3 50 64 2 16 18 32 34 48 2 0 2 0 2 0 2 0
19 25 35 41 51 57 3 9 2 1 2 1 2 1 2 1
52 58 4 10 20 26 36 42 2 0 2 0 2 0 2 0
21 27 37 43 53 59 5 11 2 1 2 1 2 1 2 1
54 60 6 12 22 28 38 44 2 0 2 0 2 0 2 0
23 29 39 45 55 61 7 13 2 1 2 1 2 1 2 1
56 62 8 14 24 30 40 46 2 0 2 0 2 0 2 0
17 31 33 47 49 63 1 15 2 1 2 1 2 1 2 1
4 17 25 33 41 49 57 1 9 3 3 3 3 3 3 3 3
50 58 2 10 18 26 34 42 1 2 1 2 1 2 1 2
19 27 35 43 51 59 3 11 3 3 3 3 3 3 3 3
52 60 4 12 20 28 36 44 1 2 1 2 1 2 1 2
21 29 37 45 53 61 5 13 3 3 3 3 3 3 3 3
54 62 6 14 22 30 38 46 1 2 1 2 1 2 1 2
23 31 39 47 55 63 7 15 3 3 3 3 3 3 3 3
56 64 8 16 24 32 40 48 1 2 1 2 1 2 1 2
5 25 33 41 49 57 1 9 17 7 7 7 7 7 7 7 7
42 50 58 2 10 18 26 34 1 2 1 2 1 2 1 2
27 35 43 51 59 3 11 19 7 7 7 7 7 7 7 7
44 52 60 4 12 20 28 36 1 2 1 2 1 2 1 2
29 37 45 53 61 5 13 21 7 7 7 7 7 7 7 7
46 54 62 6 14 22 30 38 1 2 1 2 1 2 1 2
31 39 47 55 63 7 15 23 7 7 7 7 7 7 7 7
48 56 64 8 16 24 32 40 1 2 1 2 1 2 1 2
6 33 41 49 57 1 9 17 25 0 2 0 2 0 2 0 2
34 42 50 58 2 10 18 26 0 3 0 3 0 3 0 3
35 43 51 59 3 11 19 27 0 2 0 2 0 2 0 2
36 44 52 60 4 12 20 28 0 3 0 3 0 3 0 3
37 45 53 61 5 13 21 29 0 2 0 2 0 2 0 2
38 46 54 62 6 14 22 30 0 3 0 3 0 3 0 3
39 47 55 63 7 15 23 31 0 2 0 2 0 2 0 2
40 48 56 64 8 16 24 32 0 3 0 3 0 3 0 3
Поясним добавление двух элементов 5,7 во множество флагов. Легко видеть, что после первого сдвига строк (шаг 2) матрица компонент / удовлетворяет нововведенным ус-
ловиям из второго блока ЛФП. Для всех элементов нечетных строк выполняется
(f = 1 )л(= 1 )л(f^ = 1). Соответственно их флагам состояния присваивается значение 5.
Учитывая условие ((cons = 0)л( f = 2))v((cons ф 0)л(( f = 1)v(f = 5))), а также равенство
1mod4 = 5mod4, легко видеть, что состояние элемента при f = 1 и при f = 5 изменяется одинаково. Все элементы четных срок удовлетворяют последнему условию второго блока ЛФП и не меняют значение f. Это также необходимо для повторения сдвига строк матрицы. Аналогично обосновываются другие дополнения ЛФП, и, соответственно, добавление цифры 7 в множество значений f .
3. КА-формализация схемы «1:0:1»
Цикличные сдвиги элементов здесь происходят с пропуском столбца или строки. Для лучшего понимания визуализируем для КА размера 5 х 5 слой данных (табл. 4)
Таблица 4. Проход схемы «1:0:1» для матрицы 5 х 5
Номер итерации Слой данных
0 (исходная матрица)
1 6 11 16 21
2 7 12 17 22
3 8 13 18 23
4 9 14 19 24
5 10 15 20 25
1
2 6 15 17 21
3 7 11 18 22
4 8 12 19 23
5 9 13 20 24
1 10 14 16 25
2
21 2 6 15 17
3 7 11 18 22
8 12 19 23 4
24 5 9 13 20
1 10 14 16 25
3
21 10 11 15 25
3 2 19 18 17
8 7 9 23 22
24 12 14 13 4
1 5 6 16 20
Номер итерации Слой данных
4
21 10 11 15 25
2 19 18 17 3
22 8 7 9 23
24 12 14 13 4
5 6 16 20 1
5
5 19 11 20 3
21 8 18 15 23
2 12 7 17 4
22 6 14 9 1
24 10 16 13 25
6
19 11 20 3 5
23 21 8 18 15
2 12 7 17 4
6 14 9 1 22
25 24 10 16 13
Состояния ячеек определяются четырьмя компонентами (V, g, Ь, я), представленными в табл. 5:
Таблица 5. Описание компонент клетки (V, g,Ь,я) схемы «1:0:1»
Название Флаг вертикального перемещения Флаг горизонтального перемещения Логический Флаг Регистр данных
Обозначение V g Ь я
Множество значений V е{Т,Ъ, ¿} = {1,0,-1} g ^} = {1,0, Ь е {0,1} Произвольное
Описание Отвечает за вертикальное перемещение элементов матрицы. Стрелки направление движения. Отвечает за горизонтальное перемещение элементов матрицы. Стрелки указывают направление движения. Обеспечивает чередование вертикальных и горизонтальных перемещений. Хранит элементы матрицы м=| т 11.
Флаги V и g определяют направления сдвига столбцов и строк соответственно. Значение 0 обоих флагов означает отсутствие смещения элемента на текущей итерации. Флаг Ь - бит, определяющий блок ЛФП на данном шаге алгоритма. Значение Ь , равное единице, предполагает перемещение в соответствии с ЛФП для компоненты V. Аналогично и для g при Ь = 0. Начальные конфигурации матриц вертикальных и горизонтальных компонент заданы на рис. 3. Здесь паттерн «1 0 -1» повторяется по строкам или столбцам. При п = 6 было бы двукратное повторение от левого верхнего края.
а)
1 0 -1 1 0
1 0 -1 1 0
1 0 -1 1 0
1 0 -1 1 0
1 0 -1 1 0
б)
1 1 1 1 1
0 0 0 0 0
-1 -1 -1 -1 -1
1 1 1 1 1
0 0 0 0 0
Рис. 3. Начальные состояния компонент V (а) и g (б) для матрицы 5 х5 . (начальное значение Ь -единица). В табл. 6 приведены правила перехода схемы «1:0:1».
Таблица 6. ЛФП схемы «1:0:1»
№ п/п Условие перехода Формула перехода
v = 1 s := sT ; v := 0
1 b = 1 v = 0 s = s ; v := -1
v = -1 s = ;v := 1 b := b
g = 1 s = ; g := 0
2 b = 0 g = 0 s = s ; g :=-1
g = -1 s = s^; g := 1
Легко видеть, что схема «1:0:1» состоит из шести элементарных итераций.
4. КА-формализация схемы «sin»
Определение синусоидального алгоритма в терминах клеточных автоматов также предполагает введение расширенного поля с единичной окрестностью фон-Неймана. Реализация движения элемента по функции sin х, в данном случае, состоит из четырех периодически повторяющихся итераций, представляющих собой сдвиги в условно положительном и отрицательном направлениях и сохранение текущей позиции (отсутствие перемещения). По сравнению со схемой «1:0:1» шаблон движения здесь «0 1 0 -1», что дает 8 этапов (итераций) на один проход. Эволюция слоя данных для КА размера 6 х 6 представлена в табл. 7.
Таблица 7. Проход схемы «sin» для матрицы 6 х 6
Номер итерации Слой данных
1 7 13 19 25 31
2 8 14 20 26 32
0 (исходная матрица) 3 9 15 21 27 33
4 10 16 22 28 34
5 11 17 23 29 35
6 12 18 24 30 36
1 8 13 24 25 32
2 9 14 19 26 33
1 3 10 15 20 27 34
4 11 16 21 28 35
5 12 17 22 29 36
6 7 18 23 30 31
Номер итерации Слой данных
1 8 13 24 25 32
33 2 9 14 19 26
2 3 10 15 20 27 34
11 16 21 28 35 4
5 12 17 22 29 36
31 6 7 18 23 30
33 8 7 24 19 32
3 2 13 14 27 26
3 11 10 9 20 35 34
5 16 15 28 29 4
31 12 21 22 23 36
1 6 17 18 25 30
32 33 8 7 24 19
3 2 13 14 27 26
4 10 9 20 35 34 11
5 16 15 28 29 4
36 31 12 21 22 23
1 6 17 18 25 30
32 6 8 14 24 30
3 33 13 35 27 19
5 10 2 20 28 34 26
5 9 15 21 29 11
36 16 12 18 22 4
1 31 17 7 25 23
32 6 8 14 24 30
33 13 35 27 19 3
6 10 2 20 28 34 26
11 5 9 15 21 29
36 16 12 18 22 4
31 17 7 25 23 1
31 6 35 14 23 30
32 13 20 27 24 3
7 33 2 9 28 19 26
10 5 12 15 34 29
11 16 7 18 21 4
36 17 8 25 22 1
6 35 14 23 30 31
32 13 20 27 24 3
8 26 33 2 9 28 19
10 5 12 15 34 29
16 7 18 21 4 11
36 17 8 25 22 1
В табл.8 представлены компоненты, определяющие состояния ячеек в схеме.
Название: Флаг перемещения Флаг направления движения Флаг итераций Регистр данных
Обозначение: т d и
Множество значений: т = {-1;0;1} d = {1;2;3} П = {0;1;2;3} Произвольное
Описание: Определяет направление сдвигов элементов на текущей итерации. Обпечивает сдвиги как по столбцам, так и по строкам. Определяет итерацию смены основных условий изменения двух других флагов Хранит элементы матрицы м=| |м,||.
Флаг перемещения т однозначно задает сдвиги элементов, в то время как компонента d обеспечивает работу алгоритма как по столбцам, так и по строкам. Оба флага меняют значения друг друга на каждой итерации, и для их корректного взаимоопределения вводится компонента й. ЛФП схемы приводится в табл.9.
Таблица 9. ЛФП схемы <^п»
№ п/п Условие перехода Формула перехода
d = 1 т := 1
1 d = 2 т := 0
d = 3 т := -1
(тек = -1) V (теу = 1) d := 2
¡1 = 1 теуФ -1 d := 1
теу= -1 d := 3
(теП = тоу^ ) а (те^ = тоу) Н mod2 = 0 теу= -1 d := 1
теу= 1 d := 3
и = 3 теу^ = 1 d := 1
2 тек = 0 теу^ Ф 1 d := 3
и = 1 теу^ Ф -1 d := 1
теу^ = -1 d := 3
(тек^ = теу^ ) а (теу^ = теу) ¡1 mod 2 = 0 теу^ = -1 d := 1
теу^ = 1 d := 3
и = 3 теу= 1 d := 1
теуФ 1 d := 3
3 и := (# + 1)mod4
т = -1 (т- =-1) а (т1 =-1) 5 := 5-
(т^_ = -1) а (т^ = -1) 5 :=
4 т = 0 5 := 5
т = 1 (щ = 1) а (т^ = 1) 5 := ^
(т^_ = 1) а (т^ = 1) 5 := ^
5. Временная сложность КА алгоритмов
Эмпирически установлено, что через некоторое число проходов алгоритмы переводят матрицу в исходное состояние, что принято нами за условие останова алгоритма. Это следует теоретически из конечности числа перестановок, которое для матрицы п х п для попарно различных элементов не превосходит (п2)!. Уточним, что период алгоритма N вводится не только для компоненты данных s, а также и для флагов. То есть пермутация совпадает с исходной матрицей, если обе компоненты всех её элементов принимают начальные значения. Если же для флагов совпадения нет, то такую ситуацию назовем самопересечением. В схемах «1:1» и «1:2» для матриц четных порядков период линеен, а для нечетных имеет более сложную зависимость от п. Чем длиннее период, тем лучше с точки зрения криптографии; и здесь, что несколько необычно, мы должны не сокращать, а увеличивать сложность алгоритма.
Кратко рассмотрим сначала периоды схемы «1:1». Для четных порядков матриц период равен п /2, а для нечетных порядков была доказана [19] следующая зависимость по закону наименьшего кратного: N (п ) = ЬСМ (3,5,... ,2п -1). Несмотря на быстрый рост с
увеличением п, общее число пермутаций составляет лишь ничтожную долю из возможных. Известна асимптотика:
LCM (1,2,..., m) = g (m) = exp
^ 'm(m +1) ^ m(m +1) ^
: exp(m)
2 2
V У
Легко получить, что g(2n -1) = LCM (N(n) ,2g (n)), g (m) = 0(mod (g ([m / 2]))) . И если r = [log2 n], то тогда g(2n -1) = 2rgf (n). Откуда получаем оценку для периода
N (n) : ^exp (2n) в базовой схеме. n
6. Периоды схемы «1:2»
Эмпирические зависимости N(n) в схеме «1:2» удобно показать, представив n как совокупность классов вычетов по делителю 4 (табл. 10).
Таблица 10. Периоды схемы «1:2» в зависимости от к eN
к N (4к -1) Каноническое разложение N (4к -1) N (4к) N (4к + 1) Каноническое разложение N (4к +1) N (4к + 2)
1 15 5 • 3 1 45 32 • 5 3
2 385 5 • 7-11 2 45 32 • 5 5
3 1309 7 • 11-17 3 273 3^ 7 43 7
4 1035 3-2 • 5 • 23 4 459 33 47 9
5 6061 11-19 • 29 5 231 3^ 7 41 11
6 10465 5•743^23 6 975 3^52 43 13
7 5535 32 • 5^41 7 1305 32 • 5 • 29 15
к N (4к -1) Каноническое разложение N (4к -1) N (4к) N (4к + 1) Каноническое разложение N (4к +1) N (4к + 2)
8 24769 17•31-47 8 561 3 -11-17 17
9 35245 5•7-19•53 9 2109 3-19 - 37 19
10 16107 3-7-13-59 10 2583 32 - 7 - 41 21
11 64285 5-13-23-43 11 1035 32 - 5 - 23 23
12 83425 52 - 47 - 71 12 3675 3-52-72 25
13 35343 32 - 7-11-17 13 4293 34 - 53 27
14 132385 5-11-29 - 83 14 1653 3-19 - 29 29
15 162781 . 31-59 - 89 . 15 5673 3-31-61 31
16 65835 32 -5-7-11-19 16 6435 32 - 5-11-13 33
17 236845 5 - 7 - 67-101 17 2415 3 - 5 - 7 - 23 35
18 281089 37 - 71-107 18 8103 3 - 37 - 73 37
19 110175 3-52-13-113 19 9009 32 - 7-11-13 39
20 385441 7-17 - 41-79 20 3321 34 - 41 41
21 446125 53 -43-83 21 10965 3-5-17 - 43 43
22 170955 32 - 5 - 29-131 22 12015 33 -5-89 45
После к = 22 сказывается длительность расчета на обычном компьютере (более шести часов для значения N (4к -1) при к = 23 ). Следует, однако, заметить, что речь идет о непосредственном расчете, т.е. без использования постулата о том, что общий период есть наименьшее общее кратное периодов индивидуальных элементов.
Матрицы нечетных порядков имеют более сложную структуру зависимости N(п), что можно заметить по циклическим спадам рассчитанных значений для п = 4к +1 и п = 4к — 1 . Полученные величины периодов алгоритма для нечетных п, в общем случае, зависят от к следующим образом: ^4к +1) : к2, ^4к — 1) : к3. На основе табл.9 эмпирически выведем закон, выполняющийся для начального участка натурального ряда точно:
Таблица 11. N(к) и N(n) схемы «1:2»
п(к) Условие для к N (к) N (n)
4к нет к n 4
4 к + 2 нет 2 к +1 n 2
4 к +1 кmod3 Ф 2 (4к + 1)(6к + 3) 3 n(n + 1) 2
к mod 3 = 2 (4к + 1)(2к +1) 1 n(n + 1) 2
4 к-1 к mod 3 Ф 1 (2к + 1)(4к -1)(6к -1) 1 n(3n + 1)(n + 3)
к mod 3 = 1 1(2к + 1)(4к -1)(6к -1) —n(3n + 1)(n + 3) 12
На рис. 4 представлены в двойном логарифмическом масштабе эмпирические зависимости N (п) для всех классов вычетов. Из 4 ветвей две разделяются по две подветви, отличающихся множителем при полиноме.
Рис. 4. Зависимости рассчитанных значений периодов схемы «1:2» от к на логарифмических осях. Угловой
коэффициент прямой Д ^ ^(к) равняется показателю степени d в выражении N (к) : к1.
Д lg к
7. Периоды схемы «1:0:1»
N(п) в схеме «1:0:1» проще показать, представив п как совокупность классов вычетов по делителю 3. Полученные величины периодов алгоритма зависят от к следующим образом: N(3^ : к, ^3к +1) : к2, ^3к + 2) : к3. Для классов п = 3к +1 и п = 3к+2 точная структура зависимости сложна и требует выделения подветвей для к: к(т),т Е N. В табл. 12 представлен периода схемы «1:0:1» для участка натурального ряда п е[3;182] по результатам непосредственного расчета.
Таблица 12. Рассчитанные периоды схемы «1:0:1» в зависимости от к ёИ
к N (3к) N (3к + 1) N (3к + 2)
1 2 30 40
2 4 12 126
3 6 70 360
4 8 48 110
5 10 126 156
6 12 40 2730
7 14 198 1904
8 16 120 1710
9 18 286 1680
10 20 84 8602
11 22 390 1800
к N (3к) N (3к + 1) N (3к + 2)
12 24 224 13338
13 26 510 4060
14 28 144 6510
15 30 646 1056
16 32 360 27370
17 34 798 2664
18 36 220 37050
19 38 966 21320
20 40 528 16254
21 42 1150 13860
22 44 312 62698
23 46 1350 11760
24 48 728 79050
25 50 1566 22048
26 52 420 2970
27 54 1798 54264
28 56 960 119770
29 58 2046 10980
30 60 544 144522
31 62 2310 79040
32 64 1224 57486
33 66 2590
Легко видеть, что для первого класса (п = 3к) период алгоритма линеен и равен 2
N(п) = — п . Период для класса п = 3к +1 подразделяется на три подветви в зависимости от
к (т) (рис. 5). Для данного класса вычетов алгоритм, предположительно, имеет квадратичную сложность.
ю
20
-Н(Эк*1), кГгп 1=2гп-Т
в-'
■в—
30
-10
50 г>
60
го
№
90
100
Рис. 5. Зависимости рассчитанных значений периодов схемы «101» от п = 3 к +1.
В классе п = 3к + 2 зависимость N(п) делится уже на шесть ветвей по разным законам к (т) (рис. 6). В общем случае, период алгоритма для данного класса вычетов по 3 имеет кубическую сложность.
Рис. 6. Зависимости рассчитанных значений периодов схемы «1:0:1» от п = 3к + 2 . Красными маркерами отмечены точки (табл. 13), в которых нарушаются эмпирические формулы N (п).
Таблица 13. N(т) и N(n) схемы «1:0:1»
п(к) к (т) N (т) N (п)
множитель формула
3к т 2т 2 3 п
2т -1 2(2т + 3)(2т +1) 2 9
3к +1 4 т - 2 4т (2т + 1) 1 18 (п + 5)( п +11)
4т 8(т + 1)(2т +1) 1 9
12т 2(12т + 7)(24т + 3)(12т +1)
12т + 4 2(12т - 7)(24т - 13)(12т -1) 2 27
3к + 2 12т + 6 2(12т - 5)(24т - 9)(12т +1) (п + 19)(2п + 5)(п + 1)
12т + 10 2(12т + 5)(24т - 1)(12т -1)
12т + 3 24(3т - 2)(8т - 5)(6т -1) 1
12т + 7 8(3т - 1)(24т - 7)(6т + 1) 27
n(k ) k (m) N (m) N (n)
множитель формула
12т + 2 6(4m - 3)(24m -17)(4m -1) 2 81
12m + 8 6(4m + 1)(24m - 5)(4m -1)
12m + 1 4(3m + 2)(24m + 5)(6m +1) 1 (n + 19)(2n + 5)(n + 1)
3k + 2 12m + 9 12(3m + 1)(6m - 1)(8m -1) 54
12m +11 24(2m + 3)(24m + 25)(m +1) 1 81
12m + 5 12(m + 1)(24m + 13)(2m +1) 1 162
Однако из эмпирически полученных зависимостей N(3к + 2) выпадают некоторые точки (табл. 14)
Таблица 14. Выпадающие точки в схеме «1:0:1» для п = 3к + 2
№ ветви n k (m) m1 N (n ) expect N (n ) Jact N (n ) expect N (n ) fact
1(синяя) 14 12m + 4 1 1210 110 11
2(оранжевая) 47 12m + 3 2 11616 1056 11
3(желтая) 80 12m + 2 3 32670 2970 11
4(фиолетовая) 113 12m+1 4 64372 5852 11
1(синяя) 146 12m 4 533610 48510 11
Разность порядков матриц последующей и предыдущей точек всегда равна 33. Заметно также отношение ожидаемого значение периода N(n)expect и фактического N(n)^,
равное 11 для всех особых точек. Данные точки являются членами арифметической прогрессии с начальным членом 14 и шагом 33.
8. Периоды схемы «sin»
Данный алгоритм имеет запутанную структуру зависимостей N(n), однако все n могут быть также разделены на четные и нечетные. Он как бы объединяет особенности базового алгоритма и предыдущих (табл. 15, 16). N для всех матриц с четным П в синусоидальной схеме структурно разбивается на две ветви, одна из которых состоит из трех под-ветвей.
Таблица 15. N(n) в схеме «sin» для групп n = 4k и n = 4k + 2
n( k ) k (m) N (n)
множитель формула
4k - 1 (9 / 4)n
2m +1 1
4k + 2 4m+2 1 4 9 — (n + 6)(3n + 2) 16
4m 1 2
п N (п) (в итерациях) N (п) (в проходах)
Схема «1:1» Схема ««ш» N (п) ( )!'п -100% N (п)1:1 Схема «1:1» Схема «шп»
3 60 24 40 15 3
5 1260 1,120 88.(8) 315 140
7 180180 6,336 3.516484 45,045 792
9 3063060 2,688 0.087755 765,765 336
11 58198140 15,840 0.027217 14,549,535 1,980
13 6692786100 9,853,200 0.147221 1,673,196,525 1,231,650
15 582,272,390,700 - - 145,568,097,675
17 18,050,444,111,700 2,545,920 0.0000017 4,512,611,027,925 318,240
19 667,866,432,132,900 - - 166,966,608,033,225 -
21 27,382,523,717,448,900 - - 6,845,630,929,362,225 -
23 1,177,448,519,850,302,700 - - 294,362,129,962,575,675 -
Примечание: прочерки означают крах непосредственного численного расчета из-за сверхбольших чисел
9. Индивидуальные траектории элементов
Для лучшего понимания формул для периода рассмотрим теперь траектории движения отдельных элементов матрицы (или, для краткости, индивидуальные траектории -УИТ) и их временные длины в проходах (ДИТ). Поле КА с замыканием удобно представить как бесконечную плоскую периодическую решётку (рис.7). Такое представление будем называть расширенным.
Рис. 7. Расширенное представление для матрицы 4 х 4
УИТ в схеме «1:1» напоминают путь бильярдного шара, а граница поля играет роль стенки. ДИТ зависит от начального положения элемента. Период алгоритма, если не учитывать флаги, определяется произведением ДИТ, а точнее - их наименьшим кратным.
Рис. 8 Индивидуальные траектории некоторых элементов матрицы 5 х 5 в расширенном представлении для схемы «1:1». Точками отмечены начальные положения элементов. Рядом с каждой траекторией элемента
указана ее ДИТ.
Напомним, что возможны самопересечения на пути ячейки, когда ее позиция совпадает с начальной, но направление движения отлично от исходного. Ромбовидные траектории, по сути, являются ключом к доказательству формулы периода.
10. Индивидуальные траектории элементов в схеме «1:2»
УИТ в схеме «1:2» по совокупности менее разнообразны по ДИТ, чем в схеме «1:1». Все элементы матриц четных порядков сдвигаются по одной траектории, период которой равен, соответственно периоду алгоритма. Общее количество УИТ для матриц нечетных порядков равно четырем, начиная с п = 5 (для п = 3 множество УИТ состоит из трех элементов - [1;3;5] ). На рис. 10-14 представлены некоторые УИТ матриц, порядки которых принадлежат разным классам вычетов по делителю 4. Эмпирически было установлено, что в множестве ДИТ всегда присутствует единица, а другие три - зависят от п . Длина
одной из них равна п . Элемент с такой УИТ никогда не пересекает горизонтальных границ расширенного поля (рис.12). Две другие назовем «короткой» и «длинной» УИТ, обозначив их длины как /3, /4, а траектории длины 1 и п как 1г .и 12 соответственно. В табл. 17 приведены эмпирические формулы ДИТ для нечетных порядков матриц:
1 9 1 5 1 1 5 5
15 15 15 9 15 15 15 9 15
9 9 9 9 9 9 9 9 9
9 9 9 9 9 9 9 9 9
9 9 9 9 9 9 9 9 9
9 9 9 9 9 9 9 9 9
9 15 9 9 9 15 9 9 9
9 9 9 5 9 9 9 5 9
15 1 15 1 15 1 15 15 15
Рис. 9. Таблица ДИТ для матрицы 9 х 9 в схеме «1:2». Каждому элементу матрицы соответствует своя ДИТ.
2
»
* -5Г
р г
1_
*
я
ь
Рис. 10. Траектории некоторых элементов матрицы 8 х 8 (п = 4к) в схеме «1:2» на расширенном поле.
Рис. 11. Траектории некоторых элементов матрицы 10 х10 (п = 4к + 2) в схеме «1:2» на расширенном поле.
Рис. 12. УИТ матрицы 9 х 9 (п = 4к +1) в схеме «1:2» на расширенном поле.
На рисунках отсутствуют траектории единичной длины из-за их граничного положения. При визуализации происходит пересечение и слияние с другими УИТ.
Таблица 17. Эмпирически полученные формулы ДИТ для нечетных п в схеме «1:2»
п(к) 4(п) /2(п) 1з(п) 14 (п)
4 к-1 1 п 2(п+3) 3 -(п +1) - 1 2
4 к +1 1 п 2(п+» 3 2(п +1)
Рис. 13. «Длинные» траектории элементов в схеме «1:2» для матрицы 11 х 11 (п = 4к -1) на расширенном
поле с ДИТ 1(3п +1).
Траектории длины 1 не отображены из-за особенностей их расположения в граничных столбцах и строках.
11. Индивидуальные траектории элементов в схеме «1:0:1»
Множество УИТ в схеме «1:0:1» также не обладает большим количеством элементов, однако таблица ДИТ имеет некие регулярности. Матрицы порядка п = 3к +1 имеют 4 зна-
чения ДИТ, среди которых всегда есть 1 и 2. Два оставшихся значения зависят от п линейно:
1 2 Уп = 3к +1: I(п) = ~(п + 5), /2(п) = ~(п +11)
Стоит также отметить распределение траекторий по элементам для матриц п = 3к +1, показанное на рис. 15:
5 14 14 5 14 14 5 14 14 1
14 14 2 2 14 2 2 14 2 14
14 2 14 2 2 14 2 2 14 14
5 2 2 5 2 2 5 2 2 5
14 14 2 2 14 2 2 14 2 14
14 2 14 2 2 14 2 2 14 14
5 2 2 5 2 2 5 2 2 5
14 14 2 2 14 2 2 14 2 14
14 2 14 2 2 14 2 2 14 14
5 4 14 5 14 14 5 14 14 5
Рис. 15. Таблица ДИТ для матрицы 10 х10 ( п = 3 к +1) в схеме «1:0:1»
В первых и последних столбцах и строках матрицы повторяются длины I и /2: после элемента с ДИТ I следуют два элемента с длинами /2. Данная периодичность нарушается только последним элементом первой строки, ДИТ которого для всех рассмотренных порядков матриц (до п = 100 ) равняется единице. Также можно увидеть периодичность ДИТ в столбцах, не являющихся граничными и не включающих первый и последний элементы. Столбцы ДИТ делятся на классы вычетов по делителю 3 в зависимости от номера столбца ] . То есть каждый повторяется через два.
На рис. 16 изображены все четыре УИТ для матрицы 10 х10. Данное отображение траекторий не совсем наглядно, так как невозможно показать итерации с отсутствием сдвига. Однако легко видеть как меняется направление движения элемента на границах матрицы. Из-за нарушения последовательности изменения флагов перемещения ячейка совершает колебательные движения в горизонтальном и вертикальном направлении пока не окажется в неграничном столбце или строке.
Рис. 16. УИТ для матрицы 10 х10 (п = 3 к +1) в схеме «1:0:1»
Логично предположить, что матрицы порядка п = 3к + 2 имеют другую структуру индивидуальных траекторий. Количество элементов в множестве ДИТ становится постоянным только с п = 20 и равняется пяти. До этого момента не прослеживается четкая зависимость: матрицы порядков 5, 14 17 имеют только четыре разных ДИТ. Для к > 5 множество ДИТ предположительно описываются эмпирически полученными формулами (табл. 18):
Таблица 18. Эмпирически полученные формулы ДИТ для п = 3к +1 и п = 3к + 2 в схеме «1:0:1»
п(к) А(п) /2(п) /з(п) /4 (п) /5 (п)
3 к +1 1 2 1 (п + 5) 2 (п+11)
3к + 2 2 3 (п+■> 1 (п +19) 3 ( п +1) 3 (п+5)
Примечание: формула для / (3к + 2) справедлива для к > 1, в то время как остальные зависимости / (3к + 2), / > 1 выполняются при к > 5 .
Рассмотрим распределение траекторий по элементам для матриц п = 3к+2, показанное на рис. 17:
10 10 10 10 10 10 10 10 10 9
10 9 2 2 8 2 2 8 2 2 8
10 2 9 2 2 8 2 2 8 2 10
10 2 2 4 2 2 10 2 2 10 8
10 8 2 2 9 2 2 8 2 2 8
10 2 8 2 2 9 2 2 8 2 10
10 2 2 10 2 2 4 2 2 10 8
10 8 2 2 8 2 2 9 2 2 8
10 2 8 2 2 8 2 2 9 2 10
10 2 2 10 2 2 10 2 2 4 8
9 8 8 10 8 8 10 8 8 10 9
Рис. 17. Таблица ДИТ для матрицы 11x11 ( п = Зк + 2) в схеме «1:0:1»
Распределение траекторий для класса п = Зк+2также имеет регулярности. Все элементы первого столбца и строки имеют длину траектории /4, за исключением начального
и последнего, ДИТ которых равны соответственно ^ и /3. В остальной таблице равными оказываются строки и столбцы с индексами г = j для которых выполняется условие г шоёЗ = 2. Иначе попарное равенство элементов в них нарушается последними ячейками, которые имеют ДИТ /4 или /2. Также во множестве УИТ матриц с п = Зк + 2 всегда присутствует 2.
На рис. 18 изображены все пять УИТ для матрицы 11x11.
Рис. 18. УИТ для матрицы 11x11 ( п = Зк + 2) в схеме «1:0:1». Цифры обозначают число проходов.
12. Индивидуальные траектории элементов в схеме «sin»
Все ДИТ для матриц четных порядков ( n = 4kи n = 4k+2 ) на исследованном натуральном линейно зависят от п. Также, множество длин траекторий класса n = 4k+2 всегда содержит 1 и 3.
Таблица 19. Эмпирически полученные формулы ДИТ для n = 4k и n = 4k + 2 в схеме «sin»
n(k ) li(n) n) l3(n) I4 (n) h(n) ls(n)
4k n 3 9
— —n —n - - -
4 4 4
4 k + 2 1 3 1 „ -n + 2 4 3 „ -n+2 4 3 . -n+6 4 9 . -n+6 4
Рис. 19. УИТ для матрицы 6 х 6 ( n = 4k + 2) в схеме «sin»
В классах нечетных порядков мощность множества ДИТ линейно по п. Для обнаружения закономерностей напрашиваются рекуррентные соотношения. При переходе к ^ к +1 в множество ДИТ добавляются два элемента, на 9 большие двух наибольших элементов прежнего множества ДИТ. Уже имевшиеся ДИТ увеличиваются обычно на единицу, кроме некоторых, инкремент возрастания которых равняется двум или трем. Легко видеть, что для такого нарушения характерна периодичность:
Таблица 20. Особенности изменения множества ДИТ с увеличением к в классе п = 4к +1, к е Г
к Новые Количество Число нарушений ДИТ, нарушающие Инкремент
ДИТ ДИТ возрастания ДИТ закон возрастания нарушения
1 1;4;5;7 4 0
2 14;16 6 2 6; 7 2;2
3 23;25 8 2 9;10 2;2
4 32;34 10 1 13 3
5 41;43 12 1 16 3
6 50;52 14 2 18;19 2;2
7 59;61 16 2 21;22 2;2
8 68;70 18 1 25 3
9 77;79 20 1 28 3
10 86;88 22 2 30;31 2;2
11 95;97 24 2 33;34 2;2
12 104;106 26 1 37 3
13 113;115 28 1 40 3
14 122;124 30 2 42;43 2;2
15 131;133 32 2 45;46 2;2
16 140;142 34 1 49 3
Опишем рекурсию более строго. Пусть для некоторого к дано множество всех ДИТ:
Л^ = {/ £ ЩЛк = сата{кк)
Тогда для класса п = 4к +1, к > 0 имеют место следующие свойства:
• \ = {1,4,5,7} - начальное множество
• = 2, Ьк+1 = Ьк + 2 - рекурсия числа ДИТ
• 1 е Л - всегда есть ДИТ=1
• {9к - 2,9к - 4} е Лк+1 - порождение новых элементов
• Зф е Л(к): |ф| = 1V 2 - наличие нарушителей
(к = 4ш) V (к = 4т +1) ^ (ф = {Зк +1} л (Зк + 4 е Лк+1 ))
• - правило на (к = 4т - 2) V (к = 4т -1) ^ (ф = {Зк, Зк +1} л ({Зк + 2, Зк + З} е ^ ))
нарушения
• (Ук : I е Лк ) ((/ £ (Фк и {1})) ^ ((I +1) е Лк+1)) - регулярное правило
Для класса п = 4к + З особенности изменения элементов множества ДИТ полностью идентичны классу п = 4к +1. Различие только в начальном множестве и в том, неизменны первые два ДИТ {1,3}, а не только ДИТ= 1.
Таблица 21. Особенности изменения множества ДИТ с увеличением п в классе п = 4к + 3, к е Г
k Новые Количество Число нарушений ДИТ, нарушающие Инкремент
ДИТ ДИТ возрастания ДИТ закон возрастания нарушения
1 1;3;8;9;11 5 0 - -
2 18;20 7 2 10;11 2;2
3 27;29 9 2 13;14 2;2
4 36;38 11 1 17 3
5 45;47 13 1 20 3
6 54;56 15 2 22;23 2;2
7 63;65 17 2 25;26 2;2
8 72;74 19 1 29 3
9 81;83] 21 1 32 3
10 90;92 23 2 34;35 2;2
11 99;101 25 2 37;38 2;2
12 108; 110 27 1 41 3
13 117;119 29 1 44 3
14 126;128 31 2 46;47 2;2
15 135;137 33 2 49;50 2;2
16 144;146 35 1 53 3
а) б)
Рис. 20. УИТ для матриц 5 х 5 (n = 4 k +1 )(а) и 7 х 7 (n = 4 k + 3 )(б) в схеме «sin»
13. Анализ динамики схемы «1:2» с помощью «метрик перемешанности»
Каждый ход КА дает пермутацию матрицы относительно начальной (здесь и далее считаем все элементы различными). Для описания динамики КА желательно иметь числовую характеристику «удалённости» текущей перестановки от начальной. Следуя [ 19,21], введем «метрики» тт и 1т. Пусть в некоторой матрице А порядка п записаны попарно
различные элементы са, а е 1, п2 , а В является пермутацией А . Тогда метрика тт на этих двух матрицах определяется формулой
1 п2 тт = —^ ^а ,
n
а=1
где яа =| /2 - ц | + 172 - 7 |, если (ц, 7 ) — это индексы элемента са в матрице А , а (/2, у2) — индексы са в матрице В . Выражение ъа имеет смысл расстояния между положениями элемента са в матрицах А и В по метрике Манхэттена. «Метрика» 1т вводится сначала
для одномерного массива. Считая, что начальное его состояние имеет вид ак = к, к е 1, т , /т-метрика определяется формулой:
lm =
m +1
lt + h
где Ц — максимальная длина возрастающей подпоследовательности, а ^ — максимальная
длина убывающей подпоследовательности. Чтобы применить lm к матрице, необходимо сначала сопоставить матрице массив, и, следуя соглашению в MATLAB, запишем её элементы в столбец.
Максимумом метрики mm в схеме «1:1» для четных n обладала конфигурация, соответствующая диагональной перестановкой четырех блоков исходной матрицы [19]. Для матриц нечетных порядков в схеме 1:1 были выявлены конфигурации поворота против и по часовой стрелке, а также отражения относительно центра, являющиеся максимумами и минимумом lm метрики соответственно.
Ниже приведены графики mm и lm, соответствующие схеме «1:2» для матриц различных n(к). Рассмотрим вначале четные порядки n = 4к и n = 4k + 2. Динамика обеих «метрик» показана на рис. 21 и рис. 22 (с детализацией до шага iter, а не прохода).
300 STS 2» 225 \
1BQ1
100
71
25
гт ~Т" т~ ГХ~г -тт- г~ — - — — шщ .... - -
\ тлп . г = 24 (4 к)
/ / \
> г
\ I
/
А в -В [А V /
/ \
м 1 IL _У 1
/
и 1 1 1 1 I | I „ J • \
D I 2 i 1 5 6 7 '! В 10-11 111J 14 15 1fi 17 id li 2D 21 22 23 2i 2f 2Ь 27 212В 3d 31 32 33 S4 35 П6
а)
б)
Рис. 21. Динамика «метрик» матрицы 24 х 24 в схеме «1:2» (п = 4к)
а)
г\ (т, п= 22 (4к+2}
ф Л
йЧйГ ■
---
1 ¿Л
10
30
Шаг
б)
ее
ТО
Рис. 22. Динамика «метрик» матрицы 22 х 22 в схеме «1:2» (п = 4к + 2)
На рисунках красной звездочкой и греческой буквой а отмечены значения метрик метрик mm для матриц с n = 4k и n = 4k+2. Эти конфигурации есть перестановка двух блоков, образованных вертикальным разрезом матрицы пополам. Для матриц порядков n = 4k данная конфигурация соответствует глобальному и локальному минимумам mm и lm метрик соответственно. Легко видеть, что двойная блочная перестановка достигается на середине периода (iter = 18). Матрицы порядка n = 4k + 2 также имеют минимумы обоих характеристик на шаге с блочной перестановкой, однако данная конфигурация дос-
(9 N +1) (N + 2) тигается за iter =--- шагов для четных и iter =---для нечетных k . (8,9)
Анализ матриц нечетных порядков с помощью «метрик» показал, что на всем периоде, включая максимумы и минимумы, отсутствуют какие-либо особые пермутации. Соответствующие графики представлены на рис. 23. Чуть ниже (рис. 24-27) расположены гистограммы распределения пермутаций по метрикам: «доля от общего числа перестановок -значение метрики». В верхнем левом углу рисунков, для наглядности, также приведены lm и mm на периоде. Заметим, что по сравнению со схемой «1:1», где гистограмма имеет вид, похожий на три прореженных дискретных нормальных распределения [19], гистограммы схемы «1:2» более сложны и нерегулярны. Прежде всего обнаруживается спектральный характер и запрещенные или малонаселенные зоны значений «метрик».
а)
Е00
150
1М
тт, п= =21 (4к+1)
. - '('. (и*. г ■ . . _ ■ Г * А - 1 ■ ,. "цм*.; I,. . (1! 11" ■ * •, ■ ' - . •,»!(
» ■ р * 1 т ' Г" 1 1 . ' .1 - <р 1 _» " 1 ; - 14 • .:, - * I* 1 ■ 1 ■ - * 1 л ■ * \ щ* |В ;« *;" *. ' ■ -1„ - • ■ ■ / •« . ' Г - | ■ 1 . \ .. к • * / ~ V'
" - М - > * ■ * *. г ' -1 ,■ Ч | Ч > / " м Зь «Г ;
■ *
- ' " 1 / ^ • ' * *. Л ' : ■ г, .. - -* ■ 1 * . - ч ■ * " 1
,. ' ' • 1 " 1 ■ - ^ * ь , 1 ■ 4 * * « г , , , ; ■ ■ * * \ • * ■ " ;; л -' ■4 1 * ь ■*;
1 . 1 - » . * ' * , " - * * ■ г | * ' . г - * 1 9 Л ^7777 * . - л "в
, ( ; 1 г * * , # 1 • ■ * : . - ■ 1'!".
. * р I . ■ . * г
. > '.....* - * ■ 1 * ' 1 _
200
<Ю0
ВОО
МО
10И
12НЭ
1400
Шаг
в)
Рис. 23. Динамика метрик последовательных конфигураций матрицы 11x11 в схеме «1:2» (п = 4к-1), кmod3 Ф1 и тт метрика матрицы 21х 21 (п = 4к +1), кmod3 = 2
а)
б)
Рис. 24. Распределения последовательных конфигураций матрицы 15 х 15 в схеме «1:2» (n = 4k -1), k mod 3 = 1 по значениям метрик.
а)
б)
Рис. 25. Распределения последовательных конфигураций матрицы 17 х 17 в схеме «1:2» (п = 4к +1), к шо<!3 Ф 2 по значениям метрик.
а)
б)
Рис. 26. Распределения последовательных конфигураций матрицы 21 х 21 в схеме «1:2» (п = 4к +1), к то<3 = 2 по значениям метрик.
14. Анализ динамики схемы «1:0:1» с помощью «метрик
перемешанности»
При анализе динамики пермутаций с помощью «метрик», в схеме «1:0:1» заметна некоторая периодичность значений данных характеристик. Для 1т-метрики спектральный характер усиливается. Матрицы группы п = 3к + 2 на периоде алгоритма имеют повторяющиеся значения обоих метрик (рис. 27). При анализе схемы «1:0:1» с помощью метрик перемешанности на периоде алгоритма не было обнаружено особых пермутаций.
а)
б)
Рис. 27. Динамики mm- и 1т-«метрик» в схеме «1:0:1» при и=71. Ветвь: n = 3к + 2, к = 12m + 11.
Разные ветви к (m ) данной группы имеют, соответственно, отличные друг от друга диапазоны mm - метрики, что, скорее всего, обусловлено порядком матриц. Также отличны и распределения пермутаций по метрикам, которые не могут быть примерно приближены каким-либо аналитическим выражением. Однако стоит заметить, что для всех ветвей группы n = 3к + 2 характерно преобладание на периоде пермутаций с высокими значениями mm-перемешанности (примерно 2/3 от максимального значения) (табл. 22).
Таблица 22. Плотность распределений пермутаций (гистограмма) по значениям mm-метрики для всех шести
веток класса n = 3к + 2 в схеме «1:0:1».
к (m)
12m
f (mm)
12m+1
Графики /^«метрики» имеют лишь несколько областей, где сосредоточены все последовательные конфигурации. Для данной характеристики тоже наблюдается периодичность значений, особенно заметная по выделяющимся локальным минимумам, не являющимся, однако, какими-либо особыми пермутациями (табл. 23).
Таблица 23. Плотность распределений пермутаций по значениям Im -метрики для всех шести веток класса
n = 3к + 2 в схеме «1:0:1».
к ( I )
f ( Im)
12m
к (I)
/ (1т)
12т + 3
12т + 11
Для некоторых ветвей диапазон значений /да-характеристики шире, что показано большим количеством областей значений. Анализ с помощью метрик «выпадающих точек» показал, что матрицы таких порядков имеют то же распределение последовательных конфигураций по метрикам, что и остальные матрицы в ветке (рис. 28).
Рис. 28. Гистограммы для матриц п — 47 и п - 59, принадлежащих выпадающим точкам ветви
п — 3к + 2, к — 12т + 3 в схеме «1:0:1».
15. Анализ динамики схемы «sin» с помощью «метрик перемешанности»
Характеристики перемешанности в группах (n=4k n=4k+2) довольно хаотичны и не описываются аналитическими выражениями. Для класса n=4k на половине периода алгоритма присутствует блочная перестановка как в базовом алгоритме, соответствующая максимуму /m-метрики (конфигурация а на рис.29).
25
20
15
Е Е
10
-*■ * я * « 1 1 1 * * , гпт, п=8 (4К)
* * f а > i Л в " . .. ■■ • i * Л . ' * * * * > * в Р- *
* ' " *
-М *я Mi «i
50
100
150
Шзг
а)
3.5
2 5
1.5
«я V V ■ Im, п-8 (4к)
* v i к # ■ ■ «i • • мам я Л 1* « «1 м
Я # * 1 И и 4 ■ Я * 1 • # * »Н Л * * • R I i* * * 1 Я В* # К
■ ■ <" * * V • я ■ * м * ■ # ч
*
О
50
100
150
Шаг
б)
Рис. 29. Перемешанности последовательных пермутаций матрицы 8 х 8 (n = 4k) в схеме «sin». Буквой а обозначено значение метрики для блочной перестановки из схемы «1:1».
Метрики нечетных классов имеют выделяющиеся глобальные максимумы и локальные минимумы. Однако для обоих групп на периоде отсутствуют особые конфигурации. Также для /m-метрики характерна большая дискретизация по значениям, что особенно видно на гистограммах:
а)
б)
Рис. 30. Распределения последовательных конфигураций матрицы 9 х 9 (п = 4к +1) в схеме <«1п» по
значениям метрик.
Общий вид гистограмм дада-перемешанности напоминает нормальное распределение, а на гистограммах /да-метрики заметны области значений, в которых отсутствуют какие-либо пермутации. Следует скорее говорить не о спектральности, а о прореженности распределения.
а)
б)
Рис. 31. Распределения последовательных конфигураций матрицы 11x11 (n = 4k + 3) в схеме «sin» по
значениям метрик.
Заключение
Компьютерное «наивное» исследование, результаты которого подробно рассмотрены в статье, позволяет лучше «почувствовать» свойства предложенных алгоритмов. Это открывает путь уже для теоретического исследования не только в отношении формулы периода и ДИТ, но и в отношении поведения (и эффективности) «метрик перемешанности», где требуется объяснение их абриса на гистограммах. Суммируем полученные эмпирические, а значит в отсутствии аналитического доказательства верные лишь для начального участка натурального ряда, утверждения.
Утверждение 0. Для всех схем период КА равен наименьшему общему кратному всех значений ДИТ. (Вроде просто/очевидно, но нужно еще и совпадение в слое флагов)
Утверждение 1. Для четных порядков п схема (1:1) дает блочную перестановку, а схема (1:2) дает отражение матрицы относительно вертикальной оси симметрии.
Утверждение 2. Для нечетных п в схеме (1:2) возможные ДИТ пробегают 4 значения, два из которых 1 и п (табл. 17). Неплохо заодно получить количества элементов с данными ДИТ.
Утверждение 3. Для четных п в схеме (1:2) период линейно зависит от п, а для нечетных -квадратично или кубически. Неплохо показать разницу множителя в три раза (табл. 11).
Утверждение 4. Для схемы 1(0)1 количество значений ДИТ L зависит (табл. 18) от класса вычета n по модулю 3: если остаток 0, L=1; если остаток 1, L=4; если остаток
2, L=5.
Утверждение 5. Для схемы 1(0)1 период зависит от порядка матрицы линейно, квадратично или кубически, что определяется классом вычета n по модулю три (табл.13). Внутри каждого класса есть вариации множителя перед полиномом.
Утверждение 6. Для схемы 1(0)1 существуют особые значения n, образующие арифметическую прогрессию с начальным членом 14 и разностью 33 (табл. 14), где этот коэффициент падает в 11 раз относительно ожидаемого общей формулой.
Утверждение 7. Для схемы sin период зависит от класса вычета n по модулю 4: если остаток равен 0, то закон роста линейный, если равен 2, квадратичный (табл. 15).
Утверждение 8. Для схемы sin выделяются сходные классы вычетов с остатком 1 и
3, где период подчинен сложному рекуррентному закону о числах ДИТ. Эти классы отличаются только начальным значением для рекурсии (табл. 20, 21).
Как мы полагаем, их доказательство носит технический характер. Лишь четвертая схема может доставить проблемы, связанные с общей нерегулярностью натурального ряда (в отношении простоты чисел, например). Непонятна, какая асимптотика периода у нечетных ветвей этой схемы. Идея базовой схемы изначально заключалась в конструкции наиболее простого КА, способного генерировать особые пермутации матриц. Нам не удалось за счет модификации базовой схемы повысить длину периода алгоритма, что было бы полезно с точки зрения криптографии. Но из этой неудачи возникает вопрос для по-настоящему серьезного математического исследования:
Какое наибольшее число пермутаций матрицы порядка n способен порождать до своего останова детерминированный клеточный автомат, заданный в поле nx n?
Оценка снизу следует из периода базовой схемы — exp(2n) / n . Оценка сверху тоже
( п2 Y2
очевидна - I ni!: — . Проблема, по сути, алгоритмическая и связана с ограничением
e
V e У
клеточного автомата, накладываемого на реализацию алгоритма.
Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта №20-07-00409.
Список литературы
1. Wolfram S. Random Sequence Generation by Cellular Automata // Advances in Applied Mathematics, 1986, vol. 7, Pp. 429 - 432.
2. Сухинин Б. М. Разработка генераторов псевдослучайных двоичных последовательностей на основе клеточных автоматов // Машиностроение и компьютерные технологии. 2010. №9. Режим доступа: https://cyberleninka.ru/article/n/razrabotka-generatorov-
psevdosluchaynyh-dvoichnyh-posledovatelnostey-na-osnove-kletochnyh-avtomatov (дата обращения: 25.08.2020).
3. Ярмолик В.Н., Мурашко И. А. Реализация генератора псевдослучайной последовательности на клеточных автоматах // Автоматика и вычислительная техника. 1993. №3. С. 9-13.
4. Bar dell P. Analysis of cellular automata used as pseudorandom pattern generators // Proceeings., International Test Conference. 1990. Pp. 762-768.
5. Girau B., Vlassopoulos N. Evolution of 2-dimensional cellular automata as pseudo-random number generators // Cellular automata: 10th Intern. conf. on cellular automata for research and industry: ACRI 2012 (Santorini island, Greece, Sept. 24-27, 2012): Proc. B.; Hdbl.: Springer, 2012. Pp. 611-622.
6. Мухамеджанов Д.Д., Левина А.Б. Генератор псевдослучайных чисел на основе клеточных автоматов // Научно-технический вестник информационных технологий, механики и оптики. 2018. №5. С 894-900. DOI: 10.17586/2226-1494-2018-18-5-894-900
7. S. Bilan, M. Bilan, S. Bilan, Novel pseudo-random sequence of numbers generator based cellular automata. // Information Technology and Security. - 2015. - Vol. 3, Iss. 1 (4). - Pp. 38-50. DOI: doi.org/10.20535/2411-1031.2015.3.1.57710
8. Bruno Martin, Patrick Solé. Pseudo-random Sequences Generated by Cellular Automata. International Conference on Relations, Orders and Graphs: Interactions with Computer Science, May 2008, Mahdia, Tunisia, Pp. 401-410.
9. Agrawal V, Bushnell M. Essentials of Electronic Testing for Digital, Memory, and Mixed-Signal VLSI Circuits. Springer, 2000. Pp. 712.
10. Hortensius P.D. Parallel random number generation for VLSI systems using cellular automata // IEEE Transactions on Computers. 1989. Vol. 28 (10). Pp. 1466-1473.
11. Fog A. Pseudo-Random Number Generators for Vector Processors and Multicore Processors. // Journal of modern applied statistical methods. 2015. №14(1) Pp 308 - 334.
12. Ключарёв П.Г. Построение случайных графов, предназначенных для применения в криптографических алгоритмах, основанных на обобщенных клеточных автоматах. // Математика и математическое моделирование. 2017. №3. С. 77-90. DOI: 10.24108/mathm.0317.0000076.
13. Конахович Г.Ф. , Пузыренко А.Ю. Компьютерная стеганография Теория и практика. К "МК-Пресс" 2006. Гл. 5. С. 92-97.
14. Chen Y. L., Lambooij E., Mennink B. How to Build Pseudorandom Functions From Public Random Permutations Cryptology. Cryptology ePrint Archive. Available at https://eprint.iacr.org/2019/554.pdf , accessed 09.07.2020.
15. Trevisan L. U.C. Berkeley — CS276: Cryptography: Notes for Lecture 15. Available at https://people.eecs.berkeley.edu/~luca/cs276/lecture15.pdf , accessed 09.07.2020.
16. Rubio, C. F., L. Encinas, S. H. White, Á. Rey and G. Sánchez. The Use of Linear Hybrid Cellular Automata as Pseudo Random Bit Generators in Cryptography. // Neural Parallel & Scientific Comp. 2004. №12. Pp 175-192.
17. Жуков Д. А. О порядке одной перестановки на элементах квадратной матрицы // Десятая международная научно-технической конференция «Безопасные информационные технологии»: труды. М.: Москва, Изд-во Московский государственный технический университет имени Н.Э. Баумана (национальный исследовательский университет) (Москва), 2019. С. 143-145.
18. Weichuang Guo, Junqin Zhao, Ruisong Ye. A chaos-based pseudorandom permutation and bilateral diffusion scheme for image encryption // Intern. J. of Image, Graphics and Signal Processing. 2014. Vol. 6. №11. Pp. 50-61.
19. Матюшкин И.В., Кожевников В.С Клеточно-автоматные алгоритмы пермутации матриц // ТРУДЫ МФТИ. 2019. Том 11, № 1. С.39-52.
20. Матюшкин И.В., Кожевников В.С. Клеточно-автоматный алгоритм пермутации матриц // Математика и математическое моделирование. 2019. №3. С. 1-24.
21. Деза Е.И., Деза М.М. Энциклопедический словарь расстояний: пер. с англ. М.: Наука, 2008. 444 с. [Deza E., Deza M M. Dictionary of distances. Amst.: Elsevier, 2006. 391 p.].
Mathematics & Mathematical Modelling
Electronic journal http://mathmelpub.ru ISSN 2412-5911
for Matrix I.V. Matyushkin1'2'*, P.D. Rubis2
JSC Molecular Electronics Research Institute, Moscow, Zelenograd, Russia ^National Research University of Electronic Technology, Moscow, Zelenograd, Russia
imatvuEhkin@niirneju
Keywords: cellular automata, permutation, random numbers, cryptography, metric Received: 10.07.2020, Revised: 24.07.2020
* The work was supported by Russian Foundation for Basic Research project no 20-07-00409
The work is devoted to the design and study of cellular automata (CA) for permutation (mixing of elements) of matrices. The aim of the study is to create the simplest CA-algorithm with the best mixing properties, which is important for cryptography and random number generation. CA formulations are proposed for the implementation of four cellular automata of matrix permutation. The algorithms are based on cyclical shifts of matrix elements along its rows and columns; local transition functions and cell state components are given (control is given by no more than three flags with three states). The calculations were carried out for square matrices and pairwise different elements. The success of the developed spacecraft was determined by the number of iterations before repeating the original table of numbers or by the period, as well as by the set of generated permutations.
To estimate the period, a numerical calculation was carried out, on the basis of which empirical formulas for the periods of CA-algorithms were obtained, depending on the order of the matrix. The concept of "individual trajectory of an element" (IT) is introduced - the path of an element on an expanded spacecraft field until it returns to its original position. The set of lengths of unique IT (DIT) together give the period according to the law of least common multiple; part of IT is visualized. For the simplest 1: 1 algorithm, the period is determined by the Landau function. For the other two schemes "1: 2" and "1 (0) 1" the period depends on the residue class of the order of the matrix modulo 12. For even orders it is linear, and for odd orders it is either quadratic or cubic (with different coefficients) polynomial dependence. For the "sin" scheme, similar residue classes are distinguished with a remainder of 1 and 3 modulo 4, where the period empirically obeys the recurrent law, but we failed to obtain the order explicitly.
Mathematics and Mathematical Modeling, 2020. no. 04. pp. 1 - 51.
DOI: 10.24108/mathm.0420.0000223
© I.V. Matyushkin, P.D. Rubis, 2020
Four Cellular Automata Algorithms Permutation
To assess the dynamics of the spacecraft, the distribution of permutations during the period of operation was calculated by the numerical values of the two introduced distance metrics relative to the original matrix. The algorithms "scheme 1: 2" and "scheme" 1 (0) 1 "are characterized by quasi-periodicity of the values of the mixing metrics and a distinct discreteness of the spectrum (the presence of gaps). The work presents and visualizes the results of a computational experiment in sufficient detail.
However, the evidence is known only for the period of the "1: 1" scheme. We also leave the analysis of the behavior of the mixing metrics for the future. The fundamental question also remains open: what is the largest number of permutations of a square matrix that a deterministic cellular automaton can generate before it stops, under certain constraints (synchronicity, closed boundaries, radius of a neighborhood 1 or 2)?
References
1. Wolfram S. Random Sequence Generation by Cellular Automata // Advances in Applied Mathematics, 1986, vol. 7, Pp. 429 - 432.
2. Sukhinin B. M. Development of pseudorandom binary sequences generators based on cellular automata Mashinostroenie i komp'yuternye tekhnologii. 2010;(9). Available at: https://cyberleninka.ru/article/n/razrabotka-generatorov-psevdosluchaynyh-dvoichnyh-posledovatelnostey-na-osnove-kletochnyh-avtomatov , accessed: 25.08.2020. (in Russian) (English translation: Mechanical engineering and computer technology).
3. Yarmolik V.N., Murashko I. A. Implementation of the pseudorandom sequence generator on cellular automata. Avtomatika i vychislitel'naya tekhnika., 1993;(3):9-13. (in Russian) (English translation: Automation and computer facilities).
4. Bar dell P. Analysis of cellular automata used as pseudorandom pattern generators // Proceeings., International Test Conference. 1990. Pp. 762-768.
5. Girau B., Vlassopoulos N. Evolution of 2-dimensional cellular automata as pseudo-random number generators // Cellular automata: 10th Intern. conf. on cellular automata for research and industry: ACRI 2012 (Santorini island, Greece, Sept. 24-27, 2012): Proc. B.; Hdbl.: Springer, 2012. Pp. 611-622.
6. Mukhamedjanov D.D. , Levina A.B. Pseudorandom number generator on cellular automata. Nauchno-tekhnicheskij vestnik informacionnyh tekhnologij, mekhaniki i optiki. 2018;(5). Available at: https://cyberleninka.ru/article/n/generator-psevdosluchaynyh-chisel-na-osnove-kletochnyh-avtomatov accessed: 25.08.2020. (in Russian) (Scientific and technical Bulletin of information technologies, mechanics and optics).
7. Bilan S., Bilan M., Bilan S. Novel pseudo-random sequence of numbers generator based cellular automata. // Information Technology and Security. - 2015. - Vol. 3, Iss. 1 (4). - Pp. 38-50. DOI: doi.org/10.20535/2411-1031.2015.3.1.57710
8. B. Martin, Patrick Solé. Pseudo-random Sequences Generated by Cellular Automata. International Conference on Relations, Orders and Graphs: Interactions with Computer Science, May 2008, Mahdia, Tunisia, Pp. 401-410.
9. Agrawal V, Bushnell M. Essentials of Electronic Testing for Digital, Memory, and Mixed-Signal VLSI Circuits. Springer, 2000. Pp. 712.
10. Hortensius P.D. Parallel random number generation for VLSI systems using cellular automata // IEEE Transactions on Computers. 1989. Vol. 28 (10). Pp. 1466-1473.
11. Fog A. Pseudo-Random Number Generators for Vector Processors and Multicore Processors. // Journal of modern applied statistical methods. 2015. №14(1) Pp 308 - 334.
12. Klyucharev P.G. Random Graph Construction for Cryptographic Applications. Matematika i matematiceskoe modelirovanie. 2017;(3):77-90. DOI: 10.24108/mathm.0317.0000076. (In Russ.) (English translation: Mathematics and Mathematical Modeling).
13. Konahovich G.F. , Puzyrenko A.U. Komp'yuternaya steganografiya Teoriya i praktika. [Theory and practice of computer steganography]. Kiev MK-Press Publ., 2006. Pp. 92-97. (in Russian).
14. Chen Y. L., Lambooij E., Mennink B. How to Build Pseudorandom Functions From Public Random Permutations Cryptology. Cryptology ePrint Archive. Available at https://eprint.iacr.org/2019/554.pdf , accessed 09.07.2020.
15. Trevisan L. U.C. Berkeley — CS276: Cryptography: Notes for Lecture 15. Available at https://people.eecs.berkeley.edu/~luca/cs276/lecture15.pdf , accessed 09.07.2020.
16. Rubio, C. F., L. Encinas, S. H. White, Á. Rey and G. Sánchez. The Use of Linear Hybrid Cellular Automata as Pseudo Random Bit Generators in Cryptography. // Neural Parallel & Scientific Comp. №12 (2004), Pp 175-192.
17. Zhukov D. A.O poryadke odnoj perestanovki na elementah kvadratnoj matricy [On the order of a single permutation on elements of a square matrix] Desyataya mezhdu-narodnaya nauchno-tekhnicheskoj konferenciya «Bezopasnye informacionnye tekhnologii» [10th Int. Conf. "Safety information technology"]. Moscow, 2019, pp 143-145 M.: (in Russian).
18. Weichuang Guo, Junqin Zhao, Ruisong Ye. A chaos-based pseudorandom permutation and bilateral diffusion scheme for image encryption // Intern. J. of Image, Graphics and Signal Processing. 2014. Vol. 6. No. 11. Pp. 50-61.
19. Matyushkin I.V., Kozhevnikov V.S. A Cellular Automata Algorithms for Matrix Permutations. Trudy MFTI. 2019. Vol. 11. no.1. pp. 39-52. (In Russ.) (English translation: proceedings of MIPT).
20. Matyushkin I.V., Kozhevnikov V.S. A Cellular Automata Algorithm for Matrix Permutations. Matematika i matematiceskoe modelirovanie. 2019;(3):1-24. (In Russ.) (English translation: Mathematics and Mathematical Modeling).
21. Deza E., Deza M.M. Dictionary of distances. Amst.: Elsevier, 2006. 391 p. (Russ. ed.: Deza E., Deza M.M. Entsiklopedicheskij slovar' rasstoianij. Moscow: Nauka Publ., 2008. 444 p.).
