Классические алгоритмы на графах

Обходы графов

Обходом графа называется последовательное посещение всех вершин некоторой его компоненты связности (возможно с каким-то действием в каждой вершине).

В этой главе описаны два наиболее распространённых вида обходов: в глубину и в ширину. Примеры кода приведены на трёх разных языках: C++, Javascript и Haskell.

Моделирование графа

Один из самых простых (и при этом весьма эффективных) способов моделирования графов — таблица смежности: хранение для каждой вершины графа множества её соседей или же — если рёбра графа содержат какую-то дополнительную информацию (называемую весом ребра) — таблицы, сопоставляющей каждому соседу информацию на соответствующем ребре.

В этой главе нас будут интересовать только графы без весов.

Сперва приведём пример вышеописанной модели графа на языке C++. Для простоты считаем, что вершины графа — все целые числа.

#include <map>
#include <set>

using Graph = std::map<int,std::set<int>>;

// добавление ориентированного ребра
void addEdge(Graph &g, int i, int j) {
    g[i].insert(j);
}

// добавление неориентированного ребра
void addEdge2(Graph &g, int i, int j) {
    addEdge(g,i,j);
    addEdge(g,j,i);
}

// множество соседей указанной вершины
std::set<int> neighbors(Graph &g, int v) {
    return g[v];
}

Отметим, что для языка C++ характерны мутирующие операции: при изменениях объект сохраняет свою идентичность, но изменяет своё состояние. Типичный код, создающий нетривиальный граф, выглядит так:

Graph testGraph() {
    Graph g;
    addEdge2(g,1,2);
    addEdge2(g,1,3);
    addEdge2(g,2,3);
    addEdge2(g,4,5);
    addEdge2(g,5,6);
    addEdge2(g,4,7);
    addEdge2(g,7,8);
    return g;
}

В функциональных языках предпочтение отдаётся иммутабельным (т.е. неизменяемым) структурам данных: изменение моделируется созданием новой структуры, для эффективности ссылающейся на части старой.

Например, библиотека Immutable для Javascript содержит иммутабельные разновидности наиболее распространённых структур данных. В терминах структур этой библиотеки наша модель графа может выглядеть так:

const im = require("immutable") // это для nodejs
const im = Immutable            // это для браузера
// из вышенаписанных двух строчек нужно оставить одну

function emptyGraph() {
    return im.Map()
}

function addEdge(g,i,j) {
    let neighbors = g.get(i)

    if (neighbors) {
        return g.set(i,neighbors.add(j))
    }

    return g.set(i,im.Set([j]))
}

function addEdge2(g,i,j) {
    return addEdge(addEdge(g,i,j),j,i)
}

function neighbors(g,i) {
    let neighbors = g.get(i)
    
    if (neighbors) { return neighbors }
    
    return im.Set()
}

Немного большая длина кода (по сравнению с вариантом на C++) связана с отсутствием возможности задать значение по-умолчанию для тех вершин, которые ещё не хранятся в Map. Для иммутабельных структур данных характерен следующий способ построения:

function testGraph() {
    let g = emptyGraph()
    g = addEdge2(g,1,2)
    g = addEdge2(g,1,3)
    g = addEdge2(g,2,3)
    g = addEdge2(g,4,5)
    g = addEdge2(g,5,6)
    g = addEdge2(g,4,7)
    g = addEdge2(g,7,8)
    return g;
}

Ну и в завершение приведём ту же самую модель на языке Haskell.

import Data.Map(Map)
import qualified Data.Map as Map
import Data.Set(Set)
import qualified Data.Set as Set


type Graph = Map Int (Set Int)

emptyGraph :: Graph
emptyGraph = Map.fromList []

addEdge :: Int -> Int -> Graph -> Graph
addEdge i j g = Map.alter (Just . update) i g
  where 
    update Nothing  = Set.fromList [j]
    update (Just s) = Set.insert j s

addEdge2 :: Int -> Int -> Graph -> Graph
addEdge2 i j g = addEdge i j $ addEdge j i g

neighbors :: Int -> Graph -> Set Int
neighbors i g = case Map.lookup i g of
    Nothing  -> Set.empty
    (Just s) -> s

Нетривиальный граф в Haskell строится примерно так же, как и средставми Immutable в Javascript:

-- не забыть в верху программы сделать import Data.Foldable

testGraph :: Graph
testGraph = 
    foldl' go emptyGraph [(1,2), (1,3), (2,3), (4,5), (5,6), (4,7), (7,8)]
  where
    go g (i,j) = addEdge2 i j g

Наивный обход в глубину

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

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

Реализация на C++.

#include <vector>

std::set<int> component(Graph &g, int v) {
    std::vector<int> path{v};
    std::set<int> visited{};
    
    while (path.size() > 0) {
        int current_vertex = path.back();
        visited.insert(current_vertex);
        
        bool found_new = false;
        
        auto n = neighbors(g,current_vertex);
        for (int i : n) {
            if (visited.count(i) == 0) { 
                path.push_back(i);
                found_new = true; 
                break;
            }
        }
        
        if (found_new) { continue; }
        
        path.pop_back();
    }
    
    return visited;
}

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

component :: Int -> Graph -> Set Int
component i g = go [i] Set.empty
  where
    go [] s = s
    go (v:vs) s =
      let s' = Set.insert v s
          n  = neighbors v g
          d  = Set.difference n s'
      in case Set.toList d of
          [] -> go vs s'
          (next:_) -> go (next:v:vs) s'

Рекурсивный обход в глубину

Путь, хранимый при наивной реализации, можно моделировать цепочкой вызовов рекурсивной функции. Единственная проблема — подобрать адекватное рекуррентное соотношение.

Пусть S <g> v — объединение множества S и компоненты связности вершины v в графе g без S. Для этой операции верно следующее (попытайтесь понять, почему):

1. Если S содержит v, то S <g> v = S
2. В противном случае S <g> v = (S|{v}) <g> v_1 <g> v_2 <g> ... <g> v_n, где v_i -- соседи вершины v в графе g

Приведём сперва реализацию на Javascript:

function component(g,v) {
    let op = (s,x) => {
        if (s.has(x)) { return s }
        
        s = s.add(x)
        
        for (let n of neighbors(g,x)) { s = op(s,n) }
    }
    
    return op(im.Set(), v)
}

Ну и, как всегда, на Haskell (тут она не особо короче, но зато является дословной трансляцией определения операции <g>):

component :: Int -> Graph -> Set Int
component i g = Set.empty <> i 
  where 
    s <> v
      | v `Set.member` s = s 
      | otherwise = foldl' (<>) (Set.insert v s) (neighbors v g)

Фронтовой обход в глубину

Что интересно, тот же обход в глубину можно получить, если моделировать не путь от начала обхода до текущего местоположения, а границу между теми вершинами, где алгоритм уже побывал, и теми, где ему ещё предстоит побывать.

Вернёмся к C++, чтобы показать, насколько более простым при таком подходе становится код:

#include <vector>

std::set<int> component(Graph &g, int v) {
    std::vector<int> border{v};
    std::set<int> visited{};
    
    while (border.size() > 0) {
        int current_vertex = border.back(); border.pop_back();
        
        if (visited.count(current_vertex)) { continue; }
        
        visited.insert(current_vertex);
        
        auto n = neighbors(g,current_vertex);
        for (int i : n) { 
            if (visited.count(i) == 0) { border.push_back(i); }
        }
    }
    
    return visited;
}

Ну и, конечно же, вариант на Haskell:

component :: Int -> Graph -> Set Int
component i g = go [i] Set.empty
  where
    go [] s = s
    go (v:vs) s
        | v `Set.member` s = go vs s
        | otherwise = 
          let s' = Set.insert v s
              n  = neighbors v g
              d  = Set.difference n s'
          in go (Set.toList d ++ vs) s'

Обход в ширину

Интересным свойством фронтового обхода в глубину является лёгкость его обобщения. Сам алгоритм в общем виде весьма прост: пока фронт не пуст

1. взять из него вершину
2. если она уже посещена, вернуться на предыдущий пункт
3. пометить взятую вершину как посещённую
4. добавить к фронту всех непосещённых соседей взятой вершины

Опционально можно убрать второй пункт или же проверку посещённости в четвёртом пункте (но гораздо лучше оставить и то, и другое).

В зависимости от того, как именно работает операция «взять из фронта вершину», можно получить как алгоритмы обхода «в глубину» и «в ширину», так и кучу промежуточных разновидностей обходов.

