File r38/packages/int/csolve.red from the latest check-in

module csolve;   % routines to do with the C constants.

% Author: John P. Fitch.

fluid '(!*trint ccount cmap cmatrix cval loglist neweqn);

exports backsubst4cs,createcmap,findpivot,printvecsq, % printspreadc
   spreadc,subst4eliminateds;

imports nth,interr,!*multf,printsf,printsq,quotf,putv,negf,invsq,
   negsq,addsq,multsq,mksp,addf,domainp,pnth;

symbolic procedure findpivot cvec;
% Finds first non-zero element in CVEC and returns its cell number.
% If no such element exists, result is nil.
   begin         scalar i,x;
      i:=1;
      x:=getv(cvec,i);
      while i<ccount and null x do
      << i:=i+1;
         x:=getv(cvec,i) >>;
      if null x then return nil;
      return i
   end;

symbolic procedure subst4eliminatedcs(neweqn,substorder,ceqns);
% Substitutes into NEWEQN for all the C's that have been eliminated so
% far. These are given by CEQNS. SUBSTORDER gives the order of
% substitution as well as the constant multipliers. Result is the
% transformed NEWEQN.
   if null substorder then neweqn
   else begin    scalar nxt,row,cvar,temp;
      row:=car ceqns;
      nxt:=car substorder;
      if null (cvar:=getv(neweqn,nxt)) then
         return subst4eliminatedcs(neweqn,cdr substorder,cdr ceqns);
      nxt:=getv(row,nxt);
      for i:=0 : ccount do
      << temp:=!*multf(nxt,getv(neweqn,i));
         temp:=addf(temp,negf !*multf(cvar,getv(row,i)));
         putv(neweqn,i,temp) >>;
      return subst4eliminatedcs(neweqn,cdr substorder,cdr ceqns)
   end;


symbolic procedure backsubst4cs(cs2subst,cs2solve,cmatrix);
% Solves the C-eqns and sets vector CVAL to the C-constant values
% CMATRIX is a list of matrix rows for C-eqns after Gaussian
% elimination has been performed. CS2SOLVE is a list of the remaining
% C's to evaluate and CS2SUBST are the C's we have evaluated already.
   if null cmatrix then nil
   else begin    scalar eqnn,cvar,already,substlist,temp,temp2;
      eqnn:=car cmatrix;
      cvar:=car cs2solve;
      already:=nil ./ 1; % The S.Q. nil.
      substlist:=cs2subst;
% Now substitute for previously evaluated c's:
      while not null substlist do
      << temp:=car substlist;
         if not null getv(eqnn,temp) then
            already:=addsq(already,multsq(getv(eqnn,temp) ./ 1,
                                 getv(cval,temp)));
         substlist:=cdr substlist >>;
% Now solve for the c given by cvar (any remaining c's assumed zero).
      temp:=negsq addsq(getv(eqnn,0) ./ 1,already);
      if not null (temp2:=quotf(numr temp,getv(eqnn,cvar))) then
                                       temp:=temp2 ./ denr temp
      else temp:=multsq(temp,invsq(getv(eqnn,cvar) ./ 1));
      if not null numr temp then putv(cval,cvar,
                resimp rootextractsq subs2q temp);
      backsubst4cs(reversip(cvar . reversip cs2subst),
            cdr cs2solve,cdr cmatrix)
   end;

%**********************************************************************
% Routines to deal with linear equations for the constants C.
%**********************************************************************

symbolic procedure createcmap;
%Sets LOGLIST to list of things of form (LOG C-constant f), where f is
% function linear in one of the z-variables and C-constant is in S.F.
% When creating these C-constant names, the CMAP is also set up and
% returned as the result.
   begin         scalar i,l,c;
      l:=loglist;
      i:=1;
      while not null l do <<
         c:=(int!-gensym1('c) . i) . c;
         i:=i+1;
         rplacd(car l,((mksp(caar c,1) .* 1) .+ nil) . cdar l);
         l:=cdr l >>;
      if !*trint
	then printc ("Constants Created for log and tan terms:" . c);
      return c
   end;


symbolic procedure spreadc(eqnn,cvec1,w);
% Sets a vector 'cvec1' to coefficients of c<i> in eqnn.
    if domainp eqnn then putv(cvec1,0,addf(getv(cvec1,0),
                                !*multf(eqnn,w)))
    else begin    scalar mv,t1,t2;
        spreadc(red eqnn,cvec1,w);
        mv:=mvar eqnn;
        t1:=assoc(mv,cmap); %tests if it is a c var.
        if not null t1 then return <<
            t1:=cdr t1; %loc in vector for this c.
            if not (tdeg lt eqnn=1) then interr "Not linear in c eqn";
            t2:=addf(getv(cvec1,t1),!*multf(w,lc eqnn));
            putv(cvec1,t1,t2) >>;
        t1:=((lpow eqnn) .* 1) .+ nil; %this main var as sf.
        spreadc(lc eqnn,cvec1,!*multf(w,t1))
    end;

% symbolic procedure printspreadc cvec1;
%     begin
%         for i:=0 : ccount do <<
%            prin2 i;
%            printc ":";
%            printsf(getv(cvec1,i)) >>;
%         printc "End of printspreadc output"
%     end;

% symbolic procedure printvecsq cvec;
% % Print contents of cvec which contains s.q.'s (not s.f.'s).
% % Starts from cell 1 not 0 as above routine (printspreadc).
%    begin
%       for i:=1 : ccount do <<
%         prin2 i;
%         printc ":";
%         if null getv(cvec,i) then printc "0"
%         else printsq(getv(cvec,i)) >>;
%       printc "End of printvecsq output"
%    end;

endmodule;

end;