Эксперимент 003: возможные количества чердачных сфер симплекса

(Проект эксперимента доступен здесь.)

11-16 июля 2025 года состоялась V конференция математических центров России. На ней было выступление Натальи Даурцевой, на котором была рассказана задача о чердачных сферах. Задача состоит в том, чтобы сказать, сколько вписанных в продолжения фасет какого-либо симплекса может быть. В этом эксперименте мы попробуем программно найти данные значения.

Теоретическое введение

Пусть для $n > 2$ дан какой-то $(n-1)$-мерный симплекс $\Delta$ в Евклидовом пространстве. (Формально $\Delta$ — выпуклая оболочка $n$ аффинно независимых точек.) У него есть ровно $n$ вершин и ровно $n$ фасет. Если продолжить эти фасеты до гиперплоскостей в $n-1$-мерном Евклидовом пространстве, в которое мы вложили $\Delta$, то получится разбиение этого пространства на части, единственная конечная часть которого — сам симплекс $\Delta$. Некоторые из этих частей (не все) имеют в качестве границ все $n$ гиперплоскостей. Будем называть их "чердаками".

Как мы знаем, в любой симплекс можно вписать сферу. Похожим образом в его чердаки можно также вписать сферы. Но не всегда. Например, в правильный тетраэдр вписываются ровно 5 сферы: одна внутри и ещё по одной против каждой вершины, за противоположной гранью. Несмотря на то, что в общем случае в тетраэдр можно вписать ужен 8 сфер: ещё 4 вписываются в чердаки, касаясь двух граней и продолжений двух других граней. Поэтому можно задать вопрос: а какие количества вписанных в чердаки сфер симплекса могут быть? Для $n = 3$ ответ очевиден: 4 (всегда можно вписать 1 окружность и вневписать 3 окружности). Для $n = 4$ ответ принимает все значения от $5$ до $8$. А что дальше?

Оказывается, для тетраэдров выполнены следующие утверждения:

Лемма 1. Пусть даны $\delta_1, \dots, \delta_n \in \{1; -1\}$. Обозначим за $S_i$ — объём фасеты против $i$-ой вершины симплекса. Тогда сфера касающаяся продолжения граней симплекса так, что для каждого $i \in \{1; \dots; n\}$ если она находится с той же стороны, что и $\Delta$, то $\delta_i = 1$, иначе $\delta_i = -1$, есть тогда и только тогда, когда $\sum_{i=1}^n S_i \delta_i > 0$.

Лемма 2. Симплекс с объёмами фасет $S_i$ существует тогда и только тогда, когда $2 S_i < \sum_{k=1}^n S_k$ для всех i. (Т.е. на данном набор чисел выполняются строгие "обобщённые неравенства треугольника".)

Первое доказывается очень просто. Пусть $A_i$ — $i$-ая вершина симплекса, а $C$ — какая-то точка. Тогда $V(A_1A_2\dots A_n) = \sum_{i=1}^n V(A_1\dots A_{i-1}CA_{i+1}\dots A_n)$ при правильной ориентации объёмов. Это понятно из того, что правая сторона является линейной функцией от $C$, которая равна левой константе во всех вершинах симплекса, а значит и во всём аффинном замыкании симплекса. Возьмём в качестве $C$ центр искомой сферы (если таковая имеется). Пусть также $r$ — радиус этой сферы. Тогда $V(A_1\dots A_{i-1}CA_{i+1}\dots A_n) = r S(A_1\dots A_{i-1}A_{i+1}\dots A_n) / n$. Откуда имеем

$$V(A_1A_2\dots A_n) = \frac{r}{n} \sum_{i=1}^n S(A_1\dots A_{i-1}A_{i+1}\dots A_n)$$ $$\frac{n V(A_1A_2\dots A_n)}{\sum_{i=1}^n S(A_1\dots A_{i-1}A_{i+1}\dots A_n)} = r$$

Т.е. радиус сферы определён тогда и только тогда, когда $\sum_{i=1}^n S(A_1\dots A_{i-1}A_{i+1}\dots A_n) > 0$. При этом $S(A_1\dots A_{i-1}A_{i+1}\dots A_n) = S_i \delta_i$, откуда получается искомое неравенство. Осталось понять верность в обратную сторону. Но тут достаточно действовать в обратную сторону: уравнения гиперплоскостей $(n-1)$-ой из граней линейно независимы (так как пересекаются ровно в 1 точке). Значит если их сдвинуть на $r := \frac{n V}{\sum_{i=1}^n S_i \delta_i}$ в правильные стороны, то получим центр сферы с радиусом $r$, касающуюся продолжения этих $n-1$ граней с правильных сторон. Тогда по равенству на объёмы получим, что она касается ещё и продолжения последней $n$-ой гиперплоскости с правильной стороны.

Второй факт авторы доказывать не умеют и ссылку дать тоже не могут. Извините.

Отсюда задача сводится к чисто алгебраической:

Для каждого набора $(S_1, \dots, S_n)$ положительных чисел рассмотрим количество неравенств вида $\sum_{i=1}^n \delta_i S_i > 0$, где $\delta_i \in \{1; -1\}$, которое для него выполняется. Найдите все возможные получаемые количества.

Но давайте пойдём с обратной стороны. Заметим, что все неравенства бьются на пары "противоположных": для коэффициентов неравенства $(\delta_i)_{i=1}^n$ можно рассмотреть коэффициенты $(-\delta_i)_{i=1}^n$. Если не выполнено равенство $\sum_{i=1}^n \delta_i S_i = 0$, то выполнено ровно одно из этих двух неравенств. Если же выполнено — то никакое. Таким образом можно задать вопрос, а сколько уравнений вида $\sum_{i=1}^n \delta_i S_i = 0$ (с точностью до смены знака всех слагаемых) у нас выполняется? Вычтя это число из $2^{n-1}$ мы получим один из возможных ответов.

Т.е. задача обрела следующий совсем алгебраический вид:

В пространстве $\mathbb{R}^n$ в области ограниченной неравенствами $0 < 2 x_i < \sum_{i=1}^n x_i$ рассмотрим уравнения вида $\sum_{i=1}^n \pm x_i = 0$ (отождествляя уравнения со сменой знаков). Для каждой точки области посмотрим какому числу уравнений оно удовлетворяет. Какие значения можно получить?

План компьютерного решения задачи

Несложно заметить, что все уравнения и условия-неравенства положительно однородны. Поэтому можно рассматривать только наборы, где $\sum_{i=1}^n x_i = 1$. В данном гиперпространстве область будет выглядеть как внутренность некоторого (конечного) политопа. Тогда план действий будет примерно следующий:

1. Получить список всевозможных систем уравнений вида $\sum_{i=1}^n \pm x_i = 0$, взяв только те, что имеют решения, и объединив все системы с одним и тем же множеством решений.
2. Построить политоп, полученный системой неравенств $0 < 2 x_i < \sum_{i=1}^n x_i$ в гиперплоскости $\sum_{i=1}^n x_i = 1$.
3. Для каждой системы из п. 1 определить аффинное подпространство её решений в гиперплоскости $\sum_{i=1}^n x_i = 1$ и проверить его пересечение с политопом, построенным в п. 2.
4. Вычислить все размеры систем уравнений, прошедших проверку в п.2 и вывести все полученные значения (по одному разу).

🔵 Отступление: приведение набора векторов к линейно независимому множеству определённого вида

Для дальнейших продвижений нам потребуются два вида преобразования множества векторов в ЛНМ.

Алгоритм Гаусса

Первое преобразование, которое нам понадобится — преобразование к приведённому ступенчатому виду. Но сначала определим "приведённый ступенчатый вид":

Определение. Говорят, что матрица из $m$ строк (или упорядоченный набор из $m$ векторов) имеет приведённый ступенчатый вид в определённом базисе, если для каждой $i$-ой строки матрицы (каждого $i$-ого вектора из набора) есть число $k_i$, что

$$1 \leqslant k_1 < k_2 < \dots < k_m \leqslant n$$

а также:

1. всякий коэффициент $i$-ой строки до $k_i$-ого исключительно равен $0$,
2. $k_i$-ый коэффициент $i$-ой строки равен $1$,
3. для всяких $1 \leqslant i < j \leqslant m$ $i$-ая строка имеет ноль на $k_j$-ом месте.

Иначе говоря, можно выбрать такие $k_i$, что выделяя столбцы с номерами $k_i$ получается единичная квадратная матрица размера $m$.

Для приведения данного набора векторов к приведённому ступенчатому виду есть всем известный алгоритм Гаусса.

Пусть фиксирован базис $B$. Шаг алгоритма, который расширяет ЛНМ $U$ размера $m$ в приведённом ступенчатом виде новым вектором $v$, работает следующим образом:

1. Пусть $\mathrm{firstNonZeroIndex}$ — первый индекс ненулевого коэффициента $v$ при разложении по $B$, $\mathrm{basisIndex} := 0$, $\mathrm{coefficientIndex} := 0$.
2. Пока $\mathrm{basisIndex} < m$ будем делать следующее:
   1. Если $\mathrm{firstNonZeroIndex} = n$, то добавим в новое ЛНМ все вектора старого с номера $\mathrm{basisIndex}$ включительно по номер $m$ исключительно, вернём новое ЛНМ и закончим работу алгоритма.
   2. Если $\mathrm{basisIndex} = m$, то просто разделим $v$ на его же первый ненулевой коэффициент $v[\mathrm{firstNonZeroIndex}]$, вычтем (новый) $v$ из всех добавленных элементов $u_i$ в новый ЛНМ с коэффициентом $u_i[\mathrm{firstNonZeroIndex}]$, добавим в новое ЛНМ получившийся $v$, вернём новое ЛНМ и закончим работу алгоритма.
   3. Обновим $\mathrm{coefficientIndex}$ до нового индекса первого ненулевого коэффициента $u_{\mathrm{basisIndex}}$.
   4. Если $\mathrm{coefficientIndex} > \mathrm{firstNonZeroIndex}$, то просто разделим $v$ на его первый ненулевой коэффициент ($v[\mathrm{firstNonZeroIndex}]$), каждый вектор $u_i$ для $i \in [\mathrm{basisIndex}; m)$ с индексом $k_i$ первого ненулевого коэффициента вычтем из $v$ с коэффициентом $v[k_i]$, затем вычтем $v$ из каждого $u_i$ для $i \in [0; \mathrm{basisIndex})$, после чего, наконец, добавим $v$, а за ним все оставшиеся $u_i$.
   5. Иначе добавим $u_{\mathrm{basisIndex}}$ в новое ЛНМ, вычтем его из $v$ с коэффициентом $v[\mathrm{coefficientIndex}]$, обновим $\mathrm{firstNonZeroIndex}$, инкрементируем $\mathrm{basisIndex}$ и $\mathrm{coefficientIndex}$ и продолжим цикл.

Чтобы привести всё ЛНМ к данному виду, достаточно добавлять в новое ЛНМ вектора старого по очереди с помощью данного алгоритма.

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

Лемма 3. Пусть $V_1$ и $V_2$ — два конечных множества векторов с одинаковой линейной оболочкой. Тогда два ЛНМ, получаемых при применении алгоритма Гаусса к $V_1$ и $V_2$ совпадают как упорядоченные списки векторов.

Этот факт на деле прост. Сначала докажем утверждение для случая, когда $V_1$ и $V_2$ — ЛНМ; общий факт тогда будет следовать из того, что мы выберем максимальные линейной независимые подмножества $U_1$ и $U_2$ множеств $V_1$ и $V_2$, что нет линейно независимых подмножеств с меньшими индексами элементов, отчего будет равенство линейных оболочек $\mathrm{span}(U_1) = \mathrm{span}(V_1) = \mathrm{span}(V_2) = \mathrm{span}(U_2)$, откуда будет равенство результатов преобразования $\mathrm{Gauss}(V_1) = \mathrm{Gauss}(U_1) = \mathrm{Gauss}(U_2) = \mathrm{Gauss}(V_2)$.

Тогда заметим, что приведение $V_i$ к $\mathrm{Gauss}(V_i)$ является обычным домножением на квадратную обратимую матрицу $G_i$ размера $m$ слева. При этом само по себе $V_2$ получается из $V_1$ домножением слева на квадратную обратимую матрицу $H$ размера $m$ слева. ($H$ по сути описывает смену базиса в $\mathrm{span}(V_1) = \mathrm{span}(V_2)$.) Таким образом

$$\mathrm{Gauss}(V_2) = G_2 V_2 = G_2 H V_1 = G_2 H G_1^{-1} \mathrm{Gauss}(V_1)$$

Покажем, что $M := G_2 H G_1^{-1}$ является единичной матрицей. Для этого заметим, что столбцы "ступенек" $\mathrm{Gauss}(V_1)$ — это набор из $m$ первых линейно независимых столбцов в матрице $\mathrm{Gauss}(V_1)$. А значит, что они линейно независимы и в $\mathrm{Gauss}(V_2)$, так как $M$ обратима. Аналогично наоборот: первые $m$ линейно независимых столбцов в $\mathrm{Gauss}(V_2)$ линейно независимы и в $\mathrm{Gauss}(V_1)$. Значит столбцы ступенек в $\mathrm{Gauss}(V_1)$ и в $\mathrm{Gauss}(V_1)$ совпадают и образуют в обеих матрицах единичные подматрицы. Следовательно, $M$ является единичной, откуда и сами $\mathrm{Gauss}(V_i)$ равны.

Программная реализация

Нам потребуется постепенно расширять имеющиеся ЛНМ в приведённом ступенчатом виде до новых ЛНМ. Для этого заведём простенький интерфейс:

GaussReducedRowEchelonFormBasisExtender.kt

interface BasisExtender<Vector> {
    fun KoneList<Vector>.extend(newVector: Vector): KoneList<Vector>
}

context(basisExtender: BasisExtender<Vector>)
fun <Vector> KoneList<Vector>.extend(newVector: Vector): KoneList<Vector> =
    with(basisExtender) { this@extend.extend(newVector) }

Сам интерфейс абстрагирует использование алгоритма базиса: он просто даёт фасад, для использования алгоритма, не смотря на конкретную реализацию. Функция после него позволяет использовать метод extend интерфейса BasisExtender, когда последний доступен в качестве контекстного параметра.

Сама реализация будет выглядеть так:

GaussReducedRowEchelonFormBasisExtender.kt

class GaussDiagonalBasisExtender<Number, Vector>(
    private val field: Field<Number>,
    private val vectorSpace: VectorSpace<Number, Vector>,
    private val dimension: UInt,
    private val basisDecomposition: VectorSpaceBasisDecomposition<Number, Vector, UInt>,
) : BasisExtender<Vector> {
    override fun KoneList<Vector>.extend(newVector: Vector): KoneList<Vector> = context(field, vectorSpace) {
        KoneList.build(this.size + 1u) {
            var newVector = newVector
            var firstNonZeroIndex = 0u
            while (firstNonZeroIndex < dimension && basisDecomposition.decompose(newVector)[firstNonZeroIndex].isZero()) firstNonZeroIndex++

            var basisIndex = 0u
            var coefficientIndex = 0u
            while (true) {
                if (firstNonZeroIndex == dimension) {
                    addSeveral(this@extend.size - basisIndex) { this@extend[it + basisIndex] }
                    break
                }
                if (basisIndex == this@extend.size) {
                    newVector /= basisDecomposition.decompose(newVector)[firstNonZeroIndex]
                    for (index in 0u ..< basisIndex) {
                        this[index] -= newVector * basisDecomposition.decompose(this[index])[firstNonZeroIndex]
                    }
                    add(newVector)
                    break
                }

                val basisVectorDecomposition = basisDecomposition.decompose(this@extend[basisIndex])
                while (basisVectorDecomposition[coefficientIndex].isZero()) coefficientIndex++

                if (coefficientIndex > firstNonZeroIndex) {
                    newVector /= basisDecomposition.decompose(newVector)[firstNonZeroIndex]
                    val basisSeparator = basisIndex
                    while (basisIndex < this@extend.size) {
                        val basisVectorDecomposition = basisDecomposition.decompose(this@extend[basisIndex])
                        while (basisVectorDecomposition[coefficientIndex].isZero()) coefficientIndex++
                        newVector -= this@extend[basisIndex] * basisDecomposition.decompose(newVector)[coefficientIndex]
                        basisIndex++
                        coefficientIndex++
                    }
                    for (index in 0u ..< basisSeparator) {
                        this[index] -= newVector * basisDecomposition.decompose(this[index])[firstNonZeroIndex]
                    }
                    add(newVector)
                    addSeveral(this@extend.size - basisSeparator) { this@extend[it + basisSeparator] }
                    break
                }

                val newVectorCoef = basisDecomposition.decompose(newVector)[coefficientIndex]
                add(this@extend[basisIndex])
                newVector -= this@extend[basisIndex] * newVectorCoef
                while (firstNonZeroIndex < dimension && basisDecomposition.decompose(newVector)[firstNonZeroIndex].isZero()) firstNonZeroIndex++
                basisIndex++
                coefficientIndex++
            }
        }
    }
}

