Artifact 81f3e18bd7770fc78d00a7e113b2d6438f49df1ad5dfa8d9e2b5272f6f073de3:
- Executable file
r37/packages/algint/coatesid.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: 11409) [annotate] [blame] [check-ins using] [more...]
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;