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

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

;This works for maxima 5.18. tested ok - knut gjerden
(let ()
  (setf maxima::*load-verbose* nil)
  (setf *debugger-hook* #'maxima::maxima-lisp-debugger)
  (let ((input-stream maxima::*standard-input*)
	(batch-flag nil))
    (progn
      (maxima::set-readtable-for-macsyma)
      (setf maxima::*read-default-float-format* 'double-float))
    (catch 'to-lisp
      (maxima::initialize-real-and-run-time)
      (intl::setlocale)
      (maxima::set-locale-subdir)
      (maxima::adjust-character-encoding)
      (maxima::set-pathnames)
      (when (boundp 'maxima::*maxima-prefix*)
	(push (pathname (concatenate 'string maxima::*maxima-prefix* "/share/locale/"))
	      intl::*locale-directories*))
      (setf (values input-stream batch-flag)
	    (maxima::process-maxima-args input-stream batch-flag))))
  )

; This works for maxima 5.13. not tested -knut gjerden 
;(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*))


(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 "/Users/knutgj/Code/usemaxima/infix.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)
    ,@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 (,@vars)
                ,(or declaration
                  `(declare (type double-float ,@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)


;This is the same as f1-static in normal lisp
;(loop :for i from 1 below 10000000 :sum
;      (- (* 2 (cos (* pi i))) (/ (sin (* pi i)) (* pi i))) )