Компьютерный практикум (магистратура)
Черновики лекций, осенний семестр 2021


Программа курса

1. Осенний семестр -- С++, классы, стандартная библиотека STL, 
   алгоритмы (сортировка, поиск, алгоритм Дейкстры поиска
   минимального пути в графе), структуры данных (vector, deque, set,
   map, list, priority_queue ???), работа с текстами в Unicode/UTF-8,
   графическое/оконное программирование в Qt
   
   Зачет: надо сдать по одной задаче на каждую тему в течение семестра,
   для задач будет deadline. Задачи надо присылать мне на адрес
       vladimir_borisen@mail.ru
       Subject: Магистратура ...
   http://mech.math.msu.su/~vvb/
   http://mech.math.msu.su/~vvb/Master/   -- материалы для магистратуры
   
2. Весенний семестр: Python3, SageMath -- система компьютерной математики
   (компьютерной алгебры), свободная система (оболочка для множества свободных
   математических пакетов) + надстройка над языком Python
   
Язык С++. Классы и объектно-ориентированный стиль программирования

Решение несложных геометрических задач (на плоскости R2
или в пространстве R3) с помощью классов вектор и точка
      R2Vector, R2Point на плоскости
      R3Vector, R3Point в пространстве
2 типа задач:
1) воспользоваться готовыми чужими классами для решения своей задачи;
2) разработка класса.    
   
Методы класса вектор
   R2Vector u, v, w;
   w = u + v;
   w = u - v;
   w = u*1.5;
   double s = u*v; // скалярное произведение
   R2Vector n = u.normal();
   assert(n*u == 0.);
   double alpha = u.angle(v);
   R2Point p, q;
   w = p - q;
   q = p + w;
   double d = p.distance(q);
   d = (q - p).length();
   d = (q - p).norm();
   
   u.normalize(); // Поделить вектор на его длину
   w = v.normalized();
   u.area(v);
   
Требование к первой задаче: решить ее, не использую координат векторов
и точек, используя только методы классов (чисто геометрическое решение).

Простые примеры:
1) найти расстояние от точки q до прямой.
   Прямая задается точкой p и направляющим вектором v.

    R2Vector n = v.normal();    // v = (x, y)  n = (-y, x)
    n.normalize();
    double d = (q - p)*n;

2) найти точку пересечения двух прямых (мы знаем, что они не параллельны).

    n -- нормаль к вектору v1
    уравнение второй прямой: (q - p1)*n = 0
    Так как q лежит на первой прямой, то
    q = p0 + v0*t
    Подставляем выражение для q в уравнение второй прямой:
    (p0 + v0*t - p1)*n = 0
    (p0 - p1)*n + t*(v0*n) = 0
    t = (p1 - p0)*n / (v0*n)
    
Фрагмент программы на C++
    R2Point p0; R2Vector v0;    // Первая прямая
    R2Point p1; R2Vector v1;    // Вторая прямая
    . . .
    R2Vector n = v1.normal();
    double t = (p1 - p0)*n / (v0*n);
    q = p0 + v0*t;
    
Найти угол между векторами:
    R2Vector u, v;
НЕПРАВИЛЬНОЕ РЕШЕНИЕ!!!
Воспользоваться скалярным произведением, найти косинус угла
между векторами и применить функцию arccos.
Плохо, потому что arccos(x) вблизи x = 1 ведет 
себя очень плохо, производная стремится к +- infinity
Надо использовать функцию atan или, лучше, atan2

    R2Vector n = u.normal();
    double alpha = atan2(v*n, v*u);
    
    
ПЕРЕРЫВ 10 минут до 10:48   

Как рекомендуется работать

Если Linux / Apple MAC, то используем gcc/g++,
компилируем в консоле/терминале и исполняем приложение
из командной строки.

В MS Windows 10 советую либо установить WSL2
(Windows System for Linux), устанавливаем Linux Ubuntu,
либо (более простой путь) установить CygWin
(устанавливаются все свободные программы, существующие 
в мире, плюс терминал (консоль) bash).
Надо обязательно выбрать
   -- gcc/g++
   -- make
   -- vi (vim)
   -- wget
   -- TeX/LaTeX
   -- ...
Я не советую пользоваться оболочками для разработки,
кроме QtCreator. Оболчки: CodeBlocks, Visual C++, ...
Лучше сначала просто набирать команды компиляции
в командной строке или использовать make + Makefile.  

Для редактирования файлов под MS Windows следует
использовать Notepad++  -- очень хороший текстовый редактор,
только обязательно установите замену табуляций на пробелы (4 пробела).

Для Linux есть редактора gedit, kate, mousepad,
для Apple MAC -- TextWrangler или что-другое.

Задача:
дан треугольник. Вычислить вписанную окружность
(центр и радиус). Надо решить задачу, используя только
методы классов R2Vector и R2Point.
Эти классы реализованы в файлах R2Graph.h, R2Graph.cpp.

=======================================

Для вычислений в трехмерном пространстве и 3D-графики
удобно использовать классы R3Vector и R3Point.

С помощью Google или Yandex-карты можно получить
координаты точек. Для вычисления длины дистанции
достаточно уметь вычислять расстояния между точками на
карте, точки заданы их координатами.

Как найти угол между векторами в R3?

    R3Vector u, v;
    
    R3Vector ex = u;
    ex.normalize();
    R3Vector ez = u.vectorProduct(v);
    ez.normalize();
    R3Vector ey = ez.vectorProduct(ex);
    double x = v*ex;
    double y = v*ey;
    double alpha = atan2(y, x);
    
==================================

10 Sep 2021

Объектно-ориентированный стиль программирования подходит не для
всех задач. Чисто инженерные и зачастую математические задачи
(типа решения уравнений с частными производными) удобнее решать,
оставаясь в рамках чисто процедурного программирования.

Но есть области, где без объектно-ориентированного подхода не
обойтись: прежде всего, это оконное и графическое программирование.

Пример: задачи по аналитеской геометрии (графике на плоскости 
R2 или в пространстве R3) удобно решать, используя объекты
вектор (R2Vector, R3Vector), точка (R2Point, R2Point), треугольник,
многоугольник (R2Polygon), триангуляция, ...

Пример геометрических задач, которые удобно решать с помощью этих классов:
замечательные точки треугольника.
   центры вписанной, описанной окружностей;
   перечсечения медиан, высот, разные вариции на теоремы Чевы
    точка Жергона, Нагеля, Лемуана, Ферма-Торричелли, ...
Еще больший выигрыш дает использование классов R3Vector, R3Point
в трехмерных задачах.

Пример -- задача, которая у меня встречалась на практике:
1) вычислить расстояние между двумя точками на поверхности Земли
   (точки задаются сферическими координатами: (широта, долгота))
   
   d = v1.angle(v2)*R

Как вычислить угол между векторами v1 и v2?
На плоскости:
плохое решение
   s = v1*v2
   c = s/(v1.norm()*v2.norm())
   alpha = acos(c)
правильное решение -- использовать arctg (atan в англо-язычной традиции)
   n = v1.normal()        // (x, y) -> (-y, x)
   alpha = atan2(v2*n, v2*v1)
    
В R3 есть замечательная операция векторного произведения:

   | i  j  k  |
   | x1 y1 z1 | = i*(y1*z2 - y2*z1) + j*(-x1*z2 + x2*z1) + k*(x1*y2 - x2*y1)
   | x2 y2 z2 |
   
   v1 = R3Vector(x1, y1, z1)
   v2 = R3Vector(x2, y2, z2)
   w = v1.vectorProduct(v2)
   w = R3Vector(
      y1*z2 - y2*z1, -x1*z2 + x2*z1, x1*y2 - x2*y1
   );
   
Она позволяет вычислить нормаль к плоскости.
Пусть плоскость задана тремя точка (или точкой и двумя векторами)
      p0, p1, p2
      n = (p1 - p0).vectorProduct(p2 - p0)

Уравнение плоскости:
      n*(q - p0) = 0

      a*x + b*y + c*z + d = 0
      (a, b, c) -- нормаль к плоскости

Также оно позволяет вычислить площадь параллелограмма (или треугольника),
натянутого на два вектора в пространстве.

Вычисление угла между векторами.
v1, v2
Идея: постороит ортонормированный репер в плоскости (v1, v2)
      n = v1.vectorProduct(v2)
      n.normalize()
      ez = n
      ex = v1.normalized()
      ey = ez.vectorProduct(ex)
      alpha = atan2(v2*ey, v2*ex)
      
      alpha = v1.angle(v2)
      
Расстояние между точками на поверхности Земли: earthdist.cpp

