Artifact d0b03229ff9a1331d4b97339611d14e63bed4d98495d1b3851c40f3fca01ca70:
- Executable file
r37/packages/mrvlimit/expon.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 13011) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/mrvlimit/expon.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 13011) [annotate] [blame] [check-ins using]
%--------------------------------------------------------------------------- % % these programs are written to sort the Taylor problem out; namely, the % problem of extracting the leading exponent together with its sign from % a taylor expression. % %---------------------------------------------------------------------------- module expon; load_package taylor; algebraic; expr procedure split(exp); begin scalar temp,current,ans; off mcd; ans:={}; if lisp(atom exp) then temp:={exp} else temp:=for k:=1:arglength(exp) collect part(exp,k); write "temp is ", temp; for k:=1:arglength(temp) do << current:=part(temp,k); if (lisp atom current) then << if(not freeof(current,ww)) then ans:=append(ans,{current}) else nil; >> else << if(not freeof(current,expt)) then ans:=append(ans,{current}) else nil; >>; >>; return ans; end; %load_package taylor; expr procedure collect_power(li); begin scalar ans; on rational; on exp; ans:=for k:=1:length(li) collect lpower(part(li,k),ww); return ans; end; %collect_power(split(1/2*(ww+x*ww^-1+2))); expr procedure conv(li); % converts our list of powers to exponents begin scalar ans,current; ans:={}; for k:=1:length(li) do << current:=part(li,k); %write "current is ", current; if(lisp atom current) then << if(not freeof(current,ww)) then ans:=append(ans,{1}) else nil; >> else << if(part(current,0)=expt) then ans:=append(ans,{part(current,2)}) else nil; >>; >>; return ans; end; %collect_power(split(1/2*(ww+x*ww^-1+2))); %conv(ws); %load_package assist; expr procedure find_expt(exp); begin scalar spli, coll, con, ans,ans2; %spli:=split(exp); %write "split is ", spli; coll:=collect_power(spli); write "collect is "; coll; con:=conv(coll); write "con is ", con; ans:=sortnumlist(con); write "ans is ", ans; ans2:=part(ans,1); return ans2; end; % we get something like % (expt(ww -1) %-------------------------------------------------------------------------- symbolic procedure find(u); begin off mcd; off factor; off exp; if(atom u) then << if(freeof(u,'ww)) then << if(numberp(u)) then return list('number,u) else << if(u='e) then return list('number,'e) else return list('x_exp,u) >>; >> else return list('expt,'ww,1) >> else << if(car u='expt) then return list('expt, cadr u, caddr u) else << if(car u='plus) then << if(atom cadr u and atom caddr u) then << if(length(cdr u)>2) then << if((cadr u='ww) and freeof(caddr u,'ww)) then return append({'expt,cadr u,1},find(append({'plus},cddr u))) else return append(find(cadr u), append({'plus},cddr u)) >> else << if(caddr u='ww) and freeof(cadr u,'ww) then return list('x_exp,cadr u,'expt,'ww,1) else << if(cadr u='ww and freeof(caddr u,'ww)) then return append({'expt,cadr u,1},find caddr u) else return append(find(cadr u),find(caddr u)) >> >> >> else << if(atom cadr u and pairp caddr u) then << if(length(cdr u)>2) then << if(cadr u='ww) then return append({'expt,'ww,1},find(append({'plus},cddr u))) else return append(find(cadr u),find(append({'plus},cddr u))) >> else return append(find(cadr u),find(caddr u)) >> else << if(pairp cadr u and pairp caddr u) then << if(length(cdr u)>2) then % plus has more than 2 args return append(find(cadr u),find(append({'plus},cddr u))) else return append(find(cadr u),find(caddr u)) >> else << if(pairp cadr u and atom caddr u) then << if(length(cdr u)>2) then %plus has more than two args << if(caddr u='ww) then << return find(cadr u).list('expt,'ww,1). find(append({'plus},cddr u)) >> else << if(numberp(caddr u)) then return find(cadr u).(list('number,caddr u).find(append({'plus},cdr u))) else nil >> >> else return append(find(cadr u),find(caddr u)) >> else return append(find(cadr u),{caddr u}) >> % else nil % unneccesary ? >> >> >> else << if(car u='lminus) then << if(numberp cadr u) then return list('number,u) else return find(cadr u) >> else << if(car u='quotient) then << if(numberp(cadr u) and numberp(caddr u)) then return list('number,cadr u, caddr u) else return append(find(cadr u), find(caddr u)) >> else << if(car u='minus) then << if(atom cdr u) then return find(cadr u) else << if(cadr u='expt and caddr u='ww) then return append(append({'minus},find(cadr u)),find(caddr u)) else return append({'minus},find(cadr u)) >> >> else << if(car u='times) then << if(atom cadr u and atom caddr u) then << if(not freeof(cadr u,'ww)) then return list('expt,cadr u,1) else << if(not freeof(caddr u,'ww)) then return list(nil,caddr u) else nil >> >> else << if(atom cadr u and pairp caddr u) then << if(not freeof(cadr u,'ww)) then return list('expt,cadr u,1) else return find(caddr u) >> else << if(pairp cadr u and pairp caddr u) then << if(length(cdr u))>2 then % times has +2 args return append(find(cadr u),find(append({'times},cddr u))) else return append(find(cadr u),find(caddr u)) >> else << if(pairp cadr u and atom caddr u) then << if(freeof(cadr u,'ww) and caddr u='ww) then return list('expt,'ww,1) else return append(find(cadr u),find(caddr u)) >> %else nil >> >> >> >> %else return find(cdr u) >> >> >> >> >> >>% ; end; algebraic; algebraic procedure fin(u); lisp ('list.find(u)); %---------------------------------------------------------------------------- % input to this procedure is a list % output is a list, containing all the exponents, marked expt, and any % numbers, flagged number % e.g. ww^-1+ww^-2 +2 % apply fin yields {expt,ww,-1,expt,ww,-2,number,2} % now apply find_numbers to this gives % {expt,-1,expt,-2,number,2} % the presence of the number means that there is a power of ww^0 present, % ie a constant term. If any of the expoents in the list are less than zero, % then we require the lowest one; if they are all positive, then zero is the % answer to be returned expr procedure find_numbers(li); begin scalar current,expt_list,ans,l,finished; % first, second, third; off mcd; on rational; on exp; li:=li; expt_list:={}; ans:={}; for k:=1:(length(li)-1) do << current:=part(li,k); if(current=expt) then << if(part(li,k+1)=ww) then expt_list:=append(expt_list,{expt,part(li,k+2)}); >> else << if(current=number) then expt_list:=append(expt_list,{number,part(li,k+1)}) else << if(current=x_expt) then expt_list:=append(expt_list,{x_part,part(li,k+1)}) else << if((current=lisp mk!*sq simp 'minus) and part(li,k+1)=expt) then expt_list:=append(expt_list,append({lisp reval 'minus},{expt,part(li,k+3)})); %else nil; >>; >>; >>; >>; % there is no x terms or numbers in the series exp return expt_list; end; %---------------------------------------------------------------------------- expr procedure find_least_expt(exp); begin scalar ans,find, current,result,expt_list,expt_list2,num_list, x_list; off mcd; % this causes a lot of problems when on, and some problems when off, % so I don't think I can win!!! expt_list:={}; num_list:={}; % initialisations x_list:={}; find:=fin(exp); ans:=find_numbers(find); if(lisp !*tracelimit) then write "exponent list is ", ans; %ans:=delete_all(-x,ans); if(freeof(ans,number)) then % there were no numbers in series exp, only % exponents << for k:=1:(arglength(ans)-1) do << if(part(ans,k)=lisp mk!*sq simp 'minus) then << if(numberp(part(ans,k+2)) and part(ans,k+2)<0) then expt_list:=append(expt_list,{lisp 'minus,part(ans,k+2)}); %else << % if(freeof(part(ans,k+2),x)) then % expt_list:=append(expt_list,{minus,part(ans,k+2)}); % >>; >> else << if((part(ans,k)=expt) and part(ans,k-1) neq (lisp mk!*sq simp 'minus)) then << if(numberp(part(ans,k+1))) then expt_list:=append(expt_list,{part(ans,k+1)}); >> else nil; >>; >>; %ans:=sortnumlist(ans); %result:=part(ans,1); %write "got up to here OK"; >> else << for k:=1:arglength(ans)-1 do << current:=part(ans,k); if((current=expt)) then % and part(ans,k+1)=lisp mk!*sq simp 'ww) then << if(freeof(part(ans,k+1),x)) then expt_list:=append(expt_list,{part(ans,k+1)}); >> else << if(current=number) then num_list:=append(num_list,{part(ans,k+1)}) else <<if((current=lisp mk!*sq simp 'minus) and numberp(part(ans,k+2)) and part(ans,k+2)<0) then expt_list:=append(expt_list,{lisp mk!*sq simp 'minus,part(ans,k+2)}) else nil; >>; >>; >>; >>; if(expt_list={}) then % we have only a number to deal with; ie power of % ww in series is 0 return append({number},num_list) else << if(num_list={}) then << if(freeof(expt_list,(lisp mk!*sq simp 'minus))) then << expt_list:=sortnumlist(expt_list); expt_list:={expt,part(expt_list,1)}; return expt_list; >> % else << % our list contains a power with a minus sign % want to find the least exponent, and then see if it is tagged % with a minus sign expt_list2:=expt_list; expt_list2:=delete_all((lisp mk!*sq simp 'minus),expt_list2); % list is now without minus expt_list2:=sortnumlist(expt_list2); expt_list2:=part(expt_list2,1); % smallest element, this is our expt % now want to check the sign of w with this exponent l:=0; finished:=0; while (l<=(arglength(expt_list)-1) and finished=0) do << if((part(expt_list,l)=(lisp mk!*sq simp 'minus)) and (part(expt_list,l+1)=expt_list2)) then << finished:=1; expt_list2:=append({lisp 'minus},{expt_list2}); >> else l:=l+1; >>; return expt_list2; >>; >> else << if(freeof(expt_list,lisp mk!*sq simp 'minus)) then << expt_list:=sortnumlist(expt_list); expt_list:={part(expt_list,1)}; % smallest element in the list if(part(expt_list,1)<0) then %%%%%% this is the value of e0 returned return append({expt},{part(expt_list,1)}) else return append({number},num_list); >> else << % doesn't matter what is in the number list, as minus is % present, meaning there is a negative exponent here expt_list:=delete_all(lisp mk!*sq simp 'minus, expt_list); expt_list:=sortnumlist(expt_list); return {lisp 'minus,part(expt_list,1)}; >>; >>; >>; end; endmodule; end;