Давайте разберём по частям. Сначала идёт мишура формальная обёртка для алгоритма:

GaussReducedRowEchelonFormBasisExtender.kt

class GaussDiagonalBasisExtender<Number, Vector>(
    private val field: Field<Number>,
    private val vectorSpace: VectorSpace<Number, Vector>,
    private val dimension: UInt,
    private val basisDecomposition: VectorSpaceBasisDecomposition<Number, Vector, UInt>,
) : BasisExtender<Vector> {
    override fun KoneList<Vector>.extend(newVector: Vector): KoneList<Vector> = context(field, vectorSpace) {
        KoneList.build(this.size + 1u) {
            // ...
        }
    }
}

Тут мы написали, что мы рассматриваем алгоритм над числами типа Number и векторами типа Vector, принимающий в качестве перманентных аргументов:

1. field: Field<Number> — объект, который реализует базовые операции на числах (значение нуля и единицы, операции сложения, вычитания и умножения и т.д.),
2. vectorSpace: VectorSpace<Number, Vector> — объект, который реализует базовые операции на векторах и числах (значение нулевого вектора, сложение и вычитание векторов, умножения на скаляр и т.д.),
3. dimension: UInt — размерность рассматриваемого пространства и рассматриваемого базиса basisDecomposition,
4. basisDecomposition: VectorSpaceBasisDecomposition<Number, Vector, UInt> — объект, который для всякого индекса i типа UInt (на самом деле не любого, а от 0 до размерности рассматриваемого в будущем пространства исключительно) определяет коэффициент типа Number при разложении данного вектора типа Vector по данному фиксированному базису (в будущем этим базисом будет стандартный базис из векторов вида $(0, \dots, 0, 1, 0, \dots, 0)$).

Перегрузка функции определяет то, каким образом будет реализован алгоритм. Вызов context(field, vectorSpace) { /* ... */ } просто делает field и vectorSpace доступными в качестве контекстных параметров. Вызов KoneList.build(this.size + 1u) { /* ... */ } внутри него помогает построить новый список-ЛНМ. Теперь перейдём к самому алгоритму построения.

Сначала идут строки

GaussReducedRowEchelonFormBasisExtender.kt

var newVector = newVector
var firstNonZeroIndex = 0u
while (firstNonZeroIndex < dimension && basisDecomposition.decompose(newVector)[firstNonZeroIndex].isZero()) firstNonZeroIndex++

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

GaussReducedRowEchelonFormBasisExtender.kt

var basisIndex = 0u
var coefficientIndex = 0u
while (true) {
    /* ... */
}

На начало каждой итерации цикла while basisIndex указывает на индекс первого ещё не рассмотренного вектора из нашего ЛНМ, coefficientIndex — на индекс, не больший индекса первого ненулевого коэффициента basisIndex-ого элемента имеющегося ЛНМ (если только элемент с индексом basisIndex существует; иначе coefficientIndex не имеет осмысленного значения), а firstNonZeroIndex — на индекс первого ненулевого коэффициента newVector (или dimension, если newVector нулевой). Также на начало итерации сохраняем инвариант, что мы добавили только векторы, первые ненулевые коэффициенты которых идут строго до первого ненулевого коэффициента newVector.

В начале итерации мы проверяем банальные краевые случаи:

1. если coefficientIndex == dimension, т.е. newVector — нулевой вектор,
2. если basisIndex == this@extend.size, т.е. мы прошлись по всем векторам (и вычли их из newVector с правильными коэффициентами),
3. если (после обновления coefficientIndex до первого индекса ненулевого коэффициента basisIndex-ого вектора в ЛНМ) coefficientIndex > firstNonZeroIndex, т.е. первый ненулевой коэффициент у newVector идёт строго раньше, чем у basisIndex-ого вектора.

Т.е. код имеет вид

GaussReducedRowEchelonFormBasisExtender.kt

while (true) {
    if (firstNonZeroIndex == dimension) {
        /* ... */
    }
    if (basisIndex == this@extend.size) {
        /* ... */
    }

    val basisVectorDecomposition = basisDecomposition.decompose(this@extend[basisIndex])
    while (basisVectorDecomposition[coefficientIndex].isZero()) coefficientIndex++

    if (coefficientIndex > firstNonZeroIndex) {
        /* ... */
    }

    /* ... */
}

В первом случае, достаточно просто добавить вектора, которые мы ещё не добавили и закончить работу:

GaussReducedRowEchelonFormBasisExtender.kt

if (firstNonZeroIndex == dimension) {
    addSeveral(this@extend.size - basisIndex) { this@extend[it + basisIndex] }
    break
}

Во втором случае мы получаем ненулевой вектор newVector. При этом все остальные векторы добавлены в новый список-ЛНМ. Значит, если мы добавим newVector в этот список, то мы получим новое ЛНМ в ступенчатом виде. Чтобы получить приведённый ступенчатый вид, надо отнормировать вектор по его первому ненулевому коэффициенту (сделать его единичным)4 и вычесть из всех предыдущих векторов с правильными коэффициентами:

GaussReducedRowEchelonFormBasisExtender.kt

if (basisIndex == this@extend.size) {
    newVector /= basisDecomposition.decompose(newVector)[firstNonZeroIndex]
    for (index in 0u ..< basisIndex) {
        this[index] -= newVector * basisDecomposition.decompose(this[index])[firstNonZeroIndex]
    }
    add(newVector)
    break
}

В третьем случае, у нас дан ненулевой newVector, basisIndex указывает на какой-то вектор, но первый ненулевой индекс basisIndex-ого вектора идёт строго после первого ненулевого индекса newVector. Значит, если добавить сначала newVector, а потом все остальные элементы старого ЛНМ, то получим ЛНМ в ступенчатом виде. Для приведённого ступенчатого вида нужно сначала вычесть все не добавленные вектора из newVector с правильными коэффициентами, а потом вычесть newVector с правильными коэффициентами из уже добавленных векторов:

GaussReducedRowEchelonFormBasisExtender.kt

if (coefficientIndex > firstNonZeroIndex) {
    newVector /= basisDecomposition.decompose(newVector)[firstNonZeroIndex]
    val basisSeparator = basisIndex
    while (basisIndex < this@extend.size) {
        val basisVectorDecomposition = basisDecomposition.decompose(this@extend[basisIndex])
        while (basisVectorDecomposition[coefficientIndex].isZero()) coefficientIndex++
        newVector -= this@extend[basisIndex] * basisDecomposition.decompose(newVector)[coefficientIndex]
        basisIndex++
        coefficientIndex++
    }
    for (index in 0u ..< basisSeparator) {
        this[index] -= newVector * basisDecomposition.decompose(this[index])[firstNonZeroIndex]
    }
    add(newVector)
    addSeveral(this@extend.size - basisSeparator) { this@extend[it + basisSeparator] }
    break
}

Наконец, если и третий краевой случай не выполняется, то первый ненулевой коэффициент newVector идёт не раньше, чем в basisIndex-ом векторе. Значит, второй нужно просто вычесть из первого с правильным коэффициентом, чтобы первый ненулевой коэффициент шёл строго после, а затем добавить basisIndex-ый вектор в наш ЛНМ, обновляя basisIndex, firstNonZeroIndex (и на всякий случай coefficientIndex):

GaussReducedRowEchelonFormBasisExtender.kt

val newVectorCoef = basisDecomposition.decompose(newVector)[coefficientIndex]
add(this@extend[basisIndex])
newVector -= this@extend[basisIndex] * newVectorCoef
while (firstNonZeroIndex < dimension && basisDecomposition.decompose(newVector)[firstNonZeroIndex].isZero()) firstNonZeroIndex++
basisIndex++
coefficientIndex++

Алгоритм ортогонализации Грама-Шмидта

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

Конкретно, алгоритм Грама-Шмидта заменяет набор векторов $(v_i)_{i=1}^m$ на набор векторов $(u_i)_{i=1}^m$, что $\mathrm{span}(v_1, \dots, v_i) = \mathrm{span}(u_1, \dots, u_i)$, а любые два вектора $u_i$ друг другу ортогональны (в том смысле, что нулевой вектор ортогонален любому другому). Опишем сам алгоритм.

Шаг алгоритма, который добавляет в ортогональный набор векторов $U$ размера $m$ (с выполнением условия, что $\mathrm{span}(u_1, \dots, u_i) = \mathrm{span}(v_1, \dots, v_i)$) новый вектор $v_{m+1}$, работает следующим образом:

1. Для каждого $i \in (1; \dots; m)$, что $u_i \neq \overline{0}$ просто заменим $v_{m+1}$ на $$v_{m+1} := v_{m+1} - u_i \frac{\langle v_{m+1}, u_i \rangle}{\langle u_i, u_i \rangle}$$

Т.е. фактически мы заменяем $v_{m+1}$ на

$$u_{m+1} := v_{m+1} - \sum_{\substack{i=1..m\\ u_i \neq \overline{0}}} u_i \frac{\langle v_{m+1}, u_i \rangle}{\langle u_i, u_i \rangle}$$

(Так как все $u_i$ ортогональны, то не важно, будем ли мы вычислять вычитаемые слагаемые одновременно и вычитать их одновременно, как в формуле, или будем вычислять вычитать их итеративно, как в алгоритме.)

Чтобы привести всё конечное множество векторов в ортогональный вид, достаточно в пустое множество добавлять эти вектора с помощью описанного шага. Чтобы получить ЛНМ, достаточно отбросить получаемые нулевые вектора.

Программная реализация

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

Реализация будет выглядеть так:

orthogonalizeByGramSchmidt.kt

context(_: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, *>)
fun <Number, Vector> KoneList<Vector>.orthogonalizeByGramSchmidt(): KoneList<Vector> =
    KoneList.build {
        val basis = this@orthogonalizeByGramSchmidt
        for (index in 0u ..< basis.size)
            +(0u ..< index).fold(basis[index]) { accumulator, subindex ->
                val subVector = this[subindex]
                accumulator - subVector * ((accumulator dot subVector) / (subVector dot subVector))
            }
    }

Давайте разберём по частям. Сначала идёт определение функции:

orthogonalizeByGramSchmidt.kt

context(_: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, *>)
fun <Number, Vector> KoneList<Vector>.orthogonalizeByGramSchmidt(): KoneList<Vector> =

Тут описана функция orthogonalizeByGramSchmidt с типовыми параметрами Number и Vector для чисел и векторов, с которыми мы будем работать и следующими аргументами:

1. Контекстный аргумент _: Field<Number> — объект, который реализует базовые операции на числах в поле.
2. Контекстный аргумент _: EuclideanSpaceOverField<Number, Vector, *> — объект, который реализует базовые операции на векторах и числах в евклидовом векторном поле.
3. Ресивер KoneList<Vector>. — данный базис для ортогонализации.

После = в конце этого сниппета идёт возвращаемое тело функции. А именно дальше идёт простое выражение, которое строит новый список по описанию в фигурных скобках:

orthogonalizeByGramSchmidt.kt

KoneList.build {
    /* ... */
}

Первое, что мы делаем в этих скобках — просто присвоение ресивера this@orthogonalizeByGramSchmidt нашей функции в переменную basis:

orthogonalizeByGramSchmidt.kt

val basis = this@orthogonalizeByGramSchmidt

Дальше расположен цикл for, который проходиться переменной index по всем индексам базиса basis и вносит index-ый элемент basis-а в наш конструируемый список заранее преобразовывая его алгоритмом Грама-Шмидта:

orthogonalizeByGramSchmidt.kt

for (index in 0u ..< basis.size)
    /* ... */

Разберёмся с тем, как происходит итерация алгоритма Грама-Шмидта. Взяв новый элемент базиса basis[index], нам надо вычесть из него его "компоненту" по каждому из предыдущих уже ортогонализированных элементов базиса. Если нынешнее значение нашего вектора после вычитания первых subindex компонент есть accumulator, а следующий вычитаемый вектор — subVector, то компонента accumulator по этому базисному элементу будет равна

orthogonalizeByGramSchmidt.kt

subVector * ((accumulator dot subVector) / (subVector dot subVector))

где u dot v — скалярное произведение u и v. Тогда новое значение нашего вектора будет равно

orthogonalizeByGramSchmidt.kt

accumulator - subVector * ((accumulator dot subVector) / (subVector dot subVector))

Чтобы итеративно сделать такие изменения, используем функцию fold, итерирующимся по всем предыдущим индексам (т.е. subindex из интервала 0u ..< index). При этом в качестве начального значения даём basis[index]:

orthogonalizeByGramSchmidt.kt

(0u ..< index).fold(basis[index]) { accumulator, subindex ->
    val subVector = this[subindex]
    accumulator - subVector * ((accumulator dot subVector) / (subVector dot subVector))
}

Тут this указывает на конструируемый список. Отчего val subVector = this[subindex] указывает на subindex-ый уже ортогонализированный элемент. Поэтому в блоке при fold мы сначала помещаем этот вектор в переменную subVector, а потом используем его для подмены имеющегося значения accumalator. После применения функции fold будет возвращено последнее значение accumulator, т.е. искомый вектор. Чтобы добавить его в список просто припишем + перед получаемым значением.

Таким образом список заполняется ортогонализированным базисом.

Нахождение всевозможных систем-претендентов

Теперь займёмся получением всевозможных систем из уравнений вида $\sum_{i=1}^n \pm x_i = 0$. Самый простой способ делать это — это линейным программированием перебирать всевозможные уравнения и расширять все имеющиеся системы (не удаляя и/или изменяя изначальные, а лишь добавляя новые).

Заведём для этого функцию:

computePossibleEquationSystems.kt

data class EquationCase<Covector>(
    val initialCovectors: KoneList<Covector>,
    val basis: KoneList<Covector>,
)

fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
): KoneList<KoneList<EquationCase<Covector>>> {
    TODO()
}

Здесь Covector определяет тип, описывающий уравнение, т.е. ковектор. Получаем мы (пока) размерность задачи (число $n$), а выводим все разрешимые нетривиальным способом системы искомого вида, разбитые по размерности на вложенные списки. Т.е. внутренний список с индексом $i$ во внешнем описывает все системы размерности $i$. Понятно, что системы имеют размерность не более $n$, но при размерности $n-1$ система не имеет нетривиальных решений. Поэтому такие системы мы будем сразу отбрасывать, из-за чего внешний список имеет длину $n$ (по всем размерностям от $0$ до $n-1$). При этом EquationCase описывает множество изначальные ковектора в списке initialCovectors, а также получающийся из них ЛНМ basisв приведённом ступенчатом виде. Мы специально будем по дороге вычислять эти ЛНМ, чтобы сразу отсеивать неразрешимые нетривиальным образом системы. (И чтобы не вычислять их потом.)