Мы реализовали функцию
    R3Vector radiusVector(double lat, double lon);
которая возвращает радиус-вектор единичной длины, направленный
в точку с координатами (lat, lon).

Для другой задачи нам понадобится обратная функция:
по радиус-вектору v единичной длины получить координаты (lat, lon)
Прототип функции:
    void earthCoordinates(R3Vector v, double& lat, double& lon); // Плохо 
                                                        // (но приемлемо)
    void earthCoordinates(R3Vector& v, double& lat, double& lon); // ОЧЕНЬ 
                                                                  // Плохо!
Правильно:
    void earthCoordinates(const R3Vector& v, double& lat, double& lon);

    ПЕРЕРЫВ 10 минут до 10:50
      
void earthCoordinates(const R3Vector& v, double& lat, double& lon) {
    double phi = atan2(v.y, v.x);
    lon = phi*180./PI;
    R3Vector ez(0., 0., 1.);
    double alpha = ez.angle(v);
    double theta = PI/2. - alpha;
    lat = theta*180./PI;
}

Две задачи, которые были у меня в реальном проекте.
В проекте изображались самолеты на карте вблизи аэродрома
(на расстоянии <= 200 км).
В стандарте ADBS все самолеты постоянно транслируют информацию о себе,
включая координаты (широта, долгота). Надо было изображать самолеты
на карте => отсюда вытекают две задачи:

1) по координатам центра карты и координатам точки на поверхности
   Земли получить координаты на карте в метрах (x, y);

2) обратная задача -- по координатам точки на карте получить
   координаты на поверхности Земли.
   
Рассматриваем самое простое отображение на карту: центральная проекция
с началом в центре Земли.

============================================

17 Sep 2021

Две типичные ошибки в первой задаче:

#include "R2Graph.cpp"

Большие проекты всегда разбиваются на модули.
Модули транслируются и отлаживаются отдельно.
Для связи модулей в языке C/C++ используются
header-файлы. Они НИКОГДА не содержат выполняемого кода,
а только описания.

double sqrt(double x);

Сам код функции содержится в cpp-файле (или c-файле)
math.cpp или math.c
normal_ab
    orthocenter.cpp, использует также модуль R2Graph.cpp
    
    g++ -c R2Graph.cpp
Результатом будет объектый файл *.o в Unix или *.obj в MS Windows
Код функции не должен дублироваться.

Что может быть в h-файле?

1. Прототипы функций.
2. Описание классов.
3. Задание констант.
4. Определения inline-функций.
5. Описания глобальных переменных.

inline int stringLength(const char* s) {
    const char *p = s;
    while (*p != 0) {
        ++p;
    }
    return (p - s);
}
    
    
     int l = stringLength(s); // Вставка кода функции вместо ее вызова    
    
3. Задание констант.

Болезненный момент -- я всегда всеми силами пытаюсь избежать
использование операторов препроцессора.

Например, так задается константа M_PI в стандартном файле 
/usr/include/math.h

# define M_PI           3.14159265358979323846  /* pi */

Я так никогда не пишу, потому что препроцессор --
атавизм глубокого прошлого.

Если я напишу в языке C в h-файле строку

const double PI = 3.14159265358979323846;

то, к сожалению, компилятор C создаст переменную -- это катастрофа,
потому что h-файл подключается в разные модули многократно, и переменная
с одним и тем же именем будет определена многократно => ошибка компоновщика.

Но в языке C++ строка 
const double PI = 3.14159265358979323846;
интерпретируется правильно -- как задание константы, а не как
задание ячейки памяти, содержащей данное число, которую можно
только читать.

С++ создан в 1985, а современный язык C90 был создан в 1989 под сильным
влиянием C++.

5. Описания глобальных переменных.

В h-файле категорически нельзя писать
int nMax;
В h-файле надо писать
extern int nMax;
а в каком-то ОДНОМ cpp-файле уже задавать эту константу
int nMax = 1000000;

Компиляция для отладчика

    g++ -O0 -g ...
    
Запуск отладчика
    gdb program-file
    
Команды отладчика (debugger):
    r  -- run
    b  -- breakpoint (установить точку остановки)
    n  -- next (выполнить строку)
    s  -- step (нырнуть внутрь функции)
    bt -- back trace (показать стек вызовов)
    f номер -- frame (переместиться по стеку вызовов вверх)
    
----------------------------------------

 Вторая ошибка
    double x, y;
    if (x == 0.) {
        ...
        
    }

Все действия с числами типа double выполняются приближенно.
Поэтому надо писать
const double EPS = 1e-8;

    double x, y;
    // if (x == 0.) {
    if (fabs(x) <= EPS) {
        ...
        
    }
    
    // if (x == y) {
    if (fabs(x - y) <= EPS) {
        ...
    }
    
========================================

Основная тема сегодняшнего занятия (и всего семестра):
классы в C++.

Обзор C++ (что появилось нового по сравнению с C).

Язык C
Типы данных:
   числа 
   1) вещественные числа double, float
      совет: никогда не используйте float!
      (в Python тип Float соответствует типу double в C)
   2) целые числа:
      int
      char -- 8-битовое целое число
      short -- 16-битовое
      int   -- 32-бита
      long == int, 32 бита
      long long -- 64 бита
      signed, unsigned -- модификаторы типа
      (в Java нет unsigned).
   3) в С++ появился тип bool (8 бит).
      int n;
      . . .
      if (n) {  // Не надо так писать!!!
          ...
      }
             
      if (n != 0) {
          . . .
      }
      
Конструкции языка C:
1) указатель -- переменная, которая содержит адрес (и при необходимости
   дополнительную информацию) объекта.
   int n = 123;
   int *p = &n;
   // n и *p -- это одно и то же
   n = 567;
   *p = 567;
   * -- операция de-reference (разыменование)
Арифметика указателей:
   ++p увеличивает указатель на размер одного элемента того типа,
       на который ссылается указатель.
   аналогично уменьшение, разность указателей, прибавление или
   вычитание к указателю целого числа.
   
   double *p;
   p = (double *) 1000;
   int n = (int)(p + 4);
   
Чему равно n?
Ответ: адрес увеличивается на размер 4-х элементов типа double.
       Каждый элемент имеет размер 8 байтов => адрес увеличивается на 32
       (масштабирование указателей).

Связь указателей и массивов:
    double a[10000];
    double *p = a;  // Имя массива -- это указатель
                    //                на его начальный  элемент.       
    // a[123] == *(p + 123)
    // a[123] <--> *(a + 123)
    // можно было писать 123[a]
   
---------------------------------

Что нового появилось в C++

    ПЕРЕРЫВ 10 минут до 10:50 
    
Добавлена конструкция "ссылка" (reference).
Cсылка -- это (логически) второе имя для уже существующего объекта
(псевдоним, alias).
    int n;
    int& k = n;
    // После этого n и k означает одно и то же.
Отличие от указателя: указатель можно перенаправить на другой объект.
Ссылка определяется только один раз при описании и всегда будет указывать
только на тот объект, который был задан при ее описании.
    double a[1000];
    int n;
    double& b = a[2*n + 1];
    b = 12345;  // вместо a[2*n + 1] = 12345;
Ссылки в основном используются для передачи параметров функциям и методам.
Если функция inline, то просто вместо ссылки подставляется тот
объект, на который она ссылается.
Для обычных функций ссылка передается как указатель, только при работе
со ссылкой не нужно использовать операцию разыменования *  

int f(double *p, int& n) {
    *p = 123.456;
    n = 777;    // Не нужно писать *n
    return 0;
}

    double a[100];
    int k;
    f(&(a[10]), k);
    // a[10] == 123.456, k == 777 

    int *p;
    int m, n;
    p = &m;    
    ...
    p = &n;
    
    int& r; // Ошибка! При описании ссылки надо указывать,
            // на что она ссылается.
            // При описании ссылки переменная не создается!
            
-------------------------------------------------------------

Модификатор const.

const int n = 100;
задает константу (если это переменная в памяти, то ее можно только
читать).

Константный указатель можно использовать только для чтения объекта
(объект нельзя испортить с помощью константного указателя).
Аналогично константная ссылка позволяет только читать объект
и не позволяет его менять.

Когда мы пишем прототип функции или метода, надо передавать
входные параметры только через константные указатели или ссылки.

Например, функция, которая находит длину строки, имеет прототип
что-то вроде
    int strlen(const char *s);
