coeffts.red 9.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260
  1. module coeffts;
  2. % Authors: A. C. Norman and P. M. A. Moore, 1981.
  3. fluid '(!*trfac
  4. alphalist
  5. best!-known!-factor!-list
  6. best!-known!-factors
  7. coefft!-vectors
  8. deg!-of!-unknown
  9. difference!-for!-unknown
  10. divisor!-for!-unknown
  11. factor!-level
  12. factor!-trace!-list
  13. full!-gcd
  14. hensel!-growth!-size
  15. image!-factors
  16. m!-image!-variable
  17. multivariate!-factors
  18. multivariate!-input!-poly
  19. non!-monic
  20. number!-of!-factors
  21. polyzero
  22. reconstructing!-gcd
  23. true!-leading!-coeffts
  24. unknown
  25. unknowns!-list);
  26. %**********************************************************************;
  27. % Code for trying to determine more multivariate coefficients
  28. % by inspection before using multivariate hensel construction.
  29. symbolic procedure determine!-more!-coeffts();
  30. % ...
  31. begin scalar unknowns!-list,uv,r,w,best!-known!-factor!-list;
  32. best!-known!-factors:=mkvect number!-of!-factors;
  33. uv:=mkvect number!-of!-factors;
  34. for i:=number!-of!-factors step -1 until 1 do
  35. putv(uv,i,convert!-factor!-to!-termvector(
  36. getv(image!-factors,i),getv(true!-leading!-coeffts,i)));
  37. r:=red multivariate!-input!-poly;
  38. % we know all about the leading coeffts;
  39. if not depends!-on!-var(r,m!-image!-variable)
  40. or null(w:=try!-first!-coefft(
  41. ldeg r,lc r,unknowns!-list,uv)) then <<
  42. for i:=1:number!-of!-factors do
  43. putv(best!-known!-factors,i,force!-lc(
  44. getv(image!-factors,i),getv(true!-leading!-coeffts,i)));
  45. coefft!-vectors:=uv;
  46. return nil >>;
  47. factor!-trace <<
  48. printstr
  49. "By exploiting any sparsity wrt the main variable in the";
  50. printstr "factors, we can try guessing some of the multivariate";
  51. printstr "coefficients." >>;
  52. try!-other!-coeffts(r,unknowns!-list,uv);
  53. w:=convert!-and!-trial!-divide uv;
  54. % trace!-time
  55. % if full!-gcd then prin2t "Possible gcd found"
  56. % else prin2t "Have found some coefficients";
  57. return set!-up!-globals(uv,w)
  58. end;
  59. symbolic procedure convert!-factor!-to!-termvector(u,tlc);
  60. % ...
  61. begin scalar termlist,res,n,slist;
  62. termlist:=(ldeg u . tlc) . list!-terms!-in!-factor red u;
  63. res:=mkvect (n:=length termlist);
  64. for i:=1:n do <<
  65. slist:=(caar termlist . i) . slist;
  66. putv(res,i,car termlist);
  67. termlist:=cdr termlist >>;
  68. putv(res,0,(n . (n #- 1)));
  69. unknowns!-list:=(reversip slist) . unknowns!-list;
  70. return res
  71. end;
  72. symbolic procedure try!-first!-coefft(n,c,slist,uv);
  73. % ...
  74. begin scalar combns,unknown,w,l,d,v,m;
  75. combns:=get!-term(n,slist);
  76. if (combns='no) or not null cdr combns then return nil;
  77. l:=car combns;
  78. for i:=1:number!-of!-factors do <<
  79. w:=getv(getv(uv,i),car l); % degree . coefft ;
  80. if null cdr w then <<
  81. if unknown then <<c := nil; i := number!-of!-factors + 1>>
  82. else <<unknown := i . car l; d := car w>>>>
  83. else <<
  84. c:=quotf(c,cdr w);
  85. if didntgo c then i := number!-of!-factors+1>>;
  86. l:=cdr l >>;
  87. if didntgo c then return nil;
  88. putv(v:=getv(uv,car unknown),cdr unknown,(d . c));
  89. m:=getv(v,0);
  90. putv(v,0,(car m . (cdr m #- 1)));
  91. if cdr m = 1 and factors!-complete uv then return 'complete;
  92. return c
  93. end;
  94. symbolic procedure solve!-next!-coefft(n,c,slist,uv);
  95. % ...
  96. begin scalar combns,w,unknown,deg!-of!-unknown,divisor!-for!-unknown,
  97. difference!-for!-unknown,v;
  98. difference!-for!-unknown:=polyzero;
  99. divisor!-for!-unknown:=polyzero;
  100. combns:=get!-term(n,slist);
  101. if combns='no then return 'nogood;
  102. while combns do <<
  103. w:=split!-term!-list(car combns,uv);
  104. if w='nogood then combns := nil else combns:=cdr combns >>;
  105. if w='nogood then return w;
  106. if null unknown then return;
  107. w:=quotf(addf(c,negf difference!-for!-unknown),
  108. divisor!-for!-unknown);
  109. if didntgo w then return 'nogood;
  110. putv(v:=getv(uv,car unknown),cdr unknown,(deg!-of!-unknown . w));
  111. n:=getv(v,0);
  112. putv(v,0,(car n . (cdr n #- 1)));
  113. if cdr n = 1 and factors!-complete uv then return 'complete;
  114. return w
  115. end;
  116. symbolic procedure split!-term!-list(term!-combn,uv);
  117. % ...
  118. begin scalar a,v,w;
  119. a:=1;
  120. for i:=1:number!-of!-factors do <<
  121. w:=getv(getv(uv,i),car term!-combn); % degree . coefft ;
  122. if null cdr w then
  123. if v or (unknown and not((i.car term!-combn)=unknown)) then
  124. <<v:='nogood; i := number!-of!-factors+1>>
  125. else <<
  126. unknown:=(i . car term!-combn);
  127. deg!-of!-unknown:=car w;
  128. v:=unknown >>
  129. else a:=multf(a,cdr w);
  130. if not(v eq 'nogood) then term!-combn:=cdr term!-combn >>;
  131. if v='nogood then return v;
  132. if v then divisor!-for!-unknown:=addf(divisor!-for!-unknown,a)
  133. else difference!-for!-unknown:=addf(difference!-for!-unknown,a);
  134. return 'ok
  135. end;
  136. symbolic procedure factors!-complete uv;
  137. % ...
  138. begin scalar factor!-not!-done,r;
  139. r:=t;
  140. for i:=1:number!-of!-factors do
  141. if not(cdr getv(getv(uv,i),0)=0) then
  142. if factor!-not!-done then <<r:=nil; i:=number!-of!-factors+1>>
  143. else factor!-not!-done:=t;
  144. return r
  145. end;
  146. symbolic procedure convert!-and!-trial!-divide uv;
  147. % ...
  148. begin scalar w,r,fdone!-product!-mod!-p,om;
  149. om:=set!-modulus hensel!-growth!-size;
  150. fdone!-product!-mod!-p:=1;
  151. for i:=1:number!-of!-factors do <<
  152. w:=getv(uv,i);
  153. w:= if (cdr getv(w,0))=0 then termvector2sf w
  154. else merge!-terms(getv(image!-factors,i),w);
  155. r:=quotf(multivariate!-input!-poly,w);
  156. if didntgo r then best!-known!-factor!-list:=
  157. ((i . w) . best!-known!-factor!-list)
  158. else if reconstructing!-gcd and i=1
  159. then <<full!-gcd:=if non!-monic then car primitive!.parts(
  160. list w,m!-image!-variable,nil) else w;
  161. i := number!-of!-factors+1>>
  162. else <<
  163. multivariate!-factors:=w . multivariate!-factors;
  164. fdone!-product!-mod!-p:=times!-mod!-p(
  165. reduce!-mod!-p getv(image!-factors,i),
  166. fdone!-product!-mod!-p);
  167. multivariate!-input!-poly:=r >> >>;
  168. if full!-gcd then return;
  169. if null best!-known!-factor!-list then multivariate!-factors:=
  170. primitive!.parts(multivariate!-factors,m!-image!-variable,nil)
  171. else if null cdr best!-known!-factor!-list then <<
  172. if reconstructing!-gcd then
  173. if not(caar best!-known!-factor!-list=1) then
  174. errorf("gcd is jiggered in determining other coeffts")
  175. else full!-gcd:=if non!-monic then car primitive!.parts(
  176. list multivariate!-input!-poly,
  177. m!-image!-variable,nil)
  178. else multivariate!-input!-poly
  179. else multivariate!-factors:=primitive!.parts(
  180. multivariate!-input!-poly . multivariate!-factors,
  181. m!-image!-variable,nil);
  182. best!-known!-factor!-list:=nil >>;
  183. factor!-trace <<
  184. if null best!-known!-factor!-list then
  185. printstr
  186. "We have completely determined all the factors this way"
  187. else if multivariate!-factors then <<
  188. prin2!* "We have completely determined the following factor";
  189. printstr if (length multivariate!-factors)=1 then ":" else "s:";
  190. for each ww in multivariate!-factors do printsf ww >> >>;
  191. set!-modulus om;
  192. return fdone!-product!-mod!-p
  193. end;
  194. symbolic procedure set!-up!-globals(uv,f!-product);
  195. if null best!-known!-factor!-list or full!-gcd then 'done
  196. else begin scalar i,r,n,k,flist!-mod!-p,imf,om,savek;
  197. n:=length best!-known!-factor!-list;
  198. best!-known!-factors:=mkvect n;
  199. coefft!-vectors:=mkvect n;
  200. r:=mkvect n;
  201. k:=if reconstructing!-gcd then 1 else 0;
  202. om:=set!-modulus hensel!-growth!-size;
  203. for each w in best!-known!-factor!-list do <<
  204. i:=car w; w:=cdr w;
  205. if reconstructing!-gcd and i=1 then << savek:=k; k:=1 >>
  206. else k:=k #+ 1;
  207. % in case we are reconstructing gcd we had better know
  208. % which is the gcd and which the cofactor - so don't move
  209. % move the gcd from elt one;
  210. putv(r,k,imf:=getv(image!-factors,i));
  211. flist!-mod!-p:=(reduce!-mod!-p imf) . flist!-mod!-p;
  212. putv(best!-known!-factors,k,w);
  213. putv(coefft!-vectors,k,getv(uv,i));
  214. if reconstructing!-gcd and k=1 then k:=savek;
  215. % restore k if necessary;
  216. >>;
  217. if not(n=number!-of!-factors) then <<
  218. alphalist:=for each modf in flist!-mod!-p collect
  219. (modf . remainder!-mod!-p(times!-mod!-p(f!-product,
  220. cdr get!-alpha modf),modf));
  221. number!-of!-factors:=n >>;
  222. set!-modulus om;
  223. image!-factors:=r;
  224. return 'need! to! reconstruct
  225. end;
  226. symbolic procedure get!-term(n,l);
  227. % ...
  228. if n#<0 then 'no
  229. else if null cdr l then get!-term!-n(n,car l)
  230. else begin scalar w,res;
  231. for each fterm in car l do <<
  232. w:=get!-term(n#-car fterm,cdr l);
  233. if not(w='no) then res:=
  234. append(for each v in w collect (cdr fterm . v),res) >>;
  235. return if null res then 'no else res
  236. end;
  237. symbolic procedure get!-term!-n(n,u);
  238. if null u or n #> caar u then 'no
  239. else if caar u = n then list(cdar u . nil)
  240. else get!-term!-n(n,cdr u);
  241. endmodule;
  242. end;