Мы будем, перебирая все уравнения собирать системы, начиная со случая только одной пустой системы. Поэтому начнём нашу функцию очень просто:

computePossibleEquationSystems.kt

fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
): KoneList<KoneList<EquationCase<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneArrayGrowableList<EquationCase<Covector>>()
    }

    TODO()

    return equationCases
}

Начать перебор надо с присутствия только пустой системы:

computePossibleEquationSystems.kt

fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
): KoneList<KoneList<EquationCase<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneArrayGrowableList<EquationCase<Covector>>()
    }

    equationCases[0u].add(
        EquationCase(
            initialCovectors = KoneList.empty(),
            basis = KoneList.empty(),
        )
    )

    TODO()

    return equationCases
}

Теперь научимся перебирать все уравнения. Уравнения задаются коэффициентами при каждой переменной. Так как каждый коэффициент независимо может быть равен $1$ или $-1$, то нам достаточно возвести множество $\{1; -1\}$ в декарттову степень $n$. Для этого у меня есть функция cartesianPower. Но такие уравнения бьются на пары эквивалентных, в которых одно получается из другого заменой знака всех коэффициентов. Чтобы избавиться от этой проблемы, скажем, что первый коэффициент всегда равен $1$. Таким образом у нас получается цикл

computePossibleEquationSystems.kt

fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
): KoneList<KoneList<EquationCase<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneArrayGrowableList<EquationCase<Covector>>()
    }

    equationCases[0u].add(
        EquationCase(
            initialCovectors = KoneList.empty(),
            basis = KoneList.empty(),
        )
    )

    for (coefficients in KoneList.of(1, -1).cartesianPower(dimension - 1u)) {
        TODO()
    }

    return equationCases
}

Заметим, что случай всех единиц нас тоже не устраивает, так как в таком случае решений сразу не будет в искомой области.

computePossibleEquationSystems.kt

context(_: VectorSpace<*, Covector>)
fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
    covectorBasis: VectorSpaceBasis.Finite<*, Covector>,
): KoneList<KoneList<EquationCase<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneArrayGrowableList<EquationCase<Covector>>()
    }

    equationCases[0u].add(
        EquationCase(
            initialCovectors = KoneList.empty(),
            basis = KoneList.empty(),
        )
    )

    for (coefficients in KoneList.of(1, -1).cartesianPower(dimension - 1u)) {
        if (coefficients.all { it == 1 }) continue

        TODO()
    }

    return equationCases
}

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

computePossibleEquationSystems.kt

context(_: VectorSpace<*, Covector>)
fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
    covectorBasis: VectorSpaceBasis.Finite<*, Covector>,
): KoneList<KoneList<EquationCase<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneArrayGrowableList<EquationCase<Covector>>()
    }

    equationCases[0u].add(
        EquationCase(
            initialCovectors = KoneList.empty(),
            basis = KoneList.empty(),
        )
    )

    for (coefficients in KoneList.of(1, -1).cartesianPower(dimension - 1u)) {
        if (coefficients.all { it == 1 }) continue
        val vector = coefficients.foldIndexed(covectorBasis[0u]) { index, accumulator, element ->
            accumulator + element * covectorBasis[index + 1u]
        }

        TODO()
    }

    return equationCases
}

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

computePossibleEquationSystems.kt

context(_: VectorSpace<*, Covector>)
fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
    covectorBasis: VectorSpaceBasis.Finite<*, Covector>,
): KoneList<KoneList<EquationCase<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneArrayGrowableList<EquationCase<Covector>>()
    }

    equationCases[0u].add(
        EquationCase(
            initialCovectors = KoneList.empty(),
            basis = KoneList.empty(),
        )
    )

    for (coefficients in KoneList.of(1, -1).cartesianPower(dimension - 1u)) {
        if (coefficients.all { it == 1 }) continue
        val vector = coefficients.foldIndexed(covectorBasis[0u]) { index, accumulator, element ->
            accumulator + element * covectorBasis[index + 1u]
        }

        for (dim in dimension - 1u downTo 0u) for (case in equationCases[dim].toKoneList()) {
            TODO()
        }
    }

    return equationCases
}

Подмечу, что копирование происходит в вызове .toKoneList().

Теперь нам нужно расширить ЛНМ в приведённом ступенчатом виде, который мы используем. Для этого нам потребуется аргумент типа BasisExtender<Covector>:

computePossibleEquationSystems.kt

context(_: VectorSpace<*, Covector>, _: BasisExtender<Covector>)
fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
    covectorBasis: VectorSpaceBasis.Finite<*, Covector>,
): KoneList<KoneList<EquationCase<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneArrayGrowableList<EquationCase<Covector>>()
    }

    equationCases[0u].add(
        EquationCase(
            initialCovectors = KoneList.empty(),
            basis = KoneList.empty(),
        )
    )

    for (coefficients in KoneList.of(1, -1).cartesianPower(dimension - 1u)) {
        if (coefficients.all { it == 1 }) continue
        val vector = coefficients.foldIndexed(covectorBasis[0u]) { index, accumulator, element ->
            accumulator + element * covectorBasis[index + 1u]
        }

        for (dim in dimension - 1u downTo 0u) for (case in equationCases[dim].toKoneList()) {
            val newBasis = case.basis.extend(vector)
            TODO()
        }
    }

    return equationCases
}

Получив новый базис, сформируем новый EquationCase и добавим в список с индексом newBasis.size. Но только если такой индекс существует. Получится такой итоговый код:

computePossibleEquationSystems.kt

context(_: VectorSpace<*, Covector>, _: BasisExtender<Covector>)
fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
    covectorBasis: VectorSpaceBasis.Finite<*, Covector>,
): KoneList<KoneList<EquationCase<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneArrayGrowableList<EquationCase<Covector>>()
    }

    equationCases[0u].add(
        EquationCase(
            initialCovectors = KoneList.empty(),
            basis = KoneList.empty(),
        )
    )

    for (coefficients in KoneList.of(1, -1).cartesianPower(dimension - 1u)) {
        if (coefficients.all { it == 1 }) continue
        val vector = coefficients.foldIndexed(covectorBasis[0u]) { index, accumulator, element ->
            accumulator + element * covectorBasis[index + 1u]
        }

        for (dim in dimension - 1u downTo 0u) for (case in equationCases[dim].toKoneList()) {
            val newBasis = case.basis.extend(vector)
            if (newBasis.size == dimension) continue
            equationCases[newBasis.size].add(
                EquationCase(
                    initialCovectors = KoneList.build(case.initialCovectors.size + 1u) {
                        addAllFrom(case.initialCovectors)
                        add(vector)
                    },
                    basis = newBasis,
                )
            )
        }
    }

    return equationCases
}

Построение политопа возможных объёмов фасет симплекса

По второму пункту нашего плана нужно в гиперплоскости $\sum_{i=1}^n x_i = 1$ найти все реализуемые наборы объёмов фасет симплекса. Как мы помним, это внутренность выпуклого политопа, получающегося при проведении полупространств, задаваемых уравнениями вида $x_k > 0$ и вида $\sum_{i=1}^n x_i > 2 x_k$. Значит, нам нужно построить этот политоп. Под "построить" подразумевается "найти все его вершины как положения в аффинном пространстве и выписать все его грани, указав у каждой, какие другие грани является её гранями". Для этого воспользуемся следующим планом:

1. Сделаем политоп в аффинном гиперпространстве $\sum_{i=1}^n x_i = 1$, ограниченный неравенствами $x_k > 0$. Это всего лишь симплекс натянутый на вершины вида $(0, \dots, 0, 1, 0 \dots, 0)$.
2. Итеративно для каждого неравенства вида $\sum_{i=1}^n x_i > 2 x_k$ разрубим всю политопическую конструкцию гиперплоскостью $\sum_{i=1}^n x_i - 2 x_k = 0$ и уберём все политопы, находящиеся по отрицательную сторону от этой гиперплоскости.

🔵 Отступление: разбиение политопической конструкции по линейной функции на аффинном пространстве

divideByAndPut.kt

data class PolytopeAndSign<Polytope>(
    val polytope: Polytope,
    val sign: Sign,
)

