Symbolic Integration

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

Классы функций одной переменной

Чем шире класс функции, тем сложнее брать интеграл. Например, интегрирование полиномов одной переменной

p(x) = an xn + ... + a1 x + a0

не составляет никакого труда. Уже некоторую сложность представляет интегрирование рациональной функции

r(x) = an xn + ... + a1 x + a0 bm xm + ... + b1 x + b0 .

Если степень числителя n меньше степени знаменателя m, и коэффициенты лежат в поле вещественных чисел, то интегрирование прозводится следующим образом: знаменатель раскладывается на множители вида (x - x_i)^k_i и (x² + b_j x + c_j)^l_j, где каждый x_i — корень кратности k_i, а b_j и c_j задают пару сопряжённых комплексных корней кратности l_j. Затем методом неопределённых коэффициентов дробь r(x) раскладывается на сумму конечного количества дробей вида

Ai (x-xi)ki и Bj x+Cj (x2 +bjx + cj)lj

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

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

Кажется, что интегрирование рациональных функций — уже чисто техническая задача? Это не совсем так. Если степень многочлена больше 4, то его корни в общем случае не записываются в радикалах. Это доказал Абель в 1826 г. Например, вещественный корень уравнения «x⁵ + x + a = 0» невозможно записать, пользуясь арифметическими операциями и операцией извлечения корня.¹

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

Алгебраическая функция y от переменной x задаётся неявным образом, через полином от x и y:

pn(x) yn + ... + p1(x) y + p0(x) = 0

Для примера преобразуем указанное выше уравнение 5-й степени в алгебраическую функцию y от переменной x:

y5+y+x=0

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

К классу элементарных функций относятся алгебраические функции, экспонента, логарифм и любые их конечные композиции. Подразумевается оперирование над полем комплексных чисел, что означает, что элементарны и все тригонометрические функции, т.к. они выражаются через комбинации e^ix и e^-ix. Обратные тригонометрические функции, соответственно, выражаются через логарифмы. Показательная функция a^x выражается через экспоненту и логарифм: a^x = e^{x ln(a)}.

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

Интегрирование в теории

С рациональными функциями, несмотря на указанные выше неприятности, принципиальных трудностей нет: интеграл всегда существует в виде элементарной функции, и есть процедура, которая его гарантированно найдёт.

С алгебраическими функциями дело обстоит сложнее.

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

Это даёт следующую общую запись для интегрируемых алгебраических функций:

f(x) = v0'(x) + m Σ i=1 ci vi'(x) vi(x)

(1)

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

Гипотезу Лапласа доказал Лиувилль в 1833 г., а ещё через два года он опубликовал доказательство аналогичного факта для некоторого обобщения класса алгебраических функций, т.н. элементарных функций по Лиувиллю. Они имеют форму алгебраических функций, но в роли переменных в них не только x, но ещё конечный набор логарифмов и экспонент от алгебраических функий от x:

y(x) = f(x, z₁(x), ... z_r(x)) — алгебраическая функция от x, z₁(x), ... z_r(x)

z_i(x) = ln(q_i(x)) или exp(q_i(x)), где q_i — алгебраическая функция от x.

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

⌠ ⌡ f(x, z1(x), ... zr(x)) dx = y0(x, z1(x), ... zr(x)) + m Σ i=1 Ai ln(yi(x, z1(x), ... zr(x)))

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

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

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

Итак, дифференциальная теория Галуа. Ключевое понятие в ней — это т.н. элементарное расширение поля функций², которое получается последовательным добавлением в него логарифмов, экспонент и алгебраических функций:

• ln(f(x)), где f ∈ K
• exp(f(x)), где f ∈ K
• алгебраической функции, полученной уравнением вида P(f₁(x), ... f_n(x))=0, где f₁, ... f_n ∈ K.

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

Ключевая теорема дифференциальной теории Галуа является обобщением теоремы Лиувилля и говорит о том, что если функция f из поля F имеет первообразную g в поле G — элементарном расширении F, то f записывается в форме (1), где v_i ∈ F, а g соответственно содержит функцию из F и логарифмы функций из F с весами.

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

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

Интегрирование на практике

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

Первая реализация алгоритма Риша была выполнена Джоэлом Мозесом в рамках знаменитого «Project MAC» в MIT в 1971 г. Программа под названием SIN интегрировала чисто трансцендентные функции.

