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