Алгоритм обхода в ширину выглядит так (на C++):

#include <deque>

std::set<int> component(Graph &g, int v) {
    std::deque<int> border{v};
    std::set<int> visited{};
    
    while (border.size() > 0) {
        int current_vertex = border.front(); border.pop_front();
        
        if (visited.count(current_vertex)) { continue; }
        
        visited.insert(current_vertex);
        
        auto n = neighbors(g,current_vertex);
        for (int i : n) { 
            if (visited.count(i) == 0) { border.push_back(i); }
        }
    }
    
    return visited;
}

Единственное отличие: фронт заполняется с одного конца, а опустошается — с другого. А вот картина обхода меняется значительно: теперь алгоритм посещает вершины в порядке их удалённости от начальной (сначала те, до которых можно добраться за один шаг, потом те, до которых можно добраться за два шага и так далее).

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

Кратчайший путь

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

Нами будет рассмотрен единственный алгоритм — алгоритм Дейкстры. Он представляет собой обычный обход в ширину с двумя модификациями:

• вместо множества посещённых вершин — таблица расстояний от начальной
• вместо граничной очереди — очередь с приоритетом

Напомним, что очередь с приоритетом отличается от очереди/стека тем, что:

• каждый элемент сопровождается приоритетом — элементом некоторого линейно упорядоченного множества
• операция изъятия элемента изымает элемент с наименьшим (или же наибольшим) приоритетом

Очередь с приоритетом

Эффективную очередь с приоритетом можно получить из упорядоченного мультисловаря (это когда допускается хранить множество элементов по одному ключу) или же из обычного словаря, преобразовав его в мультисловарь одним из двух способов:

• мультисловарь типа k -> v — это словарь типа k -> List v, где List v — тип списков элементов типа v
• мультисловарь типа k -> v — это словарь типа (k, int) -> v

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

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

Без лишних слов Дейкстра

Возьмём стандартный обход в ширину (из прошлой главы):

std::set<int> component(Graph &g, int v) {
    std::deque<int> border{v};
    std::set<int> visited{};
    
    while (border.size() > 0) {
        int current_vertex = border.front(); border.pop_front();
        
        if (visited.count(current_vertex)) { continue; }
        
        visited.insert(current_vertex);
        
        auto n = neighbors(g,current_vertex);
        for (int i : n) { 
            if (visited.count(i) == 0) { border.push_back(i); }
        }
    }
    
    return visited;
}

и немного модифицируем:

// всё делается в предположении, что Graph означает следующее:
using Graph = std::map<int, std::map<int, double>>;

std::map<int, double> distances(Graph &f, int v) {
    std::map<std::pair<double,int>, int> border;
    border[{0,0}] = v;
    int unique = 1;
    
    std::map<int, double> visited; // таблица расстояний
    
    while (border.size() > 0) {
        // it -- "указатель" на пару вида ((расстояние,мусор),вершина)
        auto it = border.begin();
        double d = (it->first).first;
        int current_vertex = it->second;
        border.erase(it);
        
        if (visited.count(current_vertex)) { continue; }
        
        visited[current_vertex] = d;
        
        for (auto n : g[current_vertex]) {
            int neighbor = n.first;
            double edge_length = n.second;
            
            if (visited.count(neighbor) == 0) { 
                double new_d = d + edge_length;
                border[{new_d, unique++}] = neighbor; 
            }
        }
    }
    
    return visited;
}

На каждом шагу алгоритм выбирает вершину, ближайшую к обработанному множеству. Ближайшую — по путям, все рёбра которых, кроме последнего, соединяют вершины обработанного множества.

Нетрудно заметить, что среди кратчайших путей до выбранной вершины есть путь, который имеет именно такой вид: если бы в нём было хотя бы два ребра, выходящих за пределы обработанного множества, он был бы не короче того, в котором одно такое ребро, причём кратчайшее из возможных.

Ну и для полноты картины приведём код на Haskell (без импортов и тех функций, которые не относятся к сути дела):

type Vertex = Int
type Distance = Double
type Graph = Map Vertex (Map Vertex Distance)

neighbors :: Vertex -> Graph -> [(Vertex, Distance)]
neighbors i g = case Map.lookup i g of
    Nothing -> []
    Just n  -> Map.toList n

data Front = Front !Int !(Map (Distance, Int) Vertex)

frontEmpty = Front 0 Map.empty

frontPush (v,d) (Front i m) = 
    Front (i+1) (Map.insert (d,i) v m)

frontPop (Front i m) = case Map.minViewWithKey m of
    Just ( ((d,_),v), m' ) -> Just ( (v,d), Front i m' )
    Nothing -> Nothing