(на самом деле
    size_t strlen(const char *s);
    
где size_t определяется скорее всего как
typedef unsigned long long size_t;

Грубая ошибка:
    int stringLength(char *s);
    
    int n = stringLength("abc");

Модификаторы const нужно обязательтно применять при описании
входных аргументов функции (тех аргументов, которые не изменяются
при выполнении функции).

Классы: объекты, имеющие сложное внутреннее устройство и набор методов
для работы с ними.

Давайте для примера напишем 2 класса:
class Complex;  // комплексное число, задаваемое с помощью
                // двух чисел типа double
class Matrix;   // матрица вещественных чисел

Эти классы отличаются тем, что в конструкторе Complex
мы не захватываем никаких ресурсов (в динамической памяти),
все информация содержится прямо внутри объекта класса,
нет указателей в кучу. Следовательно, не нужно явно задавать
copy-конструктор, оператор копирования, деструктор, поскольку
компилятор при их отсутствии сам их генерирует, в данном случае
правильно.

В классе Matrix

class Matrix (
private:
    int m;  // Number of rows
    int n;  // Number of columns
    double *elems;  // Array of matrix elements
                    // a_ij == elems[i*n + j] 

public:
    Matrix(int numRows = 1, int numCols = 1):
        m(numRows),
        n(numCols),
        elems(new double[m*n])
    {
        for (int i = 0; i < m*n; ++i) {
            elems[i] = 0.;
        }
    }
        
    ~Matrix() {
        delete[] elems;
    }
    
    Matrix(const Matrix& a):
        m(a.m),
        n(a.n)
        elems(new double[m*n])
    {
        for (int i = 0; i < m*n; ++i) {
            elems[i] = a.elems[i];
        }
    }
    
    Matrix& operator=(const Matrix& a) {
        . . .       // Must be defined explicitly!!!
        return *this;
    }
    
    . . .
    
};

Б.Страуструп страстно пишет, что ТАК ДЕЛАТЬ НЕЛЬЗЯ!!!
Правильное решение -- использовать стандартную библиотеку!

class Matrix (
private:
    int m;  // Number of rows
    int n;  // Number of columns
    /*... Don't do that!
    double *elems;  // Array of matrix elements
                    // a_ij == elems[i*n + j] 
    ...*/
    std::vector<double> elems;
    
public:
    Matrix(int numRows = 1, int numCols = 1):
        m(numRows),
        n(numCols),
        elems(m*n)
    {
        for (int i = 0; i < m*n; ++i) {
            elems[i] = 0.;
        }
    }
        
       
    http://mech.math.msu.su/~vvb/        
    
===============================================

24 Sep 2021

Была рассмотрена реализация класса Complex.
Реализация полностью inline
Плюсы:
1) программа не тратит время на вызов функций и методов
   и возврат из них;
2) становится возможной оптимизация кода (которая обычно
   ограничивается рамками функции).
Минус:
   раздувается код программы, поскольку одни и те же
   методы многократно повторяются.
   
Таким образом, сложные методы, которые исполняются долго
по сравнению с временем вызова функции.

Для примера, в класс Complex добавим метод нахождения
всех корней степени n из комплексного числа.

Как тест для класса Complex предлагается решить
кубическое уравнение.

Знать формулу Кардано не обязательно, поскольку она очень легко выводится.

1. Кубическое уравнение
      a'*x^3 + b'*x^2 + c'*x + d' = 0
   можно преобразовать к виду
      x^3 + a*x + b = 0
   1) делим на a';
   2) делаем замену типа сдвига x = y + s
      x^3 + b*x^2 + c*x + d = 0
      (y + s)^3 + b*(y + s)^2 + c*(y + s) + d = 0
      y^3 + 3*s*y^2 + b*y^2 + ...
      y^3 + (3*s + b)*y^2 + ...
      s = -b/3
      