@IgnorableReturnValue
context(_: Field<Number>, _: Order<Number>, _: AffineSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> AbstractPolytopicConstruction<Point>.divideByAndPut(
    linearFunction: (Point) -> Number,
    target: AbstractPolytopicConstruction<Point>,
): KoneList<KoneMap<AbstractPolytopicConstruction.Polytope<Point>, KoneList<PolytopeAndSign<AbstractPolytopicConstruction.Polytope<Point>>>>> {
    val vertexMapping = vertices.associateWith(
        keyEquality = Equality.absoluteFor(),
        keyHashing = Hashing.defaultFor(),
        keyOrder = null
    ) { PolytopeAndSign(target.addVertex(it.position), linearFunction(it.position).sign()) }
    val polytopeMapping = KoneArrayFixedCapacityList<KoneMap<AbstractPolytopicConstruction.Polytope<Point>, KoneList<PolytopeAndSign<AbstractPolytopicConstruction.Polytope<Point>>>>>(spaceDimension + 1u)
    polytopeMapping.add(
        vertexMapping.nodesView.associate(
            keyEquality = Equality.absoluteFor(),
            keyHashing = Hashing.defaultFor(),
            keyOrder = null,
        ) { it.key.asPolytope() mapsTo KoneList.of(PolytopeAndSign(it.value.polytope.asPolytope(), it.value.sign)) }
    )

    val edges = polytopesOfDimension(1u)
    val polytopeSectionMapping = KoneArrayFixedCapacityList<KoneMap<AbstractPolytopicConstruction.Polytope<Point>, Maybe<AbstractPolytopicConstruction.Polytope<Point>>>>(spaceDimension)
    polytopeSectionMapping.add(
        edges.associateWith(
            keyEquality = Equality.absoluteFor(),
            keyHashing = Hashing.defaultFor(),
            keyOrder = null,
        ) {
            val (end1, end2) = it.vertices.toKoneList()
            val end1Sign = polytopeMapping[0u][end1.asPolytope()].single().sign
            val end2Sign = polytopeMapping[0u][end2.asPolytope()].single().sign
            if (end1Sign == end2Sign || end1Sign.isZero() || end2Sign.isZero()) None
            else Some(
                target.addVertex(
                    end2.position + (end1.position - end2.position) / (linearFunction(end1.position) - linearFunction(end2.position)) * (-linearFunction(end2.position)),
                ).asPolytope()
            )
        }
    )
    polytopeMapping.add(
        edges.associateWith(
            keyEquality = Equality.absoluteFor(),
            keyHashing = Hashing.defaultFor(),
            keyOrder = null,
        ) {
            val (end1, end2) = it.vertices.toKoneList()
            when (val section = polytopeSectionMapping[0u][it]) {
                None -> KoneList.of(
                    PolytopeAndSign(
                        polytope = target.addPolytope(
                            dimension = 1u,
                            vertices = KoneReifiedSet.of(
                                vertexMapping[end1].polytope, vertexMapping[end2].polytope,
                                elementReification = Reification.defaultFor(),
                                elementEquality = Equality.absoluteFor(),
                                elementHashing = Hashing.defaultFor(),
                                elementOrder = null,
                            ),
                            faces = KoneList.of(
                                KoneReifiedSet.of(
                                    vertexMapping[end1].polytope.asPolytope(), vertexMapping[end2].polytope.asPolytope(),
                                    elementReification = Reification.defaultFor(),
                                    elementEquality = Equality.absoluteFor(),
                                    elementHashing = Hashing.defaultFor(),
                                    elementOrder = null,
                                ),
                            ),
                        ),
                        sign = vertexMapping[end1].sign.takeIf { it.isNonZero() } ?: vertexMapping[end2].sign,
                    ),
                )
                is Some<AbstractPolytopicConstruction.Polytope<Point>> -> KoneList.of(
                    PolytopeAndSign(
                        polytope = target.addPolytope(
                            dimension = 1u,
                            vertices = KoneReifiedSet.of(
                                vertexMapping[end1].polytope, section.value.vertices.single(),
                                elementReification = Reification.defaultFor(),
                                elementEquality = Equality.absoluteFor(),
                                elementHashing = Hashing.defaultFor(),
                                elementOrder = null,
                            ),
                            faces = KoneList.of(
                                KoneReifiedSet.of(
                                    vertexMapping[end1].polytope.asPolytope(), section.value,
                                    elementReification = Reification.defaultFor(),
                                    elementEquality = Equality.absoluteFor(),
                                    elementHashing = Hashing.defaultFor(),
                                    elementOrder = null,
                                ),
                            ),
                        ),
                        sign = vertexMapping[end1].sign
                    ),
                    PolytopeAndSign(
                        polytope = target.addPolytope(
                            dimension = 1u,
                            vertices = KoneReifiedSet.of(
                                vertexMapping[end2].polytope, section.value.vertices.single(),
                                elementReification = Reification.defaultFor(),
                                elementEquality = Equality.absoluteFor(),
                                elementHashing = Hashing.defaultFor(),
                                elementOrder = null,
                            ),
                            faces = KoneList.of(
                                KoneReifiedSet.of(
                                    vertexMapping[end2].polytope.asPolytope(), section.value,
                                    elementReification = Reification.defaultFor(),
                                    elementEquality = Equality.absoluteFor(),
                                    elementHashing = Hashing.defaultFor(),
                                    elementOrder = null,
                                ),
                            ),
                        ),
                        sign = vertexMapping[end2].sign
                    ),
                )
            }
        }
    )

    for (dim in 2u .. spaceDimension) {
        val polytopes = polytopesOfDimension(dim)
        polytopeSectionMapping.add(
            polytopes.associateWith(
                keyEquality = Equality.absoluteFor(),
                keyHashing = Hashing.defaultFor(),
                keyOrder = null,
            ) { polytope ->
                val signs = polytope.vertices.map { vertexMapping[it].sign }

                if (signs.all { it.isNonNegative() } || signs.all { it.isNonPositive() }) return@associateWith None

                Some(
                    target.addPolytope(
                        dimension = dim - 1u,
                        vertices = KoneReifiedSet.build(
                            elementReification = Reification.defaultFor(),
                            elementEquality = Equality.absoluteFor(),
                            elementHashing = Hashing.defaultFor(),
                            elementOrder = null,
                        ) {
                            polytope.facesOfDimension(1u).mapSomeTo(this) {
                                polytopeSectionMapping[0u][it].map { it.vertices.single() }
                            }
                            polytope.vertices.map { vertexMapping[it] }.filter { it.sign.isZero() }.mapTo(this) { it.polytope }
                        },
                        faces = KoneList(dim - 1u) { subdim ->
                            KoneReifiedSet.build(
                                elementReification = Reification.defaultFor(),
                                elementEquality = Equality.absoluteFor(),
                                elementHashing = Hashing.defaultFor(),
                                elementOrder = null,
                            ) {
                                polytope.facesOfDimension(subdim + 1u).mapSomeTo(this) {
                                    polytopeSectionMapping[subdim][it]
                                }
                                polytope.facesOfDimension(subdim).flatMapTo(this) { polytopeMapping[subdim][it].filter { it.sign.isZero() }.map { it.polytope } }
                            }
                        },
                    )
                )
            }
        )
        polytopeMapping.add(
            polytopes.associateWith(
                keyEquality = Equality.absoluteFor(),
                keyHashing = Hashing.defaultFor(),
                keyOrder = null,
            ) { polytope ->
                when (val section = polytopeSectionMapping[dim - 1u][polytope]) {
                    None -> KoneList.of(
                        PolytopeAndSign(
                            polytope = target.addPolytope(
                                dimension = dim,
                                vertices = polytope.vertices.mapTo(
                                    KoneMutableReifiedSet.of(
                                        elementReification = Reification.defaultFor(),
                                        elementEquality = Equality.absoluteFor(),
                                        elementHashing = Hashing.defaultFor(),
                                        elementOrder = null,
                                    )
                                ) { vertexMapping[it].polytope },
                                faces = KoneList(dim) { subdim ->
                                    polytope.facesOfDimension(subdim).mapTo(
                                        KoneMutableReifiedSet.of(
                                            elementReification = Reification.defaultFor(),
                                            elementEquality = Equality.absoluteFor(),
                                            elementHashing = Hashing.defaultFor(),
                                            elementOrder = null,
                                        )
                                    ) { polytopeMapping[subdim][it].single().polytope }
                                }
                            ),
                            sign = polytope.vertices.firstOfThatOrNull({ vertexMapping[it].sign }) { it.isNonZero() } ?: Sign.Zero,
                        ),
                    )
                    is Some<AbstractPolytopicConstruction.Polytope<Point>> -> KoneList.of(
                        PolytopeAndSign(
                            polytope = target.addPolytope(
                                dimension = dim,
                                vertices = KoneReifiedSet.build(
                                    elementReification = Reification.defaultFor(),
                                    elementEquality = Equality.absoluteFor(),
                                    elementHashing = Hashing.defaultFor(),
                                    elementOrder = null,
                                ) {
                                    addAllFrom(section.value.vertices)
                                    polytope.vertices.map { vertexMapping[it] }.filter { it.sign.isNonNegative() }.mapTo(this) { it.polytope }
                                },
                                faces = KoneList(dim) { subdim ->
                                    KoneReifiedSet.build(
                                        elementReification = Reification.defaultFor(),
                                        elementEquality = Equality.absoluteFor(),
                                        elementHashing = Hashing.defaultFor(),
                                        elementOrder = null,
                                    ) {
                                        if (subdim < dim - 1u) addAllFrom(section.value.facesOfDimension(subdim))
                                        else add(section.value)
                                        polytope.facesOfDimension(subdim).flatMapTo(this) { polytopeMapping[subdim][it].filter { it.sign.isNonNegative() }.map { it.polytope } }
                                    }
                                },
                            ),
                            sign = Sign.Positive,
                        ),
                        PolytopeAndSign(
                            polytope = target.addPolytope(
                                dimension = dim,
                                vertices = KoneReifiedSet.build(
                                    elementReification = Reification.defaultFor(),
                                    elementEquality = Equality.absoluteFor(),
                                    elementHashing = Hashing.defaultFor(),
                                    elementOrder = null,
                                ) {
                                    addAllFrom(section.value.vertices)
                                    polytope.vertices.map { vertexMapping[it] }.filter { it.sign.isNonPositive() }.mapTo(this) { it.polytope }
                                },
                                faces = KoneList(dim) { subdim ->
                                    KoneReifiedSet.build(
                                        elementReification = Reification.defaultFor(),
                                        elementEquality = Equality.absoluteFor(),
                                        elementHashing = Hashing.defaultFor(),
                                        elementOrder = null,
                                    ) {
                                        if (subdim < dim - 1u) addAllFrom(section.value.facesOfDimension(subdim))
                                        else add(section.value)
                                        polytope.facesOfDimension(subdim).flatMapTo(this) { polytopeMapping[subdim][it].filter { it.sign.isNonPositive() }.map { it.polytope } }
                                    }
                                },
                            ),
                            sign = Sign.Negative,
                        ),
                    )
                }
            }
        )
    }

    return polytopeMapping
}

Построение симплекса на точках вида $(0, \dots, 0, 1, 0 \dots, 0)$

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

defaultSimplex.kt

fun <Number, Vector> defaultSimplex(
    dimension: UInt,
    basis: VectorSpaceBasis.Finite<Number, Vector>,
): AbstractPolytopicConstruction<PointWrapper<Vector>> =
    AbstractPolytopicConstruction<PointWrapper<Vector>>(dimension).apply {
        val vertices = KoneList(dimension) { index -> addVertex(PointWrapper(basis[index])) }
        val faces = KoneArrayFixedCapacityList<KoneMap<KoneList<UInt>, AbstractPolytopicConstruction.Polytope<PointWrapper<Vector>>>>(dimension)
        faces.add(
            vertices.withIndex().associate(
                keyEquality = KoneList.equality(UInt.context),
                keyHashing = KoneList.hashing(UInt.context),
            ) { (index, vertex) ->
                KoneList(dimension) { if (it == index) 1u else 0u } mapsTo vertex.asPolytope()
            }
        )
        for (dim in 2u .. dimension) faces.add(
            KoneMap.build(
                keyEquality = KoneList.equality(UInt.context),
                keyHashing = KoneList.hashing(UInt.context),
            ) {
                for (flags in KoneList(dimension) { if (it < dim) 1u else 0u }.permutationsWithoutRepetitions(
                    equality = UInt.context
                ))
                    this[flags] = addPolytope(
                        dimension = dim - 1u,
                        vertices = KoneReifiedSet.build(
                            elementEquality = Equality.absoluteFor(),
                            elementHashing = Hashing.defaultFor(),
                        ) {
                            for ((index, flag) in flags.withIndex()) if (flag == 1u) add(vertices[index])
                        },
                        faces = KoneList(dim - 1u) {
                            KoneMutableReifiedSet.of<AbstractPolytopicConstruction.Polytope<PointWrapper<Vector>>>(
                                elementEquality = Equality.absoluteFor(),
                                elementHashing = Hashing.defaultFor(),
                            )
                        }.apply {
                            for (subflags in cartesianProduct(flags.map { (0u .. it).toKoneList() })) if (UInt.context { subflags.sum() } in 1u ..< dim)
                                this[UInt.context { subflags.sum() } - 1u].add(faces[UInt.context { subflags.sum() } - 1u][subflags])
                        }
                    )
            }
        )
    }

Финальный политоп

areaToCheckIntersectionWith.kt

context(_: Field<Number>, _: Order<Number>, _: AffineSpaceOverField<Number, Vector, PointWrapper<Vector>>)
fun <Number, Vector> areaToCheckIntersectionWith(
    dimension: UInt,
    basis: VectorSpaceBasis.Finite<Number, Vector>
): AbstractPolytopicConstruction.Polytope<PointWrapper<Vector>> {
    var polytopicConstruction = defaultSimplex(dimension, basis)

    for (index in 0u ..< dimension) {
        val newPolytopicConstruction = AbstractPolytopicConstruction<PointWrapper<Vector>>(dimension)

        val polytopeMapping = polytopicConstruction.divideByAndPut(
            linearFunction = {  point ->
                val vectorDecomposition = basis.decompose(point.vector)
                context(contextOf<Monoid<Number>>()) {
                    (0u ..< dimension).toKoneList().sumOf { vectorDecomposition[it] } - vectorDecomposition[index] * 2
                }
            },
            target = newPolytopicConstruction
        )

        for (vertexPolytope in polytopeMapping[0u].valuesView.flatMap { it.filter { it.sign.isNegative() }.map { it.polytope } })
            vertexPolytope.remove()

        polytopicConstruction = newPolytopicConstruction
    }

    return polytopicConstruction.polytopesOfDimension(dimension - 1u).single()
}

🔵 Отступление: построение базиса решений системы уравнений в приведённом ступенчатом виде

Чтобы проверить пересечение множества решений данной системы уравнений (в приведённом ступенчатом виде) в гиперплоскости $\sum_{i=1}^n x_i = 1$ с политопом возможных объёмов фасет симплекса, нам потребуется выразить это множество решений как аффинную оболочку некоторых точек (во всё той же гиперплоскости). Для этого стоит вычислить какой-нибудь базис этого пространства решений и "спроецировать" его на гиперплоскость. Займёмся нахождением базиса.

(Будем задавать систему уравнений множеством из $m$ векторов-строк $v_t = (v_{t, i})_{i=1}^n$, для которых ищем вектор-столбец $a$, что $v_t \cdot a = (0)$.) Поскольку наша система уравнений приведена к приведённой ступенчатой форме со ступеньками в столбцах $k_1$, ..., $k_m$, то для каждого столбца с номером $i \in (k_t, k_{t+1})$ можно рассмотреть столбец $a = (a_1, \dots, a_n)^T$, где

$$a_j = \begin{cases} v_{s,i}& \text{если } j = k_s, s \leqslant t\\ -1& \text{если } j = i\\ 0& \text{иначе} \end{cases}$$

Тогда из ступенчатости будет понятно, что для всех $s > t$ верно $v_s \cdot a = (0)$, так как все коэффициенты $v_s$ равны нулю до столбца $k_s$ исключительно, а все коэффициенты $a$ равны нулю со строки $i < k_{i+1} \leqslant k_s$ исключительно. Также из приведённой ступенчатости понятно, что для всех $s \leqslant t$ верно $v_s \cdot a = (0)$, так как $v_s$ имеет возможно ненулевые коэффициенты в столбцах $k_1$, ..., $k_t$ и коэффициент $-1$ в столбце $i$, но $v_s$ имеет во всех этих столбцах кроме $k_s$ и, быть может, $i$ нулевые коэффициенты, отчего

$$v_s \cdot a = (v_{s, k_s} \cdot a_{k_s} + v_{s, i} \cdot a_i) = (1 \cdot v_{s, i} + v_{s, i} \cdot (-1)) = (0)$$

При этом последний ненулевой коэффициент $a$ находится в строке $i$. Значит, для всех $i \notin \{k_t\}_{t=1}^m$ будет получаться по вектору $a$, и все эти векторы будут линейно независимы.

Программная реализация

Заведём для этого функцию.

reducedRowEchelonFormKernel.kt

fun <
    Vector,
    Covector,
> KoneList<Covector>.reducedRowEchelonFormKernel(
): KoneList<Vector> {
    TODO()
}

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

reducedRowEchelonFormKernel.kt

context(_: Field<Number>, vectorSpace: VectorSpace<Number, Vector>)
fun <
    Number,
    Vector,
    Covector,
> KoneList<Covector>.reducedRowEchelonFormKernel(
    dimension: UInt,
    vectorBasis: VectorSpaceBasis.Finite<Number, Vector>,
    covectorBasisDecomposition: VectorSpaceBasisDecomposition.Finite<Number, Covector>,
): KoneList<Vector> {
    TODO()
}

Для удобства переложим ковектора в переменную cobasis. И заодно сразу вернём результат функции как список векторов (который мы пока не умеем строить).

reducedRowEchelonFormKernel.kt

context(_: Field<Number>, vectorSpace: VectorSpace<Number, Vector>)
fun <
    Number,
    Vector,
    Covector,
> KoneList<Covector>.reducedRowEchelonFormKernel(
    dimension: UInt,
    vectorBasis: VectorSpaceBasis.Finite<Number, Vector>,
    covectorBasisDecomposition: VectorSpaceBasisDecomposition.Finite<Number, Covector>,
): KoneList<Vector> {
    val cobasis = this
    return KoneList.build(dimension - cobasis.size) {
        TODO()
    }
}

Для построения векторов по принципу описанном в этом разделе нам нужно пройтись по всем столбцам и, если этот столбец не является столбцом ступеньки, построить соответствующий ему вектор. Для каждого такого столбца нужно понять, какие перед этим были ступеньки. Поэтому будем одновременно на каждом новом столбце либо фиксировать новую ступеньку, либо генерить вектор. Поэтому сначала точно заведём список startIndices под найденный столбцы ступенек и индекс basisIndex рассматриваемой в данный момент строки (или псевдо-строки с индексом $m$), в которой мы ещё не дошли до ступеньки и начнём итерироваться по столбцам.

reducedRowEchelonFormKernel.kt

context(_: Field<Number>, vectorSpace: VectorSpace<Number, Vector>)
fun <
    Number,
    Vector,
    Covector,
> KoneList<Covector>.reducedRowEchelonFormKernel(
    dimension: UInt,
    vectorBasis: VectorSpaceBasis.Finite<Number, Vector>,
    covectorBasisDecomposition: VectorSpaceBasisDecomposition.Finite<Number, Covector>,
): KoneList<Vector> {
    val cobasis = this
    return KoneList.build(dimension - cobasis.size) {
        val startIndices = KoneArrayFixedCapacityList<UInt>(cobasis.size)
        var basisIndex = 0u
        for (coefficientIndex in 0u ..< dimension) {
            TODO()
        }
    }
}

Пусть мы попали на начало ступеньки. Чтобы это проверить, нужно посмотреть, что в строке с индексом basisIndex (если такая вообще есть) и в столбце с индексом coefficientIndex стоит что-то ненулевое. Попав на начало ступеньки, проверим, что там стоит единица, добавим этот столбец в startIndices и перейдём к следующим строке и столбцу.

reducedRowEchelonFormKernel.kt

context(_: Field<Number>, vectorSpace: VectorSpace<Number, Vector>)
fun <
    Number,
    Vector,
    Covector,
> KoneList<Covector>.reducedRowEchelonFormKernel(
    dimension: UInt,
    vectorBasis: VectorSpaceBasis.Finite<Number, Vector>,
    covectorBasisDecomposition: VectorSpaceBasisDecomposition.Finite<Number, Covector>,
): KoneList<Vector> {
    val cobasis = this
    return KoneList.build(dimension - cobasis.size) {
        val startIndices = KoneArrayFixedCapacityList<UInt>(cobasis.size)
        var basisIndex = 0u
        for (coefficientIndex in 0u ..< dimension) {
            if (basisIndex < cobasis.size && covectorBasisDecomposition.decompose(cobasis[basisIndex])[coefficientIndex].isNotZero()) {
                check(covectorBasisDecomposition.decompose(cobasis[basisIndex])[coefficientIndex].isOne())
                startIndices.add(coefficientIndex)
                basisIndex++
                continue
            }

            TODO()
        }
    }
}

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

reducedRowEchelonFormKernel.kt

context(_: Field<Number>, vectorSpace: VectorSpace<Number, Vector>)
fun <
    Number,
    Vector,
    Covector,
> KoneList<Covector>.reducedRowEchelonFormKernel(
    dimension: UInt,
    vectorBasis: VectorSpaceBasis.Finite<Number, Vector>,
    covectorBasisDecomposition: VectorSpaceBasisDecomposition.Finite<Number, Covector>,
): KoneList<Vector> {
    val cobasis = this
    return KoneList.build(dimension - cobasis.size) {
        val startIndices = KoneArrayFixedCapacityList<UInt>(cobasis.size)
        var basisIndex = 0u
        for (coefficientIndex in 0u ..< dimension) {
            if (basisIndex < cobasis.size && covectorBasisDecomposition.decompose(cobasis[basisIndex])[coefficientIndex].isNotZero()) {
                check(covectorBasisDecomposition.decompose(cobasis[basisIndex])[coefficientIndex].isOne())
                startIndices.add(coefficientIndex)
                basisIndex++
                continue
            }

            var vector = vectorSpace.zero
            for ((index, start) in startIndices.withIndex())
                vector += covectorBasisDecomposition.decompose(cobasis[index])[coefficientIndex] * vectorBasis[start]
            vector += valueOf(-1) * vectorBasis[coefficientIndex]

            add(vector)
        }
    }
}

🔵 Отступление: проверка пересекаемости аффинной оболочки данных точек и данного выпуклого политопа в данном аффинном подпространстве

Для проверки того, что система уравнений имеет решения в области, ограниченной фиксированным набором неравенств, мы будем проецировать всё на гиперплоскость $\sum_{i=1}^n x_i = 1$, так как и уравнения, и неравенства однородны по вектору $(x_i)_{i=1}^n$. В таком случае область, ограниченная неравенствами становится внутренностью выпуклого политопа, а множество решений системы уравнений превращается в аффинное подпространство данной гиперплоскости. Получается, что нужно проверить пересекаемость внутренности одного и того же выпуклого политопа и некоторого (каждый раз нового) аффинного подпространства.

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

Но сначала для всего этого нам понадобятся промежуточные шаги.

Шаг: найти любой флаг данного политопа

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

anyFlag.kt

fun <
    Polytope: PolytopicConstruction.Polytope<*, Polytope, *>
> Polytope.anyFlag(): KoneList<Polytope> =
    KoneList.build {
        add(this@anyFlag)
        repeat(this@anyFlag.dimension) {
            val lastPolytope = last()
            +lastPolytope.facesOfDimension(lastPolytope.dimension - 1u).first()
        }
        reverse()
    }

Шаг: для данного политопа и его грани найти в его аффинной оболочке вектор, перпендикулярный данной грани

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

Заведём для этого функцию с простой проверкой входных параметров:

innerVectorPerpendicularTo.kt

context(_: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
    Polytope: PolytopicConstruction.Polytope<Point, Polytope, *>
> Polytope.innerVectorPerpendicularTo(face: Polytope): Vector {
    require(this.dimension > face.dimension && face in this.facesOfDimension(face.dimension))
    TODO()
}

Возьмём флаг грани с добавленным после него данным политопом, чтобы обработать всё одним циклом и получить искомый базис сразу.

innerVectorPerpendicularTo.kt

context(_: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
    Polytope: PolytopicConstruction.Polytope<Point, Polytope, *>
> Polytope.innerVectorPerpendicularTo(face: Polytope): Vector {
    require(this.dimension > face.dimension && face in this.facesOfDimension(face.dimension))
    val flagAndPolytope = KoneList.build {
        addAllFrom(face.anyFlag())
        add(this@innerVectorPerpendicularTo)
    }
    TODO()
}

Возьмём начальную точку:

innerVectorPerpendicularTo.kt

context(_: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
    Polytope: PolytopicConstruction.Polytope<Point, Polytope, *>
> Polytope.innerVectorPerpendicularTo(face: Polytope): Vector {
    require(this.dimension > face.dimension && face in this.facesOfDimension(face.dimension))
    val flagAndPolytope = KoneList.build {
        addAllFrom(face.anyFlag())
        add(this@innerVectorPerpendicularTo)
    }
    val startPoint = flagAndPolytope[0u].vertices.single().position
    TODO()
}

Построим базис:

innerVectorPerpendicularTo.kt

context(_: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
    Polytope: PolytopicConstruction.Polytope<Point, Polytope, *>
> Polytope.innerVectorPerpendicularTo(face: Polytope): Vector {
    require(this.dimension > face.dimension && face in this.facesOfDimension(face.dimension))
    val flagAndPolytope = KoneList.build {
        addAllFrom(face.anyFlag())
        add(this@innerVectorPerpendicularTo)
    }
    val startPoint = flagAndPolytope[0u].vertices.single().position
    KoneList.build<Vector> {
        for (index in 0u ..< flagAndPolytope.lastIndex) {
            val innerPolytope = flagAndPolytope[index]
            val outerPolytope = flagAndPolytope[index + 1u]
            +(outerPolytope.vertices.firstThat { it !in innerPolytope.vertices }.position - startPoint)
        }
    }
    TODO()
}

Остаётся только ортогонализировать, взять последний вектор и вернуть его:

innerVectorPerpendicularTo.kt

context(_: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
    Polytope: PolytopicConstruction.Polytope<Point, Polytope, *>
> Polytope.innerVectorPerpendicularTo(face: Polytope): Vector {
    require(this.dimension > face.dimension && face in this.facesOfDimension(face.dimension))
    val flagAndPolytope = KoneList.build {
        addAllFrom(face.anyFlag())
        add(this@innerVectorPerpendicularTo)
    }
    val startPoint = flagAndPolytope[0u].vertices.single().position
    return KoneList.build<Vector> {
        for (index in 0u ..< flagAndPolytope.lastIndex) {
            val innerPolytope = flagAndPolytope[index]
            val outerPolytope = flagAndPolytope[index + 1u]
            +(outerPolytope.vertices.firstThat { it !in innerPolytope.vertices }.position - startPoint)
        }
    }.orthogonalizeByGramSchmidt().last()
}

Шаг: для данного политопа, точки вне его аффинной оболочки и некоторого вектора спроецировать политоп через точку на плоскость, перпендикулярную вектору

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

Заведём для этого функцию.

projectThroughPointInGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    TODO()
}

Сначала нужно проверить, есть ли искомая проекция уже в реестре. И если есть, то вернуть её.

projectThroughPointInGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    TODO()
}

Затем отработаем ноль-мерный случай, когда нужно спроецировать одну единственную вершину.

projectThroughPointInGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    if (area.dimension == 0u)
        return target.addVertex(
            area.vertices.single().position.let {
                val vector = (it - point)
                point + vector / (vector dot normal)
            }
        ).asPolytope().also { polytopesRegistry[0u][area] = it }

    TODO()
}

Хорошо, значит, теперь у нас случай размерности хотя бы 1. Тогда сначала спроецируем все фасеты.

projectThroughPointInGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    if (area.dimension == 0u)
        return target.addVertex(
            area.vertices.single().position.let {
                val vector = (it - point)
                point + vector / (vector dot normal)
            }
        ).asPolytope().also { polytopesRegistry[0u][area] = it }

    val newFacets = area.facesOfDimension(area.dimension - 1u).map { facet ->
        projectThroughPointInGeneralPositionOntoNormalPlane(
            point = point,
            normal = normal,
            area = facet,
            polytopesRegistry = polytopesRegistry,
            target = target,
        )
    }

    TODO()
}

Теперь соберём для проекции нашего политопа все его грани каждой размерности.

projectThroughPointInGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    if (area.dimension == 0u)
        return target.addVertex(
            area.vertices.single().position.let {
                val vector = (it - point)
                point + vector / (vector dot normal)
            }
        ).asPolytope().also { polytopesRegistry[0u][area] = it }

    val newFacets = area.facesOfDimension(area.dimension - 1u).map { facet ->
        projectThroughPointInGeneralPositionOntoNormalPlane(
            point = point,
            normal = normal,
            area = facet,
            polytopesRegistry = polytopesRegistry,
            target = target,
        )
    }

    val newFaces = KoneList(area.dimension) {
        KoneMutableReifiedSet.of<AbstractPolytopicConstruction.Polytope<Point>>(
            elementEquality = Equality.absoluteFor(),
            elementHashing = Hashing.defaultFor(),
        )
    }

    for (facet in newFacets) {
        for (dim in 0u ..< facet.dimension) newFaces[dim].addAllFrom(facet.facesOfDimension(dim))
        newFaces[facet.dimension].add(facet)
    }

    TODO()
}

Теперь только осталось построить проекцию политопа, сохранить её в реестр и вернуть.

projectThroughPointInGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    if (area.dimension == 0u)
        return target.addVertex(
            area.vertices.single().position.let {
                val vector = (it - point)
                point + vector / (vector dot normal)
            }
        ).asPolytope().also { polytopesRegistry[0u][area] = it }

    val newFacets = area.facesOfDimension(area.dimension - 1u).map { facet ->
        projectThroughPointInGeneralPositionOntoNormalPlane(
            point = point,
            normal = normal,
            area = facet,
            polytopesRegistry = polytopesRegistry,
            target = target,
        )
    }

    val newFaces = KoneList(area.dimension) {
        KoneMutableReifiedSet.of<AbstractPolytopicConstruction.Polytope<Point>>(
            elementEquality = Equality.absoluteFor(),
            elementHashing = Hashing.defaultFor(),
        )
    }

    for (facet in newFacets) {
        for (dim in 0u ..< facet.dimension) newFaces[dim].addAllFrom(facet.facesOfDimension(dim))
        newFaces[facet.dimension].add(facet)
    }

    return target.addPolytope(
        dimension = area.dimension,
        vertices = newFaces[0u].mapTo(
            KoneMutableReifiedSet.of<AbstractPolytopicConstruction.Vertex<Point>>(
                elementEquality = Equality.absoluteFor(),
                elementHashing = Hashing.defaultFor(),
            )
        ) { it.vertices.single() },
        faces = newFaces,
    ).also { polytopesRegistry[area.dimension][area] = it }
}

