Тем не менее, минимализм Лиспа очарователен. Весьма несложно разобраться в синтаксисе языка с семью операторами. Сделать собственную релизацию — тоже довольно просто. Увы, в конечном счете, именно эта простота и послужила Лиспу дурную службу. Слабая совместимость реализаций и откровенно непривычный подход к программированию, не говоря уже о скобках (с непривычки скобки действительно ломают мозг (тем более вложенные)) — это как раз то, что составляет общественное мнение о языке. Common Lisp, собственно, и есть попытка привести множество реализаций под один знаменатель, обеспечить более привычные и удобные идиомы (оставив в языке скобки (к скобкам быстро привыкаешь)). Но сейчас не про него. На примере программки для нахождения вещественных корней полиномиального уравнения я бы хотел показать, что и минималистичный Лисп тоже заслуживает право на жизнь.
В качестве реализации я возьму lis.py Питера Норвига. Это 90 строчный питоновский скрипт, который обеспечивает минималистичное подмножество диалекта Scheme на три десятка операторов. Это включая арифметику, логику, предикаты и даже излишества вроде length и list.
Ближе к задаче. Совсем несложно найти решение уравнения численно, если оно одно и находится на интервале, на котором функция монотонна. Любой способ вроде «половинного деления» или Ньютона-Рафсона подойдет. Совсем другое дело, если таких решений несколько. Полиномиальное уравнение n степени имеет не более чем n вещественных корней. Например, x^2-1=0 имеет два корня, x^2=0 — один, а x^2+1=0 — ни одного. Если n нечетное, то уравнение всегда имеет хотя бы один корень.
Мы знаем про полиномы две вещи — они непрерывны, и производная полинома — тоже полином. Ну, либо константа в экстремальном случае. Интервалы монотонности ограничены сменой знака производной. Смена знака непрерывной производной в какой либо точке означает, что в этой точке она проходит через ноль. То есть интервалы монотонности полинома ограничены корнями другого полинома — его производной. Из этого получится отличный пример рекурсивной функции.
(defun pol_roots (pol_ns) (if (= (length pol_ns) 2) (list (/ (- 0 (car pol_ns)) (car (cdr pol_ns)))) ((lambda (dp_roots) (if (null dp_roots) (find_single_root pol_ns) (pol_root_l pol_ns dp_roots))) (pol_roots (d pol_ns))))) Что тут происходит. Предполагается, что полином для решения отдается в виде его коэффициентов от младшей степени к старшей. То есть 3x^2+2x+1 = (1 2 3). Если длина списка коэффициентов полинома равна 2, то это уравнение прямой ax+b=0, и его можно решив просто разделив -b на a. Так как лисп не Common, удобностей вроде first и second не предвидится. Первый элемент возьмем в помощью (car, второй — (car (cdr. Это «голова» и «голова хвоста», соответственно. Еще одна неприятность — в этой реализации унарного минуса нет, так что вместо него придется делать (- 0.
Если же список коэффициентов полинома длиннее 2, придется решать его уже по-серьезному. Для начала найдем корни его производной и положим их в локальную переменную dp_roots. Однако let у нас нет. Можно бы сделать begin, чтобы начать блок, а в блоке использовать define, но данная реализация Лиспа не различает контекста переменной. Она не будет локальной. Поэтому мы делаем хитрость: определяем функцию с параметром dp_roots и тут же выполняем ее со значением (pol_roots (d pol_ns)), где d — функция, находящая производную полинома.
Если список корней производной пустой, это означает две вещи: степень производной — четная, степень полинома соответственно нет. То есть полином имеет один корень. Не более, потому что монотонный интервал не делится; не менее, потому что степень нечетная. Вызовем функцию для нахождения единственного корня: find_single_root.
А вот если список корней производной не пустой, то есть полином имеет несколько интервалов монотонности, мы ищем корни на каждом интервале с помощью вспомогательной функции pol_root_l.
(define pol_root_l (lambda (pol_ns xs) (append (root_a (lambda (x) (pol pol_ns x)) (car xs) -1.0) (pol_root_r pol_ns xs)))) Функция принимает список коэффициентов полинома и список корней производной. Она склеивает список корней, полученной от функции root_a для самого левого интервала, мы еще до нее дойдем, и pol_rool_r.
(define pol_root_r (lambda (pol_ns xs) (if (= (length xs) 1) (root_a (lambda (x) (pol pol_ns x)) (car xs) 1.0) (append (root_ab (lambda (x) (pol pol_ns x)) (car xs) (car (cdr xs))) (pol_root_r pol_ns (cdr xs)))))) Функция pol_root_r находит корни на всех интервалах, кроме левого. Если интервал только один, список корней находится функцией root_a. Естественно, это список из одного элемента максимум. Но он также может быть пустым, так что список тут используется для того же, для чего и монада Maybe, или шарповый нулябельный элемент.
Если интервалов больше чем один, то результат работы функции складывается из корней (одного или ни одного), найденных на первом интервале функцией root_ab, и результата работы себя же на всех остальных интервалах.
(define root_ab (lambda (fx a b) ((lambda (c) (if (= (sign (fx a)) (sign (fx b))) (list ) (if (< (abs (fx c)) EPS) (list c) (if (= (sign (fx c)) (sign (fx a))) (root_ab fx c b) (root_ab fx a c))))) (/ (+ a b) 2.0)))) Функция root_ab — это собственно и есть численная решалка. Она находит корень fx на интервале (a b). Здесь используется просто метод половинного деления, но его можно легко заменить любым другим более эффективным методом.
Если знак функции в первой точке интервала совпадает со знаком в последней, то, исходя из условия монотонности, корня нет: график функции числовую ось не пересекает. Возвращаем пустой list.
Если же знаки разные, то проверяем, может быть, середина интервала нам уже сгодится в качестве корня. Уловия могут быть самые разные, тут мы просто сравниваем результат функции в этой середине с некоторой достаточно маленькой величиной. Если абсолютное значение функции при некотором значении аргумента меньше, чем некое наперед заданное число, то этот аргумент «достаточно корень», чтобы больше ничего не искать. Если недостаточно, то мы сравниваем знак в середине интервала со знаками в начале и конце. Там, где эти знаки разнятся, там наш корень и лежит — запускаем root_ab уже для этого суженного интервала.
(define root_a (lambda (fx a d) (if (= (sign (- (fx (+ a d)) (fx a))) (sign (fx a))) (list ) (if (= (sign (fx (+ a d))) (sign (fx a))) (root_a fx (+ a d) (* d 2)) (root_ab fx a (+ a d)))))) Функция root_a выглядит замысловато, но только потому что тут много условий со знаками. Она находит корень функции fx на интервале либо от минус бесконечности до a, либо от a до плюс бесконечности в зависимости от d. Число d — это первый шаг поиска. Если оно положительное, корень будет искатсья справа от a, если отрицательное — слева.
Первое условие проверяет, есть ли на интервале корень. Мы же договорились, что функция монотонна, поэтому, например, если функция положительная, ищем мы корень справа от a, а функция возрастает, то, конечно, никакого корня на этом интервале не будет. Тут есть некоторый изъян в логике, так как для того, чтобы утверждать обратное при убывающей функции, надо еще быть уверенным, что ее предел в бесконечности меньше нуля. Одна монотонность этого не гарантирует. Но для полиномов это, конечно же так. Предел полинома в бесконечности каким-либо числом быть не может.
Для того, чтобы узнать, возрастает или убывает функция, мы просто берем разность значений в a и a+d. В общем случае так делать нельзя, надо бы посмотреть на производную, но полиномы, монотонность, бла-бла-бла — можно. Мы сравниваем это значение со значением функции в a. Если корня нет, возвращаем пустой список.
А вот если есть, делаем простую вешь: проверям, есть ли корень на интервале (a, a+d). Если есть, найдем его root_ab — и дело с концом. А если нет, то он просто где-то дальше. Мы увеличиваем d в два раза и продолжаем поиски на (a+d, a+2*d). Потом, получается, на (a+d*2, a+4*d), потом (a+4*d, a+8*d), (a+8*d, a+16*d) и так далее. Получается эдакий банрный поиск навыворот.
Имея такую замечательную функцию, вернемся к find_single_root.
(define find_single_root (lambda (pol_ns) ((lambda (fx) (append (root_a fx 0 EPS) (root_a fx 0 (- 0 EPS)))) (lambda (x) (pol pol_ns x))))) Очень простая функция. Она находит единственный корень на монотонном интервале от минус до плюс бесконечности путем склеивания корней (опять же, либо одного, либо ни одного) справа от нуля с корнями слева.
Ну, с алгоритмом вроде разобрались. Остались функции вспомогательные. Например, для нахождения производной.
Производная полинома находится очень просто. Для каждого члена, кроме a0: (ax^n)’ = anx^(n-1). На Common Lisp это можно было бы записать в одну строчку: (loop for n in (cdr pol_ns) and i from 1 collect (* n i))) То есть: берем хвост коэффициентов и каждый элемент множим на его же позицию начиная с 1. Но у нас нет удобных loop, вообще никаких нет. Более того, у нас даже нет goto, так что придется обойтись без итерации вообще.
К счастью, это совсем не сложно.
(define di (lambda (ns i) (if (null? ns) (quote ()) (append (list (* (car ns) i)) (di (cdr ns) (+ 1 i)))))) (define d (lambda (ns) (di (cdr ns) 1))) Определяем функцию со счетчиком индекса коэффициента и склеиваем элементы списка ns, помноженные на этот счетчик, один за одним. Эту функцию оборачиваем в красивую d с одним аргументом, которая режет хвост коэффициентам полинома и отдает его нашей считалке с индексом. Вообще, итерацию всегда можно превратить в рекурсию без особых сложностей, единственное, однострочник из цикла все-таки выразительнее.
Похожая ситуация и с вычислением полинома. В Scheme есть удобная функция iota — генератор списка-последовательности. (iota 3) => (0 1 2), (iota 3 1) => (1 2 3), (iota 3 1 2) => (1 3 5). Первый аргумент — сколько чисел генерировать, второй, начиная от какого, третий — с каким шагом. То есть для вычисления полинома в Scheme нам достаточно сделать что-то вроде: (apply + (map pol_ns (iota (length pol_ns) 1))).
Увы, в нашем подмножестве нет ни iota, ни даже map. Так что снова рекурсия.
(define pol (lambda (ns x) (if (null? ns) 0.0 (+ (car ns) (* x (pol (cdr ns) x)))))) Еще в нашем языке нет abs и sign.
(define sign (lambda (x) (if (< x 0) -1 1))) (define abs (lambda (x) (* (sign x) x))) Ну и EPS, конечно же.
(define EPS 0.00001) Заметьте, что если в Common Lisp пространства имен функций и переменных разделены: нам нужны defun и defvar, не говоря уж о defparameter и defconstant, то в нашем минималистичном Лиспе — нет. И то, и другое определяется в одном глобальном контексте посредством define. Это имеет смысл, так как Лиспу традиционно все равно, код перед ним или данные: он может как процессить код, так и выполнять данные без особых проблем. Но при написании больших вещей проще, наверное, все-таки использовать везде #' и funcall, чем путаться в кододанных.
Выводы
Минималистичный Лисп очень по-своему красивый язык. Очень по-своему. Самое замечательное: его синтаксические правила можно записать на листике целиком, выучить за пятнадцать минут, забыть за двадцать, и все равно писать код более-менее уверенно. Самое неприятное, код получается ну… не очень идеоматичный. И не очень лаконичный. И не очень, как следствие, понятный. Все-таки хорошо, когда под рукой есть набор примитивов не только для удобного написания кода, но еще и для удобного чтения его чужими людьми в последствии.
ссылка на оригинал статьи http://habrahabr.ru/post/190848/
Добавить комментарий