Дэвенпорт в 1981 г., основываясь на работе Риша и некоторых глубоких результатах дифференциальной алгебры и комплексного анализа, разработал алгоритм интегрирования чисто алгебраических функций и реализовал его в известной среде символьных вычислений REDUCE-2.

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

Алгоритм Дэвенпорта ещё и обладал большой вычислительной сложностью, что его автор признавал и выражал надежду, что алгоритм можно упростить. Так и произошло: Барри Трагер из MIT в 1984 г. внёс серьёзные улучшения в алгоритм Дэвенпорта. Обновлённый алгоритм обладал гораздо большей практической ценностью и был реализован в математических программах Axiom и Maple.

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

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

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

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

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

x √ x4 + 10 x2 - 96 x - 71

равна

- 1 8 ln((x6+15 x4-80 x3+27 x2-528 x+781) √ x4+10 x2-96 x-71

- x8 - 20 x6 + 128 x5 - 54 x4 + 1408 x3 - 3124 x2 - 10001),

однако в таком виде её находит только Axiom. Прочие пакеты — Mathematica, Maple, Matlab, Maxima, Reduce — в лучшем случае выдают громоздкий ответ с участием эллиптических интегралов, т.е. неспособны выдать ответ в виде элементарной функции.³

Первообразную функции e^arcsin(x), равную ½(√ 1-x² +x)e^arcsin(x), способны найти Mathematica, Axiom и Reduce.

Неопределённый интеграл ∫√arctan(x) dx невозможно найти (или определить его отсутствие) ни в одном из перечисленных пакетов; при этом Axiom выводит наиболее осмысленное сообщение: implementation incomplete (constant residues).

Нерешённые вопросы

Полная имплементация алгоритма Риша-Трагера-Бронштейна ждёт своего автора. Однако, и это ещё не всё: в теории символьного интегрирования тоже остались открытые вопросы.

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

Дело в том, что в процессе работы алгоритм Риша выстраивает «башню» элементарных расширений поля рациональных функций, с тем чтобы «верхушка» башни содержала подынтегральное выражение. При этом не всегда понятно, готова ли уже «башня» или нужно добавить в неё ещё один логарифм, экспоненту, или алгебраическую функцию. Например, если в подынтегральном выражении встретилась экспонента e^x, она добавляется в поле (происходит элементарное расширение поля функцией e^x). Если при дальнейшем анализе в подынтегральном выражении встречается e^x2, то эту функцию также следует добавить в поле, вырастив таким образом «башню» до двух этажей. Но если после этого попадётся функция e^x+x2, то алгоритм должен определить, что она принадлежит сформированному на данный момент элементарному расширению, и не наращивать «башню». Да, в описанном случае это несложно сделать; чуть сложнее определить, например, на языке элементарных функций, что arctan(x) = -arctan(1/x). Но в общем случае эта задача очень сложна⁴, и нынешние реализации алгоритма Риша используют лишь эвристики для её решения, которые по определению ненадёжны.

Как указывает Дэвенпорт и следом за ним Бронштейн, при реализации алгоритма можно «сделать вид», что нам известен т.н. базис трансцендентности для поля констант, то есть например считать e и √ 2 алгебраически независимыми. Это значит, что мы принимаем на веру, что нет никакого полинома от e и √ 2 , который был бы равен 0. Такое решение нужно принимать о каждой новой константе, появляющейся в ходе работы алгоритма, которую эвристика не позволяет связать алгебраическим выражением с ранее учтёнными константами. Если полученный базис трансцендентности окажется неверным, это может отразиться на корректности выданного ответа в том случае, если ответ гласит, что «интеграл не является элементарной функцией» (да, это касается и Axiom, о которой выше было сказано противоположное). Если же интеграл найден, его можно считать правильным.

Известным частным случаем задачи об алгебраической независимости является т.н. Constant Problem: определение того, тождественно ли нулю некоторое математическое выражение. Задача тривиально решается для многочленов, но для более широкого класса функций, как показал Даниэль Ричардсон в 1968 г., эта задача неразрешима.

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

