Меня часто спрашивают о том, как можно вывести огромную формулу в несколько тысяч слагаемых, которая является решением той или иной трудной задачи. Как вообще действовать, чтобы вывести формулы длиной свыше 10000 слагаемых, о которых рассказывается, например, здесь? На самом деле, во-первых, все эти формулы не выводятся, а «угадываются», а во-вторых, в ряде случаев это делается тривиальным образом. Хотя бывают и чрезвычайно сложные ситуации. Здесь я расскажу о самых простых примерах, так сказать, для начала. Речь пойдёт пока только о линейных рекуррентных соотношениях с постоянными коэффициентами.

Прелюдия

В детстве наверняка вам нравились задачи следующего вида: задана последовательность, угадайте следующее число. Например,

0, 1, 3, 11, 37, 124, ???

Какой ответ? На самом деле, голову ломать не нужно, ответ, очевидно, 0. И попробуйте сказать, что я не прав : ) Ведь нам ничего не сказали про природу этих чисел, не дали никаких начальных знаний, в рамках которых мы должны действовать. Поэтому любое число вместо знаков вопроса будет правильным ответом. И проще всего, мне кажется, сказать, что это 0.

Другое дело, если задача творческая, и нужно придумать максимально простое и красивое объяснение, которое охватывало бы все числа в последовательности сразу. Например, если я напишу, 1,1,2,3,5,8,13, то очевидно, что дальше будет 21, так как все узнали здесь числа Фибоначчи, и это вполне красиво и просто описывает всю последовательность. Если хотите, можете подумать над задачей выше, там есть красивое в этом смысле решение.

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

Вывод формулы из первых принципов

Иногда можно вывести формулу исходя из постановки задачи. Например, если вам дали задание посчитать количество различных последовательностей из 0 и 1, имеющих длину n, то, рассуждая сначала, вы положите n=1 и получите очевидный ответ a1=2, потом скажите, что при увеличении n на 1 мы в новую позицию ставим либо 0, либо 1, тем самым удваивая число вариантов. Получается формула

an= 2an-1    для n>1,

из которой к тому же сразу видно, что an=2n, да притом это оказывается верным для n=1 и, что совсем хорошо, получилось удобное совпадение в предельном случае n=0 (если положить по определению a0=1). Вот так мы вывели формулу.

Точно также можно вывести число способов замостить шахматную доску размером 2×n костяшками домино размером 1×2. Получим числа Фибоначчи.

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

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

Вывести нельзя угадать!

В этом подзаголовке нужно правильно поставить запятую. Если формулу можно вывести и это дело «2 минут», то ради Бога. А если вывести из первых принципов все-таки нельзя? – можно попытаться угадать, а потом доказать, что угадали правильно, то есть именно то, что нужно. И парадокс заключается в том, что в ряде случаев проще доказать, что угадали формулу правильно, чем доказать, что вы учли все возможные случаи при последовательном выводе.

Проще доказать что-либо, когда уже знаешь ответ. Но с другой стороны, угадать можно только тогда, когда знаешь, что в этой формуле может быть, а чего там быть не должно.

К делу!

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

(Все cj — константы). Это абсолютно все задачи, которые напрямую решаются методом матрицы переноса (Transfer Matrix Method, описанный мной, например, здесь) и ещё кое-какие задачи, типа 6 ферзей.

Почему метод матрицы переноса гарантирует существование такого решения? Если вы возьмете матрицу T (порядка K) и будете возводить её в степень Tn для n=0,1,2,…, то, начиная с некоторого номера N, матрица TN будет выражаться через предыдущие степени. Это прямо следует из теоремы Гамильтона – Кэли, которая гласит, что квадратная матрица T обнуляет свой характеристический полином

То есть если вместо z в этот полином подставить матрицу T, то получите матрицу из нулей:

А это означает, что старшую степень матрицы можно выразить через младшие:

Домножая обе части выражения на матрицу T сколько угодно раз, мы получим, что любая степень Tn (для n≥K) будет выражаться через предыдущие Tn-1, Tn-2, … . Иными словами, (i,j)-й элемент n-й степени матрицы переноса T (назовём его an) будет удовлетворять рекуррентному соотношению

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

Так как же его угадать-то?

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

А если не работает, то техника угадывания всё равно простая. Предположим, для начала, что порядок соотношения N вам известен. Тогда составим систему линейных уравнений из первых 2N чисел последовательности (a1, a2, …, a2N), решением которой будет искомый набор рациональных коэффициентов ci (i=1,2,…,N):

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

Как решать?