distances :: Vertex -> Graph -> Map Vertex Distance
distances i g = go (frontPush (i,0) frontEmpty) Map.empty
  where
    go front distances = case frontPop front of
      Nothing -> distances
      Just ( (v,d), front' )
        | v `Map.member` distances -> go front' distances
        | otherwise ->
          let distances' = Map.insert v d distances
              n = [(v,d) | (v,d) <- neighbors v g, 
                           not (v `Map.member` distances')]
              front'' = foldl' (flip frontPush) front' n
          in go front'' distances'

Максимальный поток

Потоком в ориентированном графе называется отображение F из множества его рёбер в какое-нибудь числовое множество такое, что для любой вершины v верно следующее: сумма F от всевозможных рёбер с концом v равна сумме F от всевозможных рёбер с началом v.

Традиционно значения отображения F называются токами, а определение потока, грубо говоря, звучит как «в каждой вершине сумма входящих токов равна сумме исходящих». Типичный пример потока — токи в электрической цепи.

Некоторые комбинаторные задачи о распределении ресурсов (типичный пример: задача о максимальном паросочетании) сводятся к так называемой задаче о максимальном потоке. Постановка её следующая:

• на каждом ребре графа задаётся пропускная способность — действительное неотрицательное число
• в граф добавляется ребро E без пропускной способности; его начало называется стоком, а конец — источником
• рассматриваются только потоки с действительными неотрицательными токами, не превосходящими пропускные способности соответствующих рёбер
• требуется найти поток, максимизирующий ток на ребре E

Часто ребро E вообще исключают из разговора, говоря вместо «ток на E» просто «ток от источника до стока». В таком случае сумма выходящих из источника токов должна быть на этот самый ток больше суммы входящих в него токов. На этот же ток сумма входящих в сток токов должна быть больше суммы токов, выходящих из него.

Паросочетание как поток

Напомним, что паросочетанием называется набор непересекающихся рёбер в двудольном графе (говоря двудольный мы имеем в виду актуально двудольный — в котором вершины уже поделены на два класса, а рёбра бывают только вида «из одного класса в другой»).

Почти всегда на практике одной долей такого графа служит некоторый набор множеств, второй долей — объединение этих множеств, а рёбрами — принадлежности элементов множествам.

Если в двудольный граф добавить вершины A и B, из A пустить по одному ребру в каждую вершину первой доли, из каждой вершины второй доли пустить ребро в B, а затем каждому ребру назначить пропускную способность 1, то целочисленный поток от A до B в таком графе эквивалентен паросочетанию в исходном графе.

Переформулировка задачи

Решать задачу о максимальном потоке в том виде, в котором она поставлена, неудобно. Преобразуем её следующим образом: для каждой пары вершин исходного графа все рёбра, их соединяющие, объеденим в одно (произвольного направления) с нижним и верхним ограничениями на ток по нему (эти два ограничения — суммы пропускных способностей в двух возможных направлениях).

Поток в таком графе нетрудно преобразовать в поток в исходном, причём такое преобразование определено не однозначно, в отличие от «обратного». Как мы увидим позже, избавление от этой неоднозначности очень сильно упрощает задачу о максимальном потоке.

Алгоритм Форда-Фалкерсона

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

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

Начальный поток — нулевой.

Имея поток, следующий строится так: ищется простой (т.е. проходящий по каждому ребру не более одного раза) путь от источника к стоку в остаточном графе. Остаточный граф определён так: ребро от вершины a до вершины b есть тогда и только тогда, когда ток от a до b хотя бы на единицу меньше пропускной способности в соответствующую сторону. Например, если рассмотреть граф из одного ребра от a до b с ограничениями тока -1 и 1, то остаточный граф устроен так:

• если по ребру идёт ток -1, в остаточном есть только ребро от a до b
• если по ребру идёт ток 1, в остаточном есть только ребро от b до a
• если по ребру идёт ток 0, в остаточном графе есть рёбра как от a до b, так и от b до a

Если путей от источника до стока в остаточном графе нет, то алгоритм завершается. Если же такой путь есть, то токи вдоль него увеличиваются на единицу.

Алгоритм Форда-Фалкерсона не постулирует конкретный способ поиска пути в остаточном графе. Например, если такой поиск происходит «в ширину», то результирующий алгоритм носит имена Эдмондса и Карпа.

В качестве оптимизации можно на каждом шаге увеличивать ток не на единицу, а на минимальную ненулевую разницу между токами потока и соответствующими пропускными способностями рёбер. Заодно такой метод действует и для нецелочисленных потоков (с одной оговоркой — можно подобрать нецелые пропускные способности, приводящие к бесконечной последовательности потоков).

Почему оно работает?

Собственно, нужно доказать, что последний поток в последовательности даёт максимально возможный ток от источника до стока.

Самая главная хитрость доказательства — взглянуть на правильный набор рёбер. Например, подойдёт такой: один конец — вершина из множества тех, до которых можно дойти от источника по рёбрам остаточного графа; другой конец — вершина не из этого множества. Все такие рёбра обладают следующим свойством: на каждом достигнута максимальная пропускная способность, причём в направлении «от источника» (это напрямую следует из определения остаточного графа).

Осталось заметить два факта (оба следуют из определения потока):

• ток от источника к стоку равен суммарной пропускной способности рассмотренных нами рёбер
• ток от источника к стоку не может превосходить эту самую суммарную пропускную способность

Значит, мы нашли поток с максимально возможным током от источника к стоку.

Практическая реализация

В главе про обходы графов мы рассматривали вариант поиска в глубину, строящий путь (а не границу множества обработанных вершин):

std::set<int> component(Graph &g, int v) {
    std::vector<int> path{v};
    std::set<int> visited{};
    
    while (path.size() > 0) {
        int current_vertex = path.back();
        visited.insert(current_vertex);
        
        bool found_new = false;
        
        auto n = neighbors(g,current_vertex);
        for (int i : n) {
            if (visited.count(i) == 0) { 
                path.push_back(i);
                found_new = true; 
                break;
            }
        }
        
        if (found_new) { continue; }
        
        path.pop_back();
    }
    
    return visited;
}

Для нужд алгоритма Форда-Фалкерсона он неплохо подходит, но его нужно модифицировать для работы с остаточным графом.

Сначала модифицируем сам тип «граф»:

using Vertex = int;
using Capacity = int;
using Graph = std::map<Vertex,std::map<Vetrex,Capacity>>;

void addEdge(Graph &g, Vertex i, Vertex j, Capacity c) {
    g[i][j] += c;
    g[j][i];  // помечаем, что i -- сосед j
}

using Current = int;
using Flow = std::map<Vertex,std::map<Vertex,Current>>;

std::set<Vertex> neighbors(Graph &g, Flow &f, int v) {
    std::set<Vertex> res;

    for (auto x : g[v]) {
        // договоримся хранить ток на ребре (a,b) сразу в 
        // двух рёбрах так, чтобы: f[a][b] == -f[a][b]
        auto current = f[v][x.first];
        
        if (current < x.second) { res.insert(x.first); }
    }
    
    return res;
}

На а теперь — поиск пути в остаточном графе:

std::vector<Vertex> 
find_residual_path(Graph &g, Flow &f, Vertex src, Vertex dst) {
    std::vector<Vertex> path{src};
    std::set<Vertex> visited{};
    
    while (path.size() > 0) {
        Vertex current_vertex = path.back();
        
        // главное -- вовремя остановиться
        if (current_vertex == dst) { break; }
        
        visited.insert(current_vertex);
        
        bool found_new = false;
        
        auto n = neighbors(g,f,current_vertex);
        for (Vertex i : n) {
            if (visited.count(i) == 0) { 
                path.push_back(i);
                found_new = true; 
                break;
            }
        }
        
        if (found_new) { continue; }
        
        path.pop_back();
    }
    
    return path;  // теперь возвращаем путь, а не связную компоненту
}

Ну и, наконец, Форд-Фалкерсон:

Flow find_max_flow(Graph &g, Vertex src, Vertex dst) {
    Flow f;

    for (;;) {
        auto path = find_residual_path(g, f, src, dst);
        
        if (path.size() == 0) { break; }
        
        for (size_t i = 1; i < path.size(); i++) {
            Vertex u = path[i-1];
            Vertex v = path[i];
        
            f[u][v] += 1;
            f[v][u] -= 1;
        }
    }
    
    return f;
}

Кодирование информации

Кодировка Хаффмана

Рассмотрим задачу кодирования данных без потерь. Это означает следующее:

• есть два множества A и B, называемые входным и выходным алфавитами
• требуется построить биекцию f между множествами LA и LB конечных последовательностей элементов A и элементов B

Кодирование называется посимвольным, если существует отображение g: A -> LB, для которого f определено как конкатенация результатов g на элементах входной последовательности.

Говоря проще, посимвольное кодирование получается как результат монадного связывания функтора моноидального освобождения. В терминах языка Haskell:

encodeBy: (a -> [b]) -> [a] -> [b]
encodeBy f xs = xs >>= f

Посимвольное кодирование (порождённое отображением g) называется равномерным, если длины всех g-образов одинаковы. Равномерное кодирование называется минимальным, если длины g-образов минимально возможны (для данной пары алфавитов A и B).

Взглянем на наиболее распространённые текстовые кодировки: ASCII, UTF8, UTF16, UTF32. Все они посимвольные; ASCII и UTF32 ещё и равномерные. При этом, если считать, что выходной алфавит состоит только из 0 и 1, то используемая в реальности разновидность ASCII (кодирующая каждый символ восемью битами) не является минимальной (зато лишний восьмой бит можно использовать по своему усмотрению без потери возможности обратного преобразования).

Эксплуатация неравномерности частот

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

А именно, если кодировать более частые символы более короткими последовательностями, можно добиться меньшего математического ожидания длины сообщения, чем при минимальном равномерном кодировании. Соответственно, определим эффективность сжатия как отношение математических ожиданий длины кодировки одного символа при равномерном и при рассматриваемом кодировании. Например, эффективность сжатия настоящего (7-битного) ASCII по сравению с 8-битным, используемым в реальности, равна 8/7.

Кодировка Хаффмана — это как раз один из способов сжатия без потерь, эксплуатирующего неравномерность частот.

Дерево кодировки

Напомним, что если выходной алфавит состоит только из 0 и 1, то посимвольные кодировки принято изображать в виде дерева, на котором 0 соответствует перемещению в левый потомок, а 1 — перемещению в правый.

Например, на языке Haskell дерево посимвольной кодировки может выглядеть так:

data Encoding a = ENode (Maybe a) (Encoding a) (Encoding a) | ENull
type Bits = [Bool]

encodingTable :: Encoding a -> [(a, Bits)]
encodingTable e = encode [] e
  where
    encodeChar _ Nothing = []
    encodeChar code (Just c) = [(c, code)]
  
    encode _ ENull = []
    encode prefix (ENode mc l r) =
        encode (prefix ++ [False]) l ++ 
        encode (prefix ++ [True]) r ++
        encodeChar prefix mc

А на C++ — вот так:

template<typename T>
struct EncodingNode {
    bool hasData;
    T data;
    size_t left, right;
};

template<typename T>
using Encoding = std::vector<EncodingNode<T>>;

using Bits = std::vector<bool>;

template<typename T>
using EncodingTable = std::map<T, Bits>;

template<typename T>
EncodingTable<T> encoding_table(Encoding<T> &e) {
    return _encode(Bits{}, e, e.size() - 1);
}

template<typename T> EncodingTable<T> 
_encode(Bits prefix, Encoding<T> &e, size_t index) {
    EncodingTable<T> res;

    if (index == (size_t)(-1)) {
        return res;
    }
    
    if (e[index].hasData) {
        res[e[index].data] = prefix;
    }
    
    Bits left_prefix = prefix;
    left_prefix.push_back(false);
    auto leftTable = _encode(left_prefix, e, e[index].left);
    
    for (auto p: leftTable) { res[p.first] = p.second; }
    
    Bits right_prefix = prefix;
    right_prefix.push_back(true);
    auto rightTable = _encode(right_prefix, e, e[index].right);
    
    for (auto p: rightTable) { res[p.first] = p.second; }
    
    return res;
}

Этот код был лишь иллюстрацией того, как в C++ можно моделировать древообразные структуры без использования указателей (указатели в C++, конечно, полезная штука, но их очень трудно корректно напрямую использовать; особенно при необходимости обеспечить возможность глубокого копирования дерева). Больше мы не будем использовать в этой главе C++, тем более что уже по приведённому куску видно, что он является дословной трансляцией соответствующего кода на Haskell, только существенно более многословной.

Префиксные кодировки

Посимвольная кодировка называется префиксной, если ни один результат кодирования символа не является префиксом другого. Другими словами, в дереве кодировки символы могут находиться только в «листьях».

Префиксные кодировки легко декодируются:

-- у префиксных кодировок проще структура:
data Encoding a = ENode (Encoding a) (Encoding a) | ELeaf a | Enull
type Bits = [Bool]


encodingTable :: Encoding a -> [(a, Bits)]
encodingTable e = encode [] e
  where
    encode _ ENull = []
    encode prefix (ELeaf x) = [(x, prefix)]
    encode prefix (ENode l r) =
        encode (prefix ++ [False]) l ++ 
        encode (prefix ++ [True]) r
        
-- без обработки ошибок декодирования, зато ленивое
decode :: Encoding a -> Bits -> [a]
decode e bits = go e bits
  where
    go pos [] = []
    go ENull _ = [] -- дальше явно нельзя декодировать
    go (ELeaf x) bits = x : go e bits
    go (ENode l _) (False:rest) = go l rest
    go (ENode _ r) (True:rest)  = go r rest

Кодировка Хаффмана

Кодировка Хаффмана строится так (алгоритм здесь короче словесного описания):

import Data.Map as Map

type Frequency = Int

huffmanEncoding :: [(a,Frequency)] -> Encoding a
huffmanEncoding freqs = go forest
  where
    forest = Map.fromList $
        ( \(i,(x,f)) -> ((f,i),ELeaf x) ) <$> zip [0..] freqs
    
    go forest = case Map.minViewWithKey forest of
        Nothing -> ENull
        Just (((f1,_),e1), forest1) -> case Map.minViewWithKey forest1 of
            Nothing -> e1
            Just (((f2,i),e2), forest2) ->
                go $ Map.insert (f1+f2,i) (ENode e1 e2) forest2

Грубо говоря, изначально у нас есть для каждого символа тривиальное дерево, помеченное частотой этого символа. На каждом шаге выбираются два дерева, помеченные наименьшими частотами, и сливаются (с суммированием частот).

Конкретно в вышеприведённой программе метки-частоты были реализованы при помощи упорядоченного мультисловаря (с ключами-частотами). Мультисловарь был сделан из обычного словаря при помощи трюка с уникальными штампами (см. главу про алгоритм Дейкстры).

Кодировка Хаффмана является самой эффективной среди всевозможных префиксных кодировок: это означает, что матожидание (по всем символам) длины кода символа в ней минимально среди всех префиксных кодировок. Доказательство — несложное упражнение на перестроение дерева с полуинвариантом.

Преобразование BW

Методы решения уравнений

Одномерные нелинейные уравнения

Рассмотрим задачу нахождения корня уравнения

\[f(x) = a\]

для более-менее произвольной функции \(f:\mathbb{R}\to\mathbb{R}\) (если функция определена не на всей числовой прямой, доопределим её чем угодно).

Бинарный поиск

Предположим, что на отрезке \(S=[L,R]\) функция \(f\) обладает свойством промежуточных значений: для любых \(x,y\in S\) и любого \(M\) между \(f(x)\) и \(f(y)\) существует точка \(m\), расположенная между \(x\) и \(y\), для которой \(f(m)=M\).

Тогда если \(f(L)\) и \(f(R)\) находятся по разные стороны от \(a\), можно организовать бинарный поиск: переходить от отрезка \([x,y]\) к отрезку \([(x+y)/2,y]\) или \([x,(x+y)/2]\) в зависимости от значения \(f((x+y)/2)\) так, чтобы \(a\) по-прежнему находилось между значениями функции на краях отрезка.

Бинарный поиск работает всегда, причём весьма стабильно: с каждым шагом абсолютная погрешность уменьшается в два раза (абсолютной погрешностью мы в данном случае называем длину отрезка, на котором мы можем гарантировать наличие корня).

Метод хорды (секущей)

Для функций, близких к линейно-неоднородным (линейно-неоднородные — функции вида \(f(x)=kx+b\)), есть существенно более быстрая разновидность бинарного поиска: идея в том, чтобы смотреть не на середину отрезка, а на пересечение «хорды», соединяющей крайние точки графика на текущем отрезке, с прямой, вертикальная координата которой равна \(a\).

Если мы сейчас находимся на отрезке \([x,y]\), то прямая, проходящая через точки \((x,f(x)))) и \((y,f(y))\) является графиком отображения

\[ t\mapsto \frac{t-y}{x-y} f(x) + \frac{t-x}{y-x} f(y) \]

Нетрудно решить уравнение

\[ \frac{t-y}{x-y} f(x) + \frac{t-x}{y-x} f(y) = a \]

Оно равносильно

\[ (t-y) f(x) - (t-x) f(y) = a(x-y) \]

Это, в свою очередь, равносильно

\[ t (f(x) - f(y)) = a(x-y) + y,f(x) - x,f(y) \]

То есть искомое значение новой «середины» отрезка равно

\[ t = \frac{x-y}{f(x)-f(y)} a + \frac{y,f(x) - x,f(y)}{f(x)-f(y)} \]

Отметим, что метод секущих способен проявлять очень плохую сходимость для функций с большим разбросом значений производной.

Метод Ньютона (касательной)

Одним из самых популярных методов решения уравнений является метод Ньютона. В отличие от вышеприведённых методов, он имеет дело не с последовательностью отрезков, а с последовательностью точек — приближений к искомому корню.

Описание

Каждая следующая точка получается как пересечение касательной, построенной к функции в текущей точке, и прямой с вертикальной координатой \(a\). Чтобы не усложнять анализ лишними выражениями, будем считать без ограничения общности \(a=0\) (то есть, формально говоря, будем вместо \(f\) рассматривать отображение \((x\mapsto f(x)-a)\)).

Производную далее будем обозначать буквой \(d\). Касательная, построенная к \(f\) в точке \(x\), является графиком отображения

\[ t \mapsto df(x)\cdot (t - x) + f(x) \]

Пересечение с нулевым значением получить нетрудно:

\[ t = x - \frac{f(x)}{df(x)} \]

Оценка

Оценим сходимость метода Ньютона. Нам понадобятся дополнительные предположения: пусть мы изначально находимся левее корня, при этом на промежутке от нас до корня функция убывает (то есть её первая производная отрицательна) и выпукла вниз (то есть её вторая производная неотрицательна).

Выпуклость обеспечивает нам то, что все последующие приближения будут левее корня. Абсолютной погрешностью теперь будем называть разность между корнем и таким приближением. Сам корень обозначим буквой \(r\), а \(n\)-е приближение к нему обозначим \(r_n\). Начинаем с \(r_0\) и действуем согласно соотношению

\[ r_{n+1} = r_n - \frac{f(r_n)}{df(r_n)} \]

Абсолютную погрешность \(n\)-го приближения обозначим \(a_n\). Нетрудно заметить (вычтя \(r\) из обеих частей соотношения), что

\[ a_{n+1} = a_n - \frac{f(r_n)}{|df(r_n)|} \]

Заметим теперь, что согласно теореме Лагранжа есть точка \(m_n\) на промежутке между \(r_n\) и \(r\), для которой

\[ f(r_n) = f(r_n) - f(r) = df(m_n)(r_n - r) = |df(m_n)| a_n \]

Отсюда получаем

\[ a_{n+1} = a_n \left( 1 - \frac{|df(m_n)|}{|df(r_n)|} \right) \]

или, что то же самое (поскольку обе производные отрицательны)

\[ a_{n+1} = a_n \left( 1 - \frac{df(m_n)}{df(r_n)} \right) \]

Линейная сходимость

Учитывая неравенства \(|df(r)|\leqslant |df(m_n)|\) и \(|df(r_0)|\geqslant |df(r_n)|\), делаем вывод, что

\[ a_{n+1} = a_n \left( 1 - \frac{|df(m_n)|}{|df(r_n)|} \right) \leqslant a_n \left( 1 - \frac{|df(r)|}{|df(r_0)|} \right) \]

Это означает, что \( a_{n} \leqslant \alpha^n a_{0} \) для некоторой положительной константы \(\alpha\lt1\).

Если ввести относительную погрешность \( \delta_n = a_n / a_0 \), то видно, что \( \delta_{n+1} \leqslant \alpha \delta_n \).

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

В более приземлённых терминах: на каждом шаге количество известных нам цифр корня увеличивается на одну и ту же величину.

Квадратичная сходимость

Теперь докажем, что метод Ньютона сходится быстрее, чем линейно.

Ещё раз применим теорему Лагранжа, только на этот раз — для производной:

\[ df(m_n) = d^2f(s_n) (m_n - r_n) + df(r_n) \]

Отметим, что точка \(s_n\) находится между \(r_n\) и \(m_n\).

С учётом этого равенства соотношение для абсолютных погрешностей принимает вид

\[ a_{n+1} = a_n \left( 1 - \frac{df(m_n)}{df(r_n)} \right) = a_n \frac{|d^2f(s_n) (r_n-m_n)|}{|df(r_n)|} \]

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

\[ a_{n+1} \leqslant a_n^2 \frac{M}{|df(r)|} \]

Чтобы эта оценка стала понятнее, «обезразмерим» её, поделив на \(a_0\) и введя обозначение для «относительной» погрешности \(\delta_n = a_n/a_0\). Получим

\[ \delta_{n+1} \leqslant \delta_n^2 \frac{Ma_0}{|df(r)|} \]

Заметим, что если мы начали с такого момента, что

\[ \frac{Ma_0}{|df(r)|} \lt 1 \]

то \(\delta_{n+1} \lt \delta_{n}^2\). Это означает, грубо говоря, что с каждым следующим приближением количество известных нам цифр удваивается.

Также заметим, что если мы начали не с такого момента, то до такого момента рано или поздно дойдём: это сразу следует из доказанной выше линейной сходимости метода Ньютона.

Линейные уравнения

Вектора

Векторным пространством называется любая алгебраическая система \(V\) с двумя операциями: сложением типа \(V\times V\to V\) и растяжением вида \(R\times V\to V\) для некоторого кольца \(R\) (далее будем считать, что \(R\) — поле действительных чисел). Эти операции должны удовлетворять некоторому списку довольно естественных свойств. А именно, сложение должно иметь структуру коммутативной группы, а растяжение должно быть связано с ним следующими свойствами:

• \( (\alpha+\beta)(a+b) = \alpha a + \beta a + \alpha b + \beta b \)
• \( 1a = a \)
• \( (\alpha\beta)a = \alpha(\beta a) \)

Говоря проще, векторное пространство — то же самое, что и морфизм кольца \(R\) в кольцо эндоморфизмов коммутативной группы.

Чаще всего на практике встречаются пространства \(R^n\): последовательности чисел из \(R\) длины \(n\), сложение и растяжение которых определяется покомпонентно. Договоримся обозначать \(\delta_i\) такую последовательность (неважно для какого именно \(n\)), в которой только \(i\)-я компонента равна 1, а остальные — нулевые. Пространства \(R^n\) называются конечномерными.

Также часто встречаются пространства отображений типа \(X\to V\), где \(V\) является векторным пространством. Сложение и растяжение таких отображений определяется поточечно. В принципе, пространства \(R^n\) являются частным случаем таких пространств для множества \(X\) мощности \(n\).

Линейные отображения

Отображение векторных пространств называется линейным, если оно коммутирует со сложением и растяжением. Говоря проще, \(f\) линейно, если

\[ f(\alpha x + \beta y) = \alpha f(x) + \beta f(y) \]

Ну или, если быть совсем кратким:

\[ f\left( \sum_i \alpha_i x_i \right) = \alpha_i \sum_i f(x_i) \]

Линейные отображения конечномерных пространств имеют очень простое описание: любой вектор конечномерного пространства допускает однозначное представление в виде суммы \(\delta_i\). Поэтому если \(x = \sum_i \alpha_i \delta_i\), то

\[ f(x) = f\left( \sum_i \alpha_i \delta_i \right) = \alpha_i \sum_i f(\delta_i) \]

Заметим, что каждый из \(f(\delta_i)\) можно также разложить по дельтам. Введём обозначение \(f_{ji}\) для коэффициентов этого разложения:

\[ f(\delta_i) = \sum_i f_{ji} \delta_j \]

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

\[ f(x)_{i} = \sum_j f_{ij} x_j \]

и для коэффициентов композиции линейных отображений

\[ (f\circ g)_{ij} = \sum_k f_{ik} g_{kj} \]

Набор коэффициентов линейного отображения называется его матрицей. Обычно её записывают в виде таблицы:

\[ \left( \begin{matrix} f_{11} & \ldots & f_{m1}\\ \vdots & & \vdots\\ f_{n1} & \ldots & f_{mn} \end{matrix} \right) \]

Вектора же часто записывают в виде матрицы-столбца: это связано с тем, что любому вектору \(x: R^m\) можно однозначно сопоставить линейное отображение \(\ulcorner x\urcorner: R\to R^m \) по правилу

\[ \ulcorner x\urcorner (r) = rx \]

Такое соглашение позволяет вообще избавиться от операции подстановки значения в отображение (она плоха тем, что неассоциативна, поэтому требует внимания к тому, как расставлены скобки) и писать все равенства в терминах композиции отображений (композиция ассоциативна). Далее мы будем придерживаться именно этого соглашения, считая, что в записи \(AB\) между \(A\) и \(B\) стоит композиция. Соответственно, числа мы будем также считать линейными отображениями \(R\to R\).

Линейные уравнения

Линейным уравнением называется уравнение вида

\[ Ax = a \]

где \(A:U\to V\) — линейное отображение, \(a:V\) — известный вектор (называемый правой частью), а \(x:U\) — неизвестный вектор, который требуется найти.

Если пространства \(U\) и \(V\) конечномерны, то можно применить метод Гаусса. С бесконечномерными такой фокус не пройдёт. Также не пройдёт он с пространствами большой размерности — алгоритмическая сложность метода Гаусса растёт как куб размерности пространства.

Поэтому мы рассмотрим так называемый итерационный способ решения линейных уравнений.

Итерация для чисел

Хорошо известно, что значение \( (1-\alpha)^{-1} \) при \( |\alpha| \lt 1 \) легко найти при помощи итерации схемы Горнера:

\[ x_0 = 0;\qquad x_{n+1} = 1 + \alpha x_n \]

Известно это из курса арифметики, в котором доказывается несложная теорема о том, что

\[ x_n = \frac{1-\alpha^n}{1-\alpha} \]

откуда следует, что при \( |\alpha| \lt 1 \) выполнено

\[ \lim_n x_n = (1-\alpha)^{-1} \]

В принципе, любое числовое уравнение вида

\[ kx = a \]

можно представить в таком виде, домножив на \( \alpha \) и сделав тривиальное преобразование

\[ (1 - (1 - \alpha k)) x = \alpha a \]

Очевидно, что выбором \( \alpha \) можно загнать \( \alpha k \) в интервал от 0 до 1.

Величина линейного отображения

С линейными отображениями всё хуже — у них нет адекватного аналога абсолютной величины числа. Самое лучшее, что есть — это точная нижняя грань множества тех чисел \(M\), для каждого из которых неравенство

\[ |Ax| \leqslant M |x| \]

выполнено для любого вектора \(x\). Она же совпадает с точной верхней гранью множества всевозможных чисел вида \(|Ax| / |x|\). Будем её в дальнейшем обозначать \( |A| \) и называть нормой отображения \( A \). Она удовлетворяет:

• неравенству треугольника \( |A+B| \leqslant |A| + |B| \)
• мультипликативному свойству \( |AB| \leqslant |A||B| \)

Именно эти свойства обеспечивают сходимость итерационного процесса, рассмотренного выше (конечно, при норме соответствующего отображения, меньшей 1).

Можно также сделать оценку снизу и рассмотреть точную верхнюю грань множества тех чисел \(N\), для каждого из которых для любого \(x\) имеет место неравенство

\[ N|x| \leqslant |Ax| \]

Эта точная грань обладает похожим мультипликавным свойством, но уже не удовлетворяет неравенству треугольника: в случае

\[ A = \left(\begin{matrix}1&0\\ 0&2\end{matrix}\right)\qquad B = \left(\begin{matrix}2&0\\ 0&1\end{matrix}\right) \]

эта величина у \(A\) и у \(B\) равна 1, а у их суммы она равна 3. Те отображения, у которых она положительна, будем называть положительными, а саму эту грань — коэффициентом положительности.

Всё плохо?

Посмотрим внимательно на уравнение

\[ (1 - (1 - \alpha A)) x = \alpha a \]

Пусть

\[ A = \left(\begin{matrix}1+a&0\\ 0&1-a\end{matrix}\right) \]

Нетрудно заметить, что как ни уменьшай \(\alpha\), норма \( (1 - \alpha A) \) всё равно равна \( 1 + |\alpha a| \), что больше единицы.

Также нетрудно заметить, что при \(U\ne V\) запись вида \(1 - \alpha A\) вообще смысла не имеет (при \(U=V\) же можно считать, что единица — это тождественное отображение).

На самом деле, всё хорошо (почти)

Вспомним, что в пространствах \( R^n \) есть стандартное скалярное произведение:

\[ x\cdot y = \sum_i x_i y_i \]

Пусть \(A: R^m\to R^n \) — линейное отображение с матрицей \( i,j\mapsto (A_{ij}) \). Отображение с матрицей \( i,j\mapsto (A_{ji}) \) будем называть транспонированным к \(A\) и обозначать \(A^*\).

Нетрудно проверить, что:

\[ x\cdot y = x^* y \]

а также, что:

\[ (AB)^* = B^* A^* \]

Отсюда (и из неравенства Коши-Буняковского) немедленно следует

\[ |Ax|^2 = x^{*} A^{*} A x = x\cdot (A^{*}A x) \leqslant |x| |A^{*}| |A| |x| = |A^{*}| |A| |x|^2 \]

Это означает, что \( |A|^2 \leqslant |A^{*}| |A| \), что эквивалентно \( |A| \leqslant |A^{*}| \). Меняя в этом рассуждении \(A\) на \(A^{*}\), получаем \( |A^{*}| \leqslant |A| \), откуда

\[ |A| = |A^{*}| \]

Использование транспонированного отображения позволяет построить сходящийся итерационный процесс для уравнения \(Ax = a\).

Что делать?

Банально перейти от уравнения

\[ Ax = a \]

к уравнению

\[ \alpha A^{*} A x = \alpha A^{*} a \]

и выбрать достаточно малый \(\alpha\). При этом итерационный процесс примет вид

\[ x_0 = b;\qquad x_{n+1} = b + Bx_n \]

где \(b = \alpha A^{*} a\) и \( B = 1 - \alpha A^{*} A \).

Скорость сходимости

У нас имеется точное равенство

\[ \alpha A^{*}A x_n = (1 - \left( \alpha A^{*} A\right)^n ) b \]

Пусть \(x\) — точное решение уравнения (если оно вообще есть). Тогда

\[ A^{*}A (x-x_n) = \alpha^n \left(A^{*} A\right)^n A^{*}a \]

Пусть \(N\) — коэффициент положительности для \(A^{*}A\). Если он не равен нулю, то

\[ |x-x_n| \leqslant \frac{|A^{*}a|}{N} \left( \alpha |A|^2 \right)^n \]

Отсюда видно, что если выбрать \( \alpha \lt |A|^{-2} \), то процесс сойдётся. Единственная проблема: непонятно, когда. Для этого нужно оценивать сверху \(N^{-1}\), что представляет собой нетривиальную задачу (по сложности сопоставимую с решением исходного уравнения). Впрочем, иногда эту оценку можно вытянуть из той модели, для исследования свойств которой это уравнение и решается.

А вот выбор \( \alpha \) осуществить легко: \( |A|^2 \) оценивается сверху суммой квадратов компонент матрицы \( A \). Таким образом получаем нижнюю оценку на \( |A|^{-2} \). Соответственно, \( \alpha \) можно выбрать, например, равным половине этой оценки.

Дифференциальные уравнения

Обыкновенное дифференциальное уравнение

Пусть \(V\) — (конечномерное) векторное пространство, а \(T\) — подмножество \(\mathbb{R}\). Обыкновенным дифференциальным уравнением (ОДУ) называется задача поиска неизвестной функции \(x: T\to V\), удовлетворяющей соотношению, в котором фигурируют её производные. Обычно считают, что это соотношение имеет вид

\[ \dot{x}(t) = f(x(t),t) \]

где \(f:V\times T \to V\). Но иногда приходится иметь дело и с неявными уравнениями вида

\[ f(\dot{x}(t), x(t), t) = 0 \]

для \(f:V\times V\times T \to V\).

Задачей Коши называется ОДУ с начальным условием: дополнительным соотношением вида \(x(0) = x_0\). Обычно задача Коши имеет единственное решение (хотя бы в некоторой окрестности начального момента времени).

Уравнения высших порядков

Если в уравнении фигурируют высшие производные (хотя бы второго порядка), то его можно свести к уравнению первого порядка путём повышения размерности. Для этого для каждой из промежуточных производных вводится новая неизвестная с очевидным уравнением.

Например, уравнение

\[ \ddot{x}(t) = f(\dot{x}(t), x(t), t) \]

может быть преобразовано к виду \[ \begin{aligned} \dot{v}(t) & = f(v(t), x(t), t) \\ \dot{x}(t) & = v(t) \end{aligned} \]

Поэтому далее будем считать, что любое ОДУ является ОДУ первого порядка.

Автономные уравенения

ОДУ называется автономным, если оно имеет вид

\[ \dot{x}(t) = f(x(t)) \]

Опять же, любое уравнение можно сделать автономным, если считать время неизвестной функцией, заданной уравнением

\[ \dot{t} = 1 \]

и начальным условием \(t(0) = 0\). Далее будем считать любое ОДУ автономным.

Явный метод Эйлера

Самым простым способом численного решения дифференциальных уравнений является явный метод Эйлера: выбирается шаг по времени \(\tau\), после чего применяется итерационный процесс

\[ x_{n+1} = x_{n} + f(x_n) \tau \]

Вектор \(x_n\) соответствует значению функции \(x\) в момент времени \(t_n = n\cdot\tau\).

К сожалению, явный метод Эйлера (как и разнообразные его модификации, в том числе часто применяемые методы Рунге-Кутта) не слишком пригодны для моделирования физических систем на больших промежутках времени: они не сохраняют энергию системы.

Полунеявный метод Эйлера

Рассмотрим типичную механическую систему

\[ \ddot{x}(t) = F(x(t)) \]

с начальными условиями \(x(0) = x_0\), \(\dot{x}(0) = v_0\).

Преобразуем её в ОДУ первого порядка:

\[ \begin{aligned} \dot{v}(t) & = F(x(t)) \\ \dot{x}(t) & = v(t) \end{aligned} \]

Явный метод Эйлера для такого ОДУ имеет вид:

\[ \begin{aligned} v_{n+1} & = v_n + F(x_n)\cdot\tau \\ x_{n+1} & = x_n + v_n\cdot\tau \end{aligned} \]

У полунеявного метода ровно одно отличие (во втором уравнении):

\[ \begin{aligned} v_{n+1} & = v_n + F(x_n)\cdot\tau \\ x_{n+1} & = x_n + v_{n+1}\cdot\tau \end{aligned} \]

То есть теперь следующее состояние вычисляется не по предыдущему, а частично — по самому себе. Засчёт такой обратной связи полученный численный метод становится существенно стабильнее.

Метод конечных элементов

Рассмотренный в прошлой главе метод Эйлера относится к классу «методов конечных разностей»: они заменяют производную на разность значений в близкие моменты времени, получая приближённое уравнение, для которого затем ищется точное (по модулю вычислительных погрешностей) решение.

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

Кусочно-линейная аппроксимация

Рассмотрим для примера уравнение вида

\[ \dot{x}(t) = k x(t) \]

с начальным условием \(x(0) = A\).

Простейшей аппроксимацией решения является кусочно-линейная: отрезок от \(0\) до \(T\) разбивается на \(N\) частей, на каждой части \(x\) приближается линейной неоднородной функцией. Такая аппроксимация почти нигде не удовлетворяет уравнению, но если потребовать, чтобы точное равенство достигалось в граничных точках (для правой односторонней производной), получится метод Эйлера.

Метод конечных элементов использует другой критерий близости: вместо того, чтобы требовать точное равенство в нескольких точках, он требует равенства взвешенных средних значений вблизи этих точек.

Обобщённая постановка

Если про две функции \(f\) и \(g\) известно, что они достаточно хорошие (например, непрерывные), то из равенства интегралов

\[ \int_a^b f = \int_a^b g \]

для всевозможных \(a\) и \(b\) следует равенство самих функций.

Более того, если функция \(f\) выражает какую-то реальную величину, не существует способа измерить её значение в какой-либо точке. Максимум, на что можно надеяться — среднее значение на некотором отрезке:

\[ \frac{\int_a^b f}{b-a} \]

Поэтому довольно естественно перейти от точного равенства функций к равенству их средних значений на произвольных отрезках. Чуть менее естественно от равенства средних значений перейти к равенству взвешенных средних на всей области определения:

\[ \int_0^T fw = \int_0^T gw \]

Такое равенство формально является более сильным — равенство средних получается, если допускать только кусочно-постоянные весовые функции \(w\). Но оно, опять же, формально более слабое, чем точное равенство. Далее будем такое равенство называть обобщённым.

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

\[ \int_0^T \dot{x}w = \int_0^T kxw \]

может быть преобразовано при помощи правила Лейбница и формулы Ньютона-Лейбница следующим образом:

\[ x(T)w(T) - Aw(0) = \int_0^T (x\dot{w} + kxw) \]

У этого уравнения есть две интересных (и полезных) особенности:

1. В него уже включено граничное условие \( x(0) = A \).
2. В нём не встречается \( \dot{x} \): оно имеет смысл даже для аппроксимаций с плохо ведущей себя производной.

Дискретизация

Любая кусочно-линейная функция может быть записана в виде линейной комбинации «зубьев». Договоримся, что \(i\)-й «зуб» — это функция \(\varphi_i\):

• равная 0 везде, кроме отрезка с центром \(iT/N\) и радиусом \(T/N\)
• линейно меняющаяся от 0 до 1 на левой половине этого отрезка
• линейно меняющаяся от 1 до 0 на правой половине этого отрезка

Всего зубьев \(N+1\) — с нулевого по \(N\)-й. Обратите внимание, что зубья с номерами \(0\) и \(N\) «половинчатые».

Пусть \(x = \sum_{i=0}^N \alpha_i \varphi_i\), где \(\alpha_i\) — неизвестные числа. Мы знаем, что \(\alpha_0 = A\). Остальные \(N\) штук нам нужно найти. Соответственно, нам нужно \(N\) уравнений.

Обычно эти уравнения получают, подставляя в качестве весовой функции те же зубья, на основе которых строится решение \(x\). В нашей ситуации отлично подойдут первые \(N\) зубьев (кроме самого последнего): для всех них верно \(w(T) = 0\), что избавляет от лишнего члена в уравнении.

Таким образом мы получаем систему из \(N\) уравнений вида

\[ -A\left(\varphi_i(0) + \int_0^T \varphi_0 (\dot{\varphi}_i + k\varphi_i) \right) = \sum_{j=1}^{N} \alpha_j \int_0^T \varphi_j(\dot{\varphi}_i + k\varphi_i) \]

для \(i\) от \(0\) до \(N-1\).

ЕГЭ

Здесь с некоторой периодичностью будут появляться статьи, посвящённые тем или иным задачам ЕГЭ, с которыми могут возникнуть проблемы при решении этого самого ЕГЭ.

Задача 25

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

1. найти все элементы, удовлетворяющие некоторому предикату P
2. посчитать композицию найденных элементов в рамках некоторой полугруппы
3. напечатать результат
4. либо найти все элементы, удовлетворяющие предикату Q
5. преобразовать их некоторой унарной операцией, построенной по результату
6. напечатать преобразованный массив

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

Вариативность предикатов

Определённую сложность представляет выражение предикатов P и Q на языке программирования. Поэтому сейчас мы разберём самые частовстречающиеся предикаты и их выражения на языке C++.

Договоримся, что массив называется a, а индекс того элемента, на который мы сейчас смотрим (в цикле) — i.

Делимость на что-нибудь

Бывает двух видов: делимость элемента и делимость его номера. Их не нужно путать. Найти чётные элементы:

a[i] % 2 == 0

Найти элементы на чётных местах:

i % 2 == 0

Неделимость на что-нибудь

Здесь нужно быть аккуратным из-за специфического поведения оператора взятия остатка в большинстве языков программирования. Например, на C++ нельзя нечётность проверять конструкцией вида

a[i] % 2 == 1

Она будет работать только для неотрицательных a[i]! Правильно же выражать высказывание «элемент имеет остаток r при делении на d» конструкцией вида

(a[i] - r) % d == 0

Цифры в записи числа

Например, задание вида «найти все числа, у которых совпадает последняя цифра в десятичной и шестнадцатиричной записи». Задание не представляет сложности, если помнить определение d-ичной системы счисления.

Напомним, что система счисления с основанием d состоит из:

• множества \(D\) цифр
• биекции \(v:D\to\mathbb{N}[0,d)\) между цифрами и натуральными числами от нуля включительно до \(d\) исключительно (эта биекция называется числовым значением цифр)
• отображения \(r:[D]\to\mathbb{N}\) из множества конечных последовательностей цифр в множество натуральных чисел
• отображения \(w:\mathbb{N}\to[D]\)

При этом должно выполняться следующее: во-первых

\[ r(x) = \sum_i v(x_i) d^i \]

Во-вторых, \(r(w(n)) = n\).

Из этого определения, например, сразу следует, что \(w(n)_0=n\bmod d\). Это как раз и есть числовое значение «последней цифры числа \(n\)». Поэтому задание «найти все числа, у которых совпадает последняя цифра в десятичной и шестнадцатиричной записи» можно решить условием вида

a[i] % 10 == a[i] % 16

Количество цифр в записи числа

Например, требуется найти все числа, которые записываются тремя цифрами в шестнадцатиричной системе счисления. Лучше всего подобные условия выражать банальными неравенствами вида:

16*16 <= a[i] and a[i] < 16*16*16

Сумма цифр в записи числа

Встречается очень редко, но встречается. Требует две дополнительных переменных и вложенный цикл:

s = 0;
k = a[i];
while (k > 0) {
    s += k % 10;  // пример для 10-чной системы счисления
    k  = k / 10;
}
// теперь в переменной s находится сумма десятичных цифр числа a[i]

Взаимоотношения с соседями

Например, нужно найти все элементы, большие обоих своих соседей. Единственная тонкость такого задания — границы цикла: вместо стандартного

for (i = 0; i < N; i++) {
    // нечто
}

нужно писать

for (i = 1; i < N-1; i++) {
    if (a[i] > a[i-1] and a[i] > a[i+1]) {  
        // что-то
    }
}

Вариативность полугруппы

Если Вы вдруг забыли, что такое полугруппа, и ещё не подсмотрели нигде до этого момента, напомним, что полугруппой с носителем \(X\) называется бинарная ассоциативная операция на множестве \(X\). А именно, отображение \(f:X\times X\to X\), для которого выполнено

\[ f(x,f(y,z)) = f(f(x,y),z) \]

или, в инфиксной записи,

\[ x f (y f z) = (x f y) f z \]

В качестве полугруппы в 25-й задаче почти всегда фигурирует одно из трёх: минимум, максимум, сумма. Причём иногда дополнительно задаётся отображение из множества целых чисел в носитель полугруппы (в дальнейшем будем его называть препроцессором).

Например, в задании «найти количество чётных элементов массива» полугруппой является сумма целых чисел, а препроцессором — отображение из множества целых чисел в себя, отображающее все целые числа в единицу.

Следует очень внимательно следить за двумя вещами:

1. часто в случае, если в массиве нет элементов, удовлетворяющих предикату P, требуется сделать нечто; это самое нечто никаких сложностей не представляет, но про него нужно не забыть
2. при поиске минимума/максимума требуется либо бинарный флаг вида «ещё не встретили ни одного хорошего числа», либо удобное начальное значение

Например, если требуется найти минимальное чётное число в массиве, следует прибегнуть к конструкции вида

encountered_good = false;
for (i = 0; i < N; i++) {
    if (a[i] % 2 == 0) {
        m = not encountered_good ? a[i] : min(m,a[i]);
        encountered_good = true;
    }
}

if (not encountered_good) {
    // сделать то самое нечто
    return 0;
}

с поправкой на ограничение по использованию переменных. Но в реальности в задании всегда указано ограничение сверху на величину чисел в массиве (например, 10000), а вот лишняя переменная под флаг недоступна. В таком случае работает конструкция вида

m = 10001; // значение на 1 больше максимально возможного значения
for (i = 0; i < N; i++) {
    if (a[i] % 2 == 0) { m = min(m,a[i]); }
}

if (m > 10000) {
    // сделать то самое нечто
    return 0;
}

Примерный вид программы

Рассмотрим, например такое задание:

1. найти максимальное число в массиве целых положительных чисел, в семиричной записи которого две или три цифры
2. уменьшить все чётные элементы массива на найденное число или на 1, если чисел указанного вида не было
3. напечатать получившийся массив

В программе дана константа N (длина массива), массив a и целочисленные переменные i и m (обычно дают на всякий случай ещё одну переменную, но не нужно на это рассчитывать).

m = 0;
for (i = 0; i < N; i++) {
    if (7 <= a[i] and a[i] < 7*7*7) {
        m = max(m,a[i]);
    }
}

if (m == 0) { m = 1; }

for (i = 0; i < N; i++) {
    if (a[i] % 2 == 0) { a[i] -= m; }
}

for (i = 0; i < N; i++) {
    cout << a[i] << endl;  // не забываем про отступы между числами
    // Плюс не забываем уточнить, чем _именно_ нужно разделять элементы
    // массива при печати.
}

Задача 27

А вот эта задача очень простая. Бывает трёх видов (сводящихся к одной и той же схеме):

• выбрать из кучи пар чисел по одному числу так, чтобы что-нибудь максимизировать или минимизировать
• найти самый частовстречающийся результат какого-нибудь отображения, применённого к каждому элементу массива
• редуцировать (в какой-нибудь полугруппе) все пары элементов массива, индексы которых в этом массиве различаются не менее чем на заданную величину

И то, и другое нужно либо сделать хоть как-то (тогда за задачу дадут 2 балла), либо за один проход по входным данным, не сохраняя их в массив (тогда за задачу дадут 4 балла).

Редукция

Самая частая разновидность.

Типичное задание: найти количество пар разных элементов входного массива, сумма которых делится на 6. Обратите внимание на то, что в конкретно этой задаче неявно присутствует минимальное ограничение на разность индексов элементов пары, равное 1.

#include <iostream>
#include <deque>
#include <map>

using namespace std;

int main() {
    unsigned N; cin >> N; // входной размер массива

    deque<unsigned> queue; // задержка между получением элементов и обработкой
    map<unsigned,unsigned> mod_count; // таблица количеств по остаткам
    unsigned answer = 0;
    for (unsigned i = 0; i < N; i++) {
        unsigned x; cin >> x; 
        queue.push_back(x);

        if (i < 1) { continue; }

        mod_count[queue.front()%6] += 1; 
        queue.pop_front();

        unsigned mod = queue.back()%6;
        answer += mod_count[(6 - mod) % 6];
    }

    cout << answer << endl;
}

Может показаться, что очередь длины 1 — весьма странное решение, но в некоторых ситуациях любые альтернативы (кроме банальных замен deque на две переменных) оказываются более сложными. Всегда есть соблазн (если в условии ничего явно не сказано при минимальную разность индексов) просто завести несколько переменных, в каждую из которых редуцировать неким образом все элементы массива, а затем написать одну большую формулу, выражающую ответ в терминах этих переменных. Такой соблазн лучше игнорировать.

Рассмотрим ещё одно похожее задание: найти максимальное произведение пары элементов входного массива (все элементы — целые положительные числа), индексы которых различаются хотя бы на 6, делящееся на 15.

#include <iostream>
#include <deque>
#include <map>

using namespace std;

int main() {
    unsigned N; cin >> N;

    deque<unsigned> queue;
    map<unsigned,unsigned> mod_max;  // таблица максимумов по остаткам
    unsigned answer = 0;
    for (unsigned i = 0; i < N; i++) {
        unsigned x; cin >> x; 
        queue.push_back(x);

        if (i < 6) { continue; }  // здесь поменялась длина очереди

        unsigned mod = queue.front() % 15;
        mod_max[mod] = max(mod_max[mod], queue.front());  // сумма сменилась на максимум
        queue.pop_front();

        mod = queue.back()%15;
        if (mod == 0) {
            for (unsigned m = 0; m < 15; m+=1) { answer = max(answer, queue.back()*mod_max[m]); }
        } else if (mod % 3 == 0) {
            for (unsigned m = 0; m < 15; m+=5) { answer = max(answer, queue.back()*mod_max[m]); }
        } else if (mod % 5 == 0) {
            for (unsigned m = 0; m < 15; m+=3) { answer = max(answer, queue.back()*mod_max[m]); }
        }
    }

    cout << answer << endl;
}

Таблица частот

Требуется, например, найти и напечатать самую частую первую цифру десятичной записи чисел. Задание весьма тривиально, если завести таблицу частот цифр.

#include <iostream>
#include <map>

using namespace std;

unsigned first_digit(unsigned x) {
    while (x > 9) { x /= 10; }
    return x;
}


int main() {
    unsigned N; cin >> N;

    map<unsigned,unsigned> digit_freq;
    for (unsigned i = 0; i < N; i++) {
        unsigned x; cin >> x; 
        digit_freq[first_digit(x)] += 1;
    }

    // дальнейшее написано в предположении, что 
    // если несколько частот совпадают, то нужно вывести 
    // наибольшую из цифр
    unsigned best_digit = 0;
    unsigned best_freq = digit_freq[0];
    for (unsigned digit = 1; digit <= 9; digit++) {
        if (digit_freq[digit] >= best_freq) {
            // если бы требовалось вывести наименьшую из цифр, 
            // достаточно было бы поменять знак на >
            best_freq = digit_freq[digit];
            best_digit = digit;
        }
    }

    cout << best_digit << endl;
}

Выбор одного из двух

Типичное задание таково. На вход подаются пары (целых положительных) чисел. Нужно выбрать из каждой пары одно число так, чтобы:

• сумма выбранных чисел не делится на 42
• она является максимальной среди возможных

Если такой выбор сделать нельзя, требуется напечатать 0. В противном случае нужно напечатать сумму.

Ключевой момент заключается в том, что:

• либо подходит сумма наибольших (в каждой паре) чисел
• либо в одной из пар можно взять вместо наибольшего числа наименьшее
• либо сумма наибольших чисел и все разности в парах делятся на 42 и, соответственно, любой выбор чисел имеет сумму, делящуюся на 42

Это приводит к следующему решению (очень напоминающему решение типичной 25 задачи):

#include <iostream>

using namespace std;

int main() {
    unsigned N; cin >> N;

    unsigned large = 10001; // ограничение сверху плюс один
    unsigned total = 0;  // сумма наибольших
    unsigned min_diff = large; // минимальная из разностей, не делящихся на 42
    for (unsigned i = 0; i < N; i++) {
        unsigned x,y; cin >> x >> y; 

        total += max(x,y);
        
        unsigned diff = max(x,y) - min(x,y);
        if (diff % 42 != 0) {
            min_diff = min(min_diff, diff);
        }
    }

    unsigned answer = total;
    if (answer % 42 == 0 and min_diff != large) {
        answer -= min_diff;
    }
    if (answer % 42 == 0) {
        answer = 0;
    }
    cout << answer << endl;
}

Общая схема

Нетрудно заметить, что вторая и третья разновидности задачи являются частными случаями первой с нулевой длиной очереди.

Собственно, общая схема решения предполагает следующие входные данные:

• начальное состояние и начальный ответ
• величину задержки (длину очереди)
• функцию обновления состояния по элементу, удалённому из очереди
• функцию обновления ответа по состоянию и элементу, пришедшему в очередь

Эту схему можно записать (на хаскеллоподобном языке) так:

solve :: s -> a -> Int -> (x -> s -> s) -> (x -> s -> a -> a) -> [x] -> a
solve s0 a0 n updateS updateA xs = go s0 a0 $ zip xs (drop n xs)
  where
    go _ a [] = a
    go s a ((old,new):rest) = 
        let s' = (updateS old s)
        in go s' (updateA new s' a) rest

По-научному такая схема (а точнее говоря, некоторое её обобщение) называется зигохистоморфизмом.

Казалось бы, причём здесь зигохистоморфные препроморфизмы?..