Шаг: для данного политопа, точки в его аффинной оболочке и некоторого вектора спроецировать политоп через точку на плоскость, перпендикулярную вектору

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

Заведём для этого функцию.

projectThroughPointInNonGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInNonGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    TODO()
}

Сначала нужно проверить, есть ли искомая проекция уже в реестре. И если есть, то вернуть её.

projectThroughPointInNonGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInNonGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    TODO()
}

Дальше обработаем 0- и 1-мерный случаи. В 0-мерном случае задача не определена, поэтому надо выдать ошибку. В 1-мерном случае надо просто взять вектор направления данной прямой и спроецировать его.

projectThroughPointInNonGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInNonGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    when (area.dimension) {
        0u -> error("Cannot have normal vector in 0-dimensional non-general position")
        1u -> {
            val vector = area.vertices.first().position - point
            return target.addVertex(point + vector / (vector dot normal)).asPolytope().also { polytopesRegistry[1u][area] = it }
        }
    }

    TODO()
}

Дальше соберём все фасеты политопа и присвоим им знаки.

projectThroughPointInNonGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInNonGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    when (area.dimension) {
        0u -> error("Cannot have normal vector in 0-dimensional non-general position")
        1u -> {
            val vector = area.vertices.first().position - point
            return target.addVertex(point + vector / (vector dot normal)).asPolytope().also { polytopesRegistry[1u][area] = it }
        }
    }

    data class FacetWithSign(
        val facet: AbstractPolytopicConstruction.Polytope<Point>,
        val sign: Sign,
    )

    val facetsSigns = area.facesOfDimension(area.dimension - 1u).map { facet ->
        val start = facet.vertices.first().position
        val facetNormal = area.innerVectorPerpendicularTo(facet)

        FacetWithSign(facet, ((point - start) dot facetNormal).sign())
    }

    TODO()
}

После чего просто выделим положительные фасеты и нулевые фасеты.

projectThroughPointInNonGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInNonGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    when (area.dimension) {
        0u -> error("Cannot have normal vector in 0-dimensional non-general position")
        1u -> {
            val vector = area.vertices.first().position - point
            return target.addVertex(point + vector / (vector dot normal)).asPolytope().also { polytopesRegistry[1u][area] = it }
        }
    }

    data class FacetWithSign(
        val facet: AbstractPolytopicConstruction.Polytope<Point>,
        val sign: Sign,
    )

    val facetsSigns = area.facesOfDimension(area.dimension - 1u).map { facet ->
        val start = facet.vertices.first().position
        val facetNormal = area.innerVectorPerpendicularTo(facet)

        FacetWithSign(facet, ((point - start) dot facetNormal).sign())
    }

    val positiveFacets = facetsSigns.filter { it.sign.isPositive() }.map { it.facet }
    val zeroFacets = facetsSigns.filter { it.sign.isZero() }.map { it.facet }

    TODO()
}

Теперь начнём собирать грани политопа-проекции и искомые хребты для проецирования.

projectThroughPointInNonGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInNonGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    when (area.dimension) {
        0u -> error("Cannot have normal vector in 0-dimensional non-general position")
        1u -> {
            val vector = area.vertices.first().position - point
            return target.addVertex(point + vector / (vector dot normal)).asPolytope().also { polytopesRegistry[1u][area] = it }
        }
    }

    data class FacetWithSign(
        val facet: AbstractPolytopicConstruction.Polytope<Point>,
        val sign: Sign,
    )

    val facetsSigns = area.facesOfDimension(area.dimension - 1u).map { facet ->
        val start = facet.vertices.first().position
        val facetNormal = area.innerVectorPerpendicularTo(facet)

        FacetWithSign(facet, ((point - start) dot facetNormal).sign())
    }

    val positiveFacets = facetsSigns.filter { it.sign.isPositive() }.map { it.facet }
    val zeroFacets = facetsSigns.filter { it.sign.isZero() }.map { it.facet }

    val faces = KoneList(area.dimension - 1u) {
        KoneMutableReifiedSet.of<AbstractPolytopicConstruction.Polytope<Point>>(
            elementReification = Reification.defaultFor(),
            elementEquality = Equality.absoluteFor(),
            elementHashing = Hashing.defaultFor(),
        )
    }

    val ridgesToProject = KoneMutableSet.of<AbstractPolytopicConstruction.Polytope<Point>>(
        elementEquality = Equality.absoluteFor(),
        elementHashing = Hashing.defaultFor(),
    )

    TODO()
}

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

projectThroughPointInNonGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInNonGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    when (area.dimension) {
        0u -> error("Cannot have normal vector in 0-dimensional non-general position")
        1u -> {
            val vector = area.vertices.first().position - point
            return target.addVertex(point + vector / (vector dot normal)).asPolytope().also { polytopesRegistry[1u][area] = it }
        }
    }

    data class FacetWithSign(
        val facet: AbstractPolytopicConstruction.Polytope<Point>,
        val sign: Sign,
    )

    val facetsSigns = area.facesOfDimension(area.dimension - 1u).map { facet ->
        val start = facet.vertices.first().position
        val facetNormal = area.innerVectorPerpendicularTo(facet)

        FacetWithSign(facet, ((point - start) dot facetNormal).sign())
    }

    val positiveFacets = facetsSigns.filter { it.sign.isPositive() }.map { it.facet }
    val zeroFacets = facetsSigns.filter { it.sign.isZero() }.map { it.facet }

    val faces = KoneList(area.dimension - 1u) {
        KoneMutableReifiedSet.of<AbstractPolytopicConstruction.Polytope<Point>>(
            elementReification = Reification.defaultFor(),
            elementEquality = Equality.absoluteFor(),
            elementHashing = Hashing.defaultFor(),
        )
    }

    val ridgesToProject = KoneMutableSet.of<AbstractPolytopicConstruction.Polytope<Point>>(
        elementEquality = Equality.absoluteFor(),
        elementHashing = Hashing.defaultFor(),
    )

    for (facet in positiveFacets) for (ridge in facet.facesOfDimension(facet.dimension - 1u))
        if (ridge in ridgesToProject) ridgesToProject.remove(ridge)
        else ridgesToProject.add(ridge)

    TODO()
}

Теперь переберём нулевые фасеты и сделаем две вещи:

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

projectThroughPointInNonGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInNonGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    when (area.dimension) {
        0u -> error("Cannot have normal vector in 0-dimensional non-general position")
        1u -> {
            val vector = area.vertices.first().position - point
            return target.addVertex(point + vector / (vector dot normal)).asPolytope().also { polytopesRegistry[1u][area] = it }
        }
    }

    data class FacetWithSign(
        val facet: AbstractPolytopicConstruction.Polytope<Point>,
        val sign: Sign,
    )

    val facetsSigns = area.facesOfDimension(area.dimension - 1u).map { facet ->
        val start = facet.vertices.first().position
        val facetNormal = area.innerVectorPerpendicularTo(facet)

        FacetWithSign(facet, ((point - start) dot facetNormal).sign())
    }

    val positiveFacets = facetsSigns.filter { it.sign.isPositive() }.map { it.facet }
    val zeroFacets = facetsSigns.filter { it.sign.isZero() }.map { it.facet }

    val faces = KoneList(area.dimension - 1u) {
        KoneMutableReifiedSet.of<AbstractPolytopicConstruction.Polytope<Point>>(
            elementReification = Reification.defaultFor(),
            elementEquality = Equality.absoluteFor(),
            elementHashing = Hashing.defaultFor(),
        )
    }

    val ridgesToProject = KoneMutableSet.of<AbstractPolytopicConstruction.Polytope<Point>>(
        elementEquality = Equality.absoluteFor(),
        elementHashing = Hashing.defaultFor(),
    )

    for (facet in positiveFacets) for (ridge in facet.facesOfDimension(facet.dimension - 1u))
        if (ridge in ridgesToProject) ridgesToProject.remove(ridge)
        else ridgesToProject.add(ridge)

    for (facet in zeroFacets) {
        ridgesToProject.removeAllFrom(facet.facesOfDimension(facet.dimension - 1u))

        val projection = projectThroughPointInNonGeneralPositionOntoNormalPlane(
            point = point,
            normal = normal,
            area = facet,
            polytopesRegistry = polytopesRegistry,
            target = target,
        )

        for (dim in 0u ..< projection.dimension) faces[dim].addAllFrom(projection.facesOfDimension(dim))
        faces[projection.dimension].add(projection)
    }

    TODO()
}

