Алгоритмы и структуры данных
http://mech.math.msu.su/~vvb/
Что требуется:
1) С++ компилятор,
2) make,
3) для оконных приложений и графики
нужно установить QT (кьют) + QtCreator
free-вариант над MinGW-64
Я советую под MS Windows установить
систему CygWin (для неграфических, консольных
приложений). В cygwin надо установить следующие пакеты:
1) gcc-g++
2) make
3) gdb
4) unzip, zip
5) wget,
6) ssh, ...
7) python3
В MS Windows есть оболочка PowerShell или
древняя CMD, но они очень неудобные.
В cygwin оболочка bash полностью такая же, как в Unix
(Linux) и в OS-X (Apple -- это тоже вариант Unix).
Под MS Windows можно установить Linux как вторую
систему, лучше всего WSL2 (Windows System for Linux),
которая эмулирует Linux Ubuntu.
Пришлите мне список студентов
на адрес vladibor2016@yandex.ru
Используем С++ со стандартной библиотекой STL
Структура данных, или ADT (Abstract Data Type) --
классы, которые организуют хранение данных и работу
с ними (поиск, добавление, удаление, циклы "для каждого"
и т.д.). Структура данных в STL С++ обычно называется
контейнером, контейнеры зависят от параметров
(тип элементов в структуре, размер, специальные функции
и т.п.), которые задаются с помощью механима template
(шаблоны).
Библиотека STL содержится в пространстве имен
namespace std;
#include <iostream>
#include <vector>
#include <string>
#include <map>
std::vector<double> a(100);
std::map<std::string, int> notebook;
1. Вспомним еще раз классы
Используя классы, мы
1) пользуемся готовыми библиотеками классов,
написанными другими людьми, для решения
наших задач;
2) создаем свои собственные классы.
Как пример, рассмотрим геометрию на плоскости R2
и классы для решения геометрических задач.
У нас есть небольшая библиотека R2Graph,
реализующая геометрию на плоскости. В ней
есть классы R2Vector и R2Point
R2Vector u, v, w;
R2Point p, q, t;
v = q - p; // вектор из точки p в точку q
t = p + w; // точка t получается из точки p
// сдвигом на вектор w
Задача 1: найти расстояние от точки t
до прямой (p, v)
(прямая задается точкой p и
направляющим вектором v).
R2Point t;
R2Point p; R2Vector v;
. . .
R2Vector n = v.normal();
n.normalize();
double d = (t - p)*n;
double dist = fabs(d);
Задача 2: найти точку пересечения двух прямых.
Дано: две прямые (p1, v1) и (p2, v2).
Вычислить точку q пересечения этих двух
прямых.
Поскольку q лежит на первой прямой, то
q = p1 + v1*t, где t <- R
Пусть n -- нормаль ко второй прямой
R2Vector n = v2.normal()
Тогда, т.к. q лежит на второй прямой, то
вектор q - p2 перпендикулярен нормали n:
(q - p2)*n = 0
Подставим q из первого уравнения
(p1 + v1*t - p2)*n = 0
(p1 - p2)*n + t*(v1*n) = 0
t*(v1*n) = (p2 - p1)*n
t = ((p2 - p1)*n) / (v1*n)
R2Point intersectLines(
const R2Point& p1, const RVector& v1,
const R2Point& p2, const R2Vector& v2
) {
R2Vector n = v2.normal();
double t = ((p2 - p1)*n) / (v1*n);
return p1 + v1*t;
}
Задача:
дан треугольник, заданный тремя точками a, b, c.
Найти вписанную окружность (центр и радиус).
Вот программа, решающая эту задачу -- файл "incircle.cpp".
Обратите внимание, что мы используем только методы
классов вектор и точка и не используем вычислений
с координатами.
#include <iostream>
#include "R2Graph.h"
void inCircle(
const R2Point& a,
const R2Point& b,
const R2Point& c,
R2Point& center,
double& radius
);
int main() {
std::cout << "Enter 3 vertices of triangle:"
<< std::endl;
R2Point a, b, c;
std::cin >> a >> b >> c;
R2Point center;
double radius;
inCircle(a, b, c, center, radius);
std::cout << "Inscribed circle:" << std::endl;
std::cout << "center = " << center
<< ", radius = " << radius << std::endl;
return 0;
}
void inCircle(
const R2Point& a,
const R2Point& b,
const R2Point& c,
R2Point& center,
double& radius
) {
R2Vector ab = b - a;
ab.normalize();
R2Vector ac = c - a;
ac.normalize();
R2Vector bisa = ab + ac;
R2Vector ba = a - b;
ba.normalize();
R2Vector bc = c - b;
bc.normalize();
R2Vector bisb = ba + bc;
intersectStraightLines(
a, bisa,
b, bisb,
center
);
radius = center.distanceToLine(a, b - a);
}
При компиляции этой программ мы подключаем также
библиотеку R2Graph (файлы "R2Graph.h" и "R2Graph.cpp"),
определяющую классы R2Vector и R2Point и некоторые функции
для решения элементарных геометрических задач на плоскости
R2 (например, вычисление точки пересечения двух прямых
или нахождение расстояния от точки до прямой).
Команда компиляции в Linux или в CygWin следующая:
g++ -o incircle incircle.cpp R2Graph.cpp
Вот пример работы этой программы:
/home/vvb/Dush2025>./incircle
Enter 3 vertices of triangle:
0 0
4 0
4 3
Inscribed circle:
center = (3, 1), radius = 1
------------------------------------
26 Feb 2025
Замечательные точки треугольника
Мы знаем по школьному курсу следующие точки тругольника:
1) центр вписанной окружности (пересечение биссектрис);
2) центр описанной окружности (circumcenter =
точка пересечения серединных перпендикуляров);
3) точка пересечения медиан треугольника (центр тяжести);
4) точка пересечения высот треугольника (ортоцентр).
Теорема Чевы (Ceva)
Три отрезка из вершин треугольника к противоположным
сторонам пересекаются в одной точке тогда и только тогда,
когда выполняется равенство:
x/y * z/t * u/w = 1
где x, y, z, t, u, w -- длины отрезков, на которое
делятся стороны треугольника, в порядке обхода против
часовой стрелки.
Другие популярные популярные точки (пересечения чевиан):
5) точка Жергона;
6) точка Нагеля;
7) точка Лемауна.
Изогональное сопряжение
Пусть x -- точка внутри треугольника.
Проведем через нее 3 чевианы.
Потом каждую из чевиан отразим симметрично относительно
биссектрисы, проведенной из той же вершины.
Три симметричных отрезка, полученные из чевиан, тоже
будут пересекалься в одной точке y. Точка y изогонально
симметрична исходной точке x.
Точка Жергона -- это пересечение чевиан, проведенных
к точкам касания вписанной окружности сторон треугольника.
Точка Нагеля -- это пересечение чевиан, проведенных
к точкам касания трех внешне-вписанных окружностей.
Точка Лемуана -- это точка, изогонально симметричная
центру тяжести треугольника, т.е. пересечение трех
симмедиан. Симмедиана -- это линия, зеркально симметричная
медиане относительно биссектрисы, выходящей из той же
вершины треугольника.
Замечательные свойства центра тяжести и точки Лемуана:
1) центр тяжести минимизирует сумму квадратов расстояний
от точки до вершин треугольника;
2) точка Лемуана минимизирует сумму квадратов расстояний
от точки до сторон треугольника.
Для примера решим задачу 2:
вычислить точку Жергона треугольника.
Функция, вычисляющая точку Жергона, использует
функцию inCircle, вычисляющую окружность, вписанную
в треугольник (ее центр и радиус), которую мы
уже реализовали ранее.
R2Point gergonne(
const R2Point& a,
const R2Point& b,
const R2Point& c
) {
R2Point inCenter;
double radius;
inCircle(a, b, c, inCenter, radius);
// std::cout << "Incircle center: " << inCenter << std::endl;
R2Point cev_a;
intersectStraightLines(
inCenter, (c - b).normal(),
b, c - b,
cev_a
);
// std::cout << "cev_a: " << cev_a << std::endl;
R2Point cev_b;
intersectStraightLines(
inCenter, (a - c).normal(),
c, a - c,
cev_b
);
// std::cout << "cev_b: " << cev_b << std::endl;
R2Point g;
intersectStraightLines(
a, cev_a - a,
b, cev_b - b,
g
);
return g;
}
Вот пример работы этой программы:
/home/vvb/Dush2025>g++ -o gergonne gergonne.cpp R2Graph.cpp
/home/vvb/Dush2025>./gergonne
Enter 3 vertices of triangle:
0 0
5 1
3 5
Gergonne point: (3.11423, 1.99903)
--------------------------------------
12 Mar 2025
Классы в C++
Прошлы раз мы рассмотрели реализацию класса Complex,
представляющего комплексные числа.
Внутреннее прелставление комплексного числа
возможно в двух вариантах:
1) объект хранит вещественную и мнимую части
комплексного числа
class Complex {
private:
double re;
double im;
. . .
};
2) объект хранит абс. величину (модуль) и аргумент
(угол вектора с осью абсцисс) комплексного числа,
такое представление называется экспоненциальным
z = r*e^{i phi}
или полярным.
class Complex {
private:
double r;
double phi;
. . .
};
Каждое из этих представлений имеет свои плюсы и минусы.
В декартовом представлении удобно реализуются операции
+ и -, в экспоненциальном -- операции * и /.
Например, в полярном представлении оператор +=
реализуется следующим образом:
class Complex {
private:
double r;
double phi;
public:
double real() const {
return r*cos(phi);
}
double imag() const {
return r*sin(phi);
}
Complex& operator+=(const Complex& z) {
double re = real() + z.real();
double im = imag() + z.imag();
r = sqrt(r*r + im*im);
phi = atan2(im, re);
return *this;
}
Complex& operator*=(const Complex& z) {
r *= z.r;
phi += z.phi;
phi = fmod(phi, 2.*M_PI);
if (phi < -M_PI)
phi += 2.*M_PI;
else if (phi > M_PI)
phi -= 2.*M_PI;
return *this;
}
. . .
};
Добавим в класс Complex несколько статических методов.
Обычные методы применяются к объектам класса, например
Complex z;
double = z.real();
Статические методы не применяются к объектам, это просто
как бы функции, лишь номинально причисленные к классу.
В классе Complex есть конструктор, который по двум
вещественным числам создает объект, у которого комплексная
и мнимая части равны двуи переданным числам.
double a, b;
. . .
Complex z = Complex(a, b);
double r, phi;
// Как создать комплексное число тоже по двум
// вещественным числам, но которые определяют не
// декартовы координаты числа, а модуль и аргумент?
// Используем статический метод polar
Complex z2 = Complex::polar(r, phi);
Реализация этого статического метода внутри класса:
class Complex {
private:
double re;
double im;
public:
. . .
static Complex polar(double r, double phi) {
return Complex(r*cos(phi), r*sin(phi));
}
};
Для примера применения класса Complex
решим кубическое уравнение:
a*x^3 + b*x^2 + c*x + d = 0
Можно считать, что a = 1, иначе просто разделим
уравнение на a.
x^3 + a*x^2 + b*x + c = 0
Можно свести уравнение к виду, когда коэффициент
при x^2 равен нулю. Иначе делаем замену типа сдвига:
x = x' + s
Получим
(x' + s)^3 + a*(x' + s)^2 + b(x' + s) + c = 0
x'^3 + 3*sx'^2 + ... + a*x'^2 + ... = 0
x'^3 + (3*s + a)*x'^2 + ... = 0
откуда получаем
s = -a/3
То есть при замене
x = x' - a/3
член при х'^2 сокращается.
Окончательно уравнение примет вид
x^3 + a*x + b = 0
Мы свели исходную задачу к уравнению такого вида.
Единственное, что нужно помнить для вывода формулы Кардано --
это замена
x = u + v
где u, v -- две независимые переменные.
(Свели простое уравнение от одной переменной к
более сложному уравнению от двух переменных.)
Подставим и раскроем скобки:
(u + v)^3 + a*(u + v) + b = 0
Считаем, что 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 + (u + v)*(3*u*v + a) = 0
============= --------------------
Приравняем отмеченные части к нулю и решим систему
из двух уравнений с двумя неизвестными:
{ u^3 + v^3 + b = 0
{ (u + v)*(3*u*v + a) = 0
Поскольку b != 0 и x = 0 не является решением, x = u + v,
то можем сократить второе уравнение на (u + v):
{ u^3 + v^3 + b = 0
{ 3*u*v + a = 0
Из второго уравнения выражаем v:
v = -a/(3*u)
и подставляем в первое уравнение:
u^3 - a^3/(27*u^3) + b = 0
Сделаем замену
u^3 = y
Получим
y - a^3/(27*y) + b = 0
Сведем к квадратному, домножив на y:
y^2 + b*y - a^3/27 = 0
Дискриминант уравнения:
d = b*b + 4*a^3/27
Достаточно получить 1 корень:
y = (-b + sqrt(d))/2
Переходя обратно, получим 3 корня исходного
уравнения, поскольку уравнение
u^3 = y
имеет ровно 3 комплексных решения (3 кубических корня
из y).
Московское время 13:45.
Перерыв до 13:55
-----------------------------------
Мы реализовали вычисление трех корней кубического уравнения
в файле cubicEq.cpp
Тест:
рассмотрим уравнение
(x - 2)*(x - 3)*(x + 5) = 0
раскроем скобки:
(x^2 - 5*x + 6)*(x + 5) = 0
x^3 - 5*x^2 + 6*x + 5*x^2 - 25*x + 30 = 0
======= =====
x^3 - 19*x + 30 = 0
Корни уравнения: 2, 3, -5.
Работа нашей программы:
Solve a cubic equation of form
x^3 + a*x + b == 0
with complex coefficients
Enter complex a, b:
-19 0
30 0
Solutions of the cubic equation:
3 + 4.44089e-16*i
-5 + 0*i
2 + 0*i
Ответ правильный!
Еще пример: решим уравнение
x^3 - x - 1 = 0
Solve a cubic equation of form
x^3 + a*x + b == 0
with complex coefficients
Enter complex a, b:
-1 0
-1 0
Solutions of the cubic equation:
1.32472 + 0*i
-0.662359 + 0.56228*i
-0.662359 + -0.56228*i
Ответ также правильный.
Проверим с помощью SageMath:
sage: res = solve(x^3 - x - 1, x)
sage: res
[x == -1/2*(1/18*sqrt(23)*sqrt(3) + 1/2)^(1/3)*(I*sqrt(3) + 1)
- 1/6*(-I*sqrt(3) + 1)/(1/18*sqrt(23)*sqrt(3) + 1/2)^(1/3),
x == -1/2*(1/18*sqrt(23)*sqrt(3) + 1/2)^(1/3)*(-I*sqrt(3) + 1)
- 1/6*(I*sqrt(3) + 1)/(1/18*sqrt(23)*sqrt(3) + 1/2)^(1/3),
x == (1/18*sqrt(23)*sqrt(3) + 1/2)^(1/3) +
1/3/(1/18*sqrt(23)*sqrt(3) + 1/2)^(1/3)]
sage: res[0].rhs()
-1/2*(1/18*sqrt(23)*sqrt(3) + 1/2)^(1/3)*(I*sqrt(3) + 1) - 1/6*(-I*sqrt(3) + 1)/(1/18*sqrt(23)*sqrt(3) + 1/2)^(1/3)
sage: res[0].rhs().n()
-0.662358978622373 - 0.562279512062301*I
sage: res[1].rhs().n()
-0.662358978622373 + 0.562279512062301*I
sage: res[2].rhs().n()
1.32471795724475
Ответы совпали.
В следующий раз мы рассмотрим еще несколько классов:
матрицы, треугольник на плоскости, прямоугольник на
плоскости, многоугольник (или контур) на плоскости.
Также рассмотрим, как надо использовать контейнеры
из стандартной библиотеки классов STL языка C++
и почему не надо использовать базовые конструкции
языка C/C++, а использовать именно контейнеры.
Также будут даны задачи на классы.
------------------------------------------------
Zoom-сессия в воскресенье 16 марта в 11:00
московского времени.
------------------------------------------------
19 Mar 2025
Классы в C++
Класс Матрица
Неправильное решение:
матрица размера m*n (m строк, n столбцов)
хранится как массив строк.
Каждая строка захватывается по-отдельности
в динамической памяти и представляет собой массив длины n.
Кроме того, хранится массив указателей на начала строк
длины m.
double **a = new double*[m];
for (int i = 0; i < m; ++i) {
a[i] = new double[n];
}
// После этого элемент матрицы с индексами i, j
// записывается как a[i][j]
// Удаление матрицы (освобождение памяти):
for (int i = 0; i < m; ++i) {
delete[] a[i];
}
delete[] a;
Это дурацкий способ работы с матрицами.
Почему он плохой?
1. Чтение/запись требует двух обращений к памяти:
a[i][j]
сначала читаем p = a[i] -- указатель на начало
i-й строки;
затем читаем сам элемент p[j].
2. Плохо используется кэш-память процессора.
Выгоднее располагать элементы матрицы в непрерывном
сегменте памяти.
Если вы будете использовать эту реализацию, то я не
буду засчитывать подобные решения.
--------------------------------------------------
Правильный способ реализации матрицы:
храним элементы матрицы в непрерывном линейном массиве
по строкам (в Фортране по столбцам, в Python
в модуле numpy возможны оба варианта, по умолчанию
используется вариант по строкам).
Пример: матрица размера 3*4 (m = 3, n = 4)
a00 a01 a02 a03
a10 a11 a12 a13
a20 a21 a22 a23
хранится в линейном массиве в следующем порядке:
a00 a01 a02 a03 a10 a11 a12 a13 a20 a21 a22 a23
--------------- =============== ---------------
строка 0 строка 1 строка 2
Соответственно элемент матрицы с индексами i, j
записывается как a[i*n + j].
Так хранятся, к примеру, изображения в видео-памяти
компьютера.
Напишем 2 варианта реализации матрицы.
Первый вариант соттветствует версии языка C++
1985-го года, до того, как была введена стандартная
библиотека STL языка C++ в 1998 г.
===============================================
26 Mar 2025
Если мы явно захватываем какие-то ресурсы при создании
объекта, то обязательно надо
1) определить copy-конструктор;
2) определить оператор копирования
Matrix& operator=(const Matrix& matr);
3) определить move-конструктор
Matrix& operator=(Matrix&& matr);
4) определить оператор перемещения (move)
Matrix& operator=(Matrix&& matr); // r-value
5) определить деструктор
~Matrix();
Это так называемое "правило 5-ти" (раньше оно называлось
"правиле трех", когда в C++ не было move-семантики).
Еще есть в C++ правило RAII --
Resouce Allocation Is Initialization:
захват ресурсов удбно выполнять в конструкторе,
освобождение в деструкторе.
Большинтсво современных не слишком квалифицированных
программистов на C++ этих правил не знают, потому что
они вместо базовых массивов языка C/C++ используют
конетейнерные классы из стандартной библиотеки.
Для этого, например, мы используем вместо базовых
и указателей на их начала мы используем класс
"динамический массов" из стандартной библиотеки STL.
Например, наш класс Матрица вместо старой реализации:
class Matrix {
private:
int m; // Number of rows
int n; // Number of columns
double *a;
. . .
};
используем правильную" реализацию с современной точки
зрения (читайте книгу Б.Страуструпа
"Язык программирования C++"). Вместо массива, который
задается указателем на его начальный элемент
double *a;
мы используем класс std::vector<double> из стандартной
библиотеки.
#include <vector>
class Matrix {
private:
int m; // Number of rows
int n; // Number of columns
// double *a;
std::vector<double> a;
. . .
};
Что нам это дает?
!!! Не нужно определять эти 5 методов, поскольку
компилятор сгенерирует их автоматически и они будут
работать правильно, поскольку для всех контейнерных
классов эти операции уже определены и работают
правильно (в отличие от копирования указателей, что
всегда приводит к катастрофе).
Рассмотрим самый важный класс из стандартной библиотеки:
динамический массив.
std::vector<double> a;
+--------------------------+-----------------+
| size() |x *** мусор **** |
+--------------------------+-----------------+
capacity()
Добавление элемента в конец динамического массива:
a.push_back(x);
Мнемоника:
начало структуры данных называется front()
конец структуры данных называется back()
Добавить в конец контейнера a:
a.push_back(x);
Удалить из конца контейнера a:
a.pop_back();
Удалить из начала контейнера a:
a.pop_front();
Методы begin() и end() возвращают вовсе не
элементы в начале и в конце контейнера!!!
Они возвращают ИТЕРАТОРЫ, указывающие на начало
контейнера и на воображаемый элемент, который следует
за последним элементом контейнера (я обычно
называею этот элемент ЗАКОНЕЦ). Таким образом,
итераторы a.begin() и a.end() задают как бы
полуинтервал, в котором находятся все элементы
контейнера.
Продолжим реализацию класса Matrix, используя
современный подход: вместо базовых массивов
языка C/C++ используем контейнерные классы из
стандартной библиотеки STL. (Есть также
библиотека Boost.)
. . .
В методе Гаусса мы на очередном шаге обрабатываем
минор матрицы, начинающийся со строки i и столбца j.
В чистой математике мы ищем ненулевой элемент в
столбце j, начиная с i-й строки. Если такого нет, то
столбец нулевой, и мы просто увеличиваем j на 1.
Если есть ненулевой элемент в строке k, то мы меняем
строки i, k так, чтобы ненулевой элемент стоял бы в
начале минора в i-й строке.
НО!!! ТАК ДЕЛАТЬ НЕЛЬЗЯ при реальных вычислениях.
ПОЧЕМУ ТАК???
Когда потом зануляем столбец, мы вычитаем из
строки k, где k > i, строку i, умноженную на коэффициент
a[k][j]/a[i][j],
так, чтобы первый элемент в строке k обнулился.
Векторная запись:
a[k] = a[k] - a[i] * (a[k][j]/a[i][j])
Если элемент a[i][j] не является максимальным по модулю
в j-м столбце, то коэффициент,
на который умножаем,
| a[k][j]/a[i][j] | > 1
и ошибки экспоненциально растут.
Следовательно, надо искать не просто ненулевой элемент,
а максимальный по модулю элемент в столбце необработанного
минора.
---------------------
9 Apr 2025
Продолжение: классы в C++. Некоторые алгоритмы
1. Мы рассмотрели класс матриц над полем R
и реализовали алгоритм Гаусса приведения матрицы
к ступенчатому виду (Row Echelon Form).
Рассмотрим более подробно
1) алгоритм решения системы линейных уравнений
в частном случае, когда матрица системы невырождена;
2) алгоритм вычисления обратной матрицы.
Основаны оба этих алгоритма на приведении произвольной
матрицы к редуцированному ступенатому виду
(Reduced Row Echelon Form):
Обычный ступенчатый вид:
0 0 ... 0 x * * * * *
0 0 ... 0 0 0 x * * *
0 0 ... 0 0 0 0 x * *
0 0 ... 0 0 0 0 0 0 x
0 0 ... 0 0 0 0 0 0 0
. . .
0 0 ... 0 0 0 0 0 0 0
где x != 0, * -- любое число.
Редуцированный ступенчатый вид:
0 0 ... 0 1 * 0 0 * 0
0 0 ... 0 0 0 1 0 * 0
0 0 ... 0 0 0 0 1 * 0
0 0 ... 0 0 0 0 0 0 1
0 0 ... 0 0 0 0 0 0 0
. . .
0 0 ... 0 0 0 0 0 0 0
Редуцированный ступенчатый вид получается
из обычного обратным проходом.
Идем снизу вверх, находим первую ненулевую
строку i, домножаем ее на элемент, обратный
к первому ненулевому элементу в строке a_ij,
после чего a_ij = 1.
Далее из всех строк с индексом k > i выше нее
вычитаем i-ю строку, умноженную на элемент a_kj
так, чтобы a_kj стал равен нулю.
Напишем это на C++ (добавим метод в классе матриц).
================================
16 Apr 2025