оптимизация простой программы сэмплера Гиббса на Common Lisp

В качестве упражнения я переписал пример программы в блоге Семплер Гиббса на разных языках (повторно)Даррена Уилкинсона.

Код приведен ниже. Этот код выполняется на моем (5-летнем) компьютере примерно за 53 секунды, используя SBCL 1.0.56, создавая образ ядра с помощью buildapp, а затем запуская его с помощью

time ./gibbs > gibbs.dat

. Так как именно так были рассчитаны тайминги для других языков в пост, я думал, что сделаю что-то сопоставимое Код C в посте выполняется примерно за 25 секунд. Я хотел бы попытаться ускорить код Lisp, если это возможно.

##############################
gibbs.lisp
##############################
(eval-when (:compile-toplevel :load-toplevel :execute)
  (require :cl-rmath) (setf *read-default-float-format* 'double-float))

(defun gibbs (N thin)
  (declare (fixnum N thin))
  (declare (optimize (speed 3) (safety 1)))
  (let ((x 0.0) (y 0.0))
    (declare (double-float x y))
    (print "Iter x y")
    (dotimes (i N)
      (dotimes (j thin)
    (declare (fixnum i j))
    (setf x (cl-rmath::rgamma 3.0 (/ 1.0 (+ (* y y) 4))))
    (setf y (cl-rmath::rnorm (/ 1.0 (+ x 1.0)) (/ 1.0 (sqrt (+ (* 2 x) 2))))))
      (format t "~a ~a ~a~%" i x y))))

(defun main (argv)
  (declare (ignore argv))
  (gibbs 50000 1000))

Затем я собрал исполняемый файл gibbsс вызовом sh gibbs.shс gibbs.shas

##################
gibbs.sh
##################
buildapp --output gibbs --asdf-tree /usr/share/common-lisp/source/ --asdf-tree /usr/local/share/common-lisp/source/ --load-system cl-rmath --load gibbs.lisp --entry main

Я получаю 6 примечаний компилятора при компиляции с помощью SBCL. 1.0.56, которые воспроизведены ниже. Я не уверен, что с ними делать, но был бы признателен за любые подсказки.

; compiling file "/home/faheem/lisp/gibbs.lisp" (written 30 MAY 2012 02:00:55 PM):

; file: /home/faheem/lisp/gibbs.lisp
; in: DEFUN GIBBS
;     (SQRT (+ (* 2 X) 2))
; 
; note: unable to
;   optimize
; due to type uncertainty:
;   The result is a (VALUES (OR (DOUBLE-FLOAT 0.0) (COMPLEX DOUBLE-FLOAT))
;                           &OPTIONAL), not a (VALUES FLOAT &REST T).

;     (/ 1.0d0 (SQRT (+ (* 2 X) 2)))
; 
; note: unable to
;   optimize
; due to type uncertainty:
;   The second argument is a (OR (DOUBLE-FLOAT 0.0)
;                                (COMPLEX DOUBLE-FLOAT)), not a (COMPLEX
;                                                                DOUBLE-FLOAT).
; 
; note: forced to do static-fun Two-arg-/ (cost 53)
;       unable to do inline float arithmetic (cost 12) because:
;       The second argument is a (OR (DOUBLE-FLOAT 0.0) (COMPLEX DOUBLE-FLOAT)), not a DOUBLE-FLOAT.
;       The result is a (VALUES (OR (COMPLEX DOUBLE-FLOAT) (DOUBLE-FLOAT 0.0))
;                               &OPTIONAL), not a (VALUES DOUBLE-FLOAT &REST T).

;     (CL-RMATH:RGAMMA 3.0d0 (/ 1.0d0 (+ (* Y Y) 4)))
; 
; note: doing float to pointer coercion (cost 13)

;     (SQRT (+ (* 2 X) 2))
; 
; note: doing float to pointer coercion (cost 13)

;     (CL-RMATH:RNORM (/ 1.0d0 (+ X 1.0d0)) (/ 1.0d0 (SQRT (+ (* 2 X) 2))))
; 
; note: doing float to pointer coercion (cost 13)
; 
; compilation unit finished
;   printed 6 notes

; /home/faheem/lisp/gibbs.fasl written
; compilation finished in 0:00:00.073

ОБНОВЛЕНИЕ 1: Ответ Райнера Йосвигауказал, что аргумент SQRT может быть отрицательным, что был источником неясных заметок компилятора, которые я видел, а именно

The result is a (VALUES (OR (DOUBLE-FLOAT 0.0) (COMPLEX DOUBLE-FLOAT))
    ;                           &OPTIONAL), not a (VALUES FLOAT &REST T).

Компилятор жаловался, что, поскольку он не знал, было ли значение аргумента положительным, результатом может быть комплексное число. Поскольку в данном примере значение xявляется выборкой, отличающейся от гамма-распределения, он всегда больше 0.На это услужливо указал Стефан в списке рассылки пользователей SBCL, (см. второе сообщение в ветке «Оптимизация простой программы сэмплера Common Lisp Gibbs» , что это можно решить, объявив x большим или равным нулю следующим образом:

(declare (type (double-float 0.0 *) x))

См. Common Lisp Hyperspec для соответствующую документацию о Типы FLOAT и Обозначения интервалов.

Кажется, это немного ускоряет код. Теперь он надежно ниже 52 секунд, но все же не так уж много. Это также оставляет примечания о

note: doing float to pointer coercion (cost 13)

. Если это невозможно исправить по какой-то причине, я хотел бы знать, почему. Кроме того, объяснение того, что означает примечание, было бы интересно, несмотря ни на что. В частности, что здесь означает слово указатель? Связано ли это с тем, что вызываются функции C? Кроме того, стоимость 13 не кажется очень полезно. Что измеряется?

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