module multihen; % Hensel construction for the multivariate case.
% (This version is highly recursive.)
% Authors: A. C. Norman and P. M. A. Moore, 1979.
fluid '(!*overshoot
!*trfac
alphavec
bad!-case
factor!-level
factor!-trace!-list
fhatvec
hensel!-growth!-size
max!-unknowns
number!-of!-factors
number!-of!-unknowns
predictions);
symbolic procedure find!-multivariate!-factors!-mod!-p(poly,
best!-factors,variable!-set);
% All arithmetic is done mod p, best-factors is overwritten.
if null variable!-set then best!-factors
else (lambda factor!-level; begin
scalar growth!-factor,b0s,res,v,
bhat0s,w,degbd,first!-time,redpoly,
predicted!-forms,number!-of!-unknowns,solve!-count,
correction!-vectors,soln!-matrices,max!-unknowns,
unknowns!-count!-list,poly!-remaining,
prediction!-results,one!-prediction!-failed;
v:=car variable!-set;
degbd:=get!-degree!-bound car v;
first!-time:=t;
growth!-factor:=make!-growth!-factor v;
poly!-remaining:=poly;
prediction!-results:=mkvect number!-of!-factors;
find!-msg1(best!-factors,growth!-factor,poly);
b0s:=reduce!-vec!-by!-one!-var!-mod!-p(best!-factors,
v,number!-of!-factors);
% The above made a copy of the vector.
for i:=1:number!-of!-factors do
putv(best!-factors,i,
difference!-mod!-p(getv(best!-factors,i),getv(b0s,i)));
redpoly:=evaluate!-mod!-p(poly,car v,cdr v);
find!-msg2(v,variable!-set);
find!-multivariate!-factors!-mod!-p(redpoly,b0s,cdr variable!-set);
% answers in b0s.
if bad!-case then return;
for i:=1:number!-of!-factors do
putv(best!-factors,i,
plus!-mod!-p(getv(b0s,i),getv(best!-factors,i)));
find!-msg3(best!-factors,v);
res:=diff!-over!-k!-mod!-p(
difference!-mod!-p(poly,
times!-vector!-mod!-p(best!-factors,number!-of!-factors)),
1,car v);
% RES is the residue and must eventually be reduced to zero.
factor!-trace << printsf res; terpri!*(nil) >>;
if not polyzerop res and
cdr variable!-set and not zerop cdr v then <<
predicted!-forms:=make!-bivariate!-vec!-mod!-p(best!-factors,
cdr variable!-set,car v,number!-of!-factors);
find!-multivariate!-factors!-mod!-p(
make!-bivariate!-mod!-p(poly,cdr variable!-set,car v),
predicted!-forms,list v);
% Answers in PREDICTED!-FORMS.
find!-msg4(predicted!-forms,v);
make!-predicted!-forms(predicted!-forms,car v);
% Sets max!-unknowns and number!-of!-unknowns.
find!-msg5();
unknowns!-count!-list:=number!-of!-unknowns;
while unknowns!-count!-list and
(car (w:=car unknowns!-count!-list))=1 do
begin scalar i,r;
unknowns!-count!-list:=cdr unknowns!-count!-list;
i:=cdr w;
w:=quotient!-mod!-p(poly!-remaining,r:=getv(best!-factors,i));
if didntgo w or
not polyzerop difference!-mod!-p(poly!-remaining,
times!-mod!-p(w,r)) then
if one!-prediction!-failed then <<
factor!-trace printstr "Predictions are no good";
max!-unknowns:=nil >>
else <<
factor!-trace <<
prin2!* "Guess for f(";
prin2!* i;
printstr ") was bad." >>;
one!-prediction!-failed:=i >>
else <<
putv(prediction!-results,i,r);
factor!-trace <<
prin2!* "Prediction for f("; prin2!* i;
prin2!* ") worked: ";
printsf r >>;
poly!-remaining:=w >>
end;
w:=length unknowns!-count!-list;
if w=1 and not one!-prediction!-failed then <<
putv(best!-factors,cdar unknowns!-count!-list,poly!-remaining);
go to exit >>
else if w=0 and one!-prediction!-failed then <<
putv(best!-factors,one!-prediction!-failed,poly!-remaining);
go to exit >>;
solve!-count:=1;
if max!-unknowns then
correction!-vectors:=
make!-correction!-vectors(best!-factors,max!-unknowns) >>;
bhat0s:=make!-multivariate!-hatvec!-mod!-p(b0s,number!-of!-factors);
return multihen1(list(res,
growth!-factor,
first!-time,
bhat0s,
b0s,
variable!-set,
solve!-count,
correction!-vectors,
unknowns!-count!-list,
best!-factors,
v,
degbd,
soln!-matrices,
predicted!-forms,
poly!-remaining,
prediction!-results,
one!-prediction!-failed),
nil);
exit:
multihen!-exit(first!-time,best!-factors,nil);
end) (factor!-level+1);
symbolic procedure multihen1(u,zz);
begin scalar res,test!-prediction,growth!-factor,first!-time,hat0s,
x0s,variable!-set,solve!-count,correction!-vectors,
unknowns!-count!-list,correction!-factor,frvec,v,
degbd,soln!-matrices,predicted!-forms,poly!-remaining,
fvec,previous!-prediction!-holds,
prediction!-results,one!-prediction!-failed,
bool,d,x1,k,kk,substres,w;
res := car u; u := cdr u;
growth!-factor := car u; u := cdr u;
first!-time := car u; u := cdr u;
hat0s := car u; u := cdr u;
x0s := car u; u := cdr u;
variable!-set := car u; u := cdr u;
solve!-count := car u; u := cdr u;
correction!-vectors := car u; u := cdr u;
unknowns!-count!-list := car u; u := cdr u;
frvec := car u; u := cdr u;
v := car u; u := cdr u;
degbd := car u; u := cdr u;
soln!-matrices := car u; u := cdr u;
predicted!-forms := car u; u := cdr u;
poly!-remaining := car u; u := cdr u;
prediction!-results := car u; u := cdr u;
if zz then <<fvec := car u; u := cdr u;
previous!-prediction!-holds := car u; u := cdr u>>;
one!-prediction!-failed := car u;
correction!-factor:=growth!-factor;
% Next power of growth-factor we are adding to the factors.
x1:=mkvect number!-of!-factors;
k:=1;
kk:=0;
temploop:
bool := nil;
while not bool and not polyzerop res and (null max!-unknowns
or null test!-prediction) do
if k>degbd then <<
factor!-trace <<
prin2!* "We have overshot the degree bound for ";
printvar car v >>;
if !*overshoot then
prin2t "Multivariate degree bound overshoot -> restart";
bad!-case:= bool := t >>
else
if polyzerop(substres:=evaluate!-mod!-p(res,car v,cdr v))
then <<
k:=iadd1 k;
res:=diff!-over!-k!-mod!-p(res,k,car v);
correction!-factor:=
times!-mod!-p(correction!-factor,growth!-factor) >>
else begin
multihen!-msg(growth!-factor,first!-time,k,kk,substres,zz);
if null zz
then <<kk := kk#+1; if first!-time then first!-time := nil>>;
solve!-for!-corrections(substres,hat0s,x0s,x1,
cdr variable!-set);
% Answers left in x1.
if bad!-case then return (bool := t);
if max!-unknowns then <<
solve!-count:=iadd1 solve!-count;
for i:=1:number!-of!-factors do
putv(getv(correction!-vectors,i),solve!-count,getv(x1,i));
if solve!-count=caar unknowns!-count!-list then
test!-prediction:=t >>;
if zz then
for i:=1:number!-of!-factors do
putv(frvec,i,plus!-mod!-p(getv(frvec,i),times!-mod!-p(
getv(x1,i),correction!-factor)));
factor!-trace <<
printstr " Giving:";
if null zz then
printvec(" f(",number!-of!-factors,",1) = ",x1)
else <<
printvec(" a(",number!-of!-factors,",1) = ",x1);
printstr " New a's are now:";
printvec(" a(",number!-of!-factors,") = ",frvec) >>>>;
d:=times!-mod!-p(correction!-factor,
if zz then form!-sum!-and!-product!-mod!-p(x1,fhatvec,
number!-of!-factors)
else terms!-done!-mod!-p(frvec,x1,correction!-factor));
if degree!-in!-variable(d,car v)>degbd then <<
factor!-trace <<
prin2!* "We have overshot the degree bound for ";
printvar car v >>;
if !*overshoot then
prin2t "Multivariate degree bound overshoot -> restart";
bad!-case:=t;
return (bool := t)>>;
d:=diff!-k!-times!-mod!-p(d,k,car v);
if null zz then
for i:=1:number!-of!-factors do
putv(frvec,i,
plus!-mod!-p(getv(frvec,i),
times!-mod!-p(getv(x1,i),correction!-factor)));
k:=iadd1 k;
res:=diff!-over!-k!-mod!-p(difference!-mod!-p(res,d),k,car v);
factor!-trace <<
if null zz then <<printstr " New factors are now:";
printvec(" f(",number!-of!-factors,") = ",frvec)>>;
prin2!* " and residue = ";
printsf res;
printstr "-------------"
>>;
correction!-factor:=
times!-mod!-p(correction!-factor,growth!-factor) end;
if not polyzerop res and not bad!-case then <<
if null zz or null soln!-matrices then
soln!-matrices
:= construct!-soln!-matrices(predicted!-forms,cdr v);
factor!-trace <<
if null zz then <<
printstr "We use the results from the Hensel growth to";
printstr "produce a set of linear equations to solve";
printstr "for coefficients in the relevant factors:" >>
else <<
printstr "The Hensel growth so far allows us to test some of";
printstr "our predictions:" >>>>;
bool := nil;
while not bool and unknowns!-count!-list and
(car (w:=car unknowns!-count!-list))=solve!-count do <<
unknowns!-count!-list:=cdr unknowns!-count!-list;
factor!-trace
print!-linear!-system(cdr w,soln!-matrices,
correction!-vectors,predicted!-forms,car v);
w:=try!-prediction(soln!-matrices,correction!-vectors,
predicted!-forms,car w,cdr w,poly!-remaining,car v,
if zz then fvec else nil,
if zz then fhatvec else nil);
if car w='singular or car w='bad!-prediction then
if one!-prediction!-failed then <<
factor!-trace printstr "Predictions were no help.";
max!-unknowns:=nil;
bool := t>>
else if null zz then one!-prediction!-failed:=cdr w
else <<
if previous!-prediction!-holds then <<
predictions:=delasc(car v,predictions);
previous!-prediction!-holds:=nil >>;
one!-prediction!-failed:=cdr w >>
else <<
putv(prediction!-results,car w,cadr w);
poly!-remaining:=caddr w >> >>;
if null max!-unknowns then <<
if zz and previous!-prediction!-holds then
predictions:=delasc(car v,predictions);
goto temploop >>;
w:=length unknowns!-count!-list;
if w>1 or (w=1 and one!-prediction!-failed) then <<
test!-prediction:=nil;
goto temploop >>;
if w=1 or one!-prediction!-failed then <<
w:=if one!-prediction!-failed then one!-prediction!-failed
else cdar unknowns!-count!-list;
putv(prediction!-results,w,
if null zz then poly!-remaining
else quotfail!-mod!-p(poly!-remaining,
getv(fhatvec,w)))>>;
for i:=1:number!-of!-factors do
putv(frvec,i,getv(prediction!-results,i));
if (not previous!-prediction!-holds or null zz)
and not one!-prediction!-failed then
predictions:=
(car v .
list(soln!-matrices,predicted!-forms,max!-unknowns,
number!-of!-unknowns))
. predictions >>;
multihen!-exit(first!-time,frvec,zz)
end;
symbolic
procedure multihen!-msg(growth!-factor,first!-time,k,kk,substres,zz);
factor!-trace <<
prin2!* "Hensel Step "; printstr (kk:=kk #+ 1);
prin2!* "-------------";
if kk>10 then printstr "-" else terpri!*(t);
prin2!* "Next corrections are for (";
prinsf growth!-factor;
if not (k=1) then <<
prin2!* ") ** ";
prin2!* k >> else prin2!* '!);
printstr ". To find these we solve:";
if zz then prin2!* " sum over i [ a(i,1)*fhat(i,0) ] = "
else prin2!* " sum over i [ f(i,1)*fhat(i,0) ] = ";
prinsf substres;
prin2!* " mod ";
prin2!* hensel!-growth!-size;
if zz then printstr " for a(i,1). "
else printstr " for f(i,1), ";
if null zz and first!-time then <<
prin2!*
" where fhat(i,0) = product over j [ f(j,0) ]";
prin2!* " / f(i,0) mod ";
printstr hensel!-growth!-size >>;
terpri!*(nil)
>>;
symbolic procedure multihen!-exit(first!-time,frvec,zz);
factor!-trace <<
if not bad!-case then
if first!-time then
if zz then printstr "But these a's are already correct."
else printstr "Therefore these factors are already correct."
else <<
if zz then <<printstr "Correct a's are:";
printvec(" a(",number!-of!-factors,") = ",frvec)>>
else <<printstr "Correct factors are:";
printvec(" f(",number!-of!-factors,") = ",frvec) >>>>;
terpri!*(nil);
printstr "**************************************************";
terpri!*(nil)>>;
symbolic procedure find!-msg1(best!-factors,growth!-factor,poly);
factor!-trace <<
printstr "Want f(i) s.t.";
prin2!* " product over i [ f(i) ] = ";
prinsf poly;
prin2!* " mod ";
printstr hensel!-growth!-size;
terpri!*(nil);
printstr "We know f(i) as follows:";
printvec(" f(",number!-of!-factors,") = ",best!-factors);
prin2!* " and we shall put in powers of ";
prinsf growth!-factor;
printstr " to find them fully."
>>;
symbolic procedure find!-msg2(v,variable!-set);
factor!-trace <<
prin2!*
"First solve the problem in one less variable by putting ";
prinvar car v; prin2!* "="; printstr cdr v;
if cdr variable!-set then <<
prin2!* "and growing wrt ";
printvar caadr variable!-set
>>;
terpri!*(nil)
>>;
symbolic procedure find!-msg3(best!-factors,v);
factor!-trace <<
prin2!* "After putting back any knowledge of ";
prinvar car v;
printstr ", we have the";
printstr "factors so far as:";
printvec(" f(",number!-of!-factors,") = ",best!-factors);
printstr "Subtracting the product of these from the polynomial";
prin2!* "and differentiating wrt "; prinvar car v;
printstr " gives a residue:"
>>;
symbolic procedure find!-msg4(predicted!-forms,v);
factor!-trace <<
printstr "To help reduce the number of Hensel steps we try";
prin2!* "predicting how many terms each factor will have wrt ";
prinvar car v; printstr ".";
printstr
"Predictions are based on the bivariate factors :";
printvec(" f(",number!-of!-factors,") = ",predicted!-forms)
>>;
symbolic procedure find!-msg5;
factor!-trace <<
terpri!*(nil);
printstr "We predict :";
for each w in number!-of!-unknowns do <<
prin2!* car w;
prin2!* " terms in f("; prin2!* cdr w; printstr '!) >>;
if (caar number!-of!-unknowns)=1 then <<
prin2!* "Since we predict only one term for f(";
prin2!* cdar number!-of!-unknowns;
printstr "), we can try";
printstr "dividing it out now:" >>
else <<
prin2!* "So we shall do at least ";
prin2!* isub1 caar number!-of!-unknowns;
prin2!* " Hensel step";
if (caar number!-of!-unknowns)=2 then printstr "."
else printstr "s." >>;
terpri!*(nil) >>;
symbolic procedure solve!-for!-corrections(c,fhatvec,fvec,resvec,vset);
% ....;
if null vset then
for i:=1:number!-of!-factors do
putv(resvec,i,
remainder!-mod!-p(
times!-mod!-p(c,getv(alphavec,i)),
getv(fvec,i)))
else (lambda factor!-level; begin
scalar residue,growth!-factor,f0s,fhat0s,v,
degbd,first!-time,redc,
predicted!-forms,max!-unknowns,solve!-count,number!-of!-unknowns,
correction!-vectors,soln!-matrices,w,previous!-prediction!-holds,
unknowns!-count!-list,poly!-remaining,
prediction!-results,one!-prediction!-failed;
v:=car vset;
degbd:=get!-degree!-bound car v;
first!-time:=t;
growth!-factor:=make!-growth!-factor v;
poly!-remaining:=c;
prediction!-results:=mkvect number!-of!-factors;
redc:=evaluate!-mod!-p(c,car v,cdr v);
solve!-msg1(c,fvec,v);
solve!-for!-corrections(redc,
fhat0s:=reduce!-vec!-by!-one!-var!-mod!-p(
fhatvec,v,number!-of!-factors),
f0s:=reduce!-vec!-by!-one!-var!-mod!-p(
fvec,v,number!-of!-factors),
resvec,
cdr vset); % Results left in RESVEC.
if bad!-case then return;
solve!-msg2(resvec,v);
residue:=diff!-over!-k!-mod!-p(difference!-mod!-p(c,
form!-sum!-and!-product!-mod!-p(resvec,fhatvec,
number!-of!-factors)),1,car v);
factor!-trace <<
printsf residue;
prin2!* " Now we shall put in the powers of ";
prinsf growth!-factor;
printstr " to find the a's fully."
>>;
if not polyzerop residue and not zerop cdr v then <<
w:=atsoc(car v,predictions);
if w then <<
previous!-prediction!-holds:=t;
factor!-trace <<
printstr
"We shall use the previous prediction for the form of";
prin2!* "polynomials wrt "; printvar car v >>;
w:=cdr w;
soln!-matrices:=car w;
predicted!-forms:=cadr w;
max!-unknowns:=caddr w;
number!-of!-unknowns:=cadr cddr w >>
else <<
factor!-trace <<
printstr
"We shall use a new prediction for the form of polynomials ";
prin2!* "wrt "; printvar car v >>;
predicted!-forms:=mkvect number!-of!-factors;
for i:=1:number!-of!-factors do
putv(predicted!-forms,i,getv(fvec,i));
% Make a copy of the factors in a vector we shall overwrite.
make!-predicted!-forms(predicted!-forms,car v);
% Sets max!-unknowns and number!-of!-unknowns.
>>;
solve!-msg3();
unknowns!-count!-list:=number!-of!-unknowns;
while unknowns!-count!-list and
(car (w:=car unknowns!-count!-list))=1 do
begin scalar i,r,wr,fi;
unknowns!-count!-list:=cdr unknowns!-count!-list;
i:=cdr w;
w:=quotient!-mod!-p(
wr:=difference!-mod!-p(poly!-remaining,
times!-mod!-p(r:=getv(resvec,i),getv(fhatvec,i))),
fi:=getv(fvec,i));
if didntgo w or not polyzerop
difference!-mod!-p(wr,times!-mod!-p(w,fi)) then
if one!-prediction!-failed then <<
factor!-trace printstr "Predictions are no good.";
max!-unknowns:=nil >>
else <<
factor!-trace <<
prin2!* "Guess for a(";
prin2!* i;
printstr ") was bad." >>;
one!-prediction!-failed:=i >>
else <<
putv(prediction!-results,i,r);
factor!-trace <<
prin2!* "Prediction for a("; prin2!* i;
prin2!* ") worked: ";
printsf r >>;
poly!-remaining:=wr >>
end;
w:=length unknowns!-count!-list;
if w=1 and not one!-prediction!-failed then <<
putv(resvec,cdar unknowns!-count!-list,
quotfail!-mod!-p(poly!-remaining,getv(fhatvec,
cdar unknowns!-count!-list)));
go to exit >>
else if w=0 and one!-prediction!-failed and max!-unknowns then <<
putv(resvec,one!-prediction!-failed,
quotfail!-mod!-p(poly!-remaining,getv(fhatvec,
one!-prediction!-failed)));
go to exit >>;
solve!-count:=1;
if max!-unknowns then
correction!-vectors:=
make!-correction!-vectors(resvec,max!-unknowns) >>;
if not polyzerop residue then first!-time:=nil;
return multihen1(list(residue,
growth!-factor,
first!-time,
fhat0s,
f0s,
vset,
solve!-count,
correction!-vectors,
unknowns!-count!-list,
resvec,
v,
degbd,
soln!-matrices,
predicted!-forms,
poly!-remaining,
prediction!-results,
fvec,
previous!-prediction!-holds,
one!-prediction!-failed),
t);
exit:
multihen!-exit(first!-time,resvec,t);
end) (factor!-level+1);
symbolic procedure solve!-msg1(c,fvec,v);
factor!-trace <<
printstr "Want a(i) s.t.";
prin2!* "(*) sum over i [ a(i)*fhat(i) ] = ";
prinsf c;
prin2!* " mod ";
printstr hensel!-growth!-size;
prin2!* " where fhat(i) = product over j [ f(j) ]";
prin2!* " / f(i) mod ";
printstr hensel!-growth!-size;
printstr " and";
printvec(" f(",number!-of!-factors,") = ",fvec);
terpri!*(nil);
prin2!*
"First solve the problem in one less variable by putting ";
prinvar car v; prin2!* '!=; printstr cdr v;
terpri!*(nil)
>>;
symbolic procedure solve!-msg2(resvec,v);
factor!-trace <<
printstr "Giving:";
printvec(" a(",number!-of!-factors,",0) = ",resvec);
printstr "Subtracting the contributions these give in (*) from";
prin2!* "the R.H.S. of (*) ";
prin2!* "and differentiating wrt "; prinvar car v;
printstr " gives a residue:"
>>;
symbolic procedure solve!-msg3;
factor!-trace <<
terpri!*(nil);
printstr "We predict :";
for each w in number!-of!-unknowns do <<
prin2!* car w;
prin2!* " terms in a("; prin2!* cdr w; printstr '!) >>;
if (caar number!-of!-unknowns)=1 then <<
prin2!* "Since we predict only one term for a(";
prin2!* cdar number!-of!-unknowns;
printstr "), we can test it right away:" >>
else <<
prin2!* "So we shall do at least ";
prin2!* isub1 caar number!-of!-unknowns;
prin2!* " Hensel step";
if (caar number!-of!-unknowns)=2 then printstr "."
else printstr "s." >>;
terpri!*(nil) >>;
endmodule;
end;