1%---------------------------------------------------------------------------- 2% 3% these programs are written to sort the Taylor problem out; namely, the 4% problem of extracting the leading exponent together with its sign from 5% a taylor expression. 6% 7%---------------------------------------------------------------------------- 8 9% Redistribution and use in source and binary forms, with or without 10% modification, are permitted provided that the following conditions are met: 11% 12% * Redistributions of source code must retain the relevant copyright 13% notice, this list of conditions and the following disclaimer. 14% * Redistributions in binary form must reproduce the above copyright 15% notice, this list of conditions and the following disclaimer in the 16% documentation and/or other materials provided with the distribution. 17% 18% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 19% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 20% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 21% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 22% CONTRIBUTORS 23% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 24% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 25% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 26% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 27% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 28% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 29% POSSIBILITY OF SUCH DAMAGE. 30% 31 32 33module expon; 34load_package taylor; 35algebraic; 36expr procedure mrv_split(exp); 37begin scalar temp,current,ans; 38 off mcd; 39 ans:={}; 40 if lisp(atom exp) then temp:={exp} 41 else temp:=for k:=1:arglength(exp) collect part(exp,k); 42 write "temp is ", temp; 43 for k:=1:arglength(temp) do 44 << current:=part(temp,k); 45 if (lisp atom current) 46 then << if(not freeof(current,ww)) 47 then ans:=append(ans,{current}) 48 else nil; 49 >> 50 else if (not freeof(current,expt)) 51 then ans:=append(ans,{current}) 52 else nil; 53 >>; 54 55 return ans; 56end; 57 58 59%load_package taylor; 60 61expr procedure collect_power(li); 62begin scalar ans; 63 on rational; on exp; 64 ans:=for k:=1:length(li) collect lpower(part(li,k),ww); 65 return ans; 66end; 67 68%collect_power(mrv_split(1/2*(ww+x*ww^-1+2))); 69 70expr procedure mrv_conv(li); % converts our list of powers to exponents 71begin scalar ans,current; 72 ans:={}; 73 for k:=1:length(li) do 74 << 75 current:=part(li,k); 76 %write "current is ", current; 77 if (lisp atom current) 78 then << if (not freeof(current,ww)) then ans:=append(ans,{1}) else nil; >> 79 else if (part(current,0)=expt) 80 then ans:=append(ans,{part(current,2)}) 81 else nil; 82 >>; 83 return ans; 84end; 85 86%collect_power(mrv_split(1/2*(ww+x*ww^-1+2))); 87 88%mrv_conv(ws); 89 90%load_package assist; 91 92expr procedure mrv_find_expt(exp); 93begin scalar spli, coll, con, ans,ans2; 94 %spli:=mrv_split(exp); %write "split is ", spli; 95 coll:=collect_power(spli); write "collect is "; coll; 96 con:=mrv_conv(coll); write "con is ", con; 97 ans:=sortnumlist(con); write "ans is ", ans; 98 ans2:=part(ans,1); 99 return ans2; 100end; 101 102% we get something like 103% (expt(ww -1) 104%-------------------------------------------------------------------------- 105symbolic procedure mrv_find(u); 106begin 107 off mcd; off factor; off exp; 108 if atom u then 109 << 110 if (freeof(u,'ww)) then 111 << if(numberp(u)) then return {'number,u} 112 else if (u='e) then return {'number,'e} 113 else return {'x_exp,u}; 114 >> 115 else return {'expt,'ww,1} 116 >> 117 else if (car u='expt) then return {'expt, cadr u, caddr u} 118 else if (car u='plus) 119 then << 120 if (atom cadr u and atom caddr u) 121 then << 122 if (length(cdr u)>2) 123 then << 124 if ((cadr u='ww) and freeof(caddr u,'ww)) 125 then return append({'expt,cadr u,1},mrv_find('plus . cddr u)) 126 else return append(mrv_find(cadr u), 'plus . cddr u) 127 >> 128 else if (caddr u='ww) and freeof(cadr u,'ww) 129 then return {'x_exp,cadr u,'expt,'ww,1} 130 else if (cadr u='ww and freeof(caddr u,'ww)) 131 then return append({'expt,cadr u,1},mrv_find caddr u) 132 else return append(mrv_find(cadr u),mrv_find(caddr u)) 133 >> 134 else if (atom cadr u and pairp caddr u) 135 then << if (length(cdr u)>2) 136 then << if (cadr u='ww) 137 then return append({'expt,'ww,1},mrv_find('plus . cddr u)) 138 else return append(mrv_find(cadr u),mrv_find('plus . cddr u)) 139 >> 140 else return append(mrv_find(cadr u),mrv_find(caddr u)) 141 >> 142 else if (pairp cadr u and pairp caddr u) 143 then << if (length(cdr u)>2) 144 % plus has more than 2 args 145 then return append(mrv_find(cadr u),mrv_find('plus . cddr u)) 146 else return append(mrv_find(cadr u),mrv_find(caddr u)) 147 >> 148 else if (pairp cadr u and atom caddr u) 149 then << if (length(cdr u)>2) 150 %plus has more than two args 151 then << if (caddr u='ww) 152 then return mrv_find(cadr u).{'expt,'ww,1}.mrv_find('plus . cddr u) 153 else if (numberp(caddr u)) 154 then return mrv_find(cadr u).({'number,caddr u}.mrv_find('plus . cdr u)) 155 else nil 156 >> 157 else return append(mrv_find(cadr u),mrv_find(caddr u)) 158 >> 159 else return append(mrv_find(cadr u),{caddr u}) 160 >> 161 else if (car u='lminus) 162 then << if (numberp cadr u) then return {'number,u} else return mrv_find(cadr u) >> 163 else if (car u='quotient) 164 then << if (numberp(cadr u) and numberp(caddr u)) 165 then return {'number,cadr u, caddr u} 166 else return append(mrv_find(cadr u), mrv_find(caddr u)) 167 >> 168 else if (car u='minus) 169 then << if (atom cdr u) then return mrv_find(cadr u) 170 else if (cadr u='expt and caddr u='ww) 171 then return append('minus . mrv_find(cadr u),mrv_find(caddr u)) 172 else return ('minus . mrv_find(cadr u)) 173 >> 174 else if (car u='times) 175 then << if (atom cadr u and atom caddr u) 176 then <<if (not freeof(cadr u,'ww)) 177 then return {'expt,cadr u,1} 178 else if (not freeof(caddr u,'ww)) 179 then return {nil,caddr u} 180 else nil 181 >> 182 else if (atom cadr u and pairp caddr u) 183 then <<if (not freeof(cadr u,'ww)) 184 then return {'expt,cadr u,1} 185 else return mrv_find(caddr u) 186 >> 187 else if (pairp cadr u and pairp caddr u) 188 then <<if (length(cdr u))>2 189 % times has +2 args 190 then return append(mrv_find(cadr u),mrv_find('times . cddr u)) 191 else return append(mrv_find(cadr u),mrv_find(caddr u)) 192 >> 193 else if (pairp cadr u and atom caddr u) 194 then <<if (freeof(cadr u,'ww) and caddr u='ww) 195 then return {'expt,'ww,1} 196 else return append(mrv_find(cadr u),mrv_find(caddr u)) 197 >> 198 %else nil 199 >> 200 %else return mrv_find(cdr u) 201end; 202 203algebraic; 204 205algebraic procedure mrv_fin(u); lisp ('list.mrv_find(u)); 206%---------------------------------------------------------------------------- 207% input to this procedure is a list 208% output is a list, containing all the exponents, marked expt, and any 209% numbers, flagged number 210% e.g. ww^-1+ww^-2 +2 211% apply fin yields {expt,ww,-1,expt,ww,-2,number,2} 212% now apply mrv_find_numbers to this gives 213% {expt,-1,expt,-2,number,2} 214% the presence of the number means that there is a power of ww^0 present, 215% ie a constant term. If any of the expoents in the list are less than zero, 216% then we require the lowest one; if they are all positive, then zero is the 217% answer to be returned 218 219expr procedure mrv_find_numbers(li); 220begin scalar current,expt_list,ans,l,finished; % first, second, third; 221 off mcd; on rational; on exp; 222 li:=li; 223 expt_list:={}; ans:={}; 224 225 for k:=1:(length(li)-1) do << 226 current:=part(li,k); 227 if(current=expt) then << 228 if (part(li,k+1)=ww) 229 then expt_list:=append(expt_list,{expt,part(li,k+2)}); 230 >> 231 else <<if (current=number) 232 then expt_list:=append(expt_list,{number,part(li,k+1)}) 233 else if (current=x_expt) 234 then expt_list:=append(expt_list,{x_part,part(li,k+1)}) 235 else if ((current=lisp mk!*sq simp 'minus) and part(li,k+1)=expt) 236 then expt_list:=append(expt_list,append({lisp reval 'minus},{expt,part(li,k+3)})); 237 %else nil; 238 >>; 239 >>; 240 % there is no x terms or numbers in the series exp 241 return expt_list; 242end; 243 244%---------------------------------------------------------------------------- 245expr procedure mrv_find_least_expt(exp); 246begin scalar ans,find, current,result,expt_list,expt_list2,num_list, x_list; 247 off mcd; % this causes a lot of problems when on, and some problems when off, 248 % so I don't think I can win!!! 249 250 expt_list:={}; 251 num_list:={}; % initialisations 252 x_list:={}; 253 254 find:=mrv_fin(exp); 255 ans:=mrv_find_numbers(find); 256 if(lisp !*tracelimit) then write "exponent list is ", ans; 257 %ans:=delete_all(-x,ans); 258 259 if(freeof(ans,number)) then % there were no numbers in series exp, only 260 % exponents 261 << 262 for k:=1:(arglength(ans)-1) do 263 <<if (part(ans,k)=lisp mk!*sq simp 'minus) then 264 <<if (numberp(part(ans,k+2)) and part(ans,k+2)<0) 265 then expt_list:=append(expt_list,{lisp 'minus,part(ans,k+2)}); 266 %else << 267 % if(freeof(part(ans,k+2),x)) then 268 % expt_list:=append(expt_list,{minus,part(ans,k+2)}); 269 % >>; 270 >> 271 else if ((part(ans,k)=expt) and part(ans,k-1) neq (lisp mk!*sq simp 'minus)) 272 then <<if (numberp(part(ans,k+1))) then expt_list:=append(expt_list,{part(ans,k+1)})>> 273 else nil; 274 >>; 275 %ans:=sortnumlist(ans); 276 %result:=part(ans,1); 277 %write "got up to here OK"; 278 >> 279 else << 280 for k:=1:arglength(ans)-1 do 281 <<current:=part(ans,k); 282 if ((current=expt)) then % and part(ans,k+1)=lisp mk!*sq simp 'ww) then 283 <<if (freeof(part(ans,k+1),x)) then expt_list:=append(expt_list,{part(ans,k+1)})>> 284 else if (current=number) then num_list:=append(num_list,{part(ans,k+1)}) 285 else if ((current=lisp mk!*sq simp 'minus) and numberp(part(ans,k+2)) and part(ans,k+2)<0) 286 then expt_list:=append(expt_list,{lisp mk!*sq simp 'minus,part(ans,k+2)}) 287 else nil; 288 >>; 289 >>; 290 if (expt_list={}) % we have only a number to deal with; ie power of 291 % ww in series is 0 292 then return append({number},num_list) 293 else if (num_list={}) 294 then <<if (freeof(expt_list,(lisp mk!*sq simp 'minus))) 295 then <<expt_list:=sortnumlist(expt_list); 296 expt_list:={expt,part(expt_list,1)}; 297 return expt_list; 298 >> 299 else << % our list contains a power with a minus sign 300 % want to find the least exponent, and then see if it is tagged 301 % with a minus sign 302 expt_list2:=expt_list; 303 expt_list2:=delete_all((lisp mk!*sq simp 'minus),expt_list2); 304 % list is now without minus 305 expt_list2:=sortnumlist(expt_list2); 306 expt_list2:=part(expt_list2,1); % smallest element, this is our expt 307 % now want to check the sign of w with this exponent 308 309 l:=0; finished:=0; 310 while (l<=(arglength(expt_list)-1) and finished=0) do 311 << if ((part(expt_list,l)=(lisp mk!*sq simp 'minus)) 312 and (part(expt_list,l+1)=expt_list2)) 313 then << finished:=1; 314 expt_list2:=append({lisp 'minus},{expt_list2}); 315 >> 316 else l:=l+1; 317 >>; 318 return expt_list2; 319 >>; 320 >> 321 else <<if (freeof(expt_list,lisp mk!*sq simp 'minus)) 322 then <<expt_list:=sortnumlist(expt_list); 323 expt_list:={part(expt_list,1)}; % smallest element in the list 324 if (part(expt_list,1)<0) 325 then %%%%%% this is the value of e0 returned 326 return append({expt},{part(expt_list,1)}) 327 else return append({number},num_list); 328 >> 329 else << % doesn't matter what is in the number list, as minus is 330 % present, meaning there is a negative exponent here 331 expt_list:=delete_all(lisp mk!*sq simp 'minus, expt_list); 332 expt_list:=sortnumlist(expt_list); 333 return {lisp 'minus,part(expt_list,1)}; 334 >>; 335 >>; 336end; 337 338endmodule; 339 340end; 341