[Maxima-commits] CVS: maxima/share/diff_form cartan_new.lisp, NONE, 1.1 coeflist.lisp, NONE, 1.1 curvture2.mac, NONE, 1.1 diff_form.mac, NONE, 1.1 example.txt, NONE, 1.1 f_star_test4.mac, NONE, 1.1 format.lisp, NONE, 1.1 frobenius.mac, NONE, 1.1 helpfunc.mac, NONE, 1.1 hodge_test3.mac, NONE, 1.1 lorentz_example.txt, NONE, 1.1 new_cartan_test4.mac, NONE, 1.1 poisson.mac, NONE, 1.1 readme_diff_form.txt, NONE, 1.1 surface_example.txt, NONE, 1.1

[Robert D. <rob...@us...>]

Update of /cvsroot/maxima/maxima/share/diff_form
In directory sc8-pr-cvs7.sourceforge.net:/tmp/cvs-serv18110

Added Files:
	cartan_new.lisp coeflist.lisp curvture2.mac diff_form.mac 
	example.txt f_star_test4.mac format.lisp frobenius.mac 
	helpfunc.mac hodge_test3.mac lorentz_example.txt 
	new_cartan_test4.mac poisson.mac readme_diff_form.txt 
	surface_example.txt 
Log Message:
Commit files submitted by Gosei Furuya in message dated Oct 29, 2006 3:16 AM
with subject line "Re: [Maxima] differencial form package".
All files are the same as submitted, except that two have been renamed:
rename cartan_init.bat --> diff_form.mac (make initialization file name same
as the name of the package, and make it a .mac file), and 
rename new_cartan_test4.bat --> new_cartan_test4.mac (make it a .mac file).

--- NEW FILE: cartan_new.lisp ---
;;this program is derived from maxima/share/calculus/cartan.lisp
;;I add clifford operator &,&2
;;I add clifford differential operator d_c
  

(DEFPROP $\| %\| VERB) 
(DEFPROP $\| &\| OP) 
(DEFPROP &\| $\| OPR) 
(ADD2LNC (QUOTE &\|) $PROPS) 
(DEFPROP $\| DIMENSION-infix DIMENSION) 
(DEFPROP $\| (32 124 32) DISSYM) 
(DEFPROP $\| 120 LBP) 
(DEFPROP $\| 180 RBP) 
(DEFPROP $\| PARSE-INFIX LED) 
(DEFPROP $\| MSIZE-INFIX GRIND) 
(DEFPROP %\| DIMENSION-infix DIMENSION) 
(DEFPROP %\|  (32 124 32) DISSYM) 
(MDEFPROP $\| ((LAMBDA) ((MLIST) $V $F) ((MPROG SIMP) ((MLIST SIMP)
 $I $J $EXT101 $EXT102 $EXT103 $EXT104) ((MSETQ SIMP) $EXT103 
(($EXPAND SIMP) $F)) ((MSETQ SIMP) $EXT102 ((MTIMES SIMP) 
(($V SIMP ARRAY) 1) (($COEFF SIMP) $EXT103 (($BASIS SIMP ARRAY) 1)))) 
((MDO SIMP) $I 2 NIL NIL $DIM NIL ((MPROGN SIMP) ((MSETQ SIMP) 
$EXT101 (($COEFF SIMP) $EXT103 (($BASIS SIMP ARRAY) $I))) 
((MCOND SIMP) ((MNOTEQUAL SIMP) $EXT101 0) ((MSETQ SIMP) $EXT101
 (($SUBSTITUTE SIMP) (($EXTSUB SIMP ARRAY) $I) $EXT101)) T $FALSE)
 ((MSETQ SIMP) $EXT102 ((MPLUS SIMP) $EXT102 ((MTIMES SIMP) $EXT101
 (($V SIMP ARRAY) $I)))))) ((MRETURN SIMP) (($EXPAND SIMP) $EXT102))))
 MEXPR) 
(ADD2LNC (QUOTE (($\|) $V $F)) $FUNCTIONS) 
(DEFPROP %\| $\| NOUN) 

