GiNaC  1.8.0
factor.cpp
Go to the documentation of this file.
1 
35 /*
36  * GiNaC Copyright (C) 1999-2020 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
inifcns.h
Interface to GiNaC's initially known functions.
GiNaC::lst
container< std::list > lst
Definition: lst.h:32
r
size_t r
Definition: factor.cpp:770
mul.h
Interface to GiNaC's products of expressions.
len
size_t len
Definition: factor.cpp:1464
GiNaC::iquo
const numeric iquo(const numeric &a, const numeric &b)
Numeric integer quotient.
Definition: numeric.cpp:2404
numeric.h
Makes the interface to the underlying bignum package available.
R
cl_modint_ring R
Definition: factor.cpp:1822
GiNaC::ex::subs
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:826
add.h
Interface to GiNaC's sums of expressions.
k
vector< int > k
Definition: factor.cpp:1466
power.h
Interface to GiNaC's symbolic exponentiation (basis^exponent).
GiNaC::sqrfree
ex sqrfree(const ex &a, const lst &l)
Compute a square-free factorization of a multivariate polynomial in Q[X].
Definition: normal.cpp:1888
GiNaC::operator/
const ex operator/(const ex &lh, const ex &rh)
Definition: operators.cpp:80
poly
upoly poly
Definition: factor.cpp:1474
GiNaC::canonicalize
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
relational.h
Interface to relations between expressions.
options
unsigned options
Definition: factor.cpp:2480
factors
upvec factors
Definition: factor.cpp:1461
m
mvec m
Definition: factor.cpp:771
GiNaC::factorial
const numeric factorial(const numeric &n)
Factorial combinatorial function.
Definition: numeric.cpp:2113
GiNaC::sqrt
const numeric sqrt(const numeric &x)
Numeric square root.
Definition: numeric.cpp:2475
GiNaC
Definition: add.cpp:38
syms
exset syms
Definition: factor.cpp:2434
GiNaC::operator<<
std::ostream & operator<<(std::ostream &os, const archive_node &n)
Write archive_node to binary data stream.
Definition: archive.cpp:200
factor.h
Polynomial factorization.
x
ex x
Definition: factor.cpp:1641
lr
umodpoly lr[2]
Definition: factor.cpp:1459
GiNaC::ex::op
ex op(size_t i) const
Definition: ex.h:136
GiNaC::numeric::subs
ex subs(const exmap &m, unsigned options=0) const override
Substitute a set of objects by arbitrary expressions.
Definition: numeric.h:110
last
size_t last
Definition: factor.cpp:1465
GiNaC::gcd
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
GiNaC::ex
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
GiNaC::normal
ex normal(const ex &thisex)
Definition: ex.h:754
value
static const bool value
Definition: factor.cpp:231
GiNaC::ex::map
ex map(map_function &f) const
Definition: ex.h:162
GiNaC::rem
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
GiNaC::expand
ex expand(const ex &thisex, unsigned options=0)
Definition: ex.h:715
ex.h
Interface to GiNaC's light-weight expression handles.
GiNaC::I
const numeric I
Imaginary unit.
Definition: numeric.cpp:1433
symbol.h
Interface to GiNaC's symbolic objects.
GiNaC::operator-
const ex operator-(const ex &lh, const ex &rh)
Definition: operators.cpp:70
GiNaC::exset
std::set< ex, ex_is_less > exset
Definition: basic.h:49
GiNaC::op
ex op(const ex &thisex, size_t i)
Definition: ex.h:811
GiNaC::irem
const numeric irem(const numeric &a, const numeric &b)
Numeric integer remainder.
Definition: numeric.cpp:2363
c
size_t c
Definition: factor.cpp:770
GiNaC::coeff
ex coeff(const ex &thisex, const ex &s, int n=1)
Definition: ex.h:742
GiNaC::factor
ex factor(const ex &poly, unsigned options)
Interface function to the outside world.
Definition: factor.cpp:2581
GiNaC::operator*
const ex operator*(const ex &lh, const ex &rh)
Definition: operators.cpp:75
GiNaC::mod
const numeric mod(const numeric &a, const numeric &b)
Modulus (in positive representation).
Definition: numeric.cpp:2328
GiNaC::degree
int degree(const ex &thisex, const ex &s)
Definition: ex.h:736
std
Definition: ex.h:972
GiNaC::numeric::degree
int degree(const ex &s) const override
Return degree of highest power in object s.
Definition: numeric.cpp:733
n
size_t n
Definition: factor.cpp:1463
GiNaC::pow
const numeric pow(const numeric &x, const numeric &y)
Definition: numeric.h:251
GiNaC::is_zero
bool is_zero(const ex &thisex)
Definition: ex.h:820
std::swap
void swap(GiNaC::ex &a, GiNaC::ex &b)
Specialization of std::swap() for ex objects.
Definition: ex.h:976
cache
vector< vector< umodpoly > > cache
Definition: factor.cpp:1460
GiNaC::abs
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2315
GiNaC::smod
const numeric smod(const numeric &a_, const numeric &b_)
Modulus (in symmetric representation).
Definition: numeric.cpp:2341
one
umodpoly one
Definition: factor.cpp:1462
evalpoint
int evalpoint
Definition: factor.cpp:1642
operators.h
Interface to GiNaC's overloaded operators.
normal.h
This file defines several functions that work on univariate and multivariate polynomials and rational...
GiNaC::operator+
const ex operator+(const ex &lh, const ex &rh)
Definition: operators.cpp:65

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