FORM  4.3
polyfact.cc
Go to the documentation of this file.
1 
5 /* #[ License : */
6 /*
7  * Copyright (C) 1984-2022 J.A.M. Vermaseren
8  * When using this file you are requested to refer to the publication
9  * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
10  * This is considered a matter of courtesy as the development was paid
11  * for by FOM the Dutch physics granting agency and we would like to
12  * be able to track its scientific use to convince FOM of its value
13  * for the community.
14  *
15  * This file is part of FORM.
16  *
17  * FORM is free software: you can redistribute it and/or modify it under the
18  * terms of the GNU General Public License as published by the Free Software
19  * Foundation, either version 3 of the License, or (at your option) any later
20  * version.
21  *
22  * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
23  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
24  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
25  * details.
26  *
27  * You should have received a copy of the GNU General Public License along
28  * with FORM. If not, see <http://www.gnu.org/licenses/>.
29  */
30 /* #] License : */
31 /*
32  #[ include :
33 */
34 
35 #include "poly.h"
36 #include "polygcd.h"
37 #include "polyfact.h"
38 
39 #include <cmath>
40 #include <vector>
41 #include <iostream>
42 #include <algorithm>
43 #include <climits>
44 
45 //#define DEBUG
46 
47 #ifdef DEBUG
48 #include "mytime.h"
49 #endif
50 
51 using namespace std;
52 
53 /*
54  #] include :
55  #[ tostring :
56 */
57 
58 // Turns a factorized_poly into a readable string
59 const string factorized_poly::tostring () const {
60 
61  // empty
62  if (factor.size()==0)
63  return "no_factors";
64 
65  string res;
66 
67  // polynomial
68  for (int i=0; i<(int)factor.size(); i++) {
69  if (i>0) res += "*";
70  res += "(";
71  res += poly(factor[i],0,1).to_string();
72  res += ")";
73  if (power[i]>1) {
74  res += "^";
75  char tmp[100];
76  snprintf (tmp,100,"%i",power[i]);
77  res += tmp;
78  }
79  }
80 
81  // modulo p^n
82  if (factor[0].modp>0) {
83  res += " (mod ";
84  char tmp[12];
85  snprintf (tmp,12,"%i",factor[0].modp);
86  res += tmp;
87  if (factor[0].modn > 1) {
88  snprintf (tmp,12,"%i",factor[0].modn);
89  res += "^";
90  res += tmp;
91  }
92  res += ")";
93  }
94 
95  return res;
96 }
97 
98 /*
99  #] tostring :
100  #[ ostream operator :
101 */
102 
103 // ostream operator for outputting a factorized_poly
104 ostream& operator<< (ostream &out, const factorized_poly &a) {
105  return out << a.tostring();
106 }
107 
108 // ostream operator for outputting a vector<T>
109 template<class T> ostream& operator<< (ostream &out, const vector<T> &v) {
110  out<<"{";
111  for (int i=0; i<(int)v.size(); i++) {
112  if (i>0) out<<",";
113  out<<v[i];
114  }
115  out<<"}";
116  return out;
117 }
118 
119 /*
120  #] ostream operator :
121  #[ add_factor :
122 */
123 
124 // adds a factor f^p to a factorization
125 void factorized_poly::add_factor(const poly &f, int p) {
126  factor.push_back(f);
127  power.push_back(p);
128 }
129 
130 /*
131  #] add_factor :
132  #[ extended_gcd_Euclidean_lifted :
133 */
134 
152 const vector<poly> polyfact::extended_gcd_Euclidean_lifted (const poly &a, const poly &b) {
153 
154 #ifdef DEBUGALL
155  cout << "*** [" << thetime() << "] CALL: extended_Euclidean_lifted("<<a<<","<<b<<")"<<endl;
156 #endif
157 
158  POLY_GETIDENTITY(a);
159 
160  // Calculate s,t,gcd (mod p) with the Extended Euclidean Algorithm.
161  poly s(a,a.modp,1);
162  poly t(b,b.modp,1);
163  poly sa(BHEAD 1,a.modp,1);
164  poly sb(BHEAD 0,b.modp,1);
165  poly ta(BHEAD 0,a.modp,1);
166  poly tb(BHEAD 1,b.modp,1);
167 
168  while (!t.is_zero()) {
169  poly x(s/t);
170  swap(s -=x*t , t);
171  swap(sa-=x*ta, ta);
172  swap(sb-=x*tb, tb);
173  }
174 
175  // Normalize the result.
176  sa /= s.integer_lcoeff();
177  sb /= s.integer_lcoeff();
178 
179  // Lift the result to modolu p^n with p-adic Newton's iteration.
180  poly samodp(sa);
181  poly sbmodp(sb);
182  poly term(BHEAD 1);
183 
184  sa.setmod(0,1);
185  sb.setmod(0,1);
186 
187  poly amodp(a,a.modp,1);
188  poly bmodp(b,a.modp,1);
189  poly error(poly(BHEAD 1) - sa*a - sb*b);
190  poly p(BHEAD a.modp);
191 
192  for (int n=1; n<a.modn && !error.is_zero(); n++) {
193  error /= p;
194  term *= p;
195  poly errormodp(error,a.modp,1);
196  poly dsa((samodp * errormodp) % bmodp);
197  // equivalent, but faster than the symmetric
198  // poly dsb((sbmodp * errormodp) % amodp);
199  poly dsb((errormodp - dsa*amodp) / bmodp);
200  sa += term * dsa;
201  sb += term * dsb;
202  error -= a*dsa + b*dsb;
203  }
204 
205  sa.setmod(a.modp,a.modn);
206  sb.setmod(a.modp,a.modn);
207 
208  // Output the result
209  vector<poly> res;
210  res.push_back(sa);
211  res.push_back(sb);
212 
213 #ifdef DEBUGALL
214  cout << "*** [" << thetime() << "] RES : extended_Euclidean_lifted("<<a<<","<<b<<") = "<<res<<endl;
215 #endif
216 
217  return res;
218 }
219 
220 /*
221  #] extended_gcd_Euclidean_lifted :
222  #[ solve_Diophantine_univariate :
223 */
224 
256 const vector<poly> polyfact::solve_Diophantine_univariate (const vector<poly> &a, const poly &b) {
257 
258 #ifdef DEBUGALL
259  cout << "*** [" << thetime() << "] CALL: solve_Diophantine_univariate(" <<a<<","<<b<<")"<<endl;
260 #endif
261 
262  POLY_GETIDENTITY(b);
263 
264  vector<poly> s(1,b);
265 
266  for (int i=0; i+1<(int)a.size(); i++) {
267  poly A(BHEAD 1,b.modp,b.modn);
268  for (int j=i+1; j<(int)a.size(); j++) A *= a[j];
269 
270  vector<poly> t(extended_gcd_Euclidean_lifted(a[i],A));
271  poly prev(s.back());
272  s.back() = t[1] * prev % a[i];
273  s.push_back(t[0] * prev % A);
274  }
275 
276 #ifdef DEBUGALL
277  cout << "*** [" << thetime() << "] RES : solve_Diophantine_univariate(" <<a<<","<<b<<") = "<<s<<endl;
278 #endif
279 
280  return s;
281 }
282 
283 /*
284  #] solve_Diophantine_univariate :
285  #[ solve_Diophantine_multivariate :
286 */
287 
319 const vector<poly> polyfact::solve_Diophantine_multivariate (const vector<poly> &a, const poly &b, const vector<int> &x, const vector<int> &c, int d) {
320 
321 #ifdef DEBUGALL
322  cout << "*** [" << thetime() << "] CALL: solve_Diophantine_multivariate(" <<a<<","<<b<<","<<x<<","<<c<<","<<d<<")"<<endl;
323 #endif
324 
325  POLY_GETIDENTITY(b);
326 
327  if (b.is_zero()) return vector<poly>(a.size(),poly(BHEAD 0));
328 
329  if (x.size() == 1) return solve_Diophantine_univariate(a,b);
330 
331  // Reduce the polynomial with the homomorphism <xm-c{m-1}>
332  poly simple(poly::simple_poly(BHEAD x.back(),c.back()));
333 
334  vector<poly> ared (a);
335  for (int i=0; i<(int)ared.size(); i++)
336  ared[i] %= simple;
337 
338  poly bred(b % simple);
339  vector<int> xred(x.begin(),x.end()-1);
340  vector<int> cred(c.begin(),c.end()-1);
341 
342  // Solve the equation in one less variable
343  vector<poly> s(solve_Diophantine_multivariate(ared,bred,xred,cred,d));
344  if (s == vector<poly>()) return vector<poly>();
345 
346  // Cache the Ai = product(aj | j!=i).
347  vector<poly> A(a.size(), poly(BHEAD 1,b.modp,b.modn));
348  for (int i=0; i<(int)a.size(); i++)
349  for (int j=0; j<(int)a.size(); j++)
350  if (i!=j) A[i] *= a[j];
351 
352  // Add the powers (xm-c{m-1})^k with ideal-adic Newton iteration.
353  poly term(BHEAD 1,b.modp,b.modn);
354 
355  poly error(b);
356  for (int i=0; i<(int)A.size(); i++)
357  error -= s[i] * A[i];
358 
359  for (int deg=1; deg<=d; deg++) {
360 
361  if (error.is_zero()) break;
362 
363  error /= simple;
364  term *= simple;
365 
366  vector<poly> ds(solve_Diophantine_multivariate(ared, error%simple, xred, cred, d));
367  if (ds == vector<poly>()) return vector<poly>();
368 
369  for (int i=0; i<(int)s.size(); i++) {
370  s[i] += ds[i] * term;
371  error -= ds[i] * A[i];
372  }
373  }
374 
375  if (!error.is_zero()) return vector<poly>();
376 
377 #ifdef DEBUGALL
378  cout << "*** [" << thetime() << "] RES : solve_Diophantine_multivariate(" <<a<<","<<b<<","<<x<<","<<c<<","<<d<<") = "<<s<<endl;
379 #endif
380 
381  return s;
382 }
383 
384 /*
385  #] solve_Diophantine_multivariate :
386  #[ lift_coefficients :
387 */
388 
421 const vector<poly> polyfact::lift_coefficients (const poly &_A, const vector<poly> &_a) {
422 
423 #ifdef DEBUG
424  cout << "*** [" << thetime() << "] CALL: lift_coefficients("<<_A<<","<<_a<<")"<<endl;
425 #endif
426 
427  POLY_GETIDENTITY(_A);
428 
429  poly A(_A);
430  vector<poly> a(_a);
431  poly term(BHEAD 1);
432 
433  int x = A.first_variable();
434 
435  // Replace the leading term of all factors with lterm(A) mod p
436  poly lead(A.integer_lcoeff());
437  for (int i=0; i<(int)a.size(); i++) {
438  a[i] *= lead / a[i].integer_lcoeff();
439  if (i>0) A*=lead;
440  }
441 
442  // Solve Diophantine equation
443  vector<poly> s(solve_Diophantine_univariate(a,poly(BHEAD 1,A.modp,1)));
444 
445  // Replace the leading term of all factors with lterm(A) mod p^n
446  for (int i=0; i<(int)a.size(); i++) {
447  a[i].setmod(A.modp,A.modn);
448  a[i] += (lead - a[i].integer_lcoeff()) * poly::simple_poly(BHEAD x,0,a[i].degree(x));
449  }
450 
451  // Calculate the error, express it in terms of ai and add corrections.
452  for (int k=2; k<=A.modn; k++) {
453  term *= poly(BHEAD A.modp);
454 
455  poly error(BHEAD -1);
456  for (int i=0; i<(int)a.size(); i++) error *= a[i];
457  error += A;
458  if (error.is_zero()) break;
459 
460  error /= term;
461  error.setmod(A.modp,1);
462 
463  for (int i=0; i<(int)a.size(); i++)
464  a[i] += term * (error * s[i] % a[i]);
465  }
466 
467  // Fix leading coefficients by dividing out integer contents.
468  for (int i=0; i<(int)a.size(); i++)
469  a[i] /= polygcd::integer_content(poly(a[i],0,1));
470 
471 #ifdef DEBUG
472  cout << "*** [" << thetime() << "] RES : lift_coefficients("<<_A<<","<<_a<<") = "<<a<<endl;
473 #endif
474 
475  return a;
476 }
477 
478 /*
479  #] lift_coefficients :
480  #[ predetermine :
481 */
482 
497 void polyfact::predetermine (int dep, const vector<vector<int> > &state, vector<vector<vector<int> > > &terms, vector<int> &term, int sumdeg) {
498  // store the term
499  if (dep == (int)state.size()) {
500  terms[sumdeg].push_back(term);
501  return;
502  }
503 
504  // recursively create new terms
505  term.push_back(0);
506 
507  for (int deg=0; sumdeg+deg<(int)state[dep].size(); deg++)
508  if (state[dep][deg] > 0) {
509  term.back() = deg;
510  predetermine(dep+1, state, terms, term, sumdeg+deg);
511  }
512 
513  term.pop_back();
514 }
515 
516 /*
517  #] predetermine :
518  #[ lift_variables :
519 */
520 
558 const vector<poly> polyfact::lift_variables (const poly &A, const vector<poly> &_a, const vector<int> &x, const vector<int> &c, const vector<poly> &lc) {
559 
560 #ifdef DEBUG
561  cout << "*** [" << thetime() << "] CALL: lift_variables("<<A<<","<<_a<<","<<x<<","<<c<<","<<lc<<")\n";
562 #endif
563 
564  // If univariate, don't lift
565  if (x.size()<=1) return _a;
566 
567  POLY_GETIDENTITY(A);
568 
569  vector<poly> a(_a);
570 
571  // First method: predetermine coefficients
572 
573  // check feasibility, otherwise it tries too many possibilities
574  int cnt = POLYFACT_MAX_PREDETERMINATION;
575  for (int i=0; i<(int)a.size(); i++) {
576  if (a[i].number_of_terms() == 0) return vector<poly>();
577  cnt /= a[i].number_of_terms();
578  }
579 
580  if (cnt>0) {
581  // state[n][d]: coefficient of x^d in a[n] is
582  // 0: non-existent, 1: undetermined, 2: determined
583  int D = A.degree(x[0]);
584  vector<vector<int> > state(a.size(), vector<int>(D+1, 0));
585 
586  for (int i=0; i<(int)a.size(); i++)
587  for (int j=1; j<a[i][0]; j+=a[i][j])
588  state[i][a[i][j+1+x[0]]] = j==1 ? 2 : 1;
589 
590  // collect all products of terms
591  vector<vector<vector<int> > > terms(D+1);
592  vector<int> term;
593  predetermine(0,state,terms,term);
594 
595  // count the number of undetermined coefficients
596  vector<int> num_undet(terms.size(),0);
597  for (int i=0; i<(int)terms.size(); i++)
598  for (int j=0; j<(int)terms[i].size(); j++)
599  for (int k=0; k<(int)terms[i][j].size(); k++)
600  if (state[k][terms[i][j][k]] == 1) num_undet[i]++;
601 
602  // replace the current leading coefficients by the correct ones
603  for (int i=0; i<(int)a.size(); i++)
604  a[i] += (lc[i] - a[i].lcoeff_univar(x[0])) * poly::simple_poly(BHEAD x[0],0,a[i].degree(x[0]));
605 
606  bool changed;
607  do {
608  changed = false;
609 
610  for (int i=0; i<(int)terms.size(); i++) {
611  // is only one coefficient in a equation is undetermined, solve
612  // the equation to determine this coefficient
613  if (num_undet[i] == 1) {
614  // generate equation
615  poly lhs(BHEAD 0), rhs(A.coefficient(x[0],i), A.modp, A.modn);
616  int which_idx=-1, which_deg=-1;
617  for (int j=0; j<(int)terms[i].size(); j++) {
618  poly coeff(BHEAD 1, A.modp, A.modn);
619  bool undet=false;
620  for (int k=0; k<(int)terms[i][j].size(); k++) {
621  if (state[k][terms[i][j][k]] == 1) {
622  undet = true;
623  which_idx=k;
624  which_deg=terms[i][j][k];
625  }
626  else
627  coeff *= a[k].coefficient(x[0], terms[i][j][k]);
628  }
629  if (undet)
630  lhs = coeff;
631  else
632  rhs -= coeff;
633  }
634  // solve equation
635  if (A.modn > 1) rhs.setmod(0,1);
636  if (lhs.is_zero() || !(rhs%lhs).is_zero()) return vector<poly>();
637  a[which_idx] += (rhs / lhs - a[which_idx].coefficient(x[0],which_deg)) * poly::simple_poly(BHEAD x[0],0,which_deg);
638  state[which_idx][which_deg] = 2;
639 
640  // update number of undetermined coefficients
641  for (int j=0; j<(int)terms.size(); j++)
642  for (int k=0; k<(int)terms[j].size(); k++)
643  if (terms[j][k][which_idx] == which_deg)
644  num_undet[j]--;
645 
646  changed = true;
647  }
648  }
649  }
650  while (changed);
651 
652  // if this is the complete result, skip lifting
653  poly check(BHEAD 1, A.modn>1?0:A.modp, 1);
654  for (int i=0; i<(int)a.size(); i++)
655  check *= a[i];
656 
657  if (check == A) return a;
658  }
659 
660  // Second method: Hensel lifting
661 
662  // Calculate A and lc's modulo Ii = <xi-c{i-1],...,xm-c{m-1}> (for i=2,...,m)
663  vector<poly> simple(x.size(), poly(BHEAD 0));
664  for (int i=(int)x.size()-2; i>=0; i--)
665  simple[i] = poly::simple_poly(BHEAD x[i+1],c[i],1);
666 
667  // Calculate the maximum degree of A in x2,...,xm
668  int maxdegA=0;
669  for (int i=1; i<(int)x.size(); i++)
670  maxdegA = MaX(maxdegA, A.degree(x[i]));
671 
672  // Iteratively add the variables x2,...,xm
673  for (int xi=1; xi<(int)x.size(); xi++) {
674  // replace the current leading coefficients by the correct ones
675  for (int i=0; i<(int)a.size(); i++)
676  a[i] += (lc[i] - a[i].lcoeff_univar(x[0])) * poly::simple_poly(BHEAD x[0],0,a[i].degree(x[0]));
677 
678  vector<poly> anew(a);
679  for (int i=0; i<(int)anew.size(); i++)
680  for (int j=xi-1; j<(int)c.size(); j++)
681  anew[i] %= simple[j];
682 
683  vector<int> xnew(x.begin(), x.begin()+xi);
684  vector<int> cnew(c.begin(), c.begin()+xi-1);
685  poly term(BHEAD 1,A.modp,A.modn);
686 
687  // Iteratively add the powers xi^k
688  for (int deg=1, maxdeg=A.degree(x[xi]); deg<=maxdeg; deg++) {
689 
690  term *= simple[xi-1];
691 
692  // Calculate the error, express it in terms of ai and add corrections.
693  poly error(BHEAD -1,A.modp,A.modn);
694  for (int i=0; i<(int)a.size(); i++) error *= a[i];
695  error += A;
696  for (int i=xi; i<(int)c.size(); i++) error %= simple[i];
697 
698  if (error.is_zero()) break;
699 
700  error /= term;
701  error %= simple[xi-1];
702 
703  vector<poly> s(solve_Diophantine_multivariate(anew,error,xnew,cnew,maxdegA));
704  if (s == vector<poly>()) return vector<poly>();
705 
706  for (int i=0; i<(int)a.size(); i++)
707  a[i] += s[i] * term;
708  }
709 
710  // check whether PRODUCT(a[i]) = A mod <xi-c{i-1},...,xm-c{m-1]> over the integers or ZZ/p
711  poly check(BHEAD -1, A.modn>1?0:A.modp, 1);
712  for (int i=0; i<(int)a.size(); i++) check *= a[i];
713  check += A;
714  for (int i=xi; i<(int)c.size(); i++) check %= simple[i];
715  if (!check.is_zero()) return vector<poly>();
716  }
717 
718 #ifdef DEBUG
719  cout << "*** [" << thetime() << "] RES : lift_variables("<<A<<","<<_a<<","<<x<<","<<c<<","<<lc<<") = " << a << endl;
720 #endif
721 
722  return a;
723 }
724 
725 /*
726  #] lift_variables :
727  #[ choose_prime :
728 */
729 
741 WORD polyfact::choose_prime (const poly &a, const vector<int> &x, WORD p) {
742 
743 #ifdef DEBUG
744  cout << "*** [" << thetime() << "] CALL: choose_prime("<<a<<","<<x<<","<<p<<")"<<endl;
745 #endif
746 
747  POLY_GETIDENTITY(a);
748 
749  poly icont_lcoeff(polygcd::integer_content(a.lcoeff_univar(x[0])));
750 
751  if (p==0) p = POLYFACT_FIRST_PRIME;
752 
753  poly icont_lcoeff_modp(BHEAD 0);
754 
755  do {
756 
757  bool is_prime;
758 
759  do {
760  p += 2;
761  is_prime = true;
762  for (int d=2; d*d<=p; d++)
763  if (p%d==0) { is_prime=false; break; }
764  }
765  while (!is_prime);
766 
767  icont_lcoeff_modp = icont_lcoeff;
768  icont_lcoeff_modp.setmod(p,1);
769  }
770  while (icont_lcoeff_modp.is_zero());
771 
772 #ifdef DEBUG
773  cout << "*** [" << thetime() << "] RES : choose_prime("<<a<<","<<x<<",?) = "<<p<<endl;
774 #endif
775 
776  return p;
777 }
778 
779 /*
780  #] choose_prime :
781  #[ choose_prime_power :
782 */
783 
798 WORD polyfact::choose_prime_power (const poly &a, WORD p) {
799 
800 #ifdef DEBUG
801  cout << "*** [" << thetime() << "] CALL: choose_prime_power("<<a<<","<<p<<")"<<endl;
802 #endif
803 
804  POLY_GETIDENTITY(a);
805 
806  // analyse the polynomial for calculating the bound
807  double maxdegree=0, maxlogcoeff=0, numterms=0;
808 
809  for (int i=1; i<a[0]; i+=a[i]) {
810  for (int j=0; j<AN.poly_num_vars; j++)
811  maxdegree = MaX(maxdegree, a[i+1+j]);
812 
813  maxlogcoeff = MaX(maxlogcoeff,
814  log(1.0+(UWORD)a[i+a[i]-2]) + // most significant digit + 1
815  BITSINWORD*log(2.0)*(ABS(a[i+a[i]-1])-1)); // number of digits
816  numterms++;
817  }
818 
819  WORD res = (WORD)ceil((log((sqrt(5.0)+1)/2)*maxdegree + maxlogcoeff + 0.5*log(numterms)) / log((double)p));
820 
821 #ifdef DEBUG
822  cout << "*** [" << thetime() << "] CALL: choose_prime_power("<<a<<","<<p<<") = "<<res<<endl;
823 #endif
824 
825  return res;
826 }
827 
828 /*
829  #] choose_prime_power :
830  #[ choose_ideal :
831 */
832 
860 const vector<int> polyfact::choose_ideal (const poly &a, int p, const factorized_poly &lc, const vector<int> &x) {
861 
862 #ifdef DEBUG
863  cout << "*** [" << thetime() << "] CALL: polyfact::choose_ideal("
864  <<a<<","<<p<<","<<lc<<","<<x<<")"<<endl;
865 #endif
866 
867  if (x.size()==1) return vector<int>();
868 
869  POLY_GETIDENTITY(a);
870 
871  vector<int> c(x.size()-1);
872 
873  int dega = a.degree(x[0]);
874  poly amodI(a);
875 
876  // choose random c
877  for (int i=0; i<(int)c.size(); i++) {
878  c[i] = 1 + wranf(BHEAD0) % ((p-1) / POLYFACT_IDEAL_FRACTION);
879  amodI %= poly::simple_poly(BHEAD x[i+1],c[i],1);
880  }
881 
882  poly amodIp(amodI);
883  amodIp.setmod(p,1);
884 
885  // check if leading coefficient is non-zero [equivalent to degree=old_degree]
886  if (amodIp.degree(x[0]) != dega)
887  return c = vector<int>();
888 
889  // check if leading coefficient is squarefree [equivalent to gcd(a,a')==const]
890  if (!polygcd::gcd_Euclidean(amodIp, amodIp.derivative(x[0])).is_integer())
891  return c = vector<int>();
892 
893  if (a.modp>0 && a.modn==1) return c;
894 
895  // check for unique prime factors in each factor lc[i] of the leading coefficient
896  vector<poly> d(1, polygcd::integer_content(amodI));
897 
898  for (int i=0; i<(int)lc.factor.size(); i++) {
899  // constant factor
900  if (i==0 && lc.factor[i].is_integer()) {
901  d[0] *= lc.factor[i];
902  continue;
903  }
904 
905  // factor modulo I
906  poly q(lc.factor[i]);
907  for (int j=0; j<(int)c.size(); j++)
908  q %= poly::simple_poly(BHEAD x[j+1],c[j]);
909  if (q.sign() == -1) q *= poly(BHEAD -1);
910 
911  // divide out common factors
912  for (int j=(int)d.size()-1; j>=0; j--) {
913  poly r(d[j]);
914  while (!r.is_one()) {
915  r = polygcd::integer_gcd(r,q);
916  q /= r;
917  }
918  }
919 
920  // check whether there is some factor left
921  if (q.is_one()) return vector<int>();
922  d.push_back(q);
923  }
924 
925 #ifdef DEBUG
926  cout << "*** [" << thetime() << "] RES : polyfact::choose_ideal("
927  <<a<<","<<p<<","<<lc<<","<<x<<") = "<<c<<endl;
928 #endif
929 
930  return c;
931 }
932 
933 /*
934  #] choose_ideal :
935  #[ squarefree_factors_Yun :
936 */
937 
944 const factorized_poly polyfact::squarefree_factors_Yun (const poly &_a) {
945 
946  factorized_poly res;
947  poly a(_a);
948 
949  int pow = 1;
950  int x = a.first_variable();
951 
952  poly b(a.derivative(x));
953  poly c(polygcd::gcd(a,b));
954 
955  while (true) {
956  a /= c;
957  b /= c;
958  b -= a.derivative(x);
959  if (b.is_zero()) break;
960  c = polygcd::gcd(a,b);
961  if (!c.is_one()) res.add_factor(c,pow);
962  pow++;
963  }
964 
965  if (!a.is_one()) res.add_factor(a,pow);
966  return res;
967 }
968 
969 /*
970  #] squarefree_factors_Yun :
971  #[ squarefree_factors_modp :
972 */
973 
980 const factorized_poly polyfact::squarefree_factors_modp (const poly &_a) {
981 
982  factorized_poly res;
983  poly a(_a);
984 
985  int pow = 1;
986  int x = a.first_variable();
987  poly b(a.derivative(x));
988 
989  // poly contains terms of the form c(x)^n (n!=c*p)
990  if (!b.is_zero()) {
991  poly c(polygcd::gcd(a,b));
992  a /= c;
993 
994  while (!a.is_one()) {
995  b = polygcd::gcd(a,c);
996  a /= b;
997  if (!a.is_one()) res.add_factor(a,pow);
998  pow++;
999  a = b;
1000  c /= a;
1001  }
1002 
1003  a = c;
1004  }
1005 
1006  // polynomial contains terms of the form c(x)^p
1007  if (!a.is_one()) {
1008  for (int i=1; i<a[1]; i+=a[i])
1009  a[i+1+x] /= a.modp;
1010  factorized_poly res2(squarefree_factors(a));
1011  for (int i=0; i<(int)res2.factor.size(); i++) {
1012  res.factor.push_back(res2.factor[i]);
1013  res.power.push_back(a.modp*res2.power[i]);
1014  }
1015  }
1016 
1017  return res;
1018 }
1019 
1020 /*
1021  #] squarefree_factors_modp :
1022  #[ squarefree_factors :
1023 */
1024 
1045 const factorized_poly polyfact::squarefree_factors (const poly &a) {
1046 
1047 #ifdef DEBUG
1048  cout << "*** [" << thetime() << "] CALL: squarefree_factors("<<a<<")\n";
1049 #endif
1050 
1051  if (a.is_one()) return factorized_poly();
1052 
1053  factorized_poly res;
1054 
1055  if (a.modp==0)
1056  res = squarefree_factors_Yun(a);
1057  else
1058  res = squarefree_factors_modp(a);
1059 
1060 #ifdef DEBUG
1061  cout << "*** [" << thetime() << "] RES : squarefree_factors("<<a<<") = "<<res<<"\n";
1062 #endif
1063 
1064  return res;
1065 }
1066 
1067 /*
1068  #] squarefree_factors :
1069  #[ Berlekamp_Qmatrix :
1070 */
1071 
1094 const vector<vector<WORD> > polyfact::Berlekamp_Qmatrix (const poly &_a) {
1095 
1096 #ifdef DEBUG
1097  cout << "*** [" << thetime() << "] CALL: Berlekamp_Qmatrix("<<_a<<")\n";
1098 #endif
1099 
1100  if (_a.all_variables() == vector<int>())
1101  return vector<vector<WORD> >(0);
1102 
1103  POLY_GETIDENTITY(_a);
1104 
1105  poly a(_a);
1106  int x = a.first_variable();
1107  int n = a.degree(x);
1108  int p = a.modp;
1109 
1110  poly lc(a.integer_lcoeff());
1111  a /= lc;
1112 
1113  vector<vector<WORD> > Q(n, vector<WORD>(n));
1114 
1115  // c is the vector of coefficients of the polynomial a
1116  vector<WORD> c(n+1,0);
1117  for (int j=1; j<a[0]; j+=a[j])
1118  c[a[j+1+x]] = a[j+1+AN.poly_num_vars] * a[j+2+AN.poly_num_vars];
1119 
1120  // d is the vector of coefficients of x^i mod a, starting with i=0
1121  vector<WORD> d(n,0);
1122  d[0]=1;
1123 
1124  for (int i=0; i<=(n-1)*p; i++) {
1125  // store the coefficients of x^(i*p) mod a
1126  if (i%p==0) Q[i/p] = d;
1127 
1128  // transform d=x^i mod a into d=x^(i+1) mod a
1129  vector<WORD> e(n);
1130  for (int j=0; j<n; j++) {
1131  e[j] = (-(LONG)d[n-1]*c[j] + (j>0?d[j-1]:0)) % p;
1132  if (e[j]<0) e[j]+=p;
1133  }
1134 
1135  d=e;
1136  }
1137 
1138  // Q = Q - I
1139  for (int i=0; i<n; i++)
1140  Q[i][i] = (Q[i][i] - 1 + p) % p;
1141 
1142  // Gaussian elimination
1143  for (int i=0; i<n; i++) {
1144  // Find pivot
1145  int ii=i; while (ii<n && Q[i][ii]==0) ii++;
1146  if (ii==n) continue;
1147 
1148  for (int k=0; k<n; k++)
1149  swap(Q[k][ii],Q[k][i]);
1150 
1151  // normalize row i, which becomes the pivot
1152  WORD mul;
1153  GetModInverses (Q[i][i], p, &mul, NULL);
1154  vector<int> idx;
1155 
1156  for (int k=0; k<n; k++) if (Q[k][i] != 0) {
1157  // store indices of non-zero elements for sparse matrices
1158  idx.push_back(k);
1159  Q[k][i] = ((LONG)Q[k][i] * mul) % p;
1160  }
1161 
1162  // reduce
1163  for (int j=0; j<n; j++)
1164  if (j!=i && Q[i][j]!=0) {
1165  LONG mul = Q[i][j];
1166  for (int k=0; k<(int)idx.size(); k++) {
1167  Q[idx[k]][j] = (Q[idx[k]][j] - mul*Q[idx[k]][i]) % p;
1168  if (Q[idx[k]][j] < 0) Q[idx[k]][j]+=p;
1169  }
1170  }
1171  }
1172 
1173  for (int i=0; i<n; i++) {
1174  // Q = Q - I
1175  Q[i][i] = Q[i][i]-1;
1176 
1177  // reduce all coefficients in the range 0,1,...,p-1
1178  for (int j=0; j<n; j++) {
1179  Q[i][j] = -Q[i][j]%p;
1180  if (Q[i][j]<0) Q[i][j]+=p;
1181  }
1182  }
1183 
1184 #ifdef DEBUG
1185  cout << "*** [" << thetime() << "] RES : Berlekamp_Qmatrix("<<_a<<") = "<<Q<<"\n";
1186 #endif
1187 
1188  return Q;
1189 }
1190 
1191 /*
1192  #] Berlekamp_Qmatrix :
1193  #[ Berlekamp_find_factors :
1194 */
1195 
1219 const vector<poly> polyfact::Berlekamp_find_factors (const poly &_a, const vector<vector<WORD> > &_q) {
1220 
1221 #ifdef DEBUG
1222  cout << "*** [" << thetime() << "] CALL: Berlekamp_find_factors("<<_a<<","<<_q<<")\n";
1223 #endif
1224 
1225  if (_a.all_variables() == vector<int>()) return vector<poly>(1,_a);
1226 
1227  POLY_GETIDENTITY(_a);
1228 
1229  vector<vector<WORD> > q=_q;
1230 
1231  int rank=0;
1232  for (int i=0; i<(int)q.size(); i++)
1233  if (q[i]!=vector<WORD>(q[i].size(),0)) rank++;
1234 
1235  poly a(_a);
1236  int x = a.first_variable();
1237  int n = a.degree(x);
1238  int p = a.modp;
1239 
1240  a /= a.integer_lcoeff();
1241 
1242  // Vector of factors, represented as dense polynomials mod p
1243  vector<vector<WORD> > fac(1, vector<WORD>(n+1,0));
1244 
1245  fac[0] = poly::to_coefficient_list(a);
1246  bool finished=false;
1247 
1248  // Loop over the columns of q + constant, i.e., an exhaustive list of possible factors
1249  for (int i=1; i<n && !finished; i++) {
1250  if (q[i] == vector<WORD>(n,0)) continue;
1251 
1252  for (int s=0; s<p && !finished; s++) {
1253  for (int j=0; j<(int)fac.size() && !finished; j++) {
1254  vector<WORD> c = polygcd::coefficient_list_gcd(fac[j], q[i], p);
1255 
1256  // If a non-trivial factor is found, add it to the list
1257  if (c.size()!=1 && c.size()!=fac[j].size()) {
1258  fac.push_back(c);
1259  fac[j] = poly::coefficient_list_divmod(fac[j], c, p, 0);
1260  if ((int)fac.size() == rank) finished=true;
1261  }
1262  }
1263 
1264  // Increase the constant term by one
1265  q[i][0] = (q[i][0]+1) % p;
1266  }
1267  }
1268 
1269  // Convert the densely represented polynomials to sparse ones
1270  vector<poly> res(fac.size(),poly(BHEAD 0, p));
1271  for (int i=0; i<(int)fac.size(); i++)
1272  res[i] = poly::from_coefficient_list(BHEAD fac[i],x,p);
1273 
1274 #ifdef DEBUG
1275  cout << "*** [" << thetime() << "] RES : Berlekamp_find_factors("<<_a<<","<<_q<<") = "<<res<<"\n";
1276 #endif
1277 
1278  return res;
1279 }
1280 
1281 /*
1282  #] Berlekamp_find_factors :
1283  #[ combine_factors :
1284 */
1285 
1307 const vector<poly> polyfact::combine_factors (const poly &a, const vector<poly> &f) {
1308 
1309 #ifdef DEBUG
1310  cout << "*** [" << thetime() << "] CALL: combine_factors("<<a<<","<<f<<")\n";
1311 #endif
1312 
1313  POLY_GETIDENTITY(a);
1314 
1315  poly a0(a,0,1);
1316  vector<poly> res;
1317 
1318  int num_used = 0;
1319  vector<bool> used(f.size(), false);
1320 
1321  // Loop over all bitmasks with num=1,2,...,size(factors)/2 bits
1322  // set, that contain only unused factors
1323  for (int num=1; num<=(int)(f.size() - num_used)/2; num++) {
1324  vector<int> next(f.size() - num_used, 0);
1325  for (int i=0; i<num; i++) next[next.size()-1-i] = 1;
1326 
1327  do {
1328  poly fac(BHEAD 1,a.modp,a.modn);
1329  for (int i=0, j=0; i<(int)f.size(); i++)
1330  if (!used[i] && next[j++]) fac *= f[i];
1331  fac /= fac.integer_lcoeff();
1332  fac *= a.integer_lcoeff();
1333  fac /= polygcd::integer_content(poly(fac,0,1));
1334 
1335  if ((a0 % fac).is_zero()) {
1336  res.push_back(fac);
1337  for (int i=0, j=0; i<(int)f.size(); i++)
1338  if (!used[i]) used[i] = next[j++];
1339  num_used += num;
1340  num--;
1341  break;
1342  }
1343  }
1344  while (next_permutation(next.begin(), next.end()));
1345  }
1346 
1347  // All unused factors together form one more factor
1348  if (num_used != (int)f.size()) {
1349  poly fac(BHEAD 1,a.modp,a.modn);
1350  for (int i=0; i<(int)f.size(); i++)
1351  if (!used[i]) fac *= f[i];
1352  fac /= fac.integer_lcoeff();
1353  fac *= a.integer_lcoeff();
1354  fac /= polygcd::integer_content(poly(fac,0,1));
1355  res.push_back(fac);
1356  }
1357 
1358 #ifdef DEBUG
1359  cout << "*** [" << thetime() << "] RES : combine_factors("<<a<<","<<f<<") = "<<res<<"\n";
1360 #endif
1361 
1362  return res;
1363 }
1364 
1365 /*
1366  #] combine_factors :
1367  #[ factorize_squarefree :
1368 */
1369 
1388 const vector<poly> polyfact::factorize_squarefree (const poly &a, const vector<int> &x) {
1389 
1390 #ifdef DEBUG
1391  cout << "*** [" << thetime() << "] CALL: factorize_squarefree("<<a<<")\n";
1392 #endif
1393 
1394  POLY_GETIDENTITY(a);
1395 
1396  WORD p=a.modp, n=a.modn;
1397 
1398  try_again:
1399 
1400  int bestp=0;
1401  int min_factors = INT_MAX;
1402  poly amodI(BHEAD 0), bestamodI(BHEAD 0);
1403  vector<int> c,d,bestc,bestd;
1404  vector<vector<WORD> > q,bestq;
1405 
1406  // Factorize leading coefficient
1407  factorized_poly lc(factorize(a.lcoeff_univar(x[0])));
1408 
1409  // Try a number of primes
1410  int prime_tries = 0;
1411 
1412  while (prime_tries<POLYFACT_NUM_CONFIRMATIONS && min_factors>1) {
1413  if (a.modp == 0) {
1414  p = choose_prime(a,x,p);
1415  n = 0;
1416  if (a.degree(x[0]) % p == 0) continue;
1417 
1418  // Univariate case: check whether the polynomial mod p is squarefree
1419  // Multivariate case: this check is done after choosing I (for efficiency)
1420  if (x.size()==1) {
1421  poly amodp(a,p,1);
1422  if (polygcd::gcd_Euclidean(amodp, amodp.derivative(x[0])).degree(x[0]) != 0)
1423  continue;
1424  }
1425  }
1426 
1427  // Try a number of ideals
1428  if (x.size()>1)
1429  for (int ideal_tries=0; ideal_tries<POLYFACT_MAX_IDEAL_TRIES; ideal_tries++) {
1430  c = choose_ideal(a,p,lc,x);
1431  if (c.size()>0) break;
1432  }
1433 
1434  if (x.size()==1 || c.size()>0) {
1435  amodI = a;
1436  for (int i=0; i<(int)c.size(); i++)
1437  amodI %= poly::simple_poly(BHEAD x[i+1],c[i]);
1438 
1439  // Determine Q-matrix and its rank. Smaller rank is better.
1440  q = Berlekamp_Qmatrix(poly(amodI,p,1));
1441  int rank=0;
1442  for (int i=0; i<(int)q.size(); i++)
1443  if (q[i]!=vector<WORD>(q[i].size(),0)) rank++;
1444 
1445  if (rank<min_factors) {
1446  bestp=p;
1447  bestc=c;
1448  bestq=q;
1449  bestamodI=amodI;
1450  min_factors = rank;
1451  prime_tries = 0;
1452  }
1453 
1454  if (rank==min_factors)
1455  prime_tries++;
1456 
1457 #ifdef DEBUG
1458  cout << "*** [" << thetime() << "] (A) : factorize_squarefree("<<a<<
1459  ") try p=" << p << " #factors=" << rank << " (min="<<min_factors<<"x"<<prime_tries<<")" << endl;
1460 #endif
1461  }
1462  }
1463 
1464  p=bestp;
1465  c=bestc;
1466  q=bestq;
1467  amodI=bestamodI;
1468 
1469  // Determine to which power of p to lift
1470  if (n==0) {
1471  n = choose_prime_power(amodI,p);
1472  n = MaX(n, choose_prime_power(a,p));
1473  }
1474 
1475  amodI.setmod(p,n);
1476 
1477 #ifdef DEBUG
1478  cout << "*** [" << thetime() << "] (B) : factorize_squarefree("<<a<<
1479  ") chosen c = " << c << ", p^n = "<<p<<"^"<<n<<endl;
1480  cout << "*** [" << thetime() << "] (C) : factorize_squarefree("<<a<<
1481  ") #factors = " << min_factors << endl;
1482 #endif
1483 
1484  // Find factors
1485  vector<poly> f(Berlekamp_find_factors(poly(amodI,p,1),q));
1486 
1487  // Lift coefficients
1488  if (f.size() > 1) {
1489  f = lift_coefficients(amodI,f);
1490  if (f==vector<poly>()) {
1491 #ifdef DEBUG
1492  cout << "factor_squarefree failed (lift_coeff step) : " << endl;
1493 #endif
1494  goto try_again;
1495  }
1496 
1497  // Combine factors
1498  if (a.modp==0) f = combine_factors(amodI,f);
1499  }
1500 
1501  // Lift variables
1502  if (f.size() == 1)
1503  f = vector<poly>(1, a);
1504 
1505  if (x.size() > 1 && f.size() > 1) {
1506 
1507  // The correct leading coefficients of the factors can be
1508  // reconstructed from prime number factors of the leading
1509  // coefficients modulo I. This is possible since all factors of
1510  // the leading coefficient have unique prime factors for the ideal
1511  // I is chosen as such.
1512 
1513  poly amodI(a);
1514  for (int i=0; i<(int)c.size(); i++)
1515  amodI %= poly::simple_poly(BHEAD x[i+1],c[i]);
1516  poly delta(polygcd::integer_content(amodI));
1517 
1518  vector<poly> lcmodI(lc.factor.size(), poly(BHEAD 0));
1519  for (int i=0; i<(int)lc.factor.size(); i++) {
1520  lcmodI[i] = lc.factor[i];
1521  for (int j=0; j<(int)c.size(); j++)
1522  lcmodI[i] %= poly::simple_poly(BHEAD x[j+1],c[j]);
1523  }
1524 
1525  vector<poly> correct_lc(f.size(), poly(BHEAD 1,p,n));
1526 
1527  for (int j=0; j<(int)f.size(); j++) {
1528  poly lc_f(f[j].integer_lcoeff() * delta);
1529  WORD nlc_f = lc_f[lc_f[1]];
1530  poly quo(BHEAD 0),rem(BHEAD 0);
1531  WORD nquo,nrem;
1532 
1533  for (int i=(int)lcmodI.size()-1; i>=0; i--) {
1534 
1535  if (i==0 && lc.factor[i].is_integer()) continue;
1536 
1537  do {
1538  DivLong((UWORD *)&lc_f[2+AN.poly_num_vars], nlc_f,
1539  (UWORD *)&lcmodI[i][2+AN.poly_num_vars], lcmodI[i][lcmodI[i][1]],
1540  (UWORD *)&quo[0], &nquo,
1541  (UWORD *)&rem[0], &nrem);
1542 
1543  if (nrem == 0) {
1544  correct_lc[j] *= lc.factor[i];
1545  lc_f.termscopy(&quo[0], 2+AN.poly_num_vars, ABS(nquo));
1546  nlc_f = nquo;
1547  }
1548  }
1549  while (nrem == 0);
1550  }
1551  }
1552 
1553  for (int i=0; i<(int)correct_lc.size(); i++) {
1554  poly correct_modI(correct_lc[i]);
1555  for (int j=0; j<(int)c.size(); j++)
1556  correct_modI %= poly::simple_poly(BHEAD x[j+1],c[j]);
1557 
1558  poly d(polygcd::integer_gcd(correct_modI, f[i].integer_lcoeff()));
1559  correct_lc[i] *= f[i].integer_lcoeff() / d;
1560  delta /= correct_modI / d;
1561  f[i] *= correct_modI / d;
1562  }
1563 
1564  // increase n, because of multiplying with delta
1565  if (!delta.is_one()) {
1566  poly deltapow(BHEAD 1);
1567  for (int i=1; i<(int)correct_lc.size(); i++)
1568  deltapow *= delta;
1569  while (!deltapow.is_zero()) {
1570  deltapow /= poly(BHEAD p);
1571  n++;
1572  }
1573 
1574  for (int i=0; i<(int)f.size(); i++) {
1575  f[i].modn = n;
1576  correct_lc[i].modn = n;
1577  }
1578  }
1579 
1580  poly aa(a,p,n);
1581 
1582  for (int i=0; i<(int)correct_lc.size(); i++) {
1583  correct_lc[i] *= delta;
1584  f[i] *= delta;
1585  if (i>0) aa *= delta;
1586  }
1587 
1588  f = lift_variables(aa,f,x,c,correct_lc);
1589 
1590  for (int i=0; i<(int)f.size(); i++)
1591  if (a.modp == 0)
1592  f[i] /= polygcd::integer_content(poly(f[i],0,1));
1593  else
1594  f[i] /= polygcd::content_univar(f[i], x[0]);
1595 
1596  if (f==vector<poly>()) {
1597 #ifdef DEBUG
1598  cout << "factor_squarefree failed (lift_var step)" << endl;
1599 #endif
1600  goto try_again;
1601  }
1602  }
1603 
1604  // set modulus of the factors correctly
1605  if (a.modp==0)
1606  for (int i=0; i<(int)f.size(); i++)
1607  f[i].setmod(0,1);
1608 
1609  // Final check (not sure if this is necessary, but it doesn't hurt)
1610  poly check(BHEAD 1,a.modp,a.modn);
1611  for (int i=0; i<(int)f.size(); i++)
1612  check *= f[i];
1613 
1614  if (check != a) {
1615 #ifdef DEBUG
1616  cout << "factor_squarefree failed (final check) : " << f << endl;
1617 #endif
1618  goto try_again;
1619  }
1620 
1621 #ifdef DEBUG
1622  cout << "*** [" << thetime() << "] RES : factorize_squarefree("<<a<<","<<x<<","<<c<<") = "<<f<<"\n";
1623 #endif
1624 
1625  return f;
1626 }
1627 
1628 /*
1629  #] factorize_squarefree :
1630  #[ factorize :
1631 */
1632 
1641 const factorized_poly polyfact::factorize (const poly &a) {
1642 
1643 #ifdef DEBUG
1644  cout << "*** [" << thetime() << "] CALL: factorize("<<a<<")\n";
1645 #endif
1646  vector<int> x = a.all_variables();
1647 
1648  // No variables, so just one factor
1649  if (x.size() == 0) {
1650  factorized_poly res;
1651  if (a.is_one()) return res;
1652  res.add_factor(a,1);
1653  return res;
1654  }
1655 
1656  // Remove content
1657  poly conta(polygcd::content_univar(a,x[0]));
1658 
1659  factorized_poly faca(factorize(conta));
1660 
1661  poly ppa(a / conta);
1662 
1663  // Find a squarefree factorization
1664  factorized_poly b(squarefree_factors(ppa));
1665 
1666 #ifdef DEBUG
1667  cout << "*** [" << thetime() << "] ... : factorize("<<a<<") : SFF = "<<b<<"\n";
1668 #endif
1669 
1670  // Factorize each squarefree factor and build the "factorized_poly"
1671  for (int i=0; i<(int)b.factor.size(); i++) {
1672 
1673  vector<poly> c(factorize_squarefree(b.factor[i], x));
1674 
1675  for (int j=0; j<(int)c.size(); j++) {
1676  faca.factor.push_back(c[j]);
1677  faca.power.push_back(b.power[i]);
1678  }
1679  }
1680 
1681 #ifdef DEBUG
1682  cout << "*** [" << thetime() << "] RES : factorize("<<a<<") = "<<faca<<"\n";
1683 #endif
1684 
1685  return faca;
1686 }
1687 
1688 /*
1689  #] factorize :
1690 */
int GetModInverses(WORD, WORD, WORD *, WORD *)
Definition: reken.c:1466
Definition: poly.h:49