Если есть функция f, получающаяся с помощью сложения, умножения и композиции из синуса, экспоненты и модуля, а коэффициенты f используют рациональные числа, π и ln(2), то задача проверки утверждения «f тождественно равна нулю» алгоритмически неразрешима.

Функция модуля существенна в доказательстве Ричардсона, и есть основания полагать (да и сам автор теоремы надеется), что без неё задача Constant Problem всё таки разрешима, что даёт надежду на создание абсолютно честного алгоритма символьного интегрирования.

И последнее: класс элементарных функций кажется несколько надуманным, не так ли? Например, функция ошибки

erf(x) = 2 √ π ⌠ ⌡ x 0 e-t2 dt

не является элементарной, как доказал Лиувилль. Но чем она хуже экспоненты и логарифма, которые мы относим к элементарным функциям? На компьютере она может быть вычислена с произвольной точностью, и даже, как экспонента и логарифм, обладает элементарной производной. Эту функцию вполне можно включить в «надэлементарные» расширения функций, так же как и эллиптические интегралы, ∫x^xdx и многие другие. Тема интегрируемости в «надэлементарных» функциях была затронута ещё автором SIN в 1970-х годах, и ныне представляет обширное поле для исследований.

Благодарности

Автор признателен А. В. Демьянову, П. А. Мозоляко и Н. Н. Васильеву за вычитку черновиков этой статьи и высказанные замечания.

Литература

J. Liouville. Premier memoire sur la determination des integrales dont la valeur est algebrique. Journal de l’Ecole Polytechnique, 14:124–148, 1833

J. Liouville. Second memoire sur la détermination des integrales dont la valeur est algebrique. Journal de l’Ecole Polytechnique, 14:149–193, 1833

J. Liouville. Mèmoire sur l’intégration d’une classe de fonctions transcendantes. J. Reine Angew. Math. Bd. 13, p. 93-118. (1835)

G. H. Hardy. The Integration of Functions of a Single Variable, 1916.

R. H. Risch (1970). The solution of the Problem of Integration in Finite Terms. Bulletin of the American Mathematical Society 76 (3): 605–608.

Joel Moses. Symbolic integration the stormy decade. Proceedings of the second ACM symposium on Symbolic and algebraic manipulation (1971) 427-440.

Дж. Дэвенпорт. Интегрирование алгебраических функций, Москва «Мир» 1985.

B. Trager. Integration of algebraic functions. PhD Thesis, MIT, 1984

M. Bronstein. Integration of elementary functions. Journal of Symbolic Computation, Volume 9, Issue 2 (1990). 117-173

M. Bronstein. Symbolic Integration Tutorial. ISSAC'98, Rostock.

Richardson, D. Some Unsolvable Problems Involving Elementary Functions of a Real Variable. J. Symbolic Logic 33, 514-520, 1968.

---

Сноски:

¹Кстати, это число называется корнем Бринга Br(a), и с его помощью можно записать все комплексные корни более общего уравнения «x⁵ + px + q = 0».

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

³Проверялись следующие версии программ: Axiom 2010, Mathematica 7, Maple 13, Matlab R2007, Maxima 5.21, Reduce 2009.

⁴Например, до сих пор никто не знает, является ли e + π рациональным числом (хотя известно, что e и π по отдельности — трансцентентные числа). Известно только равенство e^{i π} = -1 и бесконечная дробь Рамануджана, связывающая e и π.

Mark Tarver и Qi

Qi (произносится «чи») — функциональный язык программирования, с проверщиком типов полным по Тьюрингу, как и в Haskell и C++. Автор (Марк Тарвер, профессор в отставке) при этом утверждает, что его творение затмевает все существующие разработки и останется непревзойдённым навсегда, по причине использования наиболее гибкой нотации для описания типов. К тому же, оно работает поверх Лиспа и разрешает использовать функции на Лиспе.

О практическом применении языка речи не идёт; на вопрос, а что такого есть в мире написанного на Qi, автор гордо отвечает, что сам Qi на себе и написан, не считая вспомогательных библиотек к нему же. Скорее, Qi — исследовательский инструмент для разработки новых систем типов; ну и для саморазвития, конечно.

По части последнего Марк вообще увлечённая личность: в августе прошлого года помпезно попрощался с Qi и вообще с программированием, и уехал в Индию заниматься йогой. В ноябре посетовал, что-де в ашраме туго с интернетом, а в декабре (уже не с такой помпой) оповестил сообщество, что «phase of travels is over», вернулся в Англию и продолжает развивать Qi.