2. Считаем, что уравнение имеет вид
    x^3 + a*x + b = 0
    
    Основная идея: делаем замену 
    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 + 3uv(u + v) + a(u + v) = 0
    u^3 + v^3 + b + (u + v)(3uv + a) = 0
    -------------   ================
    
    { u^3 + v^3 + b = 0
    { (u + v)(3uv + a) = 0 т.к. x = u + v и x = 0 -- не корень,
                                то сократим второе уравнение на (u + v)
    
    { u^3 + v^3 + b = 0
    { 3uv + a = 0
    
    { u^3 - a^3/(27u^3) + b = 0
    { v = -a/(3u)
    
    Обозначим y = u^3
    
    y - a^3/(27*y) + b = 0
    y^2 + b*y - a^3/27 = 0

Решаем квадратное уравнение, делаем обратные замены
и получаем один из корней x.

Конкретно:
    у меня было уравнение
    x^3 - x - 1 = 0

Золотое сечение:
    x^2 - x - 1 = 0
    
Откуда оно взялось?
Мы рассматриваем контейнерный класс из STL С++
   vector<T> -- динамический массив.
В отличие от обычного массива языка C, класс "Динамический массив"
может увеличивать свой размер.

    std::vector<double> a(100);
    a.push_back(123.456);
    assert(a.size() == 101);
    a.push_back(222.);
    assert(a.size() == 102);

    double x = a[12];

    +------------------------+
    |                        |
    +------------------------+
            |  copy
            V                 not used
    +-------------------------------+
    |                         |*****|
    +-------------------------------+
           size               |
    
              capacity              |
              
    new capacity = old capasite * d
    d > 1
    
Размер памяти под массив растет как геометрическая прогрессия
co знаменателем d.

Вопрос: чему равно d?

Загадка: в книге Страуструпа об этом не написано, и повлиять на это
невозможно.

MFC, class CArray -- можно регулировать рост памяти.

В STL, по моим экспериментам:
       В Linux размер массива растет как степень двойки, d = 2
       В MS Windows d = 1.5
       В Java d = 1.125 (???)
       
Какое из чисел d более правильное?

Идея: динамический массив должен иметь возможность использовать
      повторно память, освобожденную на предудих шагах.
      
      Первоначально память размера m
      Первое увеличение:           m*d
      Второе увеличение:           m*d^2       
Память размера m*d пока не освобождена
      Третье увеличение:           m*d^3       
При этом уже освобождена память размера m и размера m*d
память m*d^2 пока не освобождена. Значит, должно выполняться
неравенство
      m*d^3 <= m + m*d
В экстремальном случае выполняется равенство
      m*d^3 = m + m*d, т.е. максимальное d должно удовлетворять
                            уравнению
                            
      d^3 - d - 1 = 0
   
    ПЕРЕРЫВ 10 минут до 11:00

     x^3 - x - 1 = 0
     
     x = u + v
     
     u^3 + v^3 - 1 + 3uv(u + v) - (u + v)  = 0
     -------------   ====================
     
 {   u^3 + v^3 - 1 = 0
 {   3uv - 1 = 0
     
Из второго уравнения v = 1/(3u)
     u^3 + 1/(27v^3) - 1 = 0
Обозначим y = u^3
    y + 1/(27y) - 1 = 0
    y^2 - y + 1/27 = 0
    
    y = ( 1 + sqrt(1 - 4/27) ) / 2
             
       
Конструирование собственных классов 

2 типа:
1) классы, не захватывающие ресурсы (в конструкторе или в процессе
   выполнения методов). Пример -- class Complex
    class Complex {
    private:
        double re;
        double im;
        
2) классы, которые захватываю ресурсы (динамическая память, ...)
   Пример:
    class Matrix (плохая реализация!)
    Под элементы матрицы мы захватываем память в куче
    с помощью оператора new
    
. . .

Внимание!
Ни в коем случае нельзя вычислять определитель матрицы с помощью
разложения по строке, решать систему линейных уравнений с 
помощью правила Крамера, вычислять обратную матрицу с помощью
алгебраических дополнений и т.п. методов из чистой математики.
Все эти методы
    АБСОЛЮТНО НЕПРИГОДНЫ ПРАКТИЧЕСКИ!!!
Например, при вычислении определителя разложение по строке в результате
получаются n! слагаемых, что на практике равно бесконечности.

ВСЕ ЭТИ АЛГОРИТМЫ НАДО РЕАЛИЗОВЫВАТЬ С ПОМОЩЬЮ МЕТОДА ГАУССА!

Обратная матрица к
1 2 3
4 5 6   Алгоритм вычисления:
7 8 0
  
1 2 3 | 1 0 0
4 5 6 | 0 1 0  --> Gauss -->
7 8 0 | 0 0 1
 
* * * | * * *   Если матрица невырождена, то
0 * * | * * *   элементы на диагонали отличны от нуля -->
0 0 * | * * *


* 0 0 | * * *   Обратный проход (снизу вверх):
0 * 0 | * * *   зануляем элементы выше диагонали -->
0 0 * | * * *

Домножаем строки на элементы, обратные к диагональным -->

1 0 0 | * * *
0 1 0 | * * *   Справа будет получена обратная матрица
0 0 1 | * * *


Важно!!!

В стандарте C++ 1985 не было шаблонов (template).
Например, мы определили матрицу с вещественными элементами.
Допустим, я хочу рассмотреть матрицы на полем Zp, мне придется заново
писать класс MatrixZp, который будет почти полностью копировать
уже написанный класс матриц над вещественными числами.

template <class ValueType> class Matrix {
    int m;
    int n;
    ValueType* elems;
public:
    // Всюда слово double заменяется на ValueType
    . . .
};

Это появилось в C++98.

Но!!!
Страуструп -- так не надо писать! Надо использовать
стандартные контейнеры из библиотеки STL языка C++

#include <vector>

template <class ValueType> class Matrix {
    int m;
    int n;
    std::vector<ValueType> elems;
public:
    // Всюда слово double заменяется на ValueType
    . . .
};

========================================================

01 Oct 2021

Типичная ошибка:
подключение *.cpp-файла с помощью

#include "R3Graph.cpp"

Подключаются только заголовочные файлы (h-файлы),
которые не содержать текста, приводящего к созданию кода
программы. Они содержат только информацию, нужную для связи модулей.

Например, в h-файле можно задать константу
const double PI = 3.1415926;
но нельзя определить глобальную переменную:
double eps = 1e-12;
Если нужна глобальная переменная, то в h-файле она описывается
как внешняя:
extern double eps;
а в каком-то cpp-файле она определяется:
double eps = 1e-12;

h-файл содержит только прототипы функций, но не их определения.
Исключение -- inline-функции.

Например,

#include <iostream>

class Complex {
    . . .
};

std::ostream operator<<(std::ostream& s, const Complex& c) {
    s << c.real() << " + " << c.imag() << "*i";
    return s;
}
Ашипка!!! Такое можно писать только в cpp-файле
(в данном случае в файле "Complex.cpp"), либо сделать эту функцию
inline-новой.

inline std::ostream operator<<(std::ostream& s, const Complex& c) {
    s << c.real() << " + " << c.imag() << "*i";
    return s;
}

==================================================================

Конструирование собственных классов

class Complex;

class Matrix;

Давайте рассмотрим более современную реализацию класса
Matrix, использующую стандартный контейнер std::vector<double>
из стандартной библиотеки STL языка C++ для хранения элементов
матрицы.

#include <vector.h>  ???? Почему так не пишем???
#include <vector>    -- подчеркивает, что vector -- не обязательно
                        файл; но это обязательная часть компилятора C++


    ПЕРЕРЫВ 10 минут до 10:50

Методы вычисления ранга, определителя, обратной матрицы и решения
системы линейных уравнений сводятся к алгоритму Гаусса приведения
матрицы к ступенчатому виду.

Напомним алгоритм вычисления обратной матрицы.

   А --> E единичная матрица
         E  = tk ... t3 t2 t1 * A    ==> домножим справа на A^(-1)
         A^(-1) = tk ... t3 t2 t1
         
   A | E
приводим расширенную матрицу так, чтобы слева была единичная матрица      
Тогда справа будет обратная матрица

     1 2 3 | 1 0 0      * * * | * * *
     4 5 6 | 0 1 0  --> 0 * * | * * * --> обратный проход
     7 8 0 | 0 0 1      0 0 * | * * *     снизу-вверх, справа налево
                   Gauss
                                                                                
     * 0 0 | * * *      1 0 0 | * * *                                        
-->  0 * 0 | * * * -->  0 1 0 | * * *
     0 0 * | * * *      0 0 1 | * * *                                           
                                Обратная                                        
                                матрица                                         
                                                                                
===========================================

8 Oct 2021

Идея: вместо самостоятельного захвата памяти в куче
Б.Страуструп предлагает всегда использовать контейнеры
из стандартной библиотеки
    vector, deque, list, set, map, ...                                                                                 
Они берут на себя всю работу по управлению динамической
памятью.

Сегодня:
главная тема -- схема построения цикла с помощью инварианта.

В 60-е -- 70-е годы было поплярно так называемое
"структурное программирование", блок-схемы, доказательства
правильности работы программ и т.п.

Доказательство правильности программ и протоколов
вноь стало актуальным в связи с ростом популярности Internet.

Правила Хоара-Флойда, позволяющие доказывать правильность
работы программы.

    assert: alpha
    if (alpha) {
        block 1;    // alpha => beta
    } else {
        block 2;    // not alpha => gamma
    }
    assert: beta or gamma
    
    
    block1: 
    для него можно доказать, что alpha => beta
    
    block2: 
    для него можно доказать, что not alpha => gamma
    
    Правило инварианта цикла
    
    assert: I(x0)
    x = x0
    assert: I(x)
    while (not Q(x)) {
        assert: I(x)
        x = T(x)    // T сохраняет инариант:
                    // I(x) => для y = T(x) выполняется I(y)
    }
    assert: I(x) and Q(x) 
    
Искусство программиста состоит  умении применить эту абстрактную
схему на практике.

Рассмотрим несколько задач на применение схемы инварианта цикла.

1. Вычисление наиболшего общего делителя двух целых чисел
   Дано: числа m, n <- Z, m != 0 или n != 0
   Надо: вычислить gcd(m, n)  (НОД)
   
Применим схему инварианта цикла
Конфигурационное пространство
    x = (a, b)
Инвариант цикла: gcd(a, b) = gcd(m, n)
Начальная точка: (a, b) = (m, n)
Условие завершения:
    b = 0
    Тогда gcd(a, 0) = a, поскольку 0 делится на любое число
(d делит x, d|x, если x можно представить в виде
                      x = d*q
0 = d*0, т.е. любое число d делить 0 => НОД(a, 0) = наибольший
                                                    делитель a
Цель: сделать вторую компоненту пары (a, b) нулевой,
т.е. надо уменьшать вторую компоненту.
                                                       
    a = q*b + r,    |r| < |b|
    
gcd(a, b) = gcd(b, r), поскольку эти две пары имеют одно и то же
множество общих делителей.
Доказательство:
    пусть d|a, d|b
    r = a - q*b,    d|a, d|q*b =>  d|r
    Обратно, пусть d|b, d|r => d|(q*b + r = a)

Завершение цикла обычно связано с функцией Ляпунова:
некоторой ограниченной снизу величиной, которая монотонно убывает
(или ограниченной сверху, монотонно возрастает).

Второй пример: алгоритм быстрого возведения в степень.

Дано: элемент полугруппы a (число, элемент кольца вычетов по модулю m,
                            матрица, полином, ... всё, что можно умножать
                            и операция умножения ассоциативна)
      n >= 0 -- целое число
Надо: вычислить a^n

Применяем схему инварианта цикла

    Инвариант:  b^k * p == a^n
    Начальные значения:
                b = a, k = n, p = 1
    Преобразование T
                        { (b*b, k/2, p),    k четное
        T(b, k, p) =    {        
                        { (b, k - 1, p*b),  k нечетное
                        
    Условие завершения:  k == 0
                         тогда ответ в переменной p 

Расширенный алгоритм Евклида

Основная теорема арифметики.
    Пусть есть два числа m, n <- Z, хотя бы одно из них отлично от 0.
    Тогда их наибольший общий делитель d = gcd(m, n) 
    выражается в виде их линейной комбинации
        d = u*m + v*n   
    для некоторых целых чисел u, v <- Z
    
Числа u, v можно найти с помощью Расширенного алгоритма Евклида.

Дано: целые числа m, n <- Z,  m != 0 или n != 0
Надо: вычислить наибольший общий делитель d = gcd(m, n)
      и два целых числа u, v такие, что
        d = u*m + v*n

        ПЕРЕРЫВ 10 минут до 10:45
        
Инвариант цикла:
  +-----------------------+  
  | gcd(a, b) = gcd(m, n) |
  | a = u1*m + v1*n       |
  | b = u2*m + v2*n       |
  +-----------------------+

Начальные значения:
    a = m, b = n,
    u1 = 1, v1 = 0,
    u2 = 0, v2 = 1
    
Условие завершения цикла: b = 0
     Тогда ответ: d = a, u = u1, v = v1
      
Преобразование в теле цикла, сохраняющее инвариант и уменьшающее |b|
    a = q*b + r, причем |r| < |b|
    (a, b) <-- (b, r)
    r = a - q*b = u1*m + v1*n - q*(u2*m + v2*n) =
                = (u1 - q*u2)*m + (v1 - q*v2)*n
                
    u1' <- u2           (новое значение u1)
    u2' <- u1 - q*u2    (новое значение u2)
    
    v1' <- v2
    v2' <- v1 - q*v2
    
Зачем нужен расширенный алгоритм Евклида?
Основная цель -- найти обратный к x элемент y в
                 кольце вычетов по модулю m
Дано: x != 0, m != 0     x, m <- Z
надо: найти y:     y*x == 1 (mod m)

(выполнить деление в кольце Z_m).

Кольце Z_m обратимы элементы, взаимно простые с m
    x обратимо в Z_m  <=> gcd(x, m) = 1
    
Найдем обратный элемент. Применив Расширенный алг. Евклида,
                         найдем u, v:
                    1 = u*x + v*m == u*x (mod m)    
                              ===
                              это 0 в Z_m
                    Тогда y = u -- искомый обратный элемент 
                               
Еще один пример на примение схемы построения цикла с помощью инварианта:
вычисление логарифма без разложения в ряд.

Дано:
    a > 1 -- основание логарифма
    x > 0,   x, y <- R
Надо:
    вычислить логарифм от x по основанию a
    y = log_a(x)      
    
Это означает, что    a^y = x
Надо, чтобы это равенство выполнялось приближенно.
Инвариант:
                    a^y * z^t == x    
у должен отличаться от точного значения не больше чем на eps.

Начальные значения:
                    y = 0, z = x, t = 1
Условие завершения:
                    z^t == 1 (примерное равенство)
Для этого достаточно, чтобы
                    1/a <= z <= a   и
                    |t| <= eps
            Тогда z^t будет близко к единице и ответ будет в y
            
Цель преобразования в теле цикла -- уменьшить t.
                    { (y-t, t, z*a),    z < 1/a
        (y, z, t) = { (y+t, t, z/a),    z > a
                    { (y, t/2, z*z),    1/a <= z <= a

=================================================================

15 Oct 2021

Структуры данных (абстрактные типы данных, контейнерные классы
в STL).

Структуры данных -- это классы, которые организуют работу с данными
(добавление, удаление, поиск и чтение данных, ...).
"Абстрактные" -- означает, что мы описываем работу с этим классом
как бы аксиоматически, отвлекаясь от внутренней реализации
класса.

Книга Кушниренко, Лебедев. "Программирование для математиков".

Более-менее стандартные структуры данных:
1) динамический массив  std::vector<E>
2) стек                 std::stack<E> (я его никогда не использую)
3) очередь              std::queue<E> (я его никогда не использую)
4) deque (дек)          std::deque<E> (я его использую
                            как сам по себе, так и вместо
                            классов stack и queue)
5) бинарная куча        std::priority_queue (я его никогда не
                            использую, потому что в нем не
                            хватает одного очень важного метода,
                            который нужен для алгоритма Дейкстры
                            нахождения минимального пути в графе;
                            мне приходится самостоятельно
                            реализовывать бинарную кучу)
6) Л2-список            std::list<E>
6') Л1-список           std::forward_list<E>
7) множество            std::set<E>
8) отображение,         std::map<KeyType, ValueType>
   словарь,
   нагруженное множество
   BinaryHeap
8') многозначное отображение std::multimap<KeyType, ValueType>

В STL предложены более-менее стандартизованные имена методов:
    push    -- добавить
                куда?   front -- начало, back -- конец
                std::vector<double> a;
                a.push_back(1.5);
    pop      -- удалить
                a.pop_back(); // Удалить конец массива
    clear()  -- очистить (удалить все элементы)
    empty()  -- проверить, пуста ли структура данных,
                возвращает bool
    size()   -- число элементов в структуре данных,
                возвращает число типа size_t, который определяется
                либо как unsigned int в 32-разрядной архитектуре,
                либо как unsigned long long в 64-разрядной архитектуре.
                
ОЧЕНЬ ВАЖНО!!!
    В STL поддерживается понятие итераторов.
    Итератор -- это вложенный класс, который как бы реализует
                указатель на элементы структуры данных.
    std::map<std::string, int> words;
    
    // Указатель на начало структуры данных
    std::map<std::string, int>::iterator i = words.begin();
      
    // Указатель на за-конец структуры данных (after-end)
    std::map<std::string, int>::iterator e = words.end();
        Массив
                0 1 2 3 4 5 6 7
                * * * * * * * *
                ^               ^
                |               |
                begin()         end() -- воображаемый элемент за концом
                                         массива
        size() == 8
        maximal index == 7
        max. index - min.index == size() - 1
        
        v.end() - v.begin() == v.size() 
            
Итераторы можно использовать точно так же, как указатели в языке C.
        ++i, --i, n = i - j, *i -- это элемент, на который ссылается
                                   указатель (разыменование, de-reference)
                              i->first, если структура данных хранит
                                        пары std::pair<A, B>
                              i->second
        i == j, i != j, i < j, ...
        i += n, i -= n
        
Кроме типа iterator, есть еще тип const_iterator, это аналог
константного указателя в языке C. Он позволяет только читать элементы,
но не позволяет их портить (записывать или изменять).
В C++ 2011 появились методы cbegin(), cend(), возвращающие
константные указатели (появилось ключевое слово auto).

Пример: есть множество целых чисел. Найти среднее арифметическое
его элементов.

C++-98:
    std::set<int> s;
    . . .
    int sum;
    std::set<int>::const_iterator i = s.begin();
    while (i != s.end()) {
        sum += *i;
        ++i; // Никогда не пишите i++ !!!
    }
    double meanValue = double(sum)/double(s.size());
    
C++-2011:
    std::set<int> s;
    . . .
    int sum;
    auto i = s.cbegin();
    while (i != s.cend()) {
        sum += *i;
        ++i; 
    }
    double meanValue = double(sum)/double(s.size());
    
==============================================

Структура данных Стек
Это как бы очередь наоборот:
представим себе очередь, в которой первым обслуживают того,
кто встал в нее последним.
Дисциплина работа очереди: FIFO -- First In First Out
Дисциплина работа стека:   LIFO -- Last In First Out

Очередь:    -------------------------------------
             A  B  C ........             X Y Z
            -------------------------------------
            front() == A                   back() == Z
            
В очередь можно добавлять элементы только в конец
            q.push_back(x);
а удалять элементы можно только из начала:
            y = q.front();  // Читаем начало очереди
            q.pop_front();  // и затем удаляем его
            
            
Stack:
        |   a   | Top             Элементы добавляются только в вершину стека
        |   b   |                 (stack -- стог сена, стопка бумаг на столе,
        |   c   |                 стопка тарелок в столовой, 
        | . . . |                 железнодорожный тупик, на который можно
        |       |                 составлять вагоны, ...)
        |   y   |
        |   x   | Bottom of stack
        +-------+
        |   A   |
        |   |   | Пружина
        |   |   |
          
Предписания стека (методы класса):
    s.push(x);      // Добавить x в вершину стека
    y = s.top();    // Прочесть элемент в вершине стека (стек не меняется)
    s.pop();        // Удалить вершину стека (глубина стека уменьшается на 1)
    
Давайте самостоятельно напишем класс стек произвольного типа
(template).

Идея реализации: храним элементы стека в динамическом массиве.
Вершина стека находится в конце массива.

Файл Stack.h (cpp-файла нет, все методы inline).

#ifndef STACK_H
#define STACK_H 1

#include <vector>
#include <cassert>
#include <stdexcept>

template <class ValueType> class Stack {
private:
    std::vector<ValueType> elements;
public:
    Stack(size_t initialSize = 64):
        elements()
    {
        elements.reserve(initialSize);
        assert(elements.size() == 0 && elements.capacity() == initialSize);
    }
    
    size_t size() const {
        return elements.size();
    }
    
    bool empty() const {
        return (size() == 0);
    }
    
    void clear() {
        elements.clear();
    }
    
    const ValueType& top() const {
        if (empty())
            throw std::underflow_error("Stack is empty");
        return elements.back();
    }
        
    ValueType& top() {
        if (empty())
            throw std::underflow_error("Stack is empty");
        return elements.back();
    }
    
    void push(const ValueType& x) {
        elements.push_back(x);
    }
    
    void pop() {
        if (empty())
            throw std::underflow_error("Stack is empty");
        elements.pop_back();
    }       
};    
    
#endif    


        ПЕРЕРЫВ 10 минут до 10:45 

Стек и обратная польская запись (RPN Reverse Polish Notation)

RPN названа в честь польского математика (мат. логика)
Яна Лукасиевича. В обычной записи арифм. выражения мы ставим
знак операции между аргументами
    2 + 7
Такая запись требует скобок для указания порядка операций:
    (2 + 7)*(9 - 3) + 5
    
Ян Лукасиевич предложил два способа записи выражений,
в которых скобки не нужны:
    1) прямая польская запись -- знак операции ставится перед
       аргументами:
       +, *, +, 2, 7, -, 9, 3, 5
       прямая польская запись формулы (2 + 7)*(9 - 3) + 5
    2) обратная польская запись -- знак операции ставится после
       аргументов:
       2, 7, +, 9, 3, -, *, 5, +
В математике мы используем прямую запись, когда рассматриваем
отображения
        z = f(x, y)
(левый стиль записи отображений, префиксная запись).
В математике также используется и правая (постфиксная) запись:
        z = (x, y)f
Такая запись удобнее, потому что композиция отображений записывается
справа налево:
    в левой записи
        (fog)(x) = f(g(x))
    сначала выполняется отображение g, затем f;
    в правой записи
        (x)(fog) = ((x)f)g
        
В программировании в объектно-ориентированных языках 
используется именно правая запись!
        x.f()
действие указывается после аргумента.

Для вычисления прямой польской записи применяется рекурсия:
       +, *, +, 2, 7, -, 9, 3, 5
             9        6
          54  
       59 
Для вычисления обратной польской записи применяется
стековый вычислитель:
       2, 7, +, 9, 3, -, *, 5, +
Читаем запись слева направо. Если это число, то добавляем его в стек.
Если это знак операции, то снимаем со стека аргументы операции,
выполняем действие и добавляем его результат обратно в стек.
    +----+ 2|  2 | 7|  7  | + | 9 | 9 | 9 | 3 | 3 | - | 6 | * | 54 | 5 | 5 |
            +----+  |  2  |   +---+   | 9 |   | 9 |   | 9 |   +----+   | 54|
                    +-----+           +---+   | 9 |   +---+            +---+
                                              +---+
 +  | 59 |                                                  
    +----+
    
RPN и стековый вычислитель используются в программировании
чрезвычайно часто!
Примеры:
   bytecode языка Java, 
   CIL от Microsoft, который используется в .Net (C#, Visual Basic)
        (Common Intermediate Language) CLR C++    _gc
   bytecode языка Python -- всё это есть обратная.
   
Выполнение байткода -- либо интерпретация, либо компиляция "на лету"
(JIT -- Just In Time compiling)

Самый красивый пример -- это язык PostScript.

===================================================

22 Oct 2021

Проект "Стековый калькулятор"

Фирма Hewlett Packard (HP) в конце 60-х выпустила
стековый калькулятор, у которого память организована в виде стека.
При выполнении бинарных операций (сложение, вычитание, умножение, ...)
аргументы операции снимаются со стека, выполняется операция,
результат заносится обратно в стек. Результат последней операции
всегда лежит на вершине стека.
    (11 - 7)*3 + 4/5
    11, 7, -, 3, *, 4, 5, /, +
    
Вычисление на стековом калькуляторе

1)
2) 11
3) 11, 7
4) 4
5) 4, 3
6) 12
7) 12, 4
8) 12, 4, 5
9) 12, 0.8
10) 12.8

