Artifact dc0d74ff35eafff4ae8936472ac97a64e2a340af27cea0a88ff1f3021414618d:
- Executable file
r38/packages/solve/modsolve.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: 5651) [annotate] [blame] [check-ins using] [more...]
module modsolve; % Solve modular. % Author: Herbert Melenk <melenk@zib.de> % Algebraic interface: m_solve(eqn/eqnlist [,variables]). % Some routines from solve and factor(modpoly) are needed. fluid'(!*trnonlnr current!-modulus); % The limit '10000000' for the current modulus has been calculated % using a 500 MHz 686 machine which needs 115 seconds for a square root % computation for that limit. For faster machines the limit may be set % a bit higher, for slower machines it should be set lower. load!-package 'solve; load!-package 'factor; put('m_solve,'psopfn,'msolve); symbolic procedure msolve(u); begin scalar s,s1,v,v1,w; s:=reval car u; s:=if eqcar(s,'list) then cdr s else {s}; if cdr u then <<v:= reval cadr u; v:=if eqcar(v,'list) then cdr v else {v}>>; % test, collect variables. s1:=for each q in s collect <<if eqcar(q,'equal) then q:='difference.cdr q; w:=numr simp q./1; v1:=union(v1,solvevars{w}); numr w>>; if null v then v:=v1; return msolve!-result if null cdr s1 then msolve!-poly(car s1,v) else msolve!-psys(s1,v) end; symbolic procedure msolve!-result u; if u='failed then u else 'list.for each v in u collect 'list.for each w in v collect{'equal,car w,cdr w}; symbolic procedure msolvesys(s1,v,tg); % Interface for the Solve package. begin scalar w,fail; if null cdr s1 then <<w:=msolve!-poly(car s1,v); go to done>>; % Reject parametric modular equation system. for each p in s1 do for each x in kernels p do if not member(x,v) then fail:=t; if fail then <<if !*trnonlnr then lprim "Cannot solve parametric modular system"; go to failed>>; w:=msolve!-psys(s1,v); done: if w='failed then go to failed; w:=for each q in w collect {for each r in q collect simp cdr r, for each r in q collect car r,1}; return if tg then t.w else w; failed: return if null cdr s1 and null cdr v and null tg then mkrootsof(car s1 ./1,car v,1) else if tg then '(failed) else 'failed end; symbolic procedure msolve!-poly1(f,x); % polynomial f(x); begin scalar w,l; if domainp f then nil else if ldeg f=1 then <<w:=safe!-modrecip lc f; erfg!*:=nil; if null w then go to enum; w:=moduntag multf(w,negf red f); if w and (w<0 or w>current!-modulus) then w:=general!-modular!-number w; w:={w}; go to done>>; enum: l:=lowestdeg(f,x,0); if l>0 then f:=quotf(f,numr simp{'expt,x,l}); f:=general!-reduce!-mod!-p moduntag f; w:=for i:=1:current!-modulus -1 join if null general!-evaluate!-mod!-p(f,x,i) then {i}; if l>0 then w:=append(w,{nil}); done: return for each q in w collect{x.prepf q} end; symbolic procedure msolve!-poly(f,l); % Solve one polynomial wrt several variables. begin scalar x,vl,limit; limit:=10000000; %%%%%%%%%%%%%%%%%%%%%%%%%%%% limit if current!-modulus>limit then <<if !*trnonlnr then lprim {"Current modulus larger than",limit}; return 'failed>>; vl:=kernels f; for each x in l do <<if not member(x,vl) then l:=delete(x,l); vl:=delete(x,vl)>>; if null l then return nil; return if vl then msolve!-polya(f,l) else msolve!-polyn(f,l) end; symbolic procedure msolve!-polyn(f,l); (if null cdr l then msolve!-poly1(f,car l) else for i:=0:current!-modulus -1 join for each s in msolve!-polyn(numr subf(f,{x.i}),cdr l) collect (x.i).s) where x=car l; symbolic procedure msolve!-polya(f,l); % 'f' is a polynomial with variables in 'l' and at least one more % formal parameter. 'f' can be solved only if 'f' is linear in one of % the variables with an invertible coefficient. Otherwise we must % return a root-of expression. begin scalar x,c,w; for each y in l do if null x then if 1=ldeg ((w:=reorder f) where kord!*={y}) then x:=y; if null x then go to none; c:=lc w; w:=red w; if not domainp c then go to none; c:=safe!-modrecip c; if null c then go to none; return{{x.prepf multf(negf w,c)}}; none: return{{car l.mk!*sq caaar mkrootsof(f./1,car l,1)}} end; symbolic procedure msolve!-psys(s,v); % Solve system 's' for variables 'v'. 's' has no additional free % parameters. begin scalar b,o,z,w; if current!-modulus * length s>1000 and primep current!-modulus then <<% Domain is a field and big problem - compute a Groebner base % first. load!-package 'groebner;load!-package 'groebnr2; o:=apply1('torder,{'list.v,'lex}); b:=groebnereval{'list.for each p in s collect prepf p}; z:=gzerodimeval{b}; % The reverse basis for increasing variable number. s:=reversip for each p in cdr b collect numr simp p; apply1('torder,cdr o)>> else <<% Rearrange system for increasing variable number. w:=for each p in s collect length(for each x in v join if smemq(x,p)then{x}).p; w:=for each p in sort(w,'lesspcar) collect cdr p>>; return msolve!-psys1(s,v)end; symbolic procedure msolve!-psys1(s,v); % Solve system by successive substitution. begin scalar w,w1,f,f1; w:={nil}; for each f in s do <<w1:=nil; for each s in w do <<f1:=general!-reduce!-mod!-p moduntag numr subf(f,s); if null f1 then w1:=s.w1 else if domainp f1 then nil else for each ns in msolve!-poly(f1,v)do w1:=append(s,ns).w1>>; w:=w1>>; return w end; endmodule; end;