Как тут не пожелать доктору Тарверу всяческих успехов.

Indigo

Мы открыли сайт с программными продуктами для органической химии, которые разрабатывались с 2004 года. Продукты реализованы на C++ и объединены в тулкит под названием Indigo.

• Bingo: картридж (плагин) к СУБД Oracle с реализацией различных видов «химического» поиска молекул и реакций.
• Dingo: библиотека для рендеринга молекул и реакций, с обёрткой на C# и command-line утилитой.
  
• Cano: библиотека для подсчёта canonical SMILES — уникального представления молекулы, с обёрткой на C#.
  
• Deco: библиотека и command-line утилита для разложения набора молекул на общую часть и заместители (R-Group deconvolution)
  
• Nucleo: библиотека для работы с нуклеотидными цепочками, с обёрткой на Java.

Поддерживаемые платформы — Windows, Linux, Mac OS X, Solaris на SPARC.

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

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

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

Мы открыли три публичные почтовые рассылки:

• http://groups.google.com/group/indigo-dev — для разработчиков
• http://groups.google.com/group/indigo-bugs — для баг-репортов
• http://groups.google.com/group/indigo-general — для всего остального
  

Ещё раз адрес сайта: http://ggasoftware.com/opensource/indigo/.

Команда разработчиков находится в Санкт-Петербурге. Будем очень рады любым комментариям, вопросам, сообщениям об ошибках.

Graph automorphisms, canonical labeling, nauty

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

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

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

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

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

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

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

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

Любой канонический код состоит из (i) определённой схемы кодирования графа и (ii) процедуры переупорядочения вершин. Нас не интересует схема кодирования; это может быть обход в глубину, или в ширину, или развёрнутая в строку матрица смежности, или просто списки вершин и рёбер... как бы то ни было, кодирование графа подразумевает какой-то порядок его вершин. Наша задача — упорядочить вершины так, чтобы результат не зависел от изначального порядка. Это называется каноническая нумерация (canonical numbering, canonical labeling).

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

Реалии, однако, мало изменились за последние 15 лет. Австралийский профессор Брендан МакКей (Brendan McKay) в 1980-х годах написал статью «Practical Graph Isomorphism» и соответствующий программный пакет под названием nauty. Долгое время никому не удалось превзойти nauty по быстродействию; хотя для некоторых частных случаев более эффективные алгоритмы всё же появляются, например вот, вот. Особого упоминания заслуживают библиотеки bliss (2007) и saucy2 (2008), которые на больших разреженных графах уверенно обгоняют nauty. Но эти алгоритмы, фактически, не более чем nauty с модификациями.

Печально другое: мало кто из программистов сумел хотя бы подойти близко к пониманию того, как работает nauty. Если мы посмотрим на современные программные пакеты, то окажется, что они либо напрямую используют эту библиотеку: Magma, MuPAD, GRAPE, Planarity, Cyberaide, igraph (использует bliss), либо используют более примитивные алгоритмы: OpenBabel, Marvin. В любом случае, об осмыслении nauty речи не идёт. Приятным исключением является InChI, в которой присутствует алгоритм, явно сделанный с оглядкой на nauty. Настоящая популярность к алгоритму МакКея придёт тогда, когда в каждом пакете он будет реализован по-своему.


Перестановки и упорядоченные разбиения

Говоря здесь о разбиениях, мы подразумеваем упорядоченные разбиения множества вершин графа, т.е. представления этого множества вершин виде упорядоченного набора его попарно непересекающихся непустых подмножеств. [a b c | d e] и [a c b | e d] — одно и то же разбиение, но [d e | a b c] от него отличается. Подмножества [a b c] и [d e] называются ячейками (cells).

Если каждая ячейка разбиения π₁ является подмножеством какой-либо ячейки разбиения π₂, то π₁ называется подразбиением π₂, а π₂ — надразбиением π₁. Также говорят, что π₁ мельче π₂, а π₂ крупнее π₁. Мы подразумеваем, что при подразбиении сохраняется порядок ячеек. Например, [a b | c | d e] является подразбиением [a b c | d e] и [a b | c d e], но не [c | a b d e].