Проект Стековый калькулятор
Директория StackCalс
Файлы Stack.h, stcalc.cpp, Makefile

        ПЕРЕРЫВ 10 минут до 11:00
        
Задача на тему "Стековый калькулятор"

Была рассмотрена схема построения цикла с помощью инварианта.
Примеры:
    1) алгоритм Евклида d = gcd(m, n);
    2) расширенный алгоритм Евклида.
       Дано: целые числа m, n
       Надо: найти d = gcd(m, n) и коэффициенты u, v представления
             d в виде линейной комбинации чисел m, n
                    d = u*m + v*n
    3) Алгоритм быстрого возведения в степень в кольце вычетов 
       по модулю m
            p = powmod(a, n, m);   // p == a^n (mod m)
    4) Вычислить обратный к x элементв кольце вычетов 
       по модулю m
            b = invmod(a, m)
       (предполагается, что элементы a, m взаимно протые gcd(a, m) = 1)
    5) китайская теорема об остатках.
       Дано: взаимно простые целые числа
                m1, m2, m3, ..., m_k
             любые целые числа
                r1, r2, r3, ..., r_k
       Надо: найти целое число x такое, что
                x == r_i (mod m_i)
       На стеке должно быть 2*k + 1 число:
            r1, m1, r2, m2, ..., r_k, m_k, k
            
