#@##############################################################################
#	
#	______________________________________________________________
#
#                       THE Forms PACKAGE
#	______________________________________________________________
#
#        Copyright (c) 1991 by Yunliang Yu. All rights reserved.
#
#   This package is for differential geometers working on exterior
#   differential systems. It can be used to set up a geometric 
#   problem using moving frames, and do systematic prolongations 
#   and other calculations in Maple. Please see the help ?Forms for 
#   the usage information.
#
#   This package may be freely distributed for non-profit purposes 
#   only. It may be modified as long as no modified version is
#   distributed and the copyright notice is retained. The author 
#   is not responsible for any damage caused by using this package.
#
#   This package was first written in Maple version 4.2 and now
#   revised in maple V release 4.  (Revised for Maple V, release 5
#   by Jeanne Clelland.)  Please send bug reports, comments
#   and ideas to yu@math.duke.edu.
#
#@##############################################################################
date_:=`$Id: forms,v 1.52 1997/07/08 18:23:03 yu Exp yu $`:
`help/text/Forms` := TEXT(
`        ______________________________________________________________`,
`   `,
`                             THE Forms PACKAGE`,
`        ______________________________________________________________`,
`   `,
`            Copyright (c) 1991 by Yunliang Yu. All rights reserved.`,
` `,
`            ~~~~~~~~~~~~~~~ Usage of this package ~~~~~~~~~~~~~~~`,
`   `,
`SYNOPSIS:   `,
`   `,
`  - Form( <n1> = <e1>, <n2> = <e2>, . . . )`,
`   `,
`        a): This procedure is used to define the basic forms used in a`,
`            calculation, (then) to define the exterior derivative of a`,
`            name or to make an assignment to a name or a function.`,
`            Each parameter must be of the form <name>=<expr>; and <name>`,
`            *may* need to be enclosed in single quotes to prevent`,
`            evaluation if <name> is a function, for example,`,
`                    Form('d(x*y)'=dxy);`,
`        b): For <name> = <expr>:`,
`              if <expr>='const' or -1, then <name> is a constant.`,
`              if <expr>='form', then  <name> is a form w/ unspecified degree.`,
`              if <expr>='scalar' or 0, then <name> is scalar-valued (all names,`,
`                 by default, are scalar-valued; and this option is mainly `,
`                 used to redefine a form to be a scalar).`,
`              if <expr> is a list: [range1, range2, .., deg], it defines`,
`                 all the names <name>.i.j.k..., with i from the range1, j `,
`                 from the range2,..., to be the forms of degree deg.`,
`              When <expr> is a positive integer, the wedge degree of <name>`,
`                 is defined to be <expr>. In addition, if <name> is a 1-form,`,
`                 then "=1" part can be omitted. `,
`              Otherwise, it assigns <expr> to <name>; and moreover, if <expr>`,
`                 is of type name or function and <name> is a form, it defines `,
`                 <expr> to be a form with the same degree as <name>. etc.`,
`        c): Assignments made through Form are remembered permanently.`,
`            i.e. they are not washed out by gc().`,
`        d): The order sequence n1, n2, n3 .... is the default `,
`            ordering of the forms in wedge products.`,
`        e): Form() w/o arguments or Form(1) can be used during a session (e.g.,`,
`            after you making any assignments to basic forms) to clear those `,
`            unnecessary data in the memory to speed up the computation.`,
`   `,
`  - Forder( n1,n2,n3,...)`,
`   `,
`        a): This procedure specifies the ordering of given forms in &^(..) to `,
`            be n1,n2,n3,.... by inserting n1,n2, ... to the appropriate `,
`            places in the ordering list.`,
`        b): Forder() w/o any arguments turns off/on the ordering.`,
`   `,
`  - Wdegree(expr)`,
`   `,
`        a): It computes the wedge degree of an expression or array.`,
`        b): It returns nonhmg when the expression is a sum of forms with`,
`            different degrees.`,
`   `,
`  - &^(expr1,expr2,...)`,
`   `,
`        It calculates the wedge product of [arrays of] differential forms.`,
`        The name of this procedure can be changed using the procedure`,
`        "WedgeProduct".`,
`   `,
`  - d(expr)`,
`   `,
`        It computes the exterior derivative of an expression or array. `,
`   `,
`  - ScalarForm(expr,'f')`,
`   `,
`        a): This function is used to separate the scalar part and the form`,
`            part of an expr, it returns a list of the coeffs of the forms. `,
`        b): If the optional 'f' is specified, it is assigned a list of `,
`            the corresponding forms.`,
`   `,
`  - Simf(expr)`,
`   `,
`        a): It simplifies an expression or array involving forms.`,
`            It will collect like terms, simplify wedge products, pull out  `,
`            scalar factors, etc.`,
`        b): If Simf is called w/ additional arguments,`,
`            Simf(expr,proc,opt), then the procedure proc with additional`,
`            arguments opt is applied on each collected scalar factor,`,
`            e.g.,Simf(expr,collect,[x,y],distributed).`,
`   `,
`  - WedgeProduct(s)`,
`   `,
`        a): To select the string ``s`` as the symbol for wedge product.`,
`        b): The default symbol is ``&^``. You can change this symbol any time`,
`            you like. No data will be lost or changed.`,
`   `,
`  - L(w,expr)`,
`   `,
`        left hook an expr ``expr`` with the dual vector field of the 1-form ``w``.`,
`        (Caution: This procedure is simply minded. It basically does the`,
`        following: write expr as  w &^ form1 + others w/o w, and return form1)`,
`   `,
`  - pick(expr,w1,w2,...,wn)`,
`   `,
`        It picks out all the wedge products in expr having the forms `,
`        w1, w2, ..., wn in them, then writes the expr as a wedge product: `,
`        form1 &^ w1&^w2&^...&^wn, and returns form1.`,
`   `,
`  - mixpar(e)`,
`   `,
`        "mixpar" takes nested calls to "diff", and sorts the sequence of`,
`         differentiations  so  that they are sorted by lexicographical`,
`         ordering on the variables that are being differentiated.`,
`   `,
`     `,
`I appreciate any comments and new ideas.  Yunliang Yu  11/18/90`,
`   `,
`VERSION: $Id: forms,v 1.52 1997/07/08 18:23:03 yu Exp yu $`,
`  `
):
`help/text/forms`:=`help/text/Forms`:
`help/text/form`:=`help/text/Forms`:
`help/text/Form`:=`help/text/Forms`:
WedgeDegree_:=table('sparse'):
FormOrder_:=table(sparse):
FormOrder_[]:=1:
FormOrder_no:=1:
Form:=proc()
	local i,j,  a,b, t, RT,  nb,wd,ax:
	global ScalarForm_1_,W_sub2_, W_sub_, Wdegree, Wdf_, Wf_, d_sub_, df_,\
		FormOrder_,FormOrder_no, WedgeDegree_, date_;
	t:=table():
	if assigned(date_) then
		lprint(`The following  functions  are  available:`);
		print(`Form, &^, d, Wdegree, Simf, Forder, ScalarForm, L, `.
 		`WedgeProduct, mixpar, pick.`);
        	date_:='date_':
	fi:
	for i to nargs do
		if type(args[i],{'name','function'}) then
		   	if CKW_(args[i]) or CKW_(EVAL(args[i])) then
				lprint(`Warning: the `.i.`th argument can NOT `.
				`be designated to be an 1-form.`):
				next;
			fi:
			WedgeDegree_[args[i]]:=1:
			FormOrder_[args[i]]:=FormOrder_no: FormOrder_no:=FormOrder_no+1:
			t[i]:=`	1-form: `,args[i]:
		elif type(args[i],`=`) then
			a:=op(1,args[i]):
			b:=EVAL(op(2,args[i])):
			if not type(a,{name,function}) then
				lprint(`Warning: The left side of the `.
				i.`th argument is not a name or function.`):
				next:
			fi:
			if b='const' or b=-1 then
				if Wdf_(a) then
					print(`Warning: the `.i.`th argument `.
					`is not allowed:`,a=b): next:
				fi:
				d(a):=0:
				WedgeDegree_[a]:=evaln(WedgeDegree_[a]):
				FormOrder_[a]:=evaln(FormOrder_[a]):
				t[i]:=`	 const: `,a:
			elif b='scalar' or b=0 then
				if Wdf_(a) then
					print(`Warning: the `.i.`th argument `.
					`is not allowed:`,a=b): next:
				fi:
				WedgeDegree_[a]:=evaln(WedgeDegree_[a]):
				FormOrder_[a]:=evaln(FormOrder_[a]):
				t[i]:=`	scalar: `,a:
			elif b='form' then  a:
				if CKW_(a) or CKW_(EVAL(a)) then
					lprint(`Warning: the `.i.`th argument `.
 					`can NOT be defined in this way.`):
					next;
				fi:
				WedgeDegree_[a]:=''Wdegree(%)'':
				FormOrder_['%%']:=FormOrder_no: FormOrder_no:=FormOrder_no+1:
				t[i]:=`	  form: `,a:
			elif type(b,'integer') and b>0 then
				if CKW_(a) or CKW_(EVAL(a)) then
					lprint(`Warning: the `.i.`th argument `.
 					`can NOT be defined in this way.`):
					next;
				fi:
				WedgeDegree_[a]:=b:
				FormOrder_[a]:=FormOrder_no: FormOrder_no:=FormOrder_no+1:
				t[i]:=`	`.b.`-form: `,a:
elif type(b,list) then
	if not type(a,name) then
		lprint(`The left side of the `.i.`th argument must be a name.`):
		next:
	fi:	nb:=nops(b):
	for j to nb do
		if not type(b[j],{range,integer}) then
		      lprint(`The right side of the `.i.`th argument must be: `.
		      `[range,range,..., (optional)deg].`):break:
		fi:
	od:
	if j<nb+1 then next fi:
	if type(b[nb],range) then wd:=1:else wd:=b[nb]:nb:=nb-1:fi:
	ax:={a}:
	for j to nb do
		ax:=map(proc(a,b) if assigned(a) then RETURN(NULL) fi:
			if type(b,integer) then a.(1..b): else a.(b) fi:
			end,	ax,b[j]):
	od:
	if wd<0 then
		for j in ax do
			d(j):=0;
			WedgeDegree_[j]:=evaln(WedgeDegree_[j]);
			FormOrder_[j]:=evaln(FormOrder_[j]):
		od:
		t[i]:=`	constant:  `,``.a.` with `.nb.` indices.`:
	elif wd=0 then
		for j in ax do
			WedgeDegree_[j] := evaln(WedgeDegree_[j]);
			FormOrder_[j]:=evaln(FormOrder_[j]):
		od:
		t[i]:=`	scalar:  `,``.a.` with `.nb.` indices.`:
	else
		for j in ax do
			WedgeDegree_[j]:=wd:
			FormOrder_[j]:=FormOrder_no: FormOrder_no:=FormOrder_no+1:
        od:
        t[i]:=`	`.wd.`-forms:  `,``.a.` with `.nb.` indices.`:
    fi:
			else
