Sep 14, 2022
Конструкция для задания списка List Comprehension
(умное задание списка, математическое задание списка)
В математике мы задаем множество элементов:
S = { x: x <- [1..29], gcd(x, 30) = 1 }
В Питоне:
s = [ x for x in range(1, 30) if gcd(x, 30) == 1 ]
В С++ цикл for моделируется с помощью while:
for (int i = 1; i < 30; ++i) {
. . .
}
эквивалентно
int i = 1;
while (i < 30) {
. . .
++i;
}
В Питоне цикл for не сводится к while
for x in s:
# Тело цикла
В Питоне3 range(a, b) -- полуинтервал целых чисел a <= x < b
(в Питоне 2 это был просто список)
Задача:
получить список всех пифагоровых троек
(a, b, c), a, b, c <- N, a^2 + b^2 = c^2, a <= b <= c <= 100,
gcd(a, b) = 1
s = [(a, b, c) for a in range(1, 101)
for b in range(a, 101)
for c in range(b, 101)
if gcd(a, b) == 1 and a**2 + b**2 == c**2]
Задача 2:
решить кубическое уравнение.
a*x^3 + b*x^2 + c*x + d = 0
Можно считать a = 1
x^3 + a*x^2 + b*x + c = 0
Можно избавиться от квадратичного члена, применив замену
типа сдвига x' = x + s
(x' + s)^3 + a*(x' + s)^2 + ... =
x'^3 + 3*x'^2*s + ... + a*x'^2 + ... =
x'^2 + (3*s + a)*x'^2 + ...
Достаточно взять
s = -a/3
Уравнение сводится к виду
x^3 + a*x + b = 0, считаем, что b != 0
Достаточно найти хотя бы 1 корень. Разделив на (x - r), где
r -- корень, мы получим квадратное уравнение.
Основная идея: сделаем замену
x = u + v
(u + v)^3 + a*(u + v) + b = 0
u^3 + 3*u^2*v + 3*u*v^2 + v^3 + a(u + v) + b = 0
u^3 + v^3 + b + 3*u*v*(u + v) + a(u + v) = 0
u^3 + v^3 + b + (3*u*v + a)(u + v) = 0
------------- ==================
Получаем систему из 2-х уравнений с 2-мя неизвестными
u^3 + v^3 + b = 0
(3*u*v + a)(u + v) = 0
Так как x = 0 -- не решение, то (u + v) = x != 0, сократим на (u + v)
u^3 + v^3 + b = 0
3*u*v + a = 0
Выразим v из 2-го уравнения и подставим в первое уравнение
v = -a/(3*u)
u^3 - a^3/(27*u^3) + b = 0
Сделаем замену
y = u^3
Получим
y - a^3/(27*y) + b = 0
27*y^2 + 27*y*b - a^3 = 0
Получили квадратное уравнение
d = (27*b)^2 + 4*27*a^3
y = (-27*b + d^(1/2))/54
Дома:
решить кубическое уравнение (получить список из 3-х решений)
http://mech.math.msu.su/~vvb/
shokurov.anton.v@yandex.ru -- лектор
21 Sep 2022
Вернемся к Питону.
Как устроен цикл "для каждого элемента структуры данных"?
for x in range(1, 21):
print((x, x**2, x**3))
Как моделируется цикл "для каждого" (цикл for) в разных
языках программирования?
В стандартной библиотеке С++ (STL -- Standard Library)
для этого используются итераторы:
class S; // структура данных
// вложенный класс class S::iterator
// class S::const_iterator
// методы begin() -- возвращает итератор,
// указывающий на начало
// cbegin() -- возвращает константный итератор
//
// end() -- возвращает итератор,
// указывающий на за-конец структуры данных
// cend()
// Итератор как бы обобщает указатель на элементы структуры данных
// Если i -- итератор, то *i -- это тот элемент, на который он ссылается
// С итератором можно выполнять действия:
// *i
// ++i
// i == j, i != j
S a;
. . .
S::iterator i = a.begin();
while (i != a.end()) {
// Действие (a)
. . .
++i;
}
Как цикл for моделируется в Питоне?
Структура данных должна иметь метод iter(s)
(на самом деле имя метода __iter__(self))
который возвращает объект i, у которого есть метод next(i)
Этот метод возвращает очередной элемент структуры данных
и переходит к следующему элементу.
Если элементы исчерпаны, то возбуждается исключение типа StopIteration
Про объекты:
Питон -- объектно ориентированный язык, переменные не имеют типов
и хранят объектные ссылки. Объектная ссылка -- это не адрес объекта,
потому что объекты могут перемещаться по памяти.
Память представляет собой "контролируемую кучу" (контролируемая
динамическая память). Постоянно работает (параллельно) процесс,
который называется "сборщик мусора" (garbage collector, или просто gc).
Он удаляет объекты, которые никому не нужны (на которые нет ссылок)
и переписывает объекты в памяти вплотную друг к другу (дефрагментация).
Идентификатор объекта, объектная ссылка, handle -- информация,
отождествляющая объект.
MAC Handle == pointer to pointer A** h;
Пример:
получить список номеров всех счастливых билетов длины n, n -- четное число.
Счастливый билет: сумма цифр первой половины номера билета ==
сумме цифр второй половины, например,
142511
Необязательная задача:
найти число счастливых билетов длины n (n -- четное число).
Помимо замечательно цикла for, в Питоне есть генераторы
(генераторные функции).
Генераторная функция -- это функция, которая выдает серию значений.
Отличается от обычной функции тем, что вместо
return value
очередное значение возвращается с помощью
yield value
При этом кадр локальных переменных запоминается, и при следующем вызове
функция продолжаетработу с прерванного места.
Пример:
как перебрать все целочисленные точки первого квадранта плоскости?
(0, 0)
(1, 0) (0, 1)
(2, 0) (1, 1) (0, 2)
(3, 0) (2, 1) (1, 2) (0, 3)
. . .
Помим генераторных функций, есть также генераторные выражения:
g = ((x, x**2) for x in range(1, 21))
Применение теории чисел в программировании:
криптография, электронная подпись, криптовалюты,
электронное голосование...
Устройство Энигмы
s1 -- перестановка первого ротора,
s2 -- второго
s3 -- третий
t -- перестановка рефлектора
t представима в виде произведения независимых транспозиций
t*t = 1
Для фиксированного положения роторов применяется перестановка
s1 s2 s3 t (s1)^-1 (s2)^-1 (s3)^-1
После каждой буквы первый ротор поворачивался на одну позицию (1/26 оборота)
Второй ротор поворачивался на одну позицию после того, как первый совершал
полный оборот,
третий -- после полного оборота второго ротора.
Польские математики (Решевский и ко.)
Теорема о том, что перестановка раскладывается в произведение
независимых циклов. После окупации Польши материалы были
переданы в Англию, где работами руководил Алан Тьюринг.
3 главных алгоритма теории целых чисел:
-- алгоритм Евклида
-- расширенный алгоритм Евклида
-- быстрое возведение в степень
Теорема:
m, n -- целые числа, m != 0 или n != 0
d = gcd(m, n) -- их наибольший общий делитель
Тогда найдутся целые числа u, v:
d = u*m + v*n
Расширенный алгоритм Евклида:
Дано: целые m, n, m != 0 или n != 0
Надо: найти d = gcd(m, n) и представление d в виде
d = u*m + v*n
Схема построения цикла с помощью инварианта.
Инвариант -- это соотношение между изменяющимися переменными,
которое сохраняется в теле цикла.
Пусть надо найти x <- I пересечение Q
Т -- отображение, сохраняющее инвариант
T: I\Q --> I
x = x0
утв: I(x)
цикл пока не Q(x)
| инвариант: I(x)
| x = T(x)
конец цикла
утв: Q(x) и I(x)
x = (a, b)
x0 = (m, n)
Inv: gcd(a, b) == gcd(m, n)
T: (a, b) --> (b, r), r = a%b
Q(a, b): b == 0 тогда gcd(a, 0) = a
Идея расширенного алг. Евклида:
выполняем обычный алгоритм Евклида, модифицируя пару (a, b),
при этом мы храним представление a в виде линейной комбинации исходных
чисел m, n, а также представление b.
Инвариант:
gcd(a, b) = gcd(m, n)
a = u1 m + v1 n
b = u2 m + v2 n
(a, b) <-- (b, r) a = q b + r, |r| < |b|
u1 <-- u2 r = a - q b = (u1 m + v1 n) - q (u2 m + v2 n) =
v1 <-- v2 = (u1 - q u2)m + (v1 - q v2)n
u2 <-- u1 - q u2
v2 <-- v1 - q v2
Задача: даны числа m, n, gcd(m, n) = 1
Найти обратный элемент к n в кольце Z_m
Применяем расширенный алгоритм Евклида, находим
1 = u*m + v*n
1 == v*n (mod m)
Искомый обратный элемент равен v.
https://vimeo.com/752054421
28 Sep 2022
Домашние задания, даны на предыдущих семинарах:
1) решить кубическое уравнение общего вида:
a*x^3 + b*x^2 + c*x + d = 0, a, b, c, d <- C (комплексные числа)
Написать функцию, которая возвращает список из трех корней уравнения;
2) (необязательная). Для четного числа n вычислить количество счастливых
билетов длины n, для хотя бы n <= 16.
3) обязательная! Задача на тему теории чисел (один из 5-ти вариантов).
Deadline -- середина октября.
Теория чисел
3 базовых алгоритма:
1) gcd(m, n) -- наибольший общий делитель (алг. Евклида);
2) extEuclidAlg(m, n) -- расширенный алг. Евклида
Найти d = gcd(m, n) и его выражение в виде линейной комбинации m, n:
d = u*m + v*n
3) алгоритм быстрого возведения в степень
powmod(x, n, m) возвращает x^n (mod m).
Кольца вычетов Z_m часто называют арифметическими кольцами.
Именно они используются для практических вычислений в криптографии.
В кольце Z числа растут очень быстро и не помещаются в памяти компьютера.
Важность расширенного алгоритма Евклида в том, что он позволяет вычислить
(быстро!) обратный элемент в кольце вычетов Z_m
y = inverse(x, m), gcd(x, m) = 1
y*x == 1 (mod m)
Применяем расш. алг. Евклида:
1 = u*x + v*m
1 == u*x (mod m)
y = u
Рассмотрим алгоритм быстрого возведения в степень в кольце Z_m
Дано: a <- Z_m, n <- Z, n >= 0
Надо: вычислить a^n (mod m)
Применим схему построения цикла с помощью инварианта
p*b^k == a^n (mod m)
В цикле меняются переменные p, b, k. Если k = 0, то ответ будет в
переменной p
p*b^k == p == a^n
Начальные значения: p = 1, b = a, k = n, для этих значений
инвариант выполняется.
Тело цикла: цель -- уменьшать k, сохраняя инвариант.
если k четное => k = k/2, b = b^2
k нечетное => k = k - 1, p = p*b
Вернемся к математике
Важнейшие теоремы
1. Малая теорема Ферма
Пусть p -- простое число, b != 0 (mod p)
Тогда b^(p-1) == 1 (mod p)
Следствие теоремы Лагранжа: порядок элемента конечной группы
делит порядок группы. Группа обратимых элементов кольца Z_p
U_p имеет порядок p-1.
Если отказаться от условия b != 0 (mod p)
(т.е. b -- любое число), то справедливо
b^p == b (mod p)
Более того, для любого b и любого целого s >= 0 выполняется
b^(s*(p-1) + 1) == b (mod p)
Нам понадобится именно эта форма Малой теоремы Ферма.
2. Китайская теорема об остатках
Пусть m = m1 m2 ... mk, где числа m_i взаимно простые:
gcd(m)i, m_j) = 1 при i != j
Тогда имеет место изоморфизм колец вычетов:
~
Z_m = Z_m1 + Z_m2 + ... + Z_mk
x == (x1, x2, ..., x_k), x_i == x (mod m_i)
Операции с этими строками производятся покомпонентно.
Пример применения Китайской теоремы об остатках.
Сколько квадратных корней из 1 в кольце Z_105?
105 = 3*5*7
x == (x1, x2, x3) x1 == x (mod 3),
x2 == x (mod 5),
x3 == x (mod 7)
x_i^2 == 1 (mod m_i)
Кольца Z_3, Z_5, Z_7 являются полями
В любом поле не больше двух корней из 1 (теорема Безу: многочлен степени n
имеет не больше чем n корней, корень из 1 -- это корень многочлена
x^2 - 1. => Квадратные корни из 1 -- это числа 1 и -1.
Поэтому всего 8 корней:
(1, 1, 1)
(-1, -1, -1)
(-1, 1, 1)
(1, -1, 1),
. . .
Как найти эти элементы?
Оказывается, имеется эффективный алгоритм для вычисления обратного
изоморфизма:
даны произвольные числа x1, x2, ..., x_k и попарно взаимно простые числа
получить x: x == x_i (mod m_i)
Доказательство Китайской теоремы об остатках.
Воспользуемся теоремой о гомоморфизме.
f: R --> H -- гомоморфизм
f: R --> Im(f) -- эпиморфизм
Тогда имеет место изоморфизм
~
R/ker f = Im(f)
причем диаграмма коммутативна:
f
R --> Im(f)
pi | /
| / phi
V V
R/ker f
phi*f = pi
Для доказательства Китайской теоремы об остатках
рассмотрим гомоморфизм
f: Z --> Z_m1 + Z_m2 + ... + Z_mk
f: x --> (x1, x2, ..., x_k) x_i == x (mod m_i)
Ядро этого гомоморфизма = множество чисел, которые делятся на все m_i
Поскольку m_i взаимно простые, то это все числа, кратные m, т.е.
это идеал mZ
~
Im(f) = Z/mZ = Z_m
Покажем, что Im(f) совпадает с Z_m1 + Z_m2 + ... + Z_mk,
т.е. f -- эпиморфизм. Это следует из того, что |Im(f)| = |Z_m| = m,
но число элементов прямой суммы также равно m:
|Z_m1 + Z_m2 + ... + Z_mk| = m
Следовательно, Im(f) = Z_m1 + Z_m2 + ... + Z_mk. Теорема доказана.
Вернемся к алгоритму Китайской теоремы об остатках.
Даны произвольные числа x1, x2, ..., x_k и
попарно взаимно простые числа m1, m2, ..., m_k.
Получить x: x == x_i (mod m_i)
Пример: найти x такой, что
x == 1 (mod 3)
x == -1 (mod 5)
x == 1 (mod 7)
(1, -1, 1) -- корень из 1 в Z_105
x == x1 (mod m1) => x = x1 + s*m1 для любого s
Хотим исправить x так, чтобы выполнялось и равенство
x == x1 (mod m1)
x == x2 (mod m2)
Подберем s так, чтобы выполнялось второе сравнение
x = x1 + s*m1 == x2 (mod m2)
Решаем уравнение
x1 + s*m1 == x2 (mod m2)
s*m1 == x2 - x1 (mod m2)
так как gcd(m1, m2) = 1, то m1 обратимо в Z_m2
Вычислим обратный элемент
u*m1 == 1 (mod m2)
Домножим уравнение на u:
s*m1*u == u*(x2 - x1) (mod m2)
s == u*(x2 - x1) (mod m2)
Для третьего уравнения ищем x' в виде
x' = x + s*(m1*m2)
x' == x3 (mod m3)
Пусть u == (m1*m2)^(-1) (mod m3) -- обратный элемент к m1*m2 в Z_m3
u*(m1*m2) == 1 (mod m3)
Уравнение записывается в виде
x + s*(m1*m2) == x3 (mod m3)
s*(m1*m2) == x3 - x (mod m3)
Домножив на u, получим
s == u*(x3 - x) (mod m3)
Подставив s в выражение x' = x + s*(m1*m2), получим
x' == x1 (mod m1)
x' == x2 (mod m2)
x' == x3 (mod m3)
x'' = x' + s*(m1*m2*m3)
x'' == x4 (mod m4)
. . .
ПЕРЕРЫВ до 11:05
Теорема Эйлера.
Обозначим через phi(m) количество обратимых элементов в
кольце вычетов Z_m, или это просто порядок группы обратимых элементов
U_m кольца Z_m.
Тогда имеет место теорема Эйлера
Пусть x -- обратимый элемент кольца Z_m, т.е. gcd(x, m) = 1.
Тогда x^phi(m) == 1 (mod m).
Давайте вычислим функцию Эйлера. Разложим число m на простые множители:
m = p1^e1 p2^e2 ... p_k^e_k
Обозначим p_i^e_i = m_i
m = m1 m2 ... m_k
Изоморфизм колец в Китайской теореме об остатках индуцирует изоморфизм
групп обратимых элементов этих колец:
~
U_m = U_m1 * U_m2 * ... * U_mk (прямое произведение групп)
x = (x1, x2, ..., x_k) x_i == x (mod m_i)
|U_m| = |U_m1| * |U_m2| * ... * |U_mk|
phi(m) = phi(m1)*phi(m2)*...*phi(m_k)
Мы свели вычисление функции Эйлера для произвольного числа m
к ее вычислению для степени простого числа p^e
То есть надо посчитать число обратимых элементов в кольце Z_p^e
Легче сосчитать число НЕОБРАТИМЫХ элементов: это все элементы из отрезка
0, 1, 2, ..., p^e-1, которые кратны p, их всего
p^e/p = p^(e-1)
Из числа всех элементов вычтем число необратимых элементов,
получим число обратимых, т.е. функцию Эйлера:
phi(p^e) = p^e - p^(e-1) = p^(e-1) * (p - 1)
Пример: phi(27) = phi(3^3) = 3^2*(3-1) = 9*2 = 18
Формула Эйлера в общем виде:
phi(p1^e1 * p2^e2 * ... * p_k^e_k) =
= (p1 - 1)*p1^(e1-1) * (p2 - 1)*p2^(e2-1) * ... * (p_k - 1)*p_k^(e_k-1)
Для вычисления функции Эйлера требуется разложить число на простые множители.
Обратное тоже верно: если известна фунция Эйлера phi(m), то мы можем
разложить число m на множители (эти две задачи эквивалентны по сложности).
Частный случай теоремы Эйлера.
Пусть число m свободно от квадратов, т.е.
m = p1 * p2 * p3 * ... * p_k, где p_i -- простые числа.
phi(m) = (p1 - 1)(p2 - 1)(p3 - 1)...(p_k - 1)
Классическая формулировка теоремы Эйлера:
для всякого основания b, взаимно простого с m, выполняется сравнение
b^phi(m) == 1 (mod m)
Для числа m, свободного от квадратов, можно отказаться от условия взаимной
простоты. Для ПРОИЗВОЛЬНОГО b и произвольного целого s >= 0
выполняется сравнение
b^() + 1) == b (mod m)
Доказательство: применим Китайскую теорему об остатках и
вторую форму Малой теоремы Ферма.
~
Z_m = Z_p1 + Z_p2 + ... Z_pk
b = (b1, b2, ..., b_k)
phi(m) = (p1 - 1)(p2 - 1) ... (pk - 1)
b_i ^ (s*phi(m) + 1) == b_i по Малой теореме Ферма, поскольку
phi(m) кратно (p_i - 1).
---------------------------------------------------------
Сзема кодирования с открытым ключом RSA
Эта схема названа в честь ее авторов, профессоров MIT (Boston):
Rumley, Shamir, Adleman
Классическая криптография:
-- шифр Цезаря -- циклический сдвиг букв в алфавите на 3 позиции
A B C D E F G ... X Y Z
A --> D
B --> E
...
X --> A
Y --> B
Z --> C
-- подстановочный шифр: биективное отображение алфавита в себя
(перестановка букв).
Любой подстановочный шифр разгадывается с помощью метода
частотного анализа. В любом естественном языке частоты разных
букв сильно различаются.
-- шифр Виженера -- шифр с перекрытием. Мы имеем ключевую фразу,
каждая буква которой задает циклический сдвиг в алфавите
в соответствии с ее позицией. Буква А задает сдвиг на 1 позицию,
буква Б на две позиции, буква Я на 33 позиции.
Пусть, например, ключевое слово ХОЛМС. Зашифруем сообщение
ИЛСИ, ПРИХОДИ
ХОЛМ СХОЛМСХ
ИЛСИ, ПРИХОДИ
ЯЫ... 1 2 3
123456789012345678901234567890123
абвгдеёжзийклмнопрстуфхцчшщъыьэюя
Х О Л М С
И --> Я (сдвиг на 23 позиции, соотв. букве Х)
Л --> Ы (сдвиг на 16 позиций, соотв. букве О)
Виженер -- придворный короля Луи XI (???)
На протяжении веков шифр Виженера считался неразгадываемым, пока
в середине XIX века майор прусской армии Фридрих Касицкий не предложил
метод разгадывания шифра Виженера на основе нахождения биграмм,
расположенных на расстояниях, кратных длине ключевой фразы.
В начале XX века Фридман нашел метод определения длины ключевой фразы,
основанный на индексе совпадения.
НапротяжениивековшифрВиженерасчиталсянеразгадываемым.
НапротяжениивековшифрВиженерасчиталсянеразгадываемым.
= = =
Сравнили 47 букв и нашли 3 совпадения.
Индекс совпадения 3/47 = 0.064
В случайном наборе символов индекс совпадения = 1/33 = 0.030
в два раза меньше!
При шифре Виженера при сдвиге, кратном длине ключевой фразы,
совпадения в исходном и зашифрованном тексте совпадают.
Следовательно, индекс совпадения резко возрастает, когда сдвиг кратен
длине ключевой фразы => так мы находим длину ключевой фразы.
Вывод: нельзя использорвать шифры, основанные на заменах букв!
Все они ненадежны.
Современные шифры используют блочное шифрование.
Сообщение записывается в двоичном виде и разбивается на блоки
фиксированной длины, каждый блок шифруется отдельно.
Плюс блок можно трактовать как двоичную запись целого числа,
или элемента кольца вычетов. Таким образом, шифрование
сводится к инъективному отображению колец вычетов:
E: Z_m --> Z_n -- инъективное отображение
или даже при m = n
E: Z_m --> Z_m -- биекция.
Дешифровка
D: Z_m --> Z_m
D E = 1 (тождественное отображение).
y = E(x)
D(y) = D(E(x)) = x.
--------------------------------------
5 Oct 2022
Применение теории чисел в криптографии
Схема кодирования с открытым ключом RSA
(названа в честь ее авторов Rumley, Shamir, Adleman, MIT).
В прошлый раз мы рассмотрели
1) китайскую теорему об остатках
m = m1 * m2 * ... * m_k, gcd(m_i, m_j) = 1 при i != j
~
Z_m = Z_m1 + Z_m2 + ... + Z_m_k
x = (x1, x2, ..., x_k), x_i == x (mod m_i)
m = p1^e1 * p2^e2 * ... * p_k^e_k
Важен частный случай, когда число m свободно от квадратов:
m = p1 * p2 * ... * p_k
2) Малая теорема Ферма:
p -- простое число
b: gcd(b, p) = 1
Тогда
b^(p-1) == 1 (mod p)
Для произвольного b и s >= 0 справедливо
b^((p-1)*s + 1) == b (mod p)
3) теорема Эйлера
phi(m) = порядок группы обратимых элементов кольца Z_m
Если m = p1^e1 * p2^e2 * ... * p_k^e_k, то
phi(m) = (p1-1)*p1^(e1-1) * ... * (p_k-1)*p_k^(e_k-1)
В частном случае, когда m свободно от квадратов
phi(m) = (p1 - 1) * (p2 - 1) * ... * (p_k-1)
Классическая формулировка теоремы Эйлера:
b: gcd(m, b) = 1 (=> b обратимо в Z_m)
Тогда b^phi(m) == 1 (mod m)
Нам нужна более общая формулировка
Пусть число m свободно от квадратов, m = p1*p2*...*p_k
Тогда для любого b и любого s >= 0 выполняется сравнение
b^(phi(m)*s + 1) == b (mod m)
----------------------------------------------------------
Схема кодирования RSA
E: Z_m --> Z_n, n >= m -- кодирующее отображение
D: Z_n --> Z_m, D E = 1 для всякого x D(E(x)) = x
E -- Encoding, D -- Decoding
Особенность схемы RSA в том, что m = n:
E: Z_m --> Z_m,
D: Z_m --> Z_m, D E = E D = 1
Построение:
Возьмем 2 больших простых числа p, q
m = p*q
Число m открыто, но разложение m на множители секретно.
phi(m) = (p - 1)*(q - 1) -- секретная
Возьмем любой элемент e, обратимый в кольце Z_phi(m)
Число e открыто.
Отображение E (кодирование) -- это просто возведение в степень e
E: Z_m --> Z_m
E: x --> x^e (mod m)
Открытый ключ -- это пара (m, e).
Для построения обратного отображения D надо вычислить обратный
к e элемент в кольце вычетов Z_phi(m):
d: d*e == 1 (mod phi(m))
Элемент d секретный.
D: Z_m --> Z_m
D: y --> y^d (mod m)
Покажем, что это два взаимно обратных отображения
D E = 1
Из сравнения
d*e == 1 (mod phi(m))
d*e = 1 + s*phi(m), s <- Z, z >= 0
Надо показать, что для произвольного x D(E(x) = x
D(E(x)) == D(x^e) == (x^e)^d (mod m) == x^(e*d) ==
== x^(1 + s*phi(m)) == x (mod m) по теореме Эйлера.
Суммируем сказанное:
m = p*q, p, q -- простые (например, 200-значные десятичные)
e -- обратимый элемент в кольце Z_phi(m)
Открытый ключ: (m, e)
d -- обратный к e элемент в кольце Z_phi(m)
d*e == 1 (mod phi(m))
Секретный ключ: (m, d)
Отображение шифрования: x --> y == x^e (mod m)
Отображение дешифровки: y --> x == y^d (mod m)
В результате мы получаем секретный канал от всех к одному.
Любой может отправить секретное сообщение, используя открытый ключ,
но расшифровать его может только обладатель секретного ключа.
PGP -- Pretty Good Privacy.
Схему RSA можно применять в обратном направлении, поскольку не только
D E = 1, но и E D = 1.
Получаем секретный канал от одного ко всем.
Это позволяет удостоверить личность отправителя сообщения.
s -- сообщение, которое я хочу отправить всем, но так,
чтобы все были уверены, что автор сообщения -- именно я.
t = D(s)
Любой может его прочитать,применив мою открытую процедуру E:
E(t) = E(D(s)) = s
Но сгенерировать сообщение t может только обладатель секретной процедуры.
Другое применение -- электронная подпись.
Я хочу послать открытый документ и к нему добавить подпись,
удостоверяющую то, что именно я -- автор этого документа.
Пусть T -- это документ, h = H(T) -- хеш-функция, вычисленная по документу.
Другие характеристики документа: l = len(T) -- его длиня в байтах,
d -- дата создания документа,
n -- номер транзакции, ...
Набор (h, l, d, n) шифруем, используя секретную процедуру.
Получаем электронную подпись:
S = D( (h, l, d, n) )
Из электронгной подписи восстанавливаются параметры документа:
(h, l, d, n) = E(S)
И сравниваем h, l с хеш-функцией, вычисленной по документу, и его длиной.
Если они совпадают, то документ считается подлинным.
Какие задачи возникают в связи со схемой RSA
1. Генерация пары ключей. Для этого надо уметь генерировать
большие простые числа p, q.
Простых чисел очень много!
Количество простых чисел, не превосходящих n, по формуле Чебышева
pi(n) ~ n/ln(n)
в отрезке длины ln(n) содержится в среднем одно простое число.
Если мы имеем тест простоты, то мы довольно быстро найдем простое
число, перебирая случайные числа нужной длины.
Можно для этого использовать вероятностный тест простоты.
Отметим, что существует способ генерации случайных простых чисел
за полиномиальное время с предъявлением сертификата простоты
(это некоторый несложный способ, удостоверяющий, что число простое --
доказательство, которое легко проверяется).
Совсем недавно (2005 ???) был найден детерминированный алгоритм
проверки простоты числа, который выдает доказательство того,
что число простое, за полиномиальное время O(n^6), где n --
число цифр в записи проверяемого числа.
2. Взлом ключа: по открытому ключу (m, e) надо вычислить
секретный ключ (m, d), т.е. вычислить число d:
d*e == 1 (mod phi(m))
Для этого надо знать функцию Эйлера phi(m)
phi(m) = (p - 1)(q - 1)
а для этого надо получить разложение m на множители
(эти две задачи эквивалентны по сложности).
ПЕРЕРЫВ 10 минут до 10:55
Пример схемы RSA
Сгенерируем
p = 1151438571896145047887447723231
q = 668203938391714894132076973151
m = p*q = 769395788557135886307479507219297725164204544725096865970881
Вычислим функцию Эйлера
phi(m) = (p-1)(q-1) =
= 769395788557135886307479507217478082653916684783077341274500
Возмем случайное e, обратимое в Z_phi(m)
e = 946665994505513
Вычислим обратный к e элемент в Z_phi(m):
d = 189336738334574291258425100105243158327302667770578542246577
Зашифруем конкретное сообщение x:
x = 123456789
y = E(x) == x^e (mod m) =
= 756333665574068956010618118370905227603762870770020946841931
Расшифруем зашифрованное сообщение y:
z == y^d (mod m) == 123456789
--------------------------------------------------------------------
Вероятностный тест простоты
1. В качестве теста простоты естественно использовать Малую теорему Ферма:
m простое ==> для всякого b != 0 (mod m)
b^(m-1) == 1 (mod m)
Это необходимое условие простоты. Будет ли оно достаточным?
Древние греки считали, что верно обращение Малой теоремы Ферма
для основания 2:
2^(m-1) == 1 (mod m) => m простое
Это неверно! П. Сарруз, 1830 привел контрпример:
m = 341 = 31*11
2^(m-1) == 1 (mod 341)
2^340 = (2^10)^34 = (1024)^34 = (341*3 + 1)^34 == 1^34 == 1 (mod 341)
То, что 341 не простое число, можно увидеть, взяв другое основание,
например 3:
3^340 == 56 (mod 341)
Но оказалось, что существуют такие "патологические" числа, которые
ведет себя как простые в Малой теореме Ферма, но не являются простыми:
для всякого основания b, взаимно простого с m (gcd(b, m) = 1),
выполняется Малая теорема Ферма:
b^(m-1) == 1 (mod m)
Такие числа теперь называются Кармайкловыми в честь американского математика
Роберта Кармайкл (Carmichael), открывшего в 1910 г. знаменитое число
561. Недавно доказано, что множество Кармайкловых чисел бесконечно:
561, 1105, 1729, ...
Для подобных чисел Малая теорема Ферма не годится в качестве теста
простоты. Но, оказывается, несложная модификация теста Ферма
дает нам вероятностный тест простоты.
Вероятностный тест использует датчик случайных чисел
и для входного числа m может выдавать один из двух ответов:
1) число m составное;
2) не знаю.
При этом в случае ответа 1 тест выдает доказателдьство того,
что число составное.
И для составного числа m вероятность ответа 2 не превосходит 1/4.
Например, если при 20 запусках тест 20 раз выдал ответ "не знаю",
то вероятность того, что число m составное, не превышает
(1/4)^20 = (1/2)^40 = (1/2^10)^4 <= (1/1000)^4 = 10^(-12)
То есть мы уверены, что число m простое, но мы не получаем доказательства
эттого факта. Тем не менее этого достаточно для использования таких
чисел в криптографии.
Тест Миллера-Рабина (начало 1970х)
Дано нечетное число m, которое мы проверяем на простоту
Генерируем случайное основание 0 < b < m
Если gcd(b, m) != 1, то m составное.
В противном случае мы как бы проверяем Малую теорему Ферма.
В ней мы возводим b в степень m-1 в кольце Z_m.
Представим число m-1 в виде
m - 1 = 2^s * t, t -- нечетное число
Сначала возводим b в нечетную степень t и затем полученное число
возводим в квадрат s раз, в конце получим b^(m-1).
Получаем последовательность
x0 == b^t (mod m)
x1 = x0^2 (mod m)
x2 = x1^2 (mod m)
. . .
x_s = x_(s-1)^2 = b^(m-1) (mod m)
Рассматриваем последовательность
x0, x1, x2, ..., x_s
В ней каждый следующий элемент есть квадрат предыдущего.
Тогда тест выдает ответ 1 (m -- составное число) в следующих случаях:
1) когда последовательность имеет вид
*, *, ..., *, 1, 1, ..., 1
2) *, *, ..., * != 1 (не выполняется Малая теорема Ферма).
Здесь звездочкой обозначены числа, не равные +-1 (mod m).
В случае 1 имеем число x_i != +- 1 (mod m), x_i^2 == 1 (mod m).
Если бы m было простым числом, то Z_m было бы полем,
а в поле могут быть только 2 квадратных корня из 1:
1 и -1
Эти корни удовлетворяют уравнению x^2 - 1 = 0,
а по теореме Безу у многочлена степени n может быть не больше n корней.
a^2 == b^2 (mod m), m делится на 2 простых числа ==>
вероятность того, что это соотношение
нетривиально, >= 1/2.
(Тривиально, когда a == +-b (mod m).)
Пример: убедимся, что число 1729 (третье кармайклово число)
не является простым.
m = 1729
m - 1 = 1728 = 27 * 2^6
Возмем случайное b = 123
и вычислим последовательность
x0 == b^27 == 1464 (mod m)
x1 = x0*x0 == 1065 (mod m)
x2 = x1*x1 == 1 (mod m)
Последнее сравнение нетривиально => m составное.
--------------------------------------------------
Алгоритмы факторизации чисел
Первый алгоритм: колесо со спицами.
Как обычно пишут алгоритм разложения на множители?
делим m сначала на 2, а потом на все нечетные пробные делители,
начиная с 3:
2,
3, 5, 7, 9, 11, 13, 15, 17, ...
Как этот алгоритм обобщить?
Сначала убедимся, что число m не делится на 2, 3, 5,
и затем будем делить на пробные делители вида
30*k + s, где s <- {1, 7, 11, 13, 17, 19, 23, 29}, k = 1, 2, 3, ...
2, 3, 5,
7, 11, 13, 17, 19, 23, 29,
31, 37, 41, 43, 47, 49, 53, 59,
61, 67, 71, 73, 77, 79, 83, 89,
91, 97, 101, 103, 107, 109, 113, 119,
. . .
----------------------------------------------------------
12 Oct 2022
Прошлый раз была рассмотрена схема кодирования с открытым ключом RSA
m = p * q,
p, q -- большие простые числа (удовлетворяющие некоторым
дополнительным условиям).
phi(m) = (p - 1)(q - 1)
e -- любой обратимый элемент в Z_phi(m)
m, e -- открыто
разложение m на множители секретно => не можем вычислить
функцию Эйлера phi(m)
Открытый ключ: (m, e)
Шифрование
E: x --> x^e (mod m)
Для дешифровки вычисляется обратный к e элемент в Z_phi(m)
d: d*e == 1 (mod phi(m))
Секретный ключ: (m, d)
Дешифровка
D: y --> y^d (mod m)
D*E = E*D = 1
В связи с этой схемой возникают следующие задачи:
1) генерация больших простых чисел,
проверка простоты произвольного числа.
Мы рассмотрели вероятностный тест простоты Миллира-Рабина;
2) для взлома схемы RSA надо уметь разлагать большое число на множители.
Эта задача интересна и с чисто математической точки зрения.
Прошлый раз был рассмотрен алгоритм "Колесо со спицами"
Алгоритмы Полларда (методы "Монте-Карло")
Это два очень простых и красивых алгоритма, которые хорошо
работают на практике. Название "Монте-Карло" происходит от того,
что они срабатывают не всегда за разумное время, требуется некоторое
везение.
rho-алгоритм Полларда основан на поиске цикла
в рекуррентной последовательности
Дано целое число m > 0, известно, что оно составное.
Надо найти нетривиальный делитель числа m: p | m.
Мы работаем в кольце Z_m.
Рассмотрим полиномиальное отображение
f: Z_m --> Z_m
где f -- произвольный многочлен степени >= 2.
Как правило, берут многочлен
f(x) == x^2 + 1 (mod m)
Берем случайный начальный элемент
x0 <- Z_m
и вычисляем бесконечную последовательность
x0,
x1 = f(x0),
x2 = f(x1),
x3 = f(x2),
...
x_i+1 = f(x_i)
...
Получаем последовательность элементов Z_m:
x0, x1, x2, x3, ...
Поскольку Z_m -- конецное множество, то она циклическая
(в начале идет отрезок апериодичности, затем она входит в цикл).
Пусть p -- нетривиальный делитель числа m. Тогда ту же последовательность
можно рассмотреть по модулю p:
__ __ __ __
x0, x1, x2, x3, ...
___
где x_i == x_i (mod p) (это образы x_i при гомоморфизме Z_m --> Z_p).
Поскольку p < m, то с большой вероятностью длина цикла
во второй последовательности меньше, чем длина цикла в первой
последовательности, т.е. для i < j
x_i != x_j (mod m)
x_i == x_j (mod p)
Это означает, что p | (x_i - x_j) и m не делит (x_i - x_j)
Поэтому d = gcd(x_i - x_j, m) нетривиален:
p | d, m не делит d.
Как найти цикл в последовательности?
Самая простая схема:
x0 <--> x1
x1 <--> x3
x2 <--> x5
x3 <--> x7
. . .
На каждом шаге расстояние между сравниваемыми элементами возрастает на
единицу. Рано или поздно оба элемента войдут в цикл и, поскольку второй
элемент движется по орбите в два раза быстрее, он рано или поздно
догонит первый элемент. Вместо сравнения элементов x, y мы вычисляем
d = gcd(x - y, m)
Алгоритм заканчивается, когда d != 1. В удачном случае d != m,
неудача когда d = m.
Примерная оценка времени работы rho-алгоритма Полларда.
Для случайного полиномиального отображения f: Z_p --> Z_p,
где deg(f) > 1, длина орбиты случайного элемента порядка
O(sqrt(p))
То есть в худшем случае m = p*q, где простые p, q примерно
одного порядка, время работы алгоритма порядка
корень четвертой степени из m.
Если m -- 20-значное десятичное число, то количество шагов алгоритма
порядка 100000.
-----------------------------------------------------------
Второй интересный и простой алгоритм факторизации --
это p-1-алгоритм Полларда. Он основан на малой теореме Ферма.
Пусть p | m -- простой делитель числа m.
Возмем случайное число b0 <- Z_m, gcd(m, b0) = 1
(если gcd(b0, m) != 1, то нам повезло, мы решили задачу факторизации).
Возьмем некоторое число N -- верхняя граница (чем больше N,
тем сильнее наш алгоритм, но при этом он работает дольше; например,
можно взять N = 10,000,000). Вычислим следующую последовательность:
b1 == b0^2 (mod m) d = gcd(b1 - 1, m)
b2 == b1^3 (mod m) d = gcd(b1 - 1, m)
b3 == b2^4 (mod m) d = gcd(b1 - 1, m)
b4 == b3^5 (mod m) d = gcd(b1 - 1, m)
. . .
b_N-1 == (b_N-2)^N d = gcd(b1 - 1, m)
Окончательно мы вычислим b0^(N!) (mod m)
Пусть p | m -- простой делитель числа m и p-1 -- гладкое число (smooth).
Гладкость определяется относительно константы n и означает, что
p - 1 раскладывается в произведение небольших взимно простых множителей:
p - 1 = q1^e1 * q2^e2 * ... * q_k^e_k,
q_i^e_i <= N
Не гладкое число -- например, когда у него имеется большой простой делитель.
Пример гладкого числа:
p = 931718707015139 -- простое число,
p - 1 = 2 * 7 * 59 * 733 * 4483 * 343267 -- гладкое число
для N = 1000000
p - 1 -- гладкое для константы N, тогда
(p - 1) | N!, N! = (p - 1)*s
b0^(N!) == (b0^(p - 1))^s == 1^s == 1 (mod p)
по Малой теореме Ферма ==>
p | (b0^(N!) - 1) = (b_N-1 - 1) ==>
p | gcd(b_N-1 - 1)
ПЕРЕРЫВ 10 минут до 10:50
Реализуем p-1-алгоритм Полларда (программа совсем простая).
. . .
Пример применения алгоритма:
p = 6328690139
q = 6458144389
m = p*q = 40871594710902480071
>>> p1Factor(m)
number of steps = 1000000
1
При N = 1000000 число m не удалось разложить, поскольку
p - 1 = 2 * 3164345069 число очень негладкое;
q - 1 = 2^2 * 3 * 463 * 1162373 число гладкое для константы
N = 1162373
Увеличив N, мы раскладываем наше число:
>>> p1Factor(m, upperBound = 1200000)
number of steps = 1162373
6458144389
В чем проблема p-1-алгоритма? Если для всех простых делителей
p числа m число p-1 не гладкое, то p-1-алгоритм Полларда бессилен.
Lenstra предложил вариант p-1-алгоритма Полларда, основанный
на работе в группе точек эллиптической кривой.
Аналогия с алгоритмом Полларда: он работает в группе
обратимых элементов кольца Z_p, p -- простое число.
Эта группа имеет фиксированный порядок p - 1, и, если оно
не гладкое, то мы ничего не можем сделать.
В алгоритме Ленстра мы используем случайную эллиптическую кривую
над полем Z_p, точки эллиптической кривой образуют группу
(по сложению), ее порядок n варьируется в пределах
p+1 - 2*sqrt(p) <= n <= p+1 + 2*sqrt(p)
Если для случайной кривой n получилось не гладким числом,
то мы можем выбрать другую случайную кривую, порядок которой уже
может быть гладким числом.
Эллиптическая кривая задается уравнением
y^2 = x^3 + a*x + b
(кривая рассматривается на проективной плоскости)
На точках кривой вводится операция сложения:
1) 0 -- это бесконечно удаленная точка (x = 0, y = 1, z = 0)
в однородных координатах;
2) любая прямая пересекает элл. кривую в трех точках
p, q, s:
p + q + s = 0
s = -(p + q)
3) чтобы сложить две точки p, q, надо провести через них прямую
и третью точку пересечения прямой (p, q) с элл. кривой надо
отразить симметрично относительно оси абсцисс.
Относительно этой операции точки кривой образуют абелеву группу
(в аддитивной нотации).
В ней умножению соответствует сложение точек, а возведению
в целую степень n > 0 соответствует умножение точки на n:
Мультипликативная нотация:
b^n = b*b*b*...*b
n раз
Аддитивная нотация:
n*t = t+t+t+...+t
n раз
Для вычисления n*t используется алгоритм, аналогичный
алгоритму бустрого возведения в степень.
Очень важно!!!
Для вычисления суммы точек на эллиптической кривой мы используем только
4 арифметические операции с координатами: +, -, *, /
Координаты суммы точек есть рациональные выражения от
координат слагаемых.
Следовательно, эллиптические кривые можно рассматривать над любыми
полями, в частности, над конечным полем Z_p, где p -- простое число.
И даже можно рассматривать эллиптические кривые над коммутиативными
корьцами, в частности, над кольцом Z_m, только сложение уже не всегда
можно выполнить, а только тогда, когда мы можем вычислить обратный
элемент к знаменателю дроби. Над Z_m обратимы элементы,
взаимно простые с m. Если мы встретили необратимый элемент,
gcd(s, m) != 1,
то мы решили нашу задачу разложения числа на множители.
Алгоритм такой же точно, как p-1-алгоритм Полларда.
Пусть p | m -- простой делитель m, n = порядок эллиптической кривой,
которую мы используем; пусть нам повезло и n -- гладкое число
относительно константы N.
Тогда берем начальную (случайную) точку t на кривой и последовательно
вычисляем
t0 = t
t1 = 2*t0 d = gcd(t1.x, m)
t2 = 3*t1 d = gcd(t2.x, m)
t3 = 4*t2 d = gcd(t3.x, m)
. . .
t_N-1 = N*t_N-2 d = gcd(t_N-1.x, m)
t_N-1 = (N!)*t0
Поскольку n гладкое относительно N, то n | (N!) и по теореме
Лагранжа t_N-1.x == 0 (mod p) ==> p | t_N-1.x ==> d = gcd(t_N-1.x, m) --
нетривиальный делитель
Как выбрат случайную кривую и начальную точку на ней?
Если выбрать случайные a, b:
y^2 = x^3 + a*x + b
то непонятно, как найти точку на кривой. Поэтому
сначала выбираем совершенно случайную точку (x, y)
и случайное число a, а затем вычисляем b из уравнения
y^2 = x^3 + a*x + b
b == y^2 - x^3 - a*x (mod m)
----------------------------------------
Алгоритм вычисления квадратного корня в поле Z_p
p -- большое простое число (например, 100-значное десятичное).
p -- нечетное число ==> число 2 необратимо в кольце Z_phi(p) = Z_(p-1)
Дано число a, которое является квадратом в Z_p.
Надо вычислить x: x*x == a (mod p)
a == x^e (mod p) e обратимо в Z_p-1
x == a^d (mod p) где d*e == 1 (mod (p-1))
Но e = 2 необратимо в Z_p-1 ==> задача нетривиальная.
Предложение.
Число a != 0 (mod p) есть квадрат по модулю p <==>
a^((p-1)/2) == 1 (mod p)
Доказательство.
==>
Пусть a -- квадрат, a == x^2 (mod p)
По Малой теореме Ферма a^((p-1)/2) == x^(p - 1) == 1 (mod p)
Обратно. Воспользуемся предложением:
Для любого конечного поля группа его ненулевых элементов циклическая.
Доказать его можно, применив, например, теорему о разложении
любой конечной абелевой группы в прямую сумму примарных
циклических групп:
Z_(p^e)
A = Z_p1^e1 + Z_p2^e2 + ... + Z_pk^ek
Обозначим p_i^e_i = m_i:
A = Z_m1 + Z_m2 + ... + Z_mk
Такая группа явлется цикической тогда и только тогда,
когда все m_i взаимно простые.
Например, группа
Z_8 + Z_27 циклическая,
а группа Z_2 + Z_4 + Z_27 не циклическая.
В первом случае экспонента группы равна 8*27 = порядку группы
во втором случае экспонента равна 4*27 и она меньше, чем порядок группы.
Итак, для конечных абелевых групп экпонента группы равна ее порядку
тогда и только тогда, когда она циклическая.
Рассмотрим группу ненулевых элементов конечного поля.
От противного. Если она не циклическая, то
ее экспонента < порядка группы,
т.е. для любого ненулевого x
x^e = 1 в поле.
Мы получили, что любой ненулевой элемент является корнем уравнения
x^e - 1 = 0
где e строго меньше, чем число ненулевых элементов поля.
Но такое невозможно по теореме Безу: число корней многочлена не
превосходит его степени.
Итак, пусть a^((p-1)/2) == 1 (mod p)
Группа ненулевых элементов поля Z_p циклическая и порождается
элементом x. Тогда
a = x^k,
Если k нечетное, k = 2*s + 1 то
a^((p-1)/2) == x^(2*s + 1)*((p-1)/2) ==
== x^(2*s*(p-1)/2) * x^((p-1)/2) ==
== 1 * x^((p-1)/2) == x^((p-1)/2) != 1,
поскольку порядок элемента x равен в точности p-1.
Проверим это утверждение на примерах:
>>> p = 552716141
>>> x = 1234567
>>> a = x*x%p
>>> a
317276752
>>> pow(a, (p-1)//2, p)
1
Критерий выполняется для элемента a, который является квадратом.
>>> pow(x, (p-1)//2, p)
552716140
Видим, что x^2 == p-1 == -1 (mod p)
то есть x не является квадратом.
Отметим, что для всякого элемента h <- Z_p,
не являющегося квадратом,
h^((p-1)/2) == -1 (mod p)
потому что h^((p-1)/2) -- это квадратный корень из 1
(по Малой теореме Ферма), а в поле может быть только
2 квадратных корня из единицы: 1 и -1.
-------------------------------------------------
21 Oct 2022
Были рассмотрены
1) вероятностный алг. проверки простоты Миллера--Рабина;
2) алгоритмы факторизации Полларда:
-- rho-алгоритм;
-- p-1 алгоритм;
вариант -- алгоритм Ленстра факторизации, основанный на
эллиптических кривых.
-- алгоритм Диксона факторизации целого числа "квадратичное
решето".
Последний алгоритм не будем рассматривать подробно. Только идея:
Эйлер: если нашли a^2 == b^2 (mod m), причем равенство нетривиальное,
a != +- b (mod m)
то мы можем разложить m на множители.
a^2 - b^2 == 0 (mod m)
(a - b)*(a + b) == 0 (mod m)
a - b != 0, a + b != 0 (mod m) ==>
gcd(a - b, m) нетривиален,
gcd(a + b, m) нетривиален.
Алгоритм "Квадратичное решето":
раскладываем на множители число m
Берем конечный набор небольших простых чисел
B = { p1, p2, p3, ..., p_k }
p_i такие, что m -- квадрат (mod p_i)
Цель -- найти два квадрата, которые совпадают (mod m)
r = [ sqrt(m) ]
a1 = r + 1 раскладываем число a1^2 - m по мультипликативной базе B
a2 = r + 2 a2^2 - m
a3 = r + 3 a3^3 - m
. . .
В-число -- это число, которое раскладывается на множители из множества B.
n = p1^e1 * p2^e2 * ... * p_k^e_k
Если среди чисел вида (r + i)^2 - m найдутся k+1 B-число, то задача решена!
(r + i)^2 - m = p1^e1 * p2^e2 * ... * p_k^e_k
(e1, e2, ..., e_k)
Получим матрицу с k+1 строкой и k столбцами.
Строки матрицы линейно зависимы над Z_2 =>
есть некоторое подмножество строк { i1, i2, ..., i_s },
сумма которых есть строка с четными координатами.
Следовательно, произведение этих чисел (r + i)^2 - m
раскладывается по мультипликативной базе в произведение простых чисел
в четных степенях => можно из этого произведения извлечь квадратный
корень. Получили число b:
b^2 == ( (r + i1)*(r + i2)* ... * (r + i_s) )^2 (mod m)
----------------------------------------------------------------------
Прошлый раз мы начали рассматривать алгоритм извлечения квадратного
корня в поле Z_p, p -- большое простое число.
Предложение:
p -- простое число, p != 2
Тогда a <- Z_p, a != 0 (mod p) является квадратом <==>
a^((p-1)/2) == 1 (mod p)
Алгоритм
Дано: число a, которое является квадратом (mod p), p -- простое число != 2
Надо: найти число x: x^2 == a (mod p)
Рассмотрим остаток от деления p на 4: возможны 2 случая {1, 3}
Простой случай: p == -1 (mod 4)
Тогда p+1 делится на 4.
Можем сразу предъявить ответ:
x = a^((p + 1)/4) (mod p)
x^2 = a^((p + 1)/2) = a^((p - 1)/2 + 1) = a^((p-1)/2) * a
Так как a -- квадрат, то a^((p-1)/2) == 1 (mod p),
x^2 == a (mod p)
Пример:
>>> p = 1211260669079337762496640110987
>>> x = 12345678901234567891234567
>>> a = x*x % p
>>> a
687649847581137696074050457379
>>> r = pow(a, (p+1)//4, p)
>>> r
12345678901234567891234567
>>> r*r % p
687649847581137696074050457379
Сложный случай: p == 1 (mod 4)
Идея алгоритма похожа на идею алгоритма проверки простоты Миллера--Рабина:
p - 1 делится на 4 ==>
(p - 1)/2 -- четное число. Извлечем из него максимальную степень двойки:
(p - 1)/2 = 2^s * q, q -- нечетное число
Так как a -- квадрат (mod p), то
a^((p - 1)/2) == 1 (mod p)
Возведем a сначала в нечетную степень q и затем s раз возведем то
что получилось в квадрат, в результате получим 1.
t == a^q (mod p)
t1 == t^2
t2 == t1^2
t3 == t2^2
. . .
t_s == t_s-1^2 == 1 (mod p)
t, t1, t2, t3, ..., t_s == 1 (mod p)
Каждый следующий элеимент в этой последовательности есть
квадрат предыдущего. Она заканчивается единицей. Поэтому
эта последовательность имеет один из двух видов:
1, 1, 1, ..., 1
либо
*, *, ..., *, -1, 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 (mod p)
Итак, в этом случае ответ r. t, t1, t2, t3, ..., t_s ==
Итак, сложный случай: t == a^q != 1 (mod p)
Пусть также
r == a^((q + 1)/2) (mod p)
r^2 == a^(q + 1) == a^q * a == t*a (mod p)
Инвариант цикла
+------------------------------------------------------
| r^2 == t*a (mod p)
| последовательность квадратов
| t, t^2, t^4, ..., t^(2^k) == 1 (mod p)
| заканчивается единицей и имеет длину не больше s+1
+------------------------------------------------------
Идея: подправим числа r, t так, чтобы инвариант сохранился,
но длина последовательности квадратов {t, t^2, t^4, ..., 1}
уменьшилась бы.
Для этого возьмем любой элемент z, который не является квадратом.
z^((p-1)/2) == -1 (mod p)
Это означает, что последовательность квадратов
h = z^q
h1 = h^2
h2 = h1^2
...
h_s = h_s-1^2 == -1 (mod p)
Последовательность для t заканчивается единицей и имеет длину не больше s+1:
*, *, *, *, -1, 1
Последовательность для h имеет длину s+1 и заканчивается -1
*, *, *, *, *, *, *, -1
Выравняем эти последовательности по элементу -1:
t
*, *, *, *, -1, 1
*, *, *, *, *, *, *, -1
h b, b^2
Исправим элементы t, r:
r' = r*b (mod p)
t' = t*b^2 (mod p)
Проверим сохранение инварианта r^2 == t*a:
r'^2 = (r*b)^2 = r^2 * b^2 = t*a * b^2 = t*b^2 * a = t' * a
Последовательность квадратов для t' будет иметь меньшую длину, т.к.
она есть произведение двух последовательностей квадратов для t и для
b^2:
t
*, *, *, *, -1, 1
*, *, *, *, -1
b^2
------------------
*, *, *, *, 1, 1
В конце-концов мы добъемся, что t == 1 (последовательность
квадратов имеет минимальную длину) и r^2 == t*a == a.
Рассмотрим пример:
>>> p = 666919534863317
>>> x = 1234567890123
>>> a = x*x%p
>>> a
607441214017907
>>> pow(a, (p-1)//2, p)
1
>>> pow(2, (p-1)//2, p)
666919534863316 # Нашли не квадрат (mod p) = 2
>>> m = (p - 1)//2
>>> q = m
>>> s = 0
>>> while q%2 == 0:
... q //= 2
... s += 1
...
>>> q
166729883715829
>>> s
1 # (p-1)/2 = 2*q, q нечетное
>>> z = 2 # не квадрат
>>> h = pow(z, q, p)
>>> h
645289532741198
>>> zSeq = squareSequence(h, p)
>>> zSeq
[645289532741198, 666919534863316, 1]
>>> p
666919534863317
>>> t = pow(a, q, p)
>>> t
666919534863316
>>> tSeq = squareSequence(t, p)
>>> tSeq
[666919534863316, 1]
>>> r = pow(a, (q+1)//2, p)
>>> r
663165187127920
>>> r*r % p == t*a % p
True
>>> r = r * 645289532741198 % p
>>> t = t * 666919534863316 % p
>>> t
1
>>> r*r % p == t*a % p
True
>>> r*r%p
607441214017907
>>> a
607441214017907
ПЕРЕРЫВ до 17:12
Классы в Питоне. Объектно-ориентированный стиль программирования
Есть две задачи:
-- научиться использовать чужие классы
-- уметь создавать свои классы и ими пользоваться.
Пример области, где классы помогают в решении задач:
геометрия на плоскости и в трехмерном пространстве.
Рассмотрим геометрию на плоскости. Используем два класса:
вектор 2-мерного линейного пространства R2Vector
точка 2-мерного аффинного пространства R2Point
Они определы в модуле R2Graph.py (модуль R3Graph для работы
в трехмерном пространстве).
Примеры задач, которые удобно решать с помощью этих классов
Что можно делать с векторами:
1) сложение, вычитание, умножение на число;
2) скалярное произведение векторов, норма (длина) вектора;
3) нормализация вектора: сделать его длину равной 1, не меняя направления
u = R2Vector(1, 1)
u.normalize() # R2Vector(0.707, 0.707)
w = R2Vector(1, 2).normalized()
s = u*v # Скалярное произведение (dot product)
4) получение нормали к вектору u: нормаль -- это вектор, который
получается из вектора u поворотом на 90 градусов против часовой стрелки
n = u.normal() # Нормаль к вектору u
# u = (x, y) n = (-y, x)
Пример: найти расстояние от точки до прямой
t -- R2Point
прямая задается точкой и направляющим вектором
(p, v)
def distanceToLine(t, p, v):
n = v.normal()
n.normalize()
return abs((t-p)*n)
Второй пример: найти точку пересечения двух прямых
Дано: две прямые (p1, v1)
(p2, v2) (точка, направляющий вектор)
Пусть известно, что они не параллельны
Надо: найти их точку пересечения q
q = p1 + v1*t поскольку q принадлежит первой прямой
n = v2.normal()
n*(q - p2) = 0 поскольку q принадлежит второй прямой
n*(p1 + v1*t - p2) = 0 найдем t
n*(p1 - p2) + t*(n*v1) = 0
t = n*(p2 - p1)/(n*v1)
def intersectLines(p1, v1, p2, v2):
n = v2.normal()
t = n*(p2 - p1)/(n*v1)
q = p1 + v1*t
return q
Задача: вычислить окружность, вписанную в треугольник
from R2Graph import *
def inCircle(a, b, c):
ab = b - a
ab.normalize()
ac = c - a
ac.normalize()
bisa = ab + ac
ba = a - b
ba.normalize()
bc = c - b
bc.normalize()
bisb = ba + bc
(_, center) = intersectLines(a, bisa, b, bisb)
radius = center.distanceToLine(a, b - a)
return (center, radius)
----------------------------------------------------------------------
28 Oct 2022
Классы в Питоне
При освоении объектно-ориентированного стиля программирования
мы сталкиваемся с двумя типами задач.
1. Научиться (привыкнуть) пользоваться чужими классами
и библиотеками; вообще применять классы при решении задач.
Классы (объектно-ориентированный стиль) не обязательны для
решения чисто вычислительных задач; но без классов невозможно
обойтись в оконном и графическом программировании.
Прошлый раз мы рассмотрели использование готовых классов
R2Vector и R2Point для решения геометрических задач на
плоскости.
2. Создание своих собственных классов.
Классы в Питоне совсем не похожи на классы в C++, Java, C# и т.п.
Отличие в том, что Питон -- язык с динамической типизацией,
а языки C++, Java, C# -- со статической.
Описание класса в C++ -- это как бы чертех объекта данного типа
(имя класса -- это тип объекта). В C++ все объекты данного класса
имеют одно и то же строение (данные-члены во всех объектах класса
одни и те же и описываются в h-файле).
В отличие от C++, у класса Питона нет никакого "чертежа".
Разные объекты одного и того же класса в принципе
могут иметь разное строение (разный набор данных-членов класса),
хотя программисты стараются этого избегать. Как правило, данные-члены
объекта создаются в его конструкторе, но могут создаваться и в любом
другом методе и даже при использовании объекта.
Существенное отличие классов Питона от С++ в том, что в Питоне нет
полиморфизма (точнее, он просто не нужен). Напомним, что в С++
могут быть несколько функций или несколько методов класса с одинаковыми
именами, но с разными наборами типов формальных параметров.
Совокупность имени метода и типов его формальных параметров (а также
константность или неконстантность) метода называется его сигнатурой.
Как правило, у любого класса в С++ есть несколько конструкторов с
разными типами аргументов. В Питоне конструктор всегда один,
но в нем могут выполняться разные действия в зависимости от типов
аргументов. В отличие от С++, в Питоне тип аргумента может быть
определен во время выполнения с помощью функции type. Например,
можно написать следующий конструктор класса R2Vector:
class R2Vector:
def __init__(self, x = 0., y = 0.):
if type(x) == list or type(x) == tuple:
self.x = float(x[0])
self.y = float(x[1])
else:
self.x = float(x)
self.y = float(y)
В результате вектор можно создавать, указывая либо список или кортеж,
содержащий 2 числа (аргумент один!), либо указывая две координаты
по отдельности:
u = R2Vector([1, 2])
v = R2vector(3, 4)
В отличие от С++, в Питоне ссылка на объект, к которому применяется метод,
указывается явно и всегда является первым аргументом метода. Обычно
ссылка на объект (на самого себя) называется self, хотя это правило
неформальное (но оно строжайше соблюдается всеми программистами!).
Специальные методы, реализующие либо стандартные операции языка для
объектов класса, либо наиболее важные стандартные функции, начинаются
и заканчиваются двумя символами подчеркивания. Например, сложение векторов
задается как метод __add__:
class R2Vector:
. . .
def __add__(self, v):
return R2Vector(self.x + v.x, self.y + v.y)
Когда в Питоне выполняется сложение двух векторов
w = u + v
то метод __add__ применяется к первому операнду сложения, а второй
операнд передается как параметр. Приведенная строка эквивалентна
следующей:
w = u.__add__(v)
(т.е. первая строка -- это синтаксический сахар для второй).
Замечательная особенность классов в Питоне состоит в том, что
для четырех арифметических операций есть методы, которые применяются
ко второму аргументу (в С++ ничего подобного нет!). Например,
нам хотелось бы реализовать умножение вектора на число как слева,
так и справа. Для этого надо реализовать два метода:
class R2Vector:
. . .
def __mul__(self, v):
'''return self*v'''
if type(v) == R2Vector:
# This is the dot product of two vectors
return self.x*v.x + self.y*v.y
else:
# Multiply the vector by the number
return R2Vector(self.x*v, self.y*v)
def __rmul__(self, c):
'''return c*self'''
return R2Vector(c*self.x, c*self.y)
Буква r в названии метода означает, что объект, к которому применяется
метод, является правым операндом умножения (right). Пример:
>>> v = R2Vector(5, 1)
>>> w = 2*v
>>> print(w)
(10, 2)
Здесь во второй строке к вектору v применяется метод __rmul__
-------------------------------------------------
11 Nov 2022
Модуль numpy для работы с массивами, матрицами, многомерными массивами
(тензорами)
Зачем нужен такой модуль? Ведь матрицы можно представлять средствами
самого Питона как список строк матрицы, каждая строка -- это список чисел.
1 2 3 4
5 6 7 8
9 10 11 12
a = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]
на самом деле, так никто не делает, такой способ представления матриц
неэффективен. Эффективная реализация матрицы предполагает, что элементы
матрицы хранятся в виде непрерывного (одномерного) массива, причем
все элементы матрицы имеют одинаковый тип.
Приведенная выше матрица при работе в языке C/C++ хранится в линейном
массиве:
1 2 3 4 5 6 7 8 9 10 11 12
0-я строка 1-я строка 2-я строка
Элемент a_i,j хранится в линейном массиве с индексом
i*n + j
для матрицы размера m x n (m строк, n столбцов).
import numpy as np
Основной объект в модуле numpy -- это многомерный массив
np.ndarray
n-dimensional array
Для освоения массивов и матриц попробуем реализовать несколько
алгоритмов линейной алгебры. Большая часть алгоритмов сводится
к методу Гаусса приведения матрицы к ступенчатому виду.
Ступенчатая матрица (row echelon form):
* ? ? ? ? ? ? ?
0 0 * ? ? ? ? ?
0 0 0 * ? ? ? ?
0 0 0 0 0 0 * ?
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
Здесь звездочками обозначены ненулевые числа,
вопросительными знаками -- любые числа.
Всякая матрица приводится к ступенчатому виду элементарными
преобразованиями двух типов:
1) обмен двух строк матрицы swapRows(a, i, j)
2) к одной строке матрицы прибавляем другую строку,
умноженную на число addRows(a, i, j, coeff)
Зачем нужна ступенчатая матрица?
Ранг ступенчатой матрицы равен рангу исходной матрицы и равен
числу ненулевых строк ступенчатой матрицы.
Определитель ступенчатой матрицы равен произведению диагональных элементов,
умноженному на (-1)^k, где k -- число произведенных обменов строк.
Ступенчатый вид позволяет решить линейную систему.
Задачи
1. Решить систему линейных уравнений с невырожденной квадратной матрицей.
2. Найти обратную матрицу.
3. Дан список векторов, векторы представляются одномерными numpy-массивами.
Найти среди них максимальную линейно-независимую систему.
maxIndependentSystem(vectorList)
Функция должна вернуть список индексов векторов из данного списка,
которые образуют базис линейной оболочки этих векторов.
Метод Гаусса:
чем отличается программистский метод Гаусса от чисто математического?
1. Мы сравниваем числа не с нулем, а с eps (по абс. величине):
считаем, что, если abs(a_ij) <= eps, то a_ij -- "нулевой" элемент.
2. В качестве разрешающего элемента надо выбирать максимальный
по модулю элемент столбца.
a[k] -= a[i]*(a[k, j]/a[i, j])
=> число (a[k, j]/a[i, j]) должно быть <= 1 по абс. величине,
иначе ошибки округления будут экспоненциально расти.
Как решаются задачи домашнего задания с помощью метода Гаусса?
1. Решение линейной системы уравнений (СЛАУ)
def solve(a, b)
'''Solve the system a@x = b for non-singular square matrix a
and an array of free terms b'''
Составим расширенную матрицу
1 2 3 2
a = 4 5 6 b = 3
7 8 0 5
расш. матрица
1 2 3 | 2 Гаусс * * * | * Обратный 1 0 0 | x0
с = 4 5 6 | 3 ==> 0 * * | * ==> 0 1 0 | x1
7 8 0 | 5 0 0 * | * проход 0 0 1 | x2
Для преобразования вектора-строки к матрице-столбцу
можно использовать функцию np.reshape:
>>> b
array([2., 3., 5.])
>>> b.reshape(3, 1)
array([[2.],
[3.],
[5.]])
Можно задать размер вдоль одной из осей как -1, тогда функция сама
вычислит правильный размер:
>>> b.reshape(-1, 1)
array([[2.],
[3.],
[5.]])
Для приписывания столбца свободных членов удобно использовать
функцию np.hstack:
>>> a
array([[1., 2., 3.],
[4., 5., 6.],
[7., 8., 0.]])
>>> c = np.hstack((a, b.reshape(-1, 1)))
>>> c
array([[1., 2., 3., 2.],
[4., 5., 6., 3.],
[7., 8., 0., 5.]])
Аргумент функции hstack -- список или кортеж матриц, которые соединяются горизонтально.
2. Вычислить обратную матрицу
Составим расширенную матрицу, приписав справа единичную матрицу:
1 2 3 1 0 0
a = 4 5 6 e = 0 1 0
7 8 0 0 0 1
расш. матрица
1 2 3 | 1 0 0 Гаусс * * * | * * * Обратный 1 0 0 | * * *
с = 4 5 6 | 0 1 0 ==> 0 * * | * * * ==> 0 1 0 | * * *
7 8 0 | 0 0 1 0 0 * | * * * проход 0 0 1 | * * *
Справа получаем обратную матрицу.
Получение единичной матрицы:
>>> e = np.zeros((3, 3))
>>> e
array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]])
>>> for i in range(e.shape[0]):
... e[i, i] = 1.
...
>>> e
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
Для приписывания справа единичной матрицы используем
функцию np.hstack:
>>> c = np.hstack((a, e))
>>> c
array([[1., 2., 3., 1., 0., 0.],
[4., 5., 6., 0., 1., 0.],
[7., 8., 0., 0., 0., 1.]])
========================================================
ВЕСЕННИЙ СЕМЕСТР 2022-23 уч. года
10 Feb 2023
Разработка оконных приложений на Питоне
с помощью модуля tkinter
Задача "построение графика интерполяционного многочлена
по заданным узлам интерполяции"
Задача линейной регрессии и ее решение методом
наименьших квадратов
Задача интерполяции
Даны узлы интерполяции
x0, x1, x2, ..., x_n x_i != x_j
y0, y1, y2, ..., y_n
Существует единственный многочлен степени n
f(x) = a0*x^n + a1*x^(n-1) + ... + a_n
такой, что
f(x_i) = y_i, i = 0, 1, ..., n
(имеем систему с матрицей Вандермонда, определитель которой
отличен от нуля)
x0^n x0^(n-1) ... 1
x1^n x1^(n-1) ... 1
. . .
xn^n xn^(n-1) ... 1
Многочлен в форме Ньютона
p_n(x) = a0 + a1(x-x0) + a2(x-x0)(x-x1) + ... + a_n(x-x0)...(x-x_n-1)
Добавим узел x_n+1, y_n+1
p_n+1(x) = p_n(x) + a_n+1(x-x0)(x-x1)...(x-x_n-1)(x-x_n)
a_n+1 = (p_n+1(x) - p_n(x))/( (x-x0)(x-x1)...(x-x_n-1)(x-x_n) )
подставим x = x+1
a_n+1 = (y_n+1 - p_n(x_n+1)/( (x_n+1-x0)(x_n+1-x1)...(x_n+1-x_n-1)(x_n+1-x_n) )
ПЕРЕРЫВ 10 минут до 18:30
Задача линейной регрессии
Есть две величины X, Y, X -- независима переменная, Y -- зависимая
Хотим построить модель зависимости между X и Y
Модель определяется набором параметров w0, w1, ..., w_n
Зависимость от параметров линейная
Y = w0*f0(X) + w1*f1(X) + ... + w_n*f_n(X)
В частности, рассмотрим задачу нахождения наилучшего интерполяционного
многочлена
Y = w0*x^n + w1*x^n-1 + ... + w_n = p(x)
пусть заранее известно, что степень многочлена deg(p) = n
и нам известны значения Y в m точках, где m >= n
Если m > n, то точного решения в общем случае нет,
попробуем найти решение, минимизирующее ошибку.
В качестве функции ошибки берется сумма квадратов разностей
(квадрат евклидова расстояния между [p(X_i)] и [Y_i]
для тренировочных данных [X_i, Y_i]
w0*X_0^n + w1*X_0^n-1 + ... + w_n = Y_0
w0*X_1^n + w1*X_1^n-1 + ... + w_n = Y_1
. . .
w0*X_m^n + w1*X_m^n-1 + ... + w_n = Y_m
Несовместная система для m > n
AW = Y
A: R^n --> R^m
AW -- это линейное подпространство в R^m
A -- матрица Вандермонда
X0^n X0^n-1 ... 1
X1^n X1^n-1 ... 1
. . .
X_m^n X_m^n-1 ... 1
Попробуем найти такое решение W = V, что квадрат расстояния
между AV и Y миимален.
Понятно, что V -- это основание перпендикуляра,
опущенного из Y на подпространство AW,
то есть вектор AV - Y перпендикулярен подпространству AW
(AV - Y, AW) = 0
(AV - Y)^T * AW = 0
Так как это выполнено для любого w <- R^n, то
(AV - Y)^T * A = 0
V'A'A - Y'A = 0 (через штрих обозначено транспонирование)
V'A'A = Y'A транспонируем обе части равенства
A' A V = A' Y
Это нормальное уравнение в методе наименьших квадратов
относительно неизвестного вектора V.
У матрицы A столбцы линейно независимы =>
матрица A'A обратима =>
V = (A' A)^(-1) A' Y
V = inverse(A' A) A' Y
Определение:
матрица
pinv(A) = inverse(A' A) A'
называется псевдообратной к матрице A
Решение задачи в методе наименьших квадратов:
V = pinv(A)*Y
10 Feb 2023
Разработка оконных приложений на Питоне
с помощью модуля tkinter
Задача "построение графика интерполяционного многочлена
по заданным узлам интерполяции"
Задача линейной регрессии и ее решение методом
наименьших квадратов
Задача интерполяции
Даны узлы интерполяции
x0, x1, x2, ..., x_n x_i != x_j
y0, y1, y2, ..., y_n
Существует единственный многочлен степени n
f(x) = a0*x^n + a1*x^(n-1) + ... + a_n
такой, что
f(x_i) = y_i, i = 0, 1, ..., n
(имеем систему с матрицей Вандермонда, определитель которой
отличен от нуля)
x0^n x0^(n-1) ... 1
x1^n x1^(n-1) ... 1
. . .
xn^n xn^(n-1) ... 1
Многочлен в форме Ньютона
p_n(x) = a0 + a1(x-x0) + a2(x-x0)(x-x1) + ... + a_n(x-x0)...(x-x_n-1)
Добавим узел x_n+1, y_n+1
p_n+1(x) = p_n(x) + a_n+1(x-x0)(x-x1)...(x-x_n-1)(x-x_n)
a_n+1 = (p_n+1(x) - p_n(x))/( (x-x0)(x-x1)...(x-x_n-1)(x-x_n) )
подставим x = x+1
a_n+1 = (y_n+1 - p_n(x_n+1)/( (x_n+1-x0)(x_n+1-x1)...(x_n+1-x_n-1)(x_n+1-x_n) )
ПЕРЕРЫВ 10 минут до 18:30
Задача линейной регрессии
Есть две величины X, Y, X -- независима переменная, Y -- зависимая
Хотим построить модель зависимости между X и Y
Модель определяется набором параметров w0, w1, ..., w_n
Зависимость от параметров линейная
Y = w0*f0(X) + w1*f1(X) + ... + w_n*f_n(X)
В частности, рассмотрим задачу нахождения наилучшего интерполяционного
многочлена
Y = w0*x^n + w1*x^n-1 + ... + w_n = p(x)
пусть заранее известно, что степень многочлена deg(p) = n
и нам известны значения Y в m точках, где m >= n
Если m > n, то точного решения в общем случае нет,
попробуем найти решение, минимизирующее ошибку.
В качестве функции ошибки берется сумма квадратов разностей
(квадрат евклидова расстояния между [p(X_i)] и [Y_i]
для тренировочных данных [X_i, Y_i]
w0*X_0^n + w1*X_0^n-1 + ... + w_n = Y_0
w0*X_1^n + w1*X_1^n-1 + ... + w_n = Y_1
. . .
w0*X_m^n + w1*X_m^n-1 + ... + w_n = Y_m
Несовместная система для m > n
AW = Y
A: R^n --> R^m
AW -- это линейное подпространство в R^m
A -- матрица Вандермонда
X0^n X0^n-1 ... 1
X1^n X1^n-1 ... 1
. . .
X_m^n X_m^n-1 ... 1
Попробуем найти такое решение W = V, что квадрат расстояния
между AV и Y минимален.
Понятно, что V -- это основание перпендикуляра,
опущенного из Y на подпространство AW,
то есть вектор AV - Y перпендикулярен подпространству AW
(AV - Y, AW) = 0
(AV - Y)^T * AW = 0
Так как это выполнено для любого w <- R^n, то
(AV - Y)^T * A = 0
V'A'A - Y'A = 0 (через штрих обозначено транспонирование)
V'A'A = Y'A транспонируем обе части равенства
A' A V = A' Y
Это нормальное уравнение в методе наименьших квадратов
относительно неизвестного вектора V.
У матрицы A столбцы линейно независимы =>
матрица A'A обратима =>
V = (A' A)^(-1) A' Y
V = inverse(A' A) A' Y
Определение:
матрица
pinv(A) = inverse(A' A) A'
называется псевдообратной к матрице A
Решение задачи в методе наименьших квадратов:
V = pinv(A)*Y