GiNaC  1.8.3
factor.cpp
Go to the documentation of this file.
1 
35 /*
36  * GiNaC Copyright (C) 1999-2022 Johannes Gutenberg University Mainz, Germany
37  *
38  * This program is free software; you can redistribute it and/or modify
39  * it under the terms of the GNU General Public License as published by
40  * the Free Software Foundation; either version 2 of the License, or
41  * (at your option) any later version.
42  *
43  * This program is distributed in the hope that it will be useful,
44  * but WITHOUT ANY WARRANTY; without even the implied warranty of
45  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
46  * GNU General Public License for more details.
47  *
48  * You should have received a copy of the GNU General Public License
49  * along with this program; if not, write to the Free Software
50  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
51  */
52 
53 //#define DEBUGFACTOR
54 
55 #include "factor.h"
56 
57 #include "ex.h"
58 #include "numeric.h"
59 #include "operators.h"
60 #include "inifcns.h"
61 #include "symbol.h"
62 #include "relational.h"
63 #include "power.h"
64 #include "mul.h"
65 #include "normal.h"
66 #include "add.h"
67 
68 #include <algorithm>
69 #include <limits>
70 #include <list>
71 #include <vector>
72 #include <stack>
73 #ifdef DEBUGFACTOR
74 #include <ostream>
75 #endif
76 using namespace std;
77 
78 #include <cln/cln.h>
79 using namespace cln;
80 
81 namespace GiNaC {
82 
83 #ifdef DEBUGFACTOR
84 #define DCOUT(str) cout << #str << endl
85 #define DCOUTVAR(var) cout << #var << ": " << var << endl
86 #define DCOUT2(str,var) cout << #str << ": " << var << endl
87 ostream& operator<<(ostream& o, const vector<int>& v)
88 {
89  auto i = v.begin(), end = v.end();
90  while ( i != end ) {
91  o << *i << " ";
92  ++i;
93  }
94  return o;
95 }
96 static ostream& operator<<(ostream& o, const vector<cl_I>& v)
97 {
98  auto i = v.begin(), end = v.end();
99  while ( i != end ) {
100  o << *i << "[" << i-v.begin() << "]" << " ";
101  ++i;
102  }
103  return o;
104 }
105 static ostream& operator<<(ostream& o, const vector<cl_MI>& v)
106 {
107  auto i = v.begin(), end = v.end();
108  while ( i != end ) {
109  o << *i << "[" << i-v.begin() << "]" << " ";
110  ++i;
111  }
112  return o;
113 }
114 ostream& operator<<(ostream& o, const vector<numeric>& v)
115 {
116  for ( size_t i=0; i<v.size(); ++i ) {
117  o << v[i] << " ";
118  }
119  return o;
120 }
121 ostream& operator<<(ostream& o, const vector<vector<cl_MI>>& v)
122 {
123  auto i = v.begin(), end = v.end();
124  while ( i != end ) {
125  o << i-v.begin() << ": " << *i << endl;
126  ++i;
127  }
128  return o;
129 }
130 #else
131 #define DCOUT(str)
132 #define DCOUTVAR(var)
133 #define DCOUT2(str,var)
134 #endif // def DEBUGFACTOR
135 
136 // anonymous namespace to hide all utility functions
137 namespace {
138 
140 // modular univariate polynomial code
141 
142 typedef std::vector<cln::cl_MI> umodpoly;
143 typedef std::vector<cln::cl_I> upoly;
144 typedef vector<umodpoly> upvec;
145 
146 // COPY FROM UPOLY.HPP
147 
148 // CHANGED size_t -> int !!!
149 template<typename T> static int degree(const T& p)
150 {
151  return p.size() - 1;
152 }
153 
154 template<typename T> static typename T::value_type lcoeff(const T& p)
155 {
156  return p[p.size() - 1];
157 }
158 
159 static bool normalize_in_field(umodpoly& a)
160 {
161  if (a.size() == 0)
162  return true;
163  if ( lcoeff(a) == a[0].ring()->one() ) {
164  return true;
165  }
166 
167  const cln::cl_MI lc_1 = recip(lcoeff(a));
168  for (std::size_t k = a.size(); k-- != 0; )
169  a[k] = a[k]*lc_1;
170  return false;
171 }
172 
173 template<typename T> static void
174 canonicalize(T& p, const typename T::size_type hint = std::numeric_limits<typename T::size_type>::max())
175 {
176  if (p.empty())
177  return;
178 
179  std::size_t i = p.size() - 1;
180  // Be fast if the polynomial is already canonicalized
181  if (!zerop(p[i]))
182  return;
183 
184  if (hint < p.size())
185  i = hint;
186 
187  bool is_zero = false;
188  do {
189  if (!zerop(p[i])) {
190  ++i;
191  break;
192  }
193  if (i == 0) {
194  is_zero = true;
195  break;
196  }
197  --i;
198  } while (true);
199 
200  if (is_zero) {
201  p.clear();
202  return;
203  }
204 
205  p.erase(p.begin() + i, p.end());
206 }
207 
208 // END COPY FROM UPOLY.HPP
209 
210 static void expt_pos(umodpoly& a, unsigned int q)
211 {
212  if ( a.empty() ) return;
213  cl_MI zero = a[0].ring()->zero();
214  int deg = degree(a);
215  a.resize(degree(a)*q+1, zero);
216  for ( int i=deg; i>0; --i ) {
217  a[i*q] = a[i];
218  a[i] = zero;
219  }
220 }
221 
222 template<bool COND, typename T = void> struct enable_if
223 {
224  typedef T type;
225 };
226 
227 template<typename T> struct enable_if<false, T> { /* empty */ };
228 
229 template<typename T> struct uvar_poly_p
230 {
231  static const bool value = false;
232 };
233 
234 template<> struct uvar_poly_p<upoly>
235 {
236  static const bool value = true;
237 };
238 
239 template<> struct uvar_poly_p<umodpoly>
240 {
241  static const bool value = true;
242 };
243 
244 template<typename T>
245 // Don't define this for anything but univariate polynomials.
246 static typename enable_if<uvar_poly_p<T>::value, T>::type
247 operator+(const T& a, const T& b)
248 {
249  int sa = a.size();
250  int sb = b.size();
251  if ( sa >= sb ) {
252  T r(sa);
253  int i = 0;
254  for ( ; i<sb; ++i ) {
255  r[i] = a[i] + b[i];
256  }
257  for ( ; i<sa; ++i ) {
258  r[i] = a[i];
259  }
260  canonicalize(r);
261  return r;
262  }
263  else {
264  T r(sb);
265  int i = 0;
266  for ( ; i<sa; ++i ) {
267  r[i] = a[i] + b[i];
268  }
269  for ( ; i<sb; ++i ) {
270  r[i] = b[i];
271  }
272  canonicalize(r);
273  return r;
274  }
275 }
276 
277 template<typename T>
278 // Don't define this for anything but univariate polynomials. Otherwise
279 // overload resolution might fail (this actually happens when compiling
280 // GiNaC with g++ 3.4).
281 static typename enable_if<uvar_poly_p<T>::value, T>::type
282 operator-(const T& a, const T& b)
283 {
284  int sa = a.size();
285  int sb = b.size();
286  if ( sa >= sb ) {
287  T r(sa);
288  int i = 0;
289  for ( ; i<sb; ++i ) {
290  r[i] = a[i] - b[i];
291  }
292  for ( ; i<sa; ++i ) {
293  r[i] = a[i];
294  }
295  canonicalize(r);
296  return r;
297  }
298  else {
299  T r(sb);
300  int i = 0;
301  for ( ; i<sa; ++i ) {
302  r[i] = a[i] - b[i];
303  }
304  for ( ; i<sb; ++i ) {
305  r[i] = -b[i];
306  }
307  canonicalize(r);
308  return r;
309  }
310 }
311 
312 static upoly operator*(const upoly& a, const upoly& b)
313 {
314  upoly c;
315  if ( a.empty() || b.empty() ) return c;
316 
317  int n = degree(a) + degree(b);
318  c.resize(n+1, 0);
319  for ( int i=0 ; i<=n; ++i ) {
320  for ( int j=0 ; j<=i; ++j ) {
321  if ( j > degree(a) || (i-j) > degree(b) ) continue;
322  c[i] = c[i] + a[j] * b[i-j];
323  }
324  }
325  canonicalize(c);
326  return c;
327 }
328 
329 static umodpoly operator*(const umodpoly& a, const umodpoly& b)
330 {
331  umodpoly c;
332  if ( a.empty() || b.empty() ) return c;
333 
334  int n = degree(a) + degree(b);
335  c.resize(n+1, a[0].ring()->zero());
336  for ( int i=0 ; i<=n; ++i ) {
337  for ( int j=0 ; j<=i; ++j ) {
338  if ( j > degree(a) || (i-j) > degree(b) ) continue;
339  c[i] = c[i] + a[j] * b[i-j];
340  }
341  }
342  canonicalize(c);
343  return c;
344 }
345 
346 static upoly operator*(const upoly& a, const cl_I& x)
347 {
348  if ( zerop(x) ) {
349  upoly r;
350  return r;
351  }
352  upoly r(a.size());
353  for ( size_t i=0; i<a.size(); ++i ) {
354  r[i] = a[i] * x;
355  }
356  return r;
357 }
358 
359 static upoly operator/(const upoly& a, const cl_I& x)
360 {
361  if ( zerop(x) ) {
362  upoly r;
363  return r;
364  }
365  upoly r(a.size());
366  for ( size_t i=0; i<a.size(); ++i ) {
367  r[i] = exquo(a[i],x);
368  }
369  return r;
370 }
371 
372 static umodpoly operator*(const umodpoly& a, const cl_MI& x)
373 {
374  umodpoly r(a.size());
375  for ( size_t i=0; i<a.size(); ++i ) {
376  r[i] = a[i] * x;
377  }
378  canonicalize(r);
379  return r;
380 }
381 
382 static void upoly_from_ex(upoly& up, const ex& e, const ex& x)
383 {
384  // assert: e is in Z[x]
385  int deg = e.degree(x);
386  up.resize(deg+1);
387  int ldeg = e.ldegree(x);
388  for ( ; deg>=ldeg; --deg ) {
389  up[deg] = the<cl_I>(ex_to<numeric>(e.coeff(x, deg)).to_cl_N());
390  }
391  for ( ; deg>=0; --deg ) {
392  up[deg] = 0;
393  }
394  canonicalize(up);
395 }
396 
397 static void umodpoly_from_upoly(umodpoly& ump, const upoly& e, const cl_modint_ring& R)
398 {
399  int deg = degree(e);
400  ump.resize(deg+1);
401  for ( ; deg>=0; --deg ) {
402  ump[deg] = R->canonhom(e[deg]);
403  }
404  canonicalize(ump);
405 }
406 
407 static void umodpoly_from_ex(umodpoly& ump, const ex& e, const ex& x, const cl_modint_ring& R)
408 {
409  // assert: e is in Z[x]
410  int deg = e.degree(x);
411  ump.resize(deg+1);
412  int ldeg = e.ldegree(x);
413  for ( ; deg>=ldeg; --deg ) {
414  cl_I coeff = the<cl_I>(ex_to<numeric>(e.coeff(x, deg)).to_cl_N());
415  ump[deg] = R->canonhom(coeff);
416  }
417  for ( ; deg>=0; --deg ) {
418  ump[deg] = R->zero();
419  }
420  canonicalize(ump);
421 }
422 
423 #ifdef DEBUGFACTOR
424 static void umodpoly_from_ex(umodpoly& ump, const ex& e, const ex& x, const cl_I& modulus)
425 {
426  umodpoly_from_ex(ump, e, x, find_modint_ring(modulus));
427 }
428 #endif
429 
430 static ex upoly_to_ex(const upoly& a, const ex& x)
431 {
432  if ( a.empty() ) return 0;
433  ex e;
434  for ( int i=degree(a); i>=0; --i ) {
435  e += numeric(a[i]) * pow(x, i);
436  }
437  return e;
438 }
439 
440 static ex umodpoly_to_ex(const umodpoly& a, const ex& x)
441 {
442  if ( a.empty() ) return 0;
443  cl_modint_ring R = a[0].ring();
444  cl_I mod = R->modulus;
445  cl_I halfmod = (mod-1) >> 1;
446  ex e;
447  for ( int i=degree(a); i>=0; --i ) {
448  cl_I n = R->retract(a[i]);
449  if ( n > halfmod ) {
450  e += numeric(n-mod) * pow(x, i);
451  } else {
452  e += numeric(n) * pow(x, i);
453  }
454  }
455  return e;
456 }
457 
458 static upoly umodpoly_to_upoly(const umodpoly& a)
459 {
460  upoly e(a.size());
461  if ( a.empty() ) return e;
462  cl_modint_ring R = a[0].ring();
463  cl_I mod = R->modulus;
464  cl_I halfmod = (mod-1) >> 1;
465  for ( int i=degree(a); i>=0; --i ) {
466  cl_I n = R->retract(a[i]);
467  if ( n > halfmod ) {
468  e[i] = n-mod;
469  } else {
470  e[i] = n;
471  }
472  }
473  return e;
474 }
475 
476 static umodpoly umodpoly_to_umodpoly(const umodpoly& a, const cl_modint_ring& R, unsigned int m)
477 {
478  umodpoly e;
479  if ( a.empty() ) return e;
480  cl_modint_ring oldR = a[0].ring();
481  size_t sa = a.size();
482  e.resize(sa+m, R->zero());
483  for ( size_t i=0; i<sa; ++i ) {
484  e[i+m] = R->canonhom(oldR->retract(a[i]));
485  }
486  canonicalize(e);
487  return e;
488 }
489 
497 static void reduce_coeff(umodpoly& a, const cl_I& x)
498 {
499  if ( a.empty() ) return;
500 
501  cl_modint_ring R = a[0].ring();
502  for (auto & i : a) {
503  // cln cannot perform this division in the modular field
504  cl_I c = R->retract(i);
505  i = cl_MI(R, the<cl_I>(c / x));
506  }
507 }
508 
516 static void rem(const umodpoly& a, const umodpoly& b, umodpoly& r)
517 {
518  int k, n;
519  n = degree(b);
520  k = degree(a) - n;
521  r = a;
522  if ( k < 0 ) return;
523 
524  do {
525  cl_MI qk = div(r[n+k], b[n]);
526  if ( !zerop(qk) ) {
527  for ( int i=0; i<n; ++i ) {
528  unsigned int j = n + k - 1 - i;
529  r[j] = r[j] - qk * b[j-k];
530  }
531  }
532  } while ( k-- );
533 
534  fill(r.begin()+n, r.end(), a[0].ring()->zero());
535  canonicalize(r);
536 }
537 
545 static void div(const umodpoly& a, const umodpoly& b, umodpoly& q)
546 {
547  int k, n;
548  n = degree(b);
549  k = degree(a) - n;
550  q.clear();
551  if ( k < 0 ) return;
552 
553  umodpoly r = a;
554  q.resize(k+1, a[0].ring()->zero());
555  do {
556  cl_MI qk = div(r[n+k], b[n]);
557  if ( !zerop(qk) ) {
558  q[k] = qk;
559  for ( int i=0; i<n; ++i ) {
560  unsigned int j = n + k - 1 - i;
561  r[j] = r[j] - qk * b[j-k];
562  }
563  }
564  } while ( k-- );
565 
566  canonicalize(q);
567 }
568 
577 static void remdiv(const umodpoly& a, const umodpoly& b, umodpoly& r, umodpoly& q)
578 {
579  int k, n;
580  n = degree(b);
581  k = degree(a) - n;
582  q.clear();
583  r = a;
584  if ( k < 0 ) return;
585 
586  q.resize(k+1, a[0].ring()->zero());
587  do {
588  cl_MI qk = div(r[n+k], b[n]);
589  if ( !zerop(qk) ) {
590  q[k] = qk;
591  for ( int i=0; i<n; ++i ) {
592  unsigned int j = n + k - 1 - i;
593  r[j] = r[j] - qk * b[j-k];
594  }
595  }
596  } while ( k-- );
597 
598  fill(r.begin()+n, r.end(), a[0].ring()->zero());
599  canonicalize(r);
600  canonicalize(q);
601 }
602 
609 static void gcd(const umodpoly& a, const umodpoly& b, umodpoly& c)
610 {
611  if ( degree(a) < degree(b) ) return gcd(b, a, c);
612 
613  c = a;
614  normalize_in_field(c);
615  umodpoly d = b;
616  normalize_in_field(d);
617  umodpoly r;
618  while ( !d.empty() ) {
619  rem(c, d, r);
620  c = d;
621  d = r;
622  }
623  normalize_in_field(c);
624 }
625 
631 static void deriv(const umodpoly& a, umodpoly& d)
632 {
633  d.clear();
634  if ( a.size() <= 1 ) return;
635 
636  d.insert(d.begin(), a.begin()+1, a.end());
637  int max = d.size();
638  for ( int i=1; i<max; ++i ) {
639  d[i] = d[i] * (i+1);
640  }
641  canonicalize(d);
642 }
643 
644 static bool unequal_one(const umodpoly& a)
645 {
646  if ( a.empty() ) return true;
647  return ( a.size() != 1 || a[0] != a[0].ring()->one() );
648 }
649 
650 static bool equal_one(const umodpoly& a)
651 {
652  return ( a.size() == 1 && a[0] == a[0].ring()->one() );
653 }
654 
660 static bool squarefree(const umodpoly& a)
661 {
662  umodpoly b;
663  deriv(a, b);
664  if ( b.empty() ) {
665  return false;
666  }
667  umodpoly c;
668  gcd(a, b, c);
669  return equal_one(c);
670 }
671 
672 // END modular univariate polynomial code
674 
676 // modular matrix
677 
678 typedef vector<cl_MI> mvec;
679 
680 class modular_matrix
681 {
682 #ifdef DEBUGFACTOR
683  friend ostream& operator<<(ostream& o, const modular_matrix& m);
684 #endif
685 public:
686  modular_matrix(size_t r_, size_t c_, const cl_MI& init) : r(r_), c(c_)
687  {
688  m.resize(c*r, init);
689  }
690  size_t rowsize() const { return r; }
691  size_t colsize() const { return c; }
692  cl_MI& operator()(size_t row, size_t col) { return m[row*c + col]; }
693  cl_MI operator()(size_t row, size_t col) const { return m[row*c + col]; }
694  void mul_col(size_t col, const cl_MI x)
695  {
696  for ( size_t rc=0; rc<r; ++rc ) {
697  std::size_t i = c*rc + col;
698  m[i] = m[i] * x;
699  }
700  }
701  void sub_col(size_t col1, size_t col2, const cl_MI fac)
702  {
703  for ( size_t rc=0; rc<r; ++rc ) {
704  std::size_t i1 = col1 + c*rc;
705  std::size_t i2 = col2 + c*rc;
706  m[i1] = m[i1] - m[i2]*fac;
707  }
708  }
709  void switch_col(size_t col1, size_t col2)
710  {
711  for ( size_t rc=0; rc<r; ++rc ) {
712  std::size_t i1 = col1 + rc*c;
713  std::size_t i2 = col2 + rc*c;
714  std::swap(m[i1], m[i2]);
715  }
716  }
717  void mul_row(size_t row, const cl_MI x)
718  {
719  for ( size_t cc=0; cc<c; ++cc ) {
720  std::size_t i = row*c + cc;
721  m[i] = m[i] * x;
722  }
723  }
724  void sub_row(size_t row1, size_t row2, const cl_MI fac)
725  {
726  for ( size_t cc=0; cc<c; ++cc ) {
727  std::size_t i1 = row1*c + cc;
728  std::size_t i2 = row2*c + cc;
729  m[i1] = m[i1] - m[i2]*fac;
730  }
731  }
732  void switch_row(size_t row1, size_t row2)
733  {
734  for ( size_t cc=0; cc<c; ++cc ) {
735  std::size_t i1 = row1*c + cc;
736  std::size_t i2 = row2*c + cc;
737  std::swap(m[i1], m[i2]);
738  }
739  }
740  bool is_col_zero(size_t col) const
741  {
742  for ( size_t rr=0; rr<r; ++rr ) {
743  std::size_t i = col + rr*c;
744  if ( !zerop(m[i]) ) {
745  return false;
746  }
747  }
748  return true;
749  }
750  bool is_row_zero(size_t row) const
751  {
752  for ( size_t cc=0; cc<c; ++cc ) {
753  std::size_t i = row*c + cc;
754  if ( !zerop(m[i]) ) {
755  return false;
756  }
757  }
758  return true;
759  }
760  void set_row(size_t row, const vector<cl_MI>& newrow)
761  {
762  for (std::size_t i2 = 0; i2 < newrow.size(); ++i2) {
763  std::size_t i1 = row*c + i2;
764  m[i1] = newrow[i2];
765  }
766  }
767  mvec::const_iterator row_begin(size_t row) const { return m.begin()+row*c; }
768  mvec::const_iterator row_end(size_t row) const { return m.begin()+row*c+r; }
769 private:
770  size_t r, c;
771  mvec m;
772 };
773 
774 #ifdef DEBUGFACTOR
775 modular_matrix operator*(const modular_matrix& m1, const modular_matrix& m2)
776 {
777  const unsigned int r = m1.rowsize();
778  const unsigned int c = m2.colsize();
779  modular_matrix o(r,c,m1(0,0));
780 
781  for ( size_t i=0; i<r; ++i ) {
782  for ( size_t j=0; j<c; ++j ) {
783  cl_MI buf;
784  buf = m1(i,0) * m2(0,j);
785  for ( size_t k=1; k<c; ++k ) {
786  buf = buf + m1(i,k)*m2(k,j);
787  }
788  o(i,j) = buf;
789  }
790  }
791  return o;
792 }
793 
794 ostream& operator<<(ostream& o, const modular_matrix& m)
795 {
796  cl_modint_ring R = m(0,0).ring();
797  o << "{";
798  for ( size_t i=0; i<m.rowsize(); ++i ) {
799  o << "{";
800  for ( size_t j=0; j<m.colsize()-1; ++j ) {
801  o << R->retract(m(i,j)) << ",";
802  }
803  o << R->retract(m(i,m.colsize()-1)) << "}";
804  if ( i != m.rowsize()-1 ) {
805  o << ",";
806  }
807  }
808  o << "}";
809  return o;
810 }
811 #endif // def DEBUGFACTOR
812 
813 // END modular matrix
815 
821 static void q_matrix(const umodpoly& a_, modular_matrix& Q)
822 {
823  umodpoly a = a_;
824  normalize_in_field(a);
825 
826  int n = degree(a);
827  unsigned int q = cl_I_to_uint(a[0].ring()->modulus);
828  umodpoly r(n, a[0].ring()->zero());
829  r[0] = a[0].ring()->one();
830  Q.set_row(0, r);
831  unsigned int max = (n-1) * q;
832  for ( size_t m=1; m<=max; ++m ) {
833  cl_MI rn_1 = r.back();
834  for ( size_t i=n-1; i>0; --i ) {
835  r[i] = r[i-1] - (rn_1 * a[i]);
836  }
837  r[0] = -rn_1 * a[0];
838  if ( (m % q) == 0 ) {
839  Q.set_row(m/q, r);
840  }
841  }
842 }
843 
849 static void nullspace(modular_matrix& M, vector<mvec>& basis)
850 {
851  const size_t n = M.rowsize();
852  const cl_MI one = M(0,0).ring()->one();
853  for ( size_t i=0; i<n; ++i ) {
854  M(i,i) = M(i,i) - one;
855  }
856  for ( size_t r=0; r<n; ++r ) {
857  size_t cc = 0;
858  for ( ; cc<n; ++cc ) {
859  if ( !zerop(M(r,cc)) ) {
860  if ( cc < r ) {
861  if ( !zerop(M(cc,cc)) ) {
862  continue;
863  }
864  M.switch_col(cc, r);
865  }
866  else if ( cc > r ) {
867  M.switch_col(cc, r);
868  }
869  break;
870  }
871  }
872  if ( cc < n ) {
873  M.mul_col(r, recip(M(r,r)));
874  for ( cc=0; cc<n; ++cc ) {
875  if ( cc != r ) {
876  M.sub_col(cc, r, M(r,cc));
877  }
878  }
879  }
880  }
881 
882  for ( size_t i=0; i<n; ++i ) {
883  M(i,i) = M(i,i) - one;
884  }
885  for ( size_t i=0; i<n; ++i ) {
886  if ( !M.is_row_zero(i) ) {
887  mvec nu(M.row_begin(i), M.row_end(i));
888  basis.push_back(nu);
889  }
890  }
891 }
892 
901 static void berlekamp(const umodpoly& a, upvec& upv)
902 {
903  cl_modint_ring R = a[0].ring();
904  umodpoly one(1, R->one());
905 
906  // find nullspace of Q matrix
907  modular_matrix Q(degree(a), degree(a), R->zero());
908  q_matrix(a, Q);
909  vector<mvec> nu;
910  nullspace(Q, nu);
911 
912  const unsigned int k = nu.size();
913  if ( k == 1 ) {
914  // irreducible
915  return;
916  }
917 
918  list<umodpoly> factors = {a};
919  unsigned int size = 1;
920  unsigned int r = 1;
921  unsigned int q = cl_I_to_uint(R->modulus);
922 
923  list<umodpoly>::iterator u = factors.begin();
924 
925  // calculate all gcd's
926  while ( true ) {
927  for ( unsigned int s=0; s<q; ++s ) {
928  umodpoly nur = nu[r];
929  nur[0] = nur[0] - cl_MI(R, s);
930  canonicalize(nur);
931  umodpoly g;
932  gcd(nur, *u, g);
933  if ( unequal_one(g) && g != *u ) {
934  umodpoly uo;
935  div(*u, g, uo);
936  if ( equal_one(uo) ) {
937  throw logic_error("berlekamp: unexpected divisor.");
938  } else {
939  *u = uo;
940  }
941  factors.push_back(g);
942  size = 0;
943  for (auto & i : factors) {
944  if (degree(i))
945  ++size;
946  }
947  if ( size == k ) {
948  for (auto & i : factors) {
949  upv.push_back(i);
950  }
951  return;
952  }
953  }
954  }
955  if ( ++r == k ) {
956  r = 1;
957  ++u;
958  }
959  }
960 }
961 
962 // modular square free factorization is not used at the moment so we deactivate
963 // the code
964 #if 0
965 
972 static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap)
973 {
974  size_t newdeg = degree(a)/prime;
975  ap.resize(newdeg+1);
976  ap[0] = a[0];
977  for ( size_t i=1; i<=newdeg; ++i ) {
978  ap[i] = a[i*prime];
979  }
980 }
981 
988 static void modsqrfree(const umodpoly& a, upvec& factors, vector<int>& mult)
989 {
990  const unsigned int prime = cl_I_to_uint(a[0].ring()->modulus);
991  int i = 1;
992  umodpoly b;
993  deriv(a, b);
994  if ( b.size() ) {
995  umodpoly c;
996  gcd(a, b, c);
997  umodpoly w;
998  div(a, c, w);
999  while ( unequal_one(w) ) {
1000  umodpoly y;
1001  gcd(w, c, y);
1002  umodpoly z;
1003  div(w, y, z);
1004  factors.push_back(z);
1005  mult.push_back(i);
1006  ++i;
1007  w = y;
1008  umodpoly buf;
1009  div(c, y, buf);
1010  c = buf;
1011  }
1012  if ( unequal_one(c) ) {
1013  umodpoly cp;
1014  expt_1_over_p(c, prime, cp);
1015  size_t previ = mult.size();
1016  modsqrfree(cp, factors, mult);
1017  for ( size_t i=previ; i<mult.size(); ++i ) {
1018  mult[i] *= prime;
1019  }
1020  }
1021  } else {
1022  umodpoly ap;
1023  expt_1_over_p(a, prime, ap);
1024  size_t previ = mult.size();
1025  modsqrfree(ap, factors, mult);
1026  for ( size_t i=previ; i<mult.size(); ++i ) {
1027  mult[i] *= prime;
1028  }
1029  }
1030 }
1031 
1032 #endif // deactivation of square free factorization
1033 
1044 static void distinct_degree_factor(const umodpoly& a_, vector<int>& degrees, upvec& ddfactors)
1045 {
1046  umodpoly a = a_;
1047 
1048  cl_modint_ring R = a[0].ring();
1049  int q = cl_I_to_int(R->modulus);
1050  int nhalf = degree(a)/2;
1051 
1052  int i = 1;
1053  umodpoly w(2);
1054  w[0] = R->zero();
1055  w[1] = R->one();
1056  umodpoly x = w;
1057 
1058  while ( i <= nhalf ) {
1059  expt_pos(w, q);
1060  umodpoly buf;
1061  rem(w, a, buf);
1062  w = buf;
1063  umodpoly wx = w - x;
1064  gcd(a, wx, buf);
1065  if ( unequal_one(buf) ) {
1066  degrees.push_back(i);
1067  ddfactors.push_back(buf);
1068  }
1069  if ( unequal_one(buf) ) {
1070  umodpoly buf2;
1071  div(a, buf, buf2);
1072  a = buf2;
1073  nhalf = degree(a)/2;
1074  rem(w, a, buf);
1075  w = buf;
1076  }
1077  ++i;
1078  }
1079  if ( unequal_one(a) ) {
1080  degrees.push_back(degree(a));
1081  ddfactors.push_back(a);
1082  }
1083 }
1084 
1095 static void same_degree_factor(const umodpoly& a, upvec& upv)
1096 {
1097  cl_modint_ring R = a[0].ring();
1098 
1099  vector<int> degrees;
1100  upvec ddfactors;
1101  distinct_degree_factor(a, degrees, ddfactors);
1102 
1103  for ( size_t i=0; i<degrees.size(); ++i ) {
1104  if ( degrees[i] == degree(ddfactors[i]) ) {
1105  upv.push_back(ddfactors[i]);
1106  } else {
1107  berlekamp(ddfactors[i], upv);
1108  }
1109  }
1110 }
1111 
1112 // Yes, we can (choose).
1113 #define USE_SAME_DEGREE_FACTOR
1114 
1125 static void factor_modular(const umodpoly& p, upvec& upv)
1126 {
1127 #ifdef USE_SAME_DEGREE_FACTOR
1128  same_degree_factor(p, upv);
1129 #else
1130  berlekamp(p, upv);
1131 #endif
1132 }
1133 
1142 static void exteuclid(const umodpoly& a, const umodpoly& b, umodpoly& s, umodpoly& t)
1143 {
1144  if ( degree(a) < degree(b) ) {
1145  exteuclid(b, a, t, s);
1146  return;
1147  }
1148 
1149  umodpoly one(1, a[0].ring()->one());
1150  umodpoly c = a; normalize_in_field(c);
1151  umodpoly d = b; normalize_in_field(d);
1152  s = one;
1153  t.clear();
1154  umodpoly d1;
1155  umodpoly d2 = one;
1156  umodpoly q;
1157  while ( true ) {
1158  div(c, d, q);
1159  umodpoly r = c - q * d;
1160  umodpoly r1 = s - q * d1;
1161  umodpoly r2 = t - q * d2;
1162  c = d;
1163  s = d1;
1164  t = d2;
1165  if ( r.empty() ) break;
1166  d = r;
1167  d1 = r1;
1168  d2 = r2;
1169  }
1170  cl_MI fac = recip(lcoeff(a) * lcoeff(c));
1171  for (auto & i : s) {
1172  i = i * fac;
1173  }
1174  canonicalize(s);
1175  fac = recip(lcoeff(b) * lcoeff(c));
1176  for (auto & i : t) {
1177  i = i * fac;
1178  }
1179  canonicalize(t);
1180 }
1181 
1188 static upoly replace_lc(const upoly& poly, const cl_I& lc)
1189 {
1190  if ( poly.empty() ) return poly;
1191  upoly r = poly;
1192  r.back() = lc;
1193  return r;
1194 }
1195 
1199 static inline cl_I calc_bound(const ex& a, const ex& x, int maxdeg)
1200 {
1201  cl_I maxcoeff = 0;
1202  cl_R coeff = 0;
1203  for ( int i=a.degree(x); i>=a.ldegree(x); --i ) {
1204  cl_I aa = abs(the<cl_I>(ex_to<numeric>(a.coeff(x, i)).to_cl_N()));
1205  if ( aa > maxcoeff ) maxcoeff = aa;
1206  coeff = coeff + square(aa);
1207  }
1208  cl_I coeffnorm = ceiling1(the<cl_R>(cln::sqrt(coeff)));
1209  cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg));
1210  return ( B > maxcoeff ) ? B : maxcoeff;
1211 }
1212 
1216 static inline cl_I calc_bound(const upoly& a, int maxdeg)
1217 {
1218  cl_I maxcoeff = 0;
1219  cl_R coeff = 0;
1220  for ( int i=degree(a); i>=0; --i ) {
1221  cl_I aa = abs(a[i]);
1222  if ( aa > maxcoeff ) maxcoeff = aa;
1223  coeff = coeff + square(aa);
1224  }
1225  cl_I coeffnorm = ceiling1(the<cl_R>(cln::sqrt(coeff)));
1226  cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg));
1227  return ( B > maxcoeff ) ? B : maxcoeff;
1228 }
1229 
1242 static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_, const umodpoly& w1_, upoly& u, upoly& w)
1243 {
1244  upoly a = a_;
1245  const cl_modint_ring& R = u1_[0].ring();
1246 
1247  // calc bound B
1248  int maxdeg = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_);
1249  cl_I maxmodulus = 2*calc_bound(a, maxdeg);
1250 
1251  // step 1
1252  cl_I alpha = lcoeff(a);
1253  a = a * alpha;
1254  umodpoly nu1 = u1_;
1255  normalize_in_field(nu1);
1256  umodpoly nw1 = w1_;
1257  normalize_in_field(nw1);
1258  upoly phi;
1259  phi = umodpoly_to_upoly(nu1) * alpha;
1260  umodpoly u1;
1261  umodpoly_from_upoly(u1, phi, R);
1262  phi = umodpoly_to_upoly(nw1) * alpha;
1263  umodpoly w1;
1264  umodpoly_from_upoly(w1, phi, R);
1265 
1266  // step 2
1267  umodpoly s;
1268  umodpoly t;
1269  exteuclid(u1, w1, s, t);
1270 
1271  // step 3
1272  u = replace_lc(umodpoly_to_upoly(u1), alpha);
1273  w = replace_lc(umodpoly_to_upoly(w1), alpha);
1274  upoly e = a - u * w;
1275  cl_I modulus = p;
1276 
1277  // step 4
1278  while ( !e.empty() && modulus < maxmodulus ) {
1279  upoly c = e / modulus;
1280  phi = umodpoly_to_upoly(s) * c;
1281  umodpoly sigmatilde;
1282  umodpoly_from_upoly(sigmatilde, phi, R);
1283  phi = umodpoly_to_upoly(t) * c;
1284  umodpoly tautilde;
1285  umodpoly_from_upoly(tautilde, phi, R);
1286  umodpoly r, q;
1287  remdiv(sigmatilde, w1, r, q);
1288  umodpoly sigma = r;
1289  phi = umodpoly_to_upoly(tautilde) + umodpoly_to_upoly(q) * umodpoly_to_upoly(u1);
1290  umodpoly tau;
1291  umodpoly_from_upoly(tau, phi, R);
1292  u = u + umodpoly_to_upoly(tau) * modulus;
1293  w = w + umodpoly_to_upoly(sigma) * modulus;
1294  e = a - u * w;
1295  modulus = modulus * p;
1296  }
1297 
1298  // step 5
1299  if ( e.empty() ) {
1300  cl_I g = u[0];
1301  for ( size_t i=1; i<u.size(); ++i ) {
1302  g = gcd(g, u[i]);
1303  if ( g == 1 ) break;
1304  }
1305  if ( g != 1 ) {
1306  u = u / g;
1307  w = w * g;
1308  }
1309  if ( alpha != 1 ) {
1310  w = w / alpha;
1311  }
1312  } else {
1313  u.clear();
1314  }
1315 }
1316 
1322 static unsigned int next_prime(unsigned int p)
1323 {
1324  static vector<unsigned int> primes;
1325  if (primes.empty()) {
1326  primes = {3, 5, 7};
1327  }
1328  if ( p >= primes.back() ) {
1329  unsigned int candidate = primes.back() + 2;
1330  while ( true ) {
1331  size_t n = primes.size()/2;
1332  for ( size_t i=0; i<n; ++i ) {
1333  if (candidate % primes[i])
1334  continue;
1335  candidate += 2;
1336  i=-1;
1337  }
1338  primes.push_back(candidate);
1339  if (candidate > p)
1340  break;
1341  }
1342  return candidate;
1343  }
1344  for (auto & it : primes) {
1345  if ( it > p ) {
1346  return it;
1347  }
1348  }
1349  throw logic_error("next_prime: should not reach this point!");
1350 }
1351 
1354 class factor_partition
1355 {
1356 public:
1358  factor_partition(const upvec& factors_) : factors(factors_)
1359  {
1360  n = factors.size();
1361  k.resize(n, 0);
1362  k[0] = 1;
1363  cache.resize(n-1);
1364  one.resize(1, factors.front()[0].ring()->one());
1365  len = 1;
1366  last = 0;
1367  split();
1368  }
1369  int operator[](size_t i) const { return k[i]; }
1370  size_t size() const { return n; }
1371  size_t size_left() const { return n-len; }
1372  size_t size_right() const { return len; }
1375  bool next()
1376  {
1377  if ( last == n-1 ) {
1378  int rem = len - 1;
1379  int p = last - 1;
1380  while ( rem ) {
1381  if ( k[p] ) {
1382  --rem;
1383  --p;
1384  continue;
1385  }
1386  last = p - 1;
1387  while ( k[last] == 0 ) { --last; }
1388  if ( last == 0 && n == 2*len ) return false;
1389  k[last++] = 0;
1390  for ( size_t i=0; i<=len-rem; ++i ) {
1391  k[last] = 1;
1392  ++last;
1393  }
1394  fill(k.begin()+last, k.end(), 0);
1395  --last;
1396  split();
1397  return true;
1398  }
1399  last = len;
1400  ++len;
1401  if ( len > n/2 ) return false;
1402  fill(k.begin(), k.begin()+len, 1);
1403  fill(k.begin()+len+1, k.end(), 0);
1404  } else {
1405  k[last++] = 0;
1406  k[last] = 1;
1407  }
1408  split();
1409  return true;
1410  }
1412  umodpoly& left() { return lr[0]; }
1414  umodpoly& right() { return lr[1]; }
1415 private:
1416  void split_cached()
1417  {
1418  size_t i = 0;
1419  do {
1420  size_t pos = i;
1421  int group = k[i++];
1422  size_t d = 0;
1423  while ( i < n && k[i] == group ) { ++d; ++i; }
1424  if ( d ) {
1425  if ( cache[pos].size() >= d ) {
1426  lr[group] = lr[group] * cache[pos][d-1];
1427  } else {
1428  if ( cache[pos].size() == 0 ) {
1429  cache[pos].push_back(factors[pos] * factors[pos+1]);
1430  }
1431  size_t j = pos + cache[pos].size() + 1;
1432  d -= cache[pos].size();
1433  while ( d ) {
1434  umodpoly buf = cache[pos].back() * factors[j];
1435  cache[pos].push_back(buf);
1436  --d;
1437  ++j;
1438  }
1439  lr[group] = lr[group] * cache[pos].back();
1440  }
1441  } else {
1442  lr[group] = lr[group] * factors[pos];
1443  }
1444  } while ( i < n );
1445  }
1446  void split()
1447  {
1448  lr[0] = one;
1449  lr[1] = one;
1450  if ( n > 6 ) {
1451  split_cached();
1452  } else {
1453  for ( size_t i=0; i<n; ++i ) {
1454  lr[k[i]] = lr[k[i]] * factors[i];
1455  }
1456  }
1457  }
1458 private:
1459  umodpoly lr[2];
1460  vector<vector<umodpoly>> cache;
1461  upvec factors;
1462  umodpoly one;
1463  size_t n;
1464  size_t len;
1465  size_t last;
1466  vector<int> k;
1467 };
1468 
1472 struct ModFactors
1473 {
1474  upoly poly;
1475  upvec factors;
1476 };
1477 
1488 static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
1489 {
1490  ex unit, cont, prim_ex;
1491  poly.unitcontprim(x, unit, cont, prim_ex);
1492  upoly prim;
1493  upoly_from_ex(prim, prim_ex, x);
1494  if (prim_ex.is_equal(1)) {
1495  return poly;
1496  }
1497 
1498  // determine proper prime and minimize number of modular factors
1499  prime = 3;
1500  unsigned int lastp = prime;
1501  cl_modint_ring R;
1502  unsigned int trials = 0;
1503  unsigned int minfactors = 0;
1504 
1505  const numeric& cont_n = ex_to<numeric>(cont);
1506  cl_I i_cont;
1507  if (cont_n.is_integer()) {
1508  i_cont = the<cl_I>(cont_n.to_cl_N());
1509  } else {
1510  // poly \in Q[x] => poly = q ipoly, ipoly \in Z[x], q \in Q
1511  // factor(poly) \equiv q factor(ipoly)
1512  i_cont = cl_I(1);
1513  }
1514  cl_I lc = lcoeff(prim)*i_cont;
1515  upvec factors;
1516  while ( trials < 2 ) {
1517  umodpoly modpoly;
1518  while ( true ) {
1519  prime = next_prime(prime);
1520  if ( !zerop(rem(lc, prime)) ) {
1521  R = find_modint_ring(prime);
1522  umodpoly_from_upoly(modpoly, prim, R);
1523  if ( squarefree(modpoly) ) break;
1524  }
1525  }
1526 
1527  // do modular factorization
1528  upvec trialfactors;
1529  factor_modular(modpoly, trialfactors);
1530  if ( trialfactors.size() <= 1 ) {
1531  // irreducible for sure
1532  return poly;
1533  }
1534 
1535  if ( minfactors == 0 || trialfactors.size() < minfactors ) {
1536  factors = trialfactors;
1537  minfactors = trialfactors.size();
1538  lastp = prime;
1539  trials = 1;
1540  } else {
1541  ++trials;
1542  }
1543  }
1544  prime = lastp;
1545  R = find_modint_ring(prime);
1546 
1547  // lift all factor combinations
1548  stack<ModFactors> tocheck;
1549  ModFactors mf;
1550  mf.poly = prim;
1551  mf.factors = factors;
1552  tocheck.push(mf);
1553  upoly f1, f2;
1554  ex result = 1;
1555  while ( tocheck.size() ) {
1556  const size_t n = tocheck.top().factors.size();
1557  factor_partition part(tocheck.top().factors);
1558  while ( true ) {
1559  // call Hensel lifting
1560  hensel_univar(tocheck.top().poly, prime, part.left(), part.right(), f1, f2);
1561  if ( !f1.empty() ) {
1562  // successful, update the stack and the result
1563  if ( part.size_left() == 1 ) {
1564  if ( part.size_right() == 1 ) {
1565  result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
1566  tocheck.pop();
1567  break;
1568  }
1569  result *= upoly_to_ex(f1, x);
1570  tocheck.top().poly = f2;
1571  for ( size_t i=0; i<n; ++i ) {
1572  if ( part[i] == 0 ) {
1573  tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1574  break;
1575  }
1576  }
1577  break;
1578  }
1579  else if ( part.size_right() == 1 ) {
1580  if ( part.size_left() == 1 ) {
1581  result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
1582  tocheck.pop();
1583  break;
1584  }
1585  result *= upoly_to_ex(f2, x);
1586  tocheck.top().poly = f1;
1587  for ( size_t i=0; i<n; ++i ) {
1588  if ( part[i] == 1 ) {
1589  tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1590  break;
1591  }
1592  }
1593  break;
1594  } else {
1595  upvec newfactors1(part.size_left()), newfactors2(part.size_right());
1596  auto i1 = newfactors1.begin(), i2 = newfactors2.begin();
1597  for ( size_t i=0; i<n; ++i ) {
1598  if ( part[i] ) {
1599  *i2++ = tocheck.top().factors[i];
1600  } else {
1601  *i1++ = tocheck.top().factors[i];
1602  }
1603  }
1604  tocheck.top().factors = newfactors1;
1605  tocheck.top().poly = f1;
1606  ModFactors mf;
1607  mf.factors = newfactors2;
1608  mf.poly = f2;
1609  tocheck.push(mf);
1610  break;
1611  }
1612  } else {
1613  // not successful
1614  if ( !part.next() ) {
1615  // if no more combinations left, return polynomial as
1616  // irreducible
1617  result *= upoly_to_ex(tocheck.top().poly, x);
1618  tocheck.pop();
1619  break;
1620  }
1621  }
1622  }
1623  }
1624 
1625  return unit * cont * result;
1626 }
1627 
1631 static inline ex factor_univariate(const ex& poly, const ex& x)
1632 {
1633  unsigned int prime;
1634  return factor_univariate(poly, x, prime);
1635 }
1636 
1639 struct EvalPoint
1640 {
1641  ex x;
1643 };
1644 
1645 #ifdef DEBUGFACTOR
1646 ostream& operator<<(ostream& o, const vector<EvalPoint>& v)
1647 {
1648  for ( size_t i=0; i<v.size(); ++i ) {
1649  o << "(" << v[i].x << "==" << v[i].evalpoint << ") ";
1650  }
1651  return o;
1652 }
1653 #endif // def DEBUGFACTOR
1654 
1655 // forward declaration
1656 static vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex& c, const vector<EvalPoint>& I, unsigned int d, unsigned int p, unsigned int k);
1657 
1673 static upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, unsigned int k)
1674 {
1675  const size_t r = a.size();
1676  cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1677  upvec q(r-1);
1678  q[r-2] = a[r-1];
1679  for ( size_t j=r-2; j>=1; --j ) {
1680  q[j-1] = a[j] * q[j];
1681  }
1682  umodpoly beta(1, R->one());
1683  upvec s;
1684  for ( size_t j=1; j<r; ++j ) {
1685  vector<ex> mdarg(2);
1686  mdarg[0] = umodpoly_to_ex(q[j-1], x);
1687  mdarg[1] = umodpoly_to_ex(a[j-1], x);
1688  vector<EvalPoint> empty;
1689  vector<ex> exsigma = multivar_diophant(mdarg, x, umodpoly_to_ex(beta, x), empty, 0, p, k);
1690  umodpoly sigma1;
1691  umodpoly_from_ex(sigma1, exsigma[0], x, R);
1692  umodpoly sigma2;
1693  umodpoly_from_ex(sigma2, exsigma[1], x, R);
1694  beta = sigma1;
1695  s.push_back(sigma2);
1696  }
1697  s.push_back(beta);
1698  return s;
1699 }
1700 
1706 static void change_modulus(const cl_modint_ring& R, umodpoly& a)
1707 {
1708  if ( a.empty() ) return;
1709  cl_modint_ring oldR = a[0].ring();
1710  for (auto & i : a) {
1711  i = R->canonhom(oldR->retract(i));
1712  }
1713  canonicalize(a);
1714 }
1715 
1730 static void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned int p, unsigned int k, umodpoly& s_, umodpoly& t_)
1731 {
1732  cl_modint_ring R = find_modint_ring(p);
1733  umodpoly amod = a;
1734  change_modulus(R, amod);
1735  umodpoly bmod = b;
1736  change_modulus(R, bmod);
1737 
1738  umodpoly smod;
1739  umodpoly tmod;
1740  exteuclid(amod, bmod, smod, tmod);
1741 
1742  cl_modint_ring Rpk = find_modint_ring(expt_pos(cl_I(p),k));
1743  umodpoly s = smod;
1744  change_modulus(Rpk, s);
1745  umodpoly t = tmod;
1746  change_modulus(Rpk, t);
1747 
1748  cl_I modulus(p);
1749  umodpoly one(1, Rpk->one());
1750  for ( size_t j=1; j<k; ++j ) {
1751  umodpoly e = one - a * s - b * t;
1752  reduce_coeff(e, modulus);
1753  umodpoly c = e;
1754  change_modulus(R, c);
1755  umodpoly sigmabar = smod * c;
1756  umodpoly taubar = tmod * c;
1757  umodpoly sigma, q;
1758  remdiv(sigmabar, bmod, sigma, q);
1759  umodpoly tau = taubar + q * amod;
1760  umodpoly sadd = sigma;
1761  change_modulus(Rpk, sadd);
1762  cl_MI modmodulus(Rpk, modulus);
1763  s = s + sadd * modmodulus;
1764  umodpoly tadd = tau;
1765  change_modulus(Rpk, tadd);
1766  t = t + tadd * modmodulus;
1767  modulus = modulus * p;
1768  }
1769 
1770  s_ = s; t_ = t;
1771 }
1772 
1788 static upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsigned int p, unsigned int k)
1789 {
1790  cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1791 
1792  const size_t r = a.size();
1793  upvec result;
1794  if ( r > 2 ) {
1795  upvec s = multiterm_eea_lift(a, x, p, k);
1796  for ( size_t j=0; j<r; ++j ) {
1797  umodpoly bmod = umodpoly_to_umodpoly(s[j], R, m);
1798  umodpoly buf;
1799  rem(bmod, a[j], buf);
1800  result.push_back(buf);
1801  }
1802  } else {
1803  umodpoly s, t;
1804  eea_lift(a[1], a[0], x, p, k, s, t);
1805  umodpoly bmod = umodpoly_to_umodpoly(s, R, m);
1806  umodpoly buf, q;
1807  remdiv(bmod, a[0], buf, q);
1808  result.push_back(buf);
1809  umodpoly t1mod = umodpoly_to_umodpoly(t, R, m);
1810  buf = t1mod + q * a[1];
1811  result.push_back(buf);
1812  }
1813 
1814  return result;
1815 }
1816 
1821 struct make_modular_map : public map_function {
1822  cl_modint_ring R;
1823  make_modular_map(const cl_modint_ring& R_) : R(R_) { }
1824  ex operator()(const ex& e) override
1825  {
1826  if ( is_a<add>(e) || is_a<mul>(e) ) {
1827  return e.map(*this);
1828  }
1829  else if ( is_a<numeric>(e) ) {
1830  numeric mod(R->modulus);
1831  numeric halfmod = (mod-1)/2;
1832  cl_MI emod = R->canonhom(the<cl_I>(ex_to<numeric>(e).to_cl_N()));
1833  numeric n(R->retract(emod));
1834  if ( n > halfmod ) {
1835  return n-mod;
1836  } else {
1837  return n;
1838  }
1839  }
1840  return e;
1841  }
1842 };
1843 
1851 static ex make_modular(const ex& e, const cl_modint_ring& R)
1852 {
1853  make_modular_map map(R);
1854  return map(e.expand());
1855 }
1856 
1874 static vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex& c, const vector<EvalPoint>& I,
1875  unsigned int d, unsigned int p, unsigned int k)
1876 {
1877  vector<ex> a = a_;
1878 
1879  const cl_I modulus = expt_pos(cl_I(p),k);
1880  const cl_modint_ring R = find_modint_ring(modulus);
1881  const size_t r = a.size();
1882  const size_t nu = I.size() + 1;
1883 
1884  vector<ex> sigma;
1885  if ( nu > 1 ) {
1886  ex xnu = I.back().x;
1887  int alphanu = I.back().evalpoint;
1888 
1889  ex A = 1;
1890  for ( size_t i=0; i<r; ++i ) {
1891  A *= a[i];
1892  }
1893  vector<ex> b(r);
1894  for ( size_t i=0; i<r; ++i ) {
1895  b[i] = normal(A / a[i]);
1896  }
1897 
1898  vector<ex> anew = a;
1899  for ( size_t i=0; i<r; ++i ) {
1900  anew[i] = anew[i].subs(xnu == alphanu);
1901  }
1902  ex cnew = c.subs(xnu == alphanu);
1903  vector<EvalPoint> Inew = I;
1904  Inew.pop_back();
1905  sigma = multivar_diophant(anew, x, cnew, Inew, d, p, k);
1906 
1907  ex buf = c;
1908  for ( size_t i=0; i<r; ++i ) {
1909  buf -= sigma[i] * b[i];
1910  }
1911  ex e = make_modular(buf, R);
1912 
1913  ex monomial = 1;
1914  for ( size_t m=1; !e.is_zero() && e.has(xnu) && m<=d; ++m ) {
1915  monomial *= (xnu - alphanu);
1916  monomial = expand(monomial);
1917  ex cm = e.diff(ex_to<symbol>(xnu), m).subs(xnu==alphanu) / factorial(m);
1918  cm = make_modular(cm, R);
1919  if ( !cm.is_zero() ) {
1920  vector<ex> delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k);
1921  ex buf = e;
1922  for ( size_t j=0; j<delta_s.size(); ++j ) {
1923  delta_s[j] *= monomial;
1924  sigma[j] += delta_s[j];
1925  buf -= delta_s[j] * b[j];
1926  }
1927  e = make_modular(buf, R);
1928  }
1929  }
1930  } else {
1931  upvec amod;
1932  for ( size_t i=0; i<a.size(); ++i ) {
1933  umodpoly up;
1934  umodpoly_from_ex(up, a[i], x, R);
1935  amod.push_back(up);
1936  }
1937 
1938  sigma.insert(sigma.begin(), r, 0);
1939  size_t nterms;
1940  ex z;
1941  if ( is_a<add>(c) ) {
1942  nterms = c.nops();
1943  z = c.op(0);
1944  } else {
1945  nterms = 1;
1946  z = c;
1947  }
1948  for ( size_t i=0; i<nterms; ++i ) {
1949  int m = z.degree(x);
1950  cl_I cm = the<cl_I>(ex_to<numeric>(z.lcoeff(x)).to_cl_N());
1951  upvec delta_s = univar_diophant(amod, x, m, p, k);
1952  cl_MI modcm;
1953  cl_I poscm = plusp(cm) ? cm : mod(cm, modulus);
1954  modcm = cl_MI(R, poscm);
1955  for ( size_t j=0; j<delta_s.size(); ++j ) {
1956  delta_s[j] = delta_s[j] * modcm;
1957  sigma[j] = sigma[j] + umodpoly_to_ex(delta_s[j], x);
1958  }
1959  if ( nterms > 1 && i+1 != nterms ) {
1960  z = c.op(i+1);
1961  }
1962  }
1963  }
1964 
1965  for ( size_t i=0; i<sigma.size(); ++i ) {
1966  sigma[i] = make_modular(sigma[i], R);
1967  }
1968 
1969  return sigma;
1970 }
1971 
1988 static ex hensel_multivar(const ex& a, const ex& x, const vector<EvalPoint>& I,
1989  unsigned int p, const cl_I& l, const upvec& u, const vector<ex>& lcU)
1990 {
1991  const size_t nu = I.size() + 1;
1992  const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),l));
1993 
1994  vector<ex> A(nu);
1995  A[nu-1] = a;
1996 
1997  for ( size_t j=nu; j>=2; --j ) {
1998  ex x = I[j-2].x;
1999  int alpha = I[j-2].evalpoint;
2000  A[j-2] = A[j-1].subs(x==alpha);
2001  A[j-2] = make_modular(A[j-2], R);
2002  }
2003 
2004  int maxdeg = a.degree(I.front().x);
2005  for ( size_t i=1; i<I.size(); ++i ) {
2006  int maxdeg2 = a.degree(I[i].x);
2007  if ( maxdeg2 > maxdeg ) maxdeg = maxdeg2;
2008  }
2009 
2010  const size_t n = u.size();
2011  vector<ex> U(n);
2012  for ( size_t i=0; i<n; ++i ) {
2013  U[i] = umodpoly_to_ex(u[i], x);
2014  }
2015 
2016  for ( size_t j=2; j<=nu; ++j ) {
2017  vector<ex> U1 = U;
2018  ex monomial = 1;
2019  for ( size_t m=0; m<n; ++m) {
2020  if ( lcU[m] != 1 ) {
2021  ex coef = lcU[m];
2022  for ( size_t i=j-1; i<nu-1; ++i ) {
2023  coef = coef.subs(I[i].x == I[i].evalpoint);
2024  }
2025  coef = make_modular(coef, R);
2026  int deg = U[m].degree(x);
2027  U[m] = U[m] - U[m].lcoeff(x) * pow(x,deg) + coef * pow(x,deg);
2028  }
2029  }
2030  ex Uprod = 1;
2031  for ( size_t i=0; i<n; ++i ) {
2032  Uprod *= U[i];
2033  }
2034  ex e = expand(A[j-1] - Uprod);
2035 
2036  vector<EvalPoint> newI;
2037  for ( size_t i=1; i<=j-2; ++i ) {
2038  newI.push_back(I[i-1]);
2039  }
2040 
2041  ex xj = I[j-2].x;
2042  int alphaj = I[j-2].evalpoint;
2043  size_t deg = A[j-1].degree(xj);
2044  for ( size_t k=1; k<=deg; ++k ) {
2045  if ( !e.is_zero() ) {
2046  monomial *= (xj - alphaj);
2047  monomial = expand(monomial);
2048  ex dif = e.diff(ex_to<symbol>(xj), k);
2049  ex c = dif.subs(xj==alphaj) / factorial(k);
2050  if ( !c.is_zero() ) {
2051  vector<ex> deltaU = multivar_diophant(U1, x, c, newI, maxdeg, p, cl_I_to_uint(l));
2052  for ( size_t i=0; i<n; ++i ) {
2053  deltaU[i] *= monomial;
2054  U[i] += deltaU[i];
2055  U[i] = make_modular(U[i], R);
2056  }
2057  ex Uprod = 1;
2058  for ( size_t i=0; i<n; ++i ) {
2059  Uprod *= U[i];
2060  }
2061  e = A[j-1] - Uprod;
2062  e = make_modular(e, R);
2063  }
2064  }
2065  }
2066  }
2067 
2068  ex acand = 1;
2069  for ( size_t i=0; i<U.size(); ++i ) {
2070  acand *= U[i];
2071  }
2072  if ( expand(a-acand).is_zero() ) {
2073  return lst(U.begin(), U.end());
2074  } else {
2075  return lst{};
2076  }
2077 }
2078 
2083 static ex put_factors_into_lst(const ex& e)
2084 {
2085  lst result;
2086  if ( is_a<numeric>(e) ) {
2087  result.append(e);
2088  return result;
2089  }
2090  if ( is_a<power>(e) ) {
2091  result.append(1);
2092  result.append(e.op(0));
2093  return result;
2094  }
2095  if ( is_a<symbol>(e) || is_a<add>(e) ) {
2096  ex icont(e.integer_content());
2097  result.append(icont);
2098  result.append(e/icont);
2099  return result;
2100  }
2101  if ( is_a<mul>(e) ) {
2102  ex nfac = 1;
2103  for ( size_t i=0; i<e.nops(); ++i ) {
2104  ex op = e.op(i);
2105  if ( is_a<numeric>(op) ) {
2106  nfac = op;
2107  }
2108  if ( is_a<power>(op) ) {
2109  result.append(op.op(0));
2110  }
2111  if ( is_a<symbol>(op) || is_a<add>(op) ) {
2112  result.append(op);
2113  }
2114  }
2115  result.prepend(nfac);
2116  return result;
2117  }
2118  throw runtime_error("put_factors_into_lst: bad term.");
2119 }
2120 
2127 static bool checkdivisors(const lst& f)
2128 {
2129  const int k = f.nops();
2130  numeric q, r;
2131  vector<numeric> d(k);
2132  d[0] = ex_to<numeric>(abs(f.op(0)));
2133  for ( int i=1; i<k; ++i ) {
2134  q = ex_to<numeric>(abs(f.op(i)));
2135  for ( int j=i-1; j>=0; --j ) {
2136  r = d[j];
2137  do {
2138  r = gcd(r, q);
2139  q = q/r;
2140  } while ( r != 1 );
2141  if ( q == 1 ) {
2142  return true;
2143  }
2144  }
2145  d[i] = q;
2146  }
2147  return false;
2148 }
2149 
2166 static void generate_set(const ex& u, const ex& vn, const exset& syms, const lst& f,
2167  numeric& modulus, ex& u0, vector<numeric>& a)
2168 {
2169  const ex& x = *syms.begin();
2170  while ( true ) {
2171  ++modulus;
2172  // generate a set of integers ...
2173  u0 = u;
2174  ex vna = vn;
2175  ex vnatry;
2176  exset::const_iterator s = syms.begin();
2177  ++s;
2178  for ( size_t i=0; i<a.size(); ++i ) {
2179  do {
2180  a[i] = mod(numeric(rand()), 2*modulus) - modulus;
2181  vnatry = vna.subs(*s == a[i]);
2182  // ... for which the leading coefficient doesn't vanish ...
2183  } while ( vnatry == 0 );
2184  vna = vnatry;
2185  u0 = u0.subs(*s == a[i]);
2186  ++s;
2187  }
2188  // ... for which u0 is square free ...
2189  ex g = gcd(u0, u0.diff(ex_to<symbol>(x)));
2190  if ( !is_a<numeric>(g) ) {
2191  continue;
2192  }
2193  if ( !is_a<numeric>(vn) ) {
2194  // ... and for which the evaluated factors have each an unique prime factor
2195  lst fnum = f;
2196  fnum.let_op(0) = fnum.op(0) * u0.content(x);
2197  for ( size_t i=1; i<fnum.nops(); ++i ) {
2198  if ( !is_a<numeric>(fnum.op(i)) ) {
2199  s = syms.begin();
2200  ++s;
2201  for ( size_t j=0; j<a.size(); ++j, ++s ) {
2202  fnum.let_op(i) = fnum.op(i).subs(*s == a[j]);
2203  }
2204  }
2205  }
2206  if ( checkdivisors(fnum) ) {
2207  continue;
2208  }
2209  }
2210  // ok, we have a valid set now
2211  return;
2212  }
2213 }
2214 
2215 // forward declaration
2216 static ex factor_sqrfree(const ex& poly);
2217 
2233 static ex factor_multivariate(const ex& poly, const exset& syms)
2234 {
2235  const ex& x = *syms.begin();
2236 
2237  // make polynomial primitive
2238  ex unit, cont, pp;
2239  poly.unitcontprim(x, unit, cont, pp);
2240  if ( !is_a<numeric>(cont) ) {
2241  return unit * factor_sqrfree(cont) * factor_sqrfree(pp);
2242  }
2243 
2244  // factor leading coefficient
2245  ex vn = pp.collect(x).lcoeff(x);
2246  ex vnlst;
2247  if ( is_a<numeric>(vn) ) {
2248  vnlst = lst{vn};
2249  }
2250  else {
2251  ex vnfactors = factor(vn);
2252  vnlst = put_factors_into_lst(vnfactors);
2253  }
2254 
2255  const unsigned int maxtrials = 3;
2256  numeric modulus = (vnlst.nops() > 3) ? vnlst.nops() : 3;
2257  vector<numeric> a(syms.size()-1, 0);
2258 
2259  // try now to factorize until we are successful
2260  while ( true ) {
2261 
2262  unsigned int trialcount = 0;
2263  unsigned int prime;
2264  int factor_count = 0;
2265  int min_factor_count = -1;
2266  ex u, delta;
2267  ex ufac, ufaclst;
2268 
2269  // try several evaluation points to reduce the number of factors
2270  while ( trialcount < maxtrials ) {
2271 
2272  // generate a set of valid evaluation points
2273  generate_set(pp, vn, syms, ex_to<lst>(vnlst), modulus, u, a);
2274 
2275  ufac = factor_univariate(u, x, prime);
2276  ufaclst = put_factors_into_lst(ufac);
2277  factor_count = ufaclst.nops()-1;
2278  delta = ufaclst.op(0);
2279 
2280  if ( factor_count <= 1 ) {
2281  // irreducible
2282  return poly;
2283  }
2284  if ( min_factor_count < 0 ) {
2285  // first time here
2286  min_factor_count = factor_count;
2287  }
2288  else if ( min_factor_count == factor_count ) {
2289  // one less to try
2290  ++trialcount;
2291  }
2292  else if ( min_factor_count > factor_count ) {
2293  // new minimum, reset trial counter
2294  min_factor_count = factor_count;
2295  trialcount = 0;
2296  }
2297  }
2298 
2299  // determine true leading coefficients for the Hensel lifting
2300  vector<ex> C(factor_count);
2301  if ( is_a<numeric>(vn) ) {
2302  // easy case
2303  for ( size_t i=1; i<ufaclst.nops(); ++i ) {
2304  C[i-1] = ufaclst.op(i).lcoeff(x);
2305  }
2306  } else {
2307  // difficult case.
2308  // we use the property of the ftilde having a unique prime factor.
2309  // details can be found in [Wan].
2310  // calculate ftilde
2311  vector<numeric> ftilde(vnlst.nops()-1);
2312  for ( size_t i=0; i<ftilde.size(); ++i ) {
2313  ex ft = vnlst.op(i+1);
2314  auto s = syms.begin();
2315  ++s;
2316  for ( size_t j=0; j<a.size(); ++j ) {
2317  ft = ft.subs(*s == a[j]);
2318  ++s;
2319  }
2320  ftilde[i] = ex_to<numeric>(ft);
2321  }
2322  // calculate D and C
2323  vector<bool> used_flag(ftilde.size(), false);
2324  vector<ex> D(factor_count, 1);
2325  if ( delta == 1 ) {
2326  for ( int i=0; i<factor_count; ++i ) {
2327  numeric prefac = ex_to<numeric>(ufaclst.op(i+1).lcoeff(x));
2328  for ( int j=ftilde.size()-1; j>=0; --j ) {
2329  int count = 0;
2330  while ( irem(prefac, ftilde[j]) == 0 ) {
2331  prefac = iquo(prefac, ftilde[j]);
2332  ++count;
2333  }
2334  if ( count ) {
2335  used_flag[j] = true;
2336  D[i] = D[i] * pow(vnlst.op(j+1), count);
2337  }
2338  }
2339  C[i] = D[i] * prefac;
2340  }
2341  } else {
2342  for ( int i=0; i<factor_count; ++i ) {
2343  numeric prefac = ex_to<numeric>(ufaclst.op(i+1).lcoeff(x));
2344  for ( int j=ftilde.size()-1; j>=0; --j ) {
2345  int count = 0;
2346  while ( irem(prefac, ftilde[j]) == 0 ) {
2347  prefac = iquo(prefac, ftilde[j]);
2348  ++count;
2349  }
2350  while ( irem(ex_to<numeric>(delta)*prefac, ftilde[j]) == 0 ) {
2351  numeric g = gcd(prefac, ex_to<numeric>(ftilde[j]));
2352  prefac = iquo(prefac, g);
2353  delta = delta / (ftilde[j]/g);
2354  ufaclst.let_op(i+1) = ufaclst.op(i+1) * (ftilde[j]/g);
2355  ++count;
2356  }
2357  if ( count ) {
2358  used_flag[j] = true;
2359  D[i] = D[i] * pow(vnlst.op(j+1), count);
2360  }
2361  }
2362  C[i] = D[i] * prefac;
2363  }
2364  }
2365  // check if something went wrong
2366  bool some_factor_unused = false;
2367  for ( size_t i=0; i<used_flag.size(); ++i ) {
2368  if ( !used_flag[i] ) {
2369  some_factor_unused = true;
2370  break;
2371  }
2372  }
2373  if ( some_factor_unused ) {
2374  continue;
2375  }
2376  }
2377 
2378  // multiply the remaining content of the univariate polynomial into the
2379  // first factor
2380  if ( delta != 1 ) {
2381  C[0] = C[0] * delta;
2382  ufaclst.let_op(1) = ufaclst.op(1) * delta;
2383  }
2384 
2385  // set up evaluation points
2386  EvalPoint ep;
2387  vector<EvalPoint> epv;
2388  auto s = syms.begin();
2389  ++s;
2390  for ( size_t i=0; i<a.size(); ++i ) {
2391  ep.x = *s++;
2392  ep.evalpoint = a[i].to_int();
2393  epv.push_back(ep);
2394  }
2395 
2396  // calc bound p^l
2397  int maxdeg = 0;
2398  for ( int i=1; i<=factor_count; ++i ) {
2399  if ( ufaclst.op(i).degree(x) > maxdeg ) {
2400  maxdeg = ufaclst[i].degree(x);
2401  }
2402  }
2403  cl_I B = 2*calc_bound(u, x, maxdeg);
2404  cl_I l = 1;
2405  cl_I pl = prime;
2406  while ( pl < B ) {
2407  l = l + 1;
2408  pl = pl * prime;
2409  }
2410 
2411  // set up modular factors (mod p^l)
2412  cl_modint_ring R = find_modint_ring(expt_pos(cl_I(prime),l));
2413  upvec modfactors(ufaclst.nops()-1);
2414  for ( size_t i=1; i<ufaclst.nops(); ++i ) {
2415  umodpoly_from_ex(modfactors[i-1], ufaclst.op(i), x, R);
2416  }
2417 
2418  // try Hensel lifting
2419  ex res = hensel_multivar(pp, x, epv, prime, l, modfactors, C);
2420  if ( res != lst{} ) {
2421  ex result = cont * unit;
2422  for ( size_t i=0; i<res.nops(); ++i ) {
2423  result *= res.op(i).content(x) * res.op(i).unit(x);
2424  result *= res.op(i).primpart(x);
2425  }
2426  return result;
2427  }
2428  }
2429 }
2430 
2433 struct find_symbols_map : public map_function {
2435  ex operator()(const ex& e) override
2436  {
2437  if ( is_a<symbol>(e) ) {
2438  syms.insert(e);
2439  return e;
2440  }
2441  return e.map(*this);
2442  }
2443 };
2444 
2448 static ex factor_sqrfree(const ex& poly)
2449 {
2450  // determine all symbols in poly
2451  find_symbols_map findsymbols;
2452  findsymbols(poly);
2453  if ( findsymbols.syms.size() == 0 ) {
2454  return poly;
2455  }
2456 
2457  if ( findsymbols.syms.size() == 1 ) {
2458  // univariate case
2459  const ex& x = *(findsymbols.syms.begin());
2460  int ld = poly.ldegree(x);
2461  if ( ld > 0 ) {
2462  // pull out direct factors
2463  ex res = factor_univariate(expand(poly/pow(x, ld)), x);
2464  return res * pow(x,ld);
2465  } else {
2466  ex res = factor_univariate(poly, x);
2467  return res;
2468  }
2469  }
2470 
2471  // multivariate case
2472  ex res = factor_multivariate(poly, findsymbols.syms);
2473  return res;
2474 }
2475 
2479 struct apply_factor_map : public map_function {
2480  unsigned options;
2481  apply_factor_map(unsigned options_) : options(options_) { }
2482  ex operator()(const ex& e) override
2483  {
2484  if ( e.info(info_flags::polynomial) ) {
2485  return factor(e, options);
2486  }
2487  if ( is_a<add>(e) ) {
2488  ex s1, s2;
2489  for ( size_t i=0; i<e.nops(); ++i ) {
2490  if ( e.op(i).info(info_flags::polynomial) ) {
2491  s1 += e.op(i);
2492  } else {
2493  s2 += e.op(i);
2494  }
2495  }
2496  return factor(s1, options) + s2.map(*this);
2497  }
2498  return e.map(*this);
2499  }
2500 };
2501 
2508 template <typename F> void
2509 factor_iter(const ex &e, F yield)
2510 {
2511  if (is_a<mul>(e)) {
2512  for (const auto &f : e) {
2513  if (is_a<power>(f)) {
2514  yield(f.op(0), f.op(1));
2515  } else {
2516  yield(f, ex(1));
2517  }
2518  }
2519  } else {
2520  if (is_a<power>(e)) {
2521  yield(e.op(0), e.op(1));
2522  } else {
2523  yield(e, ex(1));
2524  }
2525  }
2526 }
2527 
2536 static ex factor1(const ex& poly, unsigned options)
2537 {
2538  // check arguments
2539  if ( !poly.info(info_flags::polynomial) ) {
2540  if ( options & factor_options::all ) {
2541  options &= ~factor_options::all;
2542  apply_factor_map factor_map(options);
2543  return factor_map(poly);
2544  }
2545  return poly;
2546  }
2547 
2548  // determine all symbols in poly
2549  find_symbols_map findsymbols;
2550  findsymbols(poly);
2551  if ( findsymbols.syms.size() == 0 ) {
2552  return poly;
2553  }
2554  lst syms;
2555  for (auto & i : findsymbols.syms ) {
2556  syms.append(i);
2557  }
2558 
2559  // make poly square free
2560  ex sfpoly = sqrfree(poly.expand(), syms);
2561 
2562  // factorize the square free components
2563  ex res = 1;
2564  factor_iter(sfpoly,
2565  [&](const ex &f, const ex &k) {
2566  if ( is_a<add>(f) ) {
2567  res *= pow(factor_sqrfree(f), k);
2568  } else {
2569  // simple case: (monomial)^exponent
2570  res *= pow(f, k);
2571  }
2572  });
2573  return res;
2574 }
2575 
2576 } // anonymous namespace
2577 
2581 ex factor(const ex& poly, unsigned options)
2582 {
2583  ex result = 1;
2584  factor_iter(poly,
2585  [&](const ex &f1, const ex &k1) {
2586  factor_iter(factor1(f1, options),
2587  [&](const ex &f2, const ex &k2) {
2588  result *= pow(f2, k1*k2);
2589  });
2590  });
2591  return result;
2592 }
2593 
2594 } // namespace GiNaC
Interface to GiNaC's sums of expressions.
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
ex map(map_function &f) const
Definition: ex.h:162
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:826
ex op(size_t i) const
Definition: ex.h:136
ex subs(const exmap &m, unsigned options=0) const override
Substitute a set of objects by arbitrary expressions.
Definition: numeric.h:110
int degree(const ex &s) const override
Return degree of highest power in object s.
Definition: numeric.cpp:733
Interface to GiNaC's light-weight expression handles.
vector< vector< umodpoly > > cache
Definition: factor.cpp:1460
static const bool value
Definition: factor.cpp:231
upoly poly
Definition: factor.cpp:1474
upvec factors
Definition: factor.cpp:1461
int evalpoint
Definition: factor.cpp:1642
unsigned options
Definition: factor.cpp:2480
vector< int > k
Definition: factor.cpp:1466
ex x
Definition: factor.cpp:1641
size_t n
Definition: factor.cpp:1463
size_t len
Definition: factor.cpp:1464
size_t c
Definition: factor.cpp:770
size_t r
Definition: factor.cpp:770
exset syms
Definition: factor.cpp:2434
umodpoly one
Definition: factor.cpp:1462
umodpoly lr[2]
Definition: factor.cpp:1459
cl_modint_ring R
Definition: factor.cpp:1822
mvec m
Definition: factor.cpp:771
size_t last
Definition: factor.cpp:1465
Polynomial factorization.
Interface to GiNaC's initially known functions.
Interface to GiNaC's products of expressions.
Definition: add.cpp:38
bool is_zero(const ex &thisex)
Definition: ex.h:820
const numeric I
Imaginary unit.
Definition: numeric.cpp:1433
const numeric pow(const numeric &x, const numeric &y)
Definition: numeric.h:251
container< std::list > lst
Definition: lst.h:32
ex sqrfree(const ex &a, const lst &l)
Compute a square-free factorization of a multivariate polynomial in Q[X].
Definition: normal.cpp:1888
std::set< ex, ex_is_less > exset
Definition: basic.h:49
const numeric mod(const numeric &a, const numeric &b)
Modulus (in positive representation).
Definition: numeric.cpp:2328
const ex operator/(const ex &lh, const ex &rh)
Definition: operators.cpp:80
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2315
ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned options)
Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X) and b(X) in Z[X].
Definition: normal.cpp:1432
const numeric sqrt(const numeric &x)
Numeric square root.
Definition: numeric.cpp:2475
const numeric irem(const numeric &a, const numeric &b)
Numeric integer remainder.
Definition: numeric.cpp:2363
const ex operator+(const ex &lh, const ex &rh)
Definition: operators.cpp:65
const ex operator*(const ex &lh, const ex &rh)
Definition: operators.cpp:75
const numeric factorial(const numeric &n)
Factorial combinatorial function.
Definition: numeric.cpp:2113
const numeric smod(const numeric &a_, const numeric &b_)
Modulus (in symmetric representation).
Definition: numeric.cpp:2341
int degree(const ex &thisex, const ex &s)
Definition: ex.h:736
std::ostream & operator<<(std::ostream &os, const archive_node &n)
Write archive_node to binary data stream.
Definition: archive.cpp:200
const numeric iquo(const numeric &a, const numeric &b)
Numeric integer quotient.
Definition: numeric.cpp:2404
ex normal(const ex &thisex)
Definition: ex.h:754
ex op(const ex &thisex, size_t i)
Definition: ex.h:811
ex coeff(const ex &thisex, const ex &s, int n=1)
Definition: ex.h:742
ex factor(const ex &poly, unsigned options)
Interface function to the outside world.
Definition: factor.cpp:2581
int canonicalize(exvector::iterator v, const symmetry &symm)
Canonicalize the order of elements of an expression vector, according to the symmetry properties defi...
Definition: symmetry.cpp:438
ex rem(const ex &a, const ex &b, const ex &x, bool check_args)
Remainder r(x) of polynomials a(x) and b(x) in Q[x].
Definition: normal.cpp:423
const ex operator-(const ex &lh, const ex &rh)
Definition: operators.cpp:70
ex expand(const ex &thisex, unsigned options=0)
Definition: ex.h:715
Definition: ex.h:972
void swap(GiNaC::ex &a, GiNaC::ex &b)
Specialization of std::swap() for ex objects.
Definition: ex.h:976
This file defines several functions that work on univariate and multivariate polynomials and rational...
Makes the interface to the underlying bignum package available.
Interface to GiNaC's overloaded operators.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
Interface to relations between expressions.
Interface to GiNaC's symbolic objects.

This page is part of the GiNaC developer's reference. It was generated automatically by doxygen. For an introduction, see the tutorial.