Теперь, наконец, спроецируем хребты.

projectThroughPointInNonGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInNonGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    when (area.dimension) {
        0u -> error("Cannot have normal vector in 0-dimensional non-general position")
        1u -> {
            val vector = area.vertices.first().position - point
            return target.addVertex(point + vector / (vector dot normal)).asPolytope().also { polytopesRegistry[1u][area] = it }
        }
    }

    data class FacetWithSign(
        val facet: AbstractPolytopicConstruction.Polytope<Point>,
        val sign: Sign,
    )

    val facetsSigns = area.facesOfDimension(area.dimension - 1u).map { facet ->
        val start = facet.vertices.first().position
        val facetNormal = area.innerVectorPerpendicularTo(facet)

        FacetWithSign(facet, ((point - start) dot facetNormal).sign())
    }

    val positiveFacets = facetsSigns.filter { it.sign.isPositive() }.map { it.facet }
    val zeroFacets = facetsSigns.filter { it.sign.isZero() }.map { it.facet }

    val faces = KoneList(area.dimension - 1u) {
        KoneMutableReifiedSet.of<AbstractPolytopicConstruction.Polytope<Point>>(
            elementReification = Reification.defaultFor(),
            elementEquality = Equality.absoluteFor(),
            elementHashing = Hashing.defaultFor(),
        )
    }

    val ridgesToProject = KoneMutableSet.of<AbstractPolytopicConstruction.Polytope<Point>>(
        elementEquality = Equality.absoluteFor(),
        elementHashing = Hashing.defaultFor(),
    )

    for (facet in positiveFacets) for (ridge in facet.facesOfDimension(facet.dimension - 1u))
        if (ridge in ridgesToProject) ridgesToProject.remove(ridge)
        else ridgesToProject.add(ridge)

    for (facet in zeroFacets) {
        ridgesToProject.removeAllFrom(facet.facesOfDimension(facet.dimension - 1u))

        val projection = projectThroughPointInNonGeneralPositionOntoNormalPlane(
            point = point,
            normal = normal,
            area = facet,
            polytopesRegistry = polytopesRegistry,
            target = target,
        )

        for (dim in 0u ..< projection.dimension) faces[dim].addAllFrom(projection.facesOfDimension(dim))
        faces[projection.dimension].add(projection)
    }

    for (ridge in ridgesToProject) {
        val projection = projectThroughPointInGeneralPositionOntoNormalPlane(
            point = point,
            normal = normal,
            area = ridge,
            polytopesRegistry = polytopesRegistry,
            target = target,
        )

        for (dim in 0u ..< projection.dimension) faces[dim].addAllFrom(projection.facesOfDimension(dim))
        faces[projection.dimension].add(projection)
    }

    TODO()
}

Теперь осталось только собрать проекцию всего политопа, зарегистрировать её и вернуть.

projectThroughPointInNonGeneralPositionOntoNormalPlane.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
fun <
    Number,
    Vector,
    Point,
> projectThroughPointInNonGeneralPositionOntoNormalPlane(
    point: Point,
    normal: Vector,
    area: AbstractPolytopicConstruction.Polytope<Point>,
    polytopesRegistry: KoneList<KoneMutableMap<AbstractPolytopicConstruction.Polytope<Point>, AbstractPolytopicConstruction.Polytope<Point>>>,
    target: AbstractPolytopicConstruction<Point>,
): AbstractPolytopicConstruction.Polytope<Point> {
    polytopesRegistry[area.dimension].getOrNull(area)?.let { return it }

    when (area.dimension) {
        0u -> error("Cannot have normal vector in 0-dimensional non-general position")
        1u -> {
            val vector = area.vertices.first().position - point
            return target.addVertex(point + vector / (vector dot normal)).asPolytope().also { polytopesRegistry[1u][area] = it }
        }
    }

    data class FacetWithSign(
        val facet: AbstractPolytopicConstruction.Polytope<Point>,
        val sign: Sign,
    )

    val facetsSigns = area.facesOfDimension(area.dimension - 1u).map { facet ->
        val start = facet.vertices.first().position
        val facetNormal = area.innerVectorPerpendicularTo(facet)

        FacetWithSign(facet, ((point - start) dot facetNormal).sign())
    }

    val positiveFacets = facetsSigns.filter { it.sign.isPositive() }.map { it.facet }
    val zeroFacets = facetsSigns.filter { it.sign.isZero() }.map { it.facet }

    val faces = KoneList(area.dimension - 1u) {
        KoneMutableReifiedSet.of<AbstractPolytopicConstruction.Polytope<Point>>(
            elementReification = Reification.defaultFor(),
            elementEquality = Equality.absoluteFor(),
            elementHashing = Hashing.defaultFor(),
        )
    }

    val ridgesToProject = KoneMutableSet.of<AbstractPolytopicConstruction.Polytope<Point>>(
        elementEquality = Equality.absoluteFor(),
        elementHashing = Hashing.defaultFor(),
    )

    for (facet in positiveFacets) for (ridge in facet.facesOfDimension(facet.dimension - 1u))
        if (ridge in ridgesToProject) ridgesToProject.remove(ridge)
        else ridgesToProject.add(ridge)

    for (facet in zeroFacets) {
        ridgesToProject.removeAllFrom(facet.facesOfDimension(facet.dimension - 1u))

        val projection = projectThroughPointInNonGeneralPositionOntoNormalPlane(
            point = point,
            normal = normal,
            area = facet,
            polytopesRegistry = polytopesRegistry,
            target = target,
        )

        for (dim in 0u ..< projection.dimension) faces[dim].addAllFrom(projection.facesOfDimension(dim))
        faces[projection.dimension].add(projection)
    }

    for (ridge in ridgesToProject) {
        val projection = projectThroughPointInGeneralPositionOntoNormalPlane(
            point = point,
            normal = normal,
            area = ridge,
            polytopesRegistry = polytopesRegistry,
            target = target,
        )

        for (dim in 0u ..< projection.dimension) faces[dim].addAllFrom(projection.facesOfDimension(dim))
        faces[projection.dimension].add(projection)
    }

    return target.addPolytope(
        dimension = area.dimension - 1u,
        vertices = faces.first().mapTo(
            KoneMutableReifiedSet.of(
                elementReification = Reification.defaultFor(),
                elementEquality = Equality.absoluteFor(),
                elementHashing = Hashing.defaultFor(),
            )
        ) { it.vertices.single() },
        faces = faces,
    ).also { polytopesRegistry[area.dimension][area] = it }
}

Финальный шаг

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

Заведём для этого функцию.

checkIntersectionOfAreaAndSolutions.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
tailrec fun <
    Number,
    Vector,
    Point,
> checkIntersectionOfAreaAndSolutions(
    area: AbstractPolytopicConstruction.Polytope<Point>,
    solutions: KoneList<Point>,
): Boolean {
    TODO()
}

Первым делом проверим, что аффинная оболочка непустая. Если оно пустое, то можно сразу вернуть false. Иначе можно выделить первую точку из списка и все остальные. При этом если размерность выпуклой оболочки (и всего пространства) равна 0, то очевиден ответ true.

checkIntersectionOfAreaAndSolutions.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
tailrec fun <
    Number,
    Vector,
    Point,
> checkIntersectionOfAreaAndSolutions(
    area: AbstractPolytopicConstruction.Polytope<Point>,
    solutions: KoneList<Point>,
): Boolean {
    if (solutions.isEmpty()) return false

    val primarySolution = solutions.first()
    val otherSolutions = solutions.drop(1u)

    if (area.dimension == 0u) return true

    TODO()
}

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

checkIntersectionOfAreaAndSolutions.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
tailrec fun <
    Number,
    Vector,
    Point,
> checkIntersectionOfAreaAndSolutions(
    area: AbstractPolytopicConstruction.Polytope<Point>,
    solutions: KoneList<Point>,
): Boolean {
    if (solutions.isEmpty()) return false

    val primarySolution = solutions.first()
    val otherSolutions = solutions.drop(1u)

    if (area.dimension == 0u) return true

    data class FacetNormalAndIsStrictlyInside(
        val facet: AbstractPolytopicConstruction.Polytope<Point>,
        val normal: Vector,
        val isStrictlyInside: Boolean,
    )

    val facetNormalAndIsStrictlyInsides = area.facesOfDimension(area.dimension - 1u).map { facet ->
        val start = facet.vertices.first().position
        val facetNormal = area.innerVectorPerpendicularTo(facet)

        val sign = ((primarySolution - start) dot facetNormal).sign()

        if (sign.isNegative()) {
            val nonParallelSolution = otherSolutions.firstThatOrNull { ((it - primarySolution) dot facetNormal).isNotZero() } ?: return false
            val otherSolutionProjections = otherSolutions.map {
                val vector = it - primarySolution
                val vectorDotNormal = vector dot facetNormal
                if (vectorDotNormal.isNotZero()) primarySolution + vector / (vector dot facetNormal)
                else {
                    val newVector = vector + (nonParallelSolution - primarySolution)
                    primarySolution + newVector / (newVector dot facetNormal)
                }
            }
            val areaProjection = projectThroughPointInNonGeneralPositionOntoNormalPlane(
                point = primarySolution,
                normal = facetNormal,
                area = area,
                polytopesRegistry = KoneList(area.dimension + 1u) {
                    KoneMutableMap.of(
                        keyEquality = Equality.absoluteFor(),
                        keyHashing = Hashing.defaultFor(),
                    )
                },
                target = AbstractPolytopicConstruction(facet.dimension),
            )

            return checkIntersectionOfAreaAndSolutions(
                area = areaProjection,
                solutions = otherSolutionProjections,
            )
        }

        FacetNormalAndIsStrictlyInside(facet, facetNormal, sign.isPositive())
    }

    TODO()
}

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

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

checkIntersectionOfAreaAndSolutions.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
tailrec fun <
    Number,
    Vector,
    Point,
> checkIntersectionOfAreaAndSolutions(
    area: AbstractPolytopicConstruction.Polytope<Point>,
    solutions: KoneList<Point>,
): Boolean {
    if (solutions.isEmpty()) return false

    val primarySolution = solutions.first()
    val otherSolutions = solutions.drop(1u)

    if (area.dimension == 0u) return true

    data class FacetNormalAndIsStrictlyInside(
        val facet: AbstractPolytopicConstruction.Polytope<Point>,
        val normal: Vector,
        val isStrictlyInside: Boolean,
    )

    val facetNormalAndIsStrictlyInsides = area.facesOfDimension(area.dimension - 1u).map { facet ->
        val start = facet.vertices.first().position
        val facetNormal = area.innerVectorPerpendicularTo(facet)

        val sign = ((primarySolution - start) dot facetNormal).sign()

        if (sign.isNegative()) {
            val nonParallelSolution = otherSolutions.firstThatOrNull { ((it - primarySolution) dot facetNormal).isNotZero() } ?: return false
            val otherSolutionProjections = otherSolutions.map {
                val vector = it - primarySolution
                val vectorDotNormal = vector dot facetNormal
                if (vectorDotNormal.isNotZero()) primarySolution + vector / (vector dot facetNormal)
                else {
                    val newVector = vector + (nonParallelSolution - primarySolution)
                    primarySolution + newVector / (newVector dot facetNormal)
                }
            }
            val areaProjection = projectThroughPointInNonGeneralPositionOntoNormalPlane(
                point = primarySolution,
                normal = facetNormal,
                area = area,
                polytopesRegistry = KoneList(area.dimension + 1u) {
                    KoneMutableMap.of(
                        keyEquality = Equality.absoluteFor(),
                        keyHashing = Hashing.defaultFor(),
                    )
                },
                target = AbstractPolytopicConstruction(facet.dimension),
            )

            return checkIntersectionOfAreaAndSolutions(
                area = areaProjection,
                solutions = otherSolutionProjections,
            )
        }

        FacetNormalAndIsStrictlyInside(facet, facetNormal, sign.isPositive())
    }

    if (facetNormalAndIsStrictlyInsides.all { it.isStrictlyInside }) return true
    if (otherSolutions.isEmpty()) return false

    TODO()
}

Если функция всё ещё ничего не вернула, то нам не повезло находиться на границе политопа, имея не 0-мерную аффинную оболочку. Значит нужно найти новую точку, которая будет строго вне политопа. Для этого достаточно найти какой-нибудь вектор в нашей аффинной оболочке и фасету, в которой он не будет лежать. Тогда на этой прямой будет точка, лежащая по другую сторону от политопа относительно этой фасеты. Подберём такую фасету и точку вне политопа.

checkIntersectionOfAreaAndSolutions.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
tailrec fun <
    Number,
    Vector,
    Point,
> checkIntersectionOfAreaAndSolutions(
    area: AbstractPolytopicConstruction.Polytope<Point>,
    solutions: KoneList<Point>,
): Boolean {
    if (solutions.isEmpty()) return false

    val primarySolution = solutions.first()
    val otherSolutions = solutions.drop(1u)

    if (area.dimension == 0u) return true

    data class FacetNormalAndIsStrictlyInside(
        val facet: AbstractPolytopicConstruction.Polytope<Point>,
        val normal: Vector,
        val isStrictlyInside: Boolean,
    )

    val facetNormalAndIsStrictlyInsides = area.facesOfDimension(area.dimension - 1u).map { facet ->
        val start = facet.vertices.first().position
        val facetNormal = area.innerVectorPerpendicularTo(facet)

        val sign = ((primarySolution - start) dot facetNormal).sign()

        if (sign.isNegative()) {
            val nonParallelSolution = otherSolutions.firstThatOrNull { ((it - primarySolution) dot facetNormal).isNotZero() } ?: return false
            val otherSolutionProjections = otherSolutions.map {
                val vector = it - primarySolution
                val vectorDotNormal = vector dot facetNormal
                if (vectorDotNormal.isNotZero()) primarySolution + vector / (vector dot facetNormal)
                else {
                    val newVector = vector + (nonParallelSolution - primarySolution)
                    primarySolution + newVector / (newVector dot facetNormal)
                }
            }
            val areaProjection = projectThroughPointInNonGeneralPositionOntoNormalPlane(
                point = primarySolution,
                normal = facetNormal,
                area = area,
                polytopesRegistry = KoneList(area.dimension + 1u) {
                    KoneMutableMap.of(
                        keyEquality = Equality.absoluteFor(),
                        keyHashing = Hashing.defaultFor(),
                    )
                },
                target = AbstractPolytopicConstruction(facet.dimension),
            )

            return checkIntersectionOfAreaAndSolutions(
                area = areaProjection,
                solutions = otherSolutionProjections,
            )
        }

        FacetNormalAndIsStrictlyInside(facet, facetNormal, sign.isPositive())
    }

    if (facetNormalAndIsStrictlyInsides.all { it.isStrictlyInside }) return true
    if (otherSolutions.isEmpty()) return false

    val shiftVector = otherSolutions.first() - primarySolution

    val firstFacetNormalAndIsStrictlyInside = facetNormalAndIsStrictlyInsides.firstThat { (it.normal dot shiftVector).isNotZero() }
    val normal = firstFacetNormalAndIsStrictlyInside.normal
    val start = firstFacetNormalAndIsStrictlyInside.facet.vertices.first().position
    val newPrimarySolution = primarySolution + shiftVector * (((start - primarySolution) dot normal) / (shiftVector dot normal) - (shiftVector dot normal))

    TODO()
}

Теперь только и осталось всё спроецировать и повторить итерацию.

checkIntersectionOfAreaAndSolutions.kt

