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