Efficient Lisp functions -- using results of symbolic calculations numerically

Dear developers,

one of the advantages of Common Lisp is that it blurs the distinction between
compile time and runtime because in many Lisp implementations you can compile
new functions to machine language speed at runtime. Why don't we use this in
Maxima? It would set apart Maxima very much from Mathematica and Maple where
you are forced to abhor loops and iterations because the interpreted languages
make them several orders of magnitude slower than one is accustomed to. There
is a compile command in Mathematica, but it only creates byte code which is
still rather slow, not machine code.

I tend to use Maxima from Lisp (SBCL) directly: After analysing the source
code, I found the following workaround. I start sbcl normally without the
Maxima core and execute this:

------------------------------------------------------------
(defpackage :my-package (:use :common-lisp))
(in-package :my-package)

(let ()
  (load "/home/daniel/pakete/maxima-5.13.0/lisp-utils/defsystem.lisp")
  ;; This creates a lot of warnings in SBCL, but they seem to do no harm:
  (load "/home/daniel/pakete/maxima-5.13.0/src/maxima.system")
  (funcall (intern (symbol-name :operate-on-system) :mk) "maxima" 
  :load :verbose t))

(let ()
  (setf maxima::*load-verbose* nil)
  (maxima::initialize-real-and-run-time)
  (maxima::set-locale)
  (maxima::set-pathnames)
  (if (not maxima::*maxima-quiet*) (maxima::maxima-banner))
  (setq maxima::*maxima-started* t))

(defparameter *maxima-package* (find-package :maxima))

(defun maxima::intern-invert-case (string)
  (intern (maxima::maybe-invert-string-case string) *maxima-package*))
------------------------------------------------------------
  
I had to overwrite intern-invert-case in order to be able to use Maxima from
another package (The parsing of if statements didn't work with the original
version of intern-invert-case, because the keywords were interned in the wrong
package).