Разбиение, в котором все подмножества состоят из одного элемента, называется дискретным (discrete partition). Дискретное разбиение очевидным образом образует перестановку (permutation) вершин. Например, дискретное разбиение [b | c | a | d | e] образует перестановку (b, c, a, d, e).

Разбиение, состоящее из одной ячейки, называется элементарным (unit partition). Любое разбиение является подразбиением элементарного.

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

function AllPermutations(π = [V₀|...|V_k])
   if π — дискретное разбиение:
      yield π
   else:
      выбрать V_i ∈ π : |V_i| ≠ 1
      for v ∈ V_i:
         AllPermutations([V₀|...|{v}|V_i\v|...|V_k])

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


Годные разбиения

Мы называем разбиение π годным (equitable), если для него выполняется следующее свойство: для любых двух ячеек V₁, V₂ ∈ π (не обязательно различных) и для любых v₁, v₂ ∈ V₁ выполняется равенство d(v₁, V₂)=d(v₂, V₂), где d(v,V) — количество вершин в V, смежных с v. Очевидно, что все дискретные разбиения являются годнымы.

В работе МакКея доказано, что для каждого разбиения π существует единственное его наикрупнейшее подразбиение, являющееся годным. Процедура уточнения (refine) разбиения находит это годное подразбиение. В упрощённом виде она выглядит так:

function Refine(π = [V₀|...|V_k])
   for Vi ∈ π :
      for Vj ∈ π:
         подразбить V_j, сгруппировав элементы v ∈ V_j по d(v, V_i) и отсортировав по возрастанию d(v, V_i)
         заменить V_j в π на результат этого подразбиения, удлиняя цикл по V_i, но не удлиняя цикл по V_j
   return π

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

function AllEquitablePermutations(π = [V₀|...|V_k])
   Refine(π)
   if π — дискретное разбиение:
      yield π
   else:
      выбрать V_i ∈ π : |V_i| ≠ 1
      for v ∈ V_i:
         AllPermutations([V₀|...|{v}|V_i\v|...|V_k])

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


(Правая половина дерева показана не целиком. Вызовы Refine, которые не уточняют подразбиение, также опущены.)


Сравнения перестановок и автоморфизмы

Наша задача — найти каноническую перестановку множества вершин исходного графа. Для этого нужно обзавестить функцией сравнения перестановок cmp, которая бы определяла, которая из двух «лучше» (или они одинаковые). Функция должна задавать транзитивное отношение, а её результат не должен зависеть от первоначальной нумерации вершин графа. Пример такой функции: по исходному графу G и перестановкам π₁ и π₂ построить графы G₁ и G₂, вычислить их матрицы смежности, и в качестве cmp(π₁, π₂) вернуть результат лексикографического сравнения матриц, развёрнутых в строки. Таким образом, канонической перестановкой будет считаться та, для которой матрица смежности наибольшая; если таких перестановок две и более, то они считаются неразличимыми и в качестве канонической можно взять любую из них.

Какой смысл для нас имеют «неразличимые» перестановки? Пусть cmp(π₁, π₂)=0. Несложно доказать, что в этом случае перестановка γ=π₂^-1∘π₁ (композиция π₁ и перестановки, обратной π₂) является автоморфизмом вершин исходного графа.

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

Есть следующий способ найти генераторы группы автоморфизмов графа: при обработке первого встретившегося дискретного разбиения ξ оно сохраняется. Каждое последующее дискретное разибение π сравнивается с ξ: если перестановка ξ^-1∘π является автоморфизмом графа, то эта перестановка записывается в группу генераторов.


Сокращение перебора 1: автоморфизм ξ^-1∘π

После обнаружения автоморфизма ξ^-1∘π, обход можно «поднять» по дереву до общего предка π и ξ, поскольку поддерево этого предка, в котором находится π, «изоморфно» поддереву, в котором находится ξ, и не содержит новой информации. Строгое доказательство приведено в работе МакКея.


Сокращение перебора 2: автоморфизм ρ^-1∘π

