Компьютерный практикум (магистратура)
Черновики лекций, весенний семестр 2021
Весенний семестр
Цель -- овладение системой компьютерной математики
SageMath
Volfram Mathematika, MathLab, Maple, Maxima, Octava, ...
Почему SageMath?
1. Свободная бесплатная система.
Свободные системы развиваются быстрее коммерческих,
т.к. их развивает весь мир.
SageMath включает в себя практически все существующие
свободные математические пакеты.
2. SageMath -- надстройка над языком Python (с января 2020 --
над Python3). Значит, не надо учить еще один язык программирования
для тех, кто знает Python.
Начнем с изучения языка Python3. Нам интересно решать математические задачи.
1. Задачи, связанные с теорией чисел и криптографией
(кодирование с открытым ключом и т.п.)
2. Классы в Python.
3. Графические оконные приложения (с использованием модуля tkinter)
TCL/Tk.
4. Сетевые приложения на Python, параллельные приложения
(нити, объекты синхронизации)
5???. Нейросети?
SageMath
5. Графика треугольника.
6. Кубическая интерполяция, построение кубических сплайнов.
7. Системы алгебраических уравнений, базисы Грёбнера (обобщение
метода Гаусса исключения переменных в случае систем
линейных уравнений). Элементы коммутативной алгебры.
Язык программирования Python3
=============================
Гвидо Россум -- создал язык программирования, очень удобный
для него самого. Питонисты -- программисты, которые предпочитают
Питон всем другим языкам.
Python2 --> Python3 (современная версия).
Объектно ориентированный, императивный (не функциональный) язык
с динамической типизацией.
Динамическая типизация означает, что тип выражения невозможно
определить по тексту программы, он определяется только на Run Time
(при выполнении). Переменные в Питоне не имеют типов и содержат
объектные ссылки (handle) на объекты, тип определяется по объекту.
Объекты можно создавать, но нельзя самому удалять (в Питоне, если
очень хочется, то можно delete), ненужные (беспризорные) объекты
могут быть удалены сборщиком мусора GC (Garbage Collector).
Python очень популярен в тематике, связанной с нейросетями
(искусственный интеллект) Tensor Flow (Питон использует модули,
написанные на C++).
Для Python есть интерпретатор, который позволяет использовать
компьютер как калькулятор. Визуально отсутствует процесс компиляции.
Интерпретатор также сильно облегчает процесс изучения языка
(можно всё быстро попробовать).
Python интересен и тем, что включает элементы функционального стиля
программирования (списки, list comprehension и др.).
Пример (hello, World!).
Напечатать квадраты 20 первых натуральных чисел
В стиле традиционных языков
Описаний переменных в Питоне нет!
int x; double y; double a[100]; // Такого нет!
n = 1
while n <= 20:
print(n, n*n)
n += 1
print("Bye!")
Списки (на самом деле, список в Питоне -- это динамический массив,
аналог std::vector в С++)
s = [2, 3, 5, 7]
В Питоне списки могут быть неоднородными, т.е. содержать
элементы разных типов.
s = [1, "two", 3.5]
Базовые типы в Питоне:
1) int -- неограниченное целое число
2) float -- вещественное число в плавающей форме.
ВНИМАНИЕ! Это соответствует числу типа double
в языке C/C++ (8-байтовое вещественное число)
3) bool константы True, False
Логические операции: or, and, not
4) str текстовые строки "abcd", 'e2-e4; c7-c7'
Типа данных char (символ) нет!
Строки содерат символы в кодировке Unicode
Конструкция, самая популярная в Питоне:
list comprehension
Умное задание списка, математическое задание списка
Пример: числа в диапазоне от 1 до 30, которые не делятся на 2, 3, 5
s = { x : x <- [1..30], x не делится на 2, 3, 5 }
s = [x for x in range(1, 31) if x%2 != 0 and x%3 != 0 and x%5 != 0]
множество предикат (фильтр), задающий подмножество
Напечатать квадраты всех натуральных чисел <= 20
[x*x for x in range(1, 21)]
Кортеж, tuple -- это как бы список, элементы которого запрещено менять
s = [2, 3, 5, 7] # list
c = (2, 3, 5, 7) # tuple
Происхождение слова tuple
двойка, тройка, четверка, ... элементов
pair, triple, quadruple, ... tuple
Проверить простоту числа
Реализовать функцию prime(m)
def prime(m):
...
Числа Ферма
F_k = 2^(2^k) + 1
F1 = 5
F2 = 17
F3 = 257
F4 = 65537 -- простые числа
F5 = 2^32 + 1 = 4294967297 -- не простое число, Эйлер разложил его на множители
==========================================
16 Feb 2021
Рассмотрим еще несколько задач на Python3
1. Написать функцию, которая определяет, является ли число m
совершенным (perfect).
Число совершенное, если оно равно сумме своих собственных делителей.
6 = 1 + 2 + 3
28 = 1 + 2 + 4 + 7 + 14
2. Написать функцию, которая определяет, является ли число m
полусовершенным (semiperfect).
Число полусовершенное, если оно равно сумме
некоторого подмножества своих собственных делителей.
Алгоритм сводится к задаче о рюкзаке (knapsack problem).
Дан рюкзак объема V (целое число)
Дано множество S вещей, задан объем каждой вещи x_i, i=1, 2, ..., n
Можно ли рюкзак заполнить так, чтобы не осталось пустого места?
(Найти подмножество T < S такое, что сумма объемов вещей из
T равнялась бы V.)
Классы задач P и NP
Задачи класса P -- задачи, для которых существует алгоритм,
решающий ее за время t <= полином(n), n -- длина записи входа.
Задачи класса NP -- задачи, для которых существует алгоритм
проверки решения, работающий за полиномиальное время от
длины предъявленного решения.
Проблема: совпадают ли классы P и NP?
Большинство математиков верят, что нет, NP != P
но это не доказано. (Если докажут, что это так, т.е. найдут
полиномиальный алгоритм для любой NP-полной задачи, то это
будет КАТАСТРОФА.)
Задача называется NP-полной, если из того, что для нее существует
алгоритм с полиномиальным временем работы вытекает,
что для любой задачи из класса NP тоже существует полиномиальный
алгоритм. (Т.е. это самая сложная задача их класса NP, к
которой можно свести любую задачу из этого класса.)
Задача о рюкзаке NP-полная.
Булевы функции
x1, x2, ..., xn
Принимают значение 0, 1 (или False, True)
Операции or, and, not (логическое сложение, логическое умножение,
отрицание).
(x + y)*(x - z) = x*x + y*x - x*z - y*z
Любая булева функция представляется в виде ДНФ
(дизъюнктивная нормальная форма) --
сумма термов, каждый терм -- произведение либо переменных,
либо их отрицаний
x and (not y) and z or (not x) and y
True False True
Конъюнктивная нормальная форма -- произведение сумм
(x + not y + z + t) and (not x + y + not t)
Каждый сомножитель в КНФ можно интерпретировать как
некоторую аксиому. Произведение -- это набор аксиом,
которые должны одновременно выполняться.
Задача: можно ли найти набор переменных (True, False, ...)
так, чтобы КНФ приняла значение True?
Существует ли модель для данного набора аксиом?
Это задача NP-полная задача. С 2006 года проводится соревнование
между SARS, решающими данную задачу.
Задача 2
Счастливые билеты
Номер автобусного билета состоял из 6 цифр
Билет счастливый, если сумма первых трех цифр равна
сумме последних трех цифр
123456 -- не счастливый билет
123402 -- счастливый билет
Найти число счастливых билетов длины n (n четное).
Решение методом "грубой силы" (тупое):
перебрать все билеты и для каждого билета определить,
является ли он счастливым. (Получим не только число счастливых
билетов, но и список счастливых билетов.)
Цифры билета представляем в виде списка чисел от 0 до 9 длины n
Билет 123402 представляется в виде
[1, 2, 3, 4, 0, 2]
Чтобы перебрать все списки длины n из цифр 0..9,
реализуем функцию next(s), которая находит следующий список
в лексикографическом порядке
123402 -> 123403 -> 123404 -> ... -> 123409 -> 123410 -> ... -> 999999 ->
-> 000000
Задачи на дом (обязательное) -- сделать любую из трех задач
до вторника 2 марта 2021:
1) написать более быстрый алгоритм, который находит
число счастливых билетов длины n (без вычисления самого
списка всех счастливых билетов);
2) написать функцию, которая генерирует список всех перестановок
длины n
def permutations(n):
pass
>>> permutations(3)
[[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1],
[3, 1, 2], [3, 2, 1]]
(хотелось бы реализовать функцию next, которая находит
следующую перестановку в лексикографическом порядке);
3) по теореме Лагранжа, каждое целое число представимо в виде
суммы 4-х квадратов. Написать функцию, которая находит
представление данного числа в виде суммы минимального
числа квадратов (длина представления по теореме Лагранжа <= 4).
18 = 4^2 + 1^2 + 1^2 = 3^2 + 3^2
================================================================
2 Mar 2021
Замечание:
не должно быть кода, который исполняется безусловно на глобальном уровне.
def permutations(n):
res = []
...
# print(p)
res.append(p)
...
return res
if __name__ == "__main__":
n = int("Enter an order of permutations: ")
print(permutations(n))
from permutations import *
============================================================
Тема: применение теории чисел в криптографии
Enigma
Изобретение компьютеров связано с криптографией
Шифры использовались в древнем Египте (подстановочные)
Биекция алфавита
y = E(x), x -- буква алфавита
x = D(y), DE = 1 -- тождественное отображение
E -- Encryption, инъективное отображение
D -- Decryption
Шифр Цезаря -- циклический сдвиг буквы алфавита на 3 позиции
A -> D, B -> E, ..., Z -> C
Шифр Виженера -- перекрытие текста и ключевой фразы,
каждая буква ключевой фразы задает величину
сдвига в соответствии со своей позицией в алфавите.
Симметричные шифры -- ключ шифрования одновременно является и
ключом расшифровки.
Шифр Виженера был сначала взломан Фридрихом Касицким
Знаем длину ключа n
Тогда зашифрованный текст разбивается на n подтекстов
x0 x1 x2 x3 ...
Подтекст 0
x0 xn x_2n ...
Подтекст 1
x1 x_n+1 x_2n+1 ...
. . .
Подтекст i шифруется шифром Цезаря со сдвигом, соответветствующей
i-й букве ключевой фразы => применим метод частотного анализа.
Если биграмма в исходном тексте повторяется на расстоянии n,
она будет повторяться и в зашифрованном тексте.
В тексте на естественном языке вероятность повторения
биграмм намного выше, чем в случайной последовательности букв.
Метод Фридмана определения длины ключевой фразы:
вычисление индекса совпадения.
Втекстенаестественномязыкевероятностьповторения
Втекстенаестественномязыкевероятностьповторения
===
Индекс совпадения 3
3/42 = 7%
Случайная последовательность 1/32 = 3%
Если сдвиг на расстояние, кратное n, то совпадения
в зашифрованном тексте те же самые, что и в исходном тексте =>
их в 2 раза больше, чем при сдвигах на расстояние, не кратное n.
Сдвигаем текст на расстояния 2, 3, 4, 5, 6, 7, 8, 9, 10, ...
Вычисляем индекс совпадения. Как только он увеличивается скачком --
это длина ключевой фразы.
В период первой мировой войны были изобретены
роторные шифровальные машины.
Enigma использовалась Германией во второй мировой войне.
3 ротора
у каждого ротора 26 входных контактов и 26 выходных
внутри они соединялись проводами, ротору соответствовала
перестановка 26 элементов
1-й ротор поворачивался на 1 позицию после каждого нажатия
2-й ротор после полного оборота 1-го ротора
3-1 -- после полного оборота второго.
s1, s2, s3
Дальше был рефлектор:
26 контактов, разбитых на пары, пара контактов была соединена
=> перестановка, которая представляется в виде произведения
независимых транспозиций
r: r^2 = 1 (идемпотентная)
Ток проходил последовательно через 3 ротора,
подавался на рефлектор и с выходного контакта пары
подавался обратно на 3 ротора
y = u1 u2 u3 r s3 s2 s1 (x)
где u_i = s_i^(-1) u1 s1 = 1, u2 s2 = 1
sigma = u1 u2 u3 r s3 s2 s1
sigma^2 = 1
=> для расшифровки сообщения нужно было вновь пропустить
его через машину, восстановив начальное положение роторов.
Ключом шифрования служило начальное положение роторов
(оно задавалось по дню года в шифровальной книге).
В середине 30-х группа польских математиков
под руководством Мариана Решевского начали работы по
расшифровке Enigma -- в частности, развивали теорию
симметрической группы. После окупации все работы были
переданы в Англию, где были продолжены группой
под руководством Алана Тьюринга.
==========================================================
Из всех подобных историй напрашивается вывод,
что любой шифр можно взломать. Эдгар По:
человек не может придумать шифр, который
человеческая изобретательность не сможет разгадать.
Оказалось, что это неверно!
==========================================================
Блочное кодирование / шифрование
Сообщение -- последовательность двоичных разрядов
Разбиваем на блоки равной длины (например, в DES
размер блока 64 бита) и каждый блок шифруем по отдельности.
Симметричное шифрование:
применяются биективные преобразования
из опреденного набора, например, 16 преобразований;
они определяются по ключу шифрования. Их построение --
"запутывание" и "перемешивание" (confusing and diffusing).
Ключ шифрования дает и возможность дешифровки.
3DES -- ключ 128 бит (блок 128).
Схема RSA асимметричного шифрования
===================================
Rumley, Shamir, Adleman -- совершенно открытая,
очень простая.
Для взлома надо решить математическую задачу
разложения 400-значного числа на множители.
(Есть алгоритм Шора факторизации для квантового
компьютера, но пока нет квантового компьютера.)
При блочном кодировании блок можно рассматривать как
двоичную запись целого числа.
E: Z -> Z
Вместо Z рассматриваются кольца вычетов Zm = Z/mZ
Фактор-множество
Отношение эквивалентности ~ на множестве
1) x ~ x (рефлексивность)
2) x ~ y => y ~ x (симметричность)
3) x ~ y, y ~ z => x ~ z (транзитивность)
Множество M с отношением эквивалентности ~
разбивается на классы эквивалентных элементов M/~
Рассмотрим множество Z и число m <- Z, m != 0
x ~ y <=> m | (x - y) или x == y (mod m)
def
Фактор-множество Z/~ обозначается Z_m
Z_m -- это фактор-кольцо
Кольцо (Ring) -- множество с 2-мя операциями +, *
Относительно + это абелева группа
a + (b + c) = (a + b) + c
существует 0: 0 + a = a
Для любого a существует b = -a: a + b = 0
a + b = b + a
В ассоциативных кольцах умножение ассоциативно
Операции связаны дистрибутивностью:
a*(b + c) = a*b + a*c, (b + c)*a = b*a + c*a
Пример кольца: кольцо матриц над Z или над полем C
ПЕРЕРЫВ до 14:30
Если есть алгебраическая структура, то отношение
эквивалентности должно быть согласовано с операциями:
x1 ~ x2, y1 ~ y2 => x1 + y1 ~ x2 + y2
x1*y1 ~ x2*y2
Берем фактор-кольцо, рассмотрим множество элементов,
эквивалентных нулю.
Оно образуем идеал в кольце:
1) подкольцо
2) замкнут относительно умножения на
любые другие элементы кольца
x*i <- I, i*x <- I
Обратно, если I < R -- идеал в кольце R, то можно
рассмотреть отношение эквивалентности
x ~ y <=> x - y принадлежит I
def y = x + i
Все фактор-кольца -- это факторы по разным (двусторонним)
идеалам. R/I
Фактор-группы -- факторы по нормальным подгруппам
H < G
g^(-1) h g <- H
g1 ~ g2 <=> g2 = g1*h
g1 * g2^(-1) <- h
G/H
Пример: кольцо Z_5
{ ..., -15, -10, -5, 0, 5, 10, 15, ... }
{ ..., -14, -9, -4, 1, 6, 11, 16, ... }
{ ..., -13, -8, -3, 2, 7, 12, 17, ... }
{ ..., -12, -7, -2, 3, 8, 13, 18, ... }
{ ..., -11, -6, -1, 4, 9, 14, 19, ... }
Если выбрать из каждого класса по одному представителю,
получим систему остатков. Наиболее популярны
1) неотрицательная система остатков
{ 0, 1, 2, 3, 4 }
2) симметричная система
{ -2, -1, 0, 1, 2 }
Компьютер работает с кольцами Z_m, где
m = 2^32 для чисел int
m = 2^64 для чисел типа long long
неотрицательная система остатков соответствует типу int,
симметричная система остатков соответствует типу unsigned int.
Теорема. Наибольший общий делитель d чисел m и n выражается
в виде
d = u*m + v*n, где u, v -- целые числа
Числа d, u, v можно найти с помощью расширенного алгоритма Евклида.
Идея:
выполняем обычный алгоритм Евклида для (a, b),
в начале (a, b) = (m, n)
при этом храним выражение a, b в виде линейных
комбинаций исходных чисел m, n:
a = u1*m + v1*n
b = u2*m + v2*n
(a, b) <-- (b, r)
a = q*b + r => r = a - q*b
u1' = u2
v1' = v2
r = a - q*b = u1*m + v1*n - q*(u2*m + v2*n) =
= (u1 - q*u2)*m + (v1 - q*v2)*n
u2' = u1 - q*u2
v2' = v1 - q*v2
Расширенный алгоритм Евклида позволяет вычислить
обратный к x элемент в кольце вычетов по модулю m.
Теорема
Элемент x обратим в Zm <=> gcd(x, m) = 1
С помощью расширенного алгоритма Евклида получаем
1 = u*x + v*m
v*m == 0 (mod m)
1 == u*x (mod m) => u -- обратный элемент
Китайская теорема об остатках
m1, m2, ..., mk gcd(mi, mj) = 1 при i != j
m = m1*m2*...*mk
Имеет место изоморфизм колец:
Z_m = Z_m1 + Z_m2 + ... + Z_mk
x = (x1, x2, ..., xk), xi == m (mod m_i)
Операции с подобными строками выполняются покомпонентно
Z_15 = Z_3 + Z_5
7 = (1, 2)
9 = (0, 4)
7 + 9 = 16 =
= 1 = (1, 6) = (1, 1)
Поскольку любое целое число раскладывается в произведение
степеней простых, изучение колец вычетов сводится к
изучению примарных колец вычетов вида
Z_p^e
Пример: сколько квадратных корней из единицы в Z_15?
x^2 == 1 (mod 15)
Ответ: 4
(1, 1)
(-1, 1)
(1, -1)
(-1, -1)
Z_3 -- поле, в поле не больше двух корней из 1
x^2 - 1 = 0
Z_5 -- тоже поле
Пример 2: в кольце Z_105 всего 8 квадратных корней из 1.
m = p1^e1 * p2^e2 * ... * pk^ek
Пример использования китайской теоремы об остатках:
Теорема Эйлера:
Пусть
m = p1^e1 * p2^e2 * ... * pk^ek
Порядок группы обратимы элементов кольца Z_m равен
phi(m) = (p1-1)*p1^(e1-1) * (p2-1)*p2^(e2-1) * ... * (pk-1)*pk^(ek-1)
Пример: m = 100 = 2^2 * 5^2
phi(m) = (2-1)*2 * (5-1)*5 = 40
Обозначение:
U_m -- это группа обратимых элементов кольца Z_m
Из китайской теор. об ост.
phi(U_m) = |U_m| = |U_p1^e1|* ... * |U_pk^ek|
В кольце Z_p^e обратимы элементы, которые взаимно просты с p.
Сколько их? Легче посчитать дополнение к ним.
Элемент не обратим, если он делится на p, таких элементов
n/p, где n -- порядок кольца, n = p^e
p^e/p = p^(e-1)
Обратимы оставшиеся элементы
p^e - p^(e-1) = (p-1)*p^(e-1)
Имеется алгоритм для китайской теоремы об остатках:
Дано: m1, m2, ..., mk -- взаимно простые
r1, r2, ..., rk -- любые числа
Надо: найти x:
x == r_i (mod m_i)
=====================================================
Нужно сделать 1 задачу на тему Теория чисел и криптография
Deadline: 23 марта
=============================================================
Было рассмотрено кольцо вычетов по модулю m
Z_m = Z/mZ
Китайская теорема об остатках
Пусть m = p1^e1 * p2^e2 * ... * pk^ek
Обозначим pi^ei = mi
Z_m == Z_m1 + Z_m2 + ... + Z_mk -- прямая сумма
x == (x1, x2, ..., xk)
x_i == x (mod m_i)
Операции выполняются покомпонентно.
Доказательство -- из теоремы о гомоморфизме
f: Z --> Z_m1 + Z_m2 + ... + Z_mk
x |--> (x (mod m1), ..., x (mod mk))
Ker f = m1*m2*...*mk*Z, поскольку числа m1, ..., mk взаимно простые
Z/ker f = Z_m == Im(f)
|Im(f)| = |Z_m|, |Z_m1 + Z_m2 + ... + Z_mk| = | Z_m | =>
это биекция.
Это доказательство неконструктивное, но есть алгоритм,
который для любого набора остатков x1, x2, ..., xk
находит число x: x == x_i (mod m_i)
Дано: m1, m2, ..., mk -- взаимно простые числа
gcd(m_i, m_j) = 1
x1, x2, ..., xk -- любые целые числа
Надо: найти целое число x:
x == x1 (mod m1), x == x2 (mod m2), ..., x == xk (mod mk)
Идея алгоритма: индукция по k
База:
при k = 1 x = x1
Шаг индукции:
пусть уже вычислили x:
x == x1 (mod m1), x == x2 (mod m2), ..., x == x_k-1 (mod m_k-1)
Нужно исправить x так, чтобы выполнилось сравнение
x' == x_k (mod m_k)
Возьмем
x' = x + s*(m1*m2*...*m_k-1)
x + s*(m1*m2*...*m_k-1) == x_k (mod m_k)
Число m1*m2*...*m_k-1 обратимо по модулю m_k
u * (m1*m2*...*m_k-1) == 1 (mod m_k)
Домножим уравнение
s*(m1*m2*...*m_k-1) == x_k - x (mod m_k)
на u.
s*(m1*m2*...*m_k-1)*u == (x_k - x)*u (mod m_k)
s == (x_k - x)*u (mod m_k)
Из Китайской теоремы об остатках выводится теорема Эйлера
U_m -- группа обратимых элементов кольца Z_m
phi(m) = |U_m|
Из Китайской теоремы об остатках следует, что
U_m == U_m1 x U_m2 x ... x U_mk
m_i -- это степень простого числа pi^ei
phi(p^e) = p^e - p^(e-1) = (p - 1)*p^(e-1)
Частный случай: пусть m свободно от квадратов
m = p1 * p2 * ... * pk
phi(m) = (p1 - 1)*(p2 - 1)*...*(pk - 1)
Малая теорема Ферма:
Пусть p -- простое число, b -- любое целое число такое, что
b не делится на p.
Тогда b^(p-1) == 1 (mod p)
Другая форма Малой теоремы Ферма:
Пусть p -- простое число. Тогда
Для всякого целого числа x
x^p == x (mod p)
Также x^(s*(p-1) + 1) == x (mod p)
Обобщение Малой теоремы Ферма -- теорема Эйлера:
Пусть m != 0 -- произвольное целое число.
Пусть b взаимно просто с m, т.е. gcd(b, m) = 1
Тогда b^phi(m) == 1 (mod m)
(это также следствие теоремы Лагранжа).
Нам понадобится частный случай случай теоремы Эйлера:
Пусть число m свободно от крадратов:
m = p1 * p2 * ... * pk
Тогда для всякого целого x и для всякого s >= 0
x^(s*phi(m) + 1) == x (mod m)
Доказательство: следствие из Китайской теор. об остатках,
Малой теоремы Ферма во 2-й формулировке и формулы Эйлера.
x = (x1, x2, ..., xk)
phi(m) = (p1 - 1)*...*(pk - 1)
x_i^(s*phi(m) + 1) = x_i^(t*(p_i - 1) + 1) == x_i (mod p_i)
============================================================
Для схема RSA нужен случай теоремы Эйлера для k = 2
m = p*q, p, q -- простые
phi(m) = (p - 1)*(q - 1)
Теорема Эйлера:
для любого x <- Z_m и любого неотрицательного s <- Z
x^(s*phi(m) + 1) == x (mod m)
Схема кодирования с открытым ключом RSA
Основана на сложности задачи разложения числа на множители
при том, что мы можем легко находить большие простые числа.
Генерируются два больших простых числа p, q
(размер не меньше 200 десятичных цифр).
m = p*q -- открыто, но разложение m = p*q -- секретно.
Имеем phi(m) = (p - 1)*(q - 1)
Выберем любое более-менее случайное число e <- U_phi(m)
e обратимо в кольце Z_phi(m)
Отображение шифрования:
E: Z_m -> Z_m
E: x -> x^e (mod m)
Вычислим обратный к e элемент в кольце Z_phi(m)
d * e == 1 (mod phi(m))
Отображение дешифровки:
D: Z_m -> Z_m
D: y -> y^d (mod m)
Покажем, что D и E -- взаимно обратные (биективные) отображения.
Любое x
y = E(x)
Применим D к y
D(E(x)) == D(x^e) == (x^e)^d == x^(e*d)
Имеем e*d == 1 (mod phi(m))
e*d = 1 + s*phi(m)
По теореме для m, свободного от квадратов
x^(e*d) = x^(1 + s*phi(m)) == x (mod m).
Открытый ключ (public key): пара (m, e)
Секретный ключ (private key): пара (m, d)
Шифрование: сообщение x <- Z_m
отображается в y = x^e (mod m)
Расшифровка: зашифрованное сообщение y <- Z_m
отображается в x = y^d (mod m).
Иллюстрация на конкретном примере
p = 97, q = 101
m = p*q = 9797
phi(m) = (p - 1)*(q - 1) = 9600
Выберем e, взаимно простое с phi(m)
e = 1021
Обратный к e элемент в Z_phi(m)
d = 7381
Имеем e*d == 1 (mod phi(m))
Открытый ключ: (m=9797, e=1021)
Секретный ключ: (m=9797, d=7381)
Возьмем произвольное сообщение
x = 1234
Зашифрованное сообщение
y = x^e (mod m) = 1234^1021 (mod 9797) = 3868
Расшифруем y:
z = y^d (mod m) = 3868^7381 (mod 9797) = 1234.
===================================================
Какие задачи возникают в связи с этой схемой?
Прямая задача -- генерация ключа:
Надо сгенерировать 2 больших простых числа.
Достаточно иметь тест простоты числа, который
работает быстро для больших чисел. По теореме
Чебышева
pi(n) = количество простых чисел <= n
pi(n) ~ n/ln(n) при n -> infinity
В среднем в отрезке натурального ряда с числами порядка n
длины ln(n) находится одно простое число. Простых чисел
очень много. Можно проверять случайные числа,
мы очень быстро найдем простое.
Обратная задача -- взлом ключа. Для этого надо вычислить d:
d*e == 1 (mod phi(m)),
для этого надо вычислить phi(m), а для этого
надо разложить число m на множители, m = p*q, phi(m) = (p-1)*(q-1).
(Факторизация и вычисление функции Эйлера -- эквивалентные задачи.)
Две задачи
1. Тест простоты числа m.
2. Факторизация (разложение на множители) числа m.
Перерыв до 14:20
ТЕСТЫ ПРОСТОТЫ
Может ли Малая теорема Ферма служить тестом простоты?
m -- простое => для любого b, 1 < b < m,
b^(m-1) == 1 (mod m)
Необходимое условие простоты.
Является ли оно достаточным? Греки думали, что да, и даже
для b = 2
Ошибочно считалось, что, если
2^(m-1) == 1 (mod m) => m простое.
Первым показал, что это не так, франц. математик Сарруз в 1820 г.
m = 341
2^340 == 1 (mod 341)
2^340 = (2^10)^34 = 1024^34 == 1^34 == 1
1024 = 341*3 + 1 == 1 (mod 341)
341 -- не простое 341 = 11*31
Робер Каймайкл (1910) нашел такое странное число m,
не являющееся простым,
которое ведет себя как простое в Малой теореме Ферма
для для любых оснований b, взаимно простых с m:
gcd(b, m) = 1
b^(m-1) == 1 (mod m)
В честь него такое числа называются кармайкловыми.
Недавно доказано, что множество кармайкловых чисел
бесконечно
Первые 3 кармайкловых числа: 561, 1105, 1729, ...
Теорема. Число m кармайклово <=>
1) m свободно от квадратов и
m = p1*p2*...*pk, где k >= 3;
2) для всякого i
(p_i - 1) | (m - 1)
Проверим для m = 561 = 3 * 11 * 17
2 | 560,
10 | 560?
16 | 560 (16*35 = 560)
Доказательство достаточности <=
из Китайской теоремы об остатках
Пусть b: gcd(b, m) = 1
b == (b1, b2, ..., bk), b_i == b (mod p_i)
b^(m-1) == (b1^(m-1), ..., bk^(m-1))
Так как (p_i - 1) | (m - 1)
b_i ^ (m - 1) == 1 (mod p_i)
Малая теорема Ферма не подходит в качестве теста простоты
для кармайкловых чисел.
Тем не менее идея, связанная с тестом Ферма, работает в
построения вероятностного теста простота Миллера-Рабина.
Вероятностный тест простоты для числа m выдает один из двух ответов:
1) m составное (при этом предъявляется доказательство того,
что m составное);
2) не знаю.
При этом для составного m вероятность второго ответа не
больше чем 0.25.
Пример: пусть для числа m тест 100 раз выдал ответ
"не знаю". Тогда вероятность того, что оно сотавное,
не превышает (1/4)^100, т.е. оно с вероятностью почти
1 простое. Но доказательства простоты мы, увы, не получаем.
Тем не менее для практических целей такие псевдопростые числа
годятся.
Тест Миллера-Рабина основан на том, что в любом поле
не больше двух корней из единицы: 1 и -1.
Теоремы Безу: число корней многочлена в поле не
превосходит его степени
x^2 - 1 = 0
m простое <=> кольцо вычетов Z_m является полем.
(Если число m составное и нечетное, то число квадратных
корней из 1 >= 4; x^2 - 1 == 0 (mod m), то вероятность,
что x = +- 1, не меньше чем 1/2).
Тест Миллера--Рабина
Дано нечетное число m
Представляем число m - 1 = 2^s * t, где t -- нечетное число.
Берем случайное основание b <- Z_m
и вычисляем следующую последовательность:
b_0 = b^t (mod m)
b_1 = b_0^2 (mod m)
b_2 = b_1^2 (mod m)
. . .
b_s = b_(s-1)^2 = b^(m-1) (mod m)
Рассматриваем эту последовательность:
b_0, b1, b2, ..., b_s
Если
1) эта последовательность не заканчивается единицей, b_1 != 1 (mod m)
то число составное (не выполняется Малая теорема Ферма);
2) если в последовательности есть фрагмент
*, *, *, ..., *, 1, 1, ..., 1
====
где звездочками обозначены элементы, отличные от +-1,
то число составное;
3) ответ "не знаю" выдается в двух случаях
a) последовательность состоит из одних единиц
1, 1, 1, ..., 1
b) в последовательности имеется фрагмент
*, *, ..., *, -1, 1, ..., 1
Доказательство того, что тест не врёт:
в случае 1 не выполняется Малая теорема Ферма
в случае 2 имеем
b_j != +- 1, (b_j)^2 = 1 (mod m) => имеем нетривиальный
квадратный корень из 1 =>
Z_m не является полем =>
m не простое.
==========================================================
Красивые алгоритмы факторизации
1. Колесо со спицами: делим m на 2, 3, 5, если не делится,
то последовательно делим на пробные делители вида
30*k + r_i,
где r_i обратимы в Z_30
r_i <- { 1, 7, 11, 13, 17, 19, 23, 29 }
2. rho-алгоритм Полларда
3. p-1-алгоритм Полларда (методы Монте-Карло)
-------------------------------------------------
4. p-1 алгоритм Ленстра на эллиптических кривых.
5. Квадратичное решето (Диксон).
6. Вариант квадратичного решета в квадратичном расширении
конечного поля Z_p^2
-------------------------------------------------
Алгоритм Шора факторизации на квантовом компьютере.
=============================================================
=============================================================
16 Mar 2021
Было рассказано:
схема шифрования с открытым ключом RSA
m = p*q m открыто, разложение m секретно
phi(m) = (p - 1)*(q - 1)
e -- любое число, обратимое в Z_phi(m) -- открытый ключ
d -- обратный элемент к e в Z_phi(m)
d*e == 1 (mod phi(m))
Кодирующая процедура
E: x -> x^e (mod m)
Декодирующая процедура
D: y -> y^d (mod m)
D*E = E*D = 1 -- тождественное преобразование в Z_m
Открытый ключ -- пара (m, e)
Секретный ключ -- пара (m, d)
Для взлома ключа надо по е и m вычислить обратный элемент к e
в кольце Z_phi(m)
Для этого надо знать phi(m), а для этого
надо разложить m на множители.
Алгоритмы факторизации
1. Колесо со спицами -- обобщение "детского" алгоритма
деления на все нечетные пробные делители до квадратного
корня.
Сначала делим на 2, если не делится, то делим на все нечетные числа,
начиная с трех.
2, 3, 5 их произведение = 30
Сначала делим на 2, 3, 5, если не делится, то делим
на пробные делители вида 30*k + t_i, k = 0, 1, 2, ...
t_i -- числа в диапазоне 1..29, обратимые в кольце Z_30
t = [1, 7, 11, 13, 17, 19, 23, 29]
2. Алгоритмы Полларда, "методы Монте-Карло".
Первый метод -- rho-алгоритм факторизации Полларда.
Основан на поиске цикла в рекуррентной последовательности
Дано: число m
Надо: разложить число m на множители (найти нетривиальный
делитель числа m).
Работаем в кольце Z_m
Рассмотрим полиномиальное отображение
f: Z_m -> Z_m
f(x) = x^2 + 1 (mod m)
Иногда берут
f(x) = x^2 - 1 (mod m)
Берем любую начальную точку x0 <- Z_m
Вычисляется бесконечная рекуррентная последовательность
x0, x1 = f(x0), x2 = f(x1), x3 = f(x2), ...
x0, x1, x2, x3, x4, ...
Так как кольцо Z_m конечно, то последовательность периодическая.
У нее может быть начальный непериодический отрезок
и дальше повторяющийся цикл
x0, x1, x2, x3, x4, ..., x_i, ..., x_j = x_i, ...,
========> повторяющийся цикл
Пусть p | m. Делитель p нам неизвестен, но теоретически
можно рассмотреть ту же самую последовательность по модулю p
(в кольце Z_p):
x'0, x'1, x'2, ...
x'_i == x_i (mod p) Z_m -> Z_p
Это тоже периодическая последовательность и, поскольку |Z_p| < |Z_m|,
с большой долей вероятности период последовательности { x'_i }
меньше, чем период последовательности { x_i }.
Значит, найдутся два элемента x'_i == x'_j (mod p), при этом
x_i != x_j (mod m)
Это означает, что p делит (x_i - x_j)
m не делит (x_i - x_j)
p делит gcd(x_i - x_j, m)
m не делит gcd(x_i - x_j, m)
Значит, gcd(x_i - x_j, m) нетривиален
Задача сводится к нахождению цикла в последовательности { x_i }
по модулю p.
Сравниваем
x0 <--> x1
x0 <--> x2
x0 <--> x3
. . .
Это не работает, поскольку последовательность может входить в цикл
не сразу. Для преодоления этого можно сравнивать два элемента,
расстояние между которыми на каждом шаге увеличивается на 1 и
первый элемент сдвигается вперед по последовательности
x0 <--> x1
x1 <--> x3
x2 <--> x5
x3 <--> x7
. . .
В качестве сравнения используем вычисление gcd(x_i - x_j, m)
Предложение:
Пусть f: Z_n -> Z_n -- полиномиальное отображение,
где степень полинома больше 1
Тогда средняя длина цикла орбиты произвольного элемента
равна O(sqrt(n))
Оценка скорости работы алгоритма
t = O(sqrt(p)), где p -- минимальный простой делитель числа m
t <= O(m^(1/4))
Пусть m -- 20-значное десятичное число
Тогда число шагов алгоритма не превосходит 10^5 = 100,000
=====================================================
p-1-алгорить факторизации Полларда
Основан на Малой теореме Ферма
Пусть p -- простое число, b -- любое число, взаимное простое с p.
Тогда b^(p-1) == 1 (mod p)
Дано: число m
Надо: найти нетривиальный делитель числа m
Пусть p | m -- простой делитель m
Пусть также p - 1 -- это гладкое число (smooth number)
Гладкость определяется для константы N
число n является N-гладким, если оно раскладывается на небольшие
множители:
n = n1*n2*...*n_k, где n_i -- степени разных простых чисел,
n_i <= N.
Например, 100 = 2^2 * 5^2 = 4 * 25
Следовательно, число 100 является гладким для N = 30,
100 делит 30!
Поскольку p-1 -- гладкое число, то (p-1) | N!
N! = (p-1)*s
Рассмотрим случайное число b
b^(N!) == b^((p-1)*s) == (b^(p - 1))^s == (1)^s (mod p) == 1 (mod p)
=>
p | ( b^(N!) - 1 )
p | gcd( b^(N!) - 1, m )
Значит, gcd(b^(N!) - 1, m ) нетривиален
Получается, что, если p-1 гладкое, то построенная на
его основе схема RSA кодирования с открытым ключом ненадежная!
m = p*q можно разложить на множители p-1 алгоритмом Полларда
Простое число p надежное, если p-1 имеет большой простой делитель.
Есть полиномиальный алгоритм генерации больших простых чисел
с сертификатами простоты.
Сертификат -- это простой способ доказательства, что число простое,
которое любой человек может проверить.
ПЕРЕРЫВ до 14:30
Перерыв закончился
Утверждение: число p простое <=> группа обратимых элементов
кольца Z_p является циклической группой
порядка p - 1
Теорема. Группа обратимых элементов любого конечного поля является
циклической.
G^e = 1 Z_m x Z_n циклическая <=> gcd(m, n) = 1
G = Z_4 + Z_2 + Z_25 -- не циклическая |G| = 8*25
x^(4*25) = 1
Поле из p^n = m элементов
Группа обратимых элементов имеет порядок p^n - 1 = m - 1
Не циклическая => e | (m - 1), e < (m - 1)
x^e - 1 = 0
По теореме Безу уравнение имеет не больше коней, чем его степень
=> число корней <= e < m - 1
Группа ненулевых элементов поля Z_p циклическая.
Она порождается некоторым элементов a
1, a, a^2, a^3, ..., a^(p-2)
Такой элемент a называется первообразным корнем степени p.
Его можно использовать для построения датчика случайных чисел.
Берем любое начальное значение x0 <- Z_p
"Случайная последовательность":
x0, x1 = x0*a (mod p), x2 = x1*a (mod p), ...,
каждый следующий элемент получается умножением предыдущего на a.
Период такой последовательности равен p-1.
Сертификат простоты числа p
Надо предъявить a: a^(p-1) == 1 (mod p)
Для любого k < p-1 a^k != 1 (mod p)
Для проверки последнего утверждения
достаточно предъявить разложение p - 1 на простые множители
p - 1 = q1^e1 * q2^e2 * ... *q_s^e_s
достаточно проверить, что для всякого простого q_j | (p-1)
a^((p-1)/q_j) != 1 (mod p)
Сертификат простоты числа p -- это первообразный корень a
плюс разложение на множители числа
p - 1 = q1^e1 * q2^e2 *... * q_s^e_s
а также набор сертификатов для всех простых чисел q_j, j = 1, 2, ..., s.
----------------------------------------------
Задача: вычислить квадратный корень из заданного элемента a
поля Z_p вычетов по модулю p
Критерий того, что из элемента a можно извлечь квадратный корень:
ненулевой элемент a <- Z_p является квадратом <=>
a^( (p-1)/2 ) == 1 (mod p)
Доказательство.
Группа ненулевых элементов циклическая, порождается элементом b
=> квадраты -- это элементы вида b^(2*k), т.е. четные степени элемента b.
Если a = b^(2*k), то
a^((p-1)/2) = b^( 2*k*(p-1)/2 ) = b^( k*(p-1) ) =
= (b^(p-1))^k == 1 (mod p)
Обратно, если a не квадрат, то a = b^(2*k + 1),
a^((p-1)/2) = b^((2*k+1)*(p-1)/2)
Имеем (2*k + 1)*(p - 1)/2 = k*(p - 1) + (p - 1)/2
b^((2*k+1)*(p-1)/2) = b^(k*(p - 1)) * b^(p - 1)/2
По Малой теореме Ферма, b^(k*(p - 1)) == 1, значит
b^(k*(p - 1)) * b^(p - 1)/2 == b^(p - 1)/2 != 1 (mod p),
поскольку порядок элемента b равен p - 1.
Алгоритм вычисления квадратного корня в поле Z_p
================================================
Дано: число a, которое является квадратом по модулю p, где
p -- простое число
Надо: вычислить r: r^2 == a (mod p)
Рассмотрим 2 случая
1) p == -1 (mod 4) (или p == 3 (mod 4))
Это простой случай, сразу можем выдать ответ:
r == a^((p+1)/4) (mod p)
Действительно,
r^2 == a^((p+1)/2) == a^((p-1)/2 + 1) ==
== a^((p-1)/2)*a == a (mod p),
поскольку, т.к. a является квадратом, то a^((p-1)/2) == 1
2) p == 1 (mod 4)
Это сложный случай
(p-1)/2 = q*(2^s), q -- нечетное число
Рассмотрим последовательность
a^q, a^(2*g), a^(4*q), ..., a^(2^s*q) = a^((p-1)/2) = 1
Эта последовательность начинается с a^q и каждый следующий
ее элемент есть квадрат предыдущего.
Обозначим a^q = t. Рассматривается последовательность
t, t^2, t^4, t^8, ..., t^(2^s) = 1
Эта последовательность имеет один из двух видов
1, 1, 1, ..., 1
если она начинается с единицы,
или в ней обязательно есть -1:
*, *, ..., *, -1, 1, ..., 1
t, t^2, ...,
(если x != +-1, x^2 == 1 (mod p), то p не может быть простым,
поскольку
x^2 - 1 == 0, (x - 1)*(x + 1) == 0 (mod p)
p | (x - 1)*(x + 1) =>
p | (x - 1) или p | (x + 1) => x == 1
или x == -1)
Первый случай: последовательность начинается с 1
t == a^q == 1 (mod p), q нечетное
r == a^((q + 1)/2) (mod p)
r^2 = a^(q + 1) = a^q * a = a
Мы нашли искомый ответ r.
Пусть t == a^q != 1 (mod p)
Вычислим r == a^((q + 1)/2) (mod p)
r^2 == a^(q + 1) = a^q * a = t*a
Имеем ситуацию
t -- корень из 1 порядка 2^k, k >= 1 (если t последовательно
возводить в квадрат, то рано или поздно получим 1)
r такое, что
r^2 = t*a
Идея: у нас есть последовательность для t != 1
t, t^2, t^4, t^8, ..., -1, 1
Подправим t, т.е. заменим на t' = t*b для некоторог элемента b^2
так, чтобы длина этой последовательности уменьшилась.
Также при этом подправим r, т.е. заменим r на r' = r*b
так, чтобы по прежнему выполнялось равенство
r'^2 = t'*a
Будем так действовать до тех пор, пока t не станет равным 1
Возьмем какой-нибудь элемент, не являющийся квадратом.
z: z не квадрат (mod p)
Рассмотрим последовательность
z^q, z^(2*q), z^(4*q), ..., z^((p-1)/2) == -1, z^(p-1) == 1
которая начинается с z^q, а каждый следующий элемент в ней
является квадратом предыдущего. Заканчивается последовательность
единицей, предпоследний элемент в ней равен -1.
Последовательность квадратов для z^q имеет длину строго больше,
чем последовательность квадратов для a^q.
Последовательность для t = a^q:
t, t^2, ..., -1, 1
*, *, ..., *, -1, 1
Последовательность для z^q более длинная:
*, *, *, *, ..., *, -1, 1
*, b, b^2,*, ..., *, -1, 1
Через b обозначим элемент второй последовательности,
который отстоит от конца расстояние, равное
длине последовательности для a^q. Положим
t' == t*b^2 (mod p),
r' == r*b (mod p).
Соотношение между t' и r' остается верным:
r'^2 == t'*a,
поскольку из соотношения r^2 = t*a вытекает, что
r'^2 = (r*b)^2 = r^2 * b^2 = t*a * b^2 = t'*a
Последовательность квадратов для t' имеет вид:
t*b^2, t^2 * b^4, ..., (-1)*(-1) = 1
Она имеет длину, строго меньшую, чем длина
последовательности квадратов для t (поскольку она есть
произведение двух последовательностей, в которых -1 стоит
на одном и том же предпоследнем месте, при перемножении
на предпоследнем месте получим единицу).
Пример:
p = 113
a = 61 = 20^2 (mod p)
Вычислим квадратный корень из a = 61 в кольце Z_113
p - 1 == 112 == 7*2^4
t == a^q == 69
r == a^((q+1)/2) == 64
r^2 == t*a
Последовательность квадратов для t:
69, 15, 112, 1
Находим не квадрат: z = 3
Для z вычисляем последовательность квадратов:
40, 18, 98, 112, 1
Положим b = 40, b^2 = 18.
t' = t*18 = 112
r' = r*40 = 74
r'^2 == t'*a
Последовательность квадратов для t':
112, 1
40, 18, 98, 112, 1
*, *, b, b^2, 1
Повторяем шаг алгоритма:
b == 98, b^2 == 112
t'' == t'*b^2 == t'*112 == 1
r'' == r'*b == r'*98 == 20
r''^2 == t''*a == a
Ответ:
r'' = 20
========================================================
23 Mar 2021
Работа с матрицами в Питоне
В Питоне реализованы целые числа неограниченной
длины => можно реализовать рациональные числа
class Fraction в модуле fractions
Матрица в Питоне представляется массивом (списком) строк,
каждая строка -- это в свою очередь список чисел.
a = [[1, 2, 3], [4, 5, 6], [7, 8, 0]]
Элемент матрицы с индексами i, j записывается как
a[i][j]
Реализуем метод Гаусса приведения матрицы к ступенчатому
виду для
1) матрицы рациональных чисел (Fraction);
2) матрицы вещественных чисел (float).
def echelonForm(a):
a -- рациональная матрица
Функция не изменяет матрицу a и возвращает пару
(b, rank)
где b -- ступенчатый вид матрицы a,
rank -- число ненулевых строк матрицы b
b = echelonForm(a) -- НЕПРАВИЛЬНО!
(b, rank) = echelonForm(a)
ПЕРЕРЫВ до 14:15
Особенности алгоритма Гаусса в применении в вещественной матрице
1. Нельзя сравнивать элементы с нулем, т.к. вещещственные числа
представлены с ограниченной точностью. Вместо сравнения с нулем
используем сравнение абсолютной величины разности с eps
вместо x == y
используем abs(x - y) <= eps
вместо x != y
используем abs(x - y) > eps
2. В качестве разрешающего элемента (на который мы делим)
нужно выбрать максимальный по абсолютной величине
элемент в столбце, чтобы в процессе алгоритма
умножение выполнялось на число, не превосходящее
единицы по абсолютной величине.
Домашнее задание по матрицам
4 варианта
1. Вычислить обратную матрицу к рациональной матрице
def inverseOfRationalMatrix(a):
. . .
Функция должна вернуть обратную матрицы или
возбудить исключение типа ValueError для прямоугольной
или вырожденной квадратной матрицы.
2. Вычислить обратную матрицу к вещественной матрице
def inverseOfRealMatrix(a):
. . .
Функция должна вернуть обратную матрицы или
возбудить исключение типа ValueError для прямоугольной
или вырожденной квадратной матрицы.
3. Решить систему линейных уравнений над рациональными числами
с невырожденной квадратной матрицей
def solveRationalLinearSystem(a, b):
. . .
где a -- квадратная невырожденная матрица,
b -- массив свободных членов.
Функция должна вернуть массив, который представляет
решение системы.
a = [[1, 2], [3, 4]]
b = [5, 6]
x = solveRationalLinearSystem(a, b)
print(x)
[-4, 9/2]
4. Решить систему линейных уравнений над вещественными числами
с невырожденной квадратной матрицей
def solveRealLinearSystem(a, b):
. . .
Идея вычисления обратной матрицы:
1 2 3 | 1 0 0 * * * | * * *
4 5 6 | 0 1 0 ==> 0 * * | * * * ==> обратный проход
7 8 0 | 0 0 1 echelonForm 0 0 * | * * *
* * 0 | * * * * 0 0 | * * *
0 * 0 | * * * ==> 0 1 0 | * * * ==>
0 0 1 | * * * 0 0 1 | * * *
1 0 0 | * * *
0 1 0 | * * * Обратная матрица
0 0 1 | * * *
==================================================
Новая тема: классы в Питоне
Отличие классов в Питоне от классов в C++ или в Java
В С++ класс -- это как бы чертеж объектов класса,
и все объекты класса устроены одинаково (автомобили,
созданные по одному и тому же чертежу).
Это сoответствует статической типизации языков C++, Java, Haskell, ...
Python -- язык с экстремально динамической типизацией:
ничего не возможно понять, пока пока программа не заработает.
Разные объекты одного и того же класса могут иметь совсем
разное устройство.
=====================================
30 Mar 2021
Классы в Питоне
В С++
class R2Vector {
public:
double x;
double y;
public:
R2Vector(double xx = 0., double yy = 0.): // this -- неявная переменная
x(xx), // указатель на объект класса,
y(yy) // к которому применяется метод
{}
. . .
};
Python3:
class R2Vector:
def __init__(self, x = 0.0, y = 0.0):
if type(x) == list or type(x) == tuple:
self.x = float(x[0])
self.y = float(x[1])
elif type(x) == R2Vector:
self.x = x.x
self.y = x.y
else:
self.x = float(x)
self.y = float(y)
==================================
В Питоне, помимо списков и кортежей, есть еще одна встроенная
структура данных -- словарь, dict (от dictionary).
Это соответствует классу std::map в стандартной библиотеке
STL языка C++
Словарь содержит множество ключей, каждый ключ отображается
на некоторое значение. Ключи уникальны (ключ в словаре не
повторяется). Словарь можно рассматривать как ассоциативный
массив:
d = { "abcd": 1234, 77: "efgh" }
>>> d = { "abcd": 1234, 77: "efgh" }
>>> d
{'abcd': 1234, 77: 'efgh'}
>>> d["abcd"]
1234
>>> d[77]
'efgh'
>>> d[888]
Traceback (most recent call last):
File "", line 1, in
KeyError: 888
Добавление ключа и его значения в словарь:
>>> d[123] = 99
>>> d
{'abcd': 1234, 77: 'efgh', 123: 99}
Проверим, содержится ли ключ в словаре:
>>> "abcd" in d
True
>>> 99 in d
False
>>> 123 in d
True
Удаление ключа и его значения:
>>> d
{'abcd': 1234, 77: 'efgh', 123: 99}
>>> del(d[77])
>>> d
{'abcd': 1234, 123: 99}
Преимущество использования словаря заключается в том,
что в нем все операции выполняются быстро, за время
порядка O(log(n)), где n -- число ключей в словаре.
Реализация словаря основан на сбалансированном дереве
(красно-черные деревья или AVL-деревья?).
===================================================
Продолжим определение класса R2Vector
ПЕРЕРЫВ до 14:20
В модуле R2Graph реализованы классы R2Vector и R2Point,
представляющие векторы и точки на плоскости R2. С их помощью
можно решать большинство геометрических задач.
Что такое точка (элемент аффинного пространства)?
Разность двух точек -- вектор,
точка плюс/минус вектор -- точка.
p1*c1 + p2*c2 -- получаем отрезок с концами p1, p2
0 <= c1, c2 <= 1
c1 + c2 == 1
=====================================================
6 Apr 2021
Графическое и оконное программирование на Python'е
Нужно использовать какой-либо модуль, поддерживающий
оконную графику.
Модуль на основе kde
Common Desktop Environment CDE
OK All Correct GNU KDE = Window Manager + библиотека функций для
работы в оконной системе
GNOME -- конкурент KDE
XWindow -- классическая оконная система, которую используют в Unix
и в VAX VMS/Open WMS
Mac OS-X (Apple Macintosh) (система типа Unix на базе микроядра)
использует оконную систему Quartz.
В дистрибутиве Python'а для MS Windows есть только один модуль для
поддержки графики: tkinter (реализует старый toolkit TCL/TK).
tkinter -- очень простой модуль (для детского сада), программировать
простые оконные приложения на нем очень просто и удобно.
Недостаток -- меньше возможностей (по определению меньше, чем в PyQT).
В частности, не реализовано сглаживание линий (antialiasing).
Мы будем использовать tkinter.
from tkinter import *
# Создать олно верхнего уровня
root = Tk()
За размещение дочерних окон отвечает объект LayoutManager.
В tkinter есть 3 типа LayoutManager:
-- place (самый примитивный) -- просто указываются координаты окна,
окно остается на фиксированном месте
при изменении геометрии родительского окна;
-- pack (упаковщик) -- он уже может изменять расположение дочерних окон;
-- grid (решетка) -- дочерние окна располагаются в клетках таблицы.
Мы используем pack
ПЕРЕРЫВ до 14:15
====================================================
13 Apr 2021
Оконное программирование на Python3 + tkinter
Как "оживить" картинку: треугольник, возможность перетаскивания
его вершин и перерисовки катринки в реальном времени
Программа triangle4.py
Анимация в оконных программах --
находится в противоречии с идеологией оконного программирования
Классическая оконная программа занимается лишь обработкой событий,
связанных с окнами.
Оконная программа имеет очередь сообщений о событиях и последовательно
обрабатывает сообщения из начала очереди. Если очередь пуста, то программа
ничего не делает -- она ждет наступления события.
tkinter, Qt используют очень простое решение для анимации:
таймер, который посылает сигнал через заданный промежуток времени.
Обработчик сигнала может что-то поменять в данных и вызвать функцию
перерисовки окна, а затем запланировать следующий сигнал таймера.
В tkinter можно просто запланировать событие для окна и
вызов его обработчика через заданный интервал времени
(в миллисекундах)
drawArea.after(timeDelayInMilliseconds, functionToPerform)
Рассмотрим на примере цифровых часов
Файл dclock.py
ПЕРЕРЫВ до 14:30
Временная зона для Пекина Beijing -- такой нет.
Есть для Шанхая 'Asia/Shanghai'
Временная зона для Москвы: "Europe/Moscow"
Кривая Безье = аналог многочленов Бернштейна в двумерном случае
Многочлены Бернштейна использовались для доказательства
теоремы Вейерштрасса:
любая непрерывная функция на отрезке может быть равномерно
приближена многочленами
для любого eps существует многочлен
|f(x) - p(x)| <= eps на отрезке [a, b]
Кривая Безье задается точками
p0, p1, ..., pn
Рекуррентное определение
B{p0, p1, ..., pn}: [0, 1] --> R2
B{p0, p1, ..., pn}(t) -- кривая
1) B{p0}(t) = p0
2) B{p0, p1}(t) = (1-t)*p0 + t*p1
3) B{p0, p1, p3}(t) = (1-t)*B{p0, p1}(t) + t*B{p1, p2}(t)
4) B{p0,p1,p3,p4}(t) = (1-t)*B{p0,p1,p2}(t) + t*B{p1,p2,p3}(t)
. . .
4) B{p0,p1,...,pn}(t) = (1-t)*B{p0,p1,...,p_n-1}(t) + t*B{p1,p2,...,pn}(t)
Через биномиальные коэффициенты:
B{p0,p1,p3,p4}(t) = (1-t)^3*p0 + 3*(1-t)^2*t*p1 + 3*(1-t)*t^2*p2 + t^3*p3
B{p0,p1,...,pn}(t) = сумма_k=0^n C(n,k)*(1-t)^(n-k)*t^k*p_k
=========================================================
20 Apr 2021
Синтаксический разбор:
реализация сканера (лексический анализатор)
с помощью механизма регулярных выражений
и рекурсивная реализация парсера (синтаксического анализатора)
Была разобрана программа, рисующая график функции:
plotFunc.py
Исправления в plotFunc.py
Возможность рисовать график разрывной или не всюду определенной функции.
Идея реализации компилятора формул
Цель: мы преобразуем обычную infix-ную запись формулы в
в обратную польскую запись.
Reverse Polish Notation, или сокращение RPN
Значение RPN легко вычисляется с помощью стекового вычислителя.
2*(5 - 7*3) + 21*3 --> преобразуем в RPN:
2, 5, 7, 3, *, -, *, 21, 3, *, +
1 2 3 4 5 6 7 8 9 10 11
Вычисление с помощью стека:
1) 2
2) 2, 5
3) 2, 5, 7
4) 2, 5, 7, 3
5) 2, 5, 21
6) 2, -16
7) -32
8) -32, 21
9) -32, 21, 3
10) -32, 63
11) 31
Современные объектно ориентированные языки используют обратную
польскую запись для промежуточного представления программы:
Java -> bytecode,
C# -> CIL (Common Intermediate Language for .Net project)
Python -> *.pyc (что-то типа bytecode)
Программа может выполняться двумя способами:
1) интерпретация
2) компиляция "на лету" (JIT - Just In Time compiling)
Устройство компилятора
======================
Компиляция выполняется в 2 этапа:
1) представляет входной текст в виде потока лексем
Лексем == слово входного языка.
Лексика -- строение слов языка.
Синтаксис -- строение фраз языка.
Синтаксис описывается в терминах лексем.
Лексема имеет тип (число, имя функции, +, -, *, /, ^ или **,
скобки, конец строки) и значение
Лексема "число" может иметь разные численные значения
Лексема "имя функции" имеет текстовое значение
На первом проходе мы преобразуем поток символов в поток лексем:
sin(2*x + 3) - 1
Список лексем:
(NAME, "sin", 0.), (LPAR, "", 0.), (NUMBER, "", 2.),
(MUL, "", 0.), (X, "", 0.), (PLUS, "", 0.), (NUMBER, "", 3.),
(RPAR, "", 0.), (SUB, "", 0.), (NUMBER, "", 1.)
Такая программа называется сканером.
ПЕРЕРЫВ до 14:20 Закончился
2) на втором этапе компиляции работает парсер (parser, от to parse --
разбирать). Парсер выполняет синтаксический разбор и формирует
выходное представление формулы в виде обратной польской записи.
Рассмотрим реализацию сканера.
Специальная утилита Unix для создания сканеров называется LEX,
ее свободная версия flex (free lex). Используется запись лексем в виде
регулярных выражений. Мы не будем использовать lex или ее аналоги в Питоне.
Напишем сканер самостоятельно, используя модуль re для работы с
регулярными выражениями (regular expressions).
Теорема Клини:
классы автоматных и регулярных языков совпадат.
Регулярные выражения: способ построения цепочек языка.
4 операции:
1) eps образует регулярный язык ()
2) для любой буквы a язык, состоящий
из однобуквенной цепочки, регулярен a
3) произведение регулярных языков e1 e2
регулярно L1*L2 = {u1 u2 | u1 <- L1, u2 <- L2}
4) объединение двух рег. языков регулярно (e1 | e2)
L1 U L2 = { u | u <- L1 или u <- L2 }
5) итерация Клини e*
L* = { u1 u2 ... u_k, k >= 0, u_i <- L }
Как задать лексему "Вещественное число"?
1.5 100 0.5 0. 1e+30 2.5E6 1e30 1e-12
(a|b|c|d|...|z)
[a-z]
[a-zA-Z] -- любая латинская буква
[a-zA-Z_] добавили _, которое обычно трактуется как буква
Имя (идентификатор в языках программирования) -- начинается с буквы,
дальше идет любая последовательность букв и цифр
[a-zA-Z_][a-zA-Z_0-9]* -- регулярное выражение для имени.
Напришем регулярное выражение для числа:
1.5 100 0.5 0. 1e+30 2.5E6 1e30 1e-12 ---- .5
[0-9]+ любая НЕПУСТАЯ последовательность цифр, + означвает
повторение 1 или более раз.
Точка в рег. выражении означает любой символ, кроме перевода строки.
Поэтому обычную точку надо экранировать.
[0-9](\.[0-9]*)? -- числовые константы без экспоненциальной части
Вопросительный знак обозначает необязательный фрагмент
(expr)? == (expr|)
Числовые константы с необязательной экспоненциальной частью:
[0-9]+(\.[0-9]*)?((e|E)(\+|\-)?[0-9]+)? -- регулярное выражение для
числовой константы (целой или
вещественной).
Вход: "sin(1)^2 + x"
Сканер выдает: (NAME, "sin"), LPAR, (NUMBER, 1), RPAR, POW, (NUMBER 2),
PLUS, X
Парсер создает RPN:
(NUMBER, 1), (NAME, "sin"), (NUMBER 2), POW, X, PLUS
А.Пентус, М.Пентус
Математическая теория формальных языков
=================================================
27 Apr 2021
Ахо, Ульман
Формальные языки и компиляторы
Грис
Разработка компиляторов
Теория компиляции в основном основана на понятии формальной
грамматики, введенном Хомским.
Формальные грамматики были вначале придуманы для описания
естественных языков, но для сложных языков типа русского они
подходять плохо (в отличие от английского) -- плохо отражают
согласования слов, падежи, спряжения и т.п.
Язык -- подмножество множества конечных цепочек из символов алфавита
|Sigma| < inf
Элементы алфавита называюся терминалами, или обычными символами.
Для описания грамматики мы добавляем конечное множество
метасимволов, или нетерминалов
|N| < inf, Sigma не пересекается с N
Обычно элементы Sigma (терминалы) обозначаюся маленькими буквами
из начала алфавита: a, d, c, d, e, ...
Нетерминалы (метасимволы) обозначаются большими буквами:
A, B, C, R, T, ...
Любые цепочки обозначаются маленткими буквами из конца латинского
алфавита: u, v, w, x, y, z, ...
Множество нетерминалов содержит начальный нетерминал S <- N
Формальная грамматика языка L состоит из правил:
u -> v
u, v -- цепочки из терминалов и нетерминалов, причем в u содержится
хотя бы один нетерминал
Определим вывод:
один шаг вывода состоит в замене в выводимой цепочке по правилу
u -> v
подслова u на слово v. Слово v может быть пустым,
в книгах обычно пустое слово обозначается буквой eps.
Вывод начинается с начального нетерминала
и заканчивается, когда в цепочке не остается нетерминалов
(состоит только из терминалов, отсюда и происходит название
"терминал").
Пример: плохая грамматика арифметических формул
S -> a
S -> b
...
S -> z
S -> S+S
S -> S-S
S -> S*S
S -> S/S
S -> (S)
=====================================
Используем знак | (или) в правой части:
S -> a | b | c | ... | z
S -> S+S | S-S | S*S | S/S | (S)
Язык, задаваемой этой грамматикой, состоит из правильных формул,
составленных из переменных a, b, ..., z с использование скобок и
знаков арифметических операций.
Выведем формулу (a + b)*c
S => S*S => (S)*S => (S+S)*S => (a+S)*S => (a+b)*S => (a+b)*c
- - - -
Формула (a+b)*c уже не содержит нетерминалов => принадлежит языку.
Иерархия Хомского:
0) грамматики без ограничений (слева может стоять любое непустое
слово, включающее хотя бы один нетерминал).
Грамматики без ограничений обычно не рассматриваются, т.к. они
порождают класс рекурсивно перечислимых языков, для которого,
к примеру, алгоритмически неразрешима проблема принадлежности
цепочки языку.
1) контекстно зависимые грамматики (или просто контекстные),
правила имеют вид
xAy -> xuy, |u| > 0
Пример: L = { a^n b^n c^n, n >= 1 } контекстно зависимый
abc, aabbcc, aaabbbccc, aaaabbbbcccc, ...
S -> ABc | ABSc
BA -> AB
Bb -> bb
Bc -> bc
Ab -> ab
Aa -> aa
Выведем цепочку aabbcc
S => ABSc => ABABcc => AABBcc => AABbcc => AAbbcc => Aabbcc => aabbcc
- -- - -- -- --
Как избавиться от правила BA -> AB, которое не является
контекстно зависимым?
Введем еще один нетерминал R и правила
BA -> BR
BR -> AR
AR -> AB
при этим правило BA -> AB выкидываем.
Этот язык не констекстно свободный.
2) контекстно свободные грамматики: правила имеют вид
A -> u
A -- нетерминал, u -- любое слово, включая пустое.
Наша плохая грамматика арифм. формул контекстно свободная
S -> a | b | c | ... | z
S -> S+S | S-S | S*S | S/S | (S)
Почему она плохая?
Контекстно свободная грамматика называется неоднозначной,
если существую цепочки, имеющие два разных правых вывода
(или два разных левых вывода, или два разных дерева вывода).
Дерево вывода (синтаксическое дерево). Пример:
S => S*S => (S)*S => (S+S)*S => (a+S)*S => (a+b)*S => (a+b)*c
Дерево вывода
S
/|\
S * S
/|\ |
( S ) c
/|\
S + S
| |
a b
Выводимая цепочка (a+b)*s соответствует кроне дерева.
Для этой цепочки синтаксическое дерево только одно. Но рассмотрим
цепочку a + b * c, для нее можно выписать два разных дерева:
S S
/ | \ / | \
S * S S + S
/|\ | | /|\
S + S c a S * S
| | | |
a b b c
Первое дерево соответствует тому, что мы сначала складываем a и b
и затем полученную сумму умножаем на c; второе дерево соответствует
тому, что мы сначала перемножаем b и c и затем произведение
прибавляем к a (учитывает наше представление о том, что приоритет
операции * выше, чем +.
В данном случае можно изменить грамматику так, чтобы она перестала
быть неоднозначной.
Замечание: 1) все естественные языки имеют неоднозначную
грамматику:
Струи Арагвы и Куры
Узрели русские шатры.
2) существуют контекстно свободные языки, для которых не существует
однозначной грамматики (существенно неоднозначные языки).
Преобразуем грамматику арифм. формул, чтобы она стала однозначной:
учтем приоритет операций.
Рассмотрим грамматику с маркером конца строки:
(a+b)*c$, a - b*c + d$, ...
формулы заканчиваюся маркером конца строки $.
S -> F $
F -> T | F + T | F - T
T -> M | T * M | T / M
M -> E | ( F )
E -> a | b | ... | z
В этой грамматике любая формула имеет только одно дерево вывода;
например, выпишем дерево для формулы a + b * c
S
/ \
F $
/ | \
F + T
| /|\
T T * M
| | |
M M E
| | |
E E c
| |
a b
3) праволинейные языки
A -> uB
A -> v
где u, v -- терминальные цепочки.
Леволинейные языки
A -> Bu
A -> v
Теорема Клини. Следующие 4 класса языков совпадают:
1) праволинейные языки;
2) леволинейные языки;
3) автоматные языки;
4) регулярные языки.
Нетривиальная часть теоремы Клини состоит в доказательстве
импликации 3 => 4.
Лемма о разрастании (Pumping lemma).
Язык L -- контекстно свободный =>
существует целое число K такое, что:
если цепочка u из языка L достаточно длинная
|u| >= K
то u можно представить в виде
u = s x t y h
так, что
1) |x t y| <= K
2) для всякого целого i >= 1 цепочка
s x^i t y^i h
тоже принадлежит языку L.
s x t y h
s xx t yy h
s xxx t yyy h
s xxxx t yyyy h
. . .
все такие цепочки также принадлежат языку L.
Доказывается легко с помощью деревьев вывода.
Легко видеть, что язык L = { a^n b^n c^n } не удовлетворяет этому
свойству => значит, он не контекстно свободный.
ПЕРЕРЫВ 15 минут до 14:30
Задача разбора -- это задача восстановления вывода цепочки языка
для заданной терминальной цепочки (или доказательство того, что
цепочка не принадлежит языку).
Parsing, от to parse -- разбирать. Соответствующие
программы называются парсерами или синтаксическими анализаторами.
Задача разбора противоположна задаче вывода.
Дана терминальная цепочка, надо восстановить ее вывод.
Есть два класса алгоритмов разбора:
1) рекурсивный разбор, или LL(1)-разбор
L -- читаем цепочку слева направо, Left to right
L -- строим левые выводы (Left)
1 -- заглядываем на 1 символ вперед, т.е. при выборе
очередного действия
мы учитываем первый символ необработанной часто цепочки;
2) восходящий, или LR(1)-разбор, или разбор на конечном автомате
со стековой памятью (с магазинной памятью).
L -- читаем цепочку слева направо, Left to right
R -- строим правые выводы (Right)
1 -- заглядываем на 1 символ вперед, т.е. при выборе
очередного действия (shift или reduce)
мы учитываем первый символ необработанной часто цепочки;
При рекурсивном разборе дерево вывода строится сверху вниз и
слева направо. При восходящем разборе дерево вывода строится
снизу вверх. При этом второй алгоритм более сильный:
Теорема: если грамматика допускает LL(1)-разбор, то она допускает и
LR(1)-разбор, обратное, конечно же, неверно.
Утилита lex (или ее свободная версия bison) строит LR(1)-парсер
по записи контекстно свободной грамматики.
Рекурсивный парсер, в отличие от LR(1)-парсера, можно написать
вручную, практически не владея никакой теорией.
Идея рекурсивного разбора:
для каждого нетерминала реализуется функция, которая из
необработанной части входной цепочки выделяет максимальное
начало, которое выводится из этого нетерминала.
A -> aBEcD | BAca
void processA() {
t = peekToken();
if (t == a) {
skipToken(a);
processB();
processE();
skipToken(c);
processD();
} else {
processB();
processA();
skipToken(c);
skipToken(a);
}
}
Основная функция:
S -> F $
void parse() {
processF();
skipToken($);
}
Наша хорошая грамматика арифметических формул
S -> F $
F -> T | F + T | F - T
T -> M | T * M | T / M
M -> E | ( F )
E -> a | b | ... | z
не допускает рекурсивный разбор!
F -> T | F + T | F - T
Невозможно понять, какое из трех правил применимо, поскольку
необработанная часть входной цепочки во всех трех случаях
может начинаться с одного и того же терминала.
void processF() {
...
} else {
// F -> F + T
processF();
skipToken(+);
processT();
} ...
...
}
Как быть?
Есть способ преобразовать грамматику так, чтобы
она допускала левую рекурсию. Общая конструкция:
A -> A u1 | A u2 | ... | A uk | w1 | w2 | ... |wl
где все слова w_i не начинаются с нетерминала A.
Что выводится из A?
A => A u_i1 => A u_i2 u_i1 => ... => A u_is u_i_s-1... u_i1 =>
=> w_j u_is u_i_s-1... u_i1
A => Au => Auu => Auuu => ... => Auu...u => wuu...u
Вместо правил для A, перечисленных выше, используем
следующие правила:
A -> w1 ATail | w2 ATail | ... | w_l ATail
ATail -> eps | u1 ATail | u2 ATail | ... | u_k ATail
Выкидываем правила
A -> A u1 | A u2 | ... | A uk | w1 | w2 | ... |wl
заменяем их на
A -> w1 ATail | w2 ATail | ... | w_l ATail
ATail -> eps | u1 ATail | u2 ATail | ... | u_k ATail
Применим этот прием для нашей грамматики арифметических формул:
S -> F $
F -> T | F + T | F - T
T -> M | T * M | T / M
M -> E | ( F )
E -> a | b | ... | z
Получим
S -> F $
F -> T FTail
FTail -> eps | + T FTail | - T FTail
T -> M TTail
TTail -> eps | * M TTail | / M TTail
M -> E | ( F )
E -> a | b | ... | z
Рассмотрим проект lparser.py
Грамматика языка:
S -> F $
F -> T | F + T | F - T
T -> M | T * M | T / M
M -> E | M ^ E
E -> Func | NUMBER | ( F )
Func -> NAME ( F )
Добавлена операция возведения в степень ^ или **, имеющая более
высокий приоритет, чем * и /, и выполняющаяся справа
налево. Преобразуем эту грамматику к рекурсивному виду:
=========================================================
11 May 2021
Система компьютерной математики SageMath
1. Возможность выполнять символьные вычисления и получать
точный (а не приближенный!) ответ.
2. SageMath -- свободная и бесплатная система, которую развивает
весь. SageMath включает в себя большинство других математических
свободных пакетов (Maxima, Singular и др.).
3. SageMath реализован как надстройка на Python3. Объекты SageMath
реализованы как модули и классы в Python3. Практически любая
Python3-программы является корректной Sage-программой, хотя в Sage
есть и препроцессор, позволяющий записывать математические
выражения в более привычном для математиков виде. Например, в Python
операция возведения в степень записывается как **, но математики
привыкли записывать ее как ^.
x 2 - 3*x + 2 -- так удобнее, вместо
x**2 - 3*x + 2
Также можно, например, писать true или false вместо True или False.
4. Так как Sage рализован как надстройка над Python, не нужно
дополнительно учить язык для системы компьютерной математики,
плюс к тому можно использовать модули и программы Питона.
Как установить SageMath на компьютере?
В Linux устанавливается как пакет. Под MS Windows тоже можно установить
SageMath, но фактически при этом устанавливается эмулятор Linux.
Как вариант, можно на MS Windows-10 установить пакет WSL от Microsoft.
Но можно вообще ничего не устанавливать на локальном компьютере,
а работать в облаке:
http://www.sagemath.org/
Выбираем sageMathCell
переходим по ссылке, и теперь писать Sage-скрипты и испольнять их
на удаленном сервере, а результат видеть в браузере.
Пример: решим уравнение x^3 - x - 1
f(x) = x^3 - x - 1
res = solve(f(x) == 0, x)
x0 = res[0].rhs()
x1 = res[1].rhs()
x2 = res[2].rhs()
print(x0.n(), x1.n(), x2.n())
print(x2)
-0.662358978622373 - 0.562279512062301*I
-0.662358978622373 + 0.562279512062301*I
1.32471795724475
(1/18*sqrt(23)*sqrt(3) + 1/2)^(1/3) + 1/3/(1/18*sqrt(23)*sqrt(3) + 1/2)^(1/3)
Рассмотрим задачу кубической интерполяции -- построение кубического сплайна.
Построение C1-сплайна основано на следующем соображении.
Рассмотрим отрезок [x0, x1], x1 < x2 <- R
Пусть заданы значения функции в этих точках
f(x0) = y0
f(x1) = y1,
а также значения производной функции в этих точках:
f'(x0) = dy0
f'(x1) = dy1
Тогда существует единственный кубический многочлен p(x), удовлетворяющий
этим условиям:
p(x) <- R[x],
p(x0) = y0
p(x1) = y1,
p'(x0) = dy0
p'(x1) = dy1
Требуется вычислить этот многочлен.
Как это можно сделать?
Ищем p(x) в виде
y0 + c1(x-x0) + c2(x-x0)(x-x1) + c3(x-x0)^2*(x-x1)
и последовательно вычисляем коэффициенты c1, c2, c3.
Но можно просто воспользоваться SageMath.
Многочлен ищем в виде
p(x) = a0 + a1*x + a2*x^2 + a3*x^3
var("x x0 x1")
var("y0 y1 dy0 dy1")
var("a0 a1 a2 a3")
p(x) = a0 + a1*x + a2*x^2 + a3*x^3
dp(x) = derivative(p(x), x)
eq0 = (p(x=x0) == y0)
eq1 = (p(x=x1) == y1)
eq2 = (dp(x=x0) == dy0)
eq2 = (dp(x=x1) == dy1)
res = solve([eq0, eq1, eq2, eq3], [a0, a1, a2, a3])
print("a0 =", res[0].rhs())
print("a1 =", res[1].rhs())
print("a2 =", res[2].rhs())
print("a3 =", res[3].rhs())
===================================================
Система компьютерной математики SageMath
SageMath можно установить на локальном компьютере
(под Linux).
Под MS Windows 10 имеется сервис WSL -- Windows System for Linux.
В SageMath можно работать в облаке:
https://sagecell.sagemath.org/
Задача: значения неизвестной функции заданы в узлах:
(x0, y0), (x1, y1), ..., (x_n-1, y_n-1)
f(x_i) = y_i
Надо построить кубический С2-сплайн, аппроксимирующий эту
функцию
s(x):
s(x_i) = y_i,
s, s', s'' непрерывны
На каждом отрезке [x_i, x_i+1] s(x) представляет собой кубический
многочлен.
Рассмотрим упрощенный случай n = 3
Даны 3 узла
(x0, y0), (x1, y1), (x2, y2)
2 способа решения
1. Воспользоваться функций solve.
2. Решить систему линейных уравнений.
Реализуем две функции:
s = computeC2Spline(nodes)
Функция возвращает сплайн, т.е. список кубических многочленов.
Каждый многочлен задается 4-мя коэффициентами
p(x) = a0 + a1*x + a2*x^2 + a3*x^3 [a0, a1, a2, a3]
s -- список четверок из n-1 элемента, где n = len(nodes)
y = splineValue(nodes, spline, x)
Сегменты сплайна:
p0(x) = a00 + a01*x + a02*x^2 + a03*x^3
p0'(x) = a01 + 2*a02*x + 3*a03*x^2
p0''(x) = + 2*a02 + 6*a03*x
p1(x) = a10 + a11*x + a12*x^2 + a13*x^3
p1'(x) = a11 + 2*a12*x + 3*a13*x^2
p1''(x) = + 2*a12 + 6*a13*x
Система уравнений:
2 уравнения в начальном узле (x0, y0):
p0(x0) = y0
p0''(x0) = 0
4 уравнения в промежуточном узле (x1, y1):
p0(x1) = y1
p1(x1) = y1
p0'(x1) - p1'(x1) = 0
p0''(x1) - p1''(x1) = 0
2 уравнения в последнем узле (x2, y2):
p1(x2) = y2
p1''(x2) = 0
ПЕРЕРЫВ до 14:05
Алгебра матриц над полем
Какие в ней левые идеалы?
0 0 * 0
0 0 * 0 -- левый идеал (минимальный!)
0 0 * 0
0 0 * 0
Двустронних идеалов в алгебре матриц над полем нет!
Какие идеалы в Z?
Ответ -- все идеалы имею вид k*Z, k <- Z
Идеал, порожденный элементом k
I = {k*x для любых x <- R} -- главный идеал
Кольца Z и F[x] являются кольцами главных идеалов,
любой идеал порождается одним элементом.
Кольцо многочленов F[X], где |X| > 1, уже не является
кольцом главных идеалов.
Теорема о гомоморфизме
Алгебра многочленов над С (поле компл. чисел) C[X] -- свободная алгебра
A -- любая алгебра
любое отображение X --> A x_i |-> a_i
продолжается до гомоморфизма
С[X] --> A f(x1, ..., x_k) |-> f(a1, ..., a_k)
~
A = C[X]/I
f(x) <- C[x] => существует x0: f(x0) = 0
Теорема Гильберта о нулях
-------------------------
I < C[x1, x2, ..., x_n] -- собственный идеал
I не равен C[x1, x2, ..., x_n], т.е. I не содержит единицу
=> существуют числа c1, c2, ..., c_k <- C:
точка (с1, с2,..., c_k) является корнем идеала I
Пример:
x^2 + y^2 - 1
Нули этого многочлена образуют окружность радиуса 1 с центром (0, 0)
Пусть f(x) обращается в 0 на этой окружности =>
для некоторого m многочлен f(x)^m <- идеалу, порожденному многочленом
x^2 + y^2 - 1
(x+2)^3 - (x+y-1)^2 + (z-3)^2 - 5 = 0
xyz^4 - 17 x^2*y + y*z^3 - 8 = 0
(x+y+z)^3 - 1 = 0
========================================
25 May 2021
Элементы коммутативной алгебры:
теоремы Гильберта о конечности базисов и о нулях,
фактор-алгебры алгебры многочленов, идеалы в алгебре
многочленов. Базисы Грёбнера (стандартные базисы,
базисы Ширшова-Грёбнера) идеалов в алгебре многочленов
от нескольких переменных.
R = F[x0, x1, ..., x_n-1] -- алгебра многочленов над полем F.
Идеал I < R -- подалгебра, выдерживающая умножения на
на любые элементы алгебры R
для всякого u <- R, для всякого i <- I
u*i <- I
В некоммутативном случае различают двусторонние идеалы,
а также левые и правые идеалы.
В алгебре квадратных матриц порядка n двусторонних идеалов
вообще нет, а односторонним идеалом является, например,
множество матриц, у которых только элементы i-го столбца
могут быть отличны от нуля, 0 < i < n-1
0 * 0
0 * 0 -- минимальный левый идеал в алгебре матриц.
0 * 0
* * *
0 0 0 -- минимальный правый идеал
0 0 0
Теорема Гильберта о конечности базисов
Пусть I -- идеал в алгебре многочленов R = F[x0, x1, ..., x_n-1].
Тогда I порождается конечным числом элементов
I = Ideal(f1, f2, ..., fk), f_i <- R
Определение
Кольцо (алгебра) называется нётеровым, если в нем любой
идеал порождается конечным числом элементов.
Алгебра многочленов от конечного числа переменных над полем
нётерова.
Доказательство следует из утверждения:
кольцо многочленов от одной переменной над нётеровым
кольцом нётерово.
R нётерово => R[x] нётерово.
R[x, y] = R[x][y]
Теорема Гильберта о нулях.
Упрощенная формулировка.
Пусть I < C[x0, x1, ..., x_n-1] -- идеал в алгебре многочленов
над полем комплексных чисел
или любым другим алгебраически замкнутым
полем.
Пусть I != C[x0, x1, ..., x_n-1], т.е. I -- собственный идеал
Это эквивалентно тому, что I не содержит единицу.
Тогда идеал I имеет корень, т.е. существует набор элементов поля C
(c0, c1, ..., c_n-1):
для любого многочлена f(x0, x1, ..., x_n-1) <- I
f(c0, c1, ..., c_n-1) = 0.
I < M -- макс. идеал
C[x0, x1, ..., x_n-1]/M -- поле
По теореме Гаусса-Зарисского это есть конечное расширение поля C =>
это есть просто C
~
C[x0, x1, ..., x_n-1]/M = C
- - -
тогда (x0, x1, ..., x_n-1) является корнем любого многочлена из I
Стандартная формулировка Теоремы Гильберта о нулях
Пусть I < C[x0, x1, ..., x_n-1] -- идеал в алгебре многочленов
над полем комплексных чисел
или любым другим алгебраически замкнутым
полем.
f(x) <- C[x0, x1, ..., x_n-1] обращается в ноль на нулях идеала I <=>
<=> f^m <- I для некорого натурального m.
Импликация <= очевидна.
Обратная импликация доказывается с помощью так называемого "трюка Рабиновича".
Пусть f обращается в 0 на корнях идеала I.
По теореме Гильберта о конечности базисов I = Ideal(g1, g2, ..., g_k)
Рассмотрим идеал J в алгебре C[x0, x1, ..., x_n-1, y], порожденный
элементами
J = Ideal(g1, g2, ..., g_k, 1 - y*f)
У него нет корней, поскольку, если g1, g2, ..., g_k обращаются в 0,
то и f также обращается в 0 на том же наборе чисел => многочлен
1 - y*f равен 1 на этом наборе.
Идеал J совпадает со всем кольцом => J содержит единицу
1 = u1*g1 + ... + u_k*g_k + v*(1 - y*f)
Рассмотрим это равенство в поле рациональных функций от x0, ..., x_n-1, y.
Подставим вместо y = 1/f
1 = u1'*g1 + u2'*g1 + ... + u_k'*g_k
В знаменателях стоят степени f. Домножим на f^m, где m -- максимальная
степень y, входящая в выражения u_i.
f^m = w1*g1 + w2*g2 + ... + w_k*g_k
w_i = u_i'*f^m -- это многочлен от x0, x1, ..., x_n-1
Получили f^m <- I
Определение:
радикалом идеала I называется множество многочленов f таких, что
для некоторого m = m(f)
f^m <- I
Из теоремы Гильберта о нулях следует, что радикал идеала I
совпадает с множеством многочленов, обращающихся в 0 на
корнях идеала I.
Теорема:
многочлен f обращается в ноль на корнях идеала I = Ideal(g1, g2, ..., g_k) <=>
<=> Ideal(g1, g2, ..., g_k, 1 - y*f) содержит единицу.
f принадлежит радикалу идеала I = Ideal(g1, g2, ..., g_k) <=>
<=> Ideal(g1, g2, ..., g_k, 1 - y*f) содержит единицу.
Насколько конструктивна эта теория? Всегда ли можно, например,
проверить, лежит ли данный многочлен в идеале I = Ideal(g1, g2, ..., g_k)?
Другими словами, разрешима ли алгоритмически проблема равенства
элементов в фактор-алгебре F[x0, x1, ..., x_n-1]/I ?
Интересно, что в группах, в некоммутативных алгебрах, в полугруппах и т.п.
проблема равенства алгоритмически неразрешима.
Удивительно, что в коммутативных алгебрах большая часть проблем имеет
алгоритмическое решение -- например, алгоритмически разрешима проблема
равенства элементов в фактор-алгебре === проблема принадлежности
элемента идеалу в алгебре многочленов.
Алгоритмы основаны на построении базиса Грёбнера идеала.
Идея:
пусть один из многочленов, порождающих идеал I, имеет вид
u - w, где u -- старший моном многочлена (т.к. рассматриваем
многочлены над полем, можно считать, что
старший коэффициент равен 1).
Это значит, что в фактор-алгебре R/I, где R = F[x0, x1, ..., x_n-1],
справедливо равенство
u == w (mod I)
Любое вхождение монома u в некоторый моном можно заменить на w
f(x, y) = x^2 - 2*x*y
(2, 0) (1, 1)
x^2 == 2*x*y (mod I)
многочлен x^3 - y - 1 == x*(2*x*y) - y - 1 = 2*x^2*y - y - 1 ==
== 4*x*y^2 - y - 1
Пусть I = Ideal(g1, g2, ..., g_k)
g_i = u_i - w_i, где u_i -- старший моном
Редукция -- это замена монома u_i на w_i
a u_i b --> a w_i b
Система порождающих g1, g2, ..., gk называется базисом Грёбнера идеала I,
если любой элемент из идеала I приводится редукциями к нулю.
Алгоритм проверки равенства элементов в фактор-алгебре:
два элемента равны, если равны их нормальные формы.
Предложение:
нормальная форма элемента определена однозначно <=>
система порождающих идеала является базисом Грёбнера.
Теорема. Для любого идеала в алгебре многочленов над любым полем
можно построить конечный базис Грёбнера (существует алгоритм построения).
Идея:
g1 = u1 - w1 xy - f1 HOK(u1, u2) = u1*u2
g2 = u2 - w2 zt - f2
xy - f1 HOK(xy, yz) = xyz
yz - f2
xyz -> f1*z
xyz -> x*f2
s-многочлен -- это f1*z - x*f2
U = HOK(u1, u2) s-многочлен, это
(U/u1)*w1 - (U/u2)*w2
Если он не приводится к нулю, то добавляем его к системе порождающих идеала.
Можно доказать, что этот процесс рано или поздно оборвется.
ПЕРЕРЫВ до 14:30
Продолжение
Алгоритмические задачи, решаемые с помощью базисов Грёбнера
1. Определить, совместима ли система алгебраических уравнений над C.
Достаточно проверить, что 1 не содержится в идеале
<=> базис Грёбнера идеала отличен от [ 1 ]
2. Даны 2 системы алгебраических уравнений над C.
Определить, эквивалентны ли они,
т.е. имеют одинаковые множества решений (задают ли они одно и то же
аффинное многообразие)
I1 = Ideal(f1, f2, ..., f_k)
I2 = Ideal(g1, g2, ..., g_l)
Надо проверить, что g_i обращаются в ноль на корнях идеала I1 и, обратно,
f_j обращаются в ноль на нулях идеала I2.
То есть g_i <- Rad(I1)
f_j <- Rad(I2)
g <- Rad(Ideal(f1, f2, ..., f_k)) <=>
Ideal(f1, f2, ..., f_k, 1 - y*g) содержит единицу <=>
Базис Грёбнера этого идеала Ideal(f1, f2, ..., f_k, 1 - y*g)
состоит из одной единицы.
2. Дана система алгебраических уравнений над C. Проверить, имеет ли она
конечное непустое множество решений.
Отметим, что совместимость системы мы уже умеем проверять.
Теорема. Система имеет конечное множество решений <=>
Базис Грёбнера идеала левых частей уравнений при использовании
чисто лексикографического порядка имеет для каждой переменной
x_i многочлен со старшим мономом, равным степени x_i.
Чисто лексикографический порядок на мономах.
Моном = произведение переменных в некоторых степенях
x^3*y*z^4 (3, 1, 4)
y^2*z (0, 2, 1)
Получается x^2 > x*y^10*z^1000
По умолчанию Sage использует степенно-лексикографический порядок
"deglex": сначала сравниваем суммарные степени мономов и только если
они равны, сравниваем лексикографически степени по каждой переменной.
Пусть выполняется это условие. Докажем, что множество решений конечно.
Рассмотрим переменную с максимальным номером. В базисе Грёбнера есть
многочлен со старшим мономом = степени этой переменной. Этот многочлен не
содержит мономов, зависящих от других переменных => это многочлен от
одной переменной => имеет конечное множество корней.
Последовательно подставляя эти корни во все оставшиеся уравнений
и отбрасывая рассмотренное уравнение, мы получаем
систему от меньшего набора переменных, которая удовлетворяет тому же условию
=> по всем остальным переменным тоже имеется лишь конечное множество решений.
Алгоритм
Вычисляем базис Грёбнера идеала левых частей уравнений,
важно, что мы используем чисто лексикографический порядок.
Затем проверяем, что для каждой переменной в базисе Гребнера
имеется многочлен со старшим мономом, равным степени этой переменной.