Программа должна снять со стека аргументы, вычислить
требуемое значение (или несколько значений), напечатать красиво
результаты и поместить результаты обратно на стек.

Новая тема:
алгоритмы

#include <algorithm>

Алгоритмы сортировки и поиска

Очень важны, поскольку любой компьютер значительную часть
своего времени тратит на выполнение этих алгоритмов.

3-й том книги Д.Кнут "Искусство программирования"
постящен именно этим двум задачам:
Сортировка и поиск.

Задача сортировки.
Дан массив a[0], a[1], a[2], ..., a[n-1] длины n.

Элементы массива можно сравнивать между собой:
            a[i] < a[j]?

(Замечание: есть 6 команд сравнения:
    ==, !=, <, <=, >, >=
но все их можно реализовать с помощью только одной операции <
Например, два элемента равны x == y тогда и только тогда, когда
        !(x < y) && !(y < x)
Алгоритмы стандартной библиотеки STL используют только это сравнение
a < b, так что, если вы реализуете какой-то класс, для которого
хотите выполнять алгоритмы стандартной библиотеки, то надо для
этого класса реализовать хотя бы один оператор <
class A {
    . . .
public:
    . . .
    bool operator<(const A& y) const; // class method
    
};

или 

bool operator<(const A& x, const A& y); // global function

)

В задаче сортировки требуется упорядочить элементы массива по
возрастанию
    a[0] <= a[1] <= a[2] <= ... <= a[n-1]

Теорема:
для любого алгоритма сортировки найдется вход длины n,
на котором он выполнит не меньше чем
        t = log_2(n!)
сравнений.

То есть невозможно придумать алгоритм, который работает быстрее,
чем за log_2(n!) шагов. (Для оценки времени мы используем число сравнений.)

Пример:
есть 4 монеты и весы, на чашки которых можно класть только по одной монете.
Какое минимальное число взвешиваний достаточно, чтобы упорядочить монеты
по весу?

Ответ: 5 взвешиваний.

Теорема: t >= log_2(4!) = log_2(24) > 4
    => t >= 5
    
Можно предъявить конкретный алгоритм, который на любом входе
выполняет не больше 5 взвешиваний. Конкретно, это алгоритм сортировки
слиянием (Merge Sort)
                     5 монет
                   2     +     2
        1 взвешивание    +   1 взвешивание   
                     слияние:
                2 или 3 взвешивания.

Оценим ассимптотику функции log_2(n!) при n-->inf
Воспользуемся формулой Стирлинга:
            n! ~ n^n e^(-n) sqrt(2*pi*n)
        log_2(n!) = n*log_2(n) + o(n*log_2 n) ~ n*log_2(n)
        
Оптимальный алгоритм сортировки: работающий за время O(n*log(n)).

Докажем эту теорему.


Если в бинарном дереве не меньше n! вершин, то высота дерева
не меньше чем 
                h >= log_2(n!) + 1

Максимальное количество матчей в турнире и есть высота минус 1
=>  число сравнений >= log_2(n!)

Для полного бинарного дерева число терминальных вершин k = 2^(h-1)
Следовательно, h - 1 = log_2(k)

Для произвольного бинарного дерева справедливо неравенство

               h - 1 >= log_2(k)