context(field: Field<Number>, _: Order<Number>, _: EuclideanSpaceOverField<Number, Vector, Point>)
tailrec fun <
    Number,
    Vector,
    Point,
> checkIntersectionOfAreaAndSolutions(
    area: AbstractPolytopicConstruction.Polytope<Point>,
    solutions: KoneList<Point>,
): Boolean {
    if (solutions.isEmpty()) return false

    val primarySolution = solutions.first()
    val otherSolutions = solutions.drop(1u)

    if (area.dimension == 0u) return true

    data class FacetNormalAndIsStrictlyInside(
        val facet: AbstractPolytopicConstruction.Polytope<Point>,
        val normal: Vector,
        val isStrictlyInside: Boolean,
    )

    val facetNormalAndIsStrictlyInsides = area.facesOfDimension(area.dimension - 1u).map { facet ->
        val start = facet.vertices.first().position
        val facetNormal = area.innerVectorPerpendicularTo(facet)

        val sign = ((primarySolution - start) dot facetNormal).sign()

        if (sign.isNegative()) {
            val nonParallelSolution = otherSolutions.firstThatOrNull { ((it - primarySolution) dot facetNormal).isNotZero() } ?: return false
            val otherSolutionProjections = otherSolutions.map {
                val vector = it - primarySolution
                val vectorDotNormal = vector dot facetNormal
                if (vectorDotNormal.isNotZero()) primarySolution + vector / (vector dot facetNormal)
                else {
                    val newVector = vector + (nonParallelSolution - primarySolution)
                    primarySolution + newVector / (newVector dot facetNormal)
                }
            }
            val areaProjection = projectThroughPointInNonGeneralPositionOntoNormalPlane(
                point = primarySolution,
                normal = facetNormal,
                area = area,
                polytopesRegistry = KoneList(area.dimension + 1u) {
                    KoneMutableMap.of(
                        keyEquality = Equality.absoluteFor(),
                        keyHashing = Hashing.defaultFor(),
                    )
                },
                target = AbstractPolytopicConstruction(facet.dimension),
            )

            return checkIntersectionOfAreaAndSolutions(
                area = areaProjection,
                solutions = otherSolutionProjections,
            )
        }

        FacetNormalAndIsStrictlyInside(facet, facetNormal, sign.isPositive())
    }

    if (facetNormalAndIsStrictlyInsides.all { it.isStrictlyInside }) return true
    if (otherSolutions.isEmpty()) return false

    val shiftVector = otherSolutions.first() - primarySolution

    val firstFacetNormalAndIsStrictlyInside = facetNormalAndIsStrictlyInsides.firstThat { (it.normal dot shiftVector).isNotZero() }
    val normal = firstFacetNormalAndIsStrictlyInside.normal
    val start = firstFacetNormalAndIsStrictlyInside.facet.vertices.first().position
    val newPrimarySolution = primarySolution + shiftVector * (((start - primarySolution) dot normal) / (shiftVector dot normal) - (shiftVector dot normal))

    val areaProjection = projectThroughPointInNonGeneralPositionOntoNormalPlane(
        point = newPrimarySolution,
        normal = normal,
        area = area,
        polytopesRegistry = KoneList(area.dimension + 1u) {
            KoneMutableMap.of(
                keyEquality = Equality.absoluteFor(),
                keyHashing = Hashing.defaultFor(),
            )
        },
        target = AbstractPolytopicConstruction(area.dimension - 1u),
    )
    val otherSolutionProjections = otherSolutions.map {
        val vector = it - newPrimarySolution
        newPrimarySolution + (vector / (vector dot normal))
    }

    return checkIntersectionOfAreaAndSolutions(
        area = areaProjection,
        solutions = otherSolutionProjections,
    )
}

Фильтрация полученных систем линейных уравнений

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

1. Для каждого полученного ЛНМ оставить ровно одну систему с самым большим набором изначальных ковекторов. (Поскольку есть наибольшая по включению система, то достаточно взять наибольшую по размеру и не заниматься объединением систем.)
2. Для каждой системы заменить полученное ЛНМ на какой-нибудь базис ядра этого ЛНМ и проверить, что ядро пересекается с полученным политопом.

После этого можно просто у всех систем взять размеры и выделить все полученные числа (по одному разу).

Заведём для этого всего функцию:

filterEquationSystems.kt

fun <
    Covector,
> filterEquationSystems(
    equationCases: KoneList<KoneList<EquationCase<Covector>>>
): KoneList<Covector> {
    TODO()
}

Для первого пункта достаточно объединить внутренние списки систем и загрузить полученный список систем в одну мапу (которая сопоставляет каждому ЛНМ наибольшую систему):

filterEquationSystems.kt

context(
    covectorEquality: Equality<Covector>,
    covectorHashing: Hashing<Covector>,
)
fun <
    Covector,
> filterEquationSystems(
    equationCases: KoneList<KoneList<EquationCase<Covector>>>
): KoneList<Covector> {
    val equationToCases = equationCases.flatten().associateBy(
        keyEquality = KoneList.equality(covectorEquality),
        keyHashing = KoneList.hashing(covectorHashing),
        keySelector = { it.basis },
        valueTransform = { it.initialCovectors },
        resolve = { _, currentValue, newValue -> if (currentValue.size > newValue.size) currentValue else newValue }
    )

    TODO()
}

Теперь мы заменяем базис ковекторов на базис векторов-решений.

filterEquationSystems.kt

context(
    covectorEquality: Equality<Covector>,
    covectorHashing: Hashing<Covector>,
)
fun <
    Covector,
> filterEquationSystems(
    equationCases: KoneList<KoneList<EquationCase<Covector>>>
): KoneList<Covector> {
    val equationToCases = equationCases.flatten().associateBy(
        keyEquality = KoneList.equality(covectorEquality),
        keyHashing = KoneList.hashing(covectorHashing),
        keySelector = { it.basis },
        valueTransform = { it.initialCovectors },
        resolve = { _, currentValue, newValue -> if (currentValue.size > newValue.size) currentValue else newValue }
    )

    data class SolutionCase(
        val initialCovectors: KoneList<Covector>,
        val solutions: KoneList<Vector>,
    )

    val solutionCases = equationToCases.nodesView.map { (basis, initialCovectors) ->
        SolutionCase(
            initialCovectors = initialCovectors,
            solutions = basis.reducedRowEchelonFormKernel(
                dimension = dimension,
                vectorBasis = vectorBasis,
                covectorBasisDecomposition = covectorBasisDecomposition,
            ),
        )
    }

    TODO()
}

Далее фильтрация состоит из двух шагов:

1. Спроецировать базис на гиперплоскость $\sum_{i=1}^n x_i = 1$.
2. Проверить пересекаемость проекции с политопом.

Рассмотрим первый пункт. Для всякого базисного вектора $v$, если $\sum_{i=1}^n v_i \neq 0$, то $\frac{v}{\sum_{i=1}^n v_i}$ будет проекцией $v$. Если же $\sum_{i=1}^n v_i = 0$, но есть $u$, что $\sum_{i=1}^n u_i \neq 0$, то вместо $v$ можно рассмотреть $v + u$. Если же такого $u$ нет, то проекция пуста и получается сразу ответ "не имеет решений в политопе".

filterEquationSystems.kt

context(
    covectorEquality: Equality<Covector>,
    covectorHashing: Hashing<Covector>,
)
fun <
    Covector,
> filterEquationSystems(
    equationCases: KoneList<KoneList<EquationCase<Covector>>>
): KoneList<Covector> {
    val equationToCases = equationCases.flatten().associateBy(
        keyEquality = KoneList.equality(covectorEquality),
        keyHashing = KoneList.hashing(covectorHashing),
        keySelector = { it.basis },
        valueTransform = { it.initialCovectors },
        resolve = { _, currentValue, newValue -> if (currentValue.size > newValue.size) currentValue else newValue }
    )

    data class SolutionCase(
        val initialCovectors: KoneList<Covector>,
        val solutions: KoneList<Vector>,
    )

    val solutionCases = equationToCases.nodesView.map { (basis, initialCovectors) ->
        SolutionCase(
            initialCovectors = initialCovectors,
            solutions = basis.reducedRowEchelonFormKernel(
                dimension = dimension,
                vectorBasis = vectorBasis,
                covectorBasisDecomposition = covectorBasisDecomposition,
            ),
        )
    }

    fun VectorSpaceBasisDecomposition.Result<Number, UInt>.sum(): Number =
        (0u ..< dimension).toKoneList().sumOf<_, Number> { this[it] }

    val filteredSolutionCases = solutionCases.mapNotNull { solutionCase ->
        val solutionBasisDecompositions = solutionCase.solutions.map { vectorBasis.decompose(it) }
        val nonZeroSumSolutionIndex = solutionBasisDecompositions.firstIndexThat { _, decomposition ->
            decomposition.sum().isNotZero()
        }
        if (nonZeroSumSolutionIndex == solutionCase.solutions.size) return@mapNotNull null
        val nonZeroSumSolution = solutionCase.solutions[nonZeroSumSolutionIndex]
        val nonZeroSumSolutionSum = solutionBasisDecompositions[nonZeroSumSolutionIndex].sum()
        SolutionCase(
            initialCovectors = solutionCase.initialCovectors,
            solutions = solutionCase.solutions.mapIndexed { index, vector ->
                val sum = solutionBasisDecompositions[index].sum()
                if (sum.isNotZero()) vector / sum else (nonZeroSumSolution + vector) / nonZeroSumSolutionSum
            }
        )
    }

    TODO()
}

Второй же пункт получается просто применением функции checkIntersectionOfAreaAndSolutions.

filterEquationSystems.kt

context(
    covectorEquality: Equality<Covector>,
    covectorHashing: Hashing<Covector>,
)
fun <
    Covector,
> filterEquationSystems(
    equationCases: KoneList<KoneList<EquationCase<Covector>>>
): KoneList<Covector> {
    val equationToCases = equationCases.flatten().associateBy(
        keyEquality = KoneList.equality(covectorEquality),
        keyHashing = KoneList.hashing(covectorHashing),
        keySelector = { it.basis },
        valueTransform = { it.initialCovectors },
        resolve = { _, currentValue, newValue -> if (currentValue.size > newValue.size) currentValue else newValue }
    )

    data class SolutionCase(
        val initialCovectors: KoneList<Covector>,
        val solutions: KoneList<Vector>,
    )

    val solutionCases = equationToCases.nodesView.map { (basis, initialCovectors) ->
        SolutionCase(
            initialCovectors = initialCovectors,
            solutions = basis.reducedRowEchelonFormKernel(
                dimension = dimension,
                vectorBasis = vectorBasis,
                covectorBasisDecomposition = covectorBasisDecomposition,
            ),
        )
    }

    fun VectorSpaceBasisDecomposition.Result<Number, UInt>.sum(): Number =
        (0u ..< dimension).toKoneList().sumOf<_, Number> { this[it] }

    val filteredSolutionCases = solutionCases.mapNotNull { solutionCase ->
        val solutionBasisDecompositions = solutionCase.solutions.map { vectorBasis.decompose(it) }
        val nonZeroSumSolutionIndex = solutionBasisDecompositions.firstIndexThat { _, decomposition ->
            decomposition.sum().isNotZero()
        }
        if (nonZeroSumSolutionIndex == solutionCase.solutions.size) return@mapNotNull null
        val nonZeroSumSolution = solutionCase.solutions[nonZeroSumSolutionIndex]
        val nonZeroSumSolutionSum = solutionBasisDecompositions[nonZeroSumSolutionIndex].sum()
        SolutionCase(
            initialCovectors = solutionCase.initialCovectors,
            solutions = solutionCase.solutions.mapIndexed { index, vector ->
                val sum = solutionBasisDecompositions[index].sum()
                if (sum.isNotZero()) vector / sum else (nonZeroSumSolution + vector) / nonZeroSumSolutionSum
            }
        )
    }

    val areaToCheckIntersectionWith = areaToCheckIntersectionWith(dimension = dimension, basis = vectorBasis)

    val checkedSolutionCases = filteredSolutionCases.filter { solutionCase ->
        checkIntersectionOfAreaAndSolutions(
            area = areaToCheckIntersectionWith,
            solutions = solutionCase.solutions.map { PointWrapper(it) },
        )
    }

    return checkedSolutionCases.map { it.initialCovectors }
}

Финал первой версии: собираем все кусочки воедино

Теперь осталось собрать код воедино. У нас есть функция вычисления всех систем-претендентов и функция фильтрации этих систем. Остаётся только объединить их в одной функции и посчитать все возможные размеры.

possibleResultNumbers.kt

context(
    _: Field<Number>,
    _: Hashing<Number>,
    _: Order<Number>,
    _: EuclideanSpaceOverField<Number, Vector, PointWrapper<Vector>>,
    _: Hashing<Vector>,
    _: VectorSpace<*, Covector>,
    _: Hashing<Covector>,
    _: BasisExtender<Covector>,
)
fun <Number, Vector, Covector> possibleResultNumbers(
    dimension: UInt,
    vectorBasis: VectorSpaceBasis.Finite<Number, Vector>,
    covectorBasis: VectorSpaceBasis.Finite<Number, Covector>,
): KoneList<UInt> {
    val equationCases = computePossibleEquationSystems(
        dimension = dimension,
        covectorBasis = covectorBasis,
    )

    val checkedSolutionCases = equationCases.filterEquationSystems(
        dimension = dimension,
        vectorBasis = vectorBasis,
        covectorBasisDecomposition = covectorBasis,
    )

    return checkedSolutionCases
        .mapTo(
            KoneMutableSet.of(
                elementEquality = UInt.context,
                elementHashing = UInt.context,
                elementOrder = UInt.context,
            )
        ) { it.size }
        .sorted()
}

А в главной функции программы просто подставим в эту функцию конкретные типы чисел, векторов и ковекторов, а также имплементации всех аргументов.

main.kt

typealias Number = BigLongRational
typealias Vector = MDList1<Number>

fun main() {
    print("n = ")
    val dimension = readln().toUInt()

    val field: Field<Number> = Number.context
    val numberHashing: Hashing<Number> = Number.context
    val numberOrder: Order<Number> = Number.context
    val euclideanSpace: EuclideanSpaceOverField<Number, Vector, PointWrapper<Vector>> =
        EuclideanSpaceOverField.mdList1(field, dimension)
    val vectorHashing: Hashing<Vector> = MDList1.hashing(Number.context)
    val vectorBasis: VectorSpaceBasis.Finite<Number, Vector> =
        object : VectorSpaceBasis.Finite<Number, Vector> {
            override val size: UInt get() = dimension

            override fun decompose(vector: Vector): VectorSpaceBasisDecomposition.Result<Number, UInt> =
                VectorSpaceBasisDecomposition.Result { index -> vector[index] }

            override fun get(index: UInt): Vector =
                MDList1(dimension) { if (it == index) field.one else field.zero }
        }
    val basisExtender: BasisExtender<Vector> =
        GaussReducedRowEchelonFormBasisExtender(
            field = field,
            vectorSpace = euclideanSpace,
            basisDecomposition = vectorBasis
        )

    context(
        field,
        numberHashing,
        numberOrder,
        euclideanSpace,
        vectorHashing,
        basisExtender
    ) {
        println(possibleResultNumbers(dimension, vectorBasis, vectorBasis))
    }
}

Результаты

Далее описаны, какие результаты выдаёт первая версия программы и как долго она их считает.

$n$	Результат	Время
3	[0]	300-400 ms
4	[0, 1, 2, 3]	~600 ms
5	[0, 1, 2, 3, 4]	6-7 s
6	[0, 1, 2, 3, 4, 5, 6, 7, 10]	~12 m
7	???	>10 h

✨ Оптимизация первая: фильтрация систем уравнений во время их генерирования

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

