Матрица в степени: рекурренты и маршруты заданной длины

Если умеешь построить матрицу перехода и проверить ориентацию столбца состояния, листай до «Теперь сам».

Как думаешь, нужно ли вычислять все промежуточные значения

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

На что это похоже

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

Механизм: состояние достаточно для следующего шага

Для F₀ = 0, F₁ = 1, Fₙ₊₂ = Fₙ₊₁ + Fₙ возьмём столбец (Fₙ₊₁, Fₙ). Матрица [[1,1],[1,0]] переводит его в (Fₙ₊₂, Fₙ₊₁). Первая строка складывает оба элемента, вторая переносит прежний верхний вниз. Значит после n применений начальный столбец (1,0) превращается в (Fₙ₊₁,Fₙ).

Для рекурренты fₜ = c₁fₜ₋₁ + ... + cₖfₜ₋ₖ верхняя строка матрицы состоит из коэффициентов. Ниже стоят единицы, сдвигающие память, и нули в остальных клетках. Именно хранение последних k значений делает переход одинаковым для любого номера шага. Если появляется постоянное слагаемое, к состоянию добавляют координату, всегда равную единице.

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

Правило

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

Разбор: та же рекуррента двумя способами

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

MOD = 1009


def multiply(a, b):
    rows, middle, cols = len(a), len(b), len(b[0])
    result = [[0] * cols for _ in range(rows)]
    for i in range(rows):
        for t in range(middle):
            for j in range(cols):
                result[i][j] += a[i][t] * b[t][j]
                result[i][j] %= MOD
    return result


def power(a, n):
    size = len(a)
    result = [[int(i == j) for j in range(size)]
              for i in range(size)]
    while n:
        if n & 1:
            result = multiply(result, a)
        a = multiply(a, a)
        n //= 2
    return result


def fibonacci(n):
    return power([[1, 1], [1, 0]], n)[1][0]


x, y = 0, 1
for n in range(30):
    assert fibonacci(n) == x % MOD
    x, y = y, x + y
print("F10:", fibonacci(10))
print("F1000000 mod 1009:", fibonacci(1_000_000))
a = [[0, 1, 1], [1, 0, 1], [1, 1, 0]]
print("Маршрутов из 0 обратно в 0 за 2 ребра:", power(a, 2)[0][0])
print("Нулевая степень:", power(a, 0))

F10: 55
F1000000 mod 1009: 153
Маршрутов из 0 обратно в 0 за 2 ребра: 2
Нулевая степень: [[1, 0, 0], [0, 1, 0], [0, 0, 1]]

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

Разбор: маршрут допускает возвращение

Матрица смежности A[i][j] равна числу рёбер из i в j. В произведении A²[i][j] = Σ A[i][t]·A[t][j] произведение выбирает первое и второе ребро, сумма перебирает промежуточную вершину. По индукции Aᴸ считает маршруты из L рёбер. В показанном треугольнике из вершины ноль обратно в ноль за два ребра ведут два маршрута: через каждую другую вершину.

Это не подсчёт простых путей. Для запрета повторов пришлось бы помнить множество посещённых вершин; одного номера текущей вершины уже недостаточно. При L = 0 единичная матрица считает пустой маршрут из вершины в неё саму.

Что ломается без этого

Если перепутать строковое и столбцовое состояние, коэффициенты окажутся транспонированными. Программа даст числа, но не нужную последовательность. Проверка одного перехода дешевле отладки миллионной степени. Другая ошибка - назвать оценку O(log n), умолчав о размере матрицы: при растущем k кубическое умножение становится главной затратой.

Зачем это в работе

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

Теперь сам

Последовательность задана gₜ₊₂ = 2gₜ₊₁ + 3gₜ, g₀ = 1, g₁ = 2. Выпиши матрицу для столбца (gₜ₊₁,gₜ) и найди g₃.