Содержание к диссертации
Введение
1. Вычисление элементарных функций с гарантированной точностью 19
1.1. Соглашения 23
1.2. Оценка синуса на интервале [0, |] 27
1.3. Оценка синуса на интервале [0,оо] 33
1.4. Оценка косинуса 37
1.5. Оценка тангенса 38
1.6. Обратные тригонометрические функции 38
1.6.1. Оценка арксинуса
1.6.2. Оценка арккосинуса 39
1.6.3. Оценка арктангенса 39
1.7. Оценка экспоненты 39
1.8. Оценка логарифма 41
1.9. Результаты 42
2. Алгоритм коррекции решения 45
2.1. Базовые определения 46
2.2. Алгоритм М2В 47
2.3. Применение алгоритма М2В 54
2.4. Алгоритм МЗВ 56
2.5. Алгоритм коррекции решения 59
2.6. Настроечные параметры алгоритма 62
2.7. Результаты экспериментов 63
2.8. Управление алгоритмами сужения 68
3. Ньтоновские корректоры 69
3.1. Граф задачи 70
3.2. Функции-характеристики 73
3.3. Граф истории 76
3.3.1. Вершины графа истории 76
3.3.2. Создание заготовки графа истории 78
3.3.3. Эволюция графа истории 79
3.3.4. Свойства графа истории 82
3.4. Граф-индикатор несовместности 84
3.4.1. Вершины-числа 85
3.4.2. Вершины-функции 88
3.4.3. Корректность графа-индикатора 89
3.4.4. Предельное значение графа-индикатора 89
3.4.5. Теорема о предельном значении графа-индикатора 91
3.5. Алгоритм ньютоновской коррекции 92
3.5.1. Теорема о корректности алгоритма 94
3.6. Сравнение алгоритмов NSC и М2В на примере конкретной задачи 99
3.6.1. Решение задачи М алгоритмом М2В
3.6.2. Решение задачи М алгоритмом NSC 101
3.7. Результаты экспериментов 111
Заключение 123
Литература 126
- Обратные тригонометрические функции
- Алгоритм коррекции решения
- Граф-индикатор несовместности
- Сравнение алгоритмов NSC и М2В на примере конкретной задачи
Введение к работе
Большинство задач, возникающих при моделировании, составлении расписаний, оптимизации распределения ресурсов, могут рассматриваться как задачи удовлетворения ограничений (CSP - Constraint Satisfaction Problem).
Неформально задачу удовлетворения ограничений можно определить как задачу, составленную из конечного количества переменных, каждой из которых поставлена в соответствие область ее возможных значений, и конечного множества ограничений над всеми или некоторыми из этих переменных. В таких задачах требуется сопоставить каждой переменной некоторое значение из области ее возможных значений так, чтобы выполнялись одновременно все ограничения.
Задача удовлетворения ограничений впервые была формализована Хаффманом (Huffman) в 1971 г. [43].
Одним из классических примеров задачи, которая может быть представлена как задача удовлетворения ограничений, является задача расстановки N ферзей на шахматном поле размером N X N так, чтобы ни один из ферзей не "бил" другого. В такой постановке задачи мы можем построить CSP, в которой каждому ферзю соответствует переменная, область значений каждой переменной - множество полей шахматной доски, а множество ограничений состоит из n ^ ' ограничений, запрещающих каждой паре ферзей "бить" друг друга, т. е. находиться на одной строке, в одном столбце или на одной диагонали.
Рассмотрим другой пример. Пусть множество переменных задачи удовлетворения ограничений состоит из единственной переменной х, множество ее возможных значений - множество вещественных чисел, а множество ограничений состоит из ограничения х — 2х + 1 = 0. Таким образом, обычное математическое уравнение может быть представлено как задача удовлетворения ограничений. Различные системы уравнений, в том числе и с не полностью определенными коэффициентами, также могут быть рассмотрены как задачи удовлетворения ограничений.
Задачи, в которых множества возможных значений всех переменных являются подмножествами множества вещественных чисел, называют численными (numeric CSP, NCSP). Задачи, в которых множества возможных значений переменных конечны, называют конечными (finite CSP). Задачи, в которых множества возможных значений переменных являются вещественными интервалами, называют непрерывными. Задачи, все ограничения которых унарные либо бинарные, называют бинарными.
Большинство методов, посвященных решению задач удовлетворения ограничений, предназначены для решения задач с конечными областями. Среди этих методов с некоторыми допущениями можно выделить три базовых класса [75]:
Преобразование задачи. Методы, заключающиеся в замене исходной задачи эквивалентной, но предположительно более легкой для решения.
Поиск. Различные методы поиска решения.
Синтез. Методы синтеза решения.
Большая часть остальных методов является некоторой суперпозицией базовых.
Преобразование задачи - это уменьшение областей возможных значений переменных и/или изменение множества ограничений.
Приведенную выше задачу расстановки ферзей можно изменить следующим образом. Учитывая, что в одном столбце шахматной доски может находиться ровно один ферзь, сопоставим каждому ферзю столбец, в котором он будет находиться. В этом случае областью возможных значений каждой переменной будет множество строк шахматной доски. С каждой парой переменных будет связано более слабое ограничение. Оно будет запрещать этим переменным иметь такие значения, при которых ферзи, которым соответствуют эти переменные, находятся на одной строке или на одной диагонали. Решить CSP в такой постановке оказывается существенно проще.
Многие методы решения CSP применимы только к задачам специального вида, например к задачам, все ограничения которых унарные либо бинарные. Поэтому разработано множество алгоритмов приведения произвольных CSP к тому или иному специальному виду. Выполнение различных линейных подстановок, раскрытие скобок, применение формул тригонометрических преобразований, построение базисов Гребнера [18, 41, 57, 64] для систем полиномиальных уравнений также относится к этому классу методов.
Уменьшение областей значений состоит в удалении так называемых избыточных значений, т. е. таких значений, удаление которых не изменяет множества решений CSP. Методы, удаляющие избыточные значения, называются методами сужения. Среди последних выделяют методы достижения совместности, заключающиеся в удалении тех значений, которые нарушают некоторое свойство, называемое свойством совместности.
Среди различных видов совместностей для конечных задач наиболее известны совместности в узле, совместности по дугам, по путям [75]. Названия совместностей основаны на представлении CSP в виде графа в том случае, если задача бинарна, или гиперграфа в общем случае. В таком представлении каждая переменная задачи соответствует некоторой вершине, а набор ограничений, связывающих переменные хку, представляется дугой, соединяющей вершины, соответствующие переменным ж и у.
Задача называется совместной в узлах, если для каждой переменной все значения области ее возможных значений удовлетворяют унарным ограничениям над этой переменной.
Стандартный алгоритм достижения совместности в узлах называется NC (node consistency) [75] и состоит в последовательном переборе всех возможных значений каждой переменной и проверке всех унарных ограничений над этой переменной.
Бинарная CSP называется совместной по дугам, если для любой пары переменных (х,у), для каждого значения а переменной ж, удовлетворяющего всем унарным ограничениям над х, существует значение Ъ переменной у такое, что все бинарные ограничения над х и у выполняются.
Для достижения совместности по дугам было создано несколько алгоритмов [28], наиболее распространенный из них - АС-3 [53]. Алгоритм применяется к задаче, совместной в узлах. Сначала он строит очередь Q ограничений, в которую изначально ставит все ограничения решаемой задачи; далее алгоритм удаляет по очереди все ограничения из Q, выполняя следующие действия: если с(х, у) - удаляемое ограничение, то удаляет значения а из области значений переменной х, для которых не существует такого значения Ъ в области значений у, что с(а, Ь) выполняется; удаляет значения d из области значений у, для которых не существует такого значения е в области значений ж, что с(е, d) выполняется; если область возможных значений переменной х изменилась, ставит в очередь Q все ограничения, в которых фигурирует переменная х, кроме тех, что уже стоят в очереди Q, и кроме ограничения с(х,у); если область возможных значений переменной у изменилась, ставит в очередь Q все ограничения, в которых фигурирует переменная у, кроме тех, что уже стоят в очереди Q, и кроме ограничения с(х,у).
Алгоритм завершает свою работу если область значений некоторой переменной или очередь Q становятся пустыми. В первом случае исходная задача не имеет решений, во втором - полученная совместна по дугам.
Бинарная CSP называется совместной по путям, если для любой такой последовательности переменных (xi,..., хп), что значения ai,an переменных xi, хп удовлетворяют всем ограничениям над ж і, хп, существуют значения аг,...,ап-і переменных Х2,...,хп-і такие, что для всех г = 1,..., п — 1 значения аг, а,і+і переменных жг, ж,-+і удовлетворяют всем ограничениям над переменными Хі,Хі+і.
Для достижения совместности по путям также создано несколько алгоритмов [75], однако все они практически не используются, поскольку уступают в эффективности различным суперпозициям алгоритмов поиска и достижения совместности по дугам.
Для непрерывных задач также существуют понятия локальной совместности: совместности по границам (2В-, ЗВ-, ...), box-совместность и другие [25, 49]. Методы достижения перечисленных и некоторых других видов совместностей называют также методами распространения ограничений [21, 24, 27, 46, 67]. Такие методы могут оказаться единственным подходом к решению реальных задач, когда применение классических вычислительных методов или методов компьютерной алгебры затруднено наличием неточных данных или частично определенных параметров.
Второй класс методов решения CSP содержит различные методы поиска решения. Большинство существующих методов решения конечных CSP принадлежит этому классу.
Основным алгоритмом поиска для конечных CSP является бэктрэкинг (backtracking) [33, 65], известный также как хронологический бэктрэкинг (chronological backtracking). Этот алгоритм состоит в последовательном фиксировании переменных, т. е. присваива- ний им некоторых значений. Изначально все переменные незафиксированны. На каждом шаге выбирается некоторая незафиксированная переменная, после чего ищется "совместное" значение из области ее возможных значений, т. е. такое значение переменной, которое вместе со значениями фиксированных переменных не нарушает ни одного ограничения над этими переменными. Если область значений выбранной переменной не содержит совместных значений, то значение, присвоенное при фиксировании предыдущей переменной, признается несовместным, и происходит возврат на предыдущий шаг. В противном случае найденное совместное значение фиксируется, после чего происходит переход к следующему шагу, на котором выбирается и фиксируется следующая переменная. Если все переменные оказались фиксированными - решение найдено. Если все значения первой выбранной переменной оказались несовместными - задача не имеет решений.
В различных реализациях бэктрэкинга могут быть различными стратегии выбора: переменной для фиксирования; значения для присваивания; ограничения для проверки.
Удачно сделанный выбор может уменьшить пространство поиска.
Существует большое количество модификаций бэктрэкинга. Рассмотрим некоторые из них.
Алгоритм backchecking [40] разработан для приложений, в которых проверка совместности значений занимает много времени. Если при фиксировании некоторой переменной у, значение Ъ оказывается несовместным с некоторым значением а фиксированной переменной х, backchecking запоминает это и не проверяет больше значение Ь, пока значением переменной х остается а.
Алгоритм backmarking [40] является улучшением алгоритма backchecking. Как и backchecking, он уменьшает количество проверок совместности, запоминая, какие значения несовместны с ранее фиксированными переменными. Более того, он запоминает, какие значения совместны с ранее фиксированными переменными, и тоже их не проверяет.
Предположим что на текущем шаге ищется совместное значение для переменной Xj. Предположим также, что некоторое время назад все возможные значения переменной Xj были проверены и найдены несовместными. Предположим, что с тех пор backmarking произвел откат до переменной ж*, где і < j, присвоил ей новое значение, после чего нашел совместные значения для переменных Жі+і, ..., х$-\. Тогда при поиске совместного значения для переменной Xj backmarking не будет рассматривать те значения, которые "конфликтуют" с переменными х\,..., Жі-і, про них уже известно, что они несовместны. Кроме того, backmarking не будет проверять совместность между хі,... ,х^г и теми из оставшихся значений, про которые заранее известно, что они совместны.
Следующий алгоритм - backjumping [69] отличается от бэктрэкинга тем, что при необходимости возврата к предыдущему шагу он проводит анализ и возвращается сразу к той переменной, фиксирование которой привело к отсутствию совместных значений на текущем шаге.
Алгоритм forward checking [40] при фиксировании каждой переменной удаляет из областей возможных значений нефиксированных переменных значения, несовместные со значениями уже фиксированных переменных. При этом, если область значений некоторой переменной становится пустой, фиксируемое значение признается несовместным.
Существуют также алгоритмы, которые после фиксирования каждой переменной вызывают некоторый алгоритм сужения. Такие алгоритмы наряду с forward checking называются lookahead-алгоритмами [75].
Другой подкласс алгоритмов поиска - стохастические алгоритмы [50, 51, 58, 72, 74]. Они обьічно применяются в тех случаях, когда необходимо быстро найти одно произвольное решение задачи. Алгоритмы стохастического поиска - это такие алгоритмы поиска, которые включают эвристики и элементы случайности. Большинство стохастических алгоритмов не гарантируют свое завершение за конечное время, и это их самый серьезный недостаток. Их преимущество состоит в том, что они могут оказаться единственно способными найти решение некоторых задач за разумное время. Так, например, алгоритм QS4 [19] решает задачу расстановки N ферзей меньше чем за минуту, если N < 1000000, тогда как бэктрэкинг уже при N, равном 100, не находит ни одного решения за разумное время.
Существуют также эвристические алгоритмы поиска, которые гарантируют свое завершение за конечное время. Примером такого алгоритма является одна из разновидностей бэктрекинга - алгоритм iterative broadening [31]. Его идея состоит в добавлении элементов случайности при распределении вычислительных усилий по дереву поиска. Количество визитов каждого узла дерева ограничивается некоторым пределом, и при достижении этого предела все не рассмотренные ветви этого узла игнорируются. Если при заданном пределе количества посещений алгоритм не находит решение, предел увеличивается.
Некоторые алгоритмы поиска могут быть использованы для бесконечных и смешанных CSP. Так, например, если области возможных значений переменных представлены ограниченными интервалами, то в качестве "значения" переменной может быть зафиксирован некоторый достаточно узкий интервал, содержащийся в области возможных значений переменой.
Обычно алгоритмы поиска применяются вместе с алгоритмами преобразования задачи. Поскольку последние сужают пространство поиска, целесообразно использовать их перед применением алгоритмов поиска или, например, после фиксирования каждой переменной.
Третий класс алгоритмов решения CSP - алгоритмы синтеза решения [75] - применяют тогда, когда необходимо найти все решения данной задачи. Пусть {х^,..., хп} - множество переменных некоторой CSP, С - множество ее ограничений. Пусть множество С{ содержит все ограничения задачи, связывающие только переменные из подмножеств множества {х\,..., Хі}. Синтез решения - это построение сначала множества значений переменной х\, которые удовлетворяют ограничениям из множестра Сі, затем множества пар значений переменных Хі и #2, которые не нарушают ограничения из множества Сг, и т. д. В итоге получим множество кортежей длины п значений переменных xi,...,xn, которые удовлетворяют ограничениям из множества Сп, т. е. множество решений CSP. Первый алгоритм синтеза решения был представлен Фредером (Freuder) в 1978 г. [29].
Как отмечено выше, представление в виде CSP изначально использовалось для задач с конечными областями значений переменных, поэтому большая часть существующих методов решения CSP предназначена для решения конечных CSP. Параллельно с разработкой методов решения конечных CSP разрабатывались методы решения и непрерывных задач. Были обобщены понятия совместностей, для чего использовалось интервальное представление данных [1, 4, 5, 17, 39, 60].
В 1980-х гг. А. Нариньяни предложил метод недоопределенных вычислений (МНВ) [13, 14], который является методом сужения, т. е. относится к первому классу методов решения CSP. МНВ позволяет работать как с переменными, области значений которых конечны, так и с бесконечными и непрерывными областями значений. В этом методе использовано интервальное представление вещественных чисел и областей возможных значений вещественных переменных. Если множества значений всех переменных задачи конечны, алгоритм является алгоритмом достижения совместности по дугам и полностью соответствует алгоритму АС-3. На задачах, в которых области значений всех переменных вещественные интервалы, МНВ является алгоритмом достижения 2В-совместности и соответствует алгоритму М2В [49]. Для решения смешанных задач МНВ использует механизмы обоих методов.
При решении непрерывных задач применяют и специальные методы, неприменимые для конечных CSP [37, 38, 39, 42]. Обычно такие методы неприменимы и для смешанных CSP, более того, они обычно могут применяться только к квадратным системам уравнений, т. е. системам, количество уравнений которых совпадает с количеством переменных. Один из наиболее известных таких методов - интервальный метод Ньютона [39]. Для решения уравнения /(ж) = 0, где х Є [а, 6], а функция / дифференцируема в замкнутом промежутке [о, Ь], строится последовательность вложенных интервалов [a, b] = Iq С 1\ С ... С 1п, где 1І+1 — Ji І І V-t-i fl ґт \)->
Х{ t xj.
По теореме Лагранжа о среднем значении [3] все корни уравнения f(x) = 0, содержащиеся в /0, содержатся в 1п. Таким образом, интервальный метод Ньютона является методом сужения. Существуют также специальные методы для решения систем полиномиальных уравнений [23, 26, 30].
В 1988 г. в РосНИИ искусственного интеллекта создан решатель систем нелинейных уравнений с не полностью определенными данными UniCalc [2, 20]. Базовым алгоритмом решателя является метод недоопределенных вычислений. Решатель работает как с вещественными, так и с целочисленными переменными и допускает представление ограничений как в виде конечных суперпозиций элементарных функций и операций, так и в виде логических выражений. В РосНИИ ИИ были построены также пробные версии решателя с табличным представлением ограничений, с мультиинтервальным представлением областей значений переменных [16, 68] и др. [6, 45]. На достаточно большом классе задач решатель UniCalc показал существенное преимущество перед обычными методами вычислительной математики. Позже библиотека решателя UniCalc была расширена некоторыми алгоритмами поиска, разработанными для задач с конечными и непрерывными областями.
В интервальных методах для вычисления границ интервалов применяются интервальные математические библиотеки, в которых используются различные подходы к вычислению интервальных функций. Пусть значение некоторой переменной х содержится в некотором интервале [а, Ь]. Предположим, что необходимо оценить значение ех. Очевидно, что оно содержится в интервале [еа,е6], но поскольку числа еа и еь могут оказаться не машинно-представимыми, для оценки значения ех часто используют следующие интервалы: интервал, границы которого - приближенные значения чисел еа и еь. Для вычисления приближенных значений можно использовать, например, функции стандартных математических библиотек, поставляемых с компиляторами с языка С [32]; интервал, полученный так же, как и первый, расширенный на некоторое заранее посчитанное е. гарантированный интервал, т. е. интервал, нижняя граница которого - максимальное машинно-представимое число, не превосходящее еа, а верхняя граница - минимальное машинно-представимое число, не меньшее еь;
Заметим, что в первом случае истинное значение ех может не содержаться в том интервале, которым оно представляется. Однако в случае, когда вычисления носят итерационный характер, положительные погрешности одних итераций с большой вероятностью компенсируются отрицательными погрешностями других.
Второй способ применяют при оценке значений выражений, являющихся суперпозицией нескольких функций. В этом случае часто вычисляют каждую функцию первым способом, а затем расширяют полученный интервал. При этом положительные и отрицательные погрешности также компенсируют друг друга, и вероятность того, что истинное значение выражения лежит в полученном интервале, очень высока, хотя и не равна единице. В решателе UniCalc используется именно эта схема оценивания. Однако при решении больших задач даже небольшая вероятность ошибки в каждом конкретном выражении может привести к "потере" корней задачи.
Наконец, третий способ называется вычислением с гарантированной точностью. Несмотря на то что при использовании его в итеративном процессе будет происходить постепенное расширение интервала, этот способ позволяет гарантировать, что истинное значение некоторого выражения лежит в интервале, которым это значение представлено.
Анализ доступных интервальных математических библиотек показал, что они производят вычисления либо первым [32], либо вторым способом [55], причем алгоритмы из [32] используются практически во всех библиотеках, поставляемых с компиляторами. При исследовании алгоритмов из [55] было обнаружено, что они также не гарантируют наличия точного значения в получаемом интервале. В связи с этим возникла необходимость построения собственной математической библиотеки для гарантированного вычисления элементарных математических функций.
В 1990-х гг. стали появляться программные продукты, например такие, как Numerica [77], которые на некоторых задачах превосходили UniCalc по быстродействию. В рамках работ по улучшению быстродействия решателя UniCalc и расширения его библиотеки методами решения различных типов задач автор начал решать задачу разработки дополнительных методов сужения для непрерывных и смешанных CSP.
Целью диссертационной работы является создание алгоритмов сужения пространства поиска для задач с непрерывными и смешанными областями и доказательство сходимости и корректности разрабатываемых алгоритмов. При создании указанных алгоритмов автор столкнулся с необходимостью построения библиотеки для эффективного гарантированного вычисления оценок элементарных функций. Эта библиотека является отдельным результатом данной работы, который может быть использован при разработке других интервальных методов.
Диссертационная работа состоит из введения, трех глав и заключения. Общий объем работы - 134 страниц. Работа включает 3 таблицы и 6 рисунков.
Первая глава посвящена проблеме эффективного вычисления элементарных функций. После краткого обзора существующих подходов предлагаются алгоритмы нахождения гарантированных оценок значений элементарных функций, время вычисления которых уменьшено за счет увеличения размера используемой памяти. В главе приведены подробные алгоритмы вычисления верхней и нижней оценок значения синуса. Кратко описаны структуры алгоритмов для вычисления других элементарных функций.
Во второй главе представлен новый эвристический алгоритм достижения ЗВ-совмест-ности: алгоритм коррекции решения (SC-алгоритм). Он является альтернативой существующему алгоритму, который обычно приводится в литературе как стандартный алгоритм достижения ЗВ-совместности. Основой обоих алгоритмов является набор так называемых корректоров. При этом SC-алгоритм осуществляет более эффективное управление корректорами. В главе приведены доказательства корректности и сходимости алгоритма, а также результаты экспериментов.
В третьей главе представлен алгоритм ньютоновской коррекции. Построенные на его основе корректоры предназначены для за^лы стандартных корректоров, используемых в SC-алгоритме. Так же, как и в предыдущей главе, здесь приведены доказательства корректности и сходимости предлагаемого алгоритма, результаты экспериментов.
В заключении сформулированы основные результаты диссертации.
По теме диссертации опубликовано 8 работ. Основное содержание составляют резуль- таты работ [8, 11, 12, 52].
Результаты диссертации докладывались автором на XV Международной конференции по интервальной математике (г. Красноярск, 17-19 августа 1999 г.), конференции молодых ученых, посвященной 10-летию ИВТ СО РАН (г. Новосибирск, 25-26 декабря 2000 г.), конференции молодых ученых по математике, математическому моделированию и информатике (г. Новосибирск, 4-6 декабря 2001 г.), на семинарах РосНИИ ИИ, ИСИ СО РАН.
Обратные тригонометрические функции
Не все вещественные числа могут быть точно представлены в компьютере, чаще всего различные форматы представления (типы) вещественных чисел сопоставляют множество битовых последовательностей определенной длины некоторому конечному подмножеству множества вещественных чисел. Числа, входящие в такое подмножество, обычно называют числами данного типа, про остальные говорят, что они не являются числами данного типа.
Для чисел второго вида существуют различные методы представления. Среди таких методов широко распространены представление вещественного числа ближайшим числом нужного типа и представление содержащим данное число интервалом, границы которого - числа нужного типа.
Представление вещественных чисел в виде интервалов используется в методах интер вальной математики [1, 4, 5, 17, 39, 60], в методах интервального распространения огра ничений [21, 24, 27, 46, 67], в методах решения дифференциальных и обычных уравнений.
В таких методах важно корректно оперировать с интервалами так, чтобы интерваль- ный результат арифметических операций и элементарных математических функций гарантированно содержал множество значений вычисляемой операции или функции от интервалов, которыми представлены аргументы.
Поскольку результат операции или функции от машинно-представимых чисел может не являться машинно-представимым числом, необходимо уметь вычислять нижнюю и верхнюю мащинно-представимые оценки значения этого результата. Для вычисления нижних и верхних оценок результатов операций применяются так называемые флаги округления [44]. Флаг округления - это флаг состояния процессора, который влияет на результат ряда арифметических действий, выполняемых процессором, таких как, например, вычисление операций, приведение типов, округление и др. В процессорах, удовлетворяющих стандарту іеее 754 [66], существует четыре состояния флага округления: "вниз" (или "к — оо"), "вверх" (или "к +оо"), "к нулю" и "к ближайшему". Если полученное значение некоторой операции из списка операций, специфицированных в стандарте іеее 754, не является числом того типа, которым должен быть результат, то в зависимости от состояния флага округления полученное значение округляется соответственно вниз, вверх, в сторону нуля или к ближайшему числу нужного типа.
Например, вычисление левой и правой оценок значения суммы чисел а и Ь может выглядеть как с+ = а + Ъ. Здесь под set_round(down) имеется в виду установка флага округления в значение "вниз", а под set_round(up) - в значение "вверх". Установка флага округления - достаточно дорогая операция в том смысле, что она занимает достаточно много времени. Поэтому более эффективной является следующая реализация оценки сложения: tmp. Стоит отметить, что некоторые оптимизаторы могут изменить написанный код так, что скомпилированная программа будет выполнять не в точности те действия, которые от нее ожидает разработчик, в результате чего полученные значения будут неверными. Исключить подобные действия оптимизатора достаточно легко при использовании определенного набора опций, специального для каждого конкретного компилятора. Эта тема не будет обсуждаться в работе.
При вычислении оценок значений математических выражений недостаточно простой установки флага округления перед вычислением. Так, например, для вычисления оценок произведения мультипликативного ряда необходимо учитывать знаки всех множителей. Для вычисления оценок значения элементарных математических функций также недостаточно установки флага округления и вызова соответствующей "стандартной" функции. Библиотеки, поставляемые с компиляторами с языка С, как, например, библиотека Glibc [32], предназначены для приближенного вычисления элементарных функций. В этих библиотеках для вычисления элементарных функций используются мультипликативные ряды, вычисление которых с направленными округлениями неэффективно как по времени, так и по качеству.
В настоящее время создано достаточно большое количество различных алгоритмов для вычисления элементарных функций. В [22] приведены ссылки на 227 работ по этой теме. Основной недостаток опубликованных в Интернете библиотек, находящихся в свободном доступе, - это отсутствие гарантий со стороны авторов корректности их работы. Например, библиотека fdlibm [55] не строится по причине отсутствия каких-то функций, в C-XSC [78] не работают поставляемые в комплекте примеры.
Один из распространенных методов вычисления элементарных функций с направленными округлениями состоит в приближенном вычислении искомого значения с большей точностью и округлении полученного результата в необходимую сторону, т. е. для вычисления m двоичных цифр значения /(ж) вычисляется приближенное значение f(x) с относительной точностью 2 п, где п т. После чего берется т знаков полученного значения. В [63] проводится исследование зависимости значения п от количества искомых гарантированных знаков т.
Такой подход эффективен для разовых вычислений, когда время вычисления функции не критично и важно найти действительно ближайшее представимое число. В методах интервального распространения ограничений необходимым является лишь требование, чтобы нижняя оценка оцениваемого значения была меньше, а верхняя - больше самого значения. Требование минимизировать разность между верхней и нижней оценками обычно является желательным, но не необходимым.
Иными словами, быстрый алгоритм, получающий гарантированную, но иногда не оптимальную оценку значения функции, часто оказывается полезнее, чем более медленный, но всегда получающий оптимальную оценку.
В связи с изменением ситуации на компьютерном рынке, в частности, в связи с существенным удешевлением оперативной памяти, алгоритмы, использующие таблицы заранее посчитанных констант, находят все большее распространение. Одной из первых работ, в которых используются огромные таблицы констант для быстрого и точного вычисления элементарных функций, является [35], в которой предложен алгоритм АТА-М. Алгоритм использует таблицы с суммарным размером 2 Мбайта; его авторы утверждают, что время вычисления основных элементарных функций составляет Зт - 4т, где г - время выполнения умножения двух чисел с двойной точностью [44]. Однако алгоритм АТА-М предназначен для приближенного вычисления элементарных функций, а не для гарантированного вычисления оценок.
Здесь представлены алгоритмы для вычисления гарантированных оценок элементарных математических функций. Эти алгоритмы будут использованы для вычисления границ интервалов в алгоритме коррекции решения и в ньютоновских корректорах, представленных в двух следующих главах. С учетом этого в главе рассматриваются вопросы эффективности и целесообразности. Так, разность между верхней и нижней оценками значения функций должна быть как можно меньшей. Алгоритмы должны работать корректно при машинном флаге округления, установленном в значение "вниз".
Алгоритм коррекции решения
Каждый цикл внутри функции Фзв, кроме последнего, изменяет область значений хотя бы одной переменной. Так как количество FP-чисел в областях значений всех переменных конечно, время исполнения каждого корректора также конечно, то количество циклов внутри функции Фзв ограничено. Следовательно, алгоритм МЗВ сходится. Поскольку работа алгоритма МЗВ завершается, только если обнаружена несовместность или все корректоры становятся неудачными, из приведенных выше свойств ЗВ-корректоров следует, что в случае нормального завершения работы алгоритма полученная в результате применения МЗВ задача будет ЗВ-совместной. Таким образом, мы показали сходимость и корректность алгоритма МЗВ.
Алгоритм МЗВ позволяет существенно улучшить оценку решения, получаемую алгоритмом М2В. Выполнив модификацию алгоритма, можно увеличить его скорость.
Сформулируем основные цели, которые мы стремимся достичь. Во-первых, мы хотим чаще вызывать такие корректоры, от которых ожидаем наибольшего уточнения за единицу времени; во-вторых - прерывать работу корректоров, которые долго работают и не сужают свою область. Выбрав правильную стратегию определения моментов для прерывания работы корректоров и определения очередности их вызовов, мы тем самым ускорим работу алгоритма.
Для определения моментов прерывания понадобится счетчик времени. Использование действительно часов неразумно, поскольку работа алгоритма становится зависимой от загруженности процессора, его мощности и т. д. Поэтому в качестве единицы времени используется один вызов функции интерпретации внутри функции ФгвО При последующих запусках прерванных корректоров будем продолжать вычисления с того места, на котором они были прерваны, учитывая, что с момента предыдущего запуска этого корректора значения переменных могли измениться. Это означает, что обновляются значения в стеке корректора, и если прерывание произошло во время выполнения метода М2В, в очередь последнего добавляются ограничения, в которые входят изменившиеся переменные.
Если предыдущий вызов корректора не был прерван, т. е. завершил свою работу по причине удовлетворения условия 2В-совместности границы, проверим еще раз 2В-совместность, учитывая изменившиеся переменные, и в случае 2В-несовместности продолжим выполнение корректора. Для того чтобы сформулировать условия прерывания корректора, нам понадобится определение его текущей расчетной скорости. Дадим его в тех обозначениях, которые были использованы при описании функции l_corrector(і, Х, D, О): Определение 2.11. Пусть ХІ - главная переменная корректора. Пусть в настоящий момент корректором вызвана функция Фгв которая уже совершила п вызовов функций интерпретации. Обозначим Di область возможных значений переменной ХІ в той задаче, для которой был вызван корректор (см. описание функции Lcorrector). Обозначим Di область возможных значений переменной х в той задаче, для которой была вызвана функция Фгв Тогда если корректор левый, его текущей расчетной скоростью будем называть значение отношения — —«—, /,: если правый, - отношения . —, .. Теперь можно сформулировать условия прерывания работы корректоров. Условие 1. Работа корректора прерывается, если время с момента его вызова превышает некоторое значение tcau. Условие 2. Работа корректора прерывается, если время работы вызванной им функции Фгв() превышает некоторое значение 2в Условие 3. Работа корректора прерывается, если его текущая расчетная скорость станет меньше некоторого значения scorr. Значения tcau, tiB и scorr зависят от размера решаемой задачи и вычисляются некоторыми эмпирическими формулами. Далее, для определения очередности вызовов корректоров введем штрафы. Все корректоры вызываются по очереди, но если на какой-то корректор наложен штраф п, он пропускается п раз. Размер штрафа вычисляется эмпирической формулой в зависимости от результата работы корректора. Максимальный размер штрафа должен быть ограничен некоторой константой Р.
Для достижения ЗВ-совместности алгоритм должен завершать свою работу в том случае, если обнаружится несовместность или все корректоры станут неудачными (прерванный корректор не считается неудачным). Однако, как будет отмечено в следующем пункте, на практике можно прекращать работу SC-алгоритма несколько раньше.
Доказательство. Предположим, что SC-алгоритм работает бесконечно долго. Поскольку границы областей значений переменных - FP-числа, при уменьшении области значений любой переменной, она теряет как минимум одно FP-число. Так как область значений каждой переменной содержит конечное количество FP-чисел, и алгоритм завершает работу в случае, если хотя бы одна из областей стала пустой, то, начиная с некоторой итерации, области значений всех переменных не меняются. Пусть k номер этой итерации (итерациями будет считать циклы, см функцию Фзв)- Пусть Df,..., D - области значений переменных после к итераций.
Поскольку размер штрафа ограничен сверху, каждые Р итераций будет происходить как минимум один запуск каждого корректора. Так как, начиная с итерации с номером к, области значений переменных не меняются, то при запусках корректоров в их очереди не будут ставиться новые ограничения. А поскольку при каждом запуске корректора вызывается хотя бы одна функция интерпретации, по аналогии с леммой 2.1. получаем,
Граф-индикатор несовместности
Одним из ключевых понятий в описании ньютоновского корректора является понятие графа истории. Граф истории служит для сохранения информации о работе алгоритма М2В. Здесь описаны типы вершин графа истории и его возможные преобразования.
Граф истории является ориентированным графом с тремя основными типами вершин: вершинами-переменными, вершинами-числами и вершинами-функциями. Он построен на основе некоторого непустого графа задачи, который будем называть базовым для данного графа истории. Граф истории будем называть пустым, если базовый граф задачи пуст. С графом истории связан некоторый набор исторических функций.
Как будет описано далее, ньютоновский корректор сначала создает заготовку графа истории. Далее, на каждой итерации ньютоновского корректора будет выполняться преобразование, называемое исполнением вершины-отношения графа задачи с сохранением следа в графе истории. Как и в случае алгоритма М2В, исполняемая вершина-отношение называется рабочей вершиной итерации.
Если а - вершина графа истории, то time(a) - либо номер итерации корректора, на котором была создана вершина а, либо 0, если она присутствовала в заготовке. Будем говорить, что вершина а старше вершины Ь, если time(a) time(b). При этом будем говорить, что вершина b моложе вершины а. Если time(a) = time(b), будем говорить, что вершины а и Ь - ровесники. Вершины-числа делятся на константы и результаты. Вершины-константы в свою очередь делятся на границы и параметры. Среди вершин-границ выделим пару вершин, которые будем называть аргументами. Среди вершин-результатов будем выделять цели. Каждой вершине-константе соответствует ее значение. Тип значения вершины-константы - FP-число. Значение вершины-константы оговаривается при ее создании. Вершины-числа могут иметь как входящие, так и выходящие дуги. Вершины-переменные Множество вершин-переменных графа истории совпадает с множеством вершин-переменных базового графа задачи. Вершина, которая соответствует главной переменной корректора, называется главной вершиной-переменной. Вершины-переменные не имеют входящих дуг. Каждая вершина-переменная имеет две выходящие дуги, каждая из которых идет в некоторую вершину-границу. Выходящие дуги каждой вершины-переменной упорядочены.
Если v - вершина-переменная, то вершину-границу с такую, что (v,c) - первая выходящая дуга вершины v, обозначим l(v). Если v - вершина-переменная, то вершину-границу с такую, что (и, с) - вторая выходящая дуга вершины v, обозначим h(v).
По ходу пополнения графа истории некоторые дуги вида (v,c), где v - вершина-переменная, с - вершина-граница, удаляются и вместо них создаются новые. Для доказательства корректности предлагаемого алгоритма нам потребуется ссылаться на вершины-границы, в которые шли дуги из определенных вершин-переменных на определенной итерации ньютоновского корректора. В этом случае, если v - вершина-переменная, то вершину-границу с такую, что после к итераций корректора дуга (v,c) являлась первой выходящей дугой вершины v, обозначим lk(v).
Аналогично, если v - вершина-переменная, то вершину-границу с такую, что после к итераций корректора дуга (v,c) являлась второй выходящей дугой вершины и, обозначим hk{v). Каждая вершина-функция имеет связанную с ней функцию исполнения. Количество аргументов функции исполнения, связанной с данной вершиной-функцией, совпадает с количеством входящих в нее дуг. Входящие дуги каждой вершины-функции упорядочены. Если / - вершина-функция, вершину-число с такую, что дуга (с, /) существует и является г-й входящей дугой вершины /, обозначим Sj(/). Если / - вершина-функция, вершину-число с такую, что существует дуга (/,с), обозначим t(f). Вершины-функции делятся на вершины-условия, н-вершины и в-вершины. Пусть дан непустой граф задачи G, пусть Vi,...,vn - все его вершины-переменные, її, ...,/„- их значения. Пусть V{ - главная вершина-переменная, inf ij — sup Д. Тогда граф Н, содержащий вершины-переменные vi,..., vn графа G, вершины-границы Сі,... ,С2П со значениями inf Д, sup /lv .., inf /n, sup In, среди которых вершины Сгі-i и c-ii являются аргументами, дуги {vk,c2k-i), («fc,c2fc), где к = 1,...,п, при этом (t»fc,C2k_i) - первая, a {vk,c2k) -вторая выходящие дуги вершины vk, называется заготовкой графа истории, построенной на базе графа задачи G. Заготовка графа истории также является графом истории, она не имеет вершин-функций. В процессе работы корректора происходит эволюция графа истории. Изменения последнего возможны только в рамках сохранения следа в графе истории при исполнении рабочей вершины графа задачи. Для описания этого преобразования введем следующие промежуточные преобразования: создание следа границы; сохранение характеристики исполнения; сохранение условий корректности вершины-функции. Создание следа границы Пусть г - рабочая вершина текущей итерации корректора, v — t(r) - вершина-переменная, Н - граф истории, с - вершина-число. Преобразование создание следа нижней границы заключается в замене первой выходящей дуги вершины v дугой (и, с). Это преобразование будем обозначать trace(H,v, с). Преобразование создание следа верхней границы заключается в замене второй выходящей дуги вершины v дугой (г/, с). Его будем обозначать trace(H,v,c).
Сравнение алгоритмов NSC и М2В на примере конкретной задачи
Рассмотрим работу алгоритма М2В на графе G. От порядка выбора вершин для исполнения зависит время работы алгоритма М2В, поэтому в различных его реализациях применяются различные эвристики для выбора активной вершины. Сначала для простоты рассмотрим работу алгоритма при условии, что для исполнения выбирается активная вершина-отношение с минимальным индексом, а затем оценим, каким образом можно оптимизировать его работу.
Предположим, что на некоторой итерации алгоритма М2В, исключая первую, исполнилась вершина Г\. Опустим подробности, скажем лишь, что перед тем как вершина т\ будет исполнена следующий раз, будут последовательно исполнены вершины г2, г3, Г4, г3, г5, ге- Выделенным шрифтом отмечены те вершины, исполнение которых влечет изменение значения некоторой вершины-переменной.
Итак, если для исполнения выбирается активная вершина с минимальным индексом, то итерации алгоритма М2В можно разбить на группы, каждая из которых содержит итерации, последовательно исполняющие вершины гі, г2, г3, г4, г3, 7 5, /"б- При этом исполнение вершин г2, г3 и Гб не приводит к изменениям значений вершин-переменных.
Предположим, что некоторый эвристический алгоритм выбирает последовательно вершины гі, г4, Гб, и только если все они станут неактивными, выбирает одну из вершин гг, Гз, г$. В этом случае работа алгоритма М2В сведется к последовательному исполнению Несложно доказать, что при таком выборе вершины для исполнения время работы М2В будет минимальным. Предположим, что после некоторого исполнения вершины т ! значение, связанное с вершиной г і, стало [1 — 6,1], где S 1. Тогда при условии вычислений с бесконечной точностью, т.е. при условии выполнения FP = JR, после следующего исполнения вершины Г\ значение, связанное с вершиной V\, стало бы [1 — 6 + у, 1]. В нашем случае через 158915040 итераций значение вершины v\ станет равным [а, 1], где а Є FP таково, что ("а 2 ) = а а приблизительно равно 1 — 1.49 10 8. Этап 1. Создание подзадачи Итак, пусть переменная Х\ будет главной. Строим задачу М — (X, [0.5,0.5] х D2 х D3,C) и переходим ко второму этапу. Этап 2. Создание графа задачи и заготовки для графа истории По задаче М строим граф G . Он отличается от графа G, построенного в предыдущем пункте, только значением, связанным с вершиной v%. Значение вершины v\ в графе G -интервал [0.5,0.5]. На базе графа G строим заготовку графа истории Н. Она содержит три вершины- переменные г і, v2 и г»з, с которыми соответственно связаны значения [0.5,0.5], [0,оо] и [0,оо]; шесть вершин-границ blv .., Ъ, с которыми соответственно связаны значения 0.5, 0.5, 0, оо, 0 и сю; и дуги (vi, Ь2І-І), (fj, b2i), г = 1,..., 3, причем (t j, Ьгі-і) - первые выходящие дуги вершин «І, г = 1,..., 3. Вершины Ь\ и Ь2 будут аргументами. Далее активизируем все вершины-отношения графа G и переходим к третьему этапу. Этап 3. Пошаговое исполнение алгоритма М2В с сохранением информации в графе истории Здесь мы также будем выбирать для исполнения активную вершину-отношение с минимальным индексом. Итерация 1. т\ - рабочая вершина итерации. Исполняем Г\. вычисляем значение А = Р2([0.5,0.5] х [0,оо] П {(х,у)\х = у2). Получаем А = [0.25,0.25], но так как 0.25 Є FP, то новое значение вершины v2 - интервал [0.25,0.25]. В качестве нижней характеристики исполнения гх возьмем функцию - і(Уі5Ї/252/312/4) = У\- Ее упрощение - функция F\{y) = у2. Сохраняем характеристику исполнения 7 i для нижней границы значения v2: создаем вершину-результат с\, создаем н-вершину /і со связанной функцией исполнения Fx; создаем дуги (6і,Д) и (Л, ). 102 В качестве верхней характеристики исполнения гх возьмем функцию - 2(2/1,2/2,2/3,2/4) — у\- Ее упрощение - функция р2(у) = у2. Сохраняем характеристику исполнения гх для верхней границы значения г 2: создаем вершину-результат с2; создаем в-вершину /2 со связанной функцией исполнения F2; создаем дуги (62,/2) и (/г,с2). Набор условий корректности для функции F\ будет состоять из единственной функции з(2/ъ 2/2,2/3,2/4) = 2/ъ Действительно, если у Є [2/1,2/2] и уг О, то у2 27? Упрощение функции .Р3 - функция 3(2/) = у. Сохраняем условие корректности для Д: создаем вершину цель t\\ создаем вершину-условие /3 со связанной функцией исполнения F3; создаем дуги (Ьь/3), (/з, і) и ( і,еі). Набор условий корректности для функции F2 будет состоять из функций -4(2/1,2/2,2/3,2/4) = 2/2 и 5(2/1,2/2,2/3,2/4) = 2/1 +2/2- Действительно, если 2/ Є [j/i, у2], 2/2 О и 2/1 + 2/2 0, то 2/2 2/2- Упрощение F4 - функция F4(y) = у, упрощение Fs - функция - 5(2/1,2/2) =2/1 + 2/2- Сохраняем условия корректности для /2: создаем вершину-цель 2; создаем вершину-условие Д со связанной функцией исполнения F±\ создаем дуги (Ь2, Д), (/4,) и ( 2,с2). А также: создаем вершину-цель зі создаем вершину-условие /5 со связанной функцией исполнения Fs] создаем дуги (bi,/5), (62,/5), (/Б, З) И ( З,С2) Выполняем преобразования trace(H,v- _,ci) и trace(i/, г 2,с2): заменяем дугу (v2,b3) дугой (v2,cx); заменяем дугу (v2,64) дугой (v2,c2). Делаем активными вершины г і, г2, г3 и Г4 и неактивной вершину г і. Итерация 2. г2 - рабочая вершина итерации. Исполняем г2: вычисляем значение А — Рі([0.5,0.5] х [0,ею] П {(ж,у)[аг = у2). Получаем А = [0.5,0.5], значение вершины Vi не меняется. Делаем неактивной вершину г2. Итерация 3. гз - рабочая вершина итерации. Исполняем гз: вычисляем значение А — Pi([0.25,0.25] х [0, сю]П{(ж,2/)ж = у — 1). Получаем А = [0.25,0.25], значение вершины V2 не меняется. Делаем неактивной вершину г$. Итерация 4. г\ - рабочая вершина итерации. Исполняем г4: вычисляем значение А = Р2([0.25,0.25] х [0,сю] П {(х,у)\х = у - 1). Получаем А = [1.25,1.25], поскольку 1.25 Є FP, новое значение вершины v3 - интервал [1.25,1.25]. В качестве нижней характеристики исполнения г4 возьмем функцию -Рб(2/ъ2/2,2/3)2/4)2/5) = 2/1 + 2/5 с вектором параметров, состоящим из единственного элемента - числа 1. Ее упрощение - функция i 6 (2/1,2/2) = Ш + Уі- Сохраняем характеристику исполнения г4 для нижней границы значения г : создаем вершину-параметр рг со значением 1; создаем вершину-результат с3; создаем н-вершину /б со связанной функцией исполнения F ; создаем дуги (сь/6), (рі,/б) и (/в,с3). В качестве верхней характеристики исполнения г4 возьмем функцию V(2/i) 2/2,2/з, 2/4,2/5) = 2/2 + Уъ с вектором параметров, состоящим из единственного элемента - числа 1. Ее упрощение - функция 7(2/1,2/2) =2/1+2/2.