доказательство: достроим дерево до полного той же высоты,
число терминальных вершин при этом может только увеличиться 
                h - 1 = log_2(k'),   k' >= k =>
                h - 1 = log_2(k') >= log_2(k)

h - 1 -- это время работы алгоритма сортировки для самого плохого входа
                t >= log_2(n!), так как k >= n!

===========================================================

29 Oct 2021

Стандартная библиотека STL C++: итераторы,
реализация цикла "для каждого элемента структуры данных"

В языке C есть две конструкции цикла

    while (условие) {
        . . .
        тело цикла
        . . .
    }
    
Цикл for -- калька с арифметических циклов языков
Algol, Fortran
    for (int i = 0; i < n; ++i) {
        . . .
        тело цикла, зависящее от переменной i
        . . .
    }
        
        
Цикл for, вообще говоря, не нужен в C, поскольку
он целиком реализуется с помощью цикла while:

    int i = 0;
    while (i < n) {
        . . .
        тело цикла, зависящее от переменной i
        . . .
        ++i;
    }
        
Вообще говоря, в продвинутом языке программирования
обязательно нужен цикл "для каждого элемента структуры
данных"

    Set s;
    . . .
    for x in s {
        . . .
        тело цикла, зависящее от переменной x
        . . .
    }
    
Например, такие циклы реализованы в языке Python

    p = [2, 3, 5, 7, 11, 13, 17, 19]
    m = 1
    for x in p:
        m *= x
        
В Питоне имеется конструкция iterator, 
с помощью которой красиво
реализован цикл for.

Мехмат, язык Pra (Г.Лебедев, А.Кушниренко)
Кушниренко, Лебедев
Программирование для математиков

В стандартной библиотеке имеется общий подход
к реализии цикла "для каждого" и вообще алгоритмов,
которые работают с элементами структуры данных,
позволяя ничего не знать о ее внутреннем устройстве.

Пусть есть структура данных (контейнерный класс, абстрактный
тип данных) S. Тогда внутри класса S описаны классы
    class S::iterator
    class S::const_iterator
которые как бы моделируют указатели на элементы структуры данных S.
Эти "указатели" позволяют получить доступ к элементам класса.

Пример: суммирование массива в языке C

Традиционная запись:
    double a[100];
    . . .
    double s = 0.;
    int i = 0;
    while (i < 100) {
        s += a[i];
        ++i;
    }

Язык С отличался от всех существующих к концу 60-х годов языков
программирования наличием указателей. В стиле времен начала языка C
можно написать этот фрагмент, используя только указатели:
    double a[100];
    . . .
    double s = 0.;
    const double *i = a;
    const double *e = a + 100;  // Указатель на за-конец массива a
    while (i != e) {
        s += *i;
        ++i;
    }
Раньше казалось, что так писать -- это пижонство.
Но, оказалось, что такой стиль хорошо подходит для работы с контейнерными
классами в стандартной библиотеке.
    std::vector<double> a(100);
    . . .
    double s = 0.;
    std::vector<double>::const_iterator i = a.cbegin();
    std::vector<double>::const_iterator e = a.cend();
    while (i != e) {
        s += *i;
        ++i;
    }
    
Замечание:
Естественно пожаловаться на громоздкую запись типа итератора:
    std::vector<double>::const_iterator i = a.cbegin();
В современных стандартах языка C++ появилось слово auto
(на самом деле, это ключевое слово языка C, которое было в нем
с начала 70-х годов, оно означало пожелание для компилятора,
чтобы он расположил переменную не в памяти, а в регистре процессора;
но никто эти не пользовался; с С++ 2003 г. значение этого слова
изменилось, теперь оно означает, что мы предоставляем компилятору
самому определить тип выражения или переменной).

Это синтаксический сахар! Не меняется основной принцип:
язык C++ -- это язык со статической типизацией (в отличие,
например, от языка Python).

    std::vector<double> a(100);
    . . .
    double s = 0.;
    // std::vector<double>::const_iterator i = a.cbegin();
    auto i = a.cbegin();
    // std::vector<double>::const_iterator e = a.cend();
    auto e = a.cend();
    while (i != e) {
        s += *i;
        ++i;
    }

В С++ версии 1998 были только методы begin() и end(),
которые возвращали объект типа iterator,
хотя тип const_iterator существовал уже тогда.
В C++-98 нужно было писать так:
    std::vector<double> a(100);
    . . .
    double s = 0.;
    std::vector<double>::const_iterator i = a.begin();  // не cbegin()
    std::vector<double>::const_iterator e = a.end();    // не cend()
    while (i != e) {
        s += *i;
        ++i;
    }
    
Для типа std::vector итераторы как бы не очень нужны.
Но, например, для типов std::set и std::map без итераторов
вообще невозможно обойтись! Неизвестно, как внутри устроены эти
классы (на самом деле, известно: это красно-черные деревья),
вернее, программисту не нужно (и даже запрещено) полагаться
на конкретное внутреннее устройство этих классов.

Пример:
мы хотим в тексте найти самое часто встречающееся слово.
Для этого мы используем структуру данных отображение,
которое слово отображает на целое число, равное числу вхождений
слова в текст. Пусть мы заполнили эту структуру данных.
теперь найдем самое частое слово.

    std::map<std::string, int> words;
    . . .
    // Заполняем структуру words, сканируя текст, выделяя слова
    // и либо добавляя новое слово со значением 1, либо увеличивая
    // число вхождений в текст слова, уже содержащегося в words
    . . .
    std::string mostFrequentWord;
    int n = 0;
    auto i = words.cbegin();
    while (i != words.cend()) {
        // std::pair<std::string, int>   (p.first, p.second)
        if (i->second > n) {
            mostFrequentWord = i->first;
            n = i->second;
        }
        ++i;
    }
    
===========================================================

Большинство алгоритмов в стандартной библиотеке работают
с итераторами. Например, сортировка массива:

#include <algorithm>


    size_t n; 
    . . .
    std::vector<double> a(n);
    . . .
    std::sort(a.begin(), a.end());  // Сортировка всего массива
    . . .
    int i, j;
    . . .    
    std::sort(a.begin() + i, a.begin() + j);  // Сортировка отрезка массива
                                              // с i-го до j-1-го элемента
                                              
=====================================================================

Если мы самостоятельно реализуем некоторый класс,
хранящий данные, то нам нужно написать два вложенных класса:
iterator и const_iterator

template <class ValueType> class S {
    . . .
    class iterator {
        . . .
        bool operator==(const iterator&);
        bool operator!=(const iterator& i) {
            return !operator==(i);
        }
        iterator& operator++();     // Prefix increment operator ++i
        iterator operator++(int) {  // STUPID postfix operator i++
            iterator oldValue = *this;
            ++(*this);
            return oldValue;
        }
        ValueType& operator*() const;   // x = *i;
        ValueType* operator->() const;  // x0 = i->first; x1 = i->second
        . . .
     };

     class const_iterator {
        . . .   
        const ValueType& operator*() const;   // x = *i;
        const ValueType* operator->() const;  // x0 = i->first; x1 = i->second
        . . .
     };
     
     iterator begin();
     iterator end();
     
     const_iterator cbegin() const;
     const_iterator cend() const;
     . . .
};

    ПЕРЕРЫВ до 10:45
    
Вернемся к нашей теме: алгоритмы сортировки

    Дан массив размера n
        a[0], a[1], ..., a[n-1]
    Надо упорядочить его элементы по возрастанию:
        a[0] <= a[1] <= ... <= a[n-1]
        
Было доказано, что для любого алгоритма сортировки,
основанного на сравнении элементов, найдется вход, на котором он
выполнит не меньше log_2(n!) сравнений.

    log_2(n!) ~ n*log_2(n) -- почти что линейная функция,
                              поскольку функция log растет очень медленно
                              log_2(1000000) < 20
                              1000000 = 10^6 ~= (2^10)^2 = 2^20
                              10^3 < 1024 = 2^10
                              10^3 ~= 2^10
                                KB -- в тысячах, KiB -- в 1024

Время работы любого алгоритма >= n*log_2(n)
Назовем алгоритм сортировки оптимальным, если он работает за время
        O(n*log(n)) 
        
Наивные алгоритмы сортировки работают за время O(n^2):
1) пузырьковая сортировка --
   просматриваем массив от начала к концу, меняем местами
   соседние элементы, если они образуют инверсию
        a[i] > a[i+1]  =>  std::swap(a[i], a[i+1])
2) сортировка прямым выбором:
   просматриваем массив от начала к концу, находя минимальный
   элемент; ставим минимальный элемент в начало массива (выполняя
   функцию swap), повторяем то же самое для отрезка массива,
   начинающегося со второго элемента, и т.д.
   
Существует несколько оптимальных алгоритмов сортировки,
работающих за время O(n*lon(n)):
1) сортировка кучей, или пирамидальная сортировка
   HeapSort;
2) сортировка слиянием MergeSort
   (есть несколько схем такой сортировки, самая современная называется
   TimSort);
3) Radix sort (поразрядная сортировка), 
    используется лишь для сортировки ключей вида
    [s0, s1, s2, ..., s_k-1 ] фиксированной длины k, где s_i
    принадлежат дискретному конечному множеству. Например,
    идентификаторы банковских счетов в Росии -- это 20-значные
    последовательности цифр, например
        23109879871235478245
    Работае за время O(n*k)