Как было сказано, поиск канонической нумерации — только опция в алгоритме нахождения автоморфизмов. Когда эта опция включена, мы храним ещё одну перестановку ρ, которая является «лучшей» на данный момент (best so far, BSF). При обнаружении первой перестановки ξ мы копируем её в ρ. Каждую последующую перестановку π мы сравниваем с ρ. Если cmp(π, ρ) > 0, то мы копируем π в ρ. Если cmp(π, ρ) < 0, то ничего не делаем. Если же cmp(π, ρ) = 0, то ρ^-1∘π является новым автоморфизмом, и обход можно «поднять» по дереву до общего предка π и ρ.


Сокращение перебора 3: орбиты

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


В начале работы nauty каждая вершина принадлежит своей собственной орбите. По мере открывания новых автоморфизмов орбиты увеличиваются, а их количество — сокращается. Например, если найден автоморфизм (a, b, c, d, e) → (b, a, c, e, d), то список орбит из [a][b][c][d][e] становится [a b] [c] [d e]. Если после этого найден автоморфизм (a, b, c, d, e) → (a, c, b, d, e), то список орбит становится [a b c] [d e]

Рассмотрим «линию предков» первого дискретного разбиения ξ. Предположим, мы в процессе обхода дерева вернулись к какому-либо из предков ξ (назовём его ν), чтобы продолжить обход других потомков ν, перебирая другие элементы ячейки V_i. Если рассматриваемый элемент v_k ∈ V_i лежит на одной орбите с уже рассмотренным элементом v_j ∈ V_i, то существует автоморфизм, открытый в процессе обхода более ранних узлов, который переводит v_j в v_k. В работе МакКея доказано, что в этом случае нет смысла рассматривать всё поддерево, произведённое отделением v_k: все листья этого поддерева будут изоморфны ранее открытым листьям, и все автоморфизмы, которые будут найдены при обходе этого поддерева, будут повторять предыдущие.

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


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


Сокращение перебора 4: fix и mcr

С узлами, которые не являются предками ξ, предыдущее сокращение перебора делать нельзя. Для них есть менее сильный, но тоже действенный критерий отсечения потомков. Для применения этого критерия необходимо сохранять информацию о каждом автоморфизме в отдельности (вместо множества орбит, которое накапливает информацию о всех автоморфизмах в совокупности). Для каждого найденного автоморфизма γ запоминаются два множества: fix(γ) и mcr(γ). fix(γ) — это множество вершин, которое при действиии автоморфизма остаются на своих местах. mcr(γ) — это множество вершин, которые являются минимальными по номеру в своих орбитах. Имеются в виду орбиты, составленные по одному автоморфизму γ, а не те, которые используются в предыдущем способе отсечения.

Рассмотрим множество V_i, элементы которого мы перебираем для продолжения обхода дерева из какого-то узла, глубины l. В процессе пути от корня дерева до этого узла l вершин были «зафиксированы», когда мы отделяли их от их ячеек. Назовём множество этих вершин fixed.

В работе МаКея доказано следующее: если fixed ⊂ fix(γ), то можно присвоить в V_i пересечение V_i ∩ mcr(γ) без ущерба для алгоритма. Фактически, речь идёт об автоморфизмах, которые «не затрагивают» вершины, которые мы зафиксировали. Их орбиты можно эксплуатировать так же, как и в предыдущем разделе, что и делается путём пересечения V_i и mcr(γ). Означенная операция выполняется для всех γ, fix и mcr которых мы успели сохранить на данный момент. Следует отметить, что генераторов группы автоморфизмов графа может быть очень много, и сохранять все fix и mcr не всегда возможно. МакКей рекомендует ограничивать количество сохраняемых fix и mcr числом, пропорциональным количеству вершин в графе, с коэффициентом пропорции больше 1 (у него почему-то 50/32).


Алгоритм целиком

В коде имеются глобальные переменные:

• getcanon — опция получения канонической нумерации в ρ
• gca_first — уровень общего предка текущего узла и ξ
• gca_canon — уровень общего предка текущего узла и ρ
• orbits — массив индексов орбит; содержит для каждой вершины номер минимальной вершины в её орбите
• Ψ — списов посчитанных множеств fix и mcr для найденных автоморфизмов

function AutomorphismSearch(π)
   FirstNode(1, π)

