;; -*- Mode:Common-Lisp;Package:mma; Base:10 -*-o
;; Mathematica(tm)-like evaluator
;; copyright (c) 1990 Richard J. Fateman; pieces by Tak Yan, Derek Lai
;; hacked more 2/2011 to make it work in GCL, not complete. Need to put
;; two-case symbols within | |.
;; hacked more 3/10/2019 to make it work with SBCL maxima
;;(eval-when (:compile-toplevel) (load "mma"))
(in-package :mma)
;;(provide 'math-eval)(require "ucons1")(require 'math-parser "parser")
;;(require "stack1") (require "disp1")(require "pf")(require "simp1")(require "rat1")(require "match")
(in-package :mma)
(defvar COUNT 1 "line number for top-level. (same as (meval |$Line|))")
(declaim (special env *expand* *numer*)) ;; environment
;; funnyvars is a hash table containing variables which, when set,
;; cause function to be executed
(defvar funnyvars (make-hash-table :test #'eq :size 8))
(defvar emptyht (make-hash-table))
;; We entera variable in the symbol table by making a hash
;; table as its value using chash.
(defvar symtab (make-hash-table :test #'eq :size 150))
(defun today-string()
(multiple-value-bind
(sec min hr date mo year day daylight-p zone)
(get-decoded-time)
(declare (ignore sec daylight-p zone))
(let ((dayx (elt '(|Mon| |Tues| |Wednes| |Thurs| |Fri| |Satur| |Sun|) day))
(mox (elt'(|Jan| |Feb| |Mar| |Apr| |May| |Jun| |Jul| |Aug| |Sep| |Oct| |Nov| |Dec|) (1- mo))))
(format nil "~a:~2d ~aday, ~a ~s, ~s" hr min dayx mox date year))))
(defun tl () ;; top level
(let*
((*package* (find-package :mma)) h hs (hin t)
(timesofar 0)
(timeunit (/ 1.0 internal-time-units-per-second))
(env (make-stack :size 50)) ;environment for matching binding
)
(declare (special env *package*))
(initialize-mma)
(if (= COUNT 1)
(format t
"Mock-Mma ~a ~%" (today-string)))
(loop
(setq timesofar (get-internal-run-time))
(if (eql hin '|Null|) nil
(let nil(format t "~%In[~s] := " COUNT)
(finish-output t))) ;; actually In and Out are variables too.
;; get the input
(setq hin (handler-case (mma::p) (error(x)(format t "~%syntax error ~s" x)(clear-input t) '|Null|)))
(unless (eql hin '|Null|)
;; evaluate it
(setq h (handler-case (meval hin)(error(x) (clrs) ;clear stack frames
(format t "~%evaluation error ~s" x)`(|Hold| , hin))))
(|SetQQ| (ulist '|In| COUNT) hin)
(setq timesofar (- (get-internal-run-time) timesofar))
;; this is not the same as mathematica but I find it more convenient
;; this way. We've also implementing "Timing", if you prefer.
(if (eq '|True| (meval '|$Showtime|))
(format t "~%time = ~3,$ secs.~%" (* timesofar timeunit)))
;; test for exit from top level to lisp REPL
(cond ((or (member h '(EXIT |Exit|) :test 'eql)
(and (listp h)(eq (car h) '|Quit|))) ;;Quit[]
(format t"~%Exited to Lisp~%")
(return t)))
(cond((eq h '|Null|) nil)
;; don't print nil-valued answers
(t
(setq hs (list '|Set| (ulist '|Out| COUNT) h))
(disp (BuildFormat hs))))
;; regardless of printing of answer, set Out[]
(|SetQQ| (ulist '|Out| COUNT) h)
(|Set| '|$Line| (incf COUNT)))) ; loop end
))
;; this definition replaces the program in the parser file
(defun mread1()
(cond((member (pc)'( #\space #\tab #\page) :test #'eql)
(rc)(mread1))
((digit-char-p (mma::pc));; next character is a digit 0-9
(mma::collect-integer
(mma::char-to-int(read-char stream)) 10)) ;radix 10 default
;; for every alphabetic symbol, set up a hash table
(t (chash #-gcl(or(read-preserving-whitespace stream nil 'e-o-l) '|False|)
#+gcl (recase(or(read-preserving-whitespace stream nil 'e-o-l) '|False|)))
;; nil reads as False
)))
#|
We currently make each symbol table entry out of a hash table.
It's plausible to change this to use defstruct ... then
make every declared "symbol-table-entry" a structure
with (at least) the following data
(a) value for the symbol e.g. a=3 in the value cell
(b) value for expressions with the symbol as head. e.g. a[45]=x+1
we might have different "arities" e.g. a[45] has arity 1,
a[3,4] has arity=2, etc. [we don't use this now]
(c) value for the collected attributes of the symbol.
e.g. Attributes[a] ={Orderless, Protected, Listable}
(d) value for each of the attributes to make access fast with using member-test
on collected value [we don't use this now]
(e) value for function definition "built-in" e.g. numeric cosine
(f) value for user-defined function rules e.g. a[x_]:= ...
we could again use some "arity" discrimination if we expect
function definitions of different numbers of arguments. [we don't use this
now]
(g) string for symbol printing (e.g. " + " for Plus). Except that
so far the Lisp symbol-name for Plus is Plus, so we don't need this.
(h) formatting information. right now the List symbol-table entry e.g.
for Plus has the Plus formatter program stored
(i) left/right operator precedence information; display stores this as above
(k) derivative and integral info, as above
(j) messages/ documentation
(k) package? context?
If we were to revise all this we could be more specific
about possible types for the fields,
e.g. for some of these.. (c): list; (d) bit-vector;
(e) lisp-function-value; (f) list? array? (g) string; (h) program,
(k) pattern or program
|#
(defun chash(m)
(let ((*package* (find-package :mma)))
(cond
((not(symbolp m))m)
;; ((null m)nil) do we need to check for nil or t? Maybe not.
(t (cond((gethash m symtab)) ;either it's there or
(t(setf (gethash m symtab)
(make-hash-table :test #'equal :size 4)); we make a hashtable
))))
m))
;; the following stuff is make-shift.
(defun |Head| (h)(typecase h
(cons (car h))
(ratio '|Rational|)
(complex '|Complex|)
(integer '|Integer|)
;((fixnum bignum integer) '|Integer|)
;((double-float single-float) '|Real|)
(float '|Real|)
;;((rat) 'rat) ;; mockmma rat function form, e.g. Rat[x+1].
;; file descriptors? characters? graphics? arrays?
(string '|String| )
(symbol '|Symbol| )
(array '|Array|)
(rat '|Rat|) ;; special rational form
(otherwise (type-of h) )
))
;; Assignment statements treat the lhs with partial evaluation.
;; For a non-atomic Head, evaluate all the arguments, but not
;; the head. Presumably this should check attributes of the Head
;; like |HoldAll|, |Hold|First, etc. We don't do that yet,
;; 12/2010
;; but probably must do |Hold|Rest at least for RuleDelayed so Rubi can work.
;; and Release, Release|Hold|.
;; To the extent possible, I have avoiding thinking about Mma's bogus version of
;; Quote and Eval for as long as I can. for decades?
;; we evaluate the lhs partially and then the rhs.
;; we'd like to have a Quote operator, but the repeated evaluation rule
;; makes it almost impossible to work unless we check for it specially..
;; alternatively, we can set *evalonce* to t, and (vastly)
;; change the semantics. Sometimes this vast change is no change at all....
;; We evaluate args, depending on the hold-specs of the head
;; The Mathematica semantics avoids repeated evaluation by putting a time stamp
;; on each object, as well as dependencies. If the object is newer than any
;; thing it depends on, then it is already "evaluated". Otherwise it
;; has to be re-evaluated and re-time-stamped. We don't do this right now,
;; but we might have to. Mathematica is a little vague and dependencies are
;; not listed explicitly (we think), but approached approximately, e.g.
;; expression x+y+z depends on stuff on some pages (which include the
;; residences of x,y,z. But they might all be on one page. And there
;; may be other things on those pages that change, so x+y+z might be
;; unnecessarily re-evaluated. In our scheme everything is re-evaluated,
;; almost inevitably way too much.
(defun mevalargs( ha l) ;; attributes of head
;;(format t "~%mevalargs env=~s" (describe env))
(unsequence
(cond ((member '|HoldAll| ha :test 'eq)l)
((member '|HoldFirst| ha :test 'eq)
(ucons (car l)(mapcar #'meval (cdr l))))
((member '|HoldRest| ha :test 'eq) (ucons (meval (car l))(cdr l)))
(t ;;(format t "~%in mevalargs l=~s, res=~s" l (umapcar #'meval l))
(umapcar #'meval l)))))
(defun unsequence(k) ;; change ( a b (|Sequence| c d) e) to (a b c d e)
(cond ((null k) nil)
((and (consp (car k))(eq (caar k) '|Sequence|))
(append (cdar k)(unsequence (cdr k))))
(t (setf (cdr k)(unsequence (cdr k))) k)))
;; note that the name of this function conflicts with that
;; of the lisp function set, unless
;; (a) capitalization is observed OR
;; (b) the package system is protecting it..
#+ignore
(defun |Set| (lhs rhs &aux h fun);; lhs=rhs
;; the value associated with the lhs will be stored
;; in the symbol table symtab, with the key h,
;; which is either the head of the lhs,
;; or the lhs itself. That is f=45 is stored in the
;; hash table for f, under the indicator "f"
;; and f[0]=7 is stored in the hash table for f, under
;; the indicator (0).
(cond ((symbolp lhs)(setq h lhs)) ;; simple case, x=stuff
(t (setq h (car lhs)) ;; not so simple, x[i,j]= stuff, evaluate i,j
(setq lhs (mevalargs (|Attributes|(|Head| h))(cdr lhs)))))
;;(format t "Set ~s to ~s, h=~s~%" lhs rhs h)
;; this stores the material in the hash table.
;; QUESTION: M'ma doesn't do this, but we could, by storing
;; stuff on a local environment... f[x_,y_]:=Block[{},x=y+1];
;; what if (gethash h symtab) is a matrix, and this is a valid matrix
;; setting? Then we should try to store the value in the array.
;; This is insufficient error checking but...
(cond
;;else THIS happens, for a global variable. We need to fix for local var
;; on stack.
(t ;;(print 'xxxxx)
;;(chash h)
(setf (gethash lhs (gethash h symtab)) rhs)
;; (setq rhs (meval rhs)) ;; hold it, rhs is already evaluated
))
;; Next, check for special variables which, when set, cause other
;; things to happen. E.g. Precision= number means, in the
;; bigfloat system, (bigfloat-init number) gets executed.
(if (setq h (gethash h funnyvars ))(funcall h rhs))
rhs)
(defun |Set| (lhs rhs) ;; lhs=rhs
;; the value associated with the lhs will be stored
;; in the symbol table symtab, with the key h,
;; which is either the head of the lhs,
;; or the lhs itself. That is f=45 is stored in the
;; hash table for f, under the indicator "f"
;; and f[0]=7 is stored in the hash table for f, under
;; the indicator (0).
(declare (special symtab funnyvars))
(let ((h nil)(fun nil))
(cond ((symbolp lhs) ;; common case, "x=stuff"
;; 3 possibilities, x is bound on stack or or a funny-var or is global
;; (format t "~%setting symbol ~s" lhs)
(multiple-value-bind (val found)(sfind env lhs)
;; (format t "~%mvb result for symbol ~s is ~s ~s " lhs val found)
(declare (ignore val)) ; previous value irrelevant
(cond (found ;found on stack.
(schange env lhs rhs)) ;;change it there.
((setf fun (gethash lhs funnyvars)) ;one of the vars needing special work
;; (format t "~%~s is in funnyvars"~s)
(funcall fun rhs))
(t ;;global case. Not found on stack
;; (format t "~%symtab is ~s" symtab)
(setf h (gethash lhs symtab)) ;odd. should be there from parser.
;; (format t "~%hashtable for ~s is ~s" lhs h)
(cond((null h)(chash lhs) (setf h (gethash lhs symtab))))
;; set h to the hash table for "x"
(setf (gethash lhs h) rhs))))) ; put value on key "x" in "x" hashtable
;; not not so simple an assignment, e.g., x[i,j]= stuff,
((consp lhs)
(cond ((eql (|Head| lhs) '|Part|)
;; case of m[[3]]=45
;; maybe also mm[[3,4]]=45.
;;IMPORTANT: deviation from Mma spec. We return WHOLE expression
(setf rhs (uniq(setpart (append (meval (cadr lhs))nil) ;fresh copy
(cddr lhs)
rhs))) )
(t ;; it is an assignment like f[1,2]=10
(setq h (gethash (car lhs) symtab)) ;symtab for x
(cond((null h)(chash (car lhs)) (setf h (gethash (car lhs) symtab))))
;; (format t "~% Set has non-atom lhs ~s" lhs)
;; (format t "~%hashtable for ~s is ~s" (car lhs) h)
(setq lhs (mevalargs (|Attributes|(car lhs))(cdr lhs))) ; evaluate i,j
(setf (gethash lhs h) rhs)) )))
rhs)) ;return value
;; the following program has an insideous feature.
;; rr=f[a,b,c]
;; rr[[2]]=2
;; now rr is f[a,2,c]. That's good
;; but g[a,b,c]
;; comes out as g[a,2,c] because the single copy in memory of [a,b,c] has been
;; destructively altered to a,2,c.
;; How to fix this???
;; we changed the return value of rr[[2]] so that it is the new value of rr.
;; this version changes the return value of part[m,1]=45 to return a new value for all of m,
;; not just 45.
(defun setpart(m locs val)
(cond ((null (cdr locs)) ;; just one index
(setf (elt m (car locs)) val)
m)
(t(setf (elt m (car locs))(append (elt m (car locs)) nil)) ;make fresh copy
(setpart (elt m (car locs))
(cdr locs)
val))))
;;(format t "Set ~s to ~s, h=~s~%" lhs rhs h)
;; this stores the material in the hash table.
;; QUESTION: M'ma doesn't do this, but we could, by storing
;; stuff on a local environment... f[x_,y_]:=Block[{},x=y+1];
;; what if (gethash h symtab) is a matrix, and this is a valid matrix
;; setting? Then we should try to store the value in the array.
;; Next, check for special variables which, when set, cause other
;; things to happen. E.g. Precision= number means, in the
;; bigfloat system, (bigfloat-init number) gets executed.
;; there is another file (nmatrix) that defines a matrix type..
;; this should be tied in to the matrix stuff from franz, perhaps.
;; also, we have to decide which bigfloat to use... mpfun or rjf's
;; old bfstuff.
;; or MPFR. 12/2010
(defun matrix-p(x) (declare (ignore x))nil) ;;; for now, this will have to do.
(defun |SetQQ|(lhs rhs &aux h);; lhs=rhs, but don't mevaluate either.
(setq h(cond ((atom lhs) lhs)
(t (prog1 (car lhs) (setq lhs (cdr lhs))))))
(setf (gethash lhs (gethash h symtab))rhs)
(if (setq h (gethash h funnyvars ))(funcall h rhs))
rhs)
(defun |Clear|(&rest xl)
(map nil #'clear1 xl)
'|Null|)
(defun clear1(h)
;; RuleHT --clear out rule optimization hash table that seems to have bugs.
;; labels -- resets counter
(declare (special COUNT optimruleht ))
(cond ((eql h '|Labels|) (setf COUNT 0) (|Set|'|$Line| 0))
((eql h '|RuleHT|)
(setf optimruleht (make-hash-table :test 'eq)))
(t (let ((mm (gethash h symtab)))
(if mm (remhash h symtab) nil))))
nil)
;; this hash table has the optimized version of the pretty rules in it
(defparameter optimruleht (make-hash-table :test 'eq))
(defun |SetDelayed|(lhs rhs) ;; this is the function definition f[x_] := ...
(let* ((visible (ulist '|SetDelayed| lhs rhs))
(hidden (bindfix(fixopts visible))))
(setdelayed1 hidden)
(setf (gethash visible optimruleht) hidden) ;; store the optim rule/def here
'|Null|))
(defun setdelayed1(thedef)
(let* ((lhs (ucons (caadr thedef)
(umapcar #'meval (cdadr thedef)))) ;; see note re eval below
(hh (|Head| lhs)) ;; lets hope the Head is a symbol
(spot (gethash hh symtab)))
; (format t "~%lhs=~s" lhs)
;; should we evaluate the formal arguments? we need to do something, e.g.
;; f[Sqrt[x_]] should be changed to f[x_^(1/2)]
;; this could be merely simplification, not evaluation .. ugh.
;; (setf lhs (ucons (car lhs)(mevalargs nil (cdr lhs))))
;;(args (mevalargs (|Attributes|(|Head| hh)) (cdr lhs-in)))
;(args (cdr lhs-in)) ;; these should not be evaluated; they are like formal params, no?
;; push all but the 'SetDelayed on the list of definitional forms.
(if (null spot)(setf spot (setf (gethash hh symtab) (make-hash-table :test #'equal :size 4)))) ; not needed except bug?
(cond ((null (gethash '|SetDelayed| spot))
(setf (gethash '|SetDelayed| spot) (ulist (ucons lhs (cddr thedef)))))
(t (push (ucons lhs (cddr thedef)) (gethash '|SetDelayed| spot))))))
;; this assumes the value of a mathematica symbol is its lisp value
;; if it is simply a constant or non-hash-table. That means that
;; a lisp dynamically bound variable could be used to block access
;; to the globally bound variable of the same name. Better not
;; use local-variable names like Sin or Cos unless you mean to
;; inhibit access to the global names.
(defun meval-atom(h)
(declare (special env *numer*))
(if (member '|Constant| (|Attributes| h):test 'eq)
(if *numer* (numer-meval h)h) ;; works for 3, 1/2
(multiple-value-bind
(val found)
(sfind env h)
(if found val ;; return val
;; if we find it here on the env stack.
;; otherwise
;; val is a symbol, same as h see if it has a value
(let ((r (gethash val symtab)))
(if (hash-table-p r)
(gethash val r val) ;
val))))))
(defun numer-meval(h) (case h
(|Pi| pi) ;; this is just good for double float
(|E| #.(exp 1.0d0))
(otherwise (if(symbolp h) (or (get h 'numerval) h) h))))
;; uh, justa hack
;; works for 3, 1/2
;; look up the value of e on e's hashtable, key e
;; look up the value of e[10] on e's hashtable, key (10)
(defun msymbol-value (h)
(cond
((atom h) (meval-atom h))
(t (let ((tabentry (gethash (|Head| h) symtab)))
(if tabentry (gethash (cdr h) tabentry h)
h)))))
;; this guy looks on the stack; returns if there. looks on global symtab; returns if there.
;; no loop. works, evaluates "once" if that's what you want.
(defun msymbol-function(h)
;; hm, how much work should we do here?
;; it could be that h is just h, and has no binding of any sort.
;; it returns either a non-function or a (global) SetDelayed list or rules
(setf h (cond
;; if we think we should look on the stack for it, then we do this.
;; if the value is yet another name, what then? Do we use that
;; name or continue looking??
((multiple-value-bind
(val found)
(sfind env h)
(if found val nil))) ;; if the function is found on stack, it is val
;; we go on to the next clause if the function is not found
;; on the stack.
((multiple-value-bind
(val found)
(gethash h symtab)
(if found (return-from msymbol-function (gethash '|SetDelayed| val h))
;; if not found return h
h)))))
h)
;; this guy may be useful, if you want to keep looking for the value...
;;for evaluating until no change..
#+ignore
(defun msymbol-function(h)
;; hm, how much work should we do here?
;; it could be that h is just h, and has no binding of any sort.
(let ((saved h))
(loop
(setf h (cond
;; if we think we should look on the stack for it, then we do this.
;; if the value is yet another name, what then? Do we use that
;; name or continue looking??
((multiple-value-bind
(val found)
(sfind env h)
(if found val nil))) ;; if the function is found on stack, it is val
;; we go on to the next clause if the function is not found
;; on the stack.
((multiple-value-bind
(val found)
(gethash h symtab)
(if found (return-from msymbol-function (gethash '|SetDelayed| val h))
;; if not found return h
h)))))
(if (equal h saved) (return h) (setf saved h)))))
;;this version always goes for global function definition ONLY. not compatible with mma ...
#+ignore
(defun msymbol-function(h)
;; hm, how much work should we do here?
(multiple-value-bind
(val found)
(gethash h symtab)
(if found (return-from msymbol-function (gethash '|SetDelayed| val h))
;; if not found return h
h)))
;; is this going to have the right scope?
;;
;;----end of makeshift definitions
(defun mapply (hraw args expr env);; note that args= (cdr expr)
(let* ((fun nil)
(h (meval-to-function hraw))
;;(msymq (gethash h symtab))
)
;; there are 2 kinds of applications here.
;; (Function ...) which is just held.
;; (w ...) where w is (Function ...)
;; (format t "~% h=~s hraw=~s args=~s" h hraw args)
(cond ((eql h '|Function|)
(return-from mapply (ucons '|Function| (funfix (cdr expr)))))
((and (consp h)(eql (car h) '|Function|))
;; this next line should grab the attributes of this Function, if any
;; (format t "~%before args to function ~s are ~s, with env=~s" h args env)
;; (setf args (mevalargs (cdddr h) args)) will be done in mappfun
;; (format t "~%after args to function ~s are ~s" h args)
(return-from mapply (mappfun h args expr env))))
;; (format t "~%before args.. ~s are ~s, with env=~s" h args env)
;;get info on evaluating arguments and do it to args
(setf args (mevalargs (|Attributes| h) args))
(cond
;; I don't believe the comment below...
((constantp h)
;; (format t "~%h=~s" h)
; h
(setf expr (cons h args))
) ;; allows any lisp function, almost, to slip through
;; check for constant values pre-stored
((not (symbolp h)) (setq expr(ucons h args)))
;; oo... this makes lisp function go first. not right??
#+ignore
((and ;(not msymq)
(symbolp h)(fboundp h))
;(format t "~% applying a lisp function ~s to args ~s " h args )
(setq expr (apply h args)))
;;((not msymq)(setq expr (ucons h args)))
;; maybe put in a check for array here? Not now.
;; next check for user-defined function
((consp (setq fun (msymbol-function h)))
;;(setq args (mevalargs h args)) ;; not always ..
;; (format t "~%applying a user function ~s to args ~s" fun args)
(setq expr(rulesetapply h fun args)))
;; next check for built-in LISP function
;; (clearly not something that Mathematica does)
((and ;(not msymq)
(symbolp h)(fboundp h))
;(format t "~% applying a lisp function ~s to args ~s " h args )
(setq expr (apply h args )))
(t ;;(format t "~% eval fell through ~s ~s" h args)
(setq expr(ucons fun args))))
;;(format t "~% mapply result: expr=~s" expr)
expr))
;;(declaim (special phead)) ;; not needed?
;;#+ignore
(defun rulesetapply(phead rules args)
;; get attributes of phead and manipulate args
(setq args (mevalargs (|Attributes|(|Head| phead)) args))
(let*
((origfn phead)
(expr (ucons phead args)))
;; (if isflat nil (setq phead '|Sequence|))
(do ((rr rules (cdr rr)))
((null rr) ;; no more rules to try -- return original
expr)
(let* ((thisrule (car rr))
(condit #'truth)
(lhs (car thisrule))
(rhs (cadr thisrule))
(testr nil))
;; Note: if the rule was
;; f[a_,b_]:= g[a,b] /; a>b, the parsed result is
;; (SetDelayed (f ..) (Condition (g a b) (Greater a b)))
;; see if there is a Condition on the rhs of the rule
;; e.g. (Condition (foo a b) (Greater a b))
;; deal with the possibility of a Condition coming in.
(cond ((and (consp rhs)(eq (car rhs) '|Condition|))
(setf testr (caddr rhs))
;; (format t "~%try to match ~s with condition ~s" lhs testr)
(setf condit #'(lambda() ;(format t "~% evaluating condit ~s" testr)
(meval testr))) ;eg testr = (Comparison a Greater b)
(setf rhs (cadr rhs))))
;rhs = (foo a b)
;;(format t "~%lhs= ~s ~%rhs= ~s expr=~s ~%condition =~s" lhs rhs expr condit)
(if (not (eql (|Head| expr) origfn))(return-from rulesetapply expr)) ;no more rules here are relevant
;; test for matching
;; REDONE 2/2011 RJF,
;;(spushframe env phead)
(cond ((match lhs expr condit) ;; matching works
(setf expr(meval rhs)) ;; something like this..
(spopframe env) ;; pop off the match frame
(return-from rulesetapply expr)) ;; end of match, don't try more rules in do loop
(t (spopframe env)))
))))
;;(defun falsenull(h)(or (null h)(eq h '|False|))) ;; test for mockmma False or maybe nil
;;(defun notfalse(h) (not (falsenull h))) ;; anything not False or nil
;; Major evaluation function for an expression
;; see Mathematica book p 568
;;(defun meval-to-function(x) x) ;; don't evaluate function name ?
(defun meval-to-function(x) (meval x))
(defvar *evalonce* t) ;; should be t to make quote (etc etc) work
(defun meval (e)
(let ((saved e)(hd nil) (ha nil))
(if(atom e)(return-from meval (meval-atom e)))
(cond ((eql (car e) 'lispapply)
#+ignore(format t "~% lispapply found, fun= ~s, args=~s" (cadr e)
(let ((m (meval (caddr e))))
(if (consp m)(cdr m)(list m))))
(return-from meval (apply (cadr e) (let ((m (meval (caddr e))))
(if (consp m)(cdr m)(list m)))))))
(if (eq (setf ha (msymbol-value e))e) ;didn't find a value
nil (return-from meval ha)) ; DID find a value.
;; check off other constant expressions that don't evaluate.
;; perhaps Messages?
;;((patternobjectp e) e) .. What about Patterns?
;; (mapply (car foo)(cdr foo) foo env) ==> foo with no conses...
(setf e (cons (meval-to-function (|Head| e))(cdr e)))
;; (return-from meval (mapply (car e) (cdr e) e env)) ))
(setf e (mapply (car e) (cdr e) e env))
;; note the 3rd arg to mapply, just in case you want to
;; return the expression as the result without any change.
;; next step --
;;
;; do we keep evaluating until there is no change???
(setf hd (|Head| e))
(setf ha (|Attributes| hd))
;; (format t" ~%ha=~s hd=~s e=~s" ha hd e)
(cond (*evalonce* e)
((eql hd '|Function|) (funfix (cdr e)) e) ;; compute new body
((and (member '|Listable| ha) (some #'|ListQ| (cdr e)))
;; uh, (h a,b,..(List c d) ..) becomes
;; (List (h a b ..c..) (h a b .. d..) ).
(listify e))
((equal e saved) e)
((eql hd '|Hold|) e)
(t (meval e)))
))
(defun listify(e)
(let ((listlength nil)(tt 0))
(loop for i in (cdr e) do
;; scan for length
(cond ((|ListQ| i)
(setf tt (length (cdr i)))
(if listlength (cond ((= tt listlength) nil )
(t (format t "~%objects of unequal length cannot be combined ~s" e)
(signal 'error)))
;; no listlength yet. set it
(setf listlength tt)))))
;; now make a list of length listlength
(ucons '|List| (loop for i from 1 to listlength collect
(ucons (car e) (loop for j in (cdr e) collect (if (|ListQ| j)(elt j i) j)))))) )
;; Each global object X is associated with a hash table
;; and we can, for each,
;; to get the value, do (gethash X X), (gethash 'rules X) etc.
;; Local bindings hide everything.
;;Do we want to do infinite evaluation?
;;How do we keep single copies of items in memory?
;;set up initial symbol-table entries for built-in stuff
;; should also set attributes
(defvar *initialized-mma* nil)
(defun initialize-mma ()
(unless *initialized-mma*
#-:allegro (setupcases)
(mapc #'chash built-in-syms)
(globinit '|$Line| 1)
(globinit '|False| nil) ;; maybe, maybe not
(globinit '|$Showtime| nil)
(globinit '|I| #c(0 1)) ;; imaginary unit I^2 = -1.
(setattribute '|Plus| '|Flat|)
(setattribute '|Plus| '|Orderless|)
(setattribute '|Plus| '|Default| 0)
(setattribute '|Plus| '|Listable|)
(setattribute '|Times| '|Flat|)
(setattribute '|Times| '|Orderless|)
(setattribute '|Times| '|Listable|)
(setattribute '|Times| '|Default| 1)
(setattribute '|Power| '|Default| 1)
(setattribute '|And| '|HoldAll|) ; short-circuiting
(setattribute '|Or| '|HoldAll|)
(setattribute '|If| '|HoldRest|)
(setattribute '|Condition| '|HoldRest|)
(setattribute '|Set| '|HoldFirst|)
(setattribute '|Set| '|HoldSequence|) ;; whatever this does; nothing at the moment
(setattribute '|SetDelayed| '|HoldAll|)
(setattribute '|UpSet| '|HoldFirst|)
(setattribute '|UpSetDelayed| '|HoldAll|)
(setattribute '|TagSet| '|HoldFirst|)
(setattribute '|TagSetDelayed| '|HoldAll|)
(setattribute '|Pattern| '|HoldFirst|)
(setattribute '|ReplaceRepeated| '|HoldRest|)
(setattribute '|RuleDelayed| '|HoldRest|)
(setattribute '|Clear| '|HoldAll|)
(setattribute '|Do| '|HoldAll|)
(setattribute '|Table| '|HoldAll|)
(setattribute '|Every| '|HoldAll|)
(setattribute '|Some| '|HoldAll|)
(setattribute '|Function| '|HoldAll|)
(setattribute '|Timing| '|HoldAll|)
(setattribute '|AddTo| '|HoldFirst|)
(setattribute '|Increment| '|HoldFirst|)
(setattribute '|PreIncrement| '|HoldFirst|)
(setattribute '|Decrement| '|HoldFirst|)
(setattribute '|Pi| '|Constant|)
(setattribute '|E| '|Constant|)
(setattribute '|RuleDelayed| '|HoldRest|)
(setf (gethash '|SetDelayed| (gethash '|Power| symtab)) powerrules)
;; do all the ratdiff stuff
(setf (gethash '|SetDelayed| (gethash '|Int| symtab)) integraterules)
(setf (get '|Exp| 'deriv) (make1dfun '(|Exp| %))) ;; or should it be (|Power| E %)???
(setf (get '|Log| 'deriv) (make1dfun '(|Power| % -1))) ;; powersimp better than |Power|??
(setf (get '|Sin| 'deriv) (make1dfun '(|Cos| %)))
(setf (get '|Cos| 'deriv) (make1dfun '(|Times| -1 (|Sin| %))))
(setf (get '|Tan| 'deriv) (make1dfun '(|Power| (|Sec| %) 2)))
(setf (get '|Sec| 'deriv) (make1dfun '(|Times| (|Sec| %) (|Tan| %))))
(setf (get '|Sqrt| 'deriv) (make1dfun '(|Times| 1/2 (|Power| % -1/2)))) ;; not likely to be there
;;... etc rest of Trig functions
(setf (get '|Sinh| 'deriv) (make1dfun '(|Cosh| %)))
;;... etc rest of Hyperbolic functions
(setf (get '|ArcSin| 'deriv) (make1dfun '(|Power| (|Plus| -1 (|Power| % 2)) -1/2)))
(setf (get '|ArcCos| 'deriv) (make1dfun '(|Times| -1 (|Power| (|Plus| 1 (|Times| -1 (|Power| % 2))))
-1/2)))
(setf (get '|ArcTan| 'deriv) (make1dfun '(|Power| (|Plus| 1 (|Power| % 2)) -1)))
(setf (get '|ArcSec| 'deriv) (make1dfun '(|Times| (|Power| (|Plus| 1 (|Times| -1 (|Power| % -2)))
-1/2)
(|Power| % -2))))
;;... etc rest of ArcTrig functions
(setf (get '|ArcSinh| 'deriv) (make1dfun '(|Power| (|Plus| 1 (|Power| % 2)) -1/2)))
(setf (get '|ArcCosh| 'deriv)
(make1dfun '(|Times| (|Power| (|Times| (|Plus| -1 %)(|Power| (|Plus| 1 %) -1)) -1/2)
(|Power| (|Plus| 1 %) -1))))
(setf (get '|ArcTanh| 'deriv) (make1dfun '(|Power| (|Plus| 1 (|Times| -1 (|Power| % 2))) -1)))
(setf (get '|ArcSech| 'deriv)
(make1dfun '(|Times| -1 (|Power| % -1)
(|Power| (|Times| (|Plus| 1 (|Times| -1 %))
(|Power| (|Plus| 1 %) -1))
-1/2)
(|Power| (|Plus| 1 %) -1))))
;; odds and ends: Erf, ExpIntegralEi, Abs
(setf (get '|Erf| 'deriv) (make1dfun '(|Times| 2
(|Power|
(|Times| (|Exp| (|Power| % 2)) (|Power| |Pi| 1/2)) -1))))
(setf (get '|ExpIntegralEi| 'deriv) (make1dfun '(|Times| (|Exp| %) (|Power| % -1))))
;; line below is perhaps a problem when x=0..
(setf (get '|Abs| 'deriv) (make1dfun '(|Times| (|Abs| %)(|Power| % -1))))
;;;; integration properties
(setf (get '|Log| 'integ) '(|Times| % (|Plus| -1 (|Log| %))))
(setf (get '|Sin| 'integ) '(|Times| -1 (|Cos| %)))
(setf (get '|Cos| 'integ) '(|Sin| %))
(setf (get '|Tan| 'integ) '(|Times| -1 (|Log| (|Cos| %))))
(setf (get '|Sec| 'integ) '(|Times| 2 (|ArcTanh|(|Tan| (|Times| 1/2 %)))))
(setf (get '|Cot| 'integ) '(|Times| -1 (|Log| (|Sin| %))))
;;; etc
(setf (get '|ArcSin| 'integ)
'(|Plus| (|Times| % (|ArcSin| %))
(|Power| (|Plus| 1 (|Times| -1 (|Power| % 2))) 1/2)))
(setf (get '|ArcCos| 'integ)
'(|Plus| (|Times| % (|ArcCos| %))
(|Times| -1 (|Power| (|Plus| 1 (|Times| -1 (|Power| % 2))) 1/2))))
(setf (get '|ArcTan| 'integ)
'(|Plus| (|Times| % (|ArcTan| %))(|Times| -1/2 (|Log| (|Plus| 1 (|Power| % 2))))))
;;; etc
(setf (get '|ArcTanh| 'integ)
'(|Plus| (|Times| % (|ArcTanh| %))(|Times| 1/2 (|Log| (|Plus| -1 (|Power| % 2))))))
(setf (get '|Exp| 'integ) '(|Exp| %))
;;
;;Here's some more odds and ends
(setf (get '|ExpIntegralEi| 'integ) '(|Plus| (|Times| % (|ExpIntegralEi| %))
(|Times| -1 (|Exp| %))))
;;int[ erf(_X) ] := _X * erf(_X) + Pi^(-1/2) * exp(-_X^2)
(setf (get '|Erf| 'integ) '(|Plus| (|Times| % (|Erf| %))
(|Power|
(|Times| (|Exp| (|Power| % 2)) (|Power| |Pi| 1/2)) -1)))
(setf (get '|Abs| 'integ) '(|Times| 1/2 % (|Abs| %)))
;;; above computes derivative wrt 1st and only argument..
;;; how can we store the derivative wrt nth argument?
;; here's a place to start. let q=EllipticE[f,m]
;; D[q,x] is Sqrt[1-m*Sin[f]]^2]*D[f,x]
;; + (1/(2m)*(EllipticE[f,m]-EllipticF[f,m)*D[f,m]
;; in the most common situation m is not m(x) but a constant, so the 2nd term
;; disappears.
(setf (get '|EllipticE| 'deriv)
(list #'(lambda(%1 %2)`(|Power|(|Plus| 1 (|Times| -1 ,%2 (|Power| (|Sin| ,%1) 2))) 1/2))
#'(lambda(%1 %2) `(|Times| 1/2 (|Plus| (|EllipticE| ,%1 ,%2)(|Times| -1 (|EllipticF| ,%1 ,%2)))(|Power| ,%2 -1)))))
(setf (get '|EllipticE| 'arity) 2)
;; Can we use the same deal for functions of one argument??
(setf (get '|Sinc| 'deriv)
(list #'(lambda(%1)`(Times (Plus (Times ,%1 (Cos ,%1)) (Times -1 (Sin ,%1)))(Power ,%1 -2)))))
(setf (get '|Sinc| 'arity) 10) ; any number >1 doesn't matter..
(setattribute '|Module| '|HoldAll|)
;; implementing numeric evaluation, we will try to get the
;; lisp global variable *numer* set so that it is fast to check,
;; rather than going through (meval '|Numer|) all over the place
(setf (gethash '|Numer| funnyvars) #'(lambda(h) ;;(format t "~%Numer set to ~s" h)
(setf *numer* h)))
(chash '|Numeric|)
;; this is the list of Numeric functions in Mathematica 7.0
(map nil #'(lambda(r)
(if (get r symtab) nil (chash r))
(setattribute r '|Numeric|))
'(|Abs| |AiryAi|
|AiryAiPrime| |AiryAiZero| |AiryBi| |AiryBiPrime| |AiryBiZero|
|AppellF1| |ArcCos| |ArcCosh| |ArcCot| |ArcCoth| |ArcCsc| |ArcCsch|
|ArcSec| |ArcSech| |ArcSin| |ArcSinh| |ArcTan| |ArcTanh| |Arg|
|ArithmeticGeometricMean| |BellB| |BesselI| |BesselJ| |BesselJZero|
|BesselK| |BesselY| |BesselYZero| |Beta| |BetaRegularized| |Binomial|
|CatalanNumber| |Ceiling| |ChebyshevT| |ChebyshevU| |Clip|
|Conjugate| |Cos| |Cosh| |CoshIntegral| |CosIntegral| |Cot| |Coth|
|Csc| |Csch| |DedekindEta| | Divide| |EllipticE| |EllipticF|
|EllipticK| |EllipticNomeQ| |EllipticPi| |Erf| |Erfc| |Erfi| |Exp|
|ExpIntegralE| |ExpIntegralEi| |Factorial| |Factorial2| |Fibonacci|
|Floor| |FractionalPart| |FresnelC| |FresnelS| |Gamma|
|GammaRegularized| |GegenbauerC| |HankelH1| |HankelH2|
|HarmonicNumber| |HermiteH| |Hypergeometric0F1|
|Hypergeometric0F1Regularized| |Hypergeometric1F1|
|Hypergeometric1F1Regularized| |Hypergeometric2F1|
|Hypergeometric2F1Regularized| |HypergeometricU| |Im| |IntegerPart|
|InverseBetaRegularized| |InverseEllipticNomeQ| |InverseErf|
|InverseErfc| |InverseGammaRegularized| |InverseJacobiCD|
|InverseJacobiCN| |InverseJacobiCS| |InverseJacobiDC|
|InverseJacobiDN| |InverseJacobiDS| |InverseJacobiNC|
|InverseJacobiND| |InverseJacobiNS| |InverseJacobiSC|
|InverseJacobiSD| |InverseJacobiSN| |JacobiAmplitude| |JacobiCD|
|JacobiCN| |JacobiCS| |JacobiDC| |JacobiDN| |JacobiDS| |JacobiNC|
|JacobiND| |JacobiNS| |JacobiP| |JacobiSC| |JacobiSD| |JacobiSN|
|JacobiZeta| |KelvinBei| |KelvinBer| |KelvinKei| |KelvinKer|
|KleinInvariantJ| |LaguerreL| |LerchPhi| |Log| |LogGamma|
|LogIntegral| |LucasL| |MathieuC| |MathieuCharacteristicA|
|MathieuCharacteristicB| |MathieuCharacteristicExponent|
|MathieuCPrime| |MathieuS| |MathieuSPrime| |Max| |Min| |Minus| |Mod|
|ModularLambda| |Multinomial| |NevilleThetaC| |NevilleThetaD|
|NevilleThetaN| |NevilleThetaS| |ParabolicCylinderD| |Plus|
|Pochhammer| |PolyGamma| |PolyLog| |Power| |Quotient|
|QuotientRemainder| |RamanujanTau| |RamanujanTauL|
|RamanujanTauTheta| |RamanujanTauZ| |Re| |Rescale|
|RiemannSiegelTheta| |RiemannSiegelZ| |Round| |Sec| |Sech|
|SiegelTheta| |Sign| |Sin| |Sinc| |Sinh| |SinhIntegral| |SinIntegral|
|SphericalBesselJ| |SphericalBesselY| |SphericalHankelH1|
|SphericalHankelH2| |SphericalHarmonicY| |SpheroidalEigenvalue|
|SpheroidalJoiningFactor| |SpheroidalPS| |SpheroidalPSPrime|
|SpheroidalQS| |SpheroidalQSPrime| |SpheroidalRadialFactor|
|SpheroidalS1| |SpheroidalS1Prime| |SpheroidalS2| |SpheroidalS2Prime|
|Sqrt| |StruveH| |StruveL| |Subfactorial| |Subtract| |Tan| |Tanh|
|Times| |UnitStep| |WhittakerM| |WhittakerW| |Zeta| |ZetaZero|)) ;;whew
(setf *initialized-mma* t)
)
)
;; All system-level $-stuff can be initialized and stored this way
(defun globinit(name val)
(chash name); just in case it isn't already there
(setf (gethash name (gethash name symtab)) val))
;; simple debugging tool: examine contents of symtable entry
(defun showhash(x)
(maphash #'(lambda(key val) (format t "~%key=~s, val=~s" key val)) x))
; |Attributes| that evaluation routines use. maybe
; - Flat [associative, flatten out nested expressions]
; - Orderless [commutative, put args in order]
; - |Hold|First [don't evaluate first arg yet]
; - |Hold|Rest [only evaluate first arg now]
; - |HoldAll| [don't evaluate any args now]
; - Procedure [procedure to call to do actual evaluation]
; - Default [default value for optional variables]
;;[version 3]?
;; the attribute semantics are probably OK as long as only global properties
;; are being recorded.
;; define ClearAttribute similarly
;; mma does not allow setting of |Attributes| other than
;;Protected, ReadProtected, Locked, Temporary, |Hold|First, |Hold|Rest, |HoldAll|, Flat, Orderless, OneIdentity, |List|able, Constant, Stub, N|Hold|First, N|Hold|Rest, N|HoldAll|, NumericFunction, |Sequence||Hold|, and |HoldAll|Complete. [version 7]
;;bunch of additions, 1/2011 RJF
(defun setattribute(h att &optional (val '|True|)) ;; not in mathematica
(if (typep h 'hash-table) nil (setq h (gethash h symtab)));; now its a ht
(setf (gethash att h) val) ;; in h's hashtable, set the attribute att's value to true
(setf(gethash '|Attributes| h);; also in h's hashtable, set its property of |Attributes| to a list
(adjoin att (gethash '|Attributes| h)))
)
(defun |Default|(h)
(gethash '|Default| (gethash h symtab emptyht) '(|Sequence|)))
(defun |SetAttributes|(hi attlist) ;; this is in Mathematica.
(let ((h (gethash hi symtab))) ;; get the hashtable for the symbol h.
(cond ((atom attlist) ;; just one
(setattribute h attlist))
(t (map nil #'(lambda(r)(setattribute h r))(cdr attlist))
(cons '|List| (gethash '|Attributes| h))
;attlist
))
(|Attributes| hi)))
(defun |Attributes|(h)
(cond ((constantp h) '(|List| |Constant|)) ; 3, 1/2, etc
((setf h (gethash h symtab)) (ucons '|List| (gethash '|Attributes| h)))
(t '(|List|)))) ;no attributes for a non-symbol
;; NO, BAD THINGS HAPPEN (setattribute '|Attributes| '|HoldAll|)
;;(chash 'ltrace) ;; lisp trace
;;(chash 'luntrace);; lisp untrace
;;(setattribute 'ltrace '|HoldAll|)
;;(setattribute 'luntrace '|HoldAll|)
;;(defun ltrace(&rest funs) (eval `(trace ,@funs)))
;;(defun luntrace(&rest funs)(eval `(untrace ,@funs)))
;; convert all real numbers to exact rational numbers
(defun |Real|(a b) (+ a b))
;; this works only for integer x, x>0
(defun decimalsin(x)
(ceiling (* (integer-length x) 0.30102999566398114d0)))
;; handle %, %%, etc.
(defun |Out| (&rest n)
(gethash (ulist
(cond ((null n) (1- COUNT)) ; just %, previous item
((and (integerp (setf n (car n)))(minusp n)) ;%%, %%% etc
(+ COUNT n))
(t n)) ;; Out[random stuff]
)
(gethash '|Out| symtab)))
(defun |Simp|(x)(simp x)) ;; rational simplification
(defun |Rat|(x)(into-rat x)) ;; leave the answer in rational form.
(defun |UnRat|(x)(outof-rat x)) ;; convert the answer to list form.
;; convert u to a rational with a single polynomial numerator
;; and denominator. That is (x+1)^2 will be multiplied out.
;; result in rational form.
(defun |RatExpand|(u)
(let* ((*expand* t) ;; global flag to rat program
(x (into-rat u)))
(make-rat :numerator (make-fpe (fpe-expand (rat-numerator x))1)
:denominator
(make-fpe (fpe-expand (rat-denominator x)) 1))))
#+ignore (defun UnivariateDistinctDegreeFactorization (u)
(let((x (into-rat u)))
(make-rat :numerator (poly-uddf (rat-numerator x))
:denominator (poly-uddf (rat-denominator x)))))
#+ignore (defun UnivariateSquareFree etc)
#+ignore (defun UnivariateFactorizationMod(u p)
etc)
;; pick out parts of an expression. x[[y]] parses to Part[x,y].
;; (a+r+b^c) [[3,2]] is c.
;;Generalizes somewhat in that (a+b^c)[[2,r]] returns (b^c)[[r]].
;; Does not handle negative part-numbers or lists of parts as
;; done in Mathematica. Also, won't decompose defstruct items
;; unless we do something about it for each structure...
(defun |Part|(u &rest k)(part1 (meval u) k)) ;; meval??
(defun part1 (u kl)(cond ((null kl) u)
((and(integerp (car kl))
(>= (length u)(car kl)))
(part1 (nth (car kl) u)(cdr kl)))
;; leave unevaluated part if can't handle it.
(t (ucons 'Part (ucons u kl)))))
;;If is HoldRest, so (car stuff) is evaluated.
(defun |If| (&rest stuff)
; (format t "~%If stuff =~s, env=~s" stuff env)
(cond ((eql (car stuff) '|True|) ; test is True?? do we need to meval??
;;(format t "~% then ~s mevals to ~s with env ~s" (cadr stuff)(meval (cadr stuff)) env)
(meval (cadr stuff))) ;; evaluate the "then clause"
((null (car stuff)) ;test is False
; evaluate the "else" if present else return Null
(if (caddr stuff)(meval (caddr stuff)) '|Null|))
(t ;; test is neither true nor false,
(ucons '|If| stuff))))
;; basic simplification of |Times|
(defun |Times| (&rest x &aux (nums 1) oths)
(dolist (h x ;; iterate with h running over all args
;;resultform
(cond((= 1 nums)
(if (null oths) 1
(if (cdr oths) ;; more than one item in product 10/13/94
(ucons '|Times| (uniq(nreverse oths)))
(car oths))))
((= nums 0) 0)
((null oths) nums)
(t (ucons '|Times| (ucons nums (uniq(nreverse oths)))))))
;; body
(cond ((numberp h)(setq nums (* nums h))) ;; collect CL numbers
;; if you find a rat, break out !
((typep h 'rat)(return-from |Times|
(reduce-rat #'rat*
(into-rat (car x))
(cdr x))))
(t (push h oths)))))
(defun |Condition|(a test)
(let ((bool (meval test)))
(cond ((or (null bool)(eql bool '|False|)) '|Null|)
((eql bool '|True|) (meval a);a
)
(t (ulist '|Condition| a test)))))
;; f[3, 4] /. (f[x_, y_] /; x > y -> gg)
;; f[5, 4] /. (f[x_, y_] /; x > y -> gg)
;; this definition allows a=c; to take effect and return Null, so no display
(defun |CompoundExpression|(&rest x &aux result)
(do* ((i x (cdr i))
(j (car i)(car i)))
((null i) result) ;; evaluate each element in turn, return last one.
(setf result(meval j)))
#+ignore (catch :ret ;; if a return is executed somewhere inside here,
)
;; no, we want not to exit NOT from the Compound expression, but the
;; construction outside it
)
(defun |Return|(x)
;;(format t "~% throwing :ret with value ~s" x)
;;(spopframe env) ;;hm. careful here.
(throw :ret x))
(defun |Plus| (&rest x &aux (nums 0) oths)
(dolist (h x ;; iterate with h running over all args
;;resultform
(cond((zerop nums)
(cond ((null oths) 0)
((null (cdr oths)) (car oths))
(t (ucons '|Plus| (uniq(nreverse oths))))))
((null oths) nums)
(t (ucons '|Plus| (ucons nums (uniq(nreverse oths)))))))
;; body
(cond ((numberp h)(incf nums h))
;; if a rat form, break out!
((typep h 'rat)(return-from |Plus|
(reduce-rat #'rat+ (into-rat (car x)) (cdr x))))
(t (push h oths)))))
(defun powersimp (b &optional (e 1)) ;; need to handle (|Power| x) --> x
( cond;;((and (integerp b)(integerp e)) (expt b e)) ;covered by next clause
((and (integerp b)(rationalp e)) (nthroot b e))
((and (rationalp b)(rationalp e))
(|Times| (nthroot (numerator b) e)(nthroot (denominator b) (- e))))
;; maybe put in Complex for sqrt negative..
((and *numer* (numberp b)(numberp e)) (expt b e))
((and (eql '|Rat|(|Head| b))
(integerp e))
(into-rat (uniq `(|Power| ,b ,e))))
((eql e 1) b)
((and (integerp e)(eql (|Head| b) '|Power|)) ;;(x^y)^2 -> x^(2*y). (x^2)^(1/2) no change.
(ulist '|Power| (cadr b)(|Times| (caddr b) e)))
(t (ulist '|Power| b e))))
;; this illustrates how to link lisp programs like powersimp into the
;; structure of rules defined by [x_ ...]:= ..
(defparameter powerrules '(((|Power| (|Pattern| z(|BlankSequence|))) (lispapply powersimp z ))))
;; works
;;(eval-when '(load)
;; (setf (gethash '|SetDelayed| (gethash '|Power| symtab)) powerrules))
;;this makes powersimp the built-in default simplification program, but one can add additional
;; rules on top of it, e.g. zsq^u_Rational :=z^(2*u)
;;until we convert everything to this scheme, we still need
(defun |Power| (b &optional (e 1)) (powersimp b e))
#+ignore;; leave Complex of symbols alone?
(defun |Complex|(re im)
;; insufficient for bigfloats or other number types that are not lisp numbers
(cond ((numberp im)(cond ((numberp re)(complex re im)) ; both re and im are numbers
(t (uniq `(|Plus| ,re ,(complex 0 im))))))
(t (uniq `(Plus ,re (|Times| ,im ,(complex 0 1)))))))
(defun |Complex|(re im)
;; insufficient for bigfloats or other number types that are not lisp numbers
(cond ((and (numberp im)(numberp re))(complex re im)) ; both re and im are numbers
(t `(|Complex| ,re ,im))))
(defun |Rational|(n d)
;; insufficient for bigfloats or other number types that are not lisp integers
(cond ((integerp d)(cond ((integerp n)(/ n d)) ; both n and d are integers
(t (uniq `(|Times| , n ,(/ 1 d))))))
(t ;;(uniq `(|Times| ,n (|Power| ,d -1)))
`(|Rational| ,n ,d))))
#+ignore
(defun rationalsimp(n d)
;; insufficient for bigfloats or other number types that are not lisp integers
(cond ((and(integerp d)(integerp n))(/ n d))
(t (uniq `(|Times| ,n (|Power| ,d -1))))));; probably want to just leave (Rational n d)
#+ignore
(defparameter rationalrules '(((|Rational| (|Pattern| %n(|Blank|))
(|Pattern| %d(|Blank|)))
(lispapply rationalsimp %n %d )))) ;;nope
;;(setf (gethash '|Rational| (gethash '|Rational| symtab)) rationalrules)
;; note: Timing is a |HoldAll| function... otherwise the evaluation
;; of the argument would come first, and the timing would
;; be of an already evaluated expression.
(defun |Timing|(x) ; x is an expression, perhaps compound, uneval'd
(let*((timeunit (/ 1.0 internal-time-units-per-second))
(timesofar (get-internal-run-time))
(result (meval x)))
;; (format t "~% result =~s" result)
(uniq `(|List| (|Times| , (*(- (get-internal-run-time) timesofar)
timeunit)
Second) ,result))))
;; if we are stuck with all one case
;; then we are forced to do something like this for
;; every function that is already defined in lisp with the
;; same name as in mathematica (tm)
;; more notes
;; Sin[x_]/; x>3-> S[x]
;; parses to
;; (Rule (Condition (Sin (Pattern x (Blank))) (Comparison x Greater 3)) (S x))
;; we don't have evaluation of Block implemented.
;;
;; probably should add evaluation of functions and slots.
;; eg #+1&[4] should return 5.
#|
(break "t")
;;(ww (Pattern x (Blank))) Condition aha
;;(Comparison x Greater 5)
;; w[x_]:= aha /; x>5 parses into something like ...
;;(( (ww (Pattern x (Blank))) . (Condition aha (Comparison x Greater 5)) ))
;;but we need (Rule (Condition (ww (Pattern x (Blank))) (Comparison x Greater 5)) aha)
;;or something close to that
;;(trial '(Condition (ww (Pattern x (Blank))) (Comparison x Greater 5)) '(w 10))
;;(trial '(Condition (ww (Pattern x (Blank))) (Comparison x Greater 5)) '(w 4))
|#
;;just a hack ..
#|
(defun macint (x y)
(max2mma (maxima::$integrate (mma2max x) (mma2max y))))
|#
(defun |Every| (exp iter)
(every1 exp iter))
(defun |Some| (exp iter)
(some1 exp iter))
(defun every1 (exp iter) ;; kick out if anything is False (=nil)
;;Exp is an expression with a free variable itvar
;;iter looks like {i,low,hi} or (|List| i low hi) in Lisp
(case (length iter)
(1 (error "invalid iterator ~s" iter))
(2 ;; (List count)
(let ((count (second iter)))
(cond((not(integerp count)) (format t "~%expected integer iterator ~s" iter)
(signal 'error))
((< count 0)(format t "~%expected non-negative iterator ~s" iter)(signal 'error))
(t (setf exp (meval exp))
(loop for i from 1 to count do (if (eql exp '|True|) nil
(return nil)))
'|True|))))
(3
;; (List i hi); no low, assumed 1
;; or (list i (List a b c ...))
(let ((tt (third iter))
(itvar (second iter)))
(cond((and(integerp tt ) (>= tt 1))
(every1 exp (uniq `(|List| ,(second iter) 1 ,tt)))) ;; just count from 1.
((and (consp tt)(eql (car tt) '|List|))
(spush env iter nil)
(do ((i (cdr tt) (cdr i))
(res nil
(progn (schange env itvar (car i))
(if (eql '|True| (meval exp) ) nil (let()(spop env)(return nil))))))
((null i) (spop env) '|True|)))))) ;; kept on going until exhausted so must be
((4 5)
(let ((itvar (second iter)) ;; (List i low hi [step])
(hi (meval (fourth iter))) ;hi
(step (or (meval (fifth iter)) 1))) ;if missing, then 1
(spush env itvar 0) ; reserve a space
;; the case of {i, 1, 10} or {i,1,10,2} ;; set step
(do ((i (meval (third iter)) (+ step i))
(res nil
(progn (schange env itvar i)
(if (eql '|True| (meval exp) ) nil (let()(spop env)(return nil))))))
((> i hi) (spop env) '|True|))))
))
(defun some1 (exp iter) ;; kick out if anything is True
;;Exp is an expression with a free variable itvar
;;iter looks like {i,low,hi} or (|List| i low hi) in Lisp
(case (length iter)
(1 (error "invalid iterator ~s" iter))
(2 ;; (List count)
(let ((count (second iter)))
(cond((not(integerp count)) (format t "~%expected integer iterator ~s" iter)
(signal 'error))
((< count 0)(format t "~%expected non-negative iterator ~s" iter)(signal 'error))
(t (setf exp (meval exp))
(loop for i from 1 to count do (if (eql exp '|True|)
(return nil) nil))
'|True|))))
(3
;; (List i hi); no low, assumed 1
;; or (list i (List a b c ...))
(let ((tt (third iter))
(itvar (second iter)))
(cond((and(integerp tt ) (>= tt 1))
(some1 exp (uniq `(|List| ,(second iter) 1 ,tt)))) ;; just count from 1.
((and (consp tt)(eql (car tt) '|List|))
(spush env iter nil)
(do ((i (cdr tt) (cdr i))
(res nil
(progn (schange env itvar (car i))
(if (eql '|True| (meval exp)) (let()(spop env)(return '|True|)) nil))))
((null i) (spop env) nil)))))) ;; kept on going until exhausted so none..
((4 5)
(let ((itvar (second iter)) ;; (List i low hi [step])
(hi (meval (fourth iter))) ;hi
(step (or (meval (fifth iter)) 1))) ;if missing, then 1
(spush env itvar 0) ; reserve a space
;; the case of {i, 1, 10} or {i,1,10,2} ;; set step
(do ((i (meval (third iter)) (+ step i))
(res nil
(progn (schange env itvar i)
(if (eql '|True| (meval exp) ) (let()(spop env)(return nil)) nil))))
((> i hi) (spop env) nil))))
))
#|
;; This is the way to make rules / function definitions in LISP.
;; we can make these rules for, for example, Power. Then override them
;; by user stuff as necessary.
(setf h5rule '((h5 (Pattern z(BlankSequence))) (lispapply h5fun z )))) ;; works
;;(setf h5rule '((h5 (Pattern z(BlankSequence))) (lispapply (lambda(&rest r)(print r)) z )))) ;; works
(defun h5fun(&rest r)(cons 'defaulth5 r))
(setf h5rules (list h5rule))
(setf (gethash '|SetDelayed| (gethash 'h5 symtab)) h5rules)
(chash 'lispapply)
(setattribute 'lispapply '|HoldAll|)
|#
(defun iroot (a n)
;; computes a^(1/n) integer a>0 n>0 see Fitch, SIGSAM Bull Nov 74
;; second value is distance from exact root
(let ((ila (integer-length a)))
(cond ((= a 0)(values 0 0))
((< ila n) (values 1 (1- a)))
(t ;assumes integer a>0 n>=2
(do ((x (expt 2 (1+ (truncate ila n)))
(- x (truncate (+ n1 bk) n)))
(n1 (1- n)) (xn) (bk))
(nil)
(cond ((<= (setq bk (- x (truncate a (setq xn (expt x n1))))) 0)
(return (values x (- a (* x xn)))))))))))
(defun nthroot(a n) ;; eg. (nthroot 8 2/3); (nthroot 9 1/3)
(if (< a 0)(ulist '|Power| a n)
(multiple-value-bind
(r leftover)
(iroot a (denominator n))
(if (= leftover 0)(expt r (numerator n)) (ulist '|Power| a n)))))
(defun sqrtexact(a) ;; a is an integer. return sqrt(a)
(let ((complex nil)(b a))
(cond ((< b 0)(setf complex t b (- a))))
(let ((rt (nthroot b 1/2)))
(cond ((integerp rt) (if complex (setf rt (complex 0 rt))))
(t (setf rt (ulist '|Power| a 1/2))))
rt)))
;; setting LineLength=100 will change COL
(setf (gethash 'LineLength funnyvars)#'(lambda(r)(declare (special COL))
(and (integerp r)(> r 10)(setf COL r))))