Artifact e43f82357aefc85b899ae055bb9579b586ebb3b48b165ec67b543d32b0f80819:

• Executable file r38/packages/crack/crhomalg.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: 7363) [annotate] [blame] [check-ins using] [more...]

%********************************************************************
module homalgsys$
%********************************************************************
%  Routines for the efficient solution of homogeneous algebraic systems
%  Author: Thomas Wolf
%  May 2000

% extra global variables:
symbolic fluid '(tr_hom_alg)$

lisp(tr_hom_alg:=t)$

algebraic procedure bi_lin_expt(p)$
if p=1 then 0 else
if  arglength(p)<2 then 1 else
if (arglength(p)=2) and (part(p,0)=EXPT) then part(p,2)
                                         else write"error!"$

algebraic procedure find_hom_deg(p)$
% returns {dg1,dg2} where dg1 is the degree of     flin_ functions
%                         dg2 is the degree of non-flin_ functions
begin scalar tm,f$
 % take first term only
 if (arglength(p)<0) or (part(p,0) neq plus) then tm:=p 
                                             else tm:=part(p,1)$
 l1:=lisp gensym()$
 l2:=lisp gensym()$
 sb:=append(for each f in lisp(cons('LIST,              flin_ )) collect f=f*l1,
            for each f in lisp(cons('LIST,setdiff(ftem_,flin_))) collect f=f*l2)$
 tm:=sub(sb,tm)/tm$
 
 if freeof(tm,l1) then <<
  dg1:=0;
  dg2:=bi_lin_expt(tm);
 >>               else
 if freeof(tm,l2) then <<
  dg1:=bi_lin_expt(tm);
  dg2:=0;
 >>               else
 if (arglength(tm)=2) and (part(tm,0)=TIMES) then 
 if freeof(part(tm,1),l2) then <<dg1:=bi_lin_expt(part(tm,1));
                                 dg2:=bi_lin_expt(part(tm,2)) >>
                          else <<dg1:=bi_lin_expt(part(tm,2));
                                 dg2:=bi_lin_expt(part(tm,1)) >>
                                             else write"error3";
 % The next two lines are a test whether the equation is still
 % homogeneous. This is currently commented out to save time.
 % tm:=sub(sb,p)/l1**dg1/l2**dg2$
 % if tm neq p then write"INHOMOGENEOUS!"$

 return {dg1,dg2}
end$

symbolic operator make_hom_ansatz$
symbolic procedure make_hom_ansatz(f_1,f_2,d1,d2)$
begin scalar ans,ans1,ans2,h,fl,rply;

 if null f_1 then ans1:=1 else
 if null cdr f_1 then ans1:={'EXPT,car f_1,d1} else
 ans1:={'EXPT,cons('PLUS,f_1),d1}$

 if null f_2 then ans2:=1 else
 if null cdr f_2 then ans2:={'EXPT,car f_1,d2} else
 ans2:={'EXPT,cons('PLUS,f_2),d2}$

 ans:=reval {'TIMES,ans1,ans2}$

 return
 if (not pairp ans) or (car ans neq 'PLUS) then <<
  h:=gensym();
  {'LIST,{'TIMES,h,ans},{'LIST,h}}
 >>                                        else <<
  ans:=cdr ans$
  for each f in ans do <<
   h:=gensym()$
   fl:=cons(h,fl)$
   rply:=cons({'TIMES,h,f},rply)
  >>$
  {'LIST,cons('PLUS,rply),cons('LIST,fl)}
 >>
end$

symbolic procedure bi_lin_sep(p,fl_1,fl_2)$
% separation of p wrt all fl_1, fl_2
begin scalar fl_1,fl_2,f1,f2,f,su,pcp,sepli,cnd;

 for each f1 in fl_1 do <<
  if print_ then write "Separation wrt. ",f1,". "$
  su:=nil;
  for each f in fl_1 do 
  su:=cons(if f=f1 then {'EQUAL,f,1}
                   else {'EQUAL,f,0},su);
  su:=cons('LIST,su);
  pcp:=algebraic(sub(su,p));
  if print_ then <<write "Substitution done."$terpri()>>$

  for each f2 in fl_2 do <<
   su:=nil;
   for each f in fl_2 do 
   su:=cons(if f=f2 then {'EQUAL,f,1}
                    else {'EQUAL,f,0},su);
   su:=cons('LIST,su);
   cnd:=algebraic(sub(su,pcp));
   if pairp cnd and car cnd='MINUS then cnd:=cadr cnd;
   sepli:=union(list cnd,sepli);
   
  >>$
 >>$
 return cons('LIST,sepli)
end$

symbolic procedure bi_lin_eqn_lin_comb(pdes)$
% generates a linear combination of all pdes
% returns also list of unknown coefficients
begin scalar p,fl,rs,h$
 if print_ then <<
  write "Formulating a linear combination of all equations."$
  terpri()
 >>$
 for each p in pdes do
 <<h:=p$  % gensym()$
   fl:=cons(h,fl)$
   rs:=cons({'TIMES,h,get(p,'val)},rs)
 >>$
 rs:=cons('PLUS,rs)$
 return {'LIST,reval rs,cons('LIST,fl)}
end$

symbolic procedure drop_dep_bi_lin(arglist)$
begin scalar pdes,cnd,fl,f,cndcp,c,linde,again,fl_1,fl_2$
 pdes:=car arglist$
 again:=t$
 fl_1:=flin_;
 fl_2:=setdiff(ftem_,flin_);
 if pdes and cdr pdes then 
 repeat <<
  again:=nil$
  cnd:=bi_lin_eqn_lin_comb(pdes)$
  fl:=caddr cnd$
  cnd:=bi_lin_sep(cadr cnd,fl_1,fl_2)$ 
  if print_ then <<write"Now solving the linear system."$terpri()>>$
  !!arbint:=0;
  cnd:=cadr solveeval list(cnd,fl)$

  for f:=1:!!arbint do <<
   cndcp:=cnd;
   for c:=1:!!arbint do
   if c neq f then cndcp:=algebraic(sub(arbcomplex(c)=0,cndcp));
   cndcp:=cdr cndcp;

   while cndcp and 
         ((zerop caddar cndcp) or
          (not freeof(linde,cadar cndcp))) do cndcp:=cdr cndcp;
   if null cndcp then <<
    write"The computation to find redundant equations has to be done again."$
    terpri()
   >>            else linde:=cons(reval cadar cndcp,linde);
  >>$
  if null linde then write"No equations deleted."
                else write"Deleted redundant equations: ",linde$
  terpri()$
  for each f in linde do pdes:=drop_pde(f,pdes,nil)$
  if again then <<
   write"This computation has to be repeated"$
   terpri()
  >>
 >> until null again;

 return cons(pdes,cdr arglist)
end$

symbolic procedure find_factor_bi_lin(arglist)$
begin scalar ps,h,pdes,fc,rhs,lhs,lhsfl,cnd,cndcp,fl,fl_1,fl_2,
      hdg,dg1,dg2,indx1,indx2,again$

 write"Before starting to determine factorizable equations with "$terpri()$
 write"a given factor, all redundant equations have to be dropped."$terpri()$
 write"Has this already been done? (y/n) "$

 ps:=promptstring!*$
 promptstring!*:=""$ 
 repeat h:=termread() until (h='y) or (h='n)$
 If h='n then arglist:=drop_dep_bi_lin(arglist)$
 terpri()$

 fl_1:=flin_;
 fl_2:=setdiff(ftem_,flin_);
 pdes:=car arglist$
 write"Start of determining factorizable equations."$terpri()$ 
 repeat <<
  
  write"Give a factor of the equations to be found: (terminate with ; ) "$
  terpri()$
  fc:=termxread()$

  if null rhs then <<
   rhs:=bi_lin_eqn_lin_comb(pdes)$
   fl:=caddr rhs;  % alg mode list
   rhs:=cadr rhs;  % expression
  >>$

  hdg:=cdr algebraic find_hom_deg(fc)$
  dg1:=1-car  hdg$
  dg2:=1-cadr hdg$
  lhs:=make_hom_ansatz(fl_1,fl_2,dg1,dg2)$
  lhsfl:=caddr lhs;
  lhs:=cadr lhs$

  cnd:=bi_lin_sep(algebraic(fc*lhs-rhs),fl_1,fl_2)$ 
  if print_ then <<write"Now solving the linear system."$terpri()>>$
  !!arbint:=0;
  cnd:=cadr solveeval list(cnd,cons('LIST,append(cdr lhsfl,cdr fl)))$
  lhs:=algebraic(sub(cnd,lhs));

  for indx1:=1:!!arbint do <<
   cndcp:=lhs;
   for indx2:=1:!!arbint do
   if indx2 neq indx1 then cndcp:=algebraic(sub(arbcomplex(indx2)=0,cndcp))
                      else cndcp:=algebraic(sub(arbcomplex(indx2)=1,cndcp))$
   if not zerop cndcp then <<
    cndcp:={'TIMES,fc,cndcp}$

    pdes:=eqinsert(h:=mkeq(cndcp,ftem_,vl_,allflags_,t,list(0),nil,pdes),pdes)$
    if h and not freeof(pdes,h) then <<

% Drop one redundant equation

     write"New equation: ",h$ mathprint get(h,'val)$
     rhs:=nil;
    >>
   >>
  >>;
  write"Do you want to find further factorizable equations ",
       "with other factors? (y/n) "$
  repeat h:=termread() until (h='y) or (h='n)$
  If h='y then again:=t
          else again:=nil$
 >> until null again$

 promptstring!*:=ps$
 return cons(pdes,cdr arglist)
end$

endmodule$

end$