if type(b,name) or (type(b,function) and op(0,b)<>`d`) then
	if Wdegree(b)=0 then
		if Wdegree(a)=0 then
			t[i]:=`	scalar: `,b:
		else
			if not type(Wdegree(a),integer) then
				b;''Wdegree(%)'': else	 Wdegree(a):
			fi:
			WedgeDegree_[b]:='%':
			if type('%',integer) then
				t[i]:=`	`.%.`-form: `,b:
			else	t[i]:=`	  form: `,b: fi:
			FormOrder_[a]:=evaln(FormOrder_[a]):
			FormOrder_[b]:=FormOrder_no: FormOrder_no:=FormOrder_no+1:
			WedgeDegree_[a]:=evaln(WedgeDegree_[a]):
		fi:
	else
		if Wdegree(a)=0 then
			lprint(`Warning: the `.i.`th argument assigns a form `.
			`to a scalar.`): next:
		fi:
		if type(b,name) then
			if not type(Wdegree(a),integer) then
				if type(Wdegree(b),integer) then
				       t[i]:=`	`.(Wdegree(b)).`-form: `,b:
				else
					t[i]:=`	form: `,b:
				fi:
			elif not type(Wdegree(b),integer) then
				WedgeDegree_[b]:=Wdegree(a):
				t[i]:=`	`.%.`-form: `,b:
			elif Wdegree(a)<>Wdegree(b) then
				lprint(`Warning: Two sides of the `.i.`th `.
				`argument have different degrees.`): next:
			else t[i]:=`	`.(Wdegree(b)).`-form: `,b:
			fi:
		else
			if type(Wdegree(a),integer) and
				type(Wdegree(b),integer) and
				Wdegree(a)<>Wdegree(b)    then
				lprint(`Warning: Two sides of the `.i.`th `.
				`argument have different degrees.`): next:
			fi:
		fi:
		if a<>b then
			FormOrder_[a]:=evaln(FormOrder_[a]):
			WedgeDegree_[a]:=evaln(WedgeDegree_[a]):
		fi:
	fi:
elif df_(b) then
	j:=op(b):
	if Wdegree(a)=0 then lprint(`Warning: The `.i.`th argument assigns `.
		`a d(..) to a scalar.`): next:
	elif a=j then lprint(`Warning: assignment `.i.` is not allowed.`):next:
	elif Wdegree(a)=1 then
		if Wdegree(j)<>0 then
			lprint(`The right side of the `.i.`th argument `.
			`should be a d(some scalar).`): next:
		fi:
		t[i]:=`	scalar: `,j:
	elif Wdegree(j)=0 then
		if not type(Wdegree(a),integer) then
			j;''Wdegree(%)'':	 else	 Wdegree(a)-1:
		fi:
		WedgeDegree_[j]:='%':
		if type('%',integer) then
     			 t[i]:=`	`.%.`-form: `,j:
		else
			 t[i]:=`	  form: `,j:
		fi:
		WedgeDegree_[d(j)]:=evaln(WedgeDegree_[d(j)]):
		FormOrder_[d(j)]:=evaln(FormOrder_[d(j)]):
		FormOrder_[a]:=evaln(FormOrder_[a]):
		FormOrder_[j]:=FormOrder_no: FormOrder_no:=FormOrder_no+1:
	elif type(Wdegree(a),integer) then
		if type(Wdegree(j),integer) then
			if Wdegree(j)<>Wdegree(a)-1 then
				lprint(`Warning: Two sides of the `.i.`th `.
				`argument have different degrees.`): next:
			fi:
   			t[i]:=`	`.(Wdegree(j)).`-form: `,j:
		else
			WedgeDegree_[j]:=Wdegree(a)-1:
			t[i]:=`	`.%.`-form: `,j:
		fi:
	else
		t[i]:=`	form: `,j:
	fi:
	FormOrder_[a]:=evaln(FormOrder_[a]):
	WedgeDegree_[a]:=evaln(WedgeDegree_[a]):
else
	if Wdegree(a)=0 and Wdegree(b)<>0 then
		lprint(`Warning: The `.i.`th argument assigns a form to `.
 		`a scalar.`): next:
	elif type(Wdegree(a),integer) and type(Wdegree(b),integer) and
		Wdegree(a)<>Wdegree(b) then
		lprint(`Warning: Two sides of the `.i.`th assignment have `.
		`different degrees.`); next:
	fi:
	FormOrder_[a]:=evaln(FormOrder_[a]):
	WedgeDegree_[a]:=evaln(WedgeDegree_[a]):
fi:
if type(a,'name') then
	assign(a,b):
else
	j:=op(0,a):
	j(op(a)):=b:
	if j=`d` then
		a:=op(4,op(d)):
		a[b]:=0:
	fi:
fi:
			fi:
		else
		fi:
	od:
	FORGET(Wdegree,W_sub_,d_sub_,W_sub2_,ScalarForm_1_,df_,Wf_,Wdf_):
	for i in map(op,{indices(FormOrder_)}) do
		if not type(i,{'name','function'}) or Wdegree(i)=0 then
			FormOrder_[i]:=evaln(FormOrder_[i]):
		fi:
	od:
	for i in {indices(WedgeDegree_)} do op(i):
		if not type(%,{'name','function'}) or WedgeDegree_[op(i)]=0
		or WedgeDegree_[%]=0
		then
			WedgeDegree_[op(i)]:=evaln(WedgeDegree_[op(i)]):
		fi:
	od:
	RT:=proc(t)
		local i: for i in map(op,{indices(t)}) do lprint(t[i]); od:
	end:
	if nargs=0 then
	else	RT(t):
	fi:	NULL:
end:
CKW_:=proc() local n:
	global WedgeProduct_;
	for n in {args} do
		if type(n,{constant,`^`,name}) then
		elif type(n,function) then
			if op(0,n)=WedgeProduct_ then RETURN(true):
			else  RETURN(CKW_(op(n))): fi:
		else
			if member(true,map(CKW_,{op(n)})) then
				RETURN(true);
			fi:
		fi:
	od:	false