function FirstNode(level, π = [V₀|...|V_k])
   π ← Refine(π)
   if π — дискретное:
      ξ ← π
      gca_first ← level
      if getcanon:
         ρ ← π
         gca_canon ← level
      return level-1
   Vi ← первая ячейка π, в которой больше 1 вершины
   for v ∈ Vi
      if orbits[v] ≠ v:
         continue (сокращение 3)
      π'←[V₀|...|{v}|V_i\v|...|V_k]
      fixed ← fixed ∪ v
      if v — первая вершина в V_i:
         rtnlevel ← FirstNode(level+1, π')
         gca_first ← level
      else:
         rtnlevel ← OtherNode(level+1, π')
      fixed ← fixed \ v
      if rtnlevel < level:
         return rtnlevel
      if level < gca_canon:
         gca_canon ← level
   return level-1

function OtherNode(level, π = [V₀|...|V_k])
   π ← Refine(π)
   if π — дискретное:
      if cmp(π,ξ)=0:
         AddFixMcr(ξ^-1∘π)
         JoinOrbits(ξ^-1∘π)
         return gca_first (сокращение 1)
      if getcanon:
         comp ← cmp(π, ρ)
         if comp=0:
            AddFixMcr(ρ^-1∘π)
            JoinOrbits(ρ^-1∘π)
            return gca_canon (сокращение 2)
         if comp>0:
            ρ ← π
      return level-1
   Vi ← первая ячейка π, в которой больше 1 вершины
   for (fix, mcr) ∈ Ψ:
      if fixed ⊂ fix:
         V_i ← V_i ∩ mcr (сокращение 4)
   for v ∈ Vi:
      π'←[V₀|...|{v}|V_i\v|...|V_k]
      fixed ← fixed ∪ v
      rtnlevel ← OtherNode(level+1, π')
      fixed ← fixed \ v
      if rtnlevel < level:
         return rtnlevel
      if level < gca_canon:
         gca_canon ← level
   return level-1


Детали реализации

Авторская имплементация nauty содержит трюки для увеличения быстродействия и уменьшения расхода памяти:

• Для хранения графа и множеств используются битовые массивы
• Вся цепочка подразбиений от корня дерева к текущему узлу кодируется всего двумя массивами (lab и ptn). При возвращении к родительскому узлу разбиение восстанавливается функцией recover

Подробности можно узнать в nauty User’s Guide и в коде самой nauty.


Применение

Прелесть изложенного алгоритма в том, что в нём не много завязок на то, что исходное множество вершин — это вершины простого графа. Фактически, можно находить автоморфизмы и каноническую нумерацию для чего угодно, например:

• Для ориентированных графов. Никаких изменений в алгоритме не требуется. МакКей в статье упоминает ориентированные графы, и в nauty есть соответствующая опция digraph
• Для графов с раскрашенными вершинами. Всё, что требуется — начинать не с элементарного разбиения, а сгруппировать вершины по цветам. Эта опция также описана у МакКея и имеется в nauty
• Для графов с раскрашенными рёбрами. Всё, что требуется — изменить функцию cmp так, чтобы она учитывала цвета рёбер. Также, можно более эффективно написать процедуру Refine (уточнение разбиения)
• Для молекул со стереохимическими свойствами. Опять же, модификации требуются только в функции cmp

Модификации алгоритма

Не было упомянуто об огромном количестве дополнительных оптимизаций, имеющихся в коде nauty. МакКей признаёт, что у него самого не все они задокументированы. Вот только некоторые:

• «Даровые» автоморфизмы (cheap automorphisms, implicit automorphisms), для обнаружения которых не требуется доходить до дискретного разбиения. Работают в nauty только для неориентированных графов (цвета вершин допускаются)
• Проверка автоморфизма cmp(π,ξ) может делаться более эффективно без вызова cmp, с учётом того, что ξ никогда не меняется и нас интересует только то, есть ли изоморфизм
• Псевдо-канонические коды разбиений, благодаря которым можно выявить отсутствие автоморфизма, не проверяя его явно. Они же используются для предварительного определения знака функции cmp
• При выборе ячейки V_i для подразбиения можно применять различные эвристики, с целью сокращения перебора

Рекурсивные функции на Common Lisp

1. Перебор всех подписков данного списка

(defun for-all-subsets (list func)
(unless list
  (funcall func nil))
(when list
  (for-all-subsets
   (cdr list)
   (lambda (s) (funcall func (cons (car list) s))))
  (for-all-subsets (cdr list) func)))

func — функция одного параметра — вызывается для каждого подсписка.

2. Перебор всех подсписков заданного размера

(defun for-all-subsets-size (list size func)
(if (zerop size)
    (funcall func nil)
    (when list
      (for-all-subsets-size
       (cdr list) (1- size)
       (lambda (s) (funcall func (cons (car list) s))))
      (when (< size (length list-set))
         (for-all-subsets (cdr list-set) size func)))))

func — функция одного параметра — вызывается для каждого подсписка.

3. Перебор всех подсписков, полученных из данного удалением одного элемента

(defun for-all-cuts (list func)
(declare (list list) (function func))
(when list
  (funcall func (car list) (cdr list)))
(when (cdr list)
  (for-all-cuts (cdr list)
   (lambda (elem cut)
     (funcall func elem (cons (car list) cut))))))

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

4. Перебор максимальных по включению подсписков, удовлетворяющих какому-либо условию

(defun for-all-max-subsets (set func)
(if set
    (let ((saved))
      (for-all-max-subsets
       (cdr set)
       (lambda (s)
         (when (funcall func (cons (car set) s))
           (push s saved))))
      (for-all-max-subsets
       (cdr set)
       (lambda (s)
         (or (position s saved :test #'equalp)
             (funcall func s)))))
    (progn
      (funcall func nil)
      nil)))

func — функция одного параметра — получает список и должна возвращать nil, если список не годится. Если для какого-то списка func вернула ненулевое значение, все подсписки этого списка не будут проверяться.

5. Перебор всех циклических перестановок данного списка

(defun for-all-cyclic-permutations (list func)
  (funcall func list)
  (do ((iter (cdr list) (cdr iter))
             (prev list))
      ((not iter))
    (when prev
       (setf (cdr prev) nil)
       (funcall func (append iter list))
       (setf (cdr prev) iter))
       (setf prev iter)))

func — функция одного параметра — вызывается для каждой циклической перестановки.

Примеры:

> (for-all-subsets '(a b c) #'print)

(A B C)
(A B)
(A C)
(A)
(B C)
(B)
(C)
NIL

> (for-all-subsets-size '(a b c d) 2 #'print)

(A B)
(A C)
(A D)
(B C)
(B D)
(C D)

> (for-all-cuts '(a b c d) (lambda (a b) (format t "~a ~a~%" a b)))

A (B C D)
B (A C D)
C (A B D)
D (A B C)

> (for-all-max-subsets
   '(2 10 2 5 4)
   (lambda (x) (if (> (apply #'+ x) 10) nil (progn (print x) t))))

(2 2 5)
(2 2 4)
(10)
(5 4)

> (for-all-cyclic-permutations '(a b c d e) #'print)

(A B C D E) 
(B C D E A) 
(C D E A B) 
(D E A B C) 
(E A B C D) 

[Oracle] Построчное хранение данных

Минимальная единица информации в хранилище Oracle называется блок.

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

При чтении какой-либо записи Oracle берёт с диска (или из кеша) весь блок, в котором находится эта запись.

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

Пример на Oracle 11:

create table a (n number, s varchar2(4000));
create table b (n number);
begin
  for i in 1..1000000 loop
    insert into a values(i, lpad('a', 4000, 'a'));
    insert into b values(i);
  end loop;
end;

set timing on;

select sum(n) from a;

    SUM(N)
----------
5.0000E+11

Elapsed: 00:02:04.51

select sum(n) from b;

    SUM(N)
----------
5.0000E+11

Elapsed: 00:00:00.11

Хранения «по колонкам» в Oracle не предусмотрено. Для лучшего быстродействия следует избегать объединения в одной таблице полей, которые могут быть востребованы по отдельности.

practical-common-lisp.{fb2, pdf}

Перегнал Practical Common Lisp в формат fb2, не без помощи OpenOffice с ooofbtools, а также ручной правки.
Сноски, таблицы, иллюстрации, цитаты и код оформлены правильными fb2-тегами. Книга будет нормально выглядеть на LBook и прочих карманных девайсах.
Скачать можно здесь.
Алсо, pdf-файл с книгой раньше находился на сайте издательства Apress, но почему-то пропал. Здесь лежит копия.