Решать такую систему нужно методом Диксона. Дело в том, что этот метод для данной задачи хорошо подходит. Он ищёт ответ в виде P-адической аппроксимации, и точность этой аппроксимации должна быть тем выше, чем больше суммарная длина числителя и знаменателя самой длинной в этом смысле дроби в ответе. В реальных задачах очень часто длина коэффициентов значительно меньше длины чисел в исходной последовательности. Например, для задачи о количестве гамильтоновых циклов на решётке размером 10×n сами числа могут превышать по размеру 1000 десятичных знаков (для n=700), тогда как в соотношении 346-го порядка для этой задачи коэффициенты целые и максимальная длина всего 70 знаков. Метод Диксона для этой задачи сходится всего за 16 итераций самого сложного шага Lifting. Иначе говоря, работает почти мгновенно.

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

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

Кстати, тот факт, что матрица имеет Теплицеву структуру наводит на подозрение, что задача может быть решена ещё более эффективно. Я пробовал запускать алгоритм Берлекампа-Месси, который решает Теплицевы системы, записанные в форме авторегрессионного фильтра, но этот алгоритм не справился даже с системой порядка N=2086 за 5 лет процессорного времени. А других алгоритмов в этой области я не знаю… Нужны алгоритмы, которые работают с целыми / рациональными числами. А алгоритм Б-М был рассчитан на поиск приближённого решения в формате вещественных чисел, поэтому неудивительно, что при попытках формально переписать его для целых чисел, промежуточные результаты выглядели в виде дробей, числитель и знаменатель которых занимал несколько миллионов бит.

Если порядок неизвестен

Чаще всего порядок N искомого соотношения неизвестен. Но, исходя из постановки задачи, часто можно сказать, что порядок не превышает заранее известную величину K. Например, если вы решаете задачу методом матрицы переноса, то порядок этой матрицы и будет этим самым числом K. Нужно взять 2K первых чисел из последовательности: a1, a2, …, a2K и составить матрицу системы

Далее нужно отыскать её ранг. Но не честным способом (так вы его просто не найдёте), а по модулю какого-нибудь достаточно большого простого числа. Например, возьмите P=231-1. Это называется «хак»: мы не знаем, будет ли ранг найден точно, но вероятность того, что вы ошибётесь, обратно пропорциональна числу P. Так вот, что показывает ранг? Ранг – это максимальное количество линейно независимых столбцов в матрице. Поскольку наше рекуррентное соотношение является верным для всех n>N (где N — порядок, который нам ещё неизвестен, он не больше ранга матрицы), то любые подряд идущие N столбцов будут линейно независимыми. Найдя ранг, который будет равен N, мы можем из 2N первых чисел составить систему и решить её как в случае с известным порядком N. Но здесь вы должны точно знать, что взятое вами K больше, чем искомый порядок N.

Пример

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

1, 1, 2, 3, 4, 5, 6, 7, …

Предположим, что по каким-то соображениям, данным свыше, нам стало известно, что K=4 (порядок будущего соотношения не выше 4). Составляем матрицу системы:

Ранг этой матрицы равен 3. Значит предполагаемый порядок рекуррентного соотношения равен N=3. Значит нам нужно составить систему из первых 6 чисел:

Её решением будет c3=0, c2=−1, c1=2. Значит рекуррентное соотношение будет иметь вид

an=2an-1 — an-2 + 0⋅ an-3, n≥4;

Формально всё правильно, однако если убрать слагаемое с 0, то соотношение станет иметь порядок 2, а не 3. Однако «правильно» работать оно будет по-прежнему только с n=4. То, что оно уменьшило свой порядок таким хитрым образом не означает, что оно неожиданно станет работать для тех значений n, для которых мы его не выводили. Оно как работало для n≥4, так и работает только для этих n. То есть тут надо быть очень осторожным и всегда указывать, для каких n ваше рекуррентное соотношение начинает работать.

Проверка

После того, как соотношение найдено, нужно его проверить. Если вы решали задачу методом матрицы переноса, порядок которой равен K, а порядок соотношения получился равным N<K, то нужно прогнать через последовательность элементы с 2N+1 по 2K. Если они удовлетворили рекуррентному соотношению, то оно гарантированно будет правильным. Я уверен, что вы сами сможете объяснить, почему это так. На самом деле вы также сможете показать, что не обязательно брать элементы до 2K, а можно взять меньше. Но мы как бы перестраховались, это особенно полезно тогда, когда вы не хотите строго доказывать что-то на этот счёт.

Обсуждение статьи на форуме в этой теме.