Artifact 42d7e7f5e03db34ab3848f68fb4406cc279b82fddcec3bd0ae00b192551cd818:
- Executable file
r37/packages/algint/coates.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: 15245) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/algint/coates.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: 15245) [annotate] [blame] [check-ins using]
module coates; % Author: James H. Davenport. fluid '(!*tra !*trmin !*galois intvar magiclist nestedsqrts previousbasis sqrt!-intvar taylorasslist thisplace listofallsqrts listofnewsqrts basic!-listofallsqrts basic!-listofnewsqrts gaussiani !*trfield); global '(coates!-fdi); exports coates,makeinitialbasis,checkpoles,multbyallcombinations; symbolic procedure coates(places,mults,x); begin scalar u,tt; tt:=readclock(); u:=coates!-hpfsd(places,mults); if !*tra or !*trmin then printc list ('coates,'time,readclock()-tt,'milliseconds); return u end; symbolic procedure coates!-real(places,mults); begin scalar thisplace,u,v,save; if !*tra or !*trmin then << princ "Find function with zeros of order:"; printc mults; if !*tra then princ " at "; terpri!*(t); if !*tra then mapc(places,function printplace) >>; % v:=placesindiv places; % V is a list of all the substitutors in PLACES; % u:=mkunique sqrtsintree(v,intvar,nil); % if !*tra then << % princ "Sqrts on this curve:"; % terpri!*(t); % superprint u >>; % algnos:=mkunique for each j in places collect basicplace j; % if !*tra then << % printc "Algebraic numbers where residues occur:"; % superprint algnos >>; v:=mults; for each uu in places do << if (car v) < 0 then u:=(rfirstsubs uu).u; v:=cdr v >>; thisplace:=list('quotient,1,intvar); if member(thisplace,u) then << v:= finitise(places,mults); % returns list (places,mults,power of intvar to remove. u:=coates!-real(car v,cadr v); if atom u then return u; return multsq(u,!*p2q mksp(intvar,caddr v)) >>; % It is not sufficient to check the current value of U in FRACTIONAL... % as we could have zeros over infinity JHD 18/8/86; for each uu in places do if rfirstsubs uu = thisplace then u:=append(u,mapovercar cdr uu); coates!-fdi:=fractional!-degree!-at!-infinity u; % Do we need to blow everything up by a factor of two (or more) % to avoid fractional powers at infinity? if coates!-fdi iequal 1 then return coatesmodule(places,mults,intvar); if !*tra then fdi!-print(); newplace list(intvar . thisplace, list(intvar,'expt,intvar,coates!-fdi)); places := for each j in places collect fdi!-upgrade j; save:=taylorasslist; u:=coatesmodule(places, for each u in mults collect u*coates!-fdi, intvar); taylorasslist:=save; % u:=fdi!-revertsq u; % That previous line is junk, I think (JHD 22.8.86) % just because we blew up the places doesn't mean that % we should deflate the function, because that has already been done. return u end; symbolic procedure coatesmodule(places,mults,x); begin scalar pzero,mzero,u,v,basis,sqrts,magiclist,mpole,ppole; % MAGICLIST holds the list of extra unknowns created in JHDSOLVE % which must be found in CHECKPOLES (calling FINDMAGIC). sqrts:=sqrtsinplaces places; if !*tra then << princ "Sqrts on this curve:"; superprint sqrts >>; u:=places; v:=mults; while u do << if 0<car v then << mzero:=(car v).mzero; pzero:=(car u).pzero >> else << mpole:=(car v).mpole; ppole:=(car u).ppole >>; u:=cdr u; v:=cdr v >>; % ***time-hack-2***; if previousbasis then basis:=previousbasis else basis:=mkvec makeinitialbasis ppole; u:=completeplaces(ppole,mpole); basis:=integralbasis(basis,car u,cdr u,x); basis:=normalbasis(basis,x,0); u:=coatessolve(mzero,pzero,basis,nil); % The NIL is the list of special constraints needed % to force certain poles to occur in the answer. if atom u then return u; v:= checkpoles(list u,places,mults); if null v then return 'failed; if not magiclist then return u; u:=removecmsq substitutesq(u,v); % Apply the values from FINDMAGIC. if !*tra or !*trmin then << printc "These values give the function"; printsq u >>; magiclist:=nil; if checkpoles(list u,places,mults) then return u else interr "Inconsistent checkpoles" end; symbolic procedure makeinitialbasis places; begin scalar u; u:=multbyallcombinations(list(1 ./ 1), for each u in getsqrtsfromplaces places collect !*kk2q u); if !*tra then << printc "Initial basis for the space m(x)"; mapc(u,function printsq) >>; return u end; symbolic procedure multbyallcombinations(u,l); % Produces a list of all elements of u, % each multiplied by every combination of elements of l. if null l then u else multbyallcombinations(nconc(multsql(car l,u),u),cdr l); symbolic procedure multsql(u,l); % Multiplies (!*multsq) each element of l by u. if null l then nil else !*multsq(u,car l) . multsql(u,cdr l); symbolic procedure checkpoles(basis,places,mults); % Checks that the BASIS really does have all the % poles in (PLACES.MULTS). begin scalar u,v,l; go to outer2; outer: places:=cdr places; mults:=cdr mults; outer2: if null places then return if magiclist then findmagic l else t; if 0 leq car mults then go to outer; u:=basis; inner: if null u then << if !*tra then << princ "The answer from the linear equations did"; printc " not have the poles at:"; printplace car places >>; return nil >>; v:=taylorform xsubstitutesq(car u,car places); if taylorfirst v=car mults then << if magiclist then l:=taylorevaluate(v,car mults) . l; go to outer >>; if taylorfirst v < car mults then interr "Extraneous pole introduced"; u:=cdr u; go to inner end; symbolic procedure coates!-hpfsd(oplaces,omults); begin scalar mzero,pzero,mpole,ppole,fun,summzero,answer,places,mults; places:=oplaces; mults:=omults; % Keep originals in case need to use COATES!-REAL directly. summzero:=0; % holds the sum of all the mzero's. while places do << if 0<car mults then << summzero:=summzero + car mults; mzero:=(car mults).mzero; pzero:=(car places).pzero >> else << mpole:=(car mults).mpole; ppole:=(car places).ppole >>; places:=cdr places; mults:=cdr mults >>; if summzero > 2 then begin % We want to combine a zero/pole pair % so as to reduce the total index before calling coates!-real % on the remaining zeros/poles. scalar nplaces,nmults,f,multiplicity,newpole,sqrts,fz,zfound,mult1; sqrts:=getsqrtsfromplaces ppole; if !*tra or !*trmin then << princ "Operate on divisor:"; printc append(mzero,mpole); printc "at"; mapc(pzero,function printplace); mapc(ppole,function printplace) >>; iterate: nplaces:=list car pzero; multiplicity:=car mzero; nmults:=list 1; if cdr ppole then << nplaces:=(car ppole) . ( (cadr ppole) . nplaces); multiplicity:=min(multiplicity,- car mpole,- cadr mpole); nmults:=(-1) .((-1) . nmults) >> else << nplaces:=(car ppole) . nplaces; multiplicity:=min(multiplicity,(- car mpole)/2); nmults:=(-2) . nmults >>; previousbasis:=nil; f:=coates!-real(nplaces,nmults); if atom f then << if !*tra or !*trmin then printc "Failure: must try whole divisor"; return coates!-real(oplaces,omults) >>; % newpole:=removezero(findzeros(f,sqrts),car pzero). fz:=findzeros(f,sqrts,2); zfound:=assoc(car pzero,fz); if not zfound then interr list("Didn't seem to find the zero", car pzero, "we looked for"); if cdr zfound > car mzero then interr "We found too many zeros"; fz:=delete(zfound,fz); if !*tra or !*trmin then << printc "Replaced by the pole"; if fz then prettyprint fz else <<terpri(); prin2t "The zero we were already looking for">>; princ multiplicity; printc " times" >>; mult1:=car mzero - multiplicity * cdr zfound; if mult1 < 0 then <<if !*tra then printc "*** A zero has turned into a pole"; multiplicity:= car mzero / cdr zfound ; mult1:=remainder(car mzero, cdr zfound); >>; if mult1=0 then << mzero:=cdr mzero; pzero:=cdr pzero >> else rplaca(mzero,mult1); if null cdr ppole then << if (car mpole + 2*multiplicity)=0 then << ppole:=cdr ppole; mpole:=cdr mpole >> else rplaca(mpole,car mpole + 2 * multiplicity) >> else << if (cadr mpole + multiplicity)=0 then << ppole:=(car ppole) . (cddr ppole); mpole:=(car mpole) . (cddr mpole) >> else rplaca(cdr mpole,cadr mpole + multiplicity); if (car mpole + multiplicity)=0 then << ppole:=cdr ppole; mpole:=cdr mpole >> else rplaca(mpole,car mpole + multiplicity) >>; while fz do << newpole:=caar fz; mult1:=multiplicity*(cdar fz); if newpole member pzero then begin scalar m,p; while newpole neq car pzero do << m:=(car mzero).m; mzero:=cdr mzero; p:=(car pzero).p; pzero:=cdr pzero >>; if mult1 < car mzero then << mzero:=(car mzero - mult1) . cdr mzero; mzero:=nconc(m,mzero); pzero:=nconc(p,pzero); return >>; if mult1 > car mzero then << ppole:=newpole.ppole; mpole:=(car mzero - mult1) . mpole >>; mzero:=nconc(m,cdr mzero); pzero:=nconc(p,cdr pzero) end else if newpole member ppole then begin scalar m,p; m:=mpole; p:=ppole; while newpole neq car p do << p:=cdr p; m:=cdr m >>; rplaca(m,car m - mult1) end else << mpole:=nconc(mpole,list(-mult1)); ppole:=nconc(ppole,list newpole) >>; fz:=cdr fz >>; f:=mk!*sq f; if multiplicity > 1 then answer:=list('expt,f,multiplicity).answer else answer:=f.answer; summzero:=0; for each x in mzero do summzero:=summzero+x; if !*tra then << princ "Function is now: "; printc append(mzero,mpole); printc "at"; mapc(pzero,function printplace); mapc(ppole,function printplace) >>; if summzero > 2 then go to iterate; end; if pzero or ppole then << % there's something left to do fun:=coates!-real(nconc(pzero,ppole), nconc(mzero,mpole)); if null answer then return fun else answer:=(mk!*sq fun).answer >>; return !*k2q('times.answer); % This is not valid, but we hope that it will be unpicked; % (e.g. by SIMPLOG) before too much damage is caused. end; % symbolic procedure removezero(l,place); % if place member l % then (lambda u; if null cdr u then car u % else interr "Removezero") delete(place,l) % else interr "Error in removezeros"; symbolic procedure findzeros(sq,sqrts,nzeros); % NZEROS is the number of zeros known, a priori, to exist begin scalar u,potentials,answer,n,not!-answer,nz,series; u:=denr sqrt2top invsq sq; potentials:=for each v in jfactor(u,intvar) collect begin scalar w,place; w:=makemainvar(numr v,intvar); if ldeg w neq 1 then interr "Can't cope"; if red w then place:=list(intvar,'plus,intvar,prepsq(negf red w ./ lc w)) else place:=intvar . intvar; % This IF .. ELSE .. added JHD 3 Sept 1980. return place end; potentials:=list(intvar,'quotient,1,intvar).potentials; for each place in potentials do begin scalar slist,nestedsqrts; place:=list place; newplace place; u:=substitutesq(sq,place); while involvesq(u,sqrt!-intvar) do begin scalar z; z:=list list(intvar,'expt,intvar,2); place:=nconc(place,z); newplace place; u:=substitutesq(u,z); end; slist:=sqrtsinsq(u,intvar); for each v in sqrts do slist:=union(slist,sqrtsinsq(xsubstitutesq(!*kk2q v,place), intvar)); slist:=sqrtsign(slist,intvar); for each s in slist do if (n:=taylorfirst (series:=taylorform substitutesq(u,s))) > 0 then answer:=(append(place,s).n).answer else not!-answer:=list(u,place,s,series) . not!-answer; return answer; end; nz:= for each u in answer sum cdr u; if nz = nzeros then return answer; if nz > nzeros then interr "We have too many zeros"; if !*tra then printc "Couldn't find enough zeros of the function: try harder"; for each v in not!-answer do begin scalar !*galois,sqrtsavelist,sublist,s; % If we don't reset taylorasslist, then all calculations are % dubious! taylorasslist:=nil; !*galois:=t; sqrtsavelist:=listofallsqrts.listofnewsqrts; listofnewsqrts:=list mvar gaussiani; listofallsqrts:=list((argof mvar gaussiani) . gaussiani); series:=cadddr v; s:=cdr assoc(taylorfirst series,taylorlist series); for each u in sortsqrts(sqrtsinsq(s,nil),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 >>; basic!-listofnewsqrts := union(sqrtsinsf(v,nil,nil), basic!-listofnewsqrts); basic!-listofallsqrts := car listofallsqrts . basic!-listofallsqrts; v:=prepf v; sublist:= (u.v) .sublist; sqrtsavelist:=( car listofallsqrts . delete(assoc(uu,car sqrtsavelist), car sqrtsavelist)). delete(u,cdr sqrtsavelist)>>; end; listofallsqrts:=car sqrtsavelist; listofnewsqrts:=cdr sqrtsavelist; if sublist and null numr substitutesq(s,sublist) then << if !*tra or !*trfield then printc "a non-zero term has become zero"; !*galois:=nil; % We'll have done all we wanted, and mustn't confuse if (n:=taylorfirst taylorform substitutesq(car v,caddr v)) > 0 then << answer:= (append(cadr v,caddr v).n) . answer; nz:= nz+n ; if !*tra or !*trfield then printc "that found us a new zero of the function" >>>>; end; if nz = nzeros then return answer; interr "can't find enough zeros"; end; endmodule; end;