Then I wrote some functions for translating normal Lisp code into the Maxima
representation of formulas (The code is rather in a proof of concept stage, so
it doesn't cover all cases, yet):

------------------------------------------------------------
(defparameter *maxima-translation-hash* (make-hash-table :test #'eql))
(defparameter *maxima-i* #C(0d0 1d0))
(defparameter *maxima-e* (exp 1d0))

(defmacro href (hash key)
  `(gethash ,key ,hash))

(defun init-maxima-translation-hash ()
  (let ((initializations '((maxima::mplus +) (maxima::mminus -)
                           (maxima::mtimes *) (maxima::rat /)
                           (maxima::mexpt expt) (maxima::mlist list)
                           (maxima::mequal =) (maxima::mgreaterp >)
                           (maxima::mabs abs)
                           (maxima::%sin sin) (maxima::%cos cos)
                           (maxima::%tan tan) (maxima::%exp exp)
                           (maxima::%log log) (maxima::%sqrt sqrt)
                           (maxima::%atan atan)
                           (maxima::%sum sum)
                           (maxima::%derivative %diff)
                           )))
    (loop :for (a b) in initializations :do
       (setf (href *maxima-translation-hash* a) b)
       (setf (href *maxima-translation-hash* b) a))))

(init-maxima-translation-hash)

(defun concatenate-string (&rest args)
    "string concatenation"
  (apply #'concatenate 'string args))

;; This is the infix-package of page http://www.cliki.net/infix
(load "/path/to/infix2.lisp")
  
(defun translate-to-maxima (expr)
    "Translate a lisp expression to a maxima expression."
  (cond ((atom expr)
         expr)
        ((eql (car expr) 'nfx)
         (translate-to-maxima (macroexpand-1 expr)))
        ((eql (car expr) '/)
         (cond ((= (length expr) 3)
                (translate-to-maxima
                 `(* ,(second expr) (expt ,(third expr) -1))))
               ((= (length expr) 2)
                (translate-to-maxima `(expt ,(second expr) -1)))
               (t
                (error "unsupported expr: ~s" expr))))
        (t
         (let* ((head (first expr))
                (head-translated (href *maxima-translation-hash* head)))
           (when (null head-translated)
             (let ((translation
                    (intern (concatenate-string "$" (string head))
                            *maxima-package*)))
               (setf (href *maxima-translation-hash* head) translation)
               (setq head-translated translation)))
           (apply #'maxima::mfuncall
                  (cons head-translated
                        (mapcar #'translate-to-maxima (cdr expr))))))))


(defun translate-from-maxima (expr)
    "Translate a maxima expression to a lisp expression."
  (cond
    ((atom expr)
     (case expr
       (maxima::$%i *maxima-i*)
       (maxima::$%e *maxima-e*)
       (maxima::$%pi pi)
       (t expr)))
    ((consp (car expr))
     (let ((operator (caar expr)))
       (cond
         ((and (eql operator 'maxima::mexpt)
               (numberp (third expr)) (minusp (third expr)))
          (assert (= (length expr) 3))
          (if (eql (third expr) -1)
              (list '/ (translate-from-maxima (second expr)))
              (list '/ (translate-from-maxima
                        (list (first expr) (second expr) (- (third expr)))))))
         ((and (eql operator 'maxima::mexpt) (eql (third expr) 1))
          (assert (= (length expr) 3))
          (translate-from-maxima (second expr)))
         ((and (eql operator 'maxima::mexpt)
               (equal (third expr) '((maxima::rat maxima::simp) 1 2)))
          (assert (= (length expr) 3))
          (list 'sqrt (translate-from-maxima (second expr))))
         ((and (eql operator 'maxima::mtimes) (eql (second expr) -1))
          (if (= (length expr) 3)
              (list '- (translate-from-maxima (third expr)))
              (list '- (translate-from-maxima
                        (list* (first expr) (cddr expr))))))
         (t
          (let ((head-translated (href *maxima-translation-hash* operator)))
            (when (null head-translated)
              ;; {ja} Try removing % or $.
              (error "There's no backtranslation for ~s." operator))
            (cons head-translated
                  (mapcar #'translate-from-maxima (cdr expr))))))))
    (t (error "not implemented"))))
------------------------------------------------------------

In order to make the expressions in Lisp more readable for the mathematically
inclined, I use the infix package of page http://www.cliki.net/infix.

Now we can calculate a formula symbolically at compile time or runtime using
Maxima, compile it in order to boost it up to machine language speed and then
call it millions of times:

------------------------------------------------------------
(defmacro with-gensyms ((&rest vars) &body body)
  "for facilitating defmacro"
  `(let ,(mapcar #'(lambda (var)
                     `(,var (gensym ,(string var))))
                 vars)
    , at body))

(defun lambda-by-maxima-1 (vars expr &key declaration)
  "Use the result of a maxima expr to create a lambda function."
  (let* ((maxima-expr (print (translate-to-maxima expr)))
         (result-expr (translate-from-maxima maxima-expr)))
    (format t "~&The maxima result of ~%~s is ~%~s." expr result-expr)
    ;; The print statements are for debugging purposes. They pass the result on.
    (print `#'(lambda (, at vars)
                ,(or declaration
                  `(declare (type double-float , at vars) (values double-float)
                    (optimize (speed 3) (space 0) (safety 0) (debug 0))))
                ,result-expr))))

(defun lambda-by-maxima-dynamic (vars expr &key declaration)
  "Use the result of a maxima expr to create a lambda function."
  (eval (lambda-by-maxima-1 vars expr :declaration declaration)))

(defmacro lambda-by-maxima-static ((&rest vars) expr &key declaration)
  "Use the result of a maxima expr to create a lambda function."
  (lambda-by-maxima-1 vars expr :declaration declaration))
------------------------------------------------------------

Here's an example in which Maxima calculates the derivative of a function and
sums up 10 million of the result's function applications. In the dynamic
version, the exponent and the sine function in the formula are only provided
at runtime:

------------------------------------------------------------
(defun f1-static ()
  (let ((f2 (lambda-by-maxima-static (x) (nfx (diff ((sin (x ** 2)) / x) x)))))
    (loop :for i from 1 below 10000000 :sum
       (funcall f2 (nfx (sqrt (pi * i)))))))

(defun f1-dynamic (func exponent)
  (let ((f2 (lambda-by-maxima-dynamic
             `(x) `(nfx (diff ((,func (x ** ,exponent)) / x) x)))))
    (loop :for i from 1 below 10000000 :sum
      (funcall f2 (nfx (sqrt (pi * i)))))))

(f1-static)
(f1-dynamic 'sin 2)
------------------------------------------------------------

Both take about a second on my machine (using SBCL 1.0.5).

This way, we can use results of symbolic calculations in numerical algorithms
without any performance penalty caused by the computer algebra system. And we
can do this many times while the application is running, i.e. we can keep
switching back and forth from general and powerful symbolic calculations to
highly efficient numerical calculations in Lisp.

Best regards,
  Bruno Daniel