module coatesid;
% Author: James H. Davenport.
fluid '(!*tra intvar magiclist taylorasslist taylorvariable);
exports coatessolve,vecprod,coates!-lineq;
imports !*invsq,!*multsq,negsq,!*addsq,swap,check!-lineq,non!-null!-vec,
printsq,sqrt2top,mapvec,mksp,vecsort,addsq,mkilist,mkvec,mapply,
taylorformp,xsubstitutesq,taylorform,taylorevaluate,multsq,
invsq,removecmsq;
symbolic procedure coatessolve(mzero,pzero,basis,normals);
begin
scalar m,n,rightside,nnn;
% if null normals
% then normals:=list mkilist(basis,1 ./ 1);
% This provides the default normalisation,
% viz merely a de-homogenising constraint;
% No it doesn't - JHD May 1983 and August 1986.
% This may be precisely the wrong constraint, as can be seen from
% the example of SQRT(X**2-1). Fixed 19/8/86 to amend COATESMATRIX
% to insert a normalising constraint if none is provided.
nnn:=max(length normals,1);
basis:=mkvec basis;
m:=coatesmatrix(mzero,pzero,basis,normals);
n:=upbv m;
rightside:=mkvect n;
for i:=0:n do
putv(rightside,n-i,(if i < nnn
then 1
else nil) ./ 1);
n:=coates!-lineq(m,rightside);
if n eq 'failed
then return 'failed;
n:=removecmsq vecprod(n,basis);
if !*tra then <<
printc "Answer from linear equation solving is ";
printsq n >>;
return n
end;
symbolic procedure coatesmatrix(mzero,pzero,basis,normals);
% NORMALS is a list of the normalising constraints
% that we must apply. Thypically, this is NIL, and we have to
% invent one - see the code IF NULL NORMALS ...
begin
scalar ans,n1,n2,j,w,save,nextflag,save!-taylors,x!-factors,
normals!-ok,temp;
save!-taylors:=mkvect isub1 length pzero;
save:=taylorasslist;
normals!-ok:=nil;
n1:=upbv basis;
n2:=isub1 mapply(function plus2,mzero) + max(length normals,1);
% the number of constaints in all (counting from 0).
taylorvariable:=intvar;
if !*tra then <<
printc "Basis for the functions with precisely the correct poles";
mapvec(basis,function printsq) >>;
ans:=mkvect n2;
for i:=0:n2 do
putv(ans,i,mkvect n1);
for i:=0:n1 do begin
scalar xmz,xpz,k;
xmz:=mzero;
k:=j:=0;
xpz:=pzero;
while xpz do <<
newplace basicplace car xpz;
if nextflag
then w:=taylorformp list('binarytimes,
getv(save!-taylors,k),
getv(x!-factors,k))
else if not !*tra
then w:=taylorform xsubstitutesq(getv(basis,i),car xpz)
else begin
scalar flg,u,slists;
u:=xsubstitutesq(getv(basis,i),basicplace car xpz);
slists:=extenplace car xpz;
for each w in sqrtsinsq(u,intvar) do
if not assoc(w,slists)
then flg:=w.flg;
if flg then <<
printc "The following square roots were not expected";
mapc(flg,function superprint);
printc "in the substitution";
superprint car xpz;
printsq getv(basis,i) >>;
w:=taylorform xsubstitutesq(u,slists)
end;
putv(save!-taylors,k,w);
k:=iadd1 k;
for l:=0 step 1 until isub1 car xmz do <<
astore(ans,j,i,taylorevaluate(w,l));
j:=iadd1 j >>;
if null normals and j=n2 then <<
temp:=taylorevaluate(w,car xmz);
astore(ans,j,i,temp);
% The defaults normalising condition is that the coefficient
% after the last zero be a non-zero.
% Unfortunately this too may fail (JHD 21.3.87) - check for it later
normals!-ok:=normals!-ok or numr temp >>;
xpz:=cdr xpz;
xmz:=cdr xmz >>;
nextflag:=(i < n1) and
(getv(basis,i) = multsq(!*kk2q intvar,getv(basis,i+1)));
if nextflag and null x!-factors then <<
x!-factors:=mkvect upbv save!-taylors;
xpz:=pzero;
k:=0;
xmz:=invsq !*kk2q intvar;
while xpz do <<
putv(x!-factors,k,taylorform xsubstitutesq(xmz,car xpz));
xpz:=cdr xpz;
k:=iadd1 k >> >>
end;
if null normals and null normals!-ok then <<
if !*tra
then printc "Our default normalisation condition was vacuous";
astore(ans,n2,n1,1 ./ 1)>>;
while normals do <<
w:=car normals;
for k:=0:n1 do <<
astore(ans,j,k,car w);
w:=cdr w >>;
j:=iadd1 j;
normals:=cdr normals >>;
tayshorten save;
return ans
end;
symbolic procedure printmatrix(ans,n2,n1);
if !*tra
then <<
printc "Equations to be solved:";
for i:=0:n2 do begin
if null getv(ans,i)
then return;
princ "Row number ";
princ i;
for j:=0:n1 do
printsq getv(getv(ans,i),j)
end >>;
symbolic procedure vecprod(u,v);
begin
scalar w,n;
w:=nil ./ 1;
n:=upbv u;
for i:=0:n do
w:=!*addsq(w,!*multsq(getv(u,i),getv(v,i)));
return w
end;
symbolic procedure coates!-lineq(m,rightside);
begin
scalar nnn,n;
nnn:=desparse(m,rightside);
if nnn eq 'failed
then return 'failed;
m:=car nnn;
if null m
then <<
n:=cddr nnn;
goto vecprod >>;
rightside:=cadr nnn;
nnn:=cddr nnn;
n:=check!-lineq(m,rightside);
if n eq 'failed
then return n;
n:=jhdsolve(m,rightside,non!-null!-vec nnn);
if n eq 'failed
then return n;
for i:=0:upbv n do
if (m:=getv(nnn,i))
then putv(n,i,m);
vecprod:
for i:=0:upbv n do
if null getv(n,i) then putv(n,i,nil ./ 1);
return n
end;
symbolic procedure jhdsolve(m,rightside,ignore);
% Returns answer to m.answer=rightside.
% Matrix m not necessarily square.
begin
scalar ii,n1,n2,ans,u,row,swapflg,swaps;
% The SWAPFLG is true if we have changed the order of the
% columns and need later to invert this via SWAPS.
n1:=upbv m;
for i:=0:n1 do
if (u:=getv(m,i))
then (n2:=upbv u);
printmatrix(m,n1,n2);
swaps:=mkvect n2;
for i:=0:n2 do
putv(swaps,i,n2-i);
% We have the SWAPS vector, which should be a vector of indices,
% arranged like this because VECSORT sorts in decreasing order.
for i:=0:isub1 n1 do begin
scalar k,v,pivot;
tryagain:
row:=getv(m,i);
if null row
then go to interchange;
% look for a pivot in row.
k:=-1;
for j:=0:n2 do
if numr (pivot:=getv(row,j))
then <<
k:=j;
j:=n2 >>;
if k neq -1
then goto newrow;
if numr getv(rightside,i)
then <<
m:='failed;
i:=sub1 n1; %Force end of loop.
go to finished >>;
% now interchange i and last element.
interchange:
swap(m,i,n1);
swap(rightside,i,n1);
n1:=isub1 n1;
if i iequal n1
then goto finished
else goto tryagain;
newrow:
if i neq k
then <<
swapflg:=t;
swap(swaps,i,k);
% record what we have done.
for l:=0:n1 do
swap(getv(m,l),i,k) >>;
% place pivot on diagonal.
pivot:=sqrt2top negsq !*invsq pivot;
for j:=iadd1 i:n1 do begin
u:=getv(m,j);
if null u
then return;
v:=!*multsq(getv(u,i),pivot);
if numr v then <<
putv(rightside,j,
!*addsq(getv(rightside,j),!*multsq(v,getv(rightside,i))));
for l:=0:n2 do
putv(u,l,!*addsq(getv(u,l),!*multsq(v,getv(row,l)))) >>
end;
finished:
end;
if m eq 'failed
then go to failed;
% Equations were inconsistent.
while null (row:=getv(m,n1)) do
n1:=isub1 n1;
u:=nil;
for i:=0:n2 do
if numr getv(row,i)
then u:='t;
if null u
then if numr getv(rightside,n1)
then go to failed
else n1:=isub1 n1;
% Deals with a last equation which is all zero.
if n1 > n2
then go to failed;
% Too many equations to satisfy.
ans:=mkvect n2;
n2:=n2 - ignore;
ii:=n1;
for i:=0 step 1 until n1 do
if null getv(m,i)
then ii:=iadd1 ii;
if ii < n2 then <<
if !*tra then <<
printc "The equations do not completely determine the functions";
printc "Matrix:";
mapvec(m,function superprint);
printc "Right-hand side:";
superprint rightside;
printc list("Adding new symbols for ", iadd1 ii," ... ",n2) >>;
% FOR I:=IADD1 ii:N2 DO <<
% U:=GENSYM();
% MAGICLIST:=U.MAGICLIST;
% PUTV(ANS,I,!*KK2Q U) >>;
if !*tra then printc "If in doubt consult an expert">>;
% now to do the back-substitution.
% Note that the matrix is not necessarily square,
% but that we have ensured that the non-square (underdetermiined)
% parts are at the right
for i:=n1 step -1 until 0 do begin
row:=getv(m,i);
if null row
then return;
% WHILE NULL NUMR GETV(ROW,II) DO II:=ISUB1 II;
ii:=0;
while null numr getv(row,ii) do ii:=iadd1 ii;
u:=getv(rightside,i);
for j:=iadd1 ii:n2 do <<
if null getv(ans,j) then <<
magiclist:=gensym().magiclist;
putv(ans,j,!*kk2q car magiclist) >>;
u:=!*addsq(u,!*multsq(getv(row,j),negsq getv(ans,j)))
>>;
putv(ans,ii,!*multsq(u,sqrt2top !*invsq getv(row,ii)));
% II:=ISUB1 II;
end;
if swapflg
then vecsort(swaps,list ans);
return ans;
failed:
if !*tra then printc "Unable to force correct zeroes";
return 'failed
end;
symbolic procedure desparse(matrx,rightside);
begin
scalar vect,changed,n,m,zero,failed;
zero := nil ./ 1;
n:=upbv matrx;
m:=upbv getv(matrx,0);
vect := mkvect m;
% for i:=0:m do putv(vect,i,zero); %%% initialize - ach
changed:=t;
while changed and not failed do begin
changed:=nil;
for i:=0:n do
if changed or failed
then i:=n % and hence quit the loop.
else begin
scalar nzcount,row,pivot;
row:=getv(matrx,i);
if null row
then return;
nzcount:=0;
for j:=0:m do
if numr getv(row,j)
then <<
nzcount:=iadd1 nzcount;
pivot:=j >>;
if nzcount = 0
then if null numr getv(rightside,i)
then return putv(matrx,i,nil)
else return (failed:='failed);
if nzcount > 1
then return nil;
nzcount:=getv(rightside,i);
if null numr nzcount
then <<
putv(vect,pivot,zero);
go to was!-zero >>;
nzcount:=!*multsq(nzcount,!*invsq getv(row,pivot));
putv(vect,pivot,nzcount);
nzcount:=negsq nzcount;
for i:=0:n do
if (row:=getv(matrx,i))
then if numr (row:=getv(row,pivot))
then putv(rightside,i,!*addsq(getv(rightside,i),
!*multsq(row,nzcount)));
was!-zero:
for i:=0:n do
if (row:=getv(matrx,i))
then putv(row,pivot,zero);
changed:=t;
putv(matrx,i,nil);
swap(matrx,i,n);
swap(rightside,i,n);
end;
end;
if failed
then return 'failed;
changed:=t;
for i:=0:n do
if getv(matrx,i)
then changed:=nil;
if changed
then matrx:=nil;
% We have completely solved the equations by these machinations.
return matrx.(rightside . vect)
end;
symbolic procedure astore(a,i,j,val);
putv(getv(a,i),j,val);
endmodule;
end;