1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
|
(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 (special $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)))))
|
