CC BY
6Математика и математическое моделирование. 2019. № 03. С. 1 - 24.
DOI: 10.24108/mathm.0319.0000181
Математика ft Математическое
моделирование
© И.В. Матюшкин, В.С. Кожевников
Сетевое научное издание http://mathmelpub.ru ISSN 2412-5911
УДК 519.7
Клеточно-автоматный алгоритм пермутации матриц
Матюшкин И.В.1*, Кожевников В.С.2
2
Национальный исследовательский университет «Московский институт электронной техники», Зеленоград, Москва, Россия
2Московский физико-технический институт (национальный исследовательский университет), Долгопрудный, Россия
В формализации клеточных автоматов (КА) описывается алгоритм, перемешивающий элементы матрицы на основе циклических сдвигов строк и столбцов. Он показывает интересное поведение при нечетном порядке матрицы n, при котором после пермутаций матрица претерпевает поворот на ±90° и на 180° (отражение относительно центра), а рост периода N оказывается быстрее экспоненциального. Доказано, что N равно наименьшему общему кратному всех нечетных чисел, не превосходящих 2n, т.е. N = НОК(3,5,К ,2п -1). Динамика пермутаций анализируется с помощью введенных авторами двух «метрик», отражающих степень перемешан-ности. Результаты работы могут быть использованы при генерации псевдослучайных чисел.
Ключевые слова: клеточные автоматы, пермутация, перестановка, псевдослучайные числа, криптография, метрика
Представлена в редакцию: 06.05.2019, исправлена 21.05.2019
Пермутация, или перестановка элементов, матрицы — это один из способов генерации псевдослучайного числа. Простейшими пермутациями являются поворот, отражение и транспонирование матрицы [1, 2, 3]. Если элементы матрицы булевы или являются последовательностью нулей и единиц, то вся матрица отождествляется с двоичной записью некоторого числа, причем интерпретация зависит от способа преобразования двумерного индекса в одномерный.
Актуальность задачи поиска эффективных способов псевдослучайной пермутации матрицы обусловлена тем, что эта операция находит множество применений в различных областях. Генерация тестовых последовательностей в логических блоках современных
imalyuEhkin@riiimejii
Введение
сверхбольших интегральных схем (СБИС) может требовать случайной перестановки данных внутри матричного блока, например, схемы c программируемой логикой (ПЛИС).
Помимо генерации случайных чисел перестановка матричных элементов непосредственно используется в криптографии для шифрования изображений [4, 5, 6, 7]. Информация, содержащаяся в цифровых изображениях, обычно отличается высокой избыточностью, а также наличием сильной корреляции между соседними пикселями. Данное обстоятельство предъявляет особые требования к используемому методу шифрования, удовлетворить которым позволяет комбинирование диффузии и конфузии, реализованное в виде перемешивающего преобразования [8]. При составлении стойких криптографических алгоритмов по сей день руководствуются утверждением Шеннона о том, что хорошее перемешивающее преобразование часто является результатом многократного применения произведения двух простых некоммутирующих операций [7, 8]. В данной работе в качестве таких операций могут рассматриваться сдвиги строк и столбцов матрицы.
Применение клеточных автоматов (КА) для генерации случайных последовательностей [9] и в криптографии [10, 11] вообще известно достаточно давно. Исследуемый в данной работе алгоритм можно использовать в качестве источника детерминированного хаоса, если глобальную функцию перехода соответствующего КА рассматривать как дис-кретизированное нелинейное отображение, как, например, преобразование пекаря [7], отображение «кот Арнольда» [5], сдвиг Бернулли [5] и т.п.
Цель работы заключается в разработке и исследовании работы КА-алгоритма на основе применения двух «метрик перемешанности» к проходимой КА цепочке состояний. Выполненное исследование также представляет интерес для комбинаторики и метрической теории для пермутаций двумерных объектов.
1. КА-алгоритм
Квадратная матрица п х п помещена в поле КА с замыканием границ (тороидальная решётка). Предлагается алгоритм преобразования матрицы с использованием циклических сдвигов строк и столбцов. Рассмотрим его на качественном уровне (рис. 1). Алгоритм состоит из 4 периодически повторяющихся этапов. Преобразование матрицы за один этап будем называть шагом или тактом, а за все 4 этапа — проходом.
Рис. 1. Визуализация алгоритма
4 этапа алгоритма:
1. Нечётные столбцы сдвигаются вверх, чётные — вниз.
2. Нечётные строки сдвигаются вправо, чётные — влево.
3. Нечётные столбцы сдвигаются вниз, чётные — вверх.
4. Нечётные строки сдвигаются влево, чётные — вправо.
Приведём формальное описание алгоритма в терминах КА. Используется квадратная решётка п х п с замкнутыми границами (нумерация, традиционная для индексов матрицы). Окрестность клетки типа фон Неймана радиуса единица. Состояние клетки определяется двумя компонентами (/, ^ , см. табл. 1.
Таблица 1. Описание компонент клетки (/, ^
Название Флаг перемещения Регистр данных
Обозначение /
Множество значений {Т, 1, = {0,1,2,3} Произвольное
Описание Отвечает за перемещение элементов матрицы. Стрелки в основном указывают направление движения. Хранит элементы матрицы м=| |m||
Начальное состояние автомата для матрицы (рис. 2):
Vi, j е 1, n : s.. = m.
7 -У 7 11 11
j mod 2 = 1 ^ f = t j mod2 = 0 л i mod2 = 1 ^ f = i 7'mod2 = 0A7'mod2 = 0=>/>. = <Рис. 2. Начальное состояние флага f при n = 5
Для формулировки локальной функции перехода (ЛФП) введем понятие согласованности. Две клетки называются согласованными, если чётности их флагов совпадают. Число согласованности cons клетки определяется как число согласованных с ней соседей, т.е.
cons = |{ ¡л| f = f (mod 2)j|, где ¡л — локальный индекс клетки. Его значение, как и f,
можно указать стрелкой. Например, «t» — сосед сверху, « ^ » — сосед справа, а f —
значение флага клетки снизу. Саму центральную клетку индексировать не будем. Знак «:= » имеет смысл присвоения значения, т.е. операнд справа от него берётся до перехода, а слева — после перехода. Операция mod возвращает остаток от целочисленного деления (например, 9mod4 = 1, - 1mod4 = 3 )
Формальные правила перехода приведены ниже, в табл. 2. Нетрудно видеть, что по флагу f данный КА является циклическим. Блок правил 1 отвечает за перемещение элементов матрицы, а блок 2 — за изменение флага f .
№ Условие перехода Формула перехода
1 {(cons = 0) л (f = 0)) v {(cons Ф 0) л (f = 3)) s := ^
{(cons = 0) л (f = 1)) v (cons Ф 0) л (f = 0) s :=
{(cons = 0) л (f = 2)) v {(cons Ф 0) л (f = 1)) s :=
{(cons = 0) л (f = 3)) v {(cons Ф 0) л (f = 2)) s :=
2 (cons = 0) v (cons = 4) f :={ f +1) mod 4
-.((cons = 0) v (cons = 4)) f :={ f -1) mod 4
2. Результат работы алгоритма
Эмпирически установлено, что через некоторое число проходов алгоритм переводит матрицу в исходное состояние. Это позволяет назвать периодом алгоритма N минимальное (ненулевое) число проходов, через которое начальная конфигурация, содержащая попарно различные элементы, повторяется. Ясно, что такое определение не зависит от выбранной с этим условием начальной конфигурации.
Примечательным является период для матриц нечётных порядков. Для чётных порядков период линеен и равен п /2, а для нечётных имеет скорость роста, превышающую экспоненциальную. В табл. 3 представлены периоды для некоторых порядков п.
Таблица 3. Зависимость периода алгоритма N(n) от порядка матрицы
Порядок матрицы n Период N Каноническое разложение N
3 15 3-5
5 315 3-3-5-7
7 45045 з^7аыз
9 765765 3-3-5-7-11-13-17
11 14549535 з^7аыза7а9
Вычисление последней точки (п = 11) занимает около двух часов на обычном компьютере. Следующая точка (п = 13) потребовала бы около 200 часов, т.е. необходимы параллельные вычисления на суперкомпьютере. Несмотря на быстрый рост периода, число пройденных КА перестановок ничтожно мало по сравнению с их общим числом п2!. Уже при п = 3 просматривается лишь 1,6 -10-4 от общего числа. Несмотря на короткую длину полученного ряда (табл. 3), обратим внимание, что нечётные простые множители р присутствуют все и входят в каноническое разложение N в первой степени. Первоначально
3
мы предположили, что N = — ^ р, п > 3, но затем уточнили эту гипотезу: N есть наи-
2 р<2п
меньшее общее кратное чисел 3,5,К , 2п -1.
Рассмотрим теперь траектории отдельных элементов матрицы (рис. 3). Для этого поле с замкнутыми границами удобно представить не как тороидальную решётку, а как бесконечную плоскую периодическую решётку. Такое представление будем называть развёр-
нутым. Каждая из траекторий напоминает путь бильярдного шара, а граница поля играет роль стенки. Траектории коротко-периодичны, и их период зависит от начального положения элемента (в табл. 4 представлена эта зависимость для матрицы 13x13 ). Соответственно период работы алгоритма равен произведению индивидуальных периодов без учета разложения составных чисел на простые. Это в целом объясняет данные табл. 3. Обратим внимание, что возможны самопересечения траектории, когда позиция элемента совпадает с начальной, но направление его движения отлично от исходного.
12345123451234512
г
г1 п
ч - \ .—Ш
.1
г- Г
1 * | г 7
г —1 *— и
г _I 1 _ и Ц-
(_ 1 * и I
Ч ! . Г" —1
Т| Тгг
I
Рис. 3. Траектории некоторых элементов матрицы 5 х 5 в развёрнутом представлении. Исходная матрица — центральный квадрат. Точками отмечены начальные положения элементов, а для периферийных квадратов — места самопересечений траектории. Рядом с каждой точкой указано количество шагов до её достижения
элементом
Таблица 4. Периоды индивидуальных траекторий. В клетки таблицы записаны номера шагов, на которых происходит возвращение соответствующих элементов матрицы 13x13 в исходную позицию
26/100 20 84 36 68 52 27/52 68 36 84 20 100 4
12 70/100 28 92 44 76 35/60 60 76 44 92 28 100
92 36 34/100 52 84 68 43/68 84 52 100 36 92 20
28 84 44 62/100 60 92 51/76 76 92 60 100 44 84
76 52 92 68 42/100 84 59/84 100 68 92 52 76 36
44 68 60 84 76 54/100 67/92 92 100 76 84 60 68
33/60 41/68 49/76 57/84 65/92 73/100 25/50/ 75/100 25/92 25/84 25/76 25/68 25/60 25/52
60 52 76 68 92 84 27/100 46/100 84 92 68 76 52
44 84 60 100 76 92 27/92 76 58/100 60 84 44 68
76 36 92 52 100 68 27/84 84 68 38/100 52 92 36
28 100 44 92 60 76 27/76 60 92 44 66/100 28 84
92 20 100 36 84 52 27/68 68 52 84 36 30/100 20
12 92 28 76 44 60 27/60 44 76 28 92 12 74/100
Последний номер равен индивидуальному периоду и дает следующие значения (в проходах), считая по строкам и без повторов: 25, 5, 21, 9, 17, 13, 1, 3, 7, 23, 11, 19, 15. Таким образом, ожидается N(13) = 32 • 52 • 7 -11-13 -17 -19 • 23 .
Обоснование формулы N = НОК(3,5,К ,2п -1) для произвольного нечётного п > 3 проводится по следующему плану:
1. До пересечения границы поля элементы двигаются параллельно главной и побочной диагоналям поля, а точнее, лесенкой по четырём схемам: вверх-влево (^Т),
вверх-вправо (Т^), вниз-влево , вниз-вправо . Вертикальное направление
зависит от четности номера начального столбца ], а горизонтальное — начальной строки г.
Таблица 5. Схема перемещения элемента в зависимости от начальных номеров строки г и столбца ]
1—^^^ нечётное чётное
нечётное
чётное Т^
2. При пересечении горизонтальной границы поля меняется направление движения по горизонтали (например, Т^ изменится на ^Т), при пересечении вертикальной границы — по вертикали ( Т^ изменится на ^^ ).
3. Период отдельного элемента с начальными индексами (г, j) явно выражается через
положение, занимаемое данным элементом сразу после первого пересечения границы, и направление движения в этот момент.
4. Период N нечётный и меньше 2п.
5. Для каждого нечётного числа от 3 до 2п — 1 существует элемент с таким периодом.
6. N(п) = НОК N (п)}; 1 = НОК(3,5,К ,2п — 1) . Последовательно докажем каждый пункт.
1. Из определения алгоритма следует, что вертикальные и горизонтальные перемещения каждого элемента чередуются. Покажем, что вертикальное и горизонтальное направления перемещения сохраняются (до пересечения границы). Иначе говоря, если элемент переместился вверх, то через шаг он опять переместится вверх; если он переместился влево, то через шаг опять переместится влево, и т.д. В самом деле, так как направление движения столбца, в котором изначально находился элемент, через шаг изменится (рис. 1 ), то в соседних столбцах направление будет совпадать с исходным. Но в результате горизонтального перемещения на втором шаге рассматриваемый элемент и окажется в одном из соседних столбцов. Аналогично доказывается сохранение горизонтального направления. Это можно выразить формулами (1):
|/(г + 2) = /(г) + А/ 1/(г+2) = /(г)+А/
где /(г) — строка, в которой находится элемент на шаге г ; / (г) — столбец; А/, Л/ = ±1 — направления вертикального и горизонтального движения. Отсюда, с учётом очерёдности движения строк и столбцов, следуют выражения (2):
г +1"
/(г) = /о + А/о •
/(г) = /о + А/о
(2)
где (/0, к) — начальные индексы элемента; А/0, А/ — направления движения до пересечения границы; [х] — целая часть числа, х -1 < [х] < х .
Как следует из формул (2), существует шаг ^ , на котором / и / либо превышают п ,
либо становятся меньше 1, что соответствует пересечению границы поля (переходу в соседний квадрат в расширенном представлении). Так как после пересечения соотношения (1) остаются справедливыми, то при г > г1 имеют место формулы, аналогичные (2). В процессе движения может происходить множество пересечений, поэтому обозначим шаг, на котором происходит к -е пересечение, через гк . Тогда при гк < г < движение описывается формулами (3):
/(г) = /к + А/к • /(г) = /к + А/к
Г " г +1" Г ч+1" Л
V _ 2 _ _ 2 _ У
г " г' " к" Л
V _ 2 _ _ 2 _ У
, Ч < г < гк+1
(3)
где (/к, /к) — индексы элемента сразу после к -го пересечения границы, а А/,, А/к — новые направления движения. Положим по определению г0 = 0, тогда формулы (3) будут
верны при всех г > 0, обобщая формулы (2).
Из определения алгоритма следует, что на первом шаге каждый элемент переместится вертикально, причём, в случае нечётного номера столбца /0 — вверх (А/0 = -1), а в случае чётного — вниз (А/0 = 1). При этом, если номер строки / нечётный, то элемент окажется в чётной строке, а значит на втором шаге переместится влево ( А/0 = -1). Соответственно, если номер / чётный, то элемент переместится вправо ( А/о = 1). Итак,
А/о = (-1); А/о = (-1У
(4)
2. Здесь принципиальна нечётность порядка матрицы п , из которой следует, что граничные столбцы и строки нечётные, а следовательно, имеют одинаковое направление движения. Пусть на шаге ^ = ^ -1 элемент горизонтально двигался в строке /(¡к -1) е {1, п} , а на шаге I = ^ элемент пересекает горизонтальную границу поля, т.е. переходит в строку /к = п +1 -/(¡к -1) . Тогда направление горизонтального перемещения строки / (¡к -1) на шаге ^ = ^ +1 будет противоположным направлению на шаге ^ = ^ -1, но оно совпадает с направлением движения строки /к, так как и /(¡к -1), и /к нечётные. Аналогично рассматривается случай пересечения вертикальной границы поля. Итак, если была пересечена горизонтальная граница, то 41к = Д/^ , 4]к = -44-1, если же вертикальная, то 41к = ,
44 = 44-1 •
3. Мы показали, что до пересечения границы поля движение по горизонтали и по вертикали происходит в одном направлении, т.е. номера текущей строки и столбца монотонно изменяются, Сразу сделаем замечание, что формулы (2) имеют смысл лишь при 1 < /(¡), ¡(?) < п , а также в случае, когда один из индексов принимает значение 0 или п +1, что отвечает пересечению границы поля.
Рассмотрим движение элемента после пересечения горизонтальной границы на некотором шаге ^ . Сразу можно сказать, что ^ нечётно, так как горизонтальная граница может быть пересечена только при вертикальном перемещении, которое всегда происходит на нечётном шаге. Номер строки после пересечения может иметь два значения: /к е {1,п} .
4 = 1 при движении вниз, т.е. при 4^ = = 1, и ^ = п при движении вверх, т.е. при 4^ = 44_! = -1. Это можно выразить формулой
гк =1-п (1 + 4/к) + п, (5)
которая пригодится в дальнейшем. Покажем, что следующая граница поля, которую элемент пересечёт на шаге , будет вертикальная. чётно, так как горизонтальные перемещения происходят только на чётных шагах. Тогда при I = формулы (3) принимают вид
'с, - и -1 ^
/(¡к+1) = /к +4/к-I ч+' * I
> 2 Л (6)
Ж+1) = ¡к + 4к ■[ 4+1 +1 ]
Если на шаге ¡к+1 пересекается горизонтальная граница, то /(¡к+ ^ = п +1 при 4^ = 1 и / 1) = 0 при 44 = -1, откуда находим, что ¡к+1 - ^ = 2п +1. Если же пересекается верти-
кальная граница, то /(1к+^ = п +1 при А/ = 1 и /= 0 при А/ =-1. Можно написать,
ч 1+ А/, , что /(гк+1) = •(п +1), откуда
- = 2ЛЫ-А +, = п + А/,(п +1 - 2/) = (2" +1 - А ДА/" = 1 (7)
А/к [2Л - 1 А/ =-1
Каким бы ни было А/ = ±1, имеет место неравенство п + А4(" +1 - 2/) < 2п +1, т.е.
вертикальное пересечение произойдёт раньше.
Подставляя выражения (7) и (5) в формулу (6), можно найти строку пересечения на
шаге гк+1:
, = ,(/ ) = п +1+ А/,А/,(п +1 -2/) =(п +1 -Л, Ч =АЛ (8)
гк+1 = /(Ч+1) =---= \ . к- ^ к- (8)
Номер столбца /+1 = 1 при движении вправо, т.е. при А/ = 1, и /+1 = п при движении влево, т.е. при А/ = -1, что задаётся формулой (9):
1 - п
/-+1 = — (1 + / + п. (9)
Как уже было установлено, при пересечении вертикальной границы меняет знак коэффициент А,, следовательно
(А/\,, = -А/,
( -+1 к (10) 1АЛ+1 = АЛ
Итак, мы получили важные характеристики движения элемента после пересечения горизонтальной границы, такие как время (число шагов) до следующего пересечения границы, причём обязательно вертикальной, и индексы после пересечения.
Аналогично рассматривая пересечение вертикальной границы на некотором шаге ^, можно установить, что следующая пересечённая граница будет горизонтальной, найти время до пересечения - ^ и столбец пересечения /+1. Приведём только результаты:
(2п +1 - 21, А/, = 1 *к+1 - ч = п + А/, (п +1 - 2/,) = «1 к к (11)
12/к -1 Ч =-1
1 - п,л к • ч
Ч+1 = —(1 + А/к) + п
. _ п +1 + А/'кА/к (п +1 - 2/к ) _ (п +1 - Ч, А/к = А/к
ук+1
(12)
2 \Ч, А/к ^АЛ
АЧ+1 = А/к „АЛ+1 = -АЛ
(13)
<
Теперь, располагая формулами (7-13), мы имеем возможность проследить всю траекторию элемента с произвольными начальными индексами (/0, j) . Покажем, что после четырёх пересечений границы он вернётся в исходное положение. Для этого рассмотрим движение из клетки (i, j ), т.е. после первого пересечения границы.
Если первая пересечённая граница горизонтальная, то из формул (7-13) следует: а)
t2 -1 = n + Aj(n +1 -2j)
Í2 = 1 (n +1 + AiiAi(n +1 - 2ji))
1-n
j2 (1 + Aji) + n Aj2 = Aj
A/2 = -Ai
б)
t3 - t2 = n + Ai2 (n +1 - 2/2) = n - А/ (n +1 - (n +1) - Ai Aj (n +1 - 2 j )) =
= n + Aj(n +1 -2j) = t2 -t
1 - n „ A. ч 1 - n ^ ... /3 =-(1 + A/2) + n =-(1 - A/j) + n
j = - (n +1 + Ai2 Aj (n +1 - 2/2)) =
2
1 (n +1 - AiAj (n +1 - (n +1) - AiAj (n +1 - 2 j ))) :
= -(n +1 -A/Aj (n +1 -(n +1)-a/Aj,
= 1(n +1 + n +1 - 2j) = n +1 - j
Ai = Ai2 = -A/ A/3 = A/2 = -Aj
Так как Аз = -Ai, то при t e [t, t2 ) и при t e[t3, tA) горизонтальное движение происходит в разных направлениях, следовательно на шагах t2 и t4 пересекаются разные вертикальные границы (левая и правая).
в)
t4 - t3 = n + А3 (n + 1 - 2 j ) = n - Aj (n + 1 - 2(n +1 ) + 2 j ) = n + Aj (n + 1 - 2 j ) = t2 -1 /4 =1 (n +1 + А/'зА/з (n +1 - 2 j )) =1 (n +1 + AilAjl (n +1 - 2(n +1) + 2 j )) =
= (n +1)^ + A/1Aj1 • j
1 -n ,л А.ч 1 - ...
j4 = (1+А/з)+n = —(1 -АЛ)+n Ai4 = -Ai = А/ A/4 = А/з = -A/1
Так как А/4 = А/ = - А/2, то при I е [¿2, 1Ъ) и при I е [¿4, ) вертикальное движение происходит в разных направлениях, следовательно на шагах ¿3 и ^ пересекаются разные горизонтальные границы (верхняя и нижняя). Таким образом, при движении все четыре границы поля пересекаются поочерёдно.
г)
Ч - Ч = п + Д*4 (п +1 - 2/4 ) = п + А/ (и +1 - (п +1)(1 -А/'; А/) - 2А/ Ад • / ) = = п + А/1 ((п +1) А/1 А/1 - 2/ • А/1А/1) = и + А/1 (и +1 - 2/ ) = Х2 - ^
1 - п ч 1 - п ...
/5 = (1+ А/4) + п = (1+ ) + п
/ = 1 (п + 1 + А/4А/4 (п +1 - 2/4 )) =
1 (n +1 - A/jA/j (n +1 - (n +1)(1 - A/ A/) - 2A/1A/J • j )) = 1 (n +1 - A/jA/1 (A/jA/j (n +1 - 2jj))) =1 (n +1 - (n +1) + 2jj) = jj
= — (n +1 - A/ A/ (A/Aj (n +1 - 2/))) = - ( n +1 - (n +1) + 1 - n
По формуле (5) /5 = —(1 + A/1) + n = /1.
Если первой пересекается вертикальная граница, то рассуждения полностью аналогичны. В общем случае, для всех к е 1,4 верно равенство:
Ч+1 - tk = n + a(n +1 - 2/), где а = A/, / = j, если на шаге ^ была пересечена горизонтальная граница, и а = A/, / = / , если вертикальная. При этом /5 = / и j5 = j, т.е. элемент возвращается в положение (/, j) на шаге t = ^ + т, где т = 4(n + a(n +1 - 2/)).
Покажем, что т и есть период данного элемента. Для этого заметим, что алгоритм является обратимым, т.е. предыдущее состояние однозначно определяется текущим состоянием и номером шага по модулю 4, точнее числом (t mod 4), что следует из определения алгоритма (4 этапа и цикличность). Так как Tmod4 = 0, то элемент возвращается в клетку (/, j ) на том же этапе алгоритма, что и входит в неё на шаге ^ . Из обратимости алгоритма следует, что на шаге t = (t + т) -1 = т элемент находился в той же клетке, из которой он дошёл до (/, j ) за ^ шагов, т.е. в (/0, j0) . Итак, элемент действительно возвращается в (/0,j0) за т шагов, т.е. за т/4 проходов.
Докажем, что на промежутке времени (измеряемого в шагах) [1,т) нет шага, кратного 4, на котором элемент также проходил бы клетку (/0, j0) . В самом деле, если бы на некотором шаге тх, таком что тх mod 4 = 0, элемент оказался бы в (/0, j0), то первой границей, которую бы он пересёк, была бы та же, что и на шаге t , в силу цикличности алгоритма.
Но, так как пересекаемые границы чередуются, то пересечение этой границы возможно лишь на шагах ^ и ¿5 . На шаге ^ оно произойти не может, иначе тх е [1, ^) , а, как следует из формул (2), при этом г и ] изменяются монотонно от /0 и ]0, следовательно, не могут вернуться к исходным значения до шага ^ . Таким образом, тх е [¿4, ¿5) , но тогда тх = т, так как на [¿4, ¿5) каждая клетка проходится элементом лишь один раз, потому что г и ]
опять же меняются монотонно, что следует из формул (3). Итак, т — первый кратный 4 шаг, на котором происходит возвращение элемента в исходное положение. Для периода данного элемента можно написать формулу
^ = п + а(п + 1- 21), (14)
где а = А], I= ], если на шаге ^ была пересечена горизонтальная граница, и а = А/,
I = / , если вертикальная.
4. Нечётность N следует из формулы (14), как и тот факт, что < 2п — 1.
5. Найдём явную зависимость N от (/0, ]0), определив, как соотносятся (/, ]) и (/0, ]0). Для этого перепишем формулы (2) для шага ^ в виде
(i(t1) - О -A/q = (7(0 - 7 c) -Ac =
+ ^ mod 2
(15)
Вычитая одно равенство в (15) из другого, получаем
(i(ti) - ic) - A/q = (/'(ti) - 7'c) - A/c + ti mod 2
(16)
Далее следует рассмотреть различные случаи чётных и нечётных начальных индексов (/0, ]0). Например, если /0 mod2 = 0 и ]0 mod2 = 0, то, по формулам (4), А/0 = 1 и А]0 = 1.
Если ) = п +1 (случай невозможен, так как Ло = 1 > 0 и ](£) возрастает), т.е. первой пересекается правая граница, то ^ mod2 = 0, то из уравнения (16) следует, что
/1 = п +1 + ic - /с.
(17)
Так как 1 < / < п, то из полученной формулы (17) находим условие реализации данного случая:
]0 — г0 ^ 1
Если же /(^ ) = п +1, т.е. первой пересекается нижняя граница, то ^ mod2 = 1 и
7i = n + 7с - /с,
а так как 1 < 7г < п, то
7с - ic < с
Остальные случаи рассматриваются аналогичным образом. Результаты сведены в табл. 6.
Таблица 6. Зависимость , / ) и N ^ от ('0, /)
Условие '1 Щ/о
[/0 шоё 2 = 0 /о ^'о +1 п + 1 + 'о — /о 1 2('о — /о) + 2п + 1
!./отой2=о /о < 'о 1 п + /о — 'о 2(/о — 'о + п) — 1
р0 шod2 = 1 'о ^ /о + 1 'о — ./о п 2(/о — 'о) + 2п + 1
1./оШоа2 = 1 'о < /о п /о — 'о +1 2('о — /о) + 2п — 1
р0 шod2 = о 'о + /о ^ п + 2 'о + /о — п — 1 1 4п + 3 — 2('о + /о)
[./о шod 2 = 1 'о + /о < п + 1 п 'о + /о — 1 2('о + /о) — 3
[/0 mod 2 = 1 'о + /о < п 'о + /о п 2('о + /о) — 1
[./о шod2 = о 'о + /о ^ п + 1 1 'о + /о — п 4п + 1 — 2('о + /о)
Как следует из предыдущего пункта, любой индивидуальный период Щ нечётный , а значит представим в виде Щ. = 2р —1, р е Г . Нас интересует лишь случай 3 < ^ < 2п — 1, или р е 2, п . Подберём такие начальные индексы (/0, /), что = 2р —1 при заданном р е 2, п . Несложно проверить, что при нечётном р подойдут значения
к=1
и =п+1 - р
При этом реализуется случай = 2(/0 — /) + 2п — 1 = 2р — 1. При чётном же р можно взять значения
к=п
/о = п +1 — р
Тогда И/ = 4п + 1 — 2(/о + /0) = 2р — 1.
6. Период алгоритма кратен периоду каждого элемента. Так как по определению период алгоритма есть минимальное число проходов, за которое матрица возвращается в исходное состояние, то И(п) = НОК («)}. 1. Среди {Щ (п)| ^ есть лишь нечётные числа, не превышающие 2п — 1, причём все такие числа. Следовательно, N(п) = НОК(3,5,К ,2п — 1) . Ч.т.д.
3. Используемые «метрики перемешанности»
Каждый ход КА дает нам пермутацию матрицы относительно начальной (здесь и далее считаем все элементы различными). Для описания динамики КА желательно иметь числовую характеристику, назовем её mixedness, «удалённости» текущей перестановки от начальной конфигурации. Существует большое число разнообразных метрик, определённых на множестве пермутаций массива. Среди них l -метрики, Хэммингова метрика, лексикографическая и другие, которые, однако, все зависят от значений элементов перестановок [12]. Здесь же предлагаются «метрики», характеризующие лишь взаимное расположение элементов. Мы будем использовать две «метрики»: матричную mm (matrix mixedness), которая наследована от манхэттенской метрики на плоскости и является метрикой в математическом смысле, и линейную lm (linear mixedness), которая адресуется массиву данных. Метрика mm формально схожа с l -метрикой, однако адаптирована специально для матричных перестановок. «Метрика» lm не является математической метрикой, но более эффективна по сравнению с mm в поиске особых конфигураций.
Пусть в некоторой матрице A порядка n записаны попарно различные элементы са ,
а е 1, п2 , а В является пермутацией А . Тогда определим метрику тт на этих двух матрицах формулой
1 п2 тт = "Г Е 5а >
П а=1
где =| /2 - \ | +1 у2 - ^ |, если (\, ^ ) — это индексы элемента са в матрице А , а (/2, у2) — индексы са в матрице В . Выражение имеет смысл расстояния между положениями элемента са в матрицах А и В по метрике Манхэттена.
Элементы са можно отождествить (поставить во взаимно-однозначное соответствие) с парами чисел (к,1), к, I е1, п, так, чтобы элементам матрицы А = ||аи|| соответствовали пары их же индексов, например, через конструкцию аи = к + п •(/ -1) . После перестановки элемент ак1 окажется на позиции (/, ]) в матрице В = Щ , т.е. Щ = аи = (к,1) . Например:
а11 а12
а
13
a21 a22 a23
V a3i a32 a33 J
(1,1) (1,2) (1,3)' (2,1) (2,2) (2,3) (3,1) (3,2) (3,3)
(1,1) (2,2) (3,1)' (2,1) (3,2) (2,3) (1,3) (1,2) (3,3)
Г К b12 Ь13 '
b21 b22 Ь23
V Ь31 Ь32 Ь33 j
Тогда mm полученной матрицы вычисляется по формуле
mm =
1
n
i, j=1
У
где ^ =| к -11 +11 - j |, если Ь = (к,1) . Например, для приведённой выше матрицы третьего порядка тт = -1 (£31 + £12 + + ^32 + £13 ) = — (4 +1 +1 + 2 + 2) =10 .
п 9 9
Эмпирически выполнен поиск наибольшего значения тт при заданном п. В качестве начальной конфигурации были использованы либо исходная матрица, либо отраженная от центра. Алгоритм поиска простой:
1. Задаётся произвольная конфигурация матрицы.
2. Переставляются 2 элемента со случайными номерами.
3. Если при этом тт не возрастает, то перестановка элементов отменяется.
4. Повторяется пункт 2.
В диапазоне 1 < п < 40 с числом итераций 104 не найдено конфигураций со значением тт, превышающим п.
«Метрика» 1т вводится сначала для одномерного массива. Считая, что начальное его
состояние имеет вид ак = к, к е 1, т, 1т -метрика определяется формулой
, т +1 1т =-,
1т + Ь
где Ц — максимальная длина возрастающей подпоследовательности, а ^ — максимальная
длина убывающей подпоследовательности. Числитель выбран так, чтобы 1т исходного массива равнялась 1. В самом деле, для него Ц = т, а ^ = 1 (подпоследовательность из
одного элемента).
Чтобы применить 1т к матрице, необходимо сначала сопоставить матрице массив, например, записать её элементы в столбец. Это можно сделать различными способами, например, змейкой или спиралью. Мы будем использовать простую схему, заимствованную из МАТЬАБ. Столбец получается последовательной конкатенацией всех столбцов матрицы, взятых в порядке возрастания индекса. Итак, по определению 1т матрицы есть 1т соответствующего ей столбца.
Можно ожидать, что пермутации с экстремальными значениями «метрики» особенные.
4. Анализ динамики КА алгоритма с помощью «метрик»
Ниже приведены графики тт и 1т, соответствующие применению алгоритма к матрицам различных порядков. Точки, отвечающие особенным пермутациям, отмечены крестиком «х» и обозначены греческими буквами. Рассмотрим вначале четные порядки п шоё2 = 0. В качестве иллюстрации общего случая рассмотрим п = 20. Динамика обеих «метрик» показана на рис. 4 (с детализацией до такта (шага) г, а не прохода).
а)
»
10
о «
? 6
1 тт я= 20
¿жг а • 120:20) \ А
у < * 0 V А С в
г В В с А \
У \
й Аппроксимация: V = + 2Ж1-- 0,02 Л Шаг
10 15 20 25 30 35 ¿С 4&
й)
% V "" 20
7 Л|111рОКС11М44Ш1^ убм ьакмцей 1 ни \ у = {57,75}/^ ¡2,3*)
*
* * \
• * • \ > ч
• ■ * • •
>
Аппроксимация возрастающей ктьн; V = С-59,71Гх™ 2Я, Шаг
1л™— ■ ■ ■>
10 15 70 го 35 40
Рис. 4. Перемешанности последовательных конфигураций матрицы порядка п = 20, полученных применением алгоритма. В отличие от метрики тт (а) наблюдается асимметрия метрики 1т (б), по-видимому, связанная с упорядоченностью элементов матрицы по столбцам, а не строкам, и выбором флагов
первого такта
«Метрика» тт достигает максимума а на середине пути (t = 20), а 1т в той же точке близка к локальному минимуму. Соответствующая конфигурация получена диагональной перестановкой блоков (рис. 4а); в силу четности п исходная матрица легко представляется в блочном виде. То же наблюдается и при других чётных порядках п.
Уравнение аппроксимирующей тт параболы имеет вид, близкий к
2 2 X X
у =---ъ 2X ^ — + 2X . Оказывается, что на чётном шаге I метрика тт совпадает со
20 п
значением, полученным по этой формуле. На нечётном же шаге реальное значение тт
меньше точно на 0.05 = ——>1. Таким образом, для чётных п можно обобщить (равенст-
20 п
во проверено для п = 20 и п = 8 ):
Г г2
— + 2г, г шоё2 = 0
тт (г) = <! п
w ' г2 +1
--+ 2г, г шоё2 = 1
п
«Метрика» 1т в целом показывает нам две симметричные ветви, но эта симметрия нарушается в центре. Возможно, это связано с изначальной асимметрией движения столбцов в начальной конфигурации КА относительно порядка нумерации. Обе ветви аппрок-
симируются гиперболой у =
х + q
. Параметры аппроксимации МАТЬАБ (значение, ин-
тервал значений, доверительная вероятность): для возрастающей ветви -р = -59,71 (-62,46; -56,97), q =-28,16 (-28,74; -27,58), Я2 = 0,9942 ; для убывающей ветви р = 57,75 (55,16;60,33), q =-12,36 (-12,95;-11,76), Я2 = 0,9945.
Рассмотрим теперь случай нечетных порядков на двух примерах: п = 5 (рис. 5, 6) и п = 7 (рис. 7, 8). Дополнительно строим гистограммы «число пермутаций — значение метрики».
Рис. 5. Метрика тт последовательных конфигураций матрицы порядка п = 5, полученных применением
алгоритма
Рис. 6. «Метрика» 1т последовательных конфигураций матрицы порядка п = 5, полученных применением
алгоритма
Матрица в соответствует одновременно глобальному максимуму тт и глобальному минимуму 1т. Она представляет собой отражение исходной матрицы ||агу|| относительно
центра: Д = аи+1_г. и+1 . (элементы матрицы обозначены той же буквой, что и матрица, с индексами). Матрицы а и у, на которых 1т достигает глобального максимума, получены
вращением а.. против и по часовой стрелке соответственно: аг = а
Ля+1-г'
п = а
n+1-j
г. Од-
нако эти экстремумы тт -метрика «не замечает». Аналогичная ситуация наблюдается и при п = 7 .
Как и при п = 5, матрица в есть отражение исходной относительно центра, а а и у — повороты. Однако в данном случае а — поворот по часовой стрелке, а у — против часовой стрелки. То есть при п = 5 и п = 7 алгоритм проходит эти конфигурации в разном порядке. Если же заметить, что отражение Д есть ни что иное как поворот на 180°, то получится, что алгоритм поворачивает матрицу на 90° за число шагов га, равное периоду N . Напомним, что алгоритм четырехтактовый, следовательно, он возвращает матрицу в исходное состояние за 4N шагов. Для п = 5 вращение происходит против часовой стрелки, а для п = 7 — по часовой стрелке. Примерное объяснение этому факту состоит в аналогии с волной.
Обращает на себя внимание асимметрия гистограмм и наличие, во всяком случае, визуальное, для п = 7 и тт -метрики суперпозиции трех дискретных гауссовых распределений. Для больших порядков матрицы этот феномен требует более глубокого исследования.
Рис. 7. Гистограмма метрики тт для матрицы порядка п = 7, полученных применением алгоритма
а)
4,0 3,5 3,0 2,5
1.0 щ
1,0 0,5
а,о
б)
350СЮ 10000 15000
гоооо
15000 10000 5000 о
1щ X и — 7 у( 135135:3.57} и
«<45045; 3,57| * я
т » . - , - ... ■ ■ . -„ -
; 1 " = =......: '
1 ™ * X ■ « л " 4 Ь.1 ■ ■
^ {90090; 1. 00)
Шаг
■ >-
гсосо 40000 60000 ЭЮСОЗ 100000 13 СОПИ 140000 150000 120000
Гистограмма й«, п - 7
I
.....и 111! 1.
С} ^ ^ Н N Н И № '"г ^ ^ ^ 'ч Ч'Ч Ч ^ И ™ Т 14 т 'Л Г"" ° ^
Рис. 8. «Метрика» (а) 1т и её гистограмма (б) последовательных конфигураций матрицы порядка п = 7,
полученных применением алгоритма
Заключение
В результате анализа алгоритма был выявлен ряд примечательных его свойств. В частности, последовательность его периодов N для нечётных порядков матриц n задаётся формулой N(n) = НОК(3,5,K , 2n -1) . Величина периода косвенно говорит в пользу алгоритма как генератора случайных чисел. Более того, период может быть существенно увеличен небольшим усложнением алгоритма, например, введением различных длин для циклов смены направления столбцов и строк. Останавливая алгоритм в случайный момент времени, текущую перестановку матрицы можно рассматривать как массив случайных чисел или же преобразовать её каким-либо способом в одно случайное число.
Установлено, что в ходе работы алгоритма отдельные элементы (рисунок 3) перемещаются аналогично бильярдным шарам под углом 45°. Для усложнения динамики можно реализовать движение и под другим углом, например, arctg(2).
Введённые «метрики перемешанности» позволили определить, что в результате применения алгоритма матрица претерпевает вращение в случае нечётного порядка и блочную перестановку в случае чётного. В результате численного эксперимента получена примечательная гистограмма для n = 7 (рис. 7), форма которой имеет вид трёх нормальных распределений, наложенных друг на друга. Изменение «метрик» mm и lm на последовательных конфигурациях матрицы имеет псевдослучайный вид (рис. 5, 6).
Список литературы
1. Матюшкин И.В. Алгоритмы параллельных вычислений в формализации клеточных автоматов: сортировка строк и умножение чисел по схеме Атрубина // Проблемы разработки перспективных микро- и наноэлектронных систем (МЭС). 2016. № 4. С. 77-81.
2. Матюшкин И.В., Жемерикин А.В., Заплетина М.А. Клеточно-автоматные алгоритмы сортировки строк и умножения целых чисел по схеме Атрубина // Изв. высш. учеб. заведений. Электроника. 2016. Т. 21. № 6. С. 557-565.
3. Матюшкин И.В., Заплетина М.А. Отражение и транспонирование данных в матрице кле-точно-автоматного вычислителя // Изв. высш. учеб. заведений. Электроника. 2019. Т. 24. № 1. С. 51-63. DOI: 10.24151/1561-5405-2019-24-1-51-63
4. Ji Won Yoon, Hyoungshick Kim. An image encryption scheme with a pseudorandom permutation based on chaotic maps // Communications in Nonlinear Science and Numerical Simulation. 2010. Vol. 15. No. 12. Pp. 3998-4006. DOI: 10.1016/j.cnsns.2010.01.041
5. 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. DOI: 10.5815/ijigsp.2014.11.07
6. Kaur M., Kumar V. Efficient image encryption method based on improved Lorenz chaotic system // Electronics Letters. 2018. Vol. 54. No. 9. Pp. 562-564. DOI: 10.1049/el.2017.4426
7. Fridrich J. Symmetric ciphers based on two-dimensional chaotic maps // Intern. J. of Bifurcation and Chaos in Applied Sciences and Engineering. 1998. Vol. 8. No. 6. Pp. 1259-1284.
DOI: 10.1142/S021812749800098X
8. Vinod Patidar, Pareek N.K., Purohit G., Sud K.K. A robust and secure chaotic standard map based pseudorandom permutation-substitution scheme for image encryption // Optics Communications. 2011. Vol. 284. No. 19. Pp. 4331-4339. DOI: 10.1016/j.optcom.2011.05.028
9. 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. DOI: 10.1007/978-3-642-33350-7_63
10. Souyah A., Faraoun K.M. An image encryption scheme combining chaos-memory cellular automata and weighted histogram // Nonlinear Dynamics. 2016. Vol. 86. No. 1. Pp. 639-653. DOI: 10.1007/s11071-016-2912-0
11. Ключарёв П.Г. Построение случайных графов, предназначенных для применения в криптографических алгоритмах, основанных на обобщенных клеточных автоматах // Математика и математическое моделирование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2017. № 3. С. 77-90. DOI: 10.24108/mathm.0317.0000076
12. Деза Е.И., Деза М.М. Энциклопедический словарь расстояний: пер. с англ. М.: Наука, 2008. 444 с. [Deza E., Deza M M. Dictionary of distances. Amst.: Elsevier, 2006. 391 p.].
Moh^ a_ndMathematicalModelmg 2019 Mathematics & Mathematical
DOI: 10.24108/mathm.0319.0000181
[¿on ®
Electronic journal
© I-V. Matyushkin, V.S. Kozhevnikov http://mathmelpub.ru ISSN 2412-5911
A Cellular Automata Algorithm for Matrix Permutations
I.V. Matyushkin1*, V.S. Kozhevnikov2
1National Research University of Electronic Technology, Zelenograd, Moscow, Russia 2Moscow Institute of Physics and Technology, Dolgoprudny, Russia
imalyuEhkin@riiimejii
Keywords: cellular automata, permutation, rearrangement, pseudorandom numbers, cryptography, metrics Received: 06.05.2019, Revised: 21.05.2019
The article describes an algorithm to provide permutation of matrix elements through cyclic shifts of rows and columns and gives a formal description of a cellular automaton (CA) that implements this algorithm. For this, a square n x n lattice with closed boundaries and a neighbourhood of a von Neumann type cell are used.
As a result of a computational experiment for the initial orders of the matrix n, it was found that after a sufficiently large number of steps, the algorithm transfers the matrix to the original one, i.e. has period N. For odd orders of the matrix, the growth of N as a function of n is faster than exponential.
A movement of the individual elements of the matrix has been analysed to show that they move in a manner similar to billiard balls. The element moves to the matrix boundaries at an angle of 45 ° and changes direction when it reaches the bound. A defined explicit dependence of the element movement period on its initial position in the matrix allows us to prove that the global period N is equal to the least common multiple of all the odd numbers being less than 2n, i.e. N = LCM(3,5, ..., 2n-1).
To analyse the permutation dynamics, two authors-introduced "metrics" that reflect a degree of randomizing were used. One of the metrics is introduced specifically for the matrix, the other for the linear array and depends on how the matrix is transformed into a one-dimensional array. By searching for permutations with extreme metric values, it was found that as a result of permutations, even-order matrices undergo block permutation. In the case of an odd order, the matrix undergoes a rotation of ± 90 ° and 180 ° (reflection relative to the centre). Moreover, the direction of rotation depends on the order n. For example, for n = 5, counter-clockwise rotation occurs, and for n = 7, it is clockwise.
The algorithm can be used to generate pseudorandom numbers, and a value of its period indirectly argues for it. The period can be significantly increased through a slight complication of the algorithm, for example, by introducing different lengths for cycles of changing a direction of columns and rows. Stopping the algorithm at random time, one can consider the current permutation of the matrix as an array of random numbers or somehow transform it into one random number.
References
1. Matyushkin I.V. Algorithms of the parallel computations in the formalization of cellular automata: the sorting of strings and the multiplication of numbers by Atrubin's method. Problemy razrabotki perspektivnykh mikro- i nanoelektronnykh system (MES) [Problems of development of perspective micro-and nanoelectronic systems (MES], 2016, no. 4, pp. 77-81 (in Russian).
2. Matyushkin I.V., Zhemerikin A.V., Zapletina M.A. Cellular automaton algorithms for string sorting and integer multiplication by Atrubin's scheme. Izvestiia vysshikh uchebnykh zavedenij. Elektronika [J. "Proc. of Universities. Electronics"], 2016, vol. 21, no. 6, pp. 557-565 (in Russian).
3. Matyushkin I.V., Zapletina M.A. Data reflection and transposition in the matrix of a cellular automata computator. Izvestiia vysshikh uchebnykh zavedenij. Elektronika [J. "Proc. of Universities. Electronics"], 2019, vol. 24, no. 1, pp. 51-63. DOI: 10.24151/1561-5405-2019-24-1-51-63 (in Russian)
4. Ji Won Yoon, Hyoungshick Kim. An image encryption scheme with a pseudorandom permutation based on chaotic maps. Communications in Nonlinear Science and Numerical Simulation, 2010, vol. 15, no. 12, pp. 3998-4006. DOI: 10.1016/j.cnsns.2010.01.041
5. 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. DOI: 10.5815/ijigsp.2014.11.07
6. Kaur M., Kumar V. Efficient image encryption method based on improved Lorenz chaotic system. Electronics Letters, 2018, vol. 54, no. 9, pp. 562-564. DOI: 10.1049/el.2017.4426
7. Fridrich J. Symmetric ciphers based on two-dimensional chaotic maps. Intern. J. of Bifurcation and Chaos in Applied Sciences and Engineering, 1998, vol. 8, no. 6, pp. 1259-1284.
DOI: 10.1142/S021812749800098X
8. Vinod Patidar, Pareek N.K., Purohit G., Sud K.K. A robust and secure chaotic standard map based pseudorandom permutation-substitution scheme for image encryption. Optics Communications, 2011, vol. 284, no. 19, pp. 4331-4339. DOI: 10.1016/j.optcom.2011.05.028
9. 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. DOI: 10.1007/978-3-642-33350-7_63
10. Souyah A., Faraoun K.M. An image encryption scheme combining chaos-memory cellular automata and weighted histogram. Nonlinear Dynamics, 2016, vol. 86, no. 1, pp. 639-653. DOI: 10.1007/s11071-016-2912-0
11. Klyucharev P.G. Random graph construction for cryptographic applications. Matematika i matematicheskoe modelirovanie [Mathematics and Mathematical Modelling], 2017, no. 3, pp. 77-90. DOI: 10.24108/mathm.0317.0000076 (in Russian)
12. 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.).
