Собственная линейная алгебра на лиспе 1

Под векторами ниже по тексту я буду подразумевать вектора в смысле линейной алгебры, а не в смысле стандартной структуры VECTOR в лиспе.

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

Будем писать векторные операции для трехмерного вектора с SINGLE-FLOAT в качестве типа каждой компоненты. Массивы в лиспе создаются функцией MAKE-ARRAY. Таким образом, чтобы создать новый вектор с координатами (50 60 70) нужно написать


(make-array 3 :element-type 'single-float :initial-contents (list 50.0 60.0 70.0))


Код прост и понятен, но слишком длинный, да и запомнить такое не так просто. Хотелось бы писать что-то вроде (vec3f 50.0 60.0 70.0), и чтоб компилятор вместо этого выражения сам бы подставлял make-array с нужными параметрами. (На самом деле, в лиспе принято функции и макросы создания данных называть начиная с «make-», то есть в моём случае должно было бы быть «(make-vec3f 50.0 60.0 70.0)», но из соображений лаконичности я оставлю vec3f.) Напрашивается такой вот макрос


(defmacro vec3f (x y z) "Создание вектора по координатам"
  `(make-array 3 :element-type 'single-float :initial-contents (list ,x ,y ,z)))


Теперь уже можно писать (vec3f 50.0 60.0 70.0). Для работы с вектором крайне полезно уметь узнавать его координаты. В лиспе для доступа к элементам массива используется функция ELT, т.е. получение x-координаты вектора v можно осуществить вызовом (elt v 0), y-координаты — (elt v 1), z-координаты — (elt v 2). Подобные записи не всегда являются наглядными, многие люди (и я в том числе) привыкли к «x», «y» и «z», поэтому напрашивается еще тройка макросов:


;;; Привычное обращение к координатам вектора

(defmacro vec3f-x (v)   `(elt ,v 0))
(defmacro vec3f-y (v)   `(elt ,v 1))
(defmacro vec3f-z (v)   `(elt ,v 2))


База написана, приступаем к содержательной части. Так как я хочу добиться быстродействия и не сильно беспокоюсь за размер скомпилированного кода, — буду пользоваться преимущественно макросами, чтобы компилятор разворачивал выражения в код, использующий только MAKE-ARRAY, ELT и арифметические функции. Возможно, что это ошибочный подход, — поживём-увидим.

При помощи написанного можно сходу написать пачку полезных операторов:


(defmacro mul3f (v s) "Покомпонентное умножение вектора v на скаляр s"
  `(vec3f (* (vec3f-x ,v) ,s) 
          (* (vec3f-y ,v) ,s) 
          (* (vec3f-z ,v) ,s)))

(defmacro div3f (v s) "Покомпонентное деление вектора v на скаляр s"
  `(vec3f (/ (vec3f-x ,v) ,s) 
          (/ (vec3f-y ,v) ,s) 
          (/ (vec3f-z ,v) ,s)))

(defmacro dot3f (a b) "Скалярное произведение векторов"
  `(+ (* (vec3f-x ,a) (vec3f-x ,b)) 
      (* (vec3f-y ,a) (vec3f-y ,b))
      (* (vec3f-z ,a) (vec3f-z ,b))))

(defmacro qlen3f (v)  "Квадрат длины вектора"  `(dot3f ,v ,v))
(defmacro len3f  (v)  "Длина вектора"          `(sqrt (qlen3f ,v)))
(defmacro norm3f (v)  "Нормирование вектора"  
  `(div3f ,v (len3f ,v)))


Теперь вызов (len3f (vec3f 30.0 40.0 0.0)) выдаст «50.0». Осталось самое сложное в этом посте: сложение и вычитание векторов.

Сложность заключается в том, что мы хотим написать макрос не для сложения двух векторов, а для сложения сколь угодного кол-ва векторов. Например, (add3f a b c d) должно сложить четыре вектора.

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

Во-вторых, для вычисления x-координаты результата нужно собрать все x-координаты каждого вектора в список и сложить их всех. Например, если нам на вход дали список векторов (a b c d), то для получения x-координаты результата мы должны просуммировать элементы списка ((vec3f-x a) (vec3f-x b) (vec3f-x c) (vec3f-x d)).

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


(map 'list (lambda(v) (vec3f-x v)) (list (vec3f 1.0 2.0 3.0) (vec3f 4.0 5.0 6.0) (vec3f 7.0 8.0 9.0)))


вернёт список (1.0 4.0 7.0) — будет применена функция под lambda к каждому элементу под list.

Теперь вернёмся к нашим сложениям. В макросе нужно получить список выражений вида «((vec3f-x v1) (vec3f-x v2) ...)» и подставить этот список под знак суммирования в «(+ ...)». Подстановка списка в выражение осуществляется конструкцией «,@». Приходим к такому вот:


(defmacro add3f (&rest vecs) "Сложение векторов"
  `(vec3f (+ ,@(map 'list (lambda(v) `(vec3f-x ,v)) vecs)) 
          (+ ,@(map 'list (lambda(v) `(vec3f-y ,v)) vecs)) 
          (+ ,@(map 'list (lambda(v) `(vec3f-z ,v)) vecs))))


Аналогично для разности:


(defmacro sub3f (&rest vecs) "Вычитание векторов"
  `(vec3f (- ,@(map 'list (lambda(v) `(vec3f-x ,v)) vecs)) 
          (- ,@(map 'list (lambda(v) `(vec3f-y ,v)) vecs)) 
          (- ,@(map 'list (lambda(v) `(vec3f-z ,v)) vecs))))


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


(macroexpand '(add3f a b c d e f))


получим:


(MAKE-ARRAY 3 :ELEMENT-TYPE 'SINGLE-FLOAT :INITIAL-CONTENTS
    (LIST (+ (VEC3F-X A) (VEC3F-X B) (VEC3F-X C) (VEC3F-X D) 
             (VEC3F-X E) (VEC3F-X F))
          (+ (VEC3F-Y A) (VEC3F-Y B) (VEC3F-Y C) (VEC3F-Y D) 
             (VEC3F-Y E) (VEC3F-Y F))
          (+ (VEC3F-Z A) (VEC3F-Z B) (VEC3F-Z C) (VEC3F-Z D) 
             (VEC3F-Z E) (VEC3F-Z F))))


Кроме того, VEC3F-X будет раскрыто в вызов ELT, поэтому в действительности во время выполнения будет выполнено ровно столько операций, сколько нужно, и ничего лишнего не будет.

Реализовать функцию cross3f, вычисляющую векторное произведение, оставлю на домашнее задание. На этом пока всё.