Hello,
I have slightly modified the cgrind program for grinding maxima expressions
to a form acceptable by C, so that when complex numbers appear in the
expression an error message is produced. Last version at:
http://www.lpthe.jussieu.fr/~talon/cgrind.html
or below signature.
I doing that i have seen what i beleive is a bug in fortra.lisp, the lack of
a progn below unwind-protect which seems to nullify the protection.
cgrind works like that:
niobe% maxima
(%i1) load(cgrind);
(%o1) /home/michel/.maxima/cgrind.sse2f
(%i2) cgrind(integrate((x^3+7*x)/sqrt(x^2+1),x));
pow(x,2)*sqrt(pow(x,2)+1)/3.0e+0+1.9e+1*sqrt(pow(x,2)+1)/3.0e+0;
(%o2) done
(%i3) m:matrix([x, x^2], [x^(-2),9]);
[ 2 ]
[ x x ]
[ ]
(%o3) [ 1 ]
[ -- 9 ]
[ 2 ]
[ x ]
(%i4) cgrind(m);
M[0][0] = x;
M[0][1] = pow(x,2);
M[1][0] = 1/pow(x,2);
M[1][1] = 9;
(%o4) done
Now i am working on colnew, which has been ported to maxima by R. Toy. In
fact it works remarkably well. I cannot see much to be done to improve stuff
except some dumbing down pre parametrizations. But i see at least one
problem, if one wants to do continuations, the way in which colnew:colnew is
called under a lot of flet, let, etc. stuff makes it hard to modify the
call. Let me recall that continuations consist in doing further calls to
colnew keeping the old fspace and ispace, but moving slowly an optional
parameter. This way the old solution serves as a guess for the new one. This
is a very powerful feature very well exemplified in
http://www.mathworks.com/help/techdoc/math/f1-713877.html#brfhdqp-2
Colnew has provisions to do that by setting ipar(9)= 2 or 3
except for the first computation where ipar(9)=1 and a guess is provided.
Well, i dont see easy way to do that with the way the wrapper is presently
formulated. I wonder if it would not be more expedient to put the wrapper
in the colnew package, so that names don't need any more to be hidden under
flets and lets, add another parameter (the continuation parameter) to
differential equations and boundary conditions, and write a smaller wrapper
in the maxima package.
Anyaways, in doing that, i have seen two small problems.
First, there is a function "push" documented in maxima manual, but
it doesn't work:
niobe% maxima
(%i1) push(1,[]);
(%o1) push(1, [])
Second, there exists on maxima CVS a file called colnew.mac under
share/colnew, which contains
/* Compile and load the colnew routines (in lisp) */
load("load-colnew.lisp");
At least in my installation this file is missing. But it is useful for
loading colnew easily like:
niobe% maxima
(%i1) load(colnew);
;; Loading #P"/home/michel/.maxima/binary/binary-cmucl/share/colnew/binary-
cmucl/colnew-package.sse2f".
......
;; Loading #P"/home/michel/.maxima/binary/binary-cmucl/share/colnew/binary-
cmucl/colnew-if.sse2f".
(%o2) /usr/local/share/maxima/5.22.1/share/colnew/colnew.mac
With my best regards
--
Michel Talon
http://www.lpthe.jussieu.fr/~talon/cgrind.lisp
(in-package :maxima)
(macsyma-module cgrind)
;; This is heavily lifted from grind.lisp and fortra.lisp, and should have
the
;; same features as the fortran command. In order to work it needs to be
compiled and
;; then loaded, as in (for the case of cmucl):
;; :lisp (compile-file "cgrind.lisp"), copy cgrind.sse2f to ~/.maxima
;; load(cgrind)
;; Then one can run cgrind(expression) or cgrind(matrix).
;; M. Talon (2011)
(declare-top ($loadprint ;If NIL, no load message gets printed.
))
;; This function is called from Macsyma toplevel. First the arguments is
;; checked to be a single expression. Then if the argument is a
;; symbol, and the symbol is bound to a matrix, the matrix is printed
;; using an array assignment notation.
(defmspec $cgrind (l)
(setq l (fexprcheck l))
(let ((value (strmeval l)))
(cond ((msetqp l) (setq value `((mequal) ,(cadr l) ,(meval l)))))
(cond ((and (symbolp l) ($matrixp value))
($cgrindmx l value))
((and (not (atom value)) (eq (caar value) 'mequal)
(symbolp (cadr value)) ($matrixp (caddr value)))
($cgrindmx (cadr value) (caddr value)))
(t (c-print value)))))
(defun c-print (x &optional (stream *standard-output*))
;; Restructure the expression for displaying.
;; Mainly sanitizes exponentials, notably exp(2/3) becomes
;; exp(2.0/3.0)
(setq x (scanforc x))
;; Protects the modifications to mexpt from creeping out.
(unwind-protect
(progn
(defprop mexpt msz-cmexpt grind)
;; This means basic printing for atoms, grind does fancy things.
(setq *fortran-print* t)
;; Prints using the usual grind mechanisms
(mgrind x stream)(write-char #\; stream)(write-char #\Newline stream))
;; Restore usual mexpt property etc. before exiting this frame.
(defprop mexpt msz-mexpt grind)
(setq *fortran-print* nil))
'$done)
;; The only modification to grind, converts a^b to pow(a,b), but taking
;; care of appropriate bracketing. The argument l to the left of (MEXPT)
;; has to be composed backwards. Finally a^-b has special treatment.
(defun msz-cmexpt (x l r)
(setq l (msize (cadr x) (revappend '(#\p #\o #\w #\() l) (list #\,)
'mparen 'mparen)
r (if (mmminusp (setq x (nformat (caddr x))))
(msize (cadr x) (list #\-) (cons #\) r) 'mexpt rop)
(msize x nil (cons #\) r ) 'mparen 'mparen)))
(list (+ (car l) (car r)) l r))
;; Takes a name and a matrix and prints a sequence of C assignment
;; statements of the form
;; NAME[I][J] = <corresponding matrix element>
;; This requires some formatting work unnecessary for the fortran case.
(defmfun $cgrindmx (name mat &optional (stream *standard-output*) &aux
($loadprint nil))
(cond ((not (symbolp name))
(merror (intl:gettext "cgrindmx: first argument must be a symbol;
found: ~M") name))
((not ($matrixp mat))
(merror (intl:gettext "cgrindmx: second argument must be a
matrix; found: ~M") mat)))
(do ((mat (cdr mat) (cdr mat)) (i 1 (1+ i)))
((null mat))
(do ((m (cdar mat) (cdr m)) (j 1 (1+ j)))
((null m))
(format stream "~a[~a][~a] = " (string-left-trim "$" name) (1- i) (1-
j) )
(c-print (car m) stream)))
'$done)
;; This C scanning function is similar to fortscan. Prepare an expression
;; for printing by converting x^(1/2) to sqrt(x), etc. Since C has no
;; support for complex numbers, contrary to Fortran, ban them.
(defun scanforc (e)
(cond ((atom e) (cond ((eq e '$%i) ;; ban complex numbers
(merror (intl:gettext "Take real and
imaginary parts")))
(t e)))
;; %e^a -> exp(a)
((and (eq (caar e) 'mexpt) (eq (cadr e) '$%e))
(list '(%exp simp) (scanforc (caddr e))))
;; a^1/2 -> sqrt(a) 1//2 is defined as ((rat simp) 1 2)
((and (eq (caar e) 'mexpt) (alike1 (caddr e) 1//2))
(list '(%sqrt simp) (scanforc (cadr e))))
;; a^-1/2 -> 1/sqrt(a)
((and (eq (caar e) 'mexpt) (alike1 (caddr e) -1//2))
(list '(mquotient simp) 1 (list '(%sqrt simp) (scanforc (cadr
e)))))
;; (1/3)*b -> b/3.0 and (-1/3)*b -> -b/3.0
((and (eq (caar e) 'mtimes) (ratnump (cadr e))
(member (cadadr e) '(1 -1) :test #'equal))
(cond ((equal (cadadr e) 1) (scanforc-mtimes e))
(t (list '(mminus simp) (scanforc-mtimes e)))))
;; 1/3 -> 1.0/3.0
((eq (caar e) 'rat)
(list '(mquotient simp) (float (cadr e)) (float (caddr e))))
;; rat(a/b) -> a/b via ratdisrep
((eq (caar e) 'mrat) (scanforc (ratdisrep e)))
;; ban complex numbers
((and (member (caar e) '(mtimes mplus) :test #'eq)
(let ((a (simplify ($bothcoef e '$%i))))
(and (numberp (cadr a))
(numberp (caddr a))
(not (zerop1 (cadr a)))
(merror (intl:gettext "Take real and imaginary
parts"))))))
;; in general do nothing, recurse
(t (cons (car e) (mapcar 'scanforc (cdr e))))))
;; This is used above 1/3*b*c -> b*c/3.0
(defun scanforc-mtimes (e)
(list '(mquotient simp)
(cond ((null (cdddr e)) (scanforc (caddr e)))
(t (cons (car e) (mapcar 'scanforc (cddr e)))))
(float (caddr (cadr e)))))