reduce-historical/r38/packages/algint/zmodule.red

module zmodule;

% Author: James H. Davenport.

fluid '(!*galois
        !*tra
        !*trfield
        !*trmin
        basic!-listofallsqrts
        basic!-listofnewsqrts
        commonden
        gaussiani
        listofallsqrts
        listofnewsqrts
        sqrt!-places!-alist
        taylorasslist);

exports zmodule;
imports !*multf,sqrtsinsql,sortsqrts,simp,!*q2f,actualsimpsqrt,printsf;
imports prepf,substitutesq,printsq,mapply,!*multsq,mkilist;
imports mkvecf2q,mkvec,mkidenm,invsq,multsq,negsq,addsq,gcdn;
imports !*invsq,prepsq;

symbolic procedure zmodule(w);
begin
  scalar reslist,denlist,u,commonden,basis,p1,p2,hcf;
  % w is a list of elements (place.residue)=sq.
  for each v in w do <<
    u:=cdr v;
    reslist:=u.reslist;
    denlist:=(denr u).denlist >>;
  basis:=sqrtsinsql(reslist,nil);
  if null u or null cdr u or !*galois
    then go to nochange;
  reslist:=check!-sqrts!-dependence(reslist,basis);
  denlist:=for each u in reslist
             collect denr u;
nochange:
 commonden:=mapply(function(lambda u,v;
                      multf(u,quotf(v,gcdf(u,v)))),denlist)./1;
  u:=nil;
  for each v in reslist do
    u:=(numr !*multsq(v,commonden)).u;
  reslist:=u;
    % We have effectively reversed RESLIST twice, so it is in correct
    % order.
  u:=bexprn(reslist);
  basis:=car u;
  reslist:=cdr u;
  denlist:=nil;
  while basis do <<
    p1:=reslist;
    p2:=w;
    u:=nil;
    hcf:=0;
    while p1 do <<
      if 0 neq caar p1
        then  <<
          u:=((caar p2).(caar p1)).u;
          hcf:=gcdn(hcf,caar p1) >>;
      p1:=cdr p1;
      p2:=cdr p2 >>;
    if hcf neq 1
     then u:=for each uu in u collect
               (car uu). ( (cdr uu) / hcf);
    denlist:=(prepsq !*multsq(car basis,
                              multsq(!*f2q hcf,!*invsq commonden))
                  .u).denlist;
    basis:=cdr basis;
    reslist := for each j in reslist collect cdr j>>;
  return denlist
  end;


symbolic procedure bexprn(wlist);
begin
  scalar basis,replist,w,w2,w3,p1,p2;
  % wlist is a list of sf.
  w:=reverse wlist;
  replist:=nil;
  while w do <<
    w2:=sf2df car w;
    % now ensure that all elements of w2 are in the basis.
    w3:=w2;
    while w3 do <<
      % caar is the sf,cdar a its coefficient.
      if not member(caar w3,basis)
        then <<
          basis:=(caar w3).basis;
          replist:=mapcons(replist,0) >>;
          % adds car w3 to basis.
      w3:=cdr w3 >>;
    replist:=mkilist(basis,0).replist;
    % builds a new zero representation.
    w3:=w2;
    while w3 do <<
      p1:=basis;
      p2:=car replist;
      %the list for this element.
      while p1 do <<
        if caar w3 = car p1
          then rplaca(p2,cdar w3);
        p1:=cdr p1;
        p2:=cdr p2 >>;
      w3:=cdr w3 >>;
    w:=cdr w >>;
  return mkbasis(basis,replist)
  end;


symbolic procedure mkbasis(basis,reslist);
begin
  scalar row,nbasis,nreslist,u,v;
  basis:=for each u in basis collect !*f2q u;
  % basis is a list of sq's
  % reslist is a list of representations in the form
  % ( (coeff1 coeff2 ...)    ...).
  nreslist:=mkilist(reslist,nil);
  % initialise our list-of-lists.
  trynewloop:
  row := mapovercar reslist;
  reslist := for each j in reslist collect cdr j;
  if obvindep(row,nreslist)
    then u:=nil
    else u:=lindep(row,nreslist);
  if u
    then <<
      % u contains the numbers with which to add this new item into the
      % basis.
      v:=nil;
      while nbasis do <<
        v:=addsq(car nbasis,!*multsq(car basis,car u)).v;
        nbasis:=cdr nbasis;
        u:=cdr u >>;
      nbasis:=reversip v >>
    else <<
      nreslist:=pair(row,nreslist);
      nbasis:=(car basis).nbasis
      >>;
  basis:=cdr basis;
  if basis
   then go to trynewloop;
  return nbasis.nreslist
  end;


symbolic procedure obvindep(row,matrx);
  % True if row is obviously linearly independent of the
  % Rows of the matrix.
begin scalar u;
  if null car matrx
    then return t;
  % no matrix => no dependence.
nexttry:
  if null row
    then return nil;
  if 0 iequal car row
    then go to nouse;
  u:=car matrx;
testloop:
  if 0 neq car u
    then go to nouse;
  u:=cdr u;
  if u
    then go to testloop;
  return t;