Начнём с первого. Понятно, что если добавлена некоторая система уравнений с некоторым ЛНМ, то любая другая система с тем же ЛНМ тоже будет сохранена. На деле у нас перебираются $2^{n-1}$ векторов и все $2^{2^{n-1}} - 1$ систем, которые они порождают, будут сохранены. Это сильно ограничивает нашу программу по производительности (если не заставляет её вовсе закончить выполнение с исключением OutOfMemoryError). Что же делать?

На деле ответ на этот пункт простой. Мы знаем, что для каждого ЛНМ есть максимальная по включению система уравнений, порождающая это ЛНМ. И нам нужен именно размер таковой системы! Значит, вместо добавления того же ЛНМ с новой системой давайте просто расширять систему имеющегося ЛНМ.

Рассмотрим функцию computePossibleEquationSystems, которая за это отвечает. На всякий случай продублируем её нынешний код.

computePossibleEquationSystems.kt

data class EquationCase<Covector>(
    val initialCovectors: KoneList<Covector>,
    val basis: KoneList<Covector>,
)

context(_: VectorSpace<*, Covector>, _: BasisExtender<Covector>)
fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
    covectorBasis: VectorSpaceBasis.Finite<*, Covector>,
): KoneList<KoneList<EquationCase<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneArrayGrowableList<EquationCase<Covector>>()
    }

    equationCases[0u].add(
        EquationCase(
            initialCovectors = KoneList.empty(),
            basis = KoneList.empty(),
        )
    )

    for (coefficients in KoneList.of(1, -1).cartesianPower(dimension - 1u)) {
        if (coefficients.all { it == 1 }) continue
        val vector = coefficients.foldIndexed(covectorBasis[0u]) { index, accumulator, element ->
            accumulator + element * covectorBasis[index + 1u]
        }

        for (dim in dimension - 1u downTo 0u) for (case in equationCases[dim].toKoneList()) {
            val newBasis = case.basis.extend(vector)
            if (newBasis.size == dimension) continue
            equationCases[newBasis.size].add(
                EquationCase(
                    initialCovectors = KoneList.build(case.initialCovectors.size + 1u) {
                        addAllFrom(case.initialCovectors)
                        add(vector)
                    },
                    basis = newBasis,
                )
            )
        }
    }

    return equationCases
}

Во-первых, повторимся, что теперь мы хотим, чтобы equationCases (в каждом своём элементе) хранил по ккаждому ЛНМ ровно одну систему, которая будет расширяться. Значит, теперь equationCases имеет тип KoneList<KoneMutableMap<KoneList<Covector>, KoneMutableSet<Covector>>>. Здесь KoneMutableSet<Covector> хранит ковектороа системы, KoneList<Covector> описывает ковектора ЛНМ, KoneMutableMap<KoneList<Covector>, KoneMutableSet<Covector>> описывает все соответствия "ЛНМ → система уравнений" для данной размерности ЛНМ, а всё equationCases соответственно хранит все такие наборы соответствий для каждой размерности от $0$ до $n-1$. И возвращать стоит такую же структуру в неизменяемом виде, т.е. типа KoneList<KoneMap<KoneList<Covector>, KoneSet<Covector>>>.

computePossibleEquationSystems.kt

context(_: VectorSpace<*, Covector>, covectorHashing: Hashing<Covector>, _: BasisExtender<Covector>)
fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
    covectorBasis: VectorSpaceBasis.Finite<*, Covector>,
): KoneList<KoneList<EquationCase<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneMutableMap.of<KoneList<Covector>, KoneMutableSet<Covector>>(
            keyEquality = KoneList.equality(covectorSpace),
            keyHashing = KoneList.hashing(covectorHashing),
        )
    }

    equationCases[0u][KoneList.empty()] =
        KoneMutableSet.of(
            elementEquality = covectorSpace,
            elementHashing = covectorHashing,
        )

    TODO()

    return equationCases
}

Хорошо, теперь будем как обычно перебирать вектора, размерности dim и системы уравнений с ЛНМ размерности dim и расширять ЛНМ данным вектором.

computePossibleEquationSystems.kt

context(covectorSpace: VectorSpace<*, Covector>, covectorHashing: Hashing<Covector>, _: BasisExtender<Covector>)
fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
    covectorBasis: VectorSpaceBasis.Finite<*, Covector>,
): KoneList<KoneMap<KoneList<Covector>, KoneSet<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneMutableMap.of<KoneList<Covector>, KoneMutableSet<Covector>>(
            keyEquality = KoneList.equality(covectorSpace),
            keyHashing = KoneList.hashing(covectorHashing),
        )
    }

    equationCases[0u][KoneList.empty()] =
        KoneMutableSet.of(
            elementEquality = covectorSpace,
            elementHashing = covectorHashing,
        )

    for (coefficients in KoneList.of(1, -1).cartesianPower(dimension - 1u)) {
        if (coefficients.all { it == 1 }) continue
        val vector = coefficients.foldIndexed(covectorBasis[0u]) { index, accumulator, element ->
            accumulator + element * covectorBasis[index + 1u]
        }

        for (dim in 0u ..< dimension) for ((basis, initialCovectors) in equationCases[dim]) {
            if (vector in initialCovectors) continue
            val newBasis = basis.extend(vector)

            TODO()
        }
    }

    return equationCases
}

Дальше поведение отличается от поведения в версии 1. Если ЛНМ не увеличился, то просто стоит добавить перебираемый вектор в систему и идти дальше. Если получилась ЛНМ большей размерности, то нужно поискать в следующей мапе, есть ли там такая ЛНМ. Если она есть, то надо расширить сопоставляемую ей систему вектором и системой, которые мы перебираем. Иначе добавить эти же вектор и систему вместе с новой ЛНМ в мапу.

computePossibleEquationSystems.kt

context(covectorSpace: VectorSpace<*, Covector>, covectorHashing: Hashing<Covector>, _: BasisExtender<Covector>)
fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
    covectorBasis: VectorSpaceBasis.Finite<*, Covector>,
): KoneList<KoneMap<KoneList<Covector>, KoneSet<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneMutableMap.of<KoneList<Covector>, KoneMutableSet<Covector>>(
            keyEquality = KoneList.equality(covectorSpace),
            keyHashing = KoneList.hashing(covectorHashing),
        )
    }

    equationCases[0u][KoneList.empty()] =
        KoneMutableSet.of(
            elementEquality = covectorSpace,
            elementHashing = covectorHashing,
        )

    for (coefficients in KoneList.of(1, -1).cartesianPower(dimension - 1u)) {
        if (coefficients.all { it == 1 }) continue
        val vector = coefficients.foldIndexed(covectorBasis[0u]) { index, accumulator, element ->
            accumulator + element * covectorBasis[index + 1u]
        }

        for (dim in 0u ..< dimension) for ((basis, initialCovectors) in equationCases[dim]) {
            if (vector in initialCovectors) continue
            val newBasis = basis.extend(vector)
            when (newBasis.size) {
                basis.size -> {
                    initialCovectors.add(vector)
                }
                basis.size + 1u -> {
                    if (dim != dimension - 1u)
                        equationCases[dim + 1u].setOrChange(
                            newBasis,
                            {
                                KoneMutableSet.of(
                                    elementEquality = covectorSpace,
                                    elementHashing = covectorHashing,
                                ).apply {
                                    addAllFrom(initialCovectors)
                                    add(vector)
                                }
                            },
                            { oldInitialCovectors ->
                                oldInitialCovectors.add(vector)
                                oldInitialCovectors
                            }
                        )
                }
                else -> error("Unexpected size of new basis")
            }
        }
    }

    return equationCases
}

Теперь вспомним второй пункт: некоторые ЛНМ сначала добавляем, а потом отфильтровываем. Но можно ли их отфильтровать сразу? Ответ, к счастью, — да. Ведь если какое-то ЛНМ не подходит, т.е. его множество решений не пересекается с нужной областью, то и его расширения тоже будут не подходить. Значит, при расширении ЛНМ можно (и нужно) проверить, что она вообще нам подходит, и если не подходит, то просто пойти перебирать дальше. Для этого давайте не вставлять алгоритм как есть, а выразим его в виде аргумента типа следующего интерфейса:

computePossibleEquationSystems.kt

fun interface BasisValidator<Covector> {
    fun validate(basis: KoneList<Covector>): Boolean
}

context(covectorSpace: VectorSpace<*, Covector>, covectorHashing: Hashing<Covector>, _: BasisExtender<Covector>)
fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
    covectorBasis: VectorSpaceBasis.Finite<*, Covector>,
    basisValidator: BasisValidator<Covector>,
): KoneList<KoneMap<KoneList<Covector>, KoneSet<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneMutableMap.of<KoneList<Covector>, KoneMutableSet<Covector>>(
            keyEquality = KoneList.equality(covectorSpace),
            keyHashing = KoneList.hashing(covectorHashing),
        )
    }

    equationCases[0u][KoneList.empty()] =
        KoneMutableSet.of(
            elementEquality = covectorSpace,
            elementHashing = covectorHashing,
        )

    for (coefficients in KoneList.of(1, -1).cartesianPower(dimension - 1u)) {
        if (coefficients.all { it == 1 }) continue
        val vector = coefficients.foldIndexed(covectorBasis[0u]) { index, accumulator, element ->
            accumulator + element * covectorBasis[index + 1u]
        }

        for (dim in 0u ..< dimension) for ((basis, initialCovectors) in equationCases[dim]) {
            if (vector in initialCovectors) continue
            val newBasis = basis.extend(vector)
            when (newBasis.size) {
                basis.size -> {
                    initialCovectors.add(vector)
                }
                basis.size + 1u -> {
                    if (dim != dimension - 1u && basisValidator.validate(newBasis))
                        equationCases[dim + 1u].setOrChange(
                            newBasis,
                            {
                                KoneMutableSet.of(
                                    elementEquality = covectorSpace,
                                    elementHashing = covectorHashing,
                                ).apply {
                                    addAllFrom(initialCovectors)
                                    add(vector)
                                }
                            },
                            { oldInitialCovectors ->
                                oldInitialCovectors.add(vector)
                                oldInitialCovectors
                            }
                        )
                }
                else -> error("Unexpected size of new basis")
            }
        }
    }

    return equationCases
}

Теперь изменив эту функцию, удалим filterEquationSystems и исправим possibleResultNumbers.

possibleResultNumbers.kt

context(
    _: Field<Number>,
    _: Hashing<Number>,
    _: Order<Number>,
    _: EuclideanSpaceOverField<Number, Vector, PointWrapper<Vector>>,
    _: Hashing<Vector>,
    _: VectorSpace<*, Covector>,
    _: Hashing<Covector>,
    _: BasisExtender<Covector>,
)
fun <Number, Vector, Covector> possibleResultNumbers(
    dimension: UInt,
    vectorBasis: VectorSpaceBasis.Finite<Number, Vector>,
    covectorBasis: VectorSpaceBasis.Finite<Number, Covector>,
): KoneList<UInt> {
    val areaToCheckIntersectionWith = areaToCheckIntersectionWith(dimension = dimension, basis = vectorBasis)

    fun VectorSpaceBasisDecomposition.Result<Number, UInt>.sum(): Number =
        (0u ..< dimension).toKoneList().sumOf<_, Number> { this[it] }

    val equationCases = computePossibleEquationSystems(
        dimension = dimension,
        covectorBasis = covectorBasis,
        basisValidator = BasisValidator { basis ->
            val solutions = basis.reducedRowEchelonFormKernel(dimension, vectorBasis, covectorBasis)
            val solutionBasisDecompositions = solutions.map { vectorBasis.decompose(it) }
            val nonZeroSumSolutionIndex = solutionBasisDecompositions.firstIndexThat { _, decomposition ->
                (0u ..< dimension).toKoneList().sumOf<_, Number> { decomposition[it] }.isNotZero()
            }
            if (nonZeroSumSolutionIndex == solutions.size) return@BasisValidator false
            val nonZeroSumSolution = solutions[nonZeroSumSolutionIndex]
            val nonZeroSumSolutionSum = solutionBasisDecompositions[nonZeroSumSolutionIndex].sum()

            checkIntersectionOfAreaAndSolutions(
                area = areaToCheckIntersectionWith,
                solutions = solutions.mapIndexed { index, vector ->
                    val sum = solutionBasisDecompositions[index].sum()
                    PointWrapper(if (sum.isNotZero()) vector / sum else (nonZeroSumSolution + vector) / nonZeroSumSolutionSum)
                },
            )
        }
    )

    val result =
        KoneMutableSet.of(
            elementEquality = UInt.context,
            elementHashing = UInt.context,
            elementOrder = UInt.context,
        )

    for (dimensionEquations in equationCases) for (equation in dimensionEquations.valuesView) result.add(equation.size)

    return result.sorted()
}

Функция main не поменялась.

Результаты

Далее описаны, какие результаты выдаёт вторая версия программы и как долго она их считает.

$n$	Результат	Время
3	[0]	300-400 ms
4	[0, 1, 2, 3]	~500 ms
5	[0, 1, 2, 3, 4]	4-5 s
6	[0, 1, 2, 3, 4, 5, 6, 7, 10]	~5 m
7	???	>10 h

✨ Оптимизация вторая: объединение случаев ЛНМ, отличных перестановкой координат

В этом разделе заметим, что получающиеся максимальные по включению системы можно переводить друг в друга перестановками координат. Следовательно, можно попробовать не хранить все такие системы, а объединять те, которые отличаются перестановкой координат. Т.е., получая новое ЛНМ, нам нужно попробовать найти одно из старых ЛНМ, которое получается перестановкой координат и снова приведением к приведённой ступенчатой форме. При нахождении старого, добавить в его систему ковекторов новый ковектор с таким же способом переставленными координатами. Иначе, добавить новое ЛНМ с новой системой уравнений.

Рассмотрим функцию computePossibleEquationSystems, которая за это отвечает. На всякий случай продублируем её нынешний код.

computePossibleEquationSystems.kt

context(covectorSpace: VectorSpace<*, Covector>, covectorHashing: Hashing<Covector>, _: BasisExtender<Covector>)
fun <
    Covector,
> computePossibleEquationSystems(
    dimension: UInt,
    covectorBasis: VectorSpaceBasis.Finite<*, Covector>,
    basisValidator: BasisValidator<Covector>,
): KoneList<KoneMap<KoneList<Covector>, KoneSet<Covector>>> {
    val equationCases = KoneList(dimension) {
        KoneMutableMap.of<KoneList<Covector>, KoneMutableSet<Covector>>(
            keyEquality = KoneList.equality(covectorSpace),
            keyHashing = KoneList.hashing(covectorHashing),
        )
    }

    equationCases[0u][KoneList.empty()] =
        KoneMutableSet.of(
            elementEquality = covectorSpace,
            elementHashing = covectorHashing,
        )

    for (coefficients in KoneList.of(1, -1).cartesianPower(dimension - 1u)) {
        if (coefficients.all { it == 1 }) continue
        val vector = coefficients.foldIndexed(covectorBasis[0u]) { index, accumulator, element ->
            accumulator + element * covectorBasis[index + 1u]
        }

        for (dim in 0u ..< dimension) for ((basis, initialCovectors) in equationCases[dim]) {
            if (vector in initialCovectors) continue
            val newBasis = basis.extend(vector)
            when (newBasis.size) {
                basis.size -> {
                    initialCovectors.add(vector)
                }
                basis.size + 1u -> {
                    if (dim != dimension - 1u && basisValidator.validate(newBasis))
                        equationCases[dim + 1u].setOrChange(
                            newBasis,
                            {
                                KoneMutableSet.of(
                                    elementEquality = covectorSpace,
                                    elementHashing = covectorHashing,
                                ).apply {
                                    addAllFrom(initialCovectors)
                                    add(vector)
                                }
                            },
                            { oldInitialCovectors ->
                                oldInitialCovectors.add(vector)
                                oldInitialCovectors
                            }
                        )
                }
                else -> error("Unexpected size of new basis")
            }
        }
    }

    return equationCases
}