Artifact 16b1e0d36b89a73826772bfed9b1c8c467c98eedffee370484ba4fdcbdb771d6:
- Executable file
r37/packages/factor/mhensfns.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: 14442) [annotate] [blame] [check-ins using] [more...]
module mhensfns; % Authors: A. C. Norman and P. M. A. Moore, 1979. fluid '(!*trfac alphalist current!-modulus degree!-bounds delfvec factor!-level factor!-trace!-list forbidden!-primes hensel!-growth!-size image!-factors max!-unknowns multivariate!-input!-poly non!-monic number!-of!-factors number!-of!-unknowns polyzero % pt prime!-base); %**********************************************************************; % This section contains some of the functions used in % the multivariate hensel growth. (ie they are called from % section MULTIHEN or function RECONSTRUCT-MULTIVARIATE-FACTORS). ; symbolic procedure set!-degree!-bounds v; degree!-bounds:=for each var in v collect (car var . degree!-in!-variable(multivariate!-input!-poly,car var)); symbolic procedure get!-degree!-bound v; begin scalar w; w:=atsoc(v,degree!-bounds); if null w then errorf(list("Degree bound not found for ", v," in ",degree!-bounds)); return cdr w end; symbolic procedure choose!-larger!-prime n; % our prime base in the multivariate hensel must be greater than n so % this sets a new prime to be that (previous one was found to be no % good). We also set up various fluids e.g. the Alphas; % the primes we can choose are < 2**24 so if n is bigger % we collapse; if n > 2**24-1 then errorf list("CANNOT CHOOSE PRIME > GIVEN NUMBER:",n) else begin scalar p,flist!-mod!-p,k,fvec!-mod!-p,forbidden!-primes; trynewprime: if p then forbidden!-primes:=p . forbidden!-primes; p:=random!-prime(); % this chooses a word-size prime (currently 24 bits); set!-modulus p; if not(p>n) or member(p,forbidden!-primes) or polyzerop reduce!-mod!-p lc multivariate!-input!-poly then goto trynewprime; for i:=1:number!-of!-factors do flist!-mod!-p:=(reduce!-mod!-p getv(image!-factors,i) . flist!-mod!-p); alphalist:=alphas(number!-of!-factors,flist!-mod!-p,1); if alphalist='factors! not! coprime then goto trynewprime; hensel!-growth!-size:=p; prime!-base:=p; factor!-trace << prin2!* "New prime chosen: "; printstr hensel!-growth!-size >>; k:=number!-of!-factors; fvec!-mod!-p:=mkvect k; for each w in flist!-mod!-p do << putv(fvec!-mod!-p,k,w); k:=isub1 k >>; return fvec!-mod!-p end; symbolic procedure binomial!-coefft!-mod!-p(n,r); if n<r then nil else if n=r then 1 else if r=1 then !*n2f modular!-number n else begin scalar n!-c!-r,b,j; n!-c!-r:=1; b:=min(r,n-r); n:=modular!-number n; r:=modular!-number r; for i:=1:b do << j:=modular!-number i; n!-c!-r:=modular!-quotient( modular!-times(n!-c!-r, modular!-difference(n,modular!-difference(j,1))), j) >>; return !*n2f n!-c!-r end; symbolic procedure make!-multivariate!-hatvec!-mod!-p(bvec,n); % makes a vector whose ith elt is product over j [ BVEC(j) ] / BVEC(i); % NB. we must NOT actually do the division here as we are likely % to be working mod p**n (some n > 1) and the division can involve % a division by p.; begin scalar bhatvec,r; bhatvec:=mkvect n; for i:=1:n do << r:=1; for j:=1:n do if not(j=i) then r:=times!-mod!-p(r,getv(bvec,j)); putv(bhatvec,i,r) >>; return bhatvec end; symbolic procedure max!-degree!-in!-var(fvec,v); begin scalar r,d; r:=0; for i:=1:number!-of!-factors do if r<(d:=degree!-in!-variable(getv(fvec,i),v)) then r:=d; return r end; symbolic procedure make!-growth!-factor pt; % pt is of form (v . n) where v is a variable. we make the s.f. v-n; if cdr pt=0 then !*f2mod !*k2f car pt else plus!-mod!-p(!*f2mod !*k2f car pt,modular!-minus cdr pt); symbolic procedure terms!-done!-mod!-p(fvec,delfvec,delfactor); % calculate the terms introduced by the corrections in DELFVEC; begin scalar flist,delflist; for i:=1:number!-of!-factors do << flist:=getv(fvec,i) . flist; delflist:=getv(delfvec,i) . delflist >>; return terms!-done1!-mod!-p(number!-of!-factors,flist,delflist, number!-of!-factors,delfactor) end; symbolic procedure terms!-done1!-mod!-p(n,flist,delflist,r,m); if n=1 then (car flist) . (car delflist) else begin scalar k,i,f1,f2,delf1,delf2; k:=n/2; i:=1; for each f in flist do << if i>k then f2:=(f . f2) else f1:=(f . f1); i:=i+1 >>; i:=1; for each delf in delflist do << if i>k then delf2:=(delf . delf2) else delf1:=(delf . delf1); i:=i+1 >>; f1:=terms!-done1!-mod!-p(k,f1,delf1,r,m); delf1:=cdr f1; f1:=car f1; f2:=terms!-done1!-mod!-p(n-k,f2,delf2,r,m); delf2:=cdr f2; f2:=car f2; delf1:= plus!-mod!-p(plus!-mod!-p( times!-mod!-p(f1,delf2), times!-mod!-p(f2,delf1)), times!-mod!-p(times!-mod!-p(delf1,m),delf2)); if n=r then return delf1; return (times!-mod!-p(f1,f2) . delf1) end; symbolic procedure primitive!.parts(flist,var,univariate!-inputs); % finds the prim.part of each factor in flist wrt variable var; % Note that FLIST may contain univariate or multivariate S.F.s % (according to UNIVARIATE!-INPUTS) - in the former case we correct the % ALPHALIST if necessary; begin scalar c,primf; if null var then errorf "Must take primitive parts wrt some non-null variable"; if non!-monic then factor!-trace << printstr "Because we multiplied the original primitive"; printstr "polynomial by a multiple of its leading coefficient"; printstr "(see (a) above), the factors we have now are not"; printstr "necessarily primitive. However the required factors"; printstr "are merely their primitive parts." >>; return for each fw in flist collect << if not depends!-on!-var(fw,var) then errorf list("WRONG VARIABLE",var,fw); c:=comfac fw; if car c then errorf(list( "FACTOR DIVISIBLE BY MAIN VARIABLE:",fw,car c)); primf:=quotfail(fw,cdr c); if not(cdr c=1) and univariate!-inputs then multiply!-alphas(cdr c,fw,primf); primf >> end; symbolic procedure make!-predicted!-forms(pfs,v); % PFS is a vector of S.F.s which represents the sparsity of % the associated polynomials wrt V. Here PFS is adjusted to a % suitable form for handling this sparsity. ie. we record the % degrees of V in a vector for each poly in PFS. Each % monomial (in V) represents an unknown (its coefft) in the predicted % form of the associated poly. We count the maximum no of unknowns for % each poly and return the maximum of these; begin scalar l,n,pvec,j,w; max!-unknowns:=0; for i:=1:number!-of!-factors do << w:=getv(pfs,i); % get the ith poly; l:=sort(spreadvar(w,v,nil),function lessp); % Pick out the monomials in V from this poly and order % them in increasing degree; n:=iadd1 length l; % no of unknowns in predicted poly - we add % one for the constant term; number!-of!-unknowns:=(n . i) . number!-of!-unknowns; if max!-unknowns<n then max!-unknowns:=n; pvec:=mkvect isub1 n; % get space for the info on this poly; j:=0; putv(pvec,j,isub1 n); % put in the length of this vector which will vary % from poly to poly; for each m in l do putv(pvec,j:=iadd1 j,m); % put in the monomial info; putv(pfs,i,pvec); % overwrite the S.F. in PFS with the more compact vector; >>; number!-of!-unknowns:=sort(number!-of!-unknowns,function lesspcar); return max!-unknowns end; symbolic procedure make!-correction!-vectors(bfs,n); % set up space for the vector of vectors to hold the correction % terms as we generate them by the function SOLVE-FOR-CORRECTIONS. % Also put in the starting values; begin scalar cvs,cv; cvs:=mkvect number!-of!-factors; for i:=1:number!-of!-factors do << cv:=mkvect n; % each CV will hold the corrections for the ith factor; % the no of corrections we put in here depends on the % maximum no of unknowns we have in the predicted % forms, giving a set of soluble linear systems (hopefully); putv(cv,1,getv(bfs,i)); % put in the first 'corrections'; putv(cvs,i,cv) >>; return cvs end; symbolic procedure construct!-soln!-matrices(pfs,val); % Here we construct the matrices - one for each linear system % we will have to solve to see if our predicted forms of the % answer are correct. Each matrix is a vector of row-vectors % - the ijth elt is in jth slot of ith row-vector (ie zero slots % are not used here); begin scalar soln!-matrix,resvec,n,pv; resvec:=mkvect number!-of!-factors; for i:=1:number!-of!-factors do << pv:=getv(pfs,i); soln!-matrix:=mkvect(n:=iadd1 getv(pv,0)); construct!-ith!-matrix(soln!-matrix,pv,n,val); putv(resvec,i,soln!-matrix) >>; return resvec end; symbolic procedure construct!-ith!-matrix(sm,pv,n,val); begin scalar mv; mv:=mkvect n; % this will be the first row; putv(mv,1,1); % the first column represents the constant term; for j:=2:n do putv(mv,j,modular!-expt(val,getv(pv,isub1 j))); % first row is straight substitution; putv(sm,1,mv); % now for the rest of the rows: ; for j:=2:n do << mv:=mkvect n; putv(mv,1,0); construct!-matrix!-row(mv,isub1 j,pv,n,val); putv(sm,j,mv) >> end; symbolic procedure construct!-matrix!-row(mrow,j,pv,n,val); begin scalar d; for k:=2:n do << d:=getv(pv,isub1 k); % degree representing the monomial; if d<j then putv(mrow,k,0) else << d:=modular!-times(!*d2n binomial!-coefft!-mod!-p(d,j), modular!-expt(val,idifference(d,j))); % differentiate and substitute all at once; putv(mrow,k,d) >> >> end; symbolic procedure print!-linear!-systems(soln!-m,correction!-v, predicted!-f,v); << for i:=1:number!-of!-factors do print!-linear!-system(i,soln!-m,correction!-v,predicted!-f,v); terpri!*(nil) >>; symbolic procedure print!-linear!-system(i,soln!-m,correction!-v, predicted!-f,v); begin scalar pv,sm,cv,mr,n,tt; terpri!*(t); prin2!* " i = "; printstr i; terpri!*(nil); sm:=getv(soln!-m,i); cv:=getv(correction!-v,i); pv:=getv(predicted!-f,i); n:=iadd1 getv(pv,0); for j:=1:n do << % for each row in matrix ... ; prin2!* "( "; tt:=2; mr:=getv(sm,j); % matrix row; for k:=1:n do << % for each elt in row ... ; prin2!* getv(mr,k); ttab!* (tt:=tt+10) >>; prin2!* ") ( ["; if j=1 then prin2!* 1 else prinsf adjoin!-term(mksp(v,getv(pv,isub1 j)),1,polyzero); prin2!* "]"; ttab!* (tt:=tt+10); prin2!* " )"; if j=(n/2) then prin2!* " = ( " else prin2!* " ( "; prinsf getv(cv,j); ttab!* (tt:=tt+30); printstr ")"; if not(j=n) then << tt:=2; prin2!* "("; ttab!* (tt:=tt+n*10); prin2!* ") ("; ttab!* (tt:=tt+10); prin2!* " ) ("; ttab!* (tt:=tt+30); printstr ")" >> >>; terpri!*(t) end; symbolic procedure try!-prediction(sm,cv,pv,n,i,poly,v,ff,ffhat); begin scalar w,ffi,fhati; sm:=getv(sm,i); cv:=getv(cv,i); pv:=getv(pv,i); if not(n=iadd1 getv(pv,0)) then errorf list("Predicted unknowns gone wrong? ",n,iadd1 getv(pv,0)); if null getm2(sm,1,0) then << w:=lu!-factorize!-mod!-p(sm,n); if w='singular then << factor!-trace << prin2!* "Prediction for "; prin2!* if null ff then 'f else 'a; prin2!* "("; prin2!* i; printstr ") failed due to singular matrix." >>; return (w . i) >> >>; back!-substitute(sm,cv,n); w:= if null ff then try!-factor(poly,cv,pv,n,v) else << ffi := getv(ff,i); fhati := getv(ffhat,i); % The unfolding here is to get round % a bug in the PSL compiler 12/9/82. It % will be tidied back up as soon as % possible; try!-alpha(poly,cv,pv,n,v,ffi,fhati) >>; if w='bad!-prediction then << factor!-trace << prin2!* "Prediction for "; prin2!* if null ff then 'f else 'a; prin2!* "("; prin2!* i; printstr ") was an inadequate guess." >>; return (w . i) >>; factor!-trace << prin2!* "Prediction for "; prin2!* if null ff then 'f else 'a; prin2!* "("; prin2!* i; prin2!* ") worked: "; printsf car w >>; return (i . w) end; symbolic procedure try!-factor(poly,testv,predictedf,n,v); begin scalar r,w; r:=getv(testv,1); for j:=2:n do << w:=!*f2mod adjoin!-term(mksp(v,getv(predictedf,isub1 j)),1, polyzero); r:=plus!-mod!-p(r,times!-mod!-p(w,getv(testv,j))) >>; w:=quotient!-mod!-p(poly,r); if didntgo w or not polyzerop difference!-mod!-p(poly,times!-mod!-p(w,r)) then return 'bad!-prediction else return list(r,w) end; symbolic procedure try!-alpha(poly,testv,predictedf,n,v,fi,fhati); begin scalar r,w,wr; r:=getv(testv,1); for j:=2:n do << w:=!*f2mod adjoin!-term(mksp(v,getv(predictedf,isub1 j)),1, polyzero); r:=plus!-mod!-p(r,times!-mod!-p(w,getv(testv,j))) >>; if polyzerop (wr:=difference!-mod!-p(poly,times!-mod!-p(r,fhati))) then return list (r,wr); w:=quotient!-mod!-p(wr,fi); if didntgo w or not polyzerop difference!-mod!-p(wr,times!-mod!-p(w,fi)) then return 'bad!-prediction else return list(r,wr) end; endmodule; end;