share/contrib/cgrind.lisp

   1 (in-package :maxima)
   2
   3 (macsyma-module cgrind)
   4
   5 ;; This is heavily lifted from grind.lisp and fortra.lisp, and should have the
   6 ;; same features as the fortran command. In order to work it needs to be compiled and
   7 ;; then loaded, as in (for the case of cmucl):
   8 ;; :lisp (compile-file "cgrind.lisp"), copy cgrind.sse2f to ~/.maxima
   9 ;; load(cgrind)
  10 ;; Then one can run cgrind(expression) or cgrind(matrix).
  11 ;; M. Talon (2011)
  12
  13
  14 (declare-top ($loadprint ;If NIL, no load message gets printed.
  15         ))
  16
  17 ;; This function is called from Macsyma toplevel.  First the arguments is
  18 ;; checked to be a single expression. Then if the argument is a
  19 ;; symbol, and the symbol is bound to a matrix, the matrix is printed
  20 ;; using an array assignment notation.
  21
  22 (defmspec $cgrind (l)
  23   (setq l (fexprcheck l))
  24   (let ((value (strmeval l)))
  25     (cond ((msetqp l) (setq value `((mequal) ,(cadr l) ,(meval l)))))
  26     (cond ((and (symbolp l) ($matrixp value))
  27            ($cgrindmx l value))
  28           ((and (not (atom value)) (eq (caar value) 'mequal)
  29                 (symbolp (cadr value)) ($matrixp (caddr value)))
  30            ($cgrindmx (cadr value) (caddr value)))
  31           (t (c-print value)))))
  32
  33 (defun c-print (x &optional (stream *standard-output*))
  34   ;; Restructure the expression for displaying.
  35   ;; Mainly sanitizes exponentials, notably exp(2/3) becomes
  36   ;; exp(2.0/3.0)
  37
  38   (setq x (scanforc x))
  39
  40   ;; Protects the modifications to mexpt from creeping out.
  41
  42   (unwind-protect
  43
  44     (progn
  45     (defprop mexpt msz-cmexpt grind)
  46
  47     ;; This means basic printing for atoms, grind does fancy things.
  48     (setq *fortran-print* t)
  49
  50     ;; Prints using the usual grind mechanisms
  51     (mgrind x stream)(write-char #\; stream)(write-char #\Newline stream))
  52
  53     ;; Restore usual mexpt property etc. before exiting this frame.
  54     (defprop mexpt msz-mexpt grind)
  55     (setq *fortran-print* nil))
  56   '$done)
  57
  58
  59 ;; The only modification to grind, converts a^b to pow(a,b), but taking
  60 ;; care of appropriate bracketing. The argument l to the left of (MEXPT)
  61 ;; has to be composed backwards. Finally a^-b has special treatment.
  62
  63 (defun msz-cmexpt (x l r)
  64   (setq l (msize (cadr x) (revappend '(#\p #\o #\w #\() l) (list #\,) 'mparen 'mparen)
  65         r (if (mmminusp (setq x (nformat (caddr x))))
  66               (msize (cadr x) (list #\-) (cons #\) r) 'mexpt rop)
  67               (msize x nil  (cons #\) r ) 'mparen 'mparen)))
  68   (list (+ (car l) (car r)) l r))
  69
  70
  71
  72
  73
  74 ;; Takes a name and a matrix and prints a sequence of C assignment
  75 ;; statements of the form
  76 ;; NAME[I][J] = <corresponding matrix element>
  77 ;; This requires some formatting work unnecessary for the fortran case.
  78
  79 (defmfun $cgrindmx (name mat &optional (stream *standard-output*) &aux ($loadprint nil))
  80   (cond ((not (symbolp name))
  81          (merror (intl:gettext "cgrindmx: first argument must be a symbol; found: ~M") name))
  82         ((not ($matrixp mat))
  83          (merror (intl:gettext "cgrindmx: second argument must be a matrix; found: ~M") mat)))
  84   (do ((mat (cdr mat) (cdr mat)) (i 1 (1+ i)))
  85       ((null mat))
  86     (do ((m (cdar mat) (cdr m)) (j 1 (1+ j)))
  87         ((null m))
  88        (format stream "~a[~a][~a] = " (string-left-trim "$" name) (1- i) (1- j) )
  89        (c-print (car m) stream)))
  90   '$done)
  91
  92
  93
  94
  95
  96 ;; This C scanning function is similar to fortscan.  Prepare an expression
  97 ;; for printing by converting x^(1/2) to sqrt(x), etc.   Since C has no
  98 ;; support for complex numbers, contrary to Fortran, ban them.
  99
 100 (defun scanforc (e)
 101   (cond ((atom e) (cond ((eq e '$%i) ;; ban complex numbers
 102                          (merror (intl:gettext "Take real and imaginary parts")))
 103                         (t e)))
 104         ;; %e^a -> exp(a)
 105         ((and (eq (caar e) 'mexpt) (eq (cadr e) '$%e))
 106          (list '(%exp simp) (scanforc (caddr e))))
 107         ;; a^1/2 -> sqrt(a)  1//2 is defined as ((rat simp) 1 2)
 108         ((and (eq (caar e) 'mexpt) (alike1 (caddr e) 1//2))
 109          (list '(%sqrt simp) (scanforc (cadr e))))
 110         ;; a^-1/2 -> 1/sqrt(a)
 111         ((and (eq (caar e) 'mexpt) (alike1 (caddr e) -1//2))
 112          (list '(mquotient simp) 1 (list '(%sqrt simp) (scanforc (cadr e)))))
 113         ;; (1/3)*b -> b/3.0 and (-1/3)*b -> -b/3.0
 114         ((and (eq (caar e) 'mtimes) (ratnump (cadr e))
 115               (member (cadadr e) '(1 -1) :test #'equal))
 116          (cond ((equal (cadadr e) 1) (scanforc-mtimes e))
 117                (t (list '(mminus simp) (scanforc-mtimes e)))))
 118         ;; 1/3 -> 1.0/3.0
 119         ((eq (caar e) 'rat)
 120          (list '(mquotient simp) (float (cadr e)) (float (caddr e))))
 121         ;; rat(a/b) -> a/b via ratdisrep
 122         ((eq (caar e) 'mrat) (scanforc (ratdisrep e)))
 123         ;; ban complex numbers
 124         ((and (member (caar e) '(mtimes mplus) :test #'eq)
 125               (let ((a (simplify ($bothcoef e '$%i))))
 126                 (and (numberp (cadr a))
 127                      (numberp (caddr a))
 128                      (not (zerop1 (cadr a)))
 129                      (merror (intl:gettext "Take real and imaginary parts"))))))
 130         ;; in general do nothing, recurse
 131         (t (cons (car e) (mapcar 'scanforc (cdr e))))))
 132
 133 ;; This is used above 1/3*b*c -> b*c/3.0
 134 (defun scanforc-mtimes (e)
 135   (list '(mquotient simp)
 136         (cond ((null (cdddr e)) (scanforc (caddr e)))
 137               (t (cons (car e) (mapcar 'scanforc (cddr e)))))
 138         (float (caddr (cadr e)))))
 139