share/simplification/facexp.mac

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