1/* -*- MACSYMA -*- */
2eval_when(batch,ttyoff:true)$
3/*ASB;FACEXP 10
42:37pm  Wednesday, 4 November 1981
5*/
6/* commented out of DOE MACSYMA
7EVAL_WHEN(TRANSLATE,
8          IF GRAPHLOAD#TRUE THEN LOAD(['GRAPH,'FASL,'DSK,'DGVAL]),
9          TRANSCOMPILE:TRUE,
10          DECLARE([ARGDUM,ARGDUM2,OP_FCN_LIST],SPECIAL),
11          MODEDECLARE(FUNCTION(NULLLISTP,FREEOFL),BOOLEAN))$
12
13EVAL_WHEN(BATCH,IF DEBUGLOAD#TRUE THEN LOAD(['DEBUG,'FASL,'DSK,'DGVAL]))$
14*/
15
16/* Autoloads */
17/* commented out of DOE MACSYMA
18IF STATUS(FEATURE,ITS)=TRUE
19THEN (SETUP_AUTOLOAD([GENUT,FASL,DSK,DGVAL],
20                     OPMAP,INTERSECT_LIST,LDELETE,IFLOPMAP,FREEOFL,ORPARTITION,
21                     LISTOFOPS,ORPARTITIONL,RLOIEWL,IFNOTCONS,LISTTOSUM,
22                     LISTOFOPS_NONRAT,NULLLISTP,SETLIST),
23      SETUP_AUTOLOAD([INDEX,FASL,DSK,DGVAL],INDEX))$
24
25IF STATUS(FEATURE,MULTICS)=TRUE
26THEN (SETUP_AUTOLOAD(">udd>Mathlab>Brenner>genut",
27                     OPMAP,INTERSECT_LIST,LDELETE,IFLOPMAP,FREEOFL,ORPARTITION,
28                     LISTOFOPS,ORPARTITIONL,RLOIEWL,IFNOTCONS,
29                     LISTOFOPS_NONRAT,NULLLISTP,SETLIST),
30      SETUP_AUTOLOAD(">udd>Mathlab>Brenner>index",INDEX))$
31*/
32
33/*ASB;FACEX1 1
344:19pm  Monday, 7 February 1983
35  Split off from FACEXP 15
36*/
37/*
38EVAL_WHEN([TRANSLATE,BATCH,LOADFILE],
39          IF GET('FACEXP,'VERSION)=FALSE
40          THEN (LOAD('[FACEXP,FASL]),
41                LOAD('[GNDECL,FASL])))$
42*/
43
44/* GNU Maxima */
45
46/* Commented out all local SPECIAL declarations. -wj */
47
48eval_when([batch,loadfile],
49  if get('gnauto,'diageval_version)=false
50  then load("genut"))$
51
52eval_when(translate,
53          declare_translated(operator0p,multthrusplit,lopplusp,orpartition,
54                             collectterms0,collecttermsl,ldelete,listtosum,
55                             ifmultthru,intersect_list,zerosubst,facexpform1,
56                             opmap,orpartitionl,facexpform,freeofl,nextlayer,
57                             iflopmap,facexpl,argsplit,setlist,facsuml,
58                             nulllistp,autoform),
59          mode_declare(function(nulllistp,freeofl),boolean))$
60
61/* Variable definitions */
62
63define_variable(nextlayerfactor,false,boolean)$
64define_variable(facsum_combine,true,boolean)$
65
66/* Predicates */
67
68lopplusp(exp):=is(inpart(exp,0)="+")$
69
70operator0p(exp):=block(
71  [ip0dum],
72  is((ip0dum:inpart(exp,0))='operator or ip0dum=nounify('operator)))$
73
74orderlastp(exp1,exp2):=orderlessp(last(exp1),last(exp2))$
75
76/* User accessible functions */
77/*
78FACTENEXPAND(EXP,[ARGDUMLIST]):=BLOCK(
79  [INDEXEXPAND_CANONICAL:FALSE,PARTSWITCH:TRUE,INFLAG:TRUE,PIECE],
80  FACEXPTENL(CONS(INDEXEXPAND(EXP),ARGDUMLIST)))$
81*/
82factorfacsum(exp,[argdum]):=block(
83  [expdum,ip0dum,partswitch:true,inflag:true,piece],
84  if atom(exp) then return(exp),
85  if nulllistp(argdum) then return(autoform(exp)),
86  if matrixp(exp) or listp(exp) or inpart(exp,0)="="
87  then return(map(lambda([elemdum],apply('factorfacsum,cons(elemdum,argdum))),
88                  exp)),
89  expdum:autoform(exp),
90  if (ip0dum:inpart(expdum,0))="^" or ip0dum="*"
91  then map(lambda([factdum],apply('factorfacsum,cons(factdum,argdum))),
92           expdum)
93  else facsuml(cons(expdum,argdum)))$
94
95
96facsum([arglist]):=facsuml(arglist)$
97
98facsuml(arglist):=block(
99  [factorflag:false,partswitch:true,inflag:true,piece,
100   farglist:first(arglist)],
101  if matrixp(farglist)
102  then matrixmap(lambda([dum],facsuml(cons(dum,rest(arglist)))),farglist)
103  else if listp(farglist) or inpart(farglist,0)="="
104       then map(lambda([dum],facsuml(cons(dum,rest(arglist)))),farglist)
105       else block(
106  [argdum:rest(arglist),argdum2:[],ratfac:false],
107  /* DECLARE([ARGDUM,ARGDUM2],SPECIAL), */
108  setlist(argsplit(farglist,argdum),'argdum,'argdum2),
109  facexpl(cons(apply('ratsimp,cons(farglist,argdum:ratsimp(argdum))),
110               argdum))))$
111
112/* Functions mainly for internal use */
113
114facexp([arglist]):=facexpl(arglist)$
115
116nextlayer(exp):=block(
117  if not nulllistp(argdum2)
118  then iflopmap("*",
119                lambda([dum],facsuml(cons(dum,argdum2))),
120                if nextlayerfactor
121                then autoform(exp)
122                else exp)
123  else autoform(exp))$
124
125facexpl(arglist):=block(
126  [expdum:first(arglist),partitiondum,
127   argdum:rest(arglist),nextlayerfactorsave:nextlayerfactor,
128   nextlayerfactor:false,numexpdum,denexpdum],
129  modedeclare(nextlayerfactorsave,boolean),
130  /* DECLARE([NUMEXPDUM,DENEXPDUM],SPECIAL), */
131  if member('nextlayerfactor,arglist)
132  then (arglist:delete('nextlayerfactor,arglist),
133        nextlayerfactor:true)
134  else nextlayerfactor:nextlayerfactorsave,
135  if nulllistp(argdum) or length(arglist)<2 or freeofl(rest(arglist),expdum)
136  then return(nextlayer(expdum)),
137  numexpdum:facexpform(num(expdum)),
138  if (denexpdum:denom(expdum))#1 then denexpdum:facexpform(denom(expdum)),
139  if inpart(numexpdum,0)="+"
140     and not freeofl(argdum,numexpdum)
141     and not facsum_combine
142  then if denexpdum#1
143       then (partitiondum:orpartitionl(numexpdum,"+",argdum),
144             multthru(denexpdum^-1,last(partitiondum))+
145                denexpdum^-1*first(partitiondum))
146       else numexpdum
147  else numexpdum*denexpdum^-1)$
148
149facexpform(exp):=(
150  exp:opmap(exp,["+",'vplus,"*",'vstar]),
151  if inpart(exp,0)="+"
152  then facexpform1(exp)
153  else exp)$
154
155facexpform1(expdum):=block(
156  [subdum:zerosubst(argdum,expdum)],
157  /* DECLARE(SUBDUM,SPECIAL), */
158  expdum-subdum+nextlayer(subdum))$
159/*
160FACEXPTEN([ARGLIST]):=FACEXPTENL(ARGLIST)$
161
162FACTORFACEXPTEN(EXP,[ARGLIST]):=
163  IFLOPMAP("*",
164           LAMBDA([FACDUM],FACEXPTENL(CONS(FACDUM,ARGLIST))),
165           AUTOFORM(EXP))$
166
167FACEXPTENL(ARGLIST):=BLOCK(
168  [FACEXPTENFLAG:TRUE],
169  /* DECLARE(FACEXPTENFLAG,SPECIAL), */
170  modedeclare(facexptenflag,boolean),
171  facsuml(append(arglist,listoftens(first(arglist)))))$
172*/
173
174vplus(exp):=block(
175  [vpsdum:map(lambda([term],
176                     if nulllistp(intersect_list(showratvars(term),argdum))
177                     then nextlayer(term)
178                     else opmap(term,op_fcn_list)),
179              exp)],
180  /* DECLARE([OP_FCN_LIST,VPSDUM],SPECIAL), */
181  if inpart(vpsdum,0)="+"
182  then facexpform1(vpsdum)
183  else vpsdum)$
184
185vstar(exp):=block(
186  [argsexpdum:args(exp),partitiondum,expiargdum],
187  for iargdum in argsexpdum do
188    if inpart(iargdum,0)="+"
189    then if not nulllistp(intersect_list(argdum,showratvars(iargdum)))
190         then (partitiondum:orpartitionl(facexpform(iargdum),"+",argdum),
191               expiargdum:ifmultthru(1/iargdum,exp),
192               exp:ifmultthru(expiargdum,last(partitiondum))+
193                        expiargdum*nextlayer(first(partitiondum)))
194         else exp:exp/iargdum*nextlayer(iargdum),
195  if inpart(exp,0)="+"
196  then facexpform1(exp)
197  else exp)$
198
199fplus(exp):=block(
200  /* DECLARE([LIST,OP_FCN_LIST],SPECIAL), */
201  iflopmap("+",
202           lambda([dum],opmap(dum,op_fcn_list)),
203           listtosum(ldelete(list,args(exp)))))$
204
205fexpt(exp):=block(
206  [ip1exp:zerosubst(list,inpart(exp,1))],
207  /* DECLARE([LIST,IP1EXP],SPECIAL), */
208  if ip1exp=0
209  then 0
210  else ip1exp^zerosubst(list,inpart(exp,2)))$
211
212fstar(exp):=block(
213  /* DECLARE(LIST,SPECIAL), */
214  if ldelete(list,args(exp))=args(exp)
215  then map(lambda([dum],opmap(dum,op_fcn_list)),exp)
216  else 0)$
217
218zerosubst(list,exp):=
219  if member(exp,list)
220  then 0
221  else opmap(exp,["*",'fstar,"+",'fplus,"^",'fexpt])$
222
223ifmultthru(exp1,exp2):=
224  if inpart(exp2,0)="+"
225  then multthru(exp1,exp2)
226  else exp1*exp2$
227
228
229/* COLLECTTEN(EXP):=COLLECTTERMSL(EXP,LISTOFTENS(EXP))$ */
230
231collectterms(exp,[varlist]):=collecttermsl(exp,varlist)$
232
233collecttermsl(exp,varlist):=block(
234  [partswitch:true,inflag:true,piece],
235  apply('collectterms0,cons(exp,argsplit(exp,varlist))))$
236
237collectterms0(exp,thisleveldum,nextleveldum):=block(
238  [iforp:true,splitdum1,splitdum2,splitdum3,anslist:[],
239   prevdum,lsplit3,ansdum,lastdumsave,prevlastdum,
240   rthisleveldum,fthisleveldum],
241  modedeclare(lsplit3,fixnum),
242  /* DECLARE([SPLITDUM3,ANSDUM],SPECIAL), */
243  if exp=0 then return(0),
244  if nulllistp(thisleveldum) or freeofl(thisleveldum,exp)
245  then if nulllistp(nextleveldum)
246       then return(exp)
247       else (splitdum1:orpartitionl(exp,"+",nextleveldum),
248             return(collecttermsl(first(splitdum1),nextleveldum)
249                    +iflopmap("+",
250                              lambda([termdum],
251                                     collecttermsl(termdum,nextleveldum)),
252                              last(splitdum1)))),
253  rthisleveldum:rest(thisleveldum),
254  if freeof(fthisleveldum:first(thisleveldum),exp)
255  then return(collectterms0(exp,rthisleveldum,nextleveldum)),
256  splitdum1:orpartitionl(exp,"+",thisleveldum),
257  splitdum2:orpartition(last(splitdum1),"+",fthisleveldum),
258  ansdum:collecttermsl(first(splitdum1),nextleveldum)
259         +collectterms0(first(splitdum2),rthisleveldum,nextleveldum),
260  if not lopplusp(splitdum3:last(splitdum2))
261  then return(ansdum+collecttermsl(splitdum3,nextleveldum)),
262  splitdum3:sort(maplist(lambda([term],orpartition(term,"*",fthisleveldum)),
263                 splitdum3),
264                 'orderlastp),
265  lsplit3:length(splitdum3)-1,
266  prevlastdum:inpart(splitdum3,1,2),
267  prevdum:inpart(splitdum3,1,1),
268  splitdum3:rest(splitdum3),
269  for idum thru lsplit3 do
270       (if (lastdumsave:inpart(splitdum3,idum,2))=prevlastdum
271        then prevdum:prevdum+inpart(splitdum3,idum,1)
272        else (anslist:endcons([prevdum,prevlastdum],anslist),
273              prevdum:inpart(splitdum3,idum,1),
274              prevlastdum:lastdumsave),
275        if idum=lsplit3
276        then anslist:endcons([prevdum,prevlastdum],anslist)),
277  listtosum(maplist(lambda([dum],
278                           if freeofl(rthisleveldum,first(dum))
279                           then collecttermsl(first(dum),nextleveldum)
280                                    *last(dum)
281                           else multthrusplit(last(dum),
282                                              collectterms0(first(dum),rthisleveldum,
283                                                            nextleveldum),
284                                              rthisleveldum)),
285                    anslist))+ansdum)$
286
287argsplit(exp,list):=block(
288  [listargsdum:[],newlist:[]],
289  for arg in list do
290      if listp(arg)
291      then listargsdum:append(listargsdum,arg)
292      else if operator0p(arg)
293           then newlist:
294                append(newlist,apply('listofops_nonrat,cons(exp,args(arg))))
295           else newlist:cons(arg,newlist),
296  [newlist,listargsdum])$
297
298multthrusplit(factordum,sumdum,rthisleveldum):=block(
299  [splitdum1:orpartitionl(sumdum,"+",rthisleveldum)],
300  multthru(factordum,last(splitdum1))+factordum*first(splitdum1))$
301
302autoform(exp) :=
303  block([fun:get('facsum,'automatic)],
304    if not member('noun,apply('properties,[fun]))
305    then apply(fun,[exp])
306    else apply(nounify(fun),[exp]))$
307
308if get('facsum,'automatic)=false
309then put('facsum,'nonumfactor,'automatic)$
310
311sqfrfacsum([arglist]):=block(
312  [dum,autodum:get('facsum,'automatic)],
313  /* DECLARE([AUTODUM,DUM],SPECIAL), */
314  put('facsum,'sqfr,'automatic),
315  dum:facsuml(arglist),
316  put('facsum,autodum,'automatic),
317  dum)$
318
319eval_when(batch,ttyoff:false)$
320