1/* -*-Macsyma-*- */
2eval_when(batch,ttyoff:true)$
3/*ASB;STOPEX 15
42:48pm  Wednesday, 4 November 1981
57:55pm  Saturday, 29 May 1982
6  Added a DIAGEVAL_VERSION for this file.
71:48pm  Saturday, 12 June 1982
8  Changed loadflags to getversions, DEFINE_VARIABLE:'MODE.
9*/
10
11eval_when(translate,
12	  transcompile:true,
13	  define_variable:'mode,
14	  modedeclare(function(freeofl),boolean))$
15
16put('stopex,15,'diageval_version)$
17/*
18eval_when([batch,loadfile],
19	  if get('gnauto,'diageval_version)=false
20	  then load('[gnauto,fasl,dsk,dgval]))$
21*/
22
23/* GNU Maxima */
24
25/* Commented out all local SPECIAL declarations.  For other changes,
26   search for `Maxima:' below. -wj */
27
28eval_when([batch,loadfile],
29  if get('gnauto,'diageval_version)=false
30  then load("genut"))$
31
32eval_when(translate,
33	  declare_translated(exwrt_power1,varmult,distribute,exwrt_power,
34			     freeofl,stopexpandl1,orpartitionl,ldelete,
35			     stopexpandl))$
36
37/* Switches  */
38define_variable(iforp,false,boolean)$
39define_variable(expandwrt_denom,false,boolean)$
40define_variable(expandwrt_nonrat,true,boolean)$
41
42stopexpand(exp,[varlist]):=
43  if atom(exp) or mapatom(exp)
44  then exp
45  else block([partswitch:true,inflag:true,piece],
46	     stopexpandl(exp,varlist))$
47
48expandwrt(exp,[varlist]):=
49  if atom(exp) or mapatom(exp)
50  then exp
51  else block([partswitch:true,inflag:true,piece],
52	     stopexpandl(exp,varlist))$
53
54expandwrtl(exp,varlist):=stopexpandl(exp,varlist)$
55
56stopexpandl(exp,varlist):=
57  if atom(exp) or mapatom(exp)
58  then exp
59  else block([inflag:true,partswitch:true,piece,ip0dum],
60	     if (ip0dum:inpart(exp,0))="+"
61	     then map(lambda([termdum],stopexpandl(termdum,varlist)),exp)
62	     else block(
63  [nonratdum,iforp:true,dendum],
64  if expandwrt_nonrat
65  then (nonratdum:
66	ldelete(varlist,last(orpartitionl(showratvars(exp),"[",varlist))),
67	for idum in nonratdum do
68	    if not atom(idum)
69	    then exp:subst(map(lambda([dum],stopexpandl(dum,varlist)),idum),
70			   idum,exp)),
71  if expandwrt_denom and (dendum:denom(exp))#1
72  then exp:num(exp)/stopexpandl(dendum,varlist),
73  stopexpandl1(exp,varlist)))$
74
75stopexpandl1(exp,varlist):=
76  if atom(exp) or mapatom(exp)
77  then exp
78  else block([ip0dum:inpart(exp,0),dum:1,varfound:false],
79  modedeclare(varfound,boolean),
80	     if freeofl(varlist,exp)
81	     then exp
82	     else if freeof("+",exp) then return(exp),
83	     if ip0dum="+"
84	     then return(map(lambda([termdum],
85				    stopexpandl1(termdum,varlist)),exp))
86	     else if ip0dum="^"
87		  then if inpart(exp,1,0)="+"
88	               then exwrt_power(exp,varlist)
89		       else exp
90		  else if ip0dum="*"
91		       then (for idum in exp do
92			         if not freeofl(varlist,idum)
93			         then (idum:stopexpandl1(idum,varlist),
94				       if varfound
95				       then dum:distribute(dum,idum,varlist)
96				       else (varfound:true,
97					     dum:varmult(dum,idum,varlist)))
98			         else if varfound
99				      then dum:varmult(idum,dum,varlist)
100				      else dum:dum*idum,
101		             dum)
102		       else if matrixp(exp) or listp(exp)
103		            then matrixmap(lambda([dumm],
104						  stopexpandl1(dumm,varlist)),
105				           exp)
106		            else if ip0dum="." and expandwrt_nonrat
107		                 then remove_nested_dots0l(map(lambda([dum],
108							      stopexpandl1(dum,
109								     varlist)),
110							       exp),
111							   varlist)
112			         else exp)$
113
114exwrt_power(exp,varlist):=block(
115  [ip1dum,ip2dum1,exwrtlist,splitdum,fsplitdum],
116  /* DECLARE(EXWRTLIST,SPECIAL), */
117  if inpart(exp,0)#"^" then return(exp),
118  if not freeofl(varlist,ip1dum:inpart(exp,1))
119     and integerp(ip2dum1:inpart(exp,2))
120     and (mode_identity(fixnum,ip2dum1))>1
121     and inpart(ip1dum,0)="+"
122  then (splitdum:orpartitionl(ip1dum,"+",varlist),
123	if (fsplitdum:first(splitdum))#0
124	then (exwrtlist:cons(1,exwrt_power1(last(splitdum),ip2dum1,varlist)),
125	      sum(varmult(fsplitdum^kdum*ip2dum1!/(kdum!*(ip2dum1-kdum)!),
126			  first(exwrtlist:rest(exwrtlist)),
127			  varlist),
128                         /* Maxima: added MODE_IDENTITY for translator */
129		  kdum,0,mode_identity(fixnum,ip2dum1)))
130	else first(exwrt_power1(last(splitdum),ip2dum1,varlist)))
131  else exp)$
132
133exwrt_power1(exp,powerdum,varlist):=(
134  modedeclare(powerdum,fixnum),
135 block(
136  [dum:[exp,1],firstdum:stopexpandl1(exp,varlist)],
137  if powerdum=1 then return(dum),
138  if inpart(exp,0)#"+"
139  then for idum:2 thru powerdum do
140	   dum:cons(exp^idum,dum)
141  else for idum:2 thru powerdum do
142	   dum:cons(firstdum:
143		    map(lambda([dum],multthru(dum,firstdum)),exp),dum),
144  dum))$
145
146varmult(fact,exp,varlist):=block(
147  [splitdum:orpartitionl(exp,"+",varlist)],
148  fact*first(splitdum)+multthru(fact,last(splitdum)))$
149
150distribute(exp1,exp2,varlist):=block(
151  [splitexp1:orpartitionl(exp1,"+",varlist),
152   splitexp2:orpartitionl(exp2,"+",varlist),
153   fsplexp1,fsplexp2,lsplexp1,lsplexp2],
154  lsplexp1:last(splitexp1),
155  lsplexp2:last(splitexp2),
156  (fsplexp1:first(splitexp1))*(fsplexp2:first(splitexp2))
157  +(if fsplexp1#0
158    then varmult(fsplexp1,stopexpandl1(lsplexp2,varlist),varlist)
159    else 0)
160  +(if fsplexp2#0
161    then varmult(fsplexp2,stopexpandl1(lsplexp1,varlist),varlist)
162    else 0)
163  +(if inpart(lsplexp1,0)="+"
164    then map(lambda([term],stopexpandl1(term*lsplexp2,varlist)),lsplexp1)
165    else if inpart(lsplexp2,0)="+"
166	 then map(lambda([term],stopexpandl1(term*lsplexp1,varlist)),lsplexp2)
167	 else lsplexp1*lsplexp2))$
168
169expandwrt_factored(exp,[varlist]):=
170  if listp(exp) or matrixp(exp)
171  then matrixmap(lambda([dum],apply('expandwrt_factored,cons(dum,varlist))),
172		 exp)
173  else block([iforp:true,piece,partswitch:true,inflag:true,dum],
174	     dum:orpartitionl(exp,"*",varlist),
175	     first(dum)*stopexpandl(last(dum),varlist))$
176
177eval_when(batch,ttyoff:false)$
178