11 Sep 2021
Обзор курса
http://mech.math.msu.su/~vvb/ Radio-MSU.net
пройти по ссылке Algebra practicum
http://mech.math.msu.su/~vvb/AlgebraPract/
Главная цель курса -- знакомство с системой компьютерной математики
(компьютерной алгебры) SageMath.
Подобные системы: Volfram Mathematica, Maple, MatLab, Maxima, ...
Почему именно sage?
1. Это свободная и бесплатная система.
Свободные системы почти всегда лучше коммерческих.
SageMath -- объединяет все свободные математические пакеты
(GAP, Singular, Maxima, ...).
2. SageMath -- расширение языка Python3, не нужно учить еще один
язык программирования. Python, если даже его не знаешь, то обязательно
следует его выучить.
3. Python позволяет использовать компьтер как интеллектуальный
калькулятор.
Начнем с изучения Python3.
Для нас Python хорош еще и тем, что в нем используются неограниченные
целые числа, что позволяет писать программы для теории чисел, криптографии,
использовать рациональные числа Fraction.
Введение в Python
Tutorials -- вместо книг
Базовые типы: их три
int, float (= double языка C)
str, bool
Нет типа char, short, long, unsigned...
Конструкции:
список (на самом деле, динамический массив)
l = [2, 3, 5]
l[1] == 3
s = l[1:3]
tuple (кортеж)
t = (11, 13, 17, 19)
Распаковка тапла:
(x, y) = (10, 20)
Главная особенность Питона:
это язык с динамической типизацией
Переменные не имеют типов (операторы описания типов отсутствуют!)
Тип выражения невозможно определить статически,
он определяется только на Runtime.
=> не нужен полиморфизм функций и методов класса.
def length(x):
if type(x) == list:
return len(x)
elif type(x) == R2Vector:
return x.norm()
else:
raise ValueError("Non-supported type of argument")
Напишем программу, которая печатает все простые числа <= n
=======================================
25 Sep 2021
Мы познакомились с языком Python3.
Просьба в письмах с задачами писать в поле
Subject: 302 группа, задача 1
Типичные ошибки:
p - 1 = q * 2^s
q = p - 1
s = 0
while q%2 == 0:
q /= 2 #### Error!!!
s += 1
Надо писать //, если мы работаем с целыми числами (операция floordiv,
обычное деление truediv):
q //= 2
При возведении числа a в степень n по модулю m
категорически нельзя писать
b = (a**n)%m
Надо все действия выполнять в кольце Z_m:
b = pow(a, n, m)
Помните, что любой код вне функций надо выполнять условно! Только когда
ваша программа выполняется как скрипт, но не когда она используется как
модуль.
if __name__ == "__main__":
n = int(input("Enter n: "))
print(primeTest(n))
def main():
. . .
. . .
if __name__ == "__main__":
main()
Элементы теории чисел
=====================
Рассмотрим задачи, связанные со схемой RSA кодирования с открытым ключом.
s1 s2 s3 r s3^-1 s2^-1 s1^-1 r*r = 1
Первая идея криптографии: блочное кодирование -- любое сообщение
представляется в виде двоичной последовательности, эта последовательность
разбивается на блоки фиксированного размера, каждый блок шифруется
по-отдельности.
Стандарты блочного кодирования: DES, Blue Fish.
Вторая идея: блок можно трактовать как двоичную запись целого числа.
Более того, поскольку блоки имеют ограниченный размер, то их можно
рассматривать как элементы кольца вычетов по модулю m.
Таким образом, шифрование и дешифровку можно рассматривать
как взаимно обратные отображения кольца вычетов:
E: Zm --> Zm (Encryption) -- инъективно
D: Zm --> Zm (Decryption)
D E = 1 (тождественное преобразование Zm).
RSA -- придумали профессоры из MIT Rumley, Shamir, Adleman (?)
Простой алгоритм, не представляющий никакого гос. секрета,
общеизвестный, основанный на простых результатах теории чисел.
Алгоритм RSA использует асимметричное кодирование:
Процедура шифрования E открыта и общеизвестна (public key)
Процедура расшифровки D секретная (private key)
По процедуре E невозможно восстановить D.
Но схема RSA обладает еще одним замечательным свойством:
Отображения E и D биективные, то есть
не только DE = 1, но и ED = 1.
Это позволяет аутентифицировать отправителя сообщения:
исходное сообщение s
отправляется t = D(s)
Все его могут прочитать, но сформировать сообщение может только
обладатель секретной процедуры D. На этой идее основа электронная
подпись.
Третье: для взлома алгоритма RSA надо решить математическую проблему
факторизации целого числа => никакой хакер не в состоянии это сделать.
Схема RSA основа на элементарных результатах теории чисел
Малая теорема Ферма
p -- простое, b != 0 (mod p) =>
b^(p - 1) = 1 (mod p)
Форма Малой теоремы Ферма: для любого b
b^p = b (mod p)
Нам понадобится следующая форма этой теоремы: для любого b <- Zp
и любого целого k >= 0
b^(k*(p-1) + 1) = b (mod p)
Возведение в степень k*(p-1) + 1 есть тождественное отображение Zp в себя.
Китайская теорема об остатках.
Пусть m = m1*m2*...*mk, где числа m_i взаимно простые.
Тогда имеет место изоморфизм
Z_m == Z_m1 + Z_m2 + ... + Z_mk (сумма прямая)
причем изоморфизм вычисляется алгоритмически в обе стороны
(есть алгоритм для китайской теор. об остатках, вычисляющий
число по набору его остатков).
Пример применения:
сколько квадратных корней из единицы в кольце Z_105?
(105 = 3*5*7).
Ответ: 8 корней.
Решение:
Z_105 == Z_3 + Z_5 + Z_7
x == (x1, x2, x3) x1 = x (mod 3), x2 = x (mod 5), x3 = x (mod 7)
x^2 = 1 <==> x_i^2 = 1
Z_3, Z_5, Z_7 -- конечные поля, в поле тольк0 2 квадратных корня из 1:
1 и -1
Значит, любой корень из 1 в Z_105 имеет вид
x = (+-1, +-1, +-1)
(всего 8 корней).
x = (-1, 1, -1)
Частный случай теоремы Эйлера.
Пусть целое число m свободно от квадратов:
m = p1 * p2 * ... * pk
Обозначим phi(m) = (p1 - 1)*(p2 - 1)*...*(pk - 1)
Тогда для любого b <- Z_m и любого целого числа s >= 0 выолняется
b^(s*phi(m) + 1) = b (mod m)
Схема RSA
=========
Возьмем 2 больших простых числа p, q
Пусть m = p*q
m открыто, разложение m на множители секретно.
phi(m) = (p - 1)*(q - 1) секретное число.
Возьмем любое (случайное) число e, обратимое в кольце Z_phi(m)
С помощью расширенного алгоритма Евклида вычислим обратный к e
элемент d:
d*e == 1 (mod phi(m))
Отображение шифрования -- это просто возведение в степень e в
кольце Z_m
E: x --> x^e (mod m)
Обратное преобразование -- это возведение в степень d в кольце Z_m:
D: y --> y^d (mod m)
Теорема:
DE = ED = 1
Доказательство:
y = E(x) == x^e (mod m)
z = D(y) == D(x^e) == x^(d*e) (mod m)
d*e == 1 (mod phi(m))
d*e = s*phi(m) + 1 = s*(p-1)*(q-1) + 1
По частному случаю теоремы Эйлера
x^(s*(p-1)*(q-1) + 1) == x (mod m)
--------------------------------------------
Пара (m, e), которая задает отображение E: Z_m --> Z_m -- это открытый ключ.
Пара (m, d) задает отображение D: Z_m --> Z_m -- секретный ключ
(секретное число d).
Для взлома схемы надо вычислить число d => надо вычислить phi(m) =>
надо разложить m на множители.
==================================================================
9 Oct 2021
Прошлый раз было рассказано про схему кодирования с
открытым ключом RSA
m = p*q, p, q простые
m открыто, разложение m секретно
e -- любое обратимое число в Z_phi(m)
phi(m) -- функция Эйлера
d -- обратный к e элемент в Z_phi(m)
d*e == 1 (mod phi(m))
e открытое, d секретное
Шифрование E: Z_m -> Z_m
E: x -> x^e (mod m)
Открытый ключ -- это пара (m, e)
Дешифровка D: Z_m -> Z_m
D: y -> y^d (mod m)
DE = ED = 1
Секретный ключ: (m, d)
Для взлома ключа надо по паре (m, e) вычислить d:
d*e == 1 (mod phi(m))
то есть надо знать функцию Эйлера phi(m) = (p-1)*(q-1)
поскольку m = p*q.
Для вычисления функции Эйлера надо разложить m на множители
(обратное тоже верно, если phi(m) известна, то m легко
раскладывается на множители).
Задача взлома ключа по сложности эквивалентна задаче
факторизации числа m. (Факторизация == разложение на множители.)
В настоящее время не существует алгоритмов факторизации целых чисел
с полиномиальным временем работы (для полиномов такой алгоритм
существует: алгоритм Берликемпа для полиномов над Z_p, алгоритм
Ленстров для полиномов над Z). Есть алгоритм Шора факторизации
числа на квантовом компьютере, также работающий за полиномиальное
время.
В связи со схемой RSA возникают две задачи:
1) генерация больших простых чисел для создания ключа;
2) факторизация целых чисел для взлома ключа (а также для оценки
его надежности).
Эти задачи представляют и чисто математический интерес.
Задача 1 (генерация простых чисел).
Можно рассмотреть похожую задачу: проверить простоту числа.
В 70-е годы был найден вероятностный тест Миллера--Рабина.
Для составного числа он с вероятностью 1-eps выдает доказательство
того, что оно составное (eps может быть сделано сколь угодно
маленьким). Если тест выдал ответ "не знаю", то число, скорее всего,
простое, однако доказательства этого не выдается.
Уже довольно давно есть полиномиальный алгоритм генерации простых
чисел с предъявлением сертификата простоты.
Число p простое <=> когда группа обратимых элементов кольца Z_p
является циклической группой порядка p-1.
Основано на том предложении, что мультипликативная группа конечного
поля является циклической.
Порождающий элемент группы ненулевых элементов поля Z_p
называется первообразным корнем порядка p:
a <- Z_p:
a^(p-1) == 1 (mod p)
и для всякого простого q | (p-1)
a^((p-1)/q) != 1 (mod p)
Сертификат простоты -- это число a, а также разложение числа
p-1 на простые множители:
p-1 = q1^e1 * q2^e2 * ... * q_k^e_k
Плюс надо еще предъявить сертификаты простоты для всех чисел q_i
Такие деревья сертификатов можно строить за полиномиальное время.
Современные достижения: совсем недавно (2005 г. ???) был найден
полиномиальный алгоритм проверки простоты, который предъяляет
доказательство простоты числа. (Сложность O(n^6), где n -- дляна записи
числа).
Предложение: число p простое <=> все биномиальные коэффициенты С(p, k)
простые при 1 < k < p.
В качестве теста простоты можно было бы попробовать использовать
Малую теорему Ферма:
b^(p-1) = 1 (mod p) -- необходимое условие простоты.
Древние греки ошибочно считали, что если
2^(p-1) = 1 (mod p)
то p простое.
Это не так, что показал французский математик П.Саррус
в начале XIX века:
2^340 = 1 (mod 341)
но 341 = 31*11
2^340 = (2^10)^34 = 1024^34 = (341*3 + 1)^34 == 1^34 == 1 (mod 341)
Можно показать, что 341 не простое с помощью основания b = 3
3^340 = 56 (mod 341).
Удивительный результат открыл в 1910 г. американский математик
Carmichael (Кармайкл). Он нашел замечательное число m = 561,
которое ведет себя как простое в Малой теореме Ферма для всех
оснований b, взаимно простых с m:
gcd(b, m) = 1 =>
b^(m-1) = 1 (mod m)
Такие числа в его честь называются кармайкловыми.
Число m кармайклово <=> 1) m = p1*p2*...*p_k, p_i простые,
2) k >= 3,
3) (p_i - 1) | (m - 1).
Недавно было показано, что множество кармайкловых чисел бесконечно.
Для кармайкловых чисел Малая теорема Ферма не годится в качестве теста
простоты.
Кармайкловы числа:
561, 1105, 1729, ...
Например, 1729 = 7 * 13 * 19
6 | 1728, 12 | 1728, 18 | 1728
Но тест Ферма можно легко переделать в полноценный вероятностный
тест Миллера--Рабина. Проверяем простоту нечетное числа m (работаем
в кольце Z_m).pract2021.txt
Берем случайное число b <- Z_m.
m - 1 -- четное числоэ
m - 1 = q*2^s, s >= 1
Вычисляем последовательность элементов Z_m:
b^q, b^(2*q), b^(4*q), ..., b^(m-1)
каждый следующий элемент есть квадрат предыдущего,
длина последовательности равна s+1.
Утверждение:
если эта последовательность имеет вид
*, *, *, 1, 1, ...,
где * != +-1, или вид
*, *, ..., *, где последняя * отлична от 1 (mod m),
то число m составное (тест выдает ответ "число составное").
Если же последовательность имеет вид
1, 1, ..., 1
или
*, *, ..., *, -1, 1, ..., 1
то число, возможно, является простым
(тест выдает ответ "не знаю").
Теорема:
1) если тест выдал ответ "число составное", то m -- действительно
составное число;
2) вероятность ответа "не знаю" для составного числа m не превышает 1/4
(не зависит от m).
Доказательство пункта 1.
Если есть x: x != +-1 (mod m), x^2 = 1 (mod m),
то m не простое, поскольку иначе Z_m было бы полем,
а в поле нет нетривиальных квадратных корней из 1
(кроме 1 и -1), т.к. уравнение
x^2 - 1 = 0pract2021.txt
имеет не более двух решений. Если же последовательность заканчивается
не единицей, то не выполняется Малая теорема Ферма.
Доказательство пункта 2 тоже несложное (идея: любая собственная подгруппа
имеет мощность <= 1/2 можности группы).
Алгоритмы факторизации целых чисел
Есть два совсем-совсем простых алгоритма Полларда факторизации,
так называемые "методы Монте-Карло":
1) rho-алгоритм Полларда;
2) p-1 алгоритм Полларда.
==================================================
23 Oct 2021
p-1 алгоритм факторизации целых чисел Полларда
Основан на Малой теореме Ферма:
b != 0 (mod p) => b^(p-1) == 1 (mod p)
Дано составное число m, надо разложить его на множители
(достаточно найти хотя бы один нетрививальный делитель d | m).
m = p1^e1 * p2^e2 * ... * p_k^e_k
Предположим, что p1 - 1 является гладким числом:
p - 1 = q1^l1 * q2^l2 * ... * q_s^l_s
Гладкое, если для всякого i
q_i^l_i <= N
N -- параметр алгоритма, например, N = 10,000,000
Не гладкое число -- число вида n = p*t, где p простое, p не меньше
квадратного корня из n.
Рассмотрим число N!. Поскольку p-1 -- гладкое число относительно N,
то p-1 | N!
Выберем случайное число b (base -- основание степени)
Вычислим n = b^N! (mod m)
Тогда p | gcd(n - 1, m),
поскольку по Малой теореме Ферма b^N! = b^((p-1)*t) == 1 (mod p)
Если повезет и это число не делится на m
m не делит gcd(n - 1, m),
то d = gcd(n - 1, m) -- нетривиальный делитель числа m
В схеме кодирования с открытым ключом m = p*q
числа p - 1 и q - 1 не должны быть гладкими! Иначе m можно разложить на
множители.
=================================================
Ленстра (Ленстры?) предложил вариант p-1 алгоритма Полларда,
работающий с эллиптическими кривыми, он основан на теореме Лагранжа
для группы точек эллиптической кривой над полем Z_p.
Элл. кривая -- кубическая кривая
y^2 = x^3 + a*x + b
Точки элл. кривой образуют абелеву группу с операцией +.
Ее порядок зависит от констант a, b.
Важно! Сумма двух точек вычисляется с помощью четырех арифметических
операций. Ноль этой группы -- бесконечно удаленная точка, ее координата x
равна нулю.
Можно написать аналог p-1 алгоритма Полларда для точек эллиптической
кривой над кольцом Z_m. (Аналогом возведения в степень n в
мультипликативной записи является умножение на n в аддитивной записи.)
При этом выполняется деление в Z_m.
Если его нельзя выполнить
(элемент n обратим, когда gcd(n, m) = 1)
то мы нашли необратимый элемент в Z_m => m разложили на множители.
В конце алгоритма получим ноль группы == бесконечно удаленная точка,
координата x которой равна нулю в Z_p => координата x делится на p =>
p | gcd(coord.x, m)
---------------------------------------
Алгоритм вычисления квадратного корня в поле вычетов по модулю p
Дано: нечетное простое число p,
элемент a <- Z_p
a == x^2 (mod p)
Надо: найти x.
Предложение:
элемент a <- Z_p яляется квадратом <=>
a^((p-1)/2) == 1 (mod p)
Доказательство: следует из того, что группа обратимых элементов
конечного поля циклическая. Элементы виды b^2k и только они являются
квадратами, где b -- порождающий элемент группы.
Имеем p == +- 1 (mod 4) (или 1, 3 (mod 4))
Случай p == -1 (mod 4) тривиален.
В этом случае p+1 делится на 4, и
x = a^((p+1)/4) есть квадратный корень из a.
x^2 = a^((p+1)/2) = a^((p-1)/2 + 1) = a^((p-1)/2) * a == a
Сложный случай: p == 1 (mod 4)
(p - 1)/2 -- четное число.
(p - 1)/2 = q*2^s, q -- нечетное число
Возводим число a в степень (p-1)/2 следующим образом:
сначала вычисляем
t = a^q (mod p)
и затем t возводим в квадрат s раз.
Получаем последовательность квадратов:
либо
a^q = t, t^2, t^4, ..., -1, 1, ..., 1 = a^((p-1)/2)
*, *, *, *, -1, 1, ..., 1
либо t с самого начала равно единице:
a^q = 1, 1, 1, ..., 1
Второй случай также тривиален:
Имеем t = a^q == 1 (mod p)
возьмем
r = a^((q+1)/2)
Тогда
r^2 = a^(q+1) = a^q*a = t*a == a
Получили r -- квадратный корень из a.
Теперь рассмотрим общий случай:
a^q = t, t^2, t^4, ..., -1, 1, ..., 1 = a^((p-1)/2)
Обозначим
r = a^((q+1)/2)
Оно не будет корнем, но для него выполняется сравнение:
r^2 = t*a (mod p), где t = a^q (mod p)
Если t возводить в квадрат, то мы получим единицу не более чем через s шагов:
t, t^2, t^4, ..., -1, 1
*, *, ..., *, -1, 1
последовательность имеет длину не больше s+1.
Эти два условия будут инвариантом цикла.
Попробуем, сохранив инвариант, уменьшить длину последовательности квадратов.
Подправим числа r и t.
Возьмем произвольный элемент z, который не является
квадратом по модулю p.
Рассмотрим последовательность квадратов при возведении z в степень (p-1)/2
z^((p-1)/2) == -1 (mod p)
z^q, z^(2*q), z^(4*q), ..., *, -1
*, *, *, ..., *, -1
длина последовательности равна s+1, последний элемент равен -1
Звездочек в последовательности квадратов для z^q
больше, чем в последовательности для t:
t: *, *, *, ..., -1, 1, ..., 1, 1
z: *, *, *, ..., *, *, ..., *, -1
Мысленно выровняем эти последовательности так, чтобы элемент -1
оказался на одной позиции:
t: *, *, *, ..., -1, 1, ..., 1, 1
z: *, *, ...,*, *, *, ..., *, -1
^ ^
| |
b b^2
Обозначим элемент последовательности для z, соответствующий началу
последовательности для t, через b^2, предыдущий элемент через b.
Подправим с помощью элемента b элементы r, t:
r' = r*b (mod p)
t' = t*b^2 (mod p)
Инвариант сохраняется:
r'^2 = (r*b)^2 = r^2*b^2 = t*a*b^2 = t'*a
Последовательность квадратов для t' будет иметь меньшую длину,
поскольку она есть поэлементное произведение двух последовательностей:
*, *, *, ..., -1, 1
t
*, *, ..., *, -1
b^2 перемножим
*, *, +-1, 1
t'
Длина последовательности квадратов для t' уменьшится
по крайней мере на 1 => рано или позно получим t = 1.
Следующий раз рассмотрим классы в Питоне.
С++
class R2Point {
public:
double x;
double y;
public:
R2Point(double xx = 0., double yy = 0.):
x(xx),
y(yy)
{}
. . .
};
===================================================
6 Nov 2021
Классы в Питоне
1. Классы в Питоне не похожи классы в C++.
Класс в C++ -- это как бы чертеж объекта класса,
в частности, все объекты класса в C++ имеют одинаковое
устройство (набор данных-членов класса).
Язык C++ -- язык со статической типизацией: тип выражения
можно определить по тексту программы, не обязательно для
этого запускать программу.
Python -- язык с динамической типизацией.
Плюсы -- гибкость, удобство программирования.
Минусы -- меньшая надежность, трудность понимания
чужих программ.
def f(x, y):
. . .
Не знаем ни типов формальных параметров, ни типа возвращаемого
значения. Более того, функция или метод класса могут возвращать
значения разных типов в разных ситуациях.
Объект класса в Питоне не имеет фиксированного внутреннего устройства,
в отличие от C++.
В Питоне нет полиморфизма, как в C++:
в С++ может быть несколько функций или методов с одинаковым именем,
но с разными сигнатурами (сигнатура -- это число и набор типов
формальных параметров + константность метода).
В Питоне не может быть двух функций или методов с одинаковым именем.
Зато в Питоне есть оператор type, который возвращает тип выражения.
Это является великолепной заменой полиморфизму.
========================
Встроенные типы в Питоне
int, float, comlex
конструкции:
list -- список (вернее, динамический массив == std::vector)
tuple -- кортеж (1, 2, "abc")
[1, 2, "abc"] -- список.
кортеж -- это неизменяемый (immutable) объект
set -- множество
s = { 2, 3, 5, 7 }
dict -- словарь, отображение == std::map в языке C++,
хранит множество уникальных ключей (все ключи разные)
и значения этих ключей. Он как бы осуществляет отображение
множества ключей на их значения.
d = { 2: "first prime", 3: "second prime", "five": "cinque" }
d[2] == "first prime"
d["five"] == "cinque"
===========================
Рассмотрим реализацию класса на примере двух простых классов:
класс R2Vector реализует векторы на плоскости R2
R2Point -- точки на плоскости
Второй пример:
кольцо вычетов по модулю m
Здесь число m -- это общий (статический) объект для всех объектов класса.
class Zm:
m = 7
@staticmethod # Decoration
def setMod(modul):
if modul == 0:
raise ValueError("Zero modul")
Zm.m = int(modul)
if Zm.m < 0:
Zm.m = -Zm.m
. . .
=========================================================
20 Nov 2021
class GFp2 -- конечное поле из p^2 элементов
Zp[x]/(x^2 - r), где r <- Zp не является квадратом
многочлен f(x) = x^2 - r неприводим,
значит, фактор по нему поле
def __init__(self, a = 0, b = 0):
# создается элемент a*x + b (плохо! лучше наоборот,
# a + b*x, чтобы можно было, например,
# писать GFp2(1) для элемента 1, а не x).
GFp2(1) -- это многочлен x, что неудобно
удобнее, чтобы это была единица.
x = GFp2(0, 1)
y = 1/x
y = GFp2(1)/x
y = GFp2(0, 1)/x
Для этого надо было реализовать метод
def __rtruediv__(self, y):
. . .
x/y == x.__truediv__(y)
1/x == x.__rtruediv__(1)
------------------------------------------------------------------
Новая (главная) тема:
система компьютерной математики SageMath
Система компьютерной алгебры
В алгебре мы имеем дело с точными (символьными) вычислениями,
в отличие от анализа.
Volfram Mathematika, MatLab, Maple, ...
Используем слово "математики" вместо "алгебры", поскольку в SageMath
возможны не только точные (символьные) вычисления, но и приближенные
вычисления.
Почему SageMath, а не, допустим, MatLab?
1. SageMath -- некоммерческая система, совершенно свободная с
исходными кодами. Плюсы -- не надо платить денег
или воровать, но главное -- свободные системы развиваются
намного быстрее коммерческих.
Особенность SageMath состоит в том, что она включает в себя
все существующие свободные пакеты в области математики
(Singular -- многочлены, GAP -- группы, Maxima, ...)
2. SageMath использует Python в качестве языка программимрования,
для работы с SageMath не нужно учить еще один язык
программирования.
3. SageMath очень хорош именно для алгебраистов.
Как установить SageMath?
Если у вас Linux, то нет проблем (просто устанавливаете пакет SageMath).
MS Windows -- лучше всего установить WSL2 (Windows System for Linux), т.е.
установить Linux Ubuntu как сервис под MS Windows 10.
Но можно даже работать в облаке, ничего не устанавливая
на локальном компьютере -- через Jupyter Notebook, или даже
(самый простой способ!) через SageMathCell
Сайт sagemath.org
Ссылка на SageMathCell на сайте:
https://sagecell.sagemath.org/
Пример: найти центр окружности, описанной вокруг треугольника
def incenter(a, b, c):
var("cx cy")
center = vector([cx, cy])
va = vector(a)
vb = vector(b)
vc = vector(c)
dist_ca = (center - va).norm()
dist_cb = (center - vb).norm()
dist_cc = (center - vc).norm()
eq1 = (dist_ca == dist_cb)
eq2 = (dist_ca == dist_cc)
res = solve([eq1, eq2], [cx, cy])
return (res[0][0].rhs(), res[0][1].rhs())
ГДЕ БЫЛА ОШИБКА????
*******************************************************************
Ошибка была в строке
c = vector([cx, cy])
Случайно имя переменной c, обозначавшей центр вписанной окружности,
совпало с названием третьей вершины треугольника. После исправления на
"center" функция заработала.
*******************************************************************
Геометрические задачи
Пример: нарисовать треугольник, вписанную окружность и процесс ее
построения (пересечение 3-х биссектрис)
Удобно использувать
1) функцию, которая находит пересечение двух прямых;
2) функцию, которая находит нормаль к вектору.
def normal(v): # normal(x, y) = (-y, x)
return vector([-v[1], v[0])
q -- пересечение прямых (p1, v1) и (p2, v2)
n = normal(v2) # Вектор нормали ко второй прямой
q = p1 + v1*t # q пробегает все точки первой прямой при t <- R
(q - p2)*n = 0 # условие того, что q принадлежит второй прямой
(p1 + v1*t - p2)*n = 0 # подставим выражение для q в уравнение 2-й прямой
(p1 - p2)*n + v1*n*t = 0 # и найдем число t
t*(v1*n) = (p2 - p1)*n
t = (p2 - p1)*n / (v1*n)
def intersectLines(p1, v1, p2, v2):
n = normal(v2)
t = (p2 - p1)*n / (v1*n)
q = p1 + v1*t
return q
В результате получим функцию, на вход которой подаются три
вершины треугольника и которая возвращает графический объект:
изображение треугольника (синий цвет),
окружности, вписанной в него (красный цвет) и
центра окружности (точка красного цвета, которая изображается
небольшим кружочком).
def normal(v):
return vector([-v[1], v[0]])
def intersectLines(p1, v1, p2, v2):
n = normal(v2)
t = ((p2 - p1)*n) / (v1*n)
q = p1 + v1*t
return q
def inCircle(a, b, c):
va = vector(a)
vb = vector(b)
vc = vector(c)
ab = (vb - va).normalized()
ac = (vc - va).normalized()
bisa = ab + ac
ba = (va - vb).normalized()
bc = (vc - vb).normalized()
bisb = ba + bc
center = intersectLines(va, bisa, vb, bisb)
n = normal(ab).normalized()
r = (center - va)*n
t = line([va, vb, vc, va], thickness=2)
circ = circle(center, r, color="red")
centre = point(center, size=40, color="red")
return t + circ + centre
-----------------------------------------------------
4 Dec 2021
Sage: построение кубического сплайна
(как пример задачи, где помогает использование SageMath)
Замечание по домашним заданиям
Реализация класса (например, многочлены над Zp)
1) надо обязательно реализовать методы
def __repr__(self):
. . .
def __str__(self):
. . .
Методы возвращают текстовое представление объекта.
Методы, которые что-то печатают, не нужны!
Слово print вообще не должно встречаться в тексте класса;
2) для полиномов надо реализовать метод деления
def __floordiv__(self):
. . .
а не __truediv__(self):
. . .
можно реализовать также метод, которое получает сразу целую часть
частного и остаток от деления: __divmod__
(выхывается как функция divmod(p, q) и возвращает пару (частное, остаток).
Вернемся к SageMath.
Задача на построение кубического сплайна
Пример: дано изображение (матрица пикселей), требуется
его увеличить (допустим, в 10 раз).
Задача сводится к растяжению матрицы, содержащей вещественные
числа от 0 до 1.
Можно сначала растянуть матрицу по горизонтали (в 10 раз),
а затем получившуюся прямоугольную матрицу растянуть по вертикали =>
задача сводится к растяжению только в одном направлении =>
каждую строку матрицы можно рассматривать независимо как функцию
от одного переменного, заданную в дискретном множестве узлов.
Построение кубического C1-сплайна основано на следующей задаче.
Даны значение функции на концах отрезка и значения ее производных:
f(x0), f(x1), f'(x0), f'(1)
Надо приблизить эту функцию кубическим многочленом, т.е. надо найти
многочлен p(x), удовлетворяющий четырем условиям:
p(x0) = y0
p(x1) = y1
p'(x0) = d0
p'(x1) = d1
Найдем выражения для коэффициентов этого многочлена. Применим SageMath.
var("x0 x1 y0 y1 d0 d1")
var("a0 a1 a2 a3")
p(x) = a0 + a1*x + a2*x^2 + a3*x^3
dp(x) = derivative(p(x), x)
# Условия на значения полинома и его производной на концах отрезка
eq1 = (p(x=x0) == y0)
eq2 = (p(x=x1) == y1)
eq3 = (dp(x=x0) == d0)
eq4 = (dp(x=x1) == d1)
res = solve([eq1, eq2, eq3, eq4], [a0, a1, a2, a3])
coeff0 = res[0]
Задача 2
Построим кубический сплайн с непрерывными первыми и вторыми производными
в промежуточных узлах.
Надо задать краевые условия: значения вторых производных в первом
и последнем узлах, тогда сплайн определен однозначно.
Узлы (x0, y0), (x1, y1), ..., (x_n, y_n) всего n+1 узел
Сплайн состоит из n сегментов,
всего 4*n коэффициентов (кубических многочленов)
p0(x), p1(x), ..., p_n-1(n)
p_i(x) определен на отрезке (x_i, x_i+1)
В крайнем (первом) узле -- 2 уравнения
p0(x) = y0
p''0(x) = 0
В промежуточных узлах по 4 уравнения в каждом.
В узле x_i имеем 4 уравнения:
p_i-1(x_i) = y_i значение предыдущего полинома
p_i(x_i) = y_i значение i-го полинома
p'_i-1(x_i) = p'_i(x_i) значения производных слева и справа совпадают
p''_i-1(x_i) = p''_i(x_i) значения вторых производных совпадают
В крайнем (последнем) узле -- 2 уравнения
p_n(x) = y_n
p''_n(x) = 0
Всего получаем 4*n уравнений на коэффициенты полиномов, которых тоже 4*n.
Получаем линейную систему с невырожденной матрицей.
Более того, матрица получается ленточной (9-диагональной) =>
методом Гаусса система решается за линейное время O(n)
a00 a01 a02 a03
a10 a11 a12 a13
a20 a21 a22 a23
-----------------------------------------------------
18 Dec 2021
Коммутативная алгебра
Рассматривае коммутативные ассоциативные кольца, порожденные конечным числом
элементов, или коммутативные алгебры над некоторым полем.
Алгебра = кольцо + векторное пространство над полем.
Теория коммутативных ассоциативных колец и алгебр связана с
алгебраической геометрием.
Рассмотрим алгебры над полем K. Будем считать, что 1 ∈ A
(в противном случае ее можно добавить внешним образом).
Свободный объект в этом классе -- алгебра многочленов над K
K[x1, x2, ..., x_n]
Пусть A -- произвольная алгебра над K.
f: X --> A -- произвольное отображение f(x_i) = a_i
то f продолжается до гомоморфизма алгебр.
f: K[X] --> A
Пусть A -- конечно порожденная алгебра A = K<a1, a2, ..., a_n>
По теореме о гомоморфизме
A = K[X]/Ker(f)
Ker(f) = { p ∈ K[X]: f(p) = 0 } -- это идеал в алгебре многочленов
I ⊲ K[X] -- подалгебра, выдерживающая умножения
на любые элементы g ∈ K[X]
Получаем, что любая коммутативная к.п. алгебра есть фактор-алгебра
алгебры многочленов по некоторому идеалу.
A ≅ K[X]/I
Изучение алгебр сводится к изучению идеалов в алгебре многочленов.
Основная идея алгебраической геометрии -- рассмотрение аффинных
алгебраических ногообразий, соответствующих идеалам в алгебре многочленов.
I -- идеал в K[x1, ..., x_n]
Аффинное алгебраическое многообразие S --
это множество коней (или "нулей") идеала I
S = { (s1, s2, ..., s_n) ∈ K^n: для любого f(x1, ..., x_n) ∈ I
f(s1, ..., s_n) = 0 }
Две теоремы Гильберта
1. Теорема Гильберта о конечности базисов
Алгебра многочленов над произвольным полем нетерова,
т.е. любой идеал в ней порождается (как идеал) конечным числом
элементов.
Эквивалентное определение нетеровости -- любая строго возрастающая
цепочка идеалов имеет конечную длину.
Это означает, что, если у нас есть бесконечная система
алгебраических уравнений, то из нее всегда можно выбрать
конечную подсистему уравнений, имеющую то же множество решений.
2. Теорема Гильберта о нулях.
Классическая формулировка.
Пусть поле K алгебраически замкнуто
(например, K = C -- поле комплексных чисел). Тогда
многочлен f обращается в 0 на всех корнях идеала I <==>
f^m принадлежит I для некоторого целого m > 0.
Определение радикала идеала I
Rad(I) = { f: f^m ∈ I для некоторого целого m > 0 }
I ⊂ Rad(I)
Множество многочленов, обращающихся в ноль на корнях идеала I,
совпадает с радикалом идеала I.
Аффинные многообразия соответствуют идеалам J таким, что Rad(J) = J.
Более важна упрощенная формулировка теоремы Гильберта о нулях:
пусть поле K алгебраически замкнуто,
I -- собственный идеал в алгебре многочленов K[x1, ..., x_n]
(т.е. I не содержит единицу). Тогда I имеет хотя бы один корень.
Доказательство.
Если I не максимальный идеал, то погрузим I в максимальный идеал M
(это можно сделать ввиду нетеровости кольца многочленов).
Фактор-алгебра K[X]/M -- это поле, в котором элемент
(x1', x2', ..., x_n') -- корень идеала I.
Осталось показать, что это это поле совпадает c K.
Лемма (Гаусса, Зарисского, Нетер?)
Пусть расширение F поля K конечно порождено как алгебра над K.
K ⊂ F, F -- поле, F = K[a1, a2, ..., a_m]
Тогда F есть конечное расширение поля K.
Поскольку поле K алгебраически замкнуто, у него не может быть
конечных расширений, значит, K[X]/M == K. Теорема доказана.
Предложение.
Пусть K -- алгебраически замкнутое поле,
I ⊲ K[x1, ..., x_n] -- собственный идеал, I = Ideal{g1, g2, ..., g_m}
Тогда
f обращается в 0 на корнях идеала I <==>
Ideal{g1, g2, ..., g_m, 1 - y f} содержит единицу в K[x1, ..., x_n, y]
Доказательство.
==> Ideal{g1, g2, ..., g_m, 1 - y f} не может иметь корней =>
он тривиален.
<== применяется "трюк Рабиновича"
Пусть 1 ∈ Ideal{g1, g2, ..., g_m, 1 - y f}
1 = u1*g1 + ... + u_m*g_m + v*(1 - y*f), u_i, v ∈ K[X, y]
Подставим y = 1/f в поле рациональных дробей (от многочленов).
Получим
1 = w1*g1 + ... + w_m*g_m, где w_i = u_i(y = 1/f)
w_i -- дробь, в знаменателе которой стоит некоторая степень f.
Домножим на максимальную степень f.
f^s = w1'*g1 + ... + w_m'*g_m, где w_i' -- уже многочлен от x_i
Получили, что f^s ∈ Ideal{g1, ..., g_m}
Следовательно, f обращается в 0 на корнях идеала I.
Аналогично доказываетсЯ и классическая формулировка
теоремы Гильберта о нулях:
многочлен f обращается в 0 на всех корнях идеала I <==>
s
f принадлежит I для некоторого целого s > 0.
Теория базисов Грёбнера
Стандартные базисы рассматривались в различных алгебраических
структурах, например, в группах.
Ситуация: мы задаем кольцо, алгебру, группу, ...
как фактор свободного кольца, алгебры, группы по заданному
набору соотношений, в случае алгебры многочленов по набору
многочленов, обращающихся в 0.
Нам нужно научиться работать с элементами фактор-алгебры,
в часности, определить канонический вид элементов фактор-алгебры.
K[X]/I I = Ideal{f1, ..., f_m}
f_i = старший член(f_i) + cумма остальных членов
Считаем, что старший коэффициент равен 1.
f_i = lm(f_i) + u_i
lm(f_i) -- старший моном, т.е. слово от переменных
Например,
f = x^3*y*z - (xy + z)
Это означает, что по модулю идеала I
x^3 y z == xy + z
Замена x^3 y z --> xy + z называется редукцией.
Нормальная форма -- когда мономы в многочлене не делатся
на старшие члены многочленов, образующих базис идеала.
Система многочленов называется базисом Гребнера, если
нормальная форма единственна.
Теорема. Для любого идеала в кольце многочленов можно построить базис
Грёбнера.
Предложения.
1. f ∈ I <=> f редуцируется к нулю по базису Грёбнера.
2. I тривиален <=> редуцированный базис гребнера идеала I
состоит из одного элемента 1.
3. Зафиксируем чисто лексикографический порядок на мономах.
x1^s1 * x2^s2 * ... * x_n^s_n
показатель степеней (s1, s2, ..., s_n)
Мономы сравниваются по показателям степеней.
x^3 z > x^2 y^1000 z^100
(3, 0, 1) > (2, 1000, 100)
Тогда система уравнений имеет конечное дискретное множество решений
тогда и только тогда, когда
в базисе Грёбнера идеала для каждой переменной x_i найдется многочлен
со старшим мономом, равным степени этой переменной.