File r38/packages/mrvlimit/expon.red artifact d0b03229ff part of check-in f2fda60abd


 %---------------------------------------------------------------------------
%
% 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;


REDUCE Historical
REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]