nouse:
  row:=cdr row;
  matrx:=cdr matrx;
  go to nexttry
  end;


symbolic procedure sf2df sf;
if null sf
   then nil
   else if numberp sf
    then (1 . sf).nil
    else begin
      scalar a,b,c;
      a:=sf2df lc sf;
      b:=(lpow sf .* 1) .+ nil;
      while a do <<
        c:=(!*multf(caar a,b).(cdar a)).c;
        a :=cdr a >>;
      return nconc(c,sf2df red sf)
      end;





symbolic procedure check!-sqrts!-dependence(sql,sqrtl);
% Resimplifies the list of SQs SQL,
% allowing for all dependencies among the
% sqrts in SQRTl.
begin
  scalar !*galois,sublist,sqrtsavelist,changeflag;
  sqrtsavelist:=listofallsqrts.listofnewsqrts;
  listofnewsqrts:=list mvar gaussiani;
  listofallsqrts:=list((argof mvar gaussiani) . gaussiani);
  !*galois:=t;
  for each u in sortsqrts(sqrtl,nil) do begin
    scalar v,uu;
    uu:=!*q2f simp argof u;
    v:=actualsimpsqrt uu;
    listofallsqrts:=(uu.v).listofallsqrts;
    if domainp v or mvar v neq u
      then <<
        if !*tra or !*trfield then <<
           printc u;
           printc "re-expressed as";
           printsf v >>;
        v:=prepf v;
        sublist:=(u.v) . sublist;
        changeflag:=t >>
    end;
  if changeflag then <<
    sql:=for each u in sql collect
           substitutesq(u,sublist);
    taylorasslist:=nil;
    sqrt!-places!-alist:=nil;
    basic!-listofallsqrts:=listofallsqrts;
    basic!-listofnewsqrts:=listofnewsqrts;
    if !*tra or !*trmin then <<
      printc "New set of residues are";
      mapc(sql,function printsq) >> >>
    else <<
      listofallsqrts:=car sqrtsavelist;
      listofnewsqrts:=cdr sqrtsavelist >>;
  return sql
  end;



symbolic procedure lindep(row,matrx);
  begin
    scalar m,m1,n,u,v,inverse,rowsinuse,failure;
    % Inverse is the answer from the "gaussian elimination"
    % we are doing.
    % Rowsinuse has nil for rows with no "awkward" non-zero entries.
    m1:=length car matrx;
    m:=isub1 m1;
    n:=isub1 length matrx;
    % n=length row.
    row:=mkvecf2q row;
    matrx:=mkvec for each j in matrx collect mkvecf2q j;
    inverse:=mkidenm m1;
    rowsinuse:=mkvect m;
    failure:=t;
    % initialisation complete.
    for i:=0 step 1 until n do begin
      % try to kill off i'th elements in each row.
      u:=nil;
      for j:=0 step 1 until m do <<
        % try to find a  pivot element.
        if  (null u) and
            (null getv(rowsinuse,j)) and
            (numr getv(getv(matrx,i),j))
          then u:=j >>;
      if null u
        then go to nullu;
      putv(rowsinuse,u,t);
      % it is no use trying this again ---
      % u is our pivot element.
      if u iequal m
        then go to nonetokill;
      for j:=iadd1 u step 1 until m do
        if numr getv(getv(matrx,i),j)
          then <<
            v:=negsq multsq(getv(getv(matrx,i),j),
                            invsq getv(getv(matrx,i),u));
            for k:=0 step 1 until m1 do
              putv(getv(inverse,k),j,
                addsq(getv(getv(inverse,k),j),
                  multsq(v,getv(getv(inverse,k),u))));
            for k:=0 step 1 until n do
              putv(getv(matrx,k),j,
                addsq(getv(getv(matrx,k),j),
                  multsq(v,getv(getv(matrx,k),u)))) >>;
      % We have now pivoted throughout matrix.
    nonetokill:
      % now do the same in row if necessary.
      if null numr getv(row,i)
        then go to norowop;
      v:=negsq multsq(getv(row,i),
                 invsq getv(getv(matrx,i),u));
      for k:=0 step 1 until m1 do
        putv(getv(inverse,k),m1,
          addsq(getv(getv(inverse,k),m1),
            multsq(v,getv(getv(inverse,k),u))));
      for k:=0 step 1 until n do
        putv(row,k,addsq(getv(row,k),
                     multsq(v,getv(getv(matrx,k),u))));
      u:=nil;
      for k:=0 step 1 until n do
        if numr getv(row,k)
          then u:=t;
      % if u is null then row is all 0.
      if null u
        then <<
          n:=-1;
          failure:=nil >>;
    norowop:
      if !*tra then <<
        princ "At end of cycle";
        printc row;
        printc matrx;
        printc inverse >>;
      return;
    nullu:
      % there is no pivot for this u.
      if numr getv(row,i)
        then n:=-1;
        % this element cannot be killed.
      end;
    if failure
      then return nil;
    v:=nil;
    for i:=0 step 1 until m do
      v:=(negsq getv(getv(inverse,m-i),m1)).v;
    return v
    end;

endmodule;

end;