end:
Forder:=proc()
	local k,i,j,n:
	global ScalarForm_1_,W_sub2_, W_sub_, Wdegree, Wdf_, Wf_, d_sub_, df_,\
		FormOrder_, FormOrder_no;
	if nargs=0 then
		if FormOrder_[]=1 then
			FormOrder_[]:=0:
			lprint(`	The ordering is off.`):
		else
			FormOrder_[]:=1:
			lprint(`	The ordering is on, it is:`);
			k:=[]: n:=0:
			for i in map(op,{indices(FormOrder_)}) do
				for j from n by -1 to 1 do
					if FormOrder_[i]>FormOrder_[k[j]] then
						k:=[op(k[1..j]),i,op(k[j+1..n])]:
						n:=n+1:
						break:
					elif FormOrder_[i]=FormOrder_[k[j]] then
						ERROR(`The order for `.i.` and `.k[j].` are the same.`);
					fi:
				od:
				if j=0 then
					k:=[i,op(k)]:
					n:=n+1:
				fi
			od:
			print(k);
		fi:
	else
		if nops({args})<>nargs then
			ERROR(`Two or more arguments are equal.`): fi:
		k:=FormOrder_no-1:
        for i from nargs by -1 to 1 do
            if not type(args[i],{'name','function'}) or Wdegree(args[i]) = 0
                 then
                next
            fi;
            if CKW_(args[i]) or CKW_(EVAL(args[i])) then next fi;
            if FormOrder_[args[i]]<>0 and FormOrder_[args[i]]<k then
            	k:=FormOrder_[args[i]]: next:
            fi:
			for j in map(op,{indices(FormOrder_)}) do
				if FormOrder_[j]>=k then
					FormOrder_[j]:=FormOrder_[j]+1:
				fi:
			od:
            FormOrder_[args[i]]:=k: FormOrder_no:=FormOrder_no+1:
        od;
        for i in map(op,{indices(FormOrder_)}) do
			if not type(i,{'name','function'}) or Wdegree(i)=0 then
				FormOrder_[i]:=evaln(FormOrder_[i]):
			fi:
		od:
		lprint(`	The ordering of forms is:`):
		k:=[]: n:=0:
		for i in map(op,{indices(FormOrder_)}) do
			for j from n by -1 to 1 do
				if FormOrder_[i]>FormOrder_[k[j]] then
					k:=[op(k[1..j]),i,op(k[j+1..n])]:
					n:=n+1:
					break:
				elif FormOrder_[i]=FormOrder_[k[j]] then
					ERROR(`The order for `.i.` and `.k[j].` are the same.`);
				fi:
			od:
			if j=0 then
				k:=[i,op(k)]:
				n:=n+1:
			fi
		od:
		print(k);
	fi:	FORGET(Wdegree,W_sub_,d_sub_,W_sub2_,ScalarForm_1_,df_,Wf_,Wdf_):
end:
WedgeProduct:=proc(s)
	global ScalarForm_1_,W_sub2_, W_sub_, Wdegree, Wdf_, Wf_, d_sub_, df_,\
		WedgeProduct_;
	WedgeProduct_:=s:
 	s:=proc () W_sub_(args) end:
 	FORGET(Wdegree,W_sub_,d_sub_,W_sub2_,ScalarForm_1_,df_,Wf_,Wdf_):
	print(`The symbol for the wedge product is: `.WedgeProduct_); NULL:
end:
df_:=proc(x)
	options system,remember:
	if type(x,'function') and op(0,x)=`d` then RETURN(true) fi:
	false:
end:
Wf_:=proc(x)
	global WedgeProduct_;
	options system,remember:
	if type(x,'function') and op(0,x)=`WedgeProduct_` then RETURN(true) fi:
	false:
end:
Wdf_:=proc(x)
	global WedgeProduct_;
	options system,remember:
	if type(x,'function') and (op(0,x)=`d` or op(0,x)=`WedgeProduct_`) then
		RETURN(true):
	fi: false:
end:
Wdegree:=proc(e)
	global WedgeDegree_,WedgeProduct_;
	options system,remember:
	if type(e,array) then RETURN( map(Wdegree,e) ):
	elif type(e,'name') then RETURN(WedgeDegree_[e]):
	elif type(e,'function') then
		if WedgeDegree_[e]<>0 then RETURN(WedgeDegree_[e]):
		elif op(0,e)=`WedgeProduct_` then
			RETURN(  convert(map(Wdegree,[op(e)]),`+`)   ):
		elif op(0,e)=`d` then	RETURN(  Wdegree(op(e))+1  ):
		else RETURN(0):
		fi:
	elif type(e,{'constant',`^`}) then RETURN(0):
	elif type(e,`*`) then
		convert( map(Wdegree,[op(e)]), `+`):
	else
		if type(e,`+`) then
			map(Wdegree,{op(e)}):
			if nops('%')=1 then RETURN(op(%)):
			else RETURN('nonhmg') fi:
		fi:
		ERROR(`Wrong types of arguments`):
	fi:
end:
Simf:=proc(e)
	local c,f,i:
	if type(e,{equation,list,set,table,array}) then RETURN(map(Simf,args)): fi:
	c:=ScalarForm(e,f):
	if nargs>1 and type(args[2],'procedure') then
		c:=map(args[2],c,args[3..nargs]):
	else 	c:=map(normal,c):
	fi:
	SUM('c[i]*f[i]','i'=1..nops(f)):