4) "Быстрая" сортировка QuickSort, она в среднем быстрая
   (за время O(n*log(n))), но для любой стратегии выбора стержневого
   элемента существует вход, на котором она работает за время O(n^2).
   Идея: выбираем "стержневой" элемент x, меняем местами элементы
   массива следующим образом:
    элементы < x слева,
    элемента, равные x в середине,
    элементы > x справа.
        
    ****** ... ****** xxxxxxxxx *******...****
       * < x          ^         ^  x < *
                      |         |
                      idx0      idx1
       
    Раскидывание элементов налево и направо выполняется за
    время O(n) с помощью функции 
        partition(a, idx_pivot, &idx_pivot0, &idx_pivot1)   
       
    Затем применяем рекурсивно этот же алгоритм
    к левой части массива и к правой части   
    При удачном выборе стержневого элемента левая и правая части
    будут иметь размер примерно n/2 => общее число шагов при удачном
    раскладе не превышает log_2(n). При неудачном выборе стержневого
    элемента число шагов будет порядка O(n) и общее число операций n^2.
    На самом деле, есть алгоритм нахождения медианы за время O(n)
    (возможно, он требует дополнительной памяти). Достаточно найти
    даже не медиану, а элемент x:
        число элементов <= x не меньше чем n*c
        число элементов >= x не меньше чем n*c
        где 0 < c < 1 и c не зависит от входного массива.
    Медиана -- это число x такое, что
        число элементов массива, <= x, равно
        числу элементов массива, >= x
            
Важное свойство алгоритмов сортировки:
алгоритм называется стабильным (или монотонным),
если он не меняет взаимного расположения равных элементов.

Какие из перечисленных алгоритмов стабильные?
1) пузырьковая сортировки ДА
2) сортировка прямым выбором НЕТ
3) сортировка кучей НЕТ
4) сортировка слиянием ДА
5) RADIX-sort ДА
6) быстрая сортировка НЕТ.           
    
Отсюда заключаем, что среди алгоритмов сортировки общего вида
(которые основаны только на сравнении элементов)
сортировка слиянием самая лучшая (и, конкретно, TimSort).

Но мы сначала рассмотрим сортировку кучей, потому что понятие кучи
очень важно само по себе. Например, бинарная куча (или, как ее иногда
называют, очередь с приоритетами priority_queue) используется в 
ЗАМЕЧАТЕЛЬНОМ алгоритме Дейкстры нахождения минимального пути в графе
(его мы обязательно рассмотрим).

==========================================

12 Nov 2021

На прошлых лекциях была начата тема "Сортировка"

Была доказана теорема:
для любого алгоритма сортировки существует вход,
на котором он выполнит не меньше чем
    log_2(n!)
сравнений. Это дает оценку снизу для скорости работы
любого алгоритма сортировки общего вида, когда нам
ничего не известно о строении ключей, мы можем
только сравнивать ключи 
        x < y: bool

-------------------------------------

void sort(firstIterator, lastIterator);
void sort(firstIterator, lastIterator, compare);

Итераторы должны обеспечивать прямой доступ (random access)
к элементам сортируемой последовательности: 
    x = firstIterator[i];
    y = *(firstIterator + i);
    
Как упорядочить динамический массив (или deque)?
#include <algorithm>

    std::vector<double> a(1000);
    . . .
    std::sort(a.begin(), a.end());
    
Если хотим упорядочить массив по убыванию элементов, то
используем собственную функцию сравнения:

bool greaterThan(const double& x, const double& y) {
    return (x > y);
}

Сравнивать элементы можно не только с помощью функции,
но и с помощью функционального объекта. Это используется,
когда нам нужны какие-нибудь дополнительные параметры для
сравнения.

Пусть, например, мы хотим сравнивать строки произвольной длины
только по первым n символам. (Строки -- это последовательности
ASCII-символов, терминированные нулевым байтом, как в языке C.
Это на самом деле атавизм, потому что сейчас тексты представляются
символами Unicode.)

    bool compareStrings(const char* x, const char* y, size_t n); 

class StringCompare {
public:
    size_t n;
    
public:
    StringCompare(size_t nn = 0):
        n(nn)
    {}
    
    bool operator()(const char* x, const char* y) const {
        int c;
        if (n == 0) {
            c = strcmp(x, y);
        } else {
            c = strncmp(x, y);
        }
        return (c < 0);
    }
};
            
    std::vector<const char *> text(1000);
    . . .
    
    StringCompare cmp(6);
    // bool b = cmp(x, y);

    std::sort(text.begin(), text.end(), cmp);        
        
    
Для реализации всех шести операция сравнения достаточно уметь
выполнять только операцию <

    a == b  <==>   !(a < b || b < a)
    a <= b  <==>   a < b || a == b
    a > b   <==>   b < a
    a >= b  <==>   b < a || a == b
    
-------------------------------------------------------------    
    
Была доказана теорема:
для любого алгоритма сортировки существует вход,
на котором он выполнит не меньше чем
    log_2(n!)
сравнений.
Эта функция при n -> infinity эквивалентна
    log_2(n!) ~ n*log_2(n)
(доказывается с помощью формулы Стирлига
    n! ~ n^n exp(-n) sqrt(2*pi*n)

Считаем алгоритм сортировки оптимальным, если он работает за время
    O(n*log_2(n))
    
Оптимальные алгоритмы сортировки:
1) пирамидальная сортировка, или сортировка кучей
   HeapSort
2) сортировка слиянием
   MergeSort
-----------------------------------------
3) "быстрая" сортировка QuickSort -- в среднем работает
    за O(n*log(n)), но при любой (несложной) стратегии
    выбора стержневого существует вход, на котором этот алгоритм
    работает за время O(n^2). В результате от этого дурацкого
    алгоритма все постепенно отказываются. (Например, Java, Python
    используют алгоритм TimSort, который есть вариант сортировки
    слиянием)
------------------------------------------
4) RADIX Sort -- поразрядная сортировка. В ней ключи имеют специальный
   вид: последовательности фиксированной длины из элементов
   конечного дискретного множества, например,
   строки фиксированной (или ограниченной) длины из десятичных цифр,
   из ASCII-символов
   (не из символов Unicode!!! Число символов в кодировке Unicode
   больше 1000000, это уже не подходит для эффективной реализации
   RADIX Sort).
   Скорость работы O(k*n), где k -- длина ключей, n -- размер массива.

Алгоритм сортировки называется стабильным (stable), если он не меняет
взаимного порядка равных элементов. Это ОЧЕНЬ важное свойство
алгоритма сортировки. Какие алгоритмы сортировки стабильны?

1. Пузырьковая сортировка t = O(n^2)        стабильна
2. HeapSort --                              не стабильна
3. MergeSort --                             стабильна
4. QuickSort                                не стабильна
5. RADIX Sort                               стабильна

Необходимость в дополнительной памяти
(может ли алгоритм сортировки работать только за счет сравнений
и обменов std::swap(x, y) элементов массива?)

1. Пузырьковая сортировка t = O(n^2)    не требует дополнительной памяти
2. HeapSort --                          не требует дополнительной памяти
3. MergeSort --                         ТРЕБУЕТ дополнительной памяти
4. QuickSort                            не требует дополнительной памяти
5. RADIX Sort                           ТРЕБУЕТ дополнительной памяти

Возможность параллельного выполнения
лучший с этой точки зрения алгоритм -- MergeSort

HeapSort использует понятие бинарной кучи (очередь с приоритетами,
std::priority_queue), которое само по себе ценно.

        ПЕРЕРЫВ 10 минут до 10:55
    
Алгоритм сортировки слиянием.
Основная идея:
    если есть два уже отсортированных массива,
    то их можно объединить в один отсортированный массив
    за время O(n + m), где n, m -- размеры исходных массивов.
    
merge(a, n, b, m, c);
        a: x x x x x x x * * * * * * * * * 
                         ^   
                         |

        b: y y y y y * * * * * * 
                     ^   
                     |
    
    
c: x y x x y x x y y x x y  
                           ^
                           |


template <class RandomAccessIterator>
void merge(
    RandomAccessIterator firstA, RandomAccessIterator lastA,
    RandomAccessIterator firstB, RandomAccessIterator lastB,
    RandomAccessIterator firstC
) {
    while (firstA < lastA && firstB < lastB) {
        if (*firstB < firstA) {
            *firstC = *firstB; ++firstB;
        } else {
            // assert(*firstA <= *firstB);
            *firstC = *firstA; ++firstA;
        }
        ++firstC;
    }
            
    while (firstA < lastA) {
        *firstC = *firstA;
        ++firstA; ++firstC;
    }
        
    while (firstB < lastB) {
        *firstC = *firstB;
        ++firstB; ++firstC;
    }
}        
        
        
Следующий раз: обработка текстов в Unicode (формат UTF-8)
с помощью классов и функций стандартной библиотеки C++.

После будет программирование в QT.

=========================================