(DEFPROP $@ %@ VERB) 
(DEFPROP $@ &@ OP) 
(DEFPROP &@ $@ OPR) 
(ADD2LNC (QUOTE &@) $PROPS) 
(DEFPROP $@ DIMENSION-infix DIMENSION) 
(DEFPROP $@ (32 38 32) DISSYM) 
;;(DEFPROP $@ 140 LBP) 
;;(DEFPROP $@ 180 RBP)
(DEFPROP $@ 80 LBP)
(DEFPROP $@ 100 RBP) 
(DEFPROP $@ PARSE-INFIX LED) 
(DEFPROP $@ MSIZE-INFIX GRIND) 
(DEFPROP %@ DIMENSION-infix DIMENSION) 
(DEFPROP %@ (32 38 32) DISSYM)
;;exterior operator @ 
(MDEFPROP $@ 
	  ((LAMBDA) 
	   ((MLIST) $F $G) 
	   ((MPROG SIMP) 
	    ((MLIST SIMP) $I $J $EXT101 $EXT102 $EXT103 $EXT104 $EXT105) 
	    ((MSETQ SIMP) $EXT101 0) 
	    ((MSETQ SIMP) $EXT102 $TRUE) 
	    ((MSETQ SIMP) $EXT103 
	     (($EXPAND SIMP) $F)) 
	    ((MDO SIMP) $I $DIM -1 NIL 0 NIL 
	     ((MPROGN SIMP) 
	      ((MSETQ SIMP) $EXT104 
	       (($EXPAND SIMP) 
		(($BOTHCOEF SIMP) $EXT103 
		 (($BASIS SIMP ARRAY) $I)))) 
	      ((MSETQ SIMP) $EXT105 
	       (($FIRST SIMP) $EXT104)) 
	      ((MCOND SIMP) 
	       ((MNOTEQUAL SIMP) $EXT105 0) 
	       ((MPROGN SIMP) 
		((MSETQ SIMP) $EXT103 (($LAST SIMP) $EXT104))  
		((MSETQ SIMP) $EXT101 
		 ((MPLUS SIMP) $EXT101 
;;(($@ SIMP) $EXT105
;;move $EXT105 later 
		  (($@ SIMP) $EXT105
		   ((MTIMES SIMP)
		    (($BASIS SIMP ARRAY) $I) 
		    (($SUBSTITUTE SIMP) 
		     (($EXTSUBB SIMP ARRAY) $I) $G)) )))  
		((MSETQ SIMP) $EXT102 $FALSE)) T $FALSE))) 
	    ((MCOND SIMP) $EXT102 
	     ((MRETURN SIMP) 
	      (($EXPAND SIMP) ((MTIMES SIMP) $F $G))) 
	     T ((MRETURN SIMP) (($EXPAND SIMP) $EXT101))))) 
	  MEXPR)
 
(ADD2LNC (QUOTE (($@) $F $G)) $FUNCTIONS) 
(DEFPROP %@ $@ NOUN) 
;;start definition with clifford operator
(DEFPROP $& %& VERB) 
(DEFPROP $& &\& OP) 
(DEFPROP &\& $\& OPR) 
(ADD2LNC (QUOTE &\&) $PROPS) 
(DEFPROP $& DIMENSION-infix DIMENSION) 
(DEFPROP $& (32 38 32) DISSYM) 
(DEFPROP $& 140 LBP) 
(DEFPROP $& 180 RBP) 
(DEFPROP $& PARSE-INFIX LED) 
(DEFPROP $& MSIZE-INFIX GRIND) 
(DEFPROP %& DIMENSION-infix DIMENSION) 
(DEFPROP %& (32 38 32) DISSYM) 
;;clifford operator &
(MDEFPROP $& ((LAMBDA) 
	      ((MLIST) $F $G) 
	      ((MPROG SIMP)
	       ((MLIST SIMP) $I $J $EXT101 $EXT102 $EXT103 $EXT104 $EXT105)
	       ((MSETQ SIMP) $EXT101 0) 
	       ((MSETQ SIMP) $EXT102 $TRUE) 
	       ((MSETQ SIMP) $EXT103 
		(($EXPAND SIMP) $F)) 
	       ((MDO SIMP) $I $DIM -1 NIL 1 NIL 
		((MPROGN SIMP) ((MSETQ SIMP) $EXT104 
				(($EXPAND SIMP) 
				 (($BOTHCOEF SIMP) $EXT103 
				  (($BASIS SIMP ARRAY) $I))))
		 ((MSETQ SIMP) $EXT105 (($FIRST SIMP) $EXT104)) 
		 ((MCOND SIMP) ((MNOTEQUAL SIMP) $EXT105 0) 
		  ((MPROGN SIMP) ((MSETQ SIMP) $EXT103 (($LAST SIMP) $EXT104))
		   ((MSETQ SIMP) $EXT101 ((MPLUS SIMP) $EXT101 
;;(($& SIMP) $EXT105				  
					  (($& SIMP) $EXT105
					   ((MTIMES SIMP)  
					    (($BASIS SIMP ARRAY) $I) 
					    (($SUBSTITUTE SIMP) 
					     (($EXTSUBB2 SIMP ARRAY) $I) $G))
					   )))
		   ((MSETQ SIMP) $EXT102 $FALSE)) T $FALSE))) 
	       ((MCOND SIMP) $EXT102 
		((MRETURN SIMP) (($EXPAND SIMP) 
				 ((MTIMES SIMP) $F $G))) T 
				 ((MRETURN SIMP) (($EXPAND SIMP) $EXT101))))) MEXPR) 

(ADD2LNC (QUOTE (($&) $F $G)) $FUNCTIONS) 
(DEFPROP %& $& NOUN)

;;exterior differntial operator
(MDEFPROP $D ((LAMBDA) 
	      ((MLIST) $F) 
	      (($SUM SIMP) 
	       (($@ SIMP) 
		(($BASIS SIMP ARRAY) $I) 
		(($DIFF SIMP) $F 
		 (($COORDS SIMP ARRAY) $I)))
	       $I 1 $DIM)) MEXPR) 

(ADD2LNC (QUOTE (($D2) $F)) $FUNCTIONS)

;;clifford differntial operator
(MDEFPROP $D_C ((LAMBDA) 
	      ((MLIST) $F) 
	      (($SUM SIMP) 
	       (($& SIMP) 
		(($BASIS SIMP ARRAY) $I) 
		(($DIFF SIMP) $F 
		 (($COORDS SIMP ARRAY) $I)))
	       $I 1 $DIM)) MEXPR)

(ADD2LNC (QUOTE (($D_C) $F)) $FUNCTIONS)

;;another clifford operator with different basis.
;;clifford operator & and exterior operator go with.
;;but we ofen need to calc clifford algebra another basis at the same time.
;;before using &2 operator you must "infix("&2")$ "
(DEFPROP $&2 %&2 VERB) 
(DEFPROP $&2 &\&2 OP) 
(DEFPROP &\&2 $\&2 OPR) 
(ADD2LNC (QUOTE &\&2) $PROPS) 
(DEFPROP $&2 DIMENSION-infix DIMENSION) 
(DEFPROP $&2 (32 38 32) DISSYM) 
(DEFPROP $&2 140 LBP) 
(DEFPROP $&2 180 RBP) 
(DEFPROP $&2 PARSE-INFIX LED) 
(DEFPROP $&2 MSIZE-INFIX GRIND) 
(DEFPROP %&2 DIMENSION-infix DIMENSION) 
(DEFPROP %&2 (32 38 32) DISSYM) 
;;clifford operator &
(MDEFPROP $&2 ((LAMBDA) 
	      ((MLIST) $F $G) 
	      ((MPROG SIMP)
	       ((MLIST SIMP) $I $J $EXT101 $EXT102 $EXT103 $EXT104 $EXT105)
	       ((MSETQ SIMP) $EXT101 0) 
	       ((MSETQ SIMP) $EXT102 $TRUE) 
	       ((MSETQ SIMP) $EXT103 
		(($EXPAND SIMP) $F)) 
	       ((MDO SIMP) $I $N_DIM -1 NIL 1 NIL 
		((MPROGN SIMP) ((MSETQ SIMP) $EXT104 
				(($EXPAND SIMP) 
				 (($BOTHCOEF SIMP) $EXT103 
				  (($BASIS2 SIMP ARRAY) $I))))
		 ((MSETQ SIMP) $EXT105 (($FIRST SIMP) $EXT104)) 
		 ((MCOND SIMP) ((MNOTEQUAL SIMP) $EXT105 0) 
		  ((MPROGN SIMP) ((MSETQ SIMP) $EXT103 (($LAST SIMP) $EXT104))
		   ((MSETQ SIMP) $EXT101 ((MPLUS SIMP) $EXT101 
					  (($&2 SIMP) $EXT105 ((MTIMES SIMP) 
							      (($BASIS2 SIMP ARRAY) $I) 
							      (($SUBSTITUTE SIMP) 
							       (($EXTSUBB4 SIMP ARRAY) $I) $G)))))
		   ((MSETQ SIMP) $EXT102 $FALSE)) T $FALSE))) 
	       ((MCOND SIMP) $EXT102 
		((MRETURN SIMP) (($EXPAND SIMP) 
				 ((MTIMES SIMP) $F $G))) T 
				 ((MRETURN SIMP) (($EXPAND SIMP) $EXT101))))) MEXPR) 

(ADD2LNC (QUOTE (($&2) $F $G)) $FUNCTIONS) 
(DEFPROP %&2 $&2 NOUN)
 

 












--- NEW FILE: coeflist.lisp ---
;;; -*- Mode: LISP; Syntax: Common-lisp; Package: Maxima; Base: 10 -*-
;;;>******************************************************************************************
;;;>   Software developed by Bruce R. Miller
;;;>   using Symbolics Common Lisp (system 425.111, ivory revision 4)
;;;>   at NIST - Computing and Applied Mathematics Laboratory
;;;>   a part of the U.S. Government; it is therefore not subject to copyright.
;;;>******************************************************************************************

;(in-package 'climax)

;;; To run in Schelter's Maxima comment the above and uncomment these:
(in-package :maxima)
(defmacro mexp-lookup (item alist) `(assolike ,item ,alist))
(defmacro mlist* (arg1 &rest more-args) `(list* '(mlist simp) ,arg1 ,@more-args))

;;;;******************************************************************************************
;;;; Needed, but unrelated, stuff.  Possibly useful in its own right
;;;;******************************************************************************************
;;; Returns an mlist of all subexpressions of expr which `match' predicate.
;;; Predicate(expr,args,...) returns non-nil if the expression matches.
;;; [eg. a function constructed by $defmatch]
(defun $matching_parts (expr predicate &rest args)
  (let ((matches nil))
    (labels ((srch (expr)
	       (when (specrepp expr)
		 (setq expr (specdisrep expr)))
	       (when (apply #'mfuncall predicate expr args)
		 (pushnew expr matches :test #'alike1))
	       (unless (atom expr)
		 (mapc #'srch (cdr expr)))))
      (srch expr)
      (mlist* matches))))

;;; Return an  mlist of all unique function calls of the function(s) FCNS in EXPR.
;;; (FCNS can also be a single function)
(defun $function_calls (expr &rest functions)
  ;; Coerce fcns to a list of function names.
  (let ((fcns (mapcar #'(lambda (x)(or (get (setq x (getopr x)) 'verb) x)) functions)))
    ($matching_parts expr #'(lambda (p)(and (listp p)(member (caar p) fcns))))))

;;; Return an mlist of all unique arguments used by FCNS in EXPR.
(defun $function_arguments (expr &rest functions)
  (mlist* (remove-duplicates (cdr (map1 #'$args (apply #'$function_calls expr functions)))
			     :test #'alike1)))

;;; $totaldisrep only `disrep's CRE (mrat), but not POIS!
(defun totalspecdisrep (expr)
  (cond ((atom expr) expr)
	((not (or (among 'mrat expr)(among 'mpois expr))) expr)
	((eq (caar expr) 'mrat)(ratdisrep expr))
	((eq (caar expr) 'mpois) ($outofpois expr))
	(t (cons (remove 'ratsimp (car expr))(mapcar 'totalspecdisrep (cdr expr))))))

;;;;******************************************************************************************
;;;; Variable Lists
;;; A variable list consists of a list of variables, simple expressions and specs like
;;;   OPERATOR(fcn) or OPERATOR(fcn,...) represents ALL calls to fcn in the expression.
;;;   MATCH(fcn,arg..) represents subexpressions of expression which pass FCN(subexpr,args..)
;;; Instanciating the variable list involves replacing those special cases with those
;;; subexpressions of the relevant expression which pass the test.
;;;;******************************************************************************************

(defun instanciate-variable-list (vars expr caller &optional max-vars)
  (let ((ivars (mapcan #'(lambda (var)
			   (setq var (totalspecdisrep var))
			   (case (and (listp var)(caar var))
			     ($operator (cdr (apply #'$function_calls expr (cdr var))))
			     ($match    (cdr (apply #'$matching_parts expr (cdr var))))
			     (t (list var))))
		       vars)))
    (when (and max-vars (> (length ivars) max-vars))
      (merror "Too many variables for ~M: ~M" caller (mlist* ivars)))
    ivars))

;;;;******************************************************************************************
;;;; Helpers
;;;;******************************************************************************************

;;; Similar to lisp:reduce with the :key keyword.
;;; Apparently, the Lisp underneath the Sun version doesn't support it. Ugh.
;  (defmacro cl-reduce (function list key) `(lisp:reduce ,function ,list :key ,key))

(defun cl-reduce (function list key)
  (if (null list) nil
      (let ((result (funcall key (car list))))
	(dolist (item (cdr list))
	  (setq result (funcall function result (funcall key item))))
	result)))

(defun map-mlist (list) (mapcar #'(lambda (e)(mlist* e)) list))

;;;******************************************************************************************
;;; Coefficient List = Pseudo-polynomial : as list of ( (coefficient power(s) ...) ...)
;;;     coefficient: the coefficient of the monomial (anything algebraic thing)
;;;     power(s) : the power(s) of the variable(s) in the monomial (any algebraic thing)
;;; Pairs are sorted into increasing order of the power(s).   
;;; 0 is represented by NIL.
;;;******************************************************************************************
;;; NOTE on ordering of terms.  The Macsyma predicate GREAT (& friends lessthan, etc)
;;; define a total ordering, but if non-numeric elements are allowed, the ordering is not
;;; robust under addition, eg L=[1,2,m] is in order, but L+m=[m,m+2,2*m] is not.
;;; We define the ordering of A & B by determining the `sign' of A-B, where the sign is
;;; the sign of the coefficient of the leading (highest degree) term.  We can use SIGNUM1 
;;; for this. 
;;;******************************************************************************************

;;;******************************************************************************************
;;; CLIST Arithmetic

;;; Add two CLISTs
(defun clist-add (l1 l2)
  (do ((result nil))
      ((not (and l1 l2)) (if (or l1 l2)(nconc (nreverse result)(or l1 l2)) (nreverse result)))
    (do ((p1 (cdar l1) (cdr p1))
	 (p2 (cdar l2) (cdr p2)))
	((or (null p1) (not (alike1 (car p1)(car p2))))
	 (if p1
	     (push (if (plusp (signum1 (sub (car p1)(car p2))))(pop l2)(pop l1)) result)
	     (let ((c3 (add (caar l1)(caar l2))))	;If power is same, combine
	       (unless (zerop1 c3)		;And accumulate result, unless zero.
		 (push (cons c3 (cdar l1)) result))
	       (pop l1)(pop l2)))))))

;;; Multiply two CLISTs
;;; Optional ORDER is for use by series arithmetic (single variable): truncates powers>order
(defun clist-mul (l1 l2 &optional order)
  (when (and l1 l2)
    (when (> (length l1)(length l2))		; make l1 be shortest
      (psetq l1 l2 l2 l1))
    (let ((rl2 (reverse l2)))
      (flet ((mul1 (pair1)
	       (let ((c1 (car pair1)) (p1 (cdr pair1)) result)
		 (dolist (i2 rl2)
		   (let ((p (mapcar #'add p1 (cdr i2))))
		     (unless (and order (great (car p) order))
		       (push (cons (mul c1 (car i2)) p) result))))
		 result)))
 	(cl-reduce #'clist-add l1 #'mul1)))))

;;; Take the Nth power of a CLIST, using "binary expansion of the exponent" method.
;;; Built-in code to handle P^2, instead of P*P.
(defun clist-pow (l n)				; Assumes n>0
  (cond ((null l) nil)				; l=0 -> 0 (nil)
	((null (cdr l))				; single term, trivial
	 `((,(power (caar l) n) ,@(mapcar #'(lambda (p)($expand (mul p n)))(cdar l)))))
	(t (let ((l^i l) (l^n (if (logtest n 1) l)))
	     (do ((bits (lsh n -1)(lsh bits -1)))
		 ((zerop bits) l^n)
	       (do ((sq nil)			; Square l^i
		    (ll (reverse l^i) (cdr ll)))
		   ((null ll) (setq l^i sq))
		 (let* ((c1 (caar ll)) (2c1 (mul 2 c1))(p1 (cdar ll))
			(psq (list (cons (power c1 2)(mapcar #'add p1 p1)))))
		   (dolist (lll (cdr ll))
		     (push (cons (mul 2c1 (car lll))(mapcar #'add p1 (cdr lll))) psq))
		   (setq sq (if sq (clist-add sq psq) psq))))
	       (if (logtest bits 1) (setq l^n (if l^n (clist-mul l^n l^i) l^i))))))))

;;; An MBAG includes lists, arrays and equations.
;;; Given the list of {list|array|equation} elements which have been converted to CLIST's, 
;;; this function combines them into a single clist whose coefficients 
;;; are {list|array|equation}s

(defun clist-mbag (op clists)
  (let ((z (if (eq op '$MATRIX)			; the `zero' of a matrix is an mlist of 0's!!!
	       (mlist* (make-list (length (cdaaar clists)) :initial-element 0))
	       0)))
    (flet ((keylessp (l1 l2)			; does key l1 precede l2?
	     (do ((l1 l1 (cdr l1))
		  (l2 l2 (cdr l2)))
		 ((or (null l1)(not (alike1 (car l1)(car l2))))
		  (and l1 (minusp (signum1 (sub (car l1)(car l2)))))))))
      (mapcar #'(lambda (p)
		  `(((,op)
		     ,@(mapcar #'(lambda (l)(or (car (rassoc p l :test #'alike)) z)) clists))
		    ,@p))
	      (sort (cl-reduce #'union1 clists #'(lambda (e)(mapcar #'cdr e))) #'keylessp)))))

;;;;******************************************************************************************
;;;; Transform an expression into its polynomial coefficient list form.

(defun $coeffs (expr &rest vars)
  (setq expr (totalspecdisrep expr))
  (let* ((vs (instanciate-variable-list vars expr '$coeffs))
	 (zeros (make-list (length vs) :initial-element 0))
	 (cache nil))
    (dolist (v vs)				; preload the cache w/ encoded variables
      (let ((u (copy-list zeros)))
	(setf (nth (position v vs) u) 1)
	(push (cons v (list (cons 1 u))) cache)))
    (labels ((gcf (expr)			; Get coefficients.
	       (or (mexp-lookup expr cache)	; reuse cached value
		   (cdar (push (cons expr (gcf1 expr)) cache))))	; or compute & store
	     (gcf1 (expr)
	       (let ((op (and (listp expr)(caar expr))) x y)
		 (cond ((MBAGp expr)    (clist-mbag op (mapcar #'gcf (cdr expr))))
		       ((or (null op)(not (dependsall expr vs))) `((,expr . ,zeros)))
		       ((eq op 'MPLUS)  (cl-reduce #'clist-add (cdr expr) #'gcf))
		       ((eq op 'MTIMES) (cl-reduce #'clist-mul (cdr expr) #'gcf))
		       ((and (eq op 'MEXPT)	; Check that we can actually compute X^Y:
			     (setq x (gcf (second expr)) y (third expr))
			     (or (and (integerp y)(plusp y))	; Either integer y > 0
				 (and (null (cdr x))	; or x is a single monomial
				      (not (dependsall y vs))	; w/ y must be free of vars
				      (or (eql $RADEXPAND '$ALL) ; & dont care about cuts
					  (integerp y)	; or y is an integer
					  (every #'(lambda (p)(or (zerop1 p)(onep p)))
						 (cdar x))))))	; or x is linear in vars.
			(clist-pow x y))	; OK, so we CAN compute x^y (whew).
		       (t `((,expr . ,zeros)))))))
      (mlist* (mlist* '$%POLY vs)(map-mlist (gcf expr))))))

; Inverse of above: make an expression out of clist.
;;; Actually works for SERIES & Taylor too.
(defun unclist (clist vars)
  (addn (mapcar #'(lambda (e)(mul (cadr e)(muln (mapcar #'power vars (cddr e)) t))) clist) t))

;;;********************************************************************************
;;; TRIG SERIES
;;; Given an expression and a list of variables, v_i, construct the list of sine & cosine
;;; coefficients & multiples of the variables in the expression:
;;;  [[%trig, v_1, ...] sine_list, cosine_list]
;;;     sine_list: [[c,m_1,...],[c',m_1',...]....]
;;;     cosine_list: " "
;;; This version carries out `trig arithmetic' on coefficient lists (does NOT use the 
;;; enhanced poisson package (pois2m)
;;;********************************************************************************

;;;;******************************************************************************************
;;;; TLIST Arithmetic.

;;; Is L1 < (0 0 0 ...)? ie. is first non-zero element negative?
(defun list-negp (l1)
  (dolist (i1 l1)
    (unless (zerop1 i1)
      (return-from list-negp (eq -1 (signum1 i1))))))

(defun tlist-add (l1 l2) (mapcar #'clist-add l1 l2))

; multiply a cos or sin list by another.  c1p is T if L1 is cosine list, c2p ditto.
(defun tlist-mul1 (l1 l2 c1p c2p)
  (when (> (length l2)(length l1))		; Swap so that L2 is the shortest.
    (psetq l1 l2 l2 l1)
    (psetq c1p c2p c2p c1p))
  (when l2
    (let ((s1 (if (and c1p (not c2p)) -1 +1))	; cos * sines -> -1, else +1
	  (s2 (if (or c1p c2p) +1 -1))		; either are cosines -> +1 else -1
	  (s3 (if (xor c1p c2p) -1 +1))		; result is sines -> -1 else +1
	  (rl1 (reverse l1)))
      (flet ((mul1 (pr2)
	       (let ((c2 (car pr2))(m2 (cdr pr2)))
		 (if (every #'zerop1 m2)
		     (when c2p
		       (mapcar #'(lambda (pr) (cons (mul c2 (car pr))(cdr pr))) l1))
		     (let ((t1 nil)(t2 nil)(t3 nil))
		       (dolist (i1 rl1)
			 (let* ((c1 (car i1))(m1 (cdr i1))
				(cc (div (mul c1 c2) 2)))
			   (push (cons (mul s1 cc)(mapcar #'sub m1 m2)) t1)
			   (push (cons (mul s2 cc)(mapcar #'add m1 m2)) t2)))
		       (do () ((not (and t1 (list-negp (cdar t1)))))
			 (push (cons (mul s3 (caar t1))(mapcar #'neg (cdar t1))) t3)
			 (pop t1))
		       (when (and (minusp s3)(every #'zerop1 (cdar t1)))	; sin(0) ?
			 (pop t1))		; remove it.
		       (clist-add (clist-add t1 t3) t2))))))
	(cl-reduce #'clist-add l2 #'mul1)))))

(defun tlist-mul (l1 l2)
  (let ((sin1 (first l1))(cos1 (second l1))
	(sin2 (first l2))(cos2 (second l2)))
    (list (clist-add (tlist-mul1 cos1 sin2 t nil)
		     (tlist-mul1 sin1 cos2 nil t))
	  (clist-add (tlist-mul1 cos1 cos2 t t)
		     (tlist-mul1 sin1 sin2 nil nil)))))

(defun tlist-pow (l n zeros)
  (let ((sin (first l))(cos (second l)))
    (flet ((pow1 (coef m sinp)			; single {cos or sin}^n: use explicit formula
	     (let* ((s (if sinp -1 +1))		; by this point, n better be a fixnum!
		    (c (mul (if (and sinp (oddp (floor n 2))) -1 +1)
			    (div (power coef n) (power 2 (1- n)))))
		    (pow (do ((result nil)
			      (k 0)(kk n (- kk 2)))
			     ((not (plusp kk)) result)
			   (push (cons c (mapcar #'(lambda (mm)(mul kk mm)) m)) result)
			   (setq c (mul c (div (* s (- n k))(setq k (+ k 1))))))))
	       (cond ((evenp n)(list nil (cons (cons (div c 2) zeros) pow)))
		     (sinp (list pow nil))
		     (t (list nil pow))))))
      (cond ((and (null cos)(null sin))       (list nil nil))	; 0^n
	    ((zerop1 n) `( nil ((1 ,@ zeros))))
	    ((and (null (cdr cos))(null sin)) (pow1 (caar cos)(cdar cos) nil))	; cos^n
	    ((and (null (cdr sin))(null cos)) (pow1 (caar sin)(cdar sin) t))	; sin^n
	    (t (let (l^i l^n)			; Compute using "binary expansion" method.
		 (do ((bits n (lsh bits -1)))
		     ((zerop bits) l^n)
		   (setq l^i (if l^i (tlist-mul l^i l^i) l))
		   (when (oddp bits)
		     (setq l^n (if l^n (tlist-mul l^n l^i) l^i))))))))))

;;;;******************************************************************************************
;;;; Extracting Trigonometric sum coefficients.

;;; Encode a call to %sin or %cos
(defun encode-tlist (expr vs)
  (let* ((arg ($expand (cadr expr)))		; trig(arg)
	 (m (mapcar #'(lambda (v)
			(let ((mm ($coeff arg v)))
			  (setq arg (sub arg (mul mm v)))
			  mm))
		    vs))
	 (sign +1))
    (when (list-negp m)				; Make sure multiples are normalized
      (setq m (mapcar #'neg m) sign -1))
    (if (eql (caar expr) '%cos)
	(if (zerop1 arg) `(()((1 . ,m)))
	    `(((,(mul (- sign) ($sin arg)) . ,m)) ((,($cos arg) . ,m))))
	(if (zerop1 arg) `(((,sign . ,m))())
	    `(((,(mul sign ($cos arg)) . ,m)) ((,($sin arg) . ,m)))))))

;;; Transform an expression into its trigonometric coefficient list form.
(defun $trig_coeffs (expr &rest vars)
  (setq expr (totalspecdisrep expr))
  (let* ((vars (instanciate-variable-list vars expr '$alt_trig_coeffs))
	 (zeros (make-list (length vars) :initial-element 0))
	 (cache nil))
    (labels  ((gcf (expr)
		(or (mexp-lookup expr cache)
		    (cdar (push (cons expr (gcf1 expr)) cache))))
	      (gcf1 (expr)
		(let ((op (and (listp expr)(caar expr))) x y)
		  (cond ((MBAGP expr) (let ((elements (mapcar #'gcf (cdr expr))))
					(list (clist-mbag op (mapcar #'car elements))
					      (clist-mbag op (mapcar #'cadr elements)))))
			((or (null op)(not (dependsall expr vars))) `(()((,expr . ,zeros))))
			((memq op '(%SIN %COS))  (encode-tlist expr vars))
			((eq op 'MPLUS)  (cl-reduce #'tlist-add (cdr expr) #'gcf))
			((eq op 'MTIMES) (cl-reduce #'tlist-mul (cdr expr) #'gcf))
			((and (eq op 'MEXPT)	; x^y Check that we can actually compute:
			      (setq x (gcf (second expr)) y (third expr))
			      (and (integerp y)(plusp y)))	; need int y >0
			 (tlist-pow x y zeros))
			(t `(()((,expr . ,zeros))))))))
      (mlist* (mlist* '$%TRIG vars)(map-mlist (mapcar #'map-mlist (gcf expr)))))))

(defun untlist (tlist vars)
  (flet ((un1 (list trig)
	   (flet ((un2 (e)(mul (cadr e)(cons-exp trig (multl (cddr e) vars)))))
	     (addn (mapcar #'un2 list) t))))
    (addn (mapcar #'un1 tlist '(%sin %cos)) t)))

;;;********************************************************************************
;;; SERIES & TAYLOR
;;; Given an expression, a variable and an order, compute the coefficients of the
;;; expansion of the expression about variable=0 to order ORDER.
;;;  -> [[%series,variable,order],[c,p],[c',p'],...]
;;;  or [[%taylor,variable,order],[c,p],[c',p'],...]
;;; The difference is that TAYLOR computes the Taylor expansion, whereas
;;; SERIES only carries out the expansion over arithmetic functions (+,*,exp) and thus 
;;; is significantly faster.
;;;********************************************************************************

(defun $taylor_coeffs (expr var order)
  (setq expr (totalspecdisrep expr))
  (let ((var (car (instanciate-variable-list (list var) expr '$taylor_coeffs 1))))
    (labels ((make1 (expr)
	       (cond ((mbagp expr)      (clist-mbag (caar expr) (mapcar #'make1 (cdr expr))))
		     ((freeof var expr) (list (list expr 0)))
		     (t (let* ((r ($taylor expr var 0 order))
			       (ohdr (car r))
			       (hdr (list (first ohdr)(second ohdr)(third ohdr)(fourth ohdr))))
			  (if (eq (second r) 'ps)
			      (mapcar #'(lambda (p)
					  (list (specdisrep (cons hdr (cdr p)))
						(cons-exp 'rat (caar p)(cdar p))))
				      (cddddr r))
			      (list (list (specdisrep (cons hdr (cdr r))) 0))))))))
      (mlist* (mlist* '$%TAYLOR var order nil)(map-mlist (make1 expr))))))

;;;;******************************************************************************************
;;;; SLIST Arithmetic.
;;; The addition & multiplication of polynomial arithmetic are used.

;;; compute the N-th power of S through ORDER.
(defun slist-pow (s n order)
  (when s
    (let* ((m (cadar s))
	   (nm (mul n m))
	   (s_m (caar s))
	   (p (list (list (power s_m n) nm))))	; 1st term of result
      (if (null (cdr s))			; Single term
	  (or (great nm order) p)		; then trivial single term (unless high order)
	  (let* ((g (cl-reduce #'$gcd s #'(lambda (x)(sub (cadr x) m))))
		 (kmax (div (sub order nm) g)))
	    (do ((k 1 (1+ k)))
		((great k kmax) (nreverse p))
	      (let ((ff (div (add 1 n) k))
		    (trms nil))
		(dolist (s (cdr s))
		  (let ((i (div (sub (cadr s) m) g)))
		    (when (lessthan k i)(return))
		    (let ((e (member (add nm (mul (sub k i) g)) p :key #'cadr :test #'like)))
		      (when e			; multthru limits expression depth
 			(push ($multthru (mul (sub (mul i ff) 1) (car s) (caar e))) trms)))))
		(let ((pk ($multthru (div (addn trms t) s_m))))
		  (unless (zerop1 pk)
		    (push (list pk (add nm (mul k g))) p))))))))))

;;;;******************************************************************************************
;;;; Extracting Series Coefficients.

(defun $series_coeffs (expr var order)
  (setq expr (totalspecdisrep expr))
  (let ((v (car (instanciate-variable-list (list var) expr '$series_coeffs 1))))
    (setq v ($ratdisrep v))
    (labels ((mino (expr)			; Find minumum power of V in expr (for mult)
	       (let ((op (and (listp expr)(caar expr))))
		 (cond ((like expr v)   1)	; Trivial case: expr is V itself
		       ((or ($atom expr)(freeof v expr)) 0)	; `constant' case
		       ((member op '(MPLUS MLIST MEQUAL $MATRIX))
			(cl-reduce #'$min (cdr expr) #'mino))
		       ((eq op 'MTIMES)  (cl-reduce #'add (cdr expr) #'mino))
		       ((and (eq op 'MEXPT)($numberp (third expr)))	; can we compute?
			(mul (mino (second expr)) (third expr)))
		       (t 0))))			; oh, well, Treat it as constant.
	     (gcf (expr order)
	       (let ((op (and (listp expr)(caar expr))))
		 (cond ((like expr v)   `((1 1)))	; Trivial case: expr is V itself
		       ((or ($atom expr)(freeof v expr)) `((,expr 0)))	; `constant' case
		       ((MBAGp expr)
			(clist-mbag op (mapcar #'(lambda (el)(gcf el order))(cdr expr))))
		       ((eq op 'MPLUS)
			(cl-reduce #'clist-add (cdr expr) #'(lambda (el)(gcf el order))))
		       ((eq op 'MTIMES)
			(let* ((ms (mapcar #'mino (cdr expr)))
			       (mtot (addn ms t))
			       (prod '((1 0))))
			  (unless (great mtot order)
			    (do ((terms (cdr expr)(cdr terms))
				 (m ms (cdr m)))
				((null terms) prod)
			      (let ((term (gcf (car terms) (sub (add order (car m)) mtot))))
				(setq prod (clist-mul term prod order)))))))
		       ((and (eq op 'MEXPT)($numberp (third expr)))	; can we compute?
			(slist-pow (gcf (second expr) order)(third expr) order))
		       (t `((,expr 0)))))))	; just treat it as constant.
      (mlist* (mlist* '$%SERIES v order nil)(map-mlist (gcf expr order))))))

(defun unslist (clist vars)
  ($trunc (unclist clist vars)))

;;;********************************************************************************
;;; Find the coefficient associated with keys (powers or multiples) in the 
;;; coefficient list clist.
(defun $get_coef (clist &rest keys)
  (let ((sublist (case (and ($listp clist)($listp (cadr clist))(cadr (cadr clist)))
		   (($%POLY $%SERIES $%TAYLOR) (cddr clist))
		   ($%TRIG (case (car keys)
			     (($SIN %SIN) (cdr (third clist)))
			     (($COS %COS) (cdr (fourth clist)))
			     (otherwise (merror "First KEY must be SIN or COS"))))
		   (otherwise (merror "Unknown coefficient list type: ~M" clist)))))
    (or (cadar (member keys sublist :test #'alike :key #'cddr)) 0)))

;;; Reconstruct a macsyma expression from a coefficient list.
(defun $uncoef (cl)
  (let ((spec (and ($listp cl)(second cl))))
    (case (and ($listp spec)(second spec))
      ($%POLY  (unclist (cddr cl) (cddr spec)))
      (($%SERIES $%TAYLOR) (unslist (cddr cl) (cddr spec)))
      ($%TRIG   (untlist (mapcar #'cdr (cddr cl)) (cddr spec)))
      (otherwise (merror "UNCOEF: Unrecognized COEFFS form: ~M" cl)))))

;;;********************************************************************************
;;; Partition a polynomial, trig series or series into those terms whose 
;;; powers (or multiples) pass a certain test, and those who dont.
;;; Returns the pair [passed, failed].
;;; The TEST is applied to the exponents or multiples of each term.

(defun partition-clist (list test)
  (cond ((null test) (values nil list))
	((eq test T) (values list nil))
	(t  (let ((pass nil)(fail nil))
	      (dolist (item list)
		(if (is-boole-check (mapply test (cddr item) '$partition_test))
		    (push item pass)
		    (push item fail)))
	      (values (nreverse pass)(nreverse fail))))))

(defun $partition_poly (expr test &rest vars)
  (let* ((clist (apply #'$coeffs expr vars))
	 (vars (cddr (second clist))))
    (multiple-value-bind (p f)(partition-clist (cddr clist) test)
      (mlist* (unclist p vars)(unclist f vars) nil))))

(defun $partition_trig (expr sintest costest &rest vars)
  (let* ((tlist (apply #'$trig_coeffs expr vars))
	 (vars (cddr (second tlist))))
    (multiple-value-bind (sp sf)(partition-clist (cdr (third tlist)) sintest)
      (multiple-value-bind (cp cf)(partition-clist (cdr (fourth tlist)) costest)
	(mlist* (untlist (list sp cp) vars) (untlist (list sf cf) vars) nil)))))

(defun $partition_series (expr test var order)
  (let* ((clist ($series_coeffs expr var order))
	 (var (caddr (second clist))))
    (multiple-value-bind (p f)(partition-clist (cddr clist) test)
      (mlist* (unslist p var)(unslist f var) nil))))

(defun $partition_taylor (expr test var order)
  (let* ((clist ($taylor_coeffs expr var order))
	 (var (caddr (second clist))))
    (multiple-value-bind (p f)(partition-clist (cddr clist) test)
      (mlist* (unslist p var)(unslist f var) nil))))

--- NEW FILE: curvture2.mac ---
/*  written by Gosei Furuya <go....@gm...> 
#  This program is free software; you can redistribute it and/or modify        
#  it under the terms of the GNU General Public License as published by        
#  the Free Software Foundation; either version 2 of the License, or           
#  (at your option) any later version.                
*/


add_tan(_m):=block([_u:[]],_u:cross(col(_m,1),col(_m,2)),addcol(_m,_u))$
cross(_a,_b):=block([_u1,_v1],_u1:makelist(_a[j,1],j,1,3),_v1:makelist(_b[j,1],j,1,3),[_u1[2]*_v1[3]-_u1[3]*_v1[2],_u1[3]*_v1[1]-_u1[1]*_v1[3],_u1[1]*_v1[2]-_u1[2]*_v1[1]])$
 matrix_element_mult:lambda([x,y],x@y)$ 
cross2(_u1,_v1):=block([_s],_s:[_u1[2]*_v1[3]-_u1[3]*_v1[2],_u1[3]*_v1[1]-_u1[1]*_v1[3],_u1[1]*_v1[2]-_u1[2]*_v1[1]])$
make_tan():=block([_l:[],_l1,_l2,_l3,_l4,_ll,_p,_p4,_m],
	_l:map(lambda([x],diff(x,coords[1])),translist),
	_l1:map(lambda([x],x^2),_l),
	_p:ratsimp(trigsimp(apply("+",_l1))),
	_m:matrix(1/sqrt(_p)*_l),
	_ll:map(lambda([x],diff(x,coords[2])),translist),
	_l3:cross2(_l,_ll),
	_l4:map(lambda([x],x^2),_l3),
	_p4:ratsimp(trigsimp(apply("+",_l4))),
	_l2:cross2(1/sqrt(_p4)*_l3,1/sqrt(_p)*_l),
	addrow(_m,_l2,1/sqrt(_p4)*_l3))$

--- NEW FILE: diff_form.mac ---
load("cartan_new.lisp")$
infix("@")$
infix("&")$
infix("|")$

/*matrix_element_mult:lambda([x,y],x@y)$*/
load("hodge_test3.mac")$
load("f_star_test4.mac")$
load("helpfunc.mac")$
load("coeflist.lisp")$
load("format.lisp")$
load("diag")$
load("poisson.mac")$
load("frobenius.mac")$
load("curvture2.mac")$
/*
fstar_with_clf([phi,theta],[sin(phi)*cos(theta),sin(phi)*sin(theta),cos(phi)],
 	(pu1:sin(phi)*cos(theta),pu2:sin(phi)*sin(theta),pu3:cos(phi),
		h_st(pu1*d(pu2)@d(pu3)++pu2*d(pu3)@d(pu1)+pu3*d(pu1)@d(pu2))))$
*/

--- NEW FILE: example.txt ---
A. Vector analysys I 

A1.vector to form
in this package, a vector shoud be transformed to a form.
for example,[a,b,c]<---->a*Dx+b*Dy+c*Dz (so called work form)
[a,b,c]<---->a*Dy@Dz+b*Dz@Dx+c*Dx@Dy (so called fluid form)
vtof1,vtof2 need the scale_factor,so not to use f_star(),to use fstar_with_clf()
or under global environment ( see lorentz_example.txt)
 
(%i17) fstar_with_clf([x,y,z],[x,y,z],vtof1([a,b,c]));
(%o17)                        c Dz + b Dy + a Dx
(%i18) fstar_with_clf([x,y,z],[x,y,z],vtof2([a,b,c]));
(%o18)                    a Dy Dz - b Dx Dz + c Dx Dy
(%i19) fstar_with_clf([r,th,phi],[r*sin(phi)*cos(th),r*sin(phi)*sin(th),r*cos(phi)],
vtof1([a,b,c]));
(%o19)                b Dth sin(phi) r + c Dphi r + a Dr
(%i20) fstar_with_clf([r,th,phi],[r*sin(phi)*cos(th),r*sin(phi)*sin(th),r*cos(phi)],
vtof2([a,b,c]));
                                 2                       b Dphi Dr
(%o20)      a Dphi Dth sin(phi) r  + c Dr Dth sin(phi) - ---------
                                                         sin(phi)
(%i21) fstar_with_clf([r,th,phi],[r*sin(phi)*cos(th),r*sin(phi)*sin(th),r*cos(phi)],
scale_factor);
(%o21)                        [1, sin(phi) r, r]

A2. grad,rot or curle,div,laplacian
grad(f) is d(f),if A is vector,a:vtof1(A) then rot(A) is h_st(d(a)),or nest2([h_st,d],a)
h_st() is hodge star operator,div(f) is h_st(d(h_st(f))),laplacian(f) is d(h_st(d(f)))
or-1(d(antid(a))+antid(d(a)))s ok generally

1.gradiant
(%i22) fstar_with_clf([x,y,z],[x,y,z],(depends(f,[x,y,z]),d(f)));
                                df      df      df
(%o22)                       Dz -- + Dy -- + Dx --
                                dz      dy      dx
(%i23) fstar_with_clf([r,th,phi],[r*sin(phi)*cos(th),r*sin(phi)*sin(th),r*cos(phi)]
,(depends(f,[r,th,phi]),d(f)));
                              df       df         df
(%o23)                    Dth --- + Dr -- + Dphi ----
                              dth      dr        dphi 

2.rotate
(%i24) fstar_with_clf([x,y,z],[x,y,z],(depends([a,b,c],[x,y,z]),aa:vtof1([a,b,c]),
h_st(d(aa))));
                 db      da      dc      da      dc      db
(%o24)           -- Dz - -- Dz - -- Dy + -- Dy + -- Dx - -- Dx
                 dx      dy      dx      dz      dy      dz
(%i25) format(%,%poly(Dx,Dy,Dz),factor);
                   db   da        dc   da        dc   db
(%o25)            (-- - --) Dz - (-- - --) Dy + (-- - --) Dx
                   dx   dy        dx   dz        dy   dz
(%i29) fstar_with_clf([r,th,phi],[r*sin(phi)*cos(th),r*sin(phi)*sin(th),r*cos(phi)],
(depends([a,b,c],[r,th,phi]),aa:vtof1([a,b,c]),h_st(d(aa))))$

(%i30) format(%,%poly(Dr,Dth,Dphi),factor);
             db                           da
       Dphi (-- sin(phi) r + b sin(phi) - ---)
             dr                           dth                  dc          da
(%o30) --------------------------------------- - Dth sin(phi) (-- r + c - ----)
                      sin(phi)                                 dr         dphi
                                               db                          dc
                                          Dr (---- sin(phi) + b cos(phi) - ---)
                                              dphi                         dth
                                        - -------------------------------------
                                                       sin(phi) r

3.divergent
(%i31) fstar_with_clf([x,y,z],[x,y,z],(depends([a,b,c],[x,y,z]),aa:vtof1([a,b,c]),
nest2([d,h_st],aa)));
                    dc            db            da
(%o31)              -- Dx Dy Dz + -- Dx Dy Dz + -- Dx Dy Dz
                    dz            dy            dx
(%i32) fstar_with_clf([x,y,z],[x,y,z],(depends([a,b,c],[x,y,z]),aa:vtof1([a,b,c]),
nest2([h_st,d,h_st],aa)));
                                 dc   db   da
(%o32)                           -- + -- + --
                                 dz   dy   dx

(%i33) fstar_with_clf([r,th,phi],[r*sin(phi)*cos(th),r*sin(phi)*sin(th),r*cos(phi)],
(depends([a,b,c],[r,th,phi]),aa:vtof1([a,b,c]),nest2([h_st,d,h_st],aa)));
                                   db         dc
                                   ---       ----
                   c cos(phi)      dth       dphi   2 a   da
(%o33)             ---------- + ---------- + ---- + --- + --
                   sin(phi) r   sin(phi) r    r      r    dr

4.laplacian
(%i34) fstar_with_clf([x,y,z],[x,y,z],(depends(f,[x,y,z]),nest2([d,h_st,d],f)));
                            2              2              2
                           d f            d f            d f
(%o34)            Dx Dy Dz --- + Dx Dy Dz --- + Dx Dy Dz ---
                             2              2              2
                           dz             dy             dx
(%i36) fstar_with_clf([x,y,z],[x,y,z],(depends(f,[x,y,z]),nest2([h_st,d,h_st,d],f)));
                                 2     2     2
                                d f   d f   d f
(%o36)                          --- + --- + ---
                                  2     2     2
                                dz    dy    dx
(%i37) fstar_with_clf([r,th,phi],[r*sin(phi)*cos(th),r*sin(phi)*sin(th),r*cos(phi)],
(depends(f,[r,th,phi]),nest2([h_st,d,h_st,d],f)));
                                           2           2
                                          d f         d f
                 df    df                 ----       -----
               2 --   ---- cos(phi)          2           2    2
                 dr   dphi                dth        dphi    d f
(%o37)         ---- + ------------- + ------------ + ----- + ---
                r                2       2       2     2       2
                       sin(phi) r     sin (phi) r     r      dr
(%i38) fstar_with_clf([r,th,phi],[r*sin(phi)*cos(th),r*sin(phi)*sin(th),r*cos(phi)],
(depends(f,[r,th,phi]),nest2([d,h_st,d],f)));
                    2
                   d f           2                 df
(%o38) Dphi Dr Dth --- sin(phi) r  + 2 Dphi Dr Dth -- sin(phi) r
                     2                             dr
                   dr
                                                2
                                               d f
                                   Dphi Dr Dth ----
                    2                             2
                   d f                         dth                 df
    + Dphi Dr Dth ----- sin(phi) + ---------------- + Dphi Dr Dth ---- cos(phi)
                      2                sin(phi)                   dphi

A3. Example for vorticity
problem: in R3 we think velocity field v=rot(A),div(A)=0,so A is vector potential.
prove voricity w=rot(v) satisfy ¦¤A= -w
a solution: as restricted 2 spaces,A can be posed,a1*Dx+a2*Dy.
w=rot(rot(A)),then
(%i40) fstar_with_clf([x,y,z],[x,y,z],(depends([a1,a2],[x,y]),A:a1*Dx+a2*Dy,
nest2([h_st,d,h_st,d],A)));
                      2         2          2          2
                     d a2      d a1       d a2       d a1
(%o40)             - ---- Dy + ----- Dy + ----- Dx - ---- Dx
                       2       dx dy      dx dy        2
                     dx                              dy

from condition ,div(A)=0
(%i41) fstar_with_clf([x,y,z],[x,y,z],(depends([a1,a2],[x,y]),A:a1*Dx+a2*Dy,d(h_st(A))));
                          da2            da1
(%o78)                    --- Dx Dy Dz + --- Dx Dy Dz
                           dy             dx
d(a2)/dy+d(a1)/dx=0 ,then d^2(a1)/dxdy= -d^2(a2)/dx^2,d^2(a2)/dxdy=-d^2(a1)/dy^2

¦¤A is -(d(antid(A))+antid(d(A))
(%i42) fstar_with_clf([x,y,z],[x,y,z],(depends([a1,a2],[x,y]),A:a1*Dx+a2*Dy,-(nest2([d,antid],A)+nest2([antid,d],A))));
                      2         2         2         2
                     d a2      d a2      d a1      d a1
(%o42)               ---- Dy + ---- Dy + ---- Dx + ---- Dx
                       2         2         2         2
                     dy        dx        dy        dx
then 
¦¤A= -w 
but in this case,old fashion style is good with ¦¤A=grad(div(A)-rot(rot(A))

B.Changing integral variable more than two variables

B1.this is from ¦£(z)*¦£(1-z)
changing integral variable more than two variables,such like
(%i55) f_star([u,v],(s:u/(v+1),t:u*v/(v+1),exp(-(s+t))*s^(-z)*t^(z-1)*d(s)@d(t)))$
(%i57) format(%o55,%poly(u,v),factor);
                                          u v      u
                                       - ----- - -----
                               z - 1     v + 1   v + 1
                        Du Dv v      %e
(%o57)                  ------------------------------
                                    v + 1
this integral,v,[0,inf],u,[0,inf] is %pi/sin(%pi) (=beta(z,1-z))

B2. from entrance exam. of graduate course in japan
(1)integrate x+y+z=u,y+z=u*v,z=u*v*w,on 0<u<1,0<v<1,0<w<1 by pullback (x^2+y^2+z^2)dxdydz

(%i43) solve([x+y+z=u,y+z=u*v,z=u*v*w],[x,y,z]);
(%o43)           [[x = u (1 - v), y = u v (1 - w), z = u v w]]
(%i45) f_star([u,v,w],(x:u*(1-v),y:u*v*(1-w),z:u*v*w,(x^2+y^2+z^2)*d(x)@d(y)@d(z)));
                      2     2  2  2    2  2        2    2        2
(%o45)      Du Dv Dw u  v (u  v  w  + u  v  (1 - w)  + u  (1 - v) )
	
(%i46) defint(defint(defint(u^2*v*(u^2*v^2*w^2+u^2*v^2*(1-w)^2+u^2*(1-v)^2),w,0,1)
,v,0,1),u,0,1);
(%o46) 1/20
(%o47) kill(x,y,z);

slightly easy problem
(2) (d^2/dx^2+d^2/dy^2)(log(sqrt(x^2+y^2)))

(%i38) fstar_with_clf([x,y],[x,y],(d(h_st(d(log(sqrt(x^2+y^2)))))))$
(%i39) format(%,%poly(Dx,Dy),factor);
(%o39)                                 0

(3) Integrate (x*Dx@Dy) on x^2/a^2+y^2/b^2<=1,x>=0,y>=0

(%i40) fstar_with_clf([x,y],[x,y],(depends([a1,a2],[x,y]),d(a1*d(x)+a2*d(y))));
                             da2         da1
(%o40)                       --- Dx Dy - --- Dx Dy
                              dx          dy
(%i41) fstar_with_clf([x,y],[x,y],(d(1/2*x^2*Dy)));
(%o41)                             Dx Dy x
1/2*x^2*Dy x=a*cos(th),y=b*sin(th)  on boundary you may use line integrate.

B3. Integral

we must pay attention to effects of some singularity and of topology (De Rham cohomology )
this points is most difficult to calculate multivariable integral in CAS,I think.
see Frankel book "The Geometry of PHYSICS" p162 5.5 Finding Potential.

(Pontrjaguin index)
(%i74) f_star([phi,theta],(pu1:sin(phi)*cos(theta),pu2:sin(phi)*sin(theta),pu3:cos(phi),
    pu1*d(pu2)@d(pu3)+pu2*d(pu3)@d(pu1)+pu3*d(pu1)@d(pu2)))$
(%i76) trigsimp(%o74);
(%o76)                       Dphi Dtheta sin(phi)
So this map induce (x,y)-->(pu1,pu2,pu3)  IxI--->S2,boundary 1X1 -->one point.
%o76 is volume form  in S2,if we can pull back IxI,we integrate this and say
this integral value is n*4pai. n is a integer.
(%i79) f_star([x,y],(depends([p1,p2,p3],[x,y]),p1*d(p2)@d(p3)+p2*d(p3)@d(p1)+p3*d(p1)@d(p2)));

                dp2 dp3         dp2 dp3              dp1 dp3         dp1 dp3
(%o79) p1 (Dx Dy --- --- - Dx Dy --- ---) + p2 (Dx Dy --- --- - Dx Dy --- ---)
                 dx  dy          dy  dx               dy  dx          dx  dy
                                                    dp1 dp2         dp1 dp2
                                           + (Dx Dy --- --- - Dx Dy --- ---) p3
                                                    dx  dy          dx  dy

(%o79) is written by using vector,(p1,p2,p3).{(d/dx(p1,p2,p3)Xd/dy(p1,p2,p3))*Dx@Dy}
integrate with IxI,this is homotopy invariant.
generaly speaking,in exchaing coordinates, vector equations are not always invariant but differntial forms are invariant. so we can calculate forms with
some suitable or fit coordinate,then by pulling back it  with general ones. 

C.Integral factor (always possible locally,but not global)

locally necessary and  sufficient condition is d(w)@w=0. so if f^-1 is integral factor
,d(w)=d(log|f|)@w . usually we use frobenius method(or theorem) but I propose a new method
with quantization.
 
(%i72) fstar_with_clf([x,y,z],[x,y,z],trans_toexact(y*z*d(x)+x*z*d(y)+d(z)));
                         2            2
                  Dz u1 y  z + Dz u2 x  z - Dx u3 y - Dy u3 x
(%o72)            -------------------------------------------
                               2  2       2  2
                           u1 y  z  + u2 x  z  + u3
u1,u2,u3 is scale factor.
integral factor do'nt depends on metric,so u1->0,u2->0. or u2->0,u3->0 so on.
this is a quantization from clifford algebra to grassmann algebra.

(%i74) limit(limit(%o72,u1,0),u2,0);
(%o74)                           - Dx y - Dy x
(%i75) limit(limit(%o72,u3,0),u2,0);
                                      Dz
(%o75)                                --
                                      z
now w is  y*z*d(x)+x*z*d(y)+d(z),by %o74 d(w)=(- Dx y - Dy x)@w,
d(log|f|)=-( Dx y + Dy x)=-d(x*y),f^(-1)=%e^(x*y) is integral factor.
by %o75 d(w)=Dz/z @w,d(log|f|)=Dz/z,f^(-1)=1/z is also integral factor.
but %e^(x*y)is better than 1/z. see Flanders P92,93.

D. Poission bracket (symplectic form,symplectic manifold)

for example symplectic manifold dim 6,it has 2form named to symplectic forms.
 1. w is closed form 2. w^3==w@w@w is not 0 everywhere .
(if dim 2*m,w^mis not0)
(%i39) f_star([p1,q1,p2,q2,p3,q3],(omega:Dp1@Dq1+Dp2@Dq2+Dp3@Dq3,omega@omega@omega));
(%o39)                    - 6 Dp1 Dp2 Dp3 Dq1 Dq2 Dq3
(%i40) f_star([p1,q1,p2,q2,p3,q3],(omega:Dp1@Dq1+Dp2@Dq2+Dp3@Dq3,d(omega)));
(%o40)                                 0
so this omega is simplectic form.
Poission bracket is definded as coordinates free
	inner(Xf,(inner(Xg,omega))),and inner(Xg,omega) is d(g).
(%i42) f_star([p1,p2,p3,q1,q2,q3],(depends([g],[p1,p2,p3,q1,q2,q3]),omega:Dp1@Dq1+Dp2@Dq2
	+Dp3@Dq3,inner(diff(g,q1)*Dp1-diff(g,p1)*Dq1+diff(g,q2)*Dp2-diff(g,p2)*Dq2+diff(g,q3)*Dp3
	-diff(g,p3)*Dq3,omega)));
               dg        dg        dg        dg        dg        dg
(%o42)     Dq3 --- + Dq2 --- + Dq1 --- + Dp3 --- + Dp2 --- + Dp1 ---
               dq3       dq2       dq1       dp3       dp2       dp1
because poission bracket is defined
by inner(Xf,d(g)) in this package p_braket().

(%i55) f_star([p1,q1,p2,q2,p3,q3],(depends([f,g],[p1,q1,p2,q2,p3,q3]),p_braket(f,g)));
           df  dg    df  dg    df  dg    df  dg    df  dg    df  dg
(%o55)     --- --- + --- --- + --- --- - --- --- - --- --- - --- ---
           dp3 dq3   dp2 dq2   dp1 dq1   dq3 dp3   dq2 dp2   dq1 dp1

Jacobi identity

(%i63) f_star([p1,q1,p2,q2,p3,q3],(depends([f,g,h],[p1,q1,p2,q2,p3,q3]),
p_braket(p_braket(f,g),h)+p_braket(p_braket(g,h),f)+p_braket(p_braket(h,f),g)))$
(%i64) ratsimp(%);
(%o64)                                 0
(%i67) f_star([p1,q1,p2,q2,p3,q3],(depends([f,g,h],[p1,q1,p2,q2,p3,q3]),
	p_braket(f,g*h)-g*p_braket(f,h)-h*p_braket(f,g)))$
(%i68) ratsimp(%);
(%o68)                                 0

Is canonical? [p,q]--->[P,Q]
(%i69) f_star([p,q],(P:q*cot(p),Q:log(sin(p)/q),d(P)@d(Q)))$
(%i70) trigsimp(%);
(%o70)                               Dp Dq
DP@DQ=Dp@Dq show canonical coordinates each other.
if H is hamiltonian, d(q)/dt=p_braket(q,H),d(p)/dt=p_braket(p,H)
(%i71) f_star([p1,q1,p2,q2,p3,q3],(depends([f,g,h],[p1,q1,p2,q2,p3,q3]),
p_braket(f*g,h)-f*p_braket(g,h)-g*p_braket(f,h)))$ is 0
if p_braket(g,H)=0 and p_braket(f,H)=0 (f,g is first integral),then
p_braket(f*g,h)=0 so f*g is a first integral too.
from jacobi identity p_braket(f,g) is also first integral.

--- NEW FILE: f_star_test4.mac ---
/* written by Gosei Furuya <go....@gm...> 
#  This program is free software; you can redistribute it and/or modify        
#  it under the terms of the GNU General Public License as published by        
#  the Free Software Foundation; either version 2 of the License, or           
#  (at your option) any later version.                
*/


infix("@");
infix("&");

f_star(newcoords,'a_form):=
	block([dim,i:1,coords,extdim:2,basis,extsub,extsubb,pu_],
	dim:length(newcoords),array(pu_,dim),
	mode_declare([basis,extsub,extsubb],any),
	coords:newcoords,
	for i thru dim do
	(
	pu_[i]:concat(D,newcoords[i])
	),basis:makelist(pu_[i],i,1,dim),
	extsub[1]:[],
	for i thru dim do
	 (
  	 extsub[i+1]:cons(basis[i]=-basis[i],extsub[i]),
	
 	 extsubb[i]:cons(basis[i]=0,extsub[i])),ev(a_form)
	);

fstar_with_clf(newcoords,n_table,'a_form):=
	block([dim,i:1,coords,extdim:2,basis,extsub,extsubb,extsubb2,
	norm_table,scale_factor,volume,a_,b_,x_,pu_],
	mode_declare([basis,extsub,extsubb,extsubb2],any),
	dim:length(newcoords),
	coords:newcoords,array(pu_,dim),
	for i thru dim do
	 (pu_[i]:concat(D,newcoords[i])
	),
	basis:makelist(pu_[i],i,1,dim),
	extsub[1]:[],
	for i thru dim do
	 (
  		extsub[i+1]:cons(basis[i]=-basis[i],extsub[i]),
 	 	extsubb[i]:cons(basis[i]=0,extsub[i])),
	norm_table:clif_norm(n_table,coords),
	a_:solve(x_^2-apply("*",norm_table),[x_]),
	volume:rhs(a_[2]),volume:1/volume,
	scale_factor:[],
	for i:1 thru dim do 
	( a_:solve(x_^2-1/norm_table[i],[x_]),
	scale_factor:cons(rhs(a_[2]),scale_factor)
	),
	scale_factor:reverse(scale_factor),
	for i:1 thru dim do
 	( extsubb2[i]:cons(basis[i]=norm_table[i]/basis[i],extsub[i])
	),
	ev(a_form)
	);

clif_norm(list_,coords_):=
	block([dim,_p,coords,cliffordtype,ntable:[],_l:[]],
	coords:coords_,
	dim:length(coords),
	cliffordtype:makelist(1,i,1,dim),
	for i:1 thru dim do 
	(_l:map(lambda([x],diff(x,coords[i])),list_),
    		_l:map(lambda([x],x^2),_l),_p:ratsimp(trigsimp(apply("+",_l))),
		ntable:endcons(cliffordtype[i]/_p,ntable)),
	ntable
	);

/*inner[_f](_g) */
inner(_f,_g):= block([_a,_b:[],_r],
      _a:expand(_f),
      for i:1 thru dim do (_b:endcons(ratcoef(_a,basis[i]),_b)),
      _r:_b | _g);

/*Lie[_f1](_g1)    Lie differntial operator*/
Lie(_f1,_g1):=d(inner(_g1,_f1))+inner(d(_g1),_f1);




 nest2(_f,_x):=block([_a:[_x],i],if listp(_f) then (
	 _f:reverse(_f),for i:1 thru length(_f) do(_a:map(_f[i],_a))) 
	else (_a:map(_f,_a)),_a[1])$
 nest3(_f,_x,_n):=block([_a,i],_a:[_x],for i:1 thru _n do (_a:map(_f,_a)),_a)$ 

--- NEW FILE: format.lisp ---
;;; -*- Mode: LISP; Syntax: Common-lisp; Package: Maxima; Base: 10 -*-
;;;>******************************************************************************************
;;;>   Software developed by Bruce R. Miller (mi...@ca...)
;;;>   using Symbolics Common Lisp (system 425.111, ivory revision 4)
;;;>   at NIST - Computing and Applied Mathematics Laboratory
;;;>   a part of the U.S. Government; it is therefore not subject to copyright.
;;;>******************************************************************************************
;;;;******************************************************************************************
;;;; FORMAT: Package for restructuring expressions in Macsyma
;;;;******************************************************************************************

;;(in-package 'climax)				; for Macsyma Inc, Macsyma

;;; To run in Schelter's Maxima comment the above and uncomment these:
(in-package :maxima)
(defmacro mlist* (arg1 &rest more-args) `(list* '(mlist simp) ,arg1 ,@more-args))
(defun mrelationp (expr)
  (and (listp expr)(member (caar expr) '(MEQUAL MNOTEQUAL MGREATERP MLESSP MGEQP MLEQP))))

;;;;******************************************************************************************
;;; format(expr,template,...)
;;; Formats EXPR according to the TEMPLATEs given: 

(defvar *template* nil "The current template")
(defvar *templates* nil "The current template chain")
(defvar *subtemplates* nil "Current template's subtemplates")

(defun $format (expr &rest templates) (format-from-chain expr templates))

;; format according to chain in *templates*
(defun format-from-chain (expr &optional (*templates* *templates*))
  (if (null *templates*) expr
      (format-one expr (pop *templates*))))

;; format according to tmp, then pieces according to *templates*
(defun format-one (expr tmp)
  (multiple-value-bind (*template* formatter parms *subtemplates*)(parse-template tmp)
    (cond (formatter (apply #'mfuncall formatter expr parms))
	  ((or (symbolp tmp)			; Apply SPEC as function, if it CAN be
	       (and (listp tmp)(or (eq (caar tmp) 'lambda)(member 'array (cdar tmp)))))
	   (format-from-chain (let ((*templates* nil)) (mfuncall tmp expr))))
	  (t (merror "FORMAT: template ~M must be a template or function." tmp)))))

;;; Format a `piece' of an expression, accounting for any current subtemplates.
;;; If NTH is given, use NTH subtemplate for this piece, else use next subtemplate.
;;; Account for %DITTO'd templates. 
(defun $format_piece (piece &optional nth)
  (flet ((dittop (ptrn)				; If %ditto form, return repeated template
	   (and (listp ptrn)(eq (caar ptrn) '$%ditto) (cadr ptrn))))
    (let ((subtmp (cond ((null *subtemplates*) nil)	; no piecewise templates.
			((null nth)(or (dittop (car *subtemplates*))	; next one %ditto's
				       (pop *subtemplates*)))	; otherwise, remove next one.
			((setq nth (nth (1- nth) *subtemplates*))	; nth subtemplate?
			 (or (dittop nth) nth))	; strip off possible %ditto
			((dittop (car (last *subtemplates*)))))))	; last dittos, reuse it
      (if subtmp (format-one piece subtmp)(format-from-chain piece)))))

;; Format expr according to remaining chain, but disallowing subtemplates.
(defun format-w/o-subtemplates (expr)
  (when *subtemplates*
    (merror "FORMAT: Template ~M was given subtemplates ~M" *template*
	    (mlist* *subtemplates*)))
  (format-from-chain expr))

;;; given a candidate format template, return:
;;; template name, formatter function, parameters (if any) and subtemplates (if any), 
(defun parse-template (template)
  (let (op name formatter)
    (flet ((getform (symbol)
	     (and (setq formatter (or ($get symbol '$formatter)(get symbol 'formatter)))
		  (setq name symbol))))
      (cond (($numberp template) nil)
	    ((atom template) (values (getform template) formatter nil nil))
	    ((eq (caar template) 'mqapply)	; Template w/ subtemplates AND parms
	     (when (and (listp (setq op (cadr template)))
			(getform (caar op))
			(not (member 'array (cdar op))))	; but not T[...][...]
	       (values name formatter (cdr op) (cddr template))))
	    ((getform (caar template))		; Template w/ parameters OR subtemplates
	     (if (member 'array (cdar template))
		 (values name formatter nil (cdr template))
		 (values name formatter (cdr template) nil)))))))

;;;;******************************************************************************************
;;; Defining format commands.
;;; (user defined ones go on the macsyma property list)

(defmacro def-formatter (names parms &body body)
  (let* ((names (if (listp names) names (list names)))
	 (fmtr (if (atom parms) parms
		   (make-symbol (concatenate 'string (string (car names)) "-FORMATTER")))))
    `(progn 
       ,(unless (atom parms) `(defun ,fmtr ,parms ,@body))
       ,@(mapcar #'(lambda (name) `(setf (get ',name 'FORMATTER) ',fmtr)) names))))

;;;;******************************************************************************************
;;; Subtemplate aids.

(def-formatter mlist (expr &rest elements)	; merge elements w/ following chain.
  (format-from-chain expr (append elements *templates*)))

(def-formatter $%preformat (expr &rest templates)	; preformat using template chain
  (format-w/o-subtemplates (format-from-chain expr templates)))

(def-formatter $%noop format-w/o-subtemplates)	; subtemplate filler.

;;;;******************************************************************************************
;;; Arithmetic template: eg.  a*%p(x)-b

(defun template-p (expr)			; is there a template in expr?
  (if (and (listp expr)(member (caar expr) '(mplus mtimes mexpt)))
      (some #'template-p (cdr expr))		; for arithmetic, find a `real' format in args
      (parse-template expr)))

(defun partition-arithmetic-template (op args)
  ;; Find the 1 (!) term or factor with a regular template in it.
  (let ((pat (remove-if-not #'template-p args)))	; find arg with template in it
    (when (or (null pat)(cdr pat))
      (merror "FORMAT: Pattern must contain exactly 1 template ~M" (cons (list op) args)))
    (values (car pat) (simplify (cons (list op) (remove (car pat) args))))))

(def-formatter mplus (expr &rest terms)
  (multiple-value-bind (template rest)(partition-arithmetic-template 'mplus terms)
    (add (format-one (sub expr rest) template) rest)))

(def-formatter mtimes (expr &rest factors)
  (multiple-value-bind (template rest)(partition-arithmetic-template 'mtimes factors)
    (mul (format-one (div expr rest) template) rest)))

(def-formatter mexpt (expr b p)			; b^p
  (cond ((template-p b)(power (format-one (power expr (inv p)) b) p))
	((template-p p)(power b (format-one (div ($log expr)($log b)) p)))
	((merror "FORMAT: Pattern must contain exactly 1 template ~M" (power b p)))))

;;;;******************************************************************************************
;;; Control templates

;;; IF ... ELSEIF ... ELSE 
(def-formatter $%if (expr &rest predicates)
  ($format_piece expr (do ((ps predicates (cdr ps))
			   (i 1 (1+ i)))
			  ((or (null ps)(is-boole-check (mfuncall (car ps) expr))) i))))
    
(def-formatter ($%expr $%expression)(expr)	; format arguments/operands
  (when ($atom expr)
    (merror "FORMAT %EXPR: ~M doesn't have parts" expr))
  (map1 #'$format_piece expr))

;;; Convenience templates
(def-formatter $%subst (expr &rest listofeqns)
  (format-w/o-subtemplates ($substitute (mlist* listofeqns) expr)))

(def-formatter $%ratsubst (expr &rest listofeqns)
  (autoldchk '$lratsubst)
  (format-w/o-subtemplates ($lratsubst (mlist* listofeqns) expr)))

;;;;******************************************************************************************
;;; `Bag' & Relation templates.

;;; This function tries to get OPER at the top level of EXPR.
;;; OPER must be a BAG or RELATION, as must the top layers of EXPR 
;;; (down to wherever OPER is found).
;;; The interpretation is that a list of equations is equivalent to an equation 
;;; whose rhs & lhs are lists.  (and ditto for all permutations).
(defun $coerce_bag (oper expr)
  (unless (or (mbagp expr)(mrelationp expr))
    (merror "Error: ~M is not a relation, list or array: can't be made into an ~M" expr oper))
  (setq oper (or (get oper 'opr) oper))
  (flet ((swap (op x)
	   (cons (list op)
		 (mapcar #'(lambda (l)(simplify (cons (car x) l)))
			 (transpose (mapcar #'(lambda (y)(cdr ($coerce_bag op y)))(cdr x)))))))
    (cond ((eq (caar expr) oper) expr)		; oper is already at top level.
	  ((eq (caar expr) '$matrix)		; swap levels 2 & 3 (mlist & oper), then 1&2
	   (swap oper (map1 #'(lambda (x)(swap oper x)) expr)))
	  ((eq oper '$matrix)			; swap level 1 & 2 (oper & matrix), then 2 & 3.
	   (map1 #'(lambda (l)(swap 'mlist l))(swap oper expr)))
	  (t (swap oper expr)))))		; swap 1st & 2nd levels.

(defun format-bag (expr op)
  (map1 #'$format_piece ($coerce_bag op expr)))

(def-formatter ($%eq $%equation $%rel $%relation) (r &optional (op 'mequal))
  (format-bag r op))
(def-formatter $%list (expr)  (format-bag expr 'mlist))
(def-formatter $%matrix (expr)(format-bag expr '$matrix))
;;; Note: %matrix subtemplates apply to ROWS. To target elements, use %list for rows.

;;;;******************************************************************************************
;;; Targetting templates.
;;; mostly shorthand for things which can be done using subtemplates, but more concise.

(defun format-nth (expr n)
  (unless (and ($integerp n) (plusp n)(< n (length expr)))
    (merror "FORMAT ~M: ~M doesn't have an argument #~M" *template* expr n))
  (let ((new (copy-list expr)))
    (setf (nth n new)(format-w/o-subtemplates (nth n expr)))
    (simplify new)))

(def-formatter $%arg format-nth)
(def-formatter $%lhs (expr &optional (op 'mequal))
  (format-nth ($coerce_bag op expr) 1))
(def-formatter $%rhs (expr &optional (op 'mequal))
  (format-nth ($coerce_bag op expr) 2))

(def-formatter ($%el $%element)(expr &rest indices)
  (let ((array ($copymatrix ($coerce_bag '$matrix expr))))
    (apply #'marrayset ($format_piece (apply #'marrayref array indices)) array indices)
    ...
 
[truncated message content]