end:
ScalarForm:=proc(e,f)
	local ee,ff:
	e: %: while '%' <> '%%' do % od:
	ff:=op(map(ScalarForm_1_,[op(frontend(indets,['%']))])):
	frontend( ScalarForm_2_, ['%%','%']):
	if nargs=1 then RETURN('%'):
	else
		if nops('%')>nops(['%%']) then f:=['%%',1] else f:=['%%']:fi:
		RETURN('%%'):
	fi:
end:
ScalarForm_1_:=proc(x)
	global WedgeDegree_;
	options system,remember:
	if WedgeDegree_[x]<>0 or Wdf_(x) then 'x' fi:
end:
ScalarForm_2_:=proc() local r,i,c:
	c:=[seq(diff(args[1],args[i]) , i=2..nargs)]:
	if has(c,{args[2..nargs]}) then
		lprint(`Warning: invalid term of the type: `.
		`form*form.`):
	fi:
	r:=subs({seq(args[i]=0 , i=2..nargs)},args[1]):
	if r=0 and nargs>1 then c else [op(c),r] fi:
end:
W_sub_:=proc()
	local i,j,a,e,x,eq,sg,k,l:
	global W_sub2_,W_sub_, WedgeDegree_, WedgeProduct_,FormOrder_, FormOrder_no;
	options system,remember:
	if nargs<2 then
		if nargs=0 then  RETURN(1):
		else             RETURN('args'):
  		fi:
  	elif type(args[1],array) then
  		RETURN(WedgeProduct_(W_sub_array_(args[1],args[2]),args[3..nargs]));
	else
		if member(0,{args}) then RETURN(0) fi:
		for i to nargs while ( WedgeDegree_[args[i]]<>0 or df_(args[i]) ) do
			if member(args[i],{args[1..i-1]})
				and type(Wdegree(args[i]),'odd') then RETURN(0);
			fi:
		od:
		if i=nargs+1 then
			if FormOrder_[]=0 then
				`WedgeProduct_`;RETURN('%(args)'): fi:
			0:
			for i to nargs do
				FormOrder_[args[i]]:
				if '%'=0 or '%'<'%%' then break fi:
			od:
			if i=nargs+1 then
				`WedgeProduct_`; RETURN('%(args)'):
			fi:
			eq:=[args]:
			sg:=1:
			k:=i-1:
			for i from i to nargs do
				if FormOrder_[eq[i]]=0 then
					next:
				fi:
				for j from k by -1 to 1 do
					if FormOrder_[eq[i]]>FormOrder_[eq[j]] then
						sg:=sg*((-1)^Wdegree(eq[i]))^(
							convert([seq(Wdegree(eq[l]),l=j+1..i-1)],`+`)):
						eq:=[op(eq[1..j]),eq[i],op(eq[j+1..i-1]),op(eq[i+1..nargs])]:
						k:=k+1:
						break:
					fi:
				od:
				if j=0 then
					sg:=sg*((-1)^Wdegree(eq[i]))^(
						convert([seq(Wdegree(eq[l]),l=1..i-1)],`+`)):
					eq:=[eq[i],op(eq[1..i-1]),op(eq[i+1..nargs])]:
					k:=k+1:
				fi:
			od:
			if k<nargs then
				for j from k+1 to nargs do
					FormOrder_[eq[j]]:=FormOrder_no:
					FormOrder_no:=FormOrder_no+1:
				od:
			fi:
		`WedgeProduct_`;op(eq):RETURN('%%(%)'*sg):
		fi;
		for j from i to nargs do a[j]:=ScalarForm(args[j],e[j]):od:
		SUM(  convert([seq(a[l][x[l]],l=i..nargs)],`*`)*
		      eq(args[1..i-1],seq(sg(e[l][x[l]]),l=i..nargs)),
		      seq(x[l]=1..nops(e[l]),l=i..nargs)   ):
		subs([sg=W_sub2_,eq=W_sub_],'%'): %:
	fi:
end:
W_sub2_:=proc(x)
	options system,remember;
	if Wf_(x) then op(x):
	elif x=1 then
	else 'x':fi:
end:
W_sub_array_:=proc(a1,a2)
	local d1,d2,m10,m11,m20,m21,n10,n11,n20,n21,a,i,j,k:
	global WedgeProduct_;
	options system,remember:
	d1:=[op(2,eval(a1))]: d2:=[op(2,eval(a2))]:
	if nops(d1)=1 then
		m10:=NULL;
		n10:=op(1,d1[1]);n11:=op(2,d1[1]);
	else
		m10:=op(1,d1[1]): m11:=op(2,d1[1]):
		n10:=op(1,d1[2]): n11:=op(2,d1[2]):
	fi:
	if nops(d2)=1 then
		m20:=NULL;
		n20:=op(1,d2[1]);n21:=op(2,d2[1]);
	else
		m20:=op(1,d2[1]): m21:=op(2,d2[1]):
		n20:=op(1,d2[2]): n21:=op(2,d2[2]):
	fi:
	if m10=NULL then
		if m20=NULL then
			if n10 <> n20 or n11 <> n21 then
				ERROR(`Dimensions do not match.`);
			fi:
			a:=convert([seq(WedgeProduct_(a1[i],a2[i]),i=n10..n11)],`+`):
		else
			if n10<>m20 or n11<>m21 then
				ERROR(`Dimensions do not match.`);
			fi:
			a:=array(n20..n21);
			for i from n20 to n21 do
				a[i]:=convert([seq(WedgeProduct_(a1[j],a2[j,i]),
					j=n10..n11)],`+`):
			od:
		fi:
	else
		if m20=NULL then
			if n10 <> n20 or n11 <> n21 then
				ERROR(`Dimensions do not match.`);
			fi:
			a:=array(m10..m11,1..1);
			for i from m10 to m11 do
				a[i,1]:=convert([seq(WedgeProduct_(a1[i,j],a2[j]),
					j=n10..n11)],`+`):
			od:
		else
			if n10<>m20 or n11<>m21 then
				ERROR(`Dimensions do not match.`);
			fi:
			a:=array(m10..m11,n20..n21);
			for i from m10 to m11 do
				for j from n20 to n21 do
					a[i,j]:=convert([seq(WedgeProduct_(a1[i,k],a2[k,j]),
						k=n10..n11)],`+`):
				od:
			od:
		fi:
	fi:
	RETURN(eval(a));
end:
d:=proc()
	d_sub_(args):
end:
d_sub_:=proc(e)
	local i,c,f,y,t:
	global WedgeProduct_;
	options system,remember:
	if type(e,array) then
		RETURN( map(d_sub_,e) ):
	elif type(e,'name') then
		RETURN('d(e)'):
	elif type(e,'constant') then
		RETURN( 0 ):
	fi:
	if type(e,`^`) then
		if type(op(2,e),'constant') then
	     		RETURN( op(2,e)*op(1,e)^(op(2,e)-1)*d(op(1,e)) ):
		fi:
		op(1,e):
		while type('%',`^`) do
			op(1,'%'):
		od:
		if '%'=-1 then RETURN(0): fi:
 		RETURN( 'd(e)' ):
	elif type(e,'function') then
		if op(0,e)=`d` then
		    	RETURN(0)
		elif op(0,e)=`WedgeProduct_` then t:=nops(e):
			RETURN( convert([seq(WedgeProduct_(op(1..i-1,e),d(op(i,e)),
			op(i+1..t,e))*(-1)^(convert(map(Wdegree,[op(1..i-1,e)]),
			`+`)),i=1..t)],`+`) ):
		fi:
		RETURN(mixpar( SUM('diff(e,t)*d(t)','t'=indets(e,name)) )):
	elif type(e,{`+`,`*`}) then
	    	c:=ScalarForm(e,f):
			if f=[1] then c:=op(c):
				RETURN( SUM('frontend(diff,[c,y])*d(y)',
		   		'y'=frontend(indets,[c])) ):
			fi:
			RETURN( SUM(  'SUM('frontend(diff,[c[i],y])*
			WedgeProduct_(d(y),f[i])',
		   	'y'=frontend(indets,[c[i]]) )+c[i]*d(f[i])',
	           	'i'=1..nops(f)) ):
	elif type(e,{'equation','list','set'}) then
		 RETURN(map(d,e)):
	fi:
	ERROR(`Case not covered`):
end:
L:=proc(ww,ee)
	local w,e,le,ne,j,i;
	global WedgeProduct_;
	w:=EVAL(ww): e:=EVAL(ee):
	if Wdegree(w)<>1 then ERROR(`the first argument must be an 1-form.`) fi;
	if type(e,{`+`,'list','equation','set'}) then
		map( proc(e,w) L(w,e) end, e,w):
	elif e = w then 1
	elif type(e,`*`) then
		if Wdegree(e)=0 then RETURN(0) fi:
		map(proc(e,w) if Wdegree(e)=0 then 'e' else L(w,e) fi end,e,w);
	elif type(e,function) then
		if op(0,e) = WedgeProduct_ then
			le:=[op(e)]: ne:=nops(le):
  			SUM('(-1)^( SUM('Wdegree(j)','j'=le[1..i-1]) )*
			    WedgeProduct_( op(le[1..i-1]),L(w,le[i]),op(le[i+1..ne]) )',
			   	'i'=1..ne):
		else 0
		fi:
	elif type(e,{`^`,constant}) then 0:
	elif not type(e,name) then ERROR(`case not covered`)
	else 0
	fi:
end:
mixpar:=proc(e)
	local f,x,os:
	if not hastype(e,function) then RETURN('e'): fi:
	if type(e,function) then
		if op(0,e)=diff then
			x:=op(2,e): f:=op(1,e):
			while  type(f,function) and op(0,f)=diff do
				x:=op(2,f),x: f:=op(1,f):
			od:
			os:=proc(a,b)
				if lexorder(a,b) then false else true fi
			end:
			RETURN(diff(f,op(sort([x],os)))):
		else	RETURN('e'):
		fi:
	else map(mixpar,e):
	fi:
end:
pick:=proc(e,w)
	local c,f,i;
    if type(e,{'list','equation','set'}) then RETURN(map(pick,args)) fi;
	if nargs<2 then ERROR(`Usage: pick(expr,form1,...)`);fi;
    c := ScalarForm(e,f);
    f := map(pick_1_,f,args[2 .. nargs]);
    convert([seq(c[i]*f[i] , i = 1 .. nops(f))],`+`)
end:
pick_1_:=proc(e,w)
    local sg,sd,l,ee,i,wi,j;
        if Wf_(e) then
			sg:=1:sd:=0; l:= nops(e); ee:=[op(e)];
			for i from 2 to nargs do
				wi:=Wdegree(args[i]);
				if wi = 0 or not type(args[i],{name,function}) then
					ERROR(`The `.i.`-th argument must be a form.`);
				elif member(args[i],ee,'j') then
					sg:=sg*(-1)^(wi*(Wdegree(WedgeProduct_(op(ee[j+1..l])))+sd));
					ee:=subsop(j=NULL,ee); sd:=sd+wi; l:=l-1;
				else RETURN(0);
				fi;
			od;
			RETURN(sg*WedgeProduct_(op(ee)));
		else
			for i from 2 to nargs do
				if Wdegree(args[i]) = 0 or not type(args[i],{name,function})
				then
					ERROR(`The `.i.`-th argument must be a form.`);
				elif e<>args[i] then RETURN(0);fi;
			od;
			RETURN(1);
		fi;
end:
WedgeProduct(`&^`):
EVAL  :=proc(x)  x: %:  while '%'<>'%%' do %:od: RETURN('%'): end:
FORGET:=proc()  local rtab,x,i:
        for i to nargs do rtab:=op(4,op(args[i])):
                if rtab<>NULL then
                        for x in {indices(rtab)} do
                                rtab[op(x)]:=evaln(rtab[op(x)]):od:
                fi:
        od;     NULL:   end:
SUM:=proc(expr) local rel,pt,i;
    rel := expr;
    for pt in args[2 .. nargs] do
        rel := convert([seq(subs(op(1,pt) = i,rel),i = op(2,pt))],`+`)
    od;
	%;
end:
print(`Saved to the file forms.m in the current directory`);
save `forms.m`;
print(`Make the library Forms.m in the current directory`);
Forms[Form]        :=op(Form):   	 	 Form:=evaln(Form):
Forms[Forder]      :=op(Forder):               Forder:=evaln(Forder):
Forms[WedgeProduct]:=op(WedgeProduct):   WedgeProduct:=evaln(WedgeProduct):
Forms[Wdegree]     :=op(Wdegree):             Wdegree:=evaln(Wdegree):
Forms[Simf]        :=op(Simf):                   Simf:=evaln(Simf):
Forms[ScalarForm]  :=op(ScalarForm):       ScalarForm:=evaln(ScalarForm):
Forms[d]           :=op(d):                         d:=evaln(d):
Forms[`&^`]        :=op(`&^`):                   `&^`:=evaln(`&^`):
Forms[L]	   :=op(L):		            L:='L':
Forms[mixpar]  	   :=op(mixpar):	       mixpar:='mixpar':
Forms[pick]        :=op(pick):                   pick:='pick':
save  `Forms.m`; quit
#@EOF##########################################################################

--PART-BOUNDARY=.19709151104.ZM11054.duke.edu--