ginac/factor.cpp

   1 /** @file factor.cpp
   2  *
   3  *  Polynomial factorization code (implementation).
   4  *
   5  *  Algorithms used can be found in
   6  *    [W1]  An Improved Multivariate Polynomial Factoring Algorithm,
   7  *          P.S.Wang, Mathematics of Computation, Vol. 32, No. 144 (1978) 1215--1231.
   8  *    [GCL] Algorithms for Computer Algebra,
   9  *          K.O.Geddes, S.R.Czapor, G.Labahn, Springer Verlag, 1992.
  10  */
  11
  12 /*
  13  *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
  14  *
  15  *  This program is free software; you can redistribute it and/or modify
  16  *  it under the terms of the GNU General Public License as published by
  17  *  the Free Software Foundation; either version 2 of the License, or
  18  *  (at your option) any later version.
  19  *
  20  *  This program is distributed in the hope that it will be useful,
  21  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  22  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  23  *  GNU General Public License for more details.
  24  *
  25  *  You should have received a copy of the GNU General Public License
  26  *  along with this program; if not, write to the Free Software
  27  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  28  */
  29
  30 //#define DEBUGFACTOR
  31
  32 #include "factor.h"
  33
  34 #include "ex.h"
  35 #include "numeric.h"
  36 #include "operators.h"
  37 #include "inifcns.h"
  38 #include "symbol.h"
  39 #include "relational.h"
  40 #include "power.h"
  41 #include "mul.h"
  42 #include "normal.h"
  43 #include "add.h"
  44
  45 #include <algorithm>
  46 #include <cmath>
  47 #include <limits>
  48 #include <list>
  49 #include <vector>
  50 #ifdef DEBUGFACTOR
  51 #include <ostream>
  52 #endif
  53 using namespace std;
  54
  55 #include <cln/cln.h>
  56 using namespace cln;
  57
  58 namespace GiNaC {
  59
  60 #ifdef DEBUGFACTOR
  61 #define DCOUT(str) cout << #str << endl
  62 #define DCOUTVAR(var) cout << #var << ": " << var << endl
  63 #define DCOUT2(str,var) cout << #str << ": " << var << endl
  64 #else
  65 #define DCOUT(str)
  66 #define DCOUTVAR(var)
  67 #define DCOUT2(str,var)
  68 #endif
  69
  70 // anonymous namespace to hide all utility functions
  71 namespace {
  72
  73 typedef vector<cl_MI> mvec;
  74 #ifdef DEBUGFACTOR
  75 ostream& operator<<(ostream& o, const vector<int>& v)
  76 {
  77         vector<int>::const_iterator i = v.begin(), end = v.end();
  78         while ( i != end ) {
  79                 o << *i++ << " ";
  80         }
  81         return o;
  82 }
  83 ostream& operator<<(ostream& o, const vector<cl_I>& v)
  84 {
  85         vector<cl_I>::const_iterator i = v.begin(), end = v.end();
  86         while ( i != end ) {
  87                 o << *i << "[" << i-v.begin() << "]" << " ";
  88                 ++i;
  89         }
  90         return o;
  91 }
  92 ostream& operator<<(ostream& o, const vector<cl_MI>& v)
  93 {
  94         vector<cl_MI>::const_iterator i = v.begin(), end = v.end();
  95         while ( i != end ) {
  96                 o << *i << "[" << i-v.begin() << "]" << " ";
  97                 ++i;
  98         }
  99         return o;
 100 }
 101 ostream& operator<<(ostream& o, const vector< vector<cl_MI> >& v)
 102 {
 103         vector< vector<cl_MI> >::const_iterator i = v.begin(), end = v.end();
 104         while ( i != end ) {
 105                 o << i-v.begin() << ": " << *i << endl;
 106                 ++i;
 107         }
 108         return o;
 109 }
 110 #endif
 111
 112 ////////////////////////////////////////////////////////////////////////////////
 113 // modular univariate polynomial code
 114
 115 typedef std::vector<cln::cl_MI> umodpoly;
 116 typedef std::vector<cln::cl_I> upoly;
 117 typedef vector<umodpoly> upvec;
 118
 119 // COPY FROM UPOLY.HPP
 120
 121 // CHANGED size_t -> int !!!
 122 template<typename T> static int degree(const T& p)
 123 {
 124         return p.size() - 1;
 125 }
 126
 127 template<typename T> static typename T::value_type lcoeff(const T& p)
 128 {
 129         return p[p.size() - 1];
 130 }
 131
 132 static bool normalize_in_field(umodpoly& a)
 133 {
 134         if (a.size() == 0)
 135                 return true;
 136         if ( lcoeff(a) == a[0].ring()->one() ) {
 137                 return true;
 138         }
 139
 140         const cln::cl_MI lc_1 = recip(lcoeff(a));
 141         for (std::size_t k = a.size(); k-- != 0; )
 142                 a[k] = a[k]*lc_1;
 143         return false;
 144 }
 145
 146 template<typename T> static void
 147 canonicalize(T& p, const typename T::size_type hint = std::numeric_limits<typename T::size_type>::max())
 148 {
 149         if (p.empty())
 150                 return;
 151
 152         std::size_t i = p.size() - 1;
 153         // Be fast if the polynomial is already canonicalized
 154         if (!zerop(p[i]))
 155                 return;
 156
 157         if (hint < p.size())
 158                 i = hint;
 159
 160         bool is_zero = false;
 161         do {
 162                 if (!zerop(p[i])) {
 163                         ++i;
 164                         break;
 165                 }
 166                 if (i == 0) {
 167                         is_zero = true;
 168                         break;
 169                 }
 170                 --i;
 171         } while (true);
 172
 173         if (is_zero) {
 174                 p.clear();
 175                 return;
 176         }
 177
 178         p.erase(p.begin() + i, p.end());
 179 }
 180
 181 // END COPY FROM UPOLY.HPP
 182
 183 static void expt_pos(umodpoly& a, unsigned int q)
 184 {
 185         if ( a.empty() ) return;
 186         cl_MI zero = a[0].ring()->zero();
 187         int deg = degree(a);
 188         a.resize(degree(a)*q+1, zero);
 189         for ( int i=deg; i>0; --i ) {
 190                 a[i*q] = a[i];
 191                 a[i] = zero;
 192         }
 193 }
 194
 195 template<typename T>
 196 static T operator+(const T& a, const T& b)
 197 {
 198         int sa = a.size();
 199         int sb = b.size();
 200         if ( sa >= sb ) {
 201                 T r(sa);
 202                 int i = 0;
 203                 for ( ; i<sb; ++i ) {
 204                         r[i] = a[i] + b[i];
 205                 }
 206                 for ( ; i<sa; ++i ) {
 207                         r[i] = a[i];
 208                 }
 209                 canonicalize(r);
 210                 return r;
 211         }
 212         else {
 213                 T r(sb);
 214                 int i = 0;
 215                 for ( ; i<sa; ++i ) {
 216                         r[i] = a[i] + b[i];
 217                 }
 218                 for ( ; i<sb; ++i ) {
 219                         r[i] = b[i];
 220                 }
 221                 canonicalize(r);
 222                 return r;
 223         }
 224 }
 225
 226 template<typename T>
 227 static T operator-(const T& a, const T& b)
 228 {
 229         int sa = a.size();
 230         int sb = b.size();
 231         if ( sa >= sb ) {
 232                 T r(sa);
 233                 int i = 0;
 234                 for ( ; i<sb; ++i ) {
 235                         r[i] = a[i] - b[i];
 236                 }
 237                 for ( ; i<sa; ++i ) {
 238                         r[i] = a[i];
 239                 }
 240                 canonicalize(r);
 241                 return r;
 242         }
 243         else {
 244                 T r(sb);
 245                 int i = 0;
 246                 for ( ; i<sa; ++i ) {
 247                         r[i] = a[i] - b[i];
 248                 }
 249                 for ( ; i<sb; ++i ) {
 250                         r[i] = -b[i];
 251                 }
 252                 canonicalize(r);
 253                 return r;
 254         }
 255 }
 256
 257 static upoly operator*(const upoly& a, const upoly& b)
 258 {
 259         upoly c;
 260         if ( a.empty() || b.empty() ) return c;
 261
 262         int n = degree(a) + degree(b);
 263         c.resize(n+1, 0);
 264         for ( int i=0 ; i<=n; ++i ) {
 265                 for ( int j=0 ; j<=i; ++j ) {
 266                         if ( j > degree(a) || (i-j) > degree(b) ) continue;
 267                         c[i] = c[i] + a[j] * b[i-j];
 268                 }
 269         }
 270         canonicalize(c);
 271         return c;
 272 }
 273
 274 static umodpoly operator*(const umodpoly& a, const umodpoly& b)
 275 {
 276         umodpoly c;
 277         if ( a.empty() || b.empty() ) return c;
 278
 279         int n = degree(a) + degree(b);
 280         c.resize(n+1, a[0].ring()->zero());
 281         for ( int i=0 ; i<=n; ++i ) {
 282                 for ( int j=0 ; j<=i; ++j ) {
 283                         if ( j > degree(a) || (i-j) > degree(b) ) continue;
 284                         c[i] = c[i] + a[j] * b[i-j];
 285                 }
 286         }
 287         canonicalize(c);
 288         return c;
 289 }
 290
 291 static upoly operator*(const upoly& a, const cl_I& x)
 292 {
 293         if ( zerop(x) ) {
 294                 upoly r;
 295                 return r;
 296         }
 297         upoly r(a.size());
 298         for ( size_t i=0; i<a.size(); ++i ) {
 299                 r[i] = a[i] * x;
 300         }
 301         return r;
 302 }
 303
 304 static upoly operator/(const upoly& a, const cl_I& x)
 305 {
 306         if ( zerop(x) ) {
 307                 upoly r;
 308                 return r;
 309         }
 310         upoly r(a.size());
 311         for ( size_t i=0; i<a.size(); ++i ) {
 312                 r[i] = exquo(a[i],x);
 313         }
 314         return r;
 315 }
 316
 317 static umodpoly operator*(const umodpoly& a, const cl_MI& x)
 318 {
 319         umodpoly r(a.size());
 320         for ( size_t i=0; i<a.size(); ++i ) {
 321                 r[i] = a[i] * x;
 322         }
 323         canonicalize(r);
 324         return r;
 325 }
 326
 327 static void upoly_from_ex(upoly& up, const ex& e, const ex& x)
 328 {
 329         // assert: e is in Z[x]
 330         int deg = e.degree(x);
 331         up.resize(deg+1);
 332         int ldeg = e.ldegree(x);
 333         for ( ; deg>=ldeg; --deg ) {
 334                 up[deg] = the<cl_I>(ex_to<numeric>(e.coeff(x, deg)).to_cl_N());
 335         }
 336         for ( ; deg>=0; --deg ) {
 337                 up[deg] = 0;
 338         }
 339         canonicalize(up);
 340 }
 341
 342 static void umodpoly_from_upoly(umodpoly& ump, const upoly& e, const cl_modint_ring& R)
 343 {
 344         int deg = degree(e);
 345         ump.resize(deg+1);
 346         for ( ; deg>=0; --deg ) {
 347                 ump[deg] = R->canonhom(e[deg]);
 348         }
 349         canonicalize(ump);
 350 }
 351
 352 static void umodpoly_from_ex(umodpoly& ump, const ex& e, const ex& x, const cl_modint_ring& R)
 353 {
 354         // assert: e is in Z[x]
 355         int deg = e.degree(x);
 356         ump.resize(deg+1);
 357         int ldeg = e.ldegree(x);
 358         for ( ; deg>=ldeg; --deg ) {
 359                 cl_I coeff = the<cl_I>(ex_to<numeric>(e.coeff(x, deg)).to_cl_N());
 360                 ump[deg] = R->canonhom(coeff);
 361         }
 362         for ( ; deg>=0; --deg ) {
 363                 ump[deg] = R->zero();
 364         }
 365         canonicalize(ump);
 366 }
 367
 368 #ifdef DEBUGFACTOR
 369 static void umodpoly_from_ex(umodpoly& ump, const ex& e, const ex& x, const cl_I& modulus)
 370 {
 371         umodpoly_from_ex(ump, e, x, find_modint_ring(modulus));
 372 }
 373 #endif
 374
 375 static ex upoly_to_ex(const upoly& a, const ex& x)
 376 {
 377         if ( a.empty() ) return 0;
 378         ex e;
 379         for ( int i=degree(a); i>=0; --i ) {
 380                 e += numeric(a[i]) * pow(x, i);
 381         }
 382         return e;
 383 }
 384
 385 static ex umodpoly_to_ex(const umodpoly& a, const ex& x)
 386 {
 387         if ( a.empty() ) return 0;
 388         cl_modint_ring R = a[0].ring();
 389         cl_I mod = R->modulus;
 390         cl_I halfmod = (mod-1) >> 1;
 391         ex e;
 392         for ( int i=degree(a); i>=0; --i ) {
 393                 cl_I n = R->retract(a[i]);
 394                 if ( n > halfmod ) {
 395                         e += numeric(n-mod) * pow(x, i);
 396                 } else {
 397                         e += numeric(n) * pow(x, i);
 398                 }
 399         }
 400         return e;
 401 }
 402
 403 static upoly umodpoly_to_upoly(const umodpoly& a)
 404 {
 405         upoly e(a.size());
 406         if ( a.empty() ) return e;
 407         cl_modint_ring R = a[0].ring();
 408         cl_I mod = R->modulus;
 409         cl_I halfmod = (mod-1) >> 1;
 410         for ( int i=degree(a); i>=0; --i ) {
 411                 cl_I n = R->retract(a[i]);
 412                 if ( n > halfmod ) {
 413                         e[i] = n-mod;
 414                 } else {
 415                         e[i] = n;
 416                 }
 417         }
 418         return e;
 419 }
 420
 421 static umodpoly umodpoly_to_umodpoly(const umodpoly& a, const cl_modint_ring& R, unsigned int m)
 422 {
 423         umodpoly e;
 424         if ( a.empty() ) return e;
 425         cl_modint_ring oldR = a[0].ring();
 426         size_t sa = a.size();
 427         e.resize(sa+m, R->zero());
 428         for ( size_t i=0; i<sa; ++i ) {
 429                 e[i+m] = R->canonhom(oldR->retract(a[i]));
 430         }
 431         canonicalize(e);
 432         return e;
 433 }
 434
 435 /** Divides all coefficients of the polynomial a by the integer x.
 436  *  All coefficients are supposed to be divisible by x. If they are not, the
 437  *  the<cl_I> cast will raise an exception.
 438  *
 439  *  @param[in,out] a  polynomial of which the coefficients will be reduced by x
 440  *  @param[in]     x  integer that divides the coefficients
 441  */
 442 static void reduce_coeff(umodpoly& a, const cl_I& x)
 443 {
 444         if ( a.empty() ) return;
 445
 446         cl_modint_ring R = a[0].ring();
 447         umodpoly::iterator i = a.begin(), end = a.end();
 448         for ( ; i!=end; ++i ) {
 449                 // cln cannot perform this division in the modular field
 450                 cl_I c = R->retract(*i);
 451                 *i = cl_MI(R, the<cl_I>(c / x));
 452         }
 453 }
 454
 455 /** Calculates remainder of a/b.
 456  *  Assertion: a and b not empty.
 457  *
 458  *  @param[in]  a  polynomial dividend
 459  *  @param[in]  b  polynomial divisor
 460  *  @param[out] r  polynomial remainder
 461  */
 462 static void rem(const umodpoly& a, const umodpoly& b, umodpoly& r)
 463 {
 464         int k, n;
 465         n = degree(b);
 466         k = degree(a) - n;
 467         r = a;
 468         if ( k < 0 ) return;
 469
 470         do {
 471                 cl_MI qk = div(r[n+k], b[n]);
 472                 if ( !zerop(qk) ) {
 473                         for ( int i=0; i<n; ++i ) {
 474                                 unsigned int j = n + k - 1 - i;
 475                                 r[j] = r[j] - qk * b[j-k];
 476                         }
 477                 }
 478         } while ( k-- );
 479
 480         fill(r.begin()+n, r.end(), a[0].ring()->zero());
 481         canonicalize(r);
 482 }
 483
 484 /** Calculates quotient of a/b.
 485  *  Assertion: a and b not empty.
 486  *
 487  *  @param[in]  a  polynomial dividend
 488  *  @param[in]  b  polynomial divisor
 489  *  @param[out] q  polynomial quotient
 490  */
 491 static void div(const umodpoly& a, const umodpoly& b, umodpoly& q)
 492 {
 493         int k, n;
 494         n = degree(b);
 495         k = degree(a) - n;
 496         q.clear();
 497         if ( k < 0 ) return;
 498
 499         umodpoly r = a;
 500         q.resize(k+1, a[0].ring()->zero());
 501         do {
 502                 cl_MI qk = div(r[n+k], b[n]);
 503                 if ( !zerop(qk) ) {
 504                         q[k] = qk;
 505                         for ( int i=0; i<n; ++i ) {
 506                                 unsigned int j = n + k - 1 - i;
 507                                 r[j] = r[j] - qk * b[j-k];
 508                         }
 509                 }
 510         } while ( k-- );
 511
 512         canonicalize(q);
 513 }
 514
 515 /** Calculates quotient and remainder of a/b.
 516  *  Assertion: a and b not empty.
 517  *
 518  *  @param[in]  a  polynomial dividend
 519  *  @param[in]  b  polynomial divisor
 520  *  @param[out] r  polynomial remainder
 521  *  @param[out] q  polynomial quotient
 522  */
 523 static void remdiv(const umodpoly& a, const umodpoly& b, umodpoly& r, umodpoly& q)
 524 {
 525         int k, n;
 526         n = degree(b);
 527         k = degree(a) - n;
 528         q.clear();
 529         r = a;
 530         if ( k < 0 ) return;
 531
 532         q.resize(k+1, a[0].ring()->zero());
 533         do {
 534                 cl_MI qk = div(r[n+k], b[n]);
 535                 if ( !zerop(qk) ) {
 536                         q[k] = qk;
 537                         for ( int i=0; i<n; ++i ) {
 538                                 unsigned int j = n + k - 1 - i;
 539                                 r[j] = r[j] - qk * b[j-k];
 540                         }
 541                 }
 542         } while ( k-- );
 543
 544         fill(r.begin()+n, r.end(), a[0].ring()->zero());
 545         canonicalize(r);
 546         canonicalize(q);
 547 }
 548
 549 /** Calculates the GCD of polynomial a and b.
 550  *
 551  *  @param[in]  a  polynomial
 552  *  @param[in]  b  polynomial
 553  *  @param[out] c  GCD
 554  */
 555 static void gcd(const umodpoly& a, const umodpoly& b, umodpoly& c)
 556 {
 557         if ( degree(a) < degree(b) ) return gcd(b, a, c);
 558
 559         c = a;
 560         normalize_in_field(c);
 561         umodpoly d = b;
 562         normalize_in_field(d);
 563         umodpoly r;
 564         while ( !d.empty() ) {
 565                 rem(c, d, r);
 566                 c = d;
 567                 d = r;
 568         }
 569         normalize_in_field(c);
 570 }
 571
 572 /** Calculates the derivative of the polynomial a.
 573  *
 574  *  @param[in]  a  polynomial of which to take the derivative
 575  *  @param[out] d  result/derivative
 576  */
 577 static void deriv(const umodpoly& a, umodpoly& d)
 578 {
 579         d.clear();
 580         if ( a.size() <= 1 ) return;
 581
 582         d.insert(d.begin(), a.begin()+1, a.end());
 583         int max = d.size();
 584         for ( int i=1; i<max; ++i ) {
 585                 d[i] = d[i] * (i+1);
 586         }
 587         canonicalize(d);
 588 }
 589
 590 static bool unequal_one(const umodpoly& a)
 591 {
 592         if ( a.empty() ) return true;
 593         return ( a.size() != 1 || a[0] != a[0].ring()->one() );
 594 }
 595
 596 static bool equal_one(const umodpoly& a)
 597 {
 598         return ( a.size() == 1 && a[0] == a[0].ring()->one() );
 599 }
 600
 601 /** Returns true if polynomial a is square free.
 602  *
 603  *  @param[in] a  polynomial to check
 604  *  @return       true if polynomial is square free, false otherwise
 605  */
 606 static bool squarefree(const umodpoly& a)
 607 {
 608         umodpoly b;
 609         deriv(a, b);
 610         if ( b.empty() ) {
 611                 return false;
 612         }
 613         umodpoly c;
 614         gcd(a, b, c);
 615         return equal_one(c);
 616 }
 617
 618 // END modular univariate polynomial code
 619 ////////////////////////////////////////////////////////////////////////////////
 620
 621 ////////////////////////////////////////////////////////////////////////////////
 622 // modular matrix
 623
 624 class modular_matrix
 625 {
 626         friend ostream& operator<<(ostream& o, const modular_matrix& m);
 627 public:
 628         modular_matrix(size_t r_, size_t c_, const cl_MI& init) : r(r_), c(c_)
 629         {
 630                 m.resize(c*r, init);
 631         }
 632         size_t rowsize() const { return r; }
 633         size_t colsize() const { return c; }
 634         cl_MI& operator()(size_t row, size_t col) { return m[row*c + col]; }
 635         cl_MI operator()(size_t row, size_t col) const { return m[row*c + col]; }
 636         void mul_col(size_t col, const cl_MI x)
 637         {
 638                 mvec::iterator i = m.begin() + col;
 639                 for ( size_t rc=0; rc<r; ++rc ) {
 640                         *i = *i * x;
 641                         i += c;
 642                 }
 643         }
 644         void sub_col(size_t col1, size_t col2, const cl_MI fac)
 645         {
 646                 mvec::iterator i1 = m.begin() + col1;
 647                 mvec::iterator i2 = m.begin() + col2;
 648                 for ( size_t rc=0; rc<r; ++rc ) {
 649                         *i1 = *i1 - *i2 * fac;
 650                         i1 += c;
 651                         i2 += c;
 652                 }
 653         }
 654         void switch_col(size_t col1, size_t col2)
 655         {
 656                 cl_MI buf;
 657                 mvec::iterator i1 = m.begin() + col1;
 658                 mvec::iterator i2 = m.begin() + col2;
 659                 for ( size_t rc=0; rc<r; ++rc ) {
 660                         buf = *i1; *i1 = *i2; *i2 = buf;
 661                         i1 += c;
 662                         i2 += c;
 663                 }
 664         }
 665         void mul_row(size_t row, const cl_MI x)
 666         {
 667                 vector<cl_MI>::iterator i = m.begin() + row*c;
 668                 for ( size_t cc=0; cc<c; ++cc ) {
 669                         *i = *i * x;
 670                         ++i;
 671                 }
 672         }
 673         void sub_row(size_t row1, size_t row2, const cl_MI fac)
 674         {
 675                 vector<cl_MI>::iterator i1 = m.begin() + row1*c;
 676                 vector<cl_MI>::iterator i2 = m.begin() + row2*c;
 677                 for ( size_t cc=0; cc<c; ++cc ) {
 678                         *i1 = *i1 - *i2 * fac;
 679                         ++i1;
 680                         ++i2;
 681                 }
 682         }
 683         void switch_row(size_t row1, size_t row2)
 684         {
 685                 cl_MI buf;
 686                 vector<cl_MI>::iterator i1 = m.begin() + row1*c;
 687                 vector<cl_MI>::iterator i2 = m.begin() + row2*c;
 688                 for ( size_t cc=0; cc<c; ++cc ) {
 689                         buf = *i1; *i1 = *i2; *i2 = buf;
 690                         ++i1;
 691                         ++i2;
 692                 }
 693         }
 694         bool is_col_zero(size_t col) const
 695         {
 696                 mvec::const_iterator i = m.begin() + col;
 697                 for ( size_t rr=0; rr<r; ++rr ) {
 698                         if ( !zerop(*i) ) {
 699                                 return false;
 700                         }
 701                         i += c;
 702                 }
 703                 return true;
 704         }
 705         bool is_row_zero(size_t row) const
 706         {
 707                 mvec::const_iterator i = m.begin() + row*c;
 708                 for ( size_t cc=0; cc<c; ++cc ) {
 709                         if ( !zerop(*i) ) {
 710                                 return false;
 711                         }
 712                         ++i;
 713                 }
 714                 return true;
 715         }
 716         void set_row(size_t row, const vector<cl_MI>& newrow)
 717         {
 718                 mvec::iterator i1 = m.begin() + row*c;
 719                 mvec::const_iterator i2 = newrow.begin(), end = newrow.end();
 720                 for ( ; i2 != end; ++i1, ++i2 ) {
 721                         *i1 = *i2;
 722                 }
 723         }
 724         mvec::const_iterator row_begin(size_t row) const { return m.begin()+row*c; }
 725         mvec::const_iterator row_end(size_t row) const { return m.begin()+row*c+r; }
 726 private:
 727         size_t r, c;
 728         mvec m;
 729 };
 730
 731 #ifdef DEBUGFACTOR
 732 modular_matrix operator*(const modular_matrix& m1, const modular_matrix& m2)
 733 {
 734         const unsigned int r = m1.rowsize();
 735         const unsigned int c = m2.colsize();
 736         modular_matrix o(r,c,m1(0,0));
 737
 738         for ( size_t i=0; i<r; ++i ) {
 739                 for ( size_t j=0; j<c; ++j ) {
 740                         cl_MI buf;
 741                         buf = m1(i,0) * m2(0,j);
 742                         for ( size_t k=1; k<c; ++k ) {
 743                                 buf = buf + m1(i,k)*m2(k,j);
 744                         }
 745                         o(i,j) = buf;
 746                 }
 747         }
 748         return o;
 749 }
 750
 751 ostream& operator<<(ostream& o, const modular_matrix& m)
 752 {
 753         cl_modint_ring R = m(0,0).ring();
 754         o << "{";
 755         for ( size_t i=0; i<m.rowsize(); ++i ) {
 756                 o << "{";
 757                 for ( size_t j=0; j<m.colsize()-1; ++j ) {
 758                         o << R->retract(m(i,j)) << ",";
 759                 }
 760                 o << R->retract(m(i,m.colsize()-1)) << "}";
 761                 if ( i != m.rowsize()-1 ) {
 762                         o << ",";
 763                 }
 764         }
 765         o << "}";
 766         return o;
 767 }
 768 #endif // def DEBUGFACTOR
 769
 770 // END modular matrix
 771 ////////////////////////////////////////////////////////////////////////////////
 772
 773 static void q_matrix(const umodpoly& a_, modular_matrix& Q)
 774 {
 775         umodpoly a = a_;
 776         normalize_in_field(a);
 777
 778         int n = degree(a);
 779         unsigned int q = cl_I_to_uint(a[0].ring()->modulus);
 780         umodpoly r(n, a[0].ring()->zero());
 781         r[0] = a[0].ring()->one();
 782         Q.set_row(0, r);
 783         unsigned int max = (n-1) * q;
 784         for ( size_t m=1; m<=max; ++m ) {
 785                 cl_MI rn_1 = r.back();
 786                 for ( size_t i=n-1; i>0; --i ) {
 787                         r[i] = r[i-1] - (rn_1 * a[i]);
 788                 }
 789                 r[0] = -rn_1 * a[0];
 790                 if ( (m % q) == 0 ) {
 791                         Q.set_row(m/q, r);
 792                 }
 793         }
 794 }
 795
 796 static void nullspace(modular_matrix& M, vector<mvec>& basis)
 797 {
 798         const size_t n = M.rowsize();
 799         const cl_MI one = M(0,0).ring()->one();
 800         for ( size_t i=0; i<n; ++i ) {
 801                 M(i,i) = M(i,i) - one;
 802         }
 803         for ( size_t r=0; r<n; ++r ) {
 804                 size_t cc = 0;
 805                 for ( ; cc<n; ++cc ) {
 806                         if ( !zerop(M(r,cc)) ) {
 807                                 if ( cc < r ) {
 808                                         if ( !zerop(M(cc,cc)) ) {
 809                                                 continue;
 810                                         }
 811                                         M.switch_col(cc, r);
 812                                 }
 813                                 else if ( cc > r ) {
 814                                         M.switch_col(cc, r);
 815                                 }
 816                                 break;
 817                         }
 818                 }
 819                 if ( cc < n ) {
 820                         M.mul_col(r, recip(M(r,r)));
 821                         for ( cc=0; cc<n; ++cc ) {
 822                                 if ( cc != r ) {
 823                                         M.sub_col(cc, r, M(r,cc));
 824                                 }
 825                         }
 826                 }
 827         }
 828
 829         for ( size_t i=0; i<n; ++i ) {
 830                 M(i,i) = M(i,i) - one;
 831         }
 832         for ( size_t i=0; i<n; ++i ) {
 833                 if ( !M.is_row_zero(i) ) {
 834                         mvec nu(M.row_begin(i), M.row_end(i));
 835                         basis.push_back(nu);
 836                 }
 837         }
 838 }
 839
 840 static void berlekamp(const umodpoly& a, upvec& upv)
 841 {
 842         cl_modint_ring R = a[0].ring();
 843         umodpoly one(1, R->one());
 844
 845         modular_matrix Q(degree(a), degree(a), R->zero());
 846         q_matrix(a, Q);
 847         vector<mvec> nu;
 848         nullspace(Q, nu);
 849
 850         const unsigned int k = nu.size();
 851         if ( k == 1 ) {
 852                 return;
 853         }
 854
 855         list<umodpoly> factors;
 856         factors.push_back(a);
 857         unsigned int size = 1;
 858         unsigned int r = 1;
 859         unsigned int q = cl_I_to_uint(R->modulus);
 860
 861         list<umodpoly>::iterator u = factors.begin();
 862
 863         while ( true ) {
 864                 for ( unsigned int s=0; s<q; ++s ) {
 865                         umodpoly nur = nu[r];
 866                         nur[0] = nur[0] - cl_MI(R, s);
 867                         canonicalize(nur);
 868                         umodpoly g;
 869                         gcd(nur, *u, g);
 870                         if ( unequal_one(g) && g != *u ) {
 871                                 umodpoly uo;
 872                                 div(*u, g, uo);
 873                                 if ( equal_one(uo) ) {
 874                                         throw logic_error("berlekamp: unexpected divisor.");
 875                                 }
 876                                 else {
 877                                         *u = uo;
 878                                 }
 879                                 factors.push_back(g);
 880                                 size = 0;
 881                                 list<umodpoly>::const_iterator i = factors.begin(), end = factors.end();
 882                                 while ( i != end ) {
 883                                         if ( degree(*i) ) ++size;
 884                                         ++i;
 885                                 }
 886                                 if ( size == k ) {
 887                                         list<umodpoly>::const_iterator i = factors.begin(), end = factors.end();
 888                                         while ( i != end ) {
 889                                                 upv.push_back(*i++);
 890                                         }
 891                                         return;
 892                                 }
 893                         }
 894                 }
 895                 if ( ++r == k ) {
 896                         r = 1;
 897                         ++u;
 898                 }
 899         }
 900 }
 901
 902 static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap)
 903 {
 904         size_t newdeg = degree(a)/prime;
 905         ap.resize(newdeg+1);
 906         ap[0] = a[0];
 907         for ( size_t i=1; i<=newdeg; ++i ) {
 908                 ap[i] = a[i*prime];
 909         }
 910 }
 911
 912 static void modsqrfree(const umodpoly& a, upvec& factors, vector<int>& mult)
 913 {
 914         const unsigned int prime = cl_I_to_uint(a[0].ring()->modulus);
 915         int i = 1;
 916         umodpoly b;
 917         deriv(a, b);
 918         if ( b.size() ) {
 919                 umodpoly c;
 920                 gcd(a, b, c);
 921                 umodpoly w;
 922                 div(a, c, w);
 923                 while ( unequal_one(w) ) {
 924                         umodpoly y;
 925                         gcd(w, c, y);
 926                         umodpoly z;
 927                         div(w, y, z);
 928                         factors.push_back(z);
 929                         mult.push_back(i);
 930                         ++i;
 931                         w = y;
 932                         umodpoly buf;
 933                         div(c, y, buf);
 934                         c = buf;
 935                 }
 936                 if ( unequal_one(c) ) {
 937                         umodpoly cp;
 938                         expt_1_over_p(c, prime, cp);
 939                         size_t previ = mult.size();
 940                         modsqrfree(cp, factors, mult);
 941                         for ( size_t i=previ; i<mult.size(); ++i ) {
 942                                 mult[i] *= prime;
 943                         }
 944                 }
 945         }
 946         else {
 947                 umodpoly ap;
 948                 expt_1_over_p(a, prime, ap);
 949                 size_t previ = mult.size();
 950                 modsqrfree(ap, factors, mult);
 951                 for ( size_t i=previ; i<mult.size(); ++i ) {
 952                         mult[i] *= prime;
 953                 }
 954         }
 955 }
 956
 957 static void distinct_degree_factor(const umodpoly& a_, vector<int>& degrees, upvec& ddfactors)
 958 {
 959         umodpoly a = a_;
 960
 961         cl_modint_ring R = a[0].ring();
 962         int q = cl_I_to_int(R->modulus);
 963         int nhalf = degree(a)/2;
 964
 965         int i = 1;
 966         umodpoly w(2);
 967         w[0] = R->zero();
 968         w[1] = R->one();
 969         umodpoly x = w;
 970
 971         while ( i <= nhalf ) {
 972                 expt_pos(w, q);
 973                 umodpoly buf;
 974                 rem(w, a, buf);
 975                 w = buf;
 976                 umodpoly wx = w - x;
 977                 gcd(a, wx, buf);
 978                 if ( unequal_one(buf) ) {
 979                         degrees.push_back(i);
 980                         ddfactors.push_back(buf);
 981                 }
 982                 if ( unequal_one(buf) ) {
 983                         umodpoly buf2;
 984                         div(a, buf, buf2);
 985                         a = buf2;
 986                         nhalf = degree(a)/2;
 987                         rem(w, a, buf);
 988                         w = buf;
 989                 }
 990                 ++i;
 991         }
 992         if ( unequal_one(a) ) {
 993                 degrees.push_back(degree(a));
 994                 ddfactors.push_back(a);
 995         }
 996 }
 997
 998 static void same_degree_factor(const umodpoly& a, upvec& upv)
 999 {
1000         cl_modint_ring R = a[0].ring();
1001
1002         vector<int> degrees;
1003         upvec ddfactors;
1004         distinct_degree_factor(a, degrees, ddfactors);
1005
1006         for ( size_t i=0; i<degrees.size(); ++i ) {
1007                 if ( degrees[i] == degree(ddfactors[i]) ) {
1008                         upv.push_back(ddfactors[i]);
1009                 }
1010                 else {
1011                         berlekamp(ddfactors[i], upv);
1012                 }
1013         }
1014 }
1015
1016 #define USE_SAME_DEGREE_FACTOR
1017
1018 static void factor_modular(const umodpoly& p, upvec& upv)
1019 {
1020 #ifdef USE_SAME_DEGREE_FACTOR
1021         same_degree_factor(p, upv);
1022 #else
1023         berlekamp(p, upv);
1024 #endif
1025 }
1026
1027 /** Calculates polynomials s and t such that a*s+b*t==1.
1028  *  Assertion: a and b are relatively prime and not zero.
1029  *
1030  *  @param[in]  a  polynomial
1031  *  @param[in]  b  polynomial
1032  *  @param[out] s  polynomial
1033  *  @param[out] t  polynomial
1034  */
1035 static void exteuclid(const umodpoly& a, const umodpoly& b, umodpoly& s, umodpoly& t)
1036 {
1037         if ( degree(a) < degree(b) ) {
1038                 exteuclid(b, a, t, s);
1039                 return;
1040         }
1041
1042         umodpoly one(1, a[0].ring()->one());
1043         umodpoly c = a; normalize_in_field(c);
1044         umodpoly d = b; normalize_in_field(d);
1045         s = one;
1046         t.clear();
1047         umodpoly d1;
1048         umodpoly d2 = one;
1049         umodpoly q;
1050         while ( true ) {
1051                 div(c, d, q);
1052                 umodpoly r = c - q * d;
1053                 umodpoly r1 = s - q * d1;
1054                 umodpoly r2 = t - q * d2;
1055                 c = d;
1056                 s = d1;
1057                 t = d2;
1058                 if ( r.empty() ) break;
1059                 d = r;
1060                 d1 = r1;
1061                 d2 = r2;
1062         }
1063         cl_MI fac = recip(lcoeff(a) * lcoeff(c));
1064         umodpoly::iterator i = s.begin(), end = s.end();
1065         for ( ; i!=end; ++i ) {
1066                 *i = *i * fac;
1067         }
1068         canonicalize(s);
1069         fac = recip(lcoeff(b) * lcoeff(c));
1070         i = t.begin(), end = t.end();
1071         for ( ; i!=end; ++i ) {
1072                 *i = *i * fac;
1073         }
1074         canonicalize(t);
1075 }
1076
1077 static upoly replace_lc(const upoly& poly, const cl_I& lc)
1078 {
1079         if ( poly.empty() ) return poly;
1080         upoly r = poly;
1081         r.back() = lc;
1082         return r;
1083 }
1084
1085 static inline cl_I calc_bound(const ex& a, const ex& x, int maxdeg)
1086 {
1087         cl_I maxcoeff = 0;
1088         cl_R coeff = 0;
1089         for ( int i=a.degree(x); i>=a.ldegree(x); --i ) {
1090                 cl_I aa = abs(the<cl_I>(ex_to<numeric>(a.coeff(x, i)).to_cl_N()));
1091                 if ( aa > maxcoeff ) maxcoeff = aa;
1092                 coeff = coeff + square(aa);
1093         }
1094         cl_I coeffnorm = ceiling1(the<cl_R>(cln::sqrt(coeff)));
1095         cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg));
1096         return ( B > maxcoeff ) ? B : maxcoeff;
1097 }
1098
1099 static inline cl_I calc_bound(const upoly& a, int maxdeg)
1100 {
1101         cl_I maxcoeff = 0;
1102         cl_R coeff = 0;
1103         for ( int i=degree(a); i>=0; --i ) {
1104                 cl_I aa = abs(a[i]);
1105                 if ( aa > maxcoeff ) maxcoeff = aa;
1106                 coeff = coeff + square(aa);
1107         }
1108         cl_I coeffnorm = ceiling1(the<cl_R>(cln::sqrt(coeff)));
1109         cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg));
1110         return ( B > maxcoeff ) ? B : maxcoeff;
1111 }
1112
1113 static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_, const umodpoly& w1_, upoly& u, upoly& w)
1114 {
1115         upoly a = a_;
1116         const cl_modint_ring& R = u1_[0].ring();
1117
1118         // calc bound B
1119         int maxdeg = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_);
1120         cl_I maxmodulus = 2*calc_bound(a, maxdeg);
1121
1122         // step 1
1123         cl_I alpha = lcoeff(a);
1124         a = a * alpha;
1125         umodpoly nu1 = u1_;
1126         normalize_in_field(nu1);
1127         umodpoly nw1 = w1_;
1128         normalize_in_field(nw1);
1129         upoly phi;
1130         phi = umodpoly_to_upoly(nu1) * alpha;
1131         umodpoly u1;
1132         umodpoly_from_upoly(u1, phi, R);
1133         phi = umodpoly_to_upoly(nw1) * alpha;
1134         umodpoly w1;
1135         umodpoly_from_upoly(w1, phi, R);
1136
1137         // step 2
1138         umodpoly s;
1139         umodpoly t;
1140         exteuclid(u1, w1, s, t);
1141
1142         // step 3
1143         u = replace_lc(umodpoly_to_upoly(u1), alpha);
1144         w = replace_lc(umodpoly_to_upoly(w1), alpha);
1145         upoly e = a - u * w;
1146         cl_I modulus = p;
1147
1148         // step 4
1149         while ( !e.empty() && modulus < maxmodulus ) {
1150                 upoly c = e / modulus;
1151                 phi = umodpoly_to_upoly(s) * c;
1152                 umodpoly sigmatilde;
1153                 umodpoly_from_upoly(sigmatilde, phi, R);
1154                 phi = umodpoly_to_upoly(t) * c;
1155                 umodpoly tautilde;
1156                 umodpoly_from_upoly(tautilde, phi, R);
1157                 umodpoly r, q;
1158                 remdiv(sigmatilde, w1, r, q);
1159                 umodpoly sigma = r;
1160                 phi = umodpoly_to_upoly(tautilde) + umodpoly_to_upoly(q) * umodpoly_to_upoly(u1);
1161                 umodpoly tau;
1162                 umodpoly_from_upoly(tau, phi, R);
1163                 u = u + umodpoly_to_upoly(tau) * modulus;
1164                 w = w + umodpoly_to_upoly(sigma) * modulus;
1165                 e = a - u * w;
1166                 modulus = modulus * p;
1167         }
1168
1169         // step 5
1170         if ( e.empty() ) {
1171                 cl_I g = u[0];
1172                 for ( size_t i=1; i<u.size(); ++i ) {
1173                         g = gcd(g, u[i]);
1174                         if ( g == 1 ) break;
1175                 }
1176                 if ( g != 1 ) {
1177                         u = u / g;
1178                         w = w * g;
1179                 }
1180                 if ( alpha != 1 ) {
1181                         w = w / alpha;
1182                 }
1183         }
1184         else {
1185                 u.clear();
1186         }
1187 }
1188
1189 static unsigned int next_prime(unsigned int p)
1190 {
1191         static vector<unsigned int> primes;
1192         if ( primes.size() == 0 ) {
1193                 primes.push_back(3); primes.push_back(5); primes.push_back(7);
1194         }
1195         vector<unsigned int>::const_iterator it = primes.begin();
1196         if ( p >= primes.back() ) {
1197                 unsigned int candidate = primes.back() + 2;
1198                 while ( true ) {
1199                         size_t n = primes.size()/2;
1200                         for ( size_t i=0; i<n; ++i ) {
1201                                 if ( candidate % primes[i] ) continue;
1202                                 candidate += 2;
1203                                 i=-1;
1204                         }
1205                         primes.push_back(candidate);
1206                         if ( candidate > p ) break;
1207                 }
1208                 return candidate;
1209         }
1210         vector<unsigned int>::const_iterator end = primes.end();
1211         for ( ; it!=end; ++it ) {
1212                 if ( *it > p ) {
1213                         return *it;
1214                 }
1215         }
1216         throw logic_error("next_prime: should not reach this point!");
1217 }
1218
1219 class factor_partition
1220 {
1221 public:
1222         factor_partition(const upvec& factors_) : factors(factors_)
1223         {
1224                 n = factors.size();
1225                 k.resize(n, 0);
1226                 k[0] = 1;
1227                 cache.resize(n-1);
1228                 one.resize(1, factors.front()[0].ring()->one());
1229                 len = 1;
1230                 last = 0;
1231                 split();
1232         }
1233         int operator[](size_t i) const { return k[i]; }
1234         size_t size() const { return n; }
1235         size_t size_left() const { return n-len; }
1236         size_t size_right() const { return len; }
1237 #ifdef DEBUGFACTOR
1238         void get() const { DCOUTVAR(k); }
1239 #endif
1240         bool next()
1241         {
1242                 if ( last == n-1 ) {
1243                         int rem = len - 1;
1244                         int p = last - 1;
1245                         while ( rem ) {
1246                                 if ( k[p] ) {
1247                                         --rem;
1248                                         --p;
1249                                         continue;
1250                                 }
1251                                 last = p - 1;
1252                                 while ( k[last] == 0 ) { --last; }
1253                                 if ( last == 0 && n == 2*len ) return false;
1254                                 k[last++] = 0;
1255                                 for ( size_t i=0; i<=len-rem; ++i ) {
1256                                         k[last] = 1;
1257                                         ++last;
1258                                 }
1259                                 fill(k.begin()+last, k.end(), 0);
1260                                 --last;
1261                                 split();
1262                                 return true;
1263                         }
1264                         last = len;
1265                         ++len;
1266                         if ( len > n/2 ) return false;
1267                         fill(k.begin(), k.begin()+len, 1);
1268                         fill(k.begin()+len+1, k.end(), 0);
1269                 }
1270                 else {
1271                         k[last++] = 0;
1272                         k[last] = 1;
1273                 }
1274                 split();
1275                 return true;
1276         }
1277         umodpoly& left() { return lr[0]; }
1278         umodpoly& right() { return lr[1]; }
1279 private:
1280         void split_cached()
1281         {
1282                 size_t i = 0;
1283                 do {
1284                         size_t pos = i;
1285                         int group = k[i++];
1286                         size_t d = 0;
1287                         while ( i < n && k[i] == group ) { ++d; ++i; }
1288                         if ( d ) {
1289                                 if ( cache[pos].size() >= d ) {
1290                                         lr[group] = lr[group] * cache[pos][d-1];
1291                                 }
1292                                 else {
1293                                         if ( cache[pos].size() == 0 ) {
1294                                                 cache[pos].push_back(factors[pos] * factors[pos+1]);
1295                                         }
1296                                         size_t j = pos + cache[pos].size() + 1;
1297                                         d -= cache[pos].size();
1298                                         while ( d ) {
1299                                                 umodpoly buf = cache[pos].back() * factors[j];
1300                                                 cache[pos].push_back(buf);
1301                                                 --d;
1302                                                 ++j;
1303                                         }
1304                                         lr[group] = lr[group] * cache[pos].back();
1305                                 }
1306                         }
1307                         else {
1308                                 lr[group] = lr[group] * factors[pos];
1309                         }
1310                 } while ( i < n );
1311         }
1312         void split()
1313         {
1314                 lr[0] = one;
1315                 lr[1] = one;
1316                 if ( n > 6 ) {
1317                         split_cached();
1318                 }
1319                 else {
1320                         for ( size_t i=0; i<n; ++i ) {
1321                                 lr[k[i]] = lr[k[i]] * factors[i];
1322                         }
1323                 }
1324         }
1325 private:
1326         umodpoly lr[2];
1327         vector< vector<umodpoly> > cache;
1328         upvec factors;
1329         umodpoly one;
1330         size_t n;
1331         size_t len;
1332         size_t last;
1333         vector<int> k;
1334 };
1335
1336 struct ModFactors
1337 {
1338         upoly poly;
1339         upvec factors;
1340 };
1341
1342 static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
1343 {
1344         ex unit, cont, prim_ex;
1345         poly.unitcontprim(x, unit, cont, prim_ex);
1346         upoly prim;
1347         upoly_from_ex(prim, prim_ex, x);
1348
1349         // determine proper prime and minimize number of modular factors
1350         prime = 3;
1351         unsigned int lastp = prime;
1352         cl_modint_ring R;
1353         unsigned int trials = 0;
1354         unsigned int minfactors = 0;
1355         cl_I lc = lcoeff(prim) * the<cl_I>(ex_to<numeric>(cont).to_cl_N());
1356         upvec factors;
1357         while ( trials < 2 ) {
1358                 umodpoly modpoly;
1359                 while ( true ) {
1360                         prime = next_prime(prime);
1361                         if ( !zerop(rem(lc, prime)) ) {
1362                                 R = find_modint_ring(prime);
1363                                 umodpoly_from_upoly(modpoly, prim, R);
1364                                 if ( squarefree(modpoly) ) break;
1365                         }
1366                 }
1367
1368                 // do modular factorization
1369                 upvec trialfactors;
1370                 factor_modular(modpoly, trialfactors);
1371                 if ( trialfactors.size() <= 1 ) {
1372                         // irreducible for sure
1373                         return poly;
1374                 }
1375
1376                 if ( minfactors == 0 || trialfactors.size() < minfactors ) {
1377                         factors = trialfactors;
1378                         minfactors = trialfactors.size();
1379                         lastp = prime;
1380                         trials = 1;
1381                 }
1382                 else {
1383                         ++trials;
1384                 }
1385         }
1386         prime = lastp;
1387         R = find_modint_ring(prime);
1388
1389         // lift all factor combinations
1390         stack<ModFactors> tocheck;
1391         ModFactors mf;
1392         mf.poly = prim;
1393         mf.factors = factors;
1394         tocheck.push(mf);
1395         upoly f1, f2;
1396         ex result = 1;
1397         while ( tocheck.size() ) {
1398                 const size_t n = tocheck.top().factors.size();
1399                 factor_partition part(tocheck.top().factors);
1400                 while ( true ) {
1401                         hensel_univar(tocheck.top().poly, prime, part.left(), part.right(), f1, f2);
1402                         if ( !f1.empty() ) {
1403                                 if ( part.size_left() == 1 ) {
1404                                         if ( part.size_right() == 1 ) {
1405                                                 result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
1406                                                 tocheck.pop();
1407                                                 break;
1408                                         }
1409                                         result *= upoly_to_ex(f1, x);
1410                                         tocheck.top().poly = f2;
1411                                         for ( size_t i=0; i<n; ++i ) {
1412                                                 if ( part[i] == 0 ) {
1413                                                         tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1414                                                         break;
1415                                                 }
1416                                         }
1417                                         break;
1418                                 }
1419                                 else if ( part.size_right() == 1 ) {
1420                                         if ( part.size_left() == 1 ) {
1421                                                 result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
1422                                                 tocheck.pop();
1423                                                 break;
1424                                         }
1425                                         result *= upoly_to_ex(f2, x);
1426                                         tocheck.top().poly = f1;
1427                                         for ( size_t i=0; i<n; ++i ) {
1428                                                 if ( part[i] == 1 ) {
1429                                                         tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1430                                                         break;
1431                                                 }
1432                                         }
1433                                         break;
1434                                 }
1435                                 else {
1436                                         upvec newfactors1(part.size_left()), newfactors2(part.size_right());
1437                                         upvec::iterator i1 = newfactors1.begin(), i2 = newfactors2.begin();
1438                                         for ( size_t i=0; i<n; ++i ) {
1439                                                 if ( part[i] ) {
1440                                                         *i2++ = tocheck.top().factors[i];
1441                                                 }
1442                                                 else {
1443                                                         *i1++ = tocheck.top().factors[i];
1444                                                 }
1445                                         }
1446                                         tocheck.top().factors = newfactors1;
1447                                         tocheck.top().poly = f1;
1448                                         ModFactors mf;
1449                                         mf.factors = newfactors2;
1450                                         mf.poly = f2;
1451                                         tocheck.push(mf);
1452                                         break;
1453                                 }
1454                         }
1455                         else {
1456                                 if ( !part.next() ) {
1457                                         result *= upoly_to_ex(tocheck.top().poly, x);
1458                                         tocheck.pop();
1459                                         break;
1460                                 }
1461                         }
1462                 }
1463         }
1464
1465         return unit * cont * result;
1466 }
1467
1468 static inline ex factor_univariate(const ex& poly, const ex& x)
1469 {
1470         unsigned int prime;
1471         return factor_univariate(poly, x, prime);
1472 }
1473
1474 struct EvalPoint
1475 {
1476         ex x;
1477         int evalpoint;
1478 };
1479
1480 // forward declaration
1481 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);
1482
1483 static upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, unsigned int k)
1484 {
1485         const size_t r = a.size();
1486         cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1487         upvec q(r-1);
1488         q[r-2] = a[r-1];
1489         for ( size_t j=r-2; j>=1; --j ) {
1490                 q[j-1] = a[j] * q[j];
1491         }
1492         umodpoly beta(1, R->one());
1493         upvec s;
1494         for ( size_t j=1; j<r; ++j ) {
1495                 vector<ex> mdarg(2);
1496                 mdarg[0] = umodpoly_to_ex(q[j-1], x);
1497                 mdarg[1] = umodpoly_to_ex(a[j-1], x);
1498                 vector<EvalPoint> empty;
1499                 vector<ex> exsigma = multivar_diophant(mdarg, x, umodpoly_to_ex(beta, x), empty, 0, p, k);
1500                 umodpoly sigma1;
1501                 umodpoly_from_ex(sigma1, exsigma[0], x, R);
1502                 umodpoly sigma2;
1503                 umodpoly_from_ex(sigma2, exsigma[1], x, R);
1504                 beta = sigma1;
1505                 s.push_back(sigma2);
1506         }
1507         s.push_back(beta);
1508         return s;
1509 }
1510
1511 /**
1512  *  Assert: a not empty.
1513  */
1514 static void change_modulus(const cl_modint_ring& R, umodpoly& a)
1515 {
1516         if ( a.empty() ) return;
1517         cl_modint_ring oldR = a[0].ring();
1518         umodpoly::iterator i = a.begin(), end = a.end();
1519         for ( ; i!=end; ++i ) {
1520                 *i = R->canonhom(oldR->retract(*i));
1521         }
1522         canonicalize(a);
1523 }
1524
1525 static void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned int p, unsigned int k, umodpoly& s_, umodpoly& t_)
1526 {
1527         cl_modint_ring R = find_modint_ring(p);
1528         umodpoly amod = a;
1529         change_modulus(R, amod);
1530         umodpoly bmod = b;
1531         change_modulus(R, bmod);
1532
1533         umodpoly smod;
1534         umodpoly tmod;
1535         exteuclid(amod, bmod, smod, tmod);
1536
1537         cl_modint_ring Rpk = find_modint_ring(expt_pos(cl_I(p),k));
1538         umodpoly s = smod;
1539         change_modulus(Rpk, s);
1540         umodpoly t = tmod;
1541         change_modulus(Rpk, t);
1542
1543         cl_I modulus(p);
1544         umodpoly one(1, Rpk->one());
1545         for ( size_t j=1; j<k; ++j ) {
1546                 umodpoly e = one - a * s - b * t;
1547                 reduce_coeff(e, modulus);
1548                 umodpoly c = e;
1549                 change_modulus(R, c);
1550                 umodpoly sigmabar = smod * c;
1551                 umodpoly taubar = tmod * c;
1552                 umodpoly sigma, q;
1553                 remdiv(sigmabar, bmod, sigma, q);
1554                 umodpoly tau = taubar + q * amod;
1555                 umodpoly sadd = sigma;
1556                 change_modulus(Rpk, sadd);
1557                 cl_MI modmodulus(Rpk, modulus);
1558                 s = s + sadd * modmodulus;
1559                 umodpoly tadd = tau;
1560                 change_modulus(Rpk, tadd);
1561                 t = t + tadd * modmodulus;
1562                 modulus = modulus * p;
1563         }
1564
1565         s_ = s; t_ = t;
1566 }
1567
1568 static upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsigned int p, unsigned int k)
1569 {
1570         cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1571
1572         const size_t r = a.size();
1573         upvec result;
1574         if ( r > 2 ) {
1575                 upvec s = multiterm_eea_lift(a, x, p, k);
1576                 for ( size_t j=0; j<r; ++j ) {
1577                         umodpoly bmod = umodpoly_to_umodpoly(s[j], R, m);
1578                         umodpoly buf;
1579                         rem(bmod, a[j], buf);
1580                         result.push_back(buf);
1581                 }
1582         }
1583         else {
1584                 umodpoly s, t;
1585                 eea_lift(a[1], a[0], x, p, k, s, t);
1586                 umodpoly bmod = umodpoly_to_umodpoly(s, R, m);
1587                 umodpoly buf, q;
1588                 remdiv(bmod, a[0], buf, q);
1589                 result.push_back(buf);
1590                 umodpoly t1mod = umodpoly_to_umodpoly(t, R, m);
1591                 buf = t1mod + q * a[1];
1592                 result.push_back(buf);
1593         }
1594
1595         return result;
1596 }
1597
1598 struct make_modular_map : public map_function {
1599         cl_modint_ring R;
1600         make_modular_map(const cl_modint_ring& R_) : R(R_) { }
1601         ex operator()(const ex& e)
1602         {
1603                 if ( is_a<add>(e) || is_a<mul>(e) ) {
1604                         return e.map(*this);
1605                 }
1606                 else if ( is_a<numeric>(e) ) {
1607                         numeric mod(R->modulus);
1608                         numeric halfmod = (mod-1)/2;
1609                         cl_MI emod = R->canonhom(the<cl_I>(ex_to<numeric>(e).to_cl_N()));
1610                         numeric n(R->retract(emod));
1611                         if ( n > halfmod ) {
1612                                 return n-mod;
1613                         }
1614                         else {
1615                                 return n;
1616                         }
1617                 }
1618                 return e;
1619         }
1620 };
1621
1622 static ex make_modular(const ex& e, const cl_modint_ring& R)
1623 {
1624         make_modular_map map(R);
1625         return map(e.expand());
1626 }
1627
1628 static vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex& c, const vector<EvalPoint>& I,
1629                                     unsigned int d, unsigned int p, unsigned int k)
1630 {
1631         vector<ex> a = a_;
1632
1633         const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1634         const size_t r = a.size();
1635         const size_t nu = I.size() + 1;
1636
1637         vector<ex> sigma;
1638         if ( nu > 1 ) {
1639                 ex xnu = I.back().x;
1640                 int alphanu = I.back().evalpoint;
1641
1642                 ex A = 1;
1643                 for ( size_t i=0; i<r; ++i ) {
1644                         A *= a[i];
1645                 }
1646                 vector<ex> b(r);
1647                 for ( size_t i=0; i<r; ++i ) {
1648                         b[i] = normal(A / a[i]);
1649                 }
1650
1651                 vector<ex> anew = a;
1652                 for ( size_t i=0; i<r; ++i ) {
1653                         anew[i] = anew[i].subs(xnu == alphanu);
1654                 }
1655                 ex cnew = c.subs(xnu == alphanu);
1656                 vector<EvalPoint> Inew = I;
1657                 Inew.pop_back();
1658                 sigma = multivar_diophant(anew, x, cnew, Inew, d, p, k);
1659
1660                 ex buf = c;
1661                 for ( size_t i=0; i<r; ++i ) {
1662                         buf -= sigma[i] * b[i];
1663                 }
1664                 ex e = make_modular(buf, R);
1665
1666                 ex monomial = 1;
1667                 for ( size_t m=1; !e.is_zero() && e.has(xnu) && m<=d; ++m ) {
1668                         monomial *= (xnu - alphanu);
1669                         monomial = expand(monomial);
1670                         ex cm = e.diff(ex_to<symbol>(xnu), m).subs(xnu==alphanu) / factorial(m);
1671                         cm = make_modular(cm, R);
1672                         if ( !cm.is_zero() ) {
1673                                 vector<ex> delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k);
1674                                 ex buf = e;
1675                                 for ( size_t j=0; j<delta_s.size(); ++j ) {
1676                                         delta_s[j] *= monomial;
1677                                         sigma[j] += delta_s[j];
1678                                         buf -= delta_s[j] * b[j];
1679                                 }
1680                                 e = make_modular(buf, R);
1681                         }
1682                 }
1683         }
1684         else {
1685                 upvec amod;
1686                 for ( size_t i=0; i<a.size(); ++i ) {
1687                         umodpoly up;
1688                         umodpoly_from_ex(up, a[i], x, R);
1689                         amod.push_back(up);
1690                 }
1691
1692                 sigma.insert(sigma.begin(), r, 0);
1693                 size_t nterms;
1694                 ex z;
1695                 if ( is_a<add>(c) ) {
1696                         nterms = c.nops();
1697                         z = c.op(0);
1698                 }
1699                 else {
1700                         nterms = 1;
1701                         z = c;
1702                 }
1703                 for ( size_t i=0; i<nterms; ++i ) {
1704                         int m = z.degree(x);
1705                         cl_I cm = the<cl_I>(ex_to<numeric>(z.lcoeff(x)).to_cl_N());
1706                         upvec delta_s = univar_diophant(amod, x, m, p, k);
1707                         cl_MI modcm;
1708                         cl_I poscm = cm;
1709                         while ( poscm < 0 ) {
1710                                 poscm = poscm + expt_pos(cl_I(p),k);
1711                         }
1712                         modcm = cl_MI(R, poscm);
1713                         for ( size_t j=0; j<delta_s.size(); ++j ) {
1714                                 delta_s[j] = delta_s[j] * modcm;
1715                                 sigma[j] = sigma[j] + umodpoly_to_ex(delta_s[j], x);
1716                         }
1717                         if ( nterms > 1 ) {
1718                                 z = c.op(i+1);
1719                         }
1720                 }
1721         }
1722
1723         for ( size_t i=0; i<sigma.size(); ++i ) {
1724                 sigma[i] = make_modular(sigma[i], R);
1725         }
1726
1727         return sigma;
1728 }
1729
1730 #ifdef DEBUGFACTOR
1731 ostream& operator<<(ostream& o, const vector<EvalPoint>& v)
1732 {
1733         for ( size_t i=0; i<v.size(); ++i ) {
1734                 o << "(" << v[i].x << "==" << v[i].evalpoint << ") ";
1735         }
1736         return o;
1737 }
1738 #endif // def DEBUGFACTOR
1739
1740 static ex hensel_multivar(const ex& a, const ex& x, const vector<EvalPoint>& I, unsigned int p, const cl_I& l, const upvec& u, const vector<ex>& lcU)
1741 {
1742         const size_t nu = I.size() + 1;
1743         const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),l));
1744
1745         vector<ex> A(nu);
1746         A[nu-1] = a;
1747
1748         for ( size_t j=nu; j>=2; --j ) {
1749                 ex x = I[j-2].x;
1750                 int alpha = I[j-2].evalpoint;
1751                 A[j-2] = A[j-1].subs(x==alpha);
1752                 A[j-2] = make_modular(A[j-2], R);
1753         }
1754
1755         int maxdeg = a.degree(I.front().x);
1756         for ( size_t i=1; i<I.size(); ++i ) {
1757                 int maxdeg2 = a.degree(I[i].x);
1758                 if ( maxdeg2 > maxdeg ) maxdeg = maxdeg2;
1759         }
1760
1761         const size_t n = u.size();
1762         vector<ex> U(n);
1763         for ( size_t i=0; i<n; ++i ) {
1764                 U[i] = umodpoly_to_ex(u[i], x);
1765         }
1766
1767         for ( size_t j=2; j<=nu; ++j ) {
1768                 vector<ex> U1 = U;
1769                 ex monomial = 1;
1770                 for ( size_t m=0; m<n; ++m) {
1771                         if ( lcU[m] != 1 ) {
1772                                 ex coef = lcU[m];
1773                                 for ( size_t i=j-1; i<nu-1; ++i ) {
1774                                         coef = coef.subs(I[i].x == I[i].evalpoint);
1775                                 }
1776                                 coef = make_modular(coef, R);
1777                                 int deg = U[m].degree(x);
1778                                 U[m] = U[m] - U[m].lcoeff(x) * pow(x,deg) + coef * pow(x,deg);
1779                         }
1780                 }
1781                 ex Uprod = 1;
1782                 for ( size_t i=0; i<n; ++i ) {
1783                         Uprod *= U[i];
1784                 }
1785                 ex e = expand(A[j-1] - Uprod);
1786
1787                 vector<EvalPoint> newI;
1788                 for ( size_t i=1; i<=j-2; ++i ) {
1789                         newI.push_back(I[i-1]);
1790                 }
1791
1792                 ex xj = I[j-2].x;
1793                 int alphaj = I[j-2].evalpoint;
1794                 size_t deg = A[j-1].degree(xj);
1795                 for ( size_t k=1; k<=deg; ++k ) {
1796                         if ( !e.is_zero() ) {
1797                                 monomial *= (xj - alphaj);
1798                                 monomial = expand(monomial);
1799                                 ex dif = e.diff(ex_to<symbol>(xj), k);
1800                                 ex c = dif.subs(xj==alphaj) / factorial(k);
1801                                 if ( !c.is_zero() ) {
1802                                         vector<ex> deltaU = multivar_diophant(U1, x, c, newI, maxdeg, p, cl_I_to_uint(l));
1803                                         for ( size_t i=0; i<n; ++i ) {
1804                                                 deltaU[i] *= monomial;
1805                                                 U[i] += deltaU[i];
1806                                                 U[i] = make_modular(U[i], R);
1807                                         }
1808                                         ex Uprod = 1;
1809                                         for ( size_t i=0; i<n; ++i ) {
1810                                                 Uprod *= U[i];
1811                                         }
1812                                         e = A[j-1] - Uprod;
1813                                         e = make_modular(e, R);
1814                                 }
1815                         }
1816                 }
1817         }
1818
1819         ex acand = 1;
1820         for ( size_t i=0; i<U.size(); ++i ) {
1821                 acand *= U[i];
1822         }
1823         if ( expand(a-acand).is_zero() ) {
1824                 lst res;
1825                 for ( size_t i=0; i<U.size(); ++i ) {
1826                         res.append(U[i]);
1827                 }
1828                 return res;
1829         }
1830         else {
1831                 lst res;
1832                 return lst();
1833         }
1834 }
1835
1836 static ex put_factors_into_lst(const ex& e)
1837 {
1838         lst result;
1839
1840         if ( is_a<numeric>(e) ) {
1841                 result.append(e);
1842                 return result;
1843         }
1844         if ( is_a<power>(e) ) {
1845                 result.append(1);
1846                 result.append(e.op(0));
1847                 return result;
1848         }
1849         if ( is_a<symbol>(e) || is_a<add>(e) ) {
1850                 result.append(1);
1851                 result.append(e);
1852                 return result;
1853         }
1854         if ( is_a<mul>(e) ) {
1855                 ex nfac = 1;
1856                 for ( size_t i=0; i<e.nops(); ++i ) {
1857                         ex op = e.op(i);
1858                         if ( is_a<numeric>(op) ) {
1859                                 nfac = op;
1860                         }
1861                         if ( is_a<power>(op) ) {
1862                                 result.append(op.op(0));
1863                         }
1864                         if ( is_a<symbol>(op) || is_a<add>(op) ) {
1865                                 result.append(op);
1866                         }
1867                 }
1868                 result.prepend(nfac);
1869                 return result;
1870         }
1871         throw runtime_error("put_factors_into_lst: bad term.");
1872 }
1873
1874 #ifdef DEBUGFACTOR
1875 ostream& operator<<(ostream& o, const vector<numeric>& v)
1876 {
1877         for ( size_t i=0; i<v.size(); ++i ) {
1878                 o << v[i] << " ";
1879         }
1880         return o;
1881 }
1882 #endif // def DEBUGFACTOR
1883
1884 /** Checks whether in a set of numbers each has a unique prime factor.
1885  *
1886  *  @param[in]  f  list of numbers to check
1887  *  @return        true: if number set is bad, false: otherwise
1888  */
1889 static bool checkdivisors(const lst& f)
1890 {
1891         const int k = f.nops();
1892         numeric q, r;
1893         vector<numeric> d(k);
1894         d[0] = ex_to<numeric>(abs(f.op(0)));
1895         for ( int i=1; i<k; ++i ) {
1896                 q = ex_to<numeric>(abs(f.op(i)));
1897                 for ( int j=i-1; j>=0; --j ) {
1898                         r = d[j];
1899                         do {
1900                                 r = gcd(r, q);
1901                                 q = q/r;
1902                         } while ( r != 1 );
1903                         if ( q == 1 ) {
1904                                 return true;
1905                         }
1906                 }
1907                 d[i] = q;
1908         }
1909         return false;
1910 }
1911
1912 /** Generates a set of evaluation points for a multivariate polynomial.
1913  *  The set fulfills the following conditions:
1914  *  1. lcoeff(evaluated_polynomial) does not vanish
1915  *  2. factors of lcoeff(evaluated_polynomial) have each a unique prime factor
1916  *  3. evaluated_polynomial is square free
1917  *  See [W1] for more details.
1918  *
1919  *  @param[in]     u        multivariate polynomial to be factored
1920  *  @param[in]     vn       leading coefficient of u in x (x==first symbol in syms)
1921  *  @param[in]     syms     set of symbols that appear in u
1922  *  @param[in]     f        lst containing the factors of the leading coefficient vn
1923  *  @param[in,out] modulus  integer modulus for random number generation (i.e. |a_i| < modulus)
1924  *  @param[out]    u0       returns the evaluated (univariate) polynomial
1925  *  @param[out]    a        returns the valid evaluation points. must have initial size equal
1926  *                          number of symbols-1 before calling generate_set
1927  */
1928 static void generate_set(const ex& u, const ex& vn, const exset& syms, const lst& f,
1929                          numeric& modulus, ex& u0, vector<numeric>& a)
1930 {
1931         const ex& x = *syms.begin();
1932         while ( true ) {
1933                 ++modulus;
1934                 /* generate a set of integers ... */
1935                 u0 = u;
1936                 ex vna = vn;
1937                 ex vnatry;
1938                 exset::const_iterator s = syms.begin();
1939                 ++s;
1940                 for ( size_t i=0; i<a.size(); ++i ) {
1941                         do {
1942                                 a[i] = mod(numeric(rand()), 2*modulus) - modulus;
1943                                 vnatry = vna.subs(*s == a[i]);
1944                                 /* ... for which the leading coefficient doesn't vanish ... */
1945                         } while ( vnatry == 0 );
1946                         vna = vnatry;
1947                         u0 = u0.subs(*s == a[i]);
1948                         ++s;
1949                 }
1950                 /* ... for which u0 is square free ... */
1951                 ex g = gcd(u0, u0.diff(ex_to<symbol>(x)));
1952                 if ( !is_a<numeric>(g) ) {
1953                         continue;
1954                 }
1955                 if ( !is_a<numeric>(vn) ) {
1956                         /* ... and for which the evaluated factors have each an unique prime factor */
1957                         lst fnum = f;
1958                         fnum.let_op(0) = fnum.op(0) * u0.content(x);
1959                         for ( size_t i=1; i<fnum.nops(); ++i ) {
1960                                 if ( !is_a<numeric>(fnum.op(i)) ) {
1961                                         s = syms.begin();
1962                                         ++s;
1963                                         for ( size_t j=0; j<a.size(); ++j, ++s ) {
1964                                                 fnum.let_op(i) = fnum.op(i).subs(*s == a[j]);
1965                                         }
1966                                 }
1967                         }
1968                         if ( checkdivisors(fnum) ) {
1969                                 continue;
1970                         }
1971                 }
1972                 /* ok, we have a valid set now */
1973                 return;
1974         }
1975 }
1976
1977 // forward declaration
1978 static ex factor_sqrfree(const ex& poly);
1979
1980 /**
1981  *  ASSERT: poly is expanded
1982  */
1983 static ex factor_multivariate(const ex& poly, const exset& syms)
1984 {
1985         exset::const_iterator s;
1986         const ex& x = *syms.begin();
1987
1988         /* make polynomial primitive */
1989         ex p = poly.collect(x);
1990         ex cont = p.lcoeff(x);
1991         for ( int i=p.degree(x)-1; i>=p.ldegree(x); --i ) {
1992                 cont = gcd(cont, p.coeff(x,i));
1993                 if ( cont == 1 ) break;
1994         }
1995         ex pp = expand(normal(p / cont));
1996         if ( !is_a<numeric>(cont) ) {
1997                 return factor_sqrfree(cont) * factor_sqrfree(pp);
1998         }
1999
2000         /* factor leading coefficient */
2001         ex vn = pp.collect(x).lcoeff(x);
2002         ex vnlst;
2003         if ( is_a<numeric>(vn) ) {
2004                 vnlst = lst(vn);
2005         }
2006         else {
2007                 ex vnfactors = factor(vn);
2008                 vnlst = put_factors_into_lst(vnfactors);
2009         }
2010
2011         const unsigned int maxtrials = 3;
2012         numeric modulus = (vnlst.nops() > 3) ? vnlst.nops() : 3;
2013         vector<numeric> a(syms.size()-1, 0);
2014
2015         /* try now to factorize until we are successful */
2016         while ( true ) {
2017
2018                 unsigned int trialcount = 0;
2019                 unsigned int prime;
2020                 int factor_count = 0;
2021                 int min_factor_count = -1;
2022                 ex u, delta;
2023                 ex ufac, ufaclst;
2024
2025                 /* try several evaluation points to reduce the number of modular factors */
2026                 while ( trialcount < maxtrials ) {
2027
2028                         /* generate a set of valid evaluation points */
2029                         generate_set(pp, vn, syms, ex_to<lst>(vnlst), modulus, u, a);
2030
2031                         ufac = factor_univariate(u, x, prime);
2032                         ufaclst = put_factors_into_lst(ufac);
2033                         factor_count = ufaclst.nops()-1;
2034                         delta = ufaclst.op(0);
2035
2036                         if ( factor_count <= 1 ) {
2037                                 /* irreducible */
2038                                 return poly;
2039                         }
2040                         if ( min_factor_count < 0 ) {
2041                                 /* first time here */
2042                                 min_factor_count = factor_count;
2043                         }
2044                         else if ( min_factor_count == factor_count ) {
2045                                 /* one less to try */
2046                                 ++trialcount;
2047                         }
2048                         else if ( min_factor_count > factor_count ) {
2049                                 /* new minimum, reset trial counter */
2050                                 min_factor_count = factor_count;
2051                                 trialcount = 0;
2052                         }
2053                 }
2054
2055                 // determine true leading coefficients for the Hensel lifting
2056                 vector<ex> C(factor_count);
2057                 if ( is_a<numeric>(vn) ) {
2058                         for ( size_t i=1; i<ufaclst.nops(); ++i ) {
2059                                 C[i-1] = ufaclst.op(i).lcoeff(x);
2060                         }
2061                 }
2062                 else {
2063                         vector<numeric> ftilde(vnlst.nops()-1);
2064                         for ( size_t i=0; i<ftilde.size(); ++i ) {
2065                                 ex ft = vnlst.op(i+1);
2066                                 s = syms.begin();
2067                                 ++s;
2068                                 for ( size_t j=0; j<a.size(); ++j ) {
2069                                         ft = ft.subs(*s == a[j]);
2070                                         ++s;
2071                                 }
2072                                 ftilde[i] = ex_to<numeric>(ft);
2073                         }
2074
2075                         vector<bool> used_flag(ftilde.size(), false);
2076                         vector<ex> D(factor_count, 1);
2077                         if ( delta == 1 ) {
2078                                 for ( int i=0; i<factor_count; ++i ) {
2079                                         numeric prefac = ex_to<numeric>(ufaclst.op(i+1).lcoeff(x));
2080                                         for ( int j=ftilde.size()-1; j>=0; --j ) {
2081                                                 int count = 0;
2082                                                 while ( irem(prefac, ftilde[j]) == 0 ) {
2083                                                         prefac = iquo(prefac, ftilde[j]);
2084                                                         ++count;
2085                                                 }
2086                                                 if ( count ) {
2087                                                         used_flag[j] = true;
2088                                                         D[i] = D[i] * pow(vnlst.op(j+1), count);
2089                                                 }
2090                                         }
2091                                         C[i] = D[i] * prefac;
2092                                 }
2093                         }
2094                         else {
2095                                 for ( int i=0; i<factor_count; ++i ) {
2096                                         numeric prefac = ex_to<numeric>(ufaclst.op(i+1).lcoeff(x));
2097                                         for ( int j=ftilde.size()-1; j>=0; --j ) {
2098                                                 int count = 0;
2099                                                 while ( irem(prefac, ftilde[j]) == 0 ) {
2100                                                         prefac = iquo(prefac, ftilde[j]);
2101                                                         ++count;
2102                                                 }
2103                                                 while ( irem(ex_to<numeric>(delta)*prefac, ftilde[j]) == 0 ) {
2104                                                         numeric g = gcd(prefac, ex_to<numeric>(ftilde[j]));
2105                                                         prefac = iquo(prefac, g);
2106                                                         delta = delta / (ftilde[j]/g);
2107                                                         ufaclst.let_op(i+1) = ufaclst.op(i+1) * (ftilde[j]/g);
2108                                                         ++count;
2109                                                 }
2110                                                 if ( count ) {
2111                                                         used_flag[j] = true;
2112                                                         D[i] = D[i] * pow(vnlst.op(j+1), count);
2113                                                 }
2114                                         }
2115                                         C[i] = D[i] * prefac;
2116                                 }
2117                         }
2118
2119                         bool some_factor_unused = false;
2120                         for ( size_t i=0; i<used_flag.size(); ++i ) {
2121                                 if ( !used_flag[i] ) {
2122                                         some_factor_unused = true;
2123                                         break;
2124                                 }
2125                         }
2126                         if ( some_factor_unused ) {
2127                                 continue;
2128                         }
2129                 }
2130
2131                 if ( delta != 1 ) {
2132                         C[0] = C[0] * delta;
2133                         ufaclst.let_op(1) = ufaclst.op(1) * delta;
2134                 }
2135
2136                 EvalPoint ep;
2137                 vector<EvalPoint> epv;
2138                 s = syms.begin();
2139                 ++s;
2140                 for ( size_t i=0; i<a.size(); ++i ) {
2141                         ep.x = *s++;
2142                         ep.evalpoint = a[i].to_int();
2143                         epv.push_back(ep);
2144                 }
2145
2146                 // calc bound p^l
2147                 int maxdeg = 0;
2148                 for ( int i=1; i<=factor_count; ++i ) {
2149                         if ( ufaclst.op(i).degree(x) > maxdeg ) {
2150                                 maxdeg = ufaclst[i].degree(x);
2151                         }
2152                 }
2153                 cl_I B = 2*calc_bound(u, x, maxdeg);
2154                 cl_I l = 1;
2155                 cl_I pl = prime;
2156                 while ( pl < B ) {
2157                         l = l + 1;
2158                         pl = pl * prime;
2159                 }
2160                 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(prime),l));
2161                 upvec modfactors(ufaclst.nops()-1);
2162                 for ( size_t i=1; i<ufaclst.nops(); ++i ) {
2163                         umodpoly_from_ex(modfactors[i-1], ufaclst.op(i), x, R);
2164                 }
2165
2166                 ex res = hensel_multivar(pp, x, epv, prime, l, modfactors, C);
2167                 if ( res != lst() ) {
2168                         ex result = cont;
2169                         for ( size_t i=0; i<res.nops(); ++i ) {
2170                                 result *= res.op(i).content(x) * res.op(i).unit(x);
2171                                 result *= res.op(i).primpart(x);
2172                         }
2173                         return result;
2174                 }
2175         }
2176 }
2177
2178 struct find_symbols_map : public map_function {
2179         exset syms;
2180         ex operator()(const ex& e)
2181         {
2182                 if ( is_a<symbol>(e) ) {
2183                         syms.insert(e);
2184                         return e;
2185                 }
2186                 return e.map(*this);
2187         }
2188 };
2189
2190 static ex factor_sqrfree(const ex& poly)
2191 {
2192         // determine all symbols in poly
2193         find_symbols_map findsymbols;
2194         findsymbols(poly);
2195         if ( findsymbols.syms.size() == 0 ) {
2196                 return poly;
2197         }
2198
2199         if ( findsymbols.syms.size() == 1 ) {
2200                 // univariate case
2201                 const ex& x = *(findsymbols.syms.begin());
2202                 if ( poly.ldegree(x) > 0 ) {
2203                         int ld = poly.ldegree(x);
2204                         ex res = factor_univariate(expand(poly/pow(x, ld)), x);
2205                         return res * pow(x,ld);
2206                 }
2207                 else {
2208                         ex res = factor_univariate(poly, x);
2209                         return res;
2210                 }
2211         }
2212
2213         // multivariate case
2214         ex res = factor_multivariate(poly, findsymbols.syms);
2215         return res;
2216 }
2217
2218 struct apply_factor_map : public map_function {
2219         unsigned options;
2220         apply_factor_map(unsigned options_) : options(options_) { }
2221         ex operator()(const ex& e)
2222         {
2223                 if ( e.info(info_flags::polynomial) ) {
2224                         return factor(e, options);
2225                 }
2226                 if ( is_a<add>(e) ) {
2227                         ex s1, s2;
2228                         for ( size_t i=0; i<e.nops(); ++i ) {
2229                                 if ( e.op(i).info(info_flags::polynomial) ) {
2230                                         s1 += e.op(i);
2231                                 }
2232                                 else {
2233                                         s2 += e.op(i);
2234                                 }
2235                         }
2236                         s1 = s1.eval();
2237                         s2 = s2.eval();
2238                         return factor(s1, options) + s2.map(*this);
2239                 }
2240                 return e.map(*this);
2241         }
2242 };
2243
2244 } // anonymous namespace
2245
2246 ex factor(const ex& poly, unsigned options)
2247 {
2248         // check arguments
2249         if ( !poly.info(info_flags::polynomial) ) {
2250                 if ( options & factor_options::all ) {
2251                         options &= ~factor_options::all;
2252                         apply_factor_map factor_map(options);
2253                         return factor_map(poly);
2254                 }
2255                 return poly;
2256         }
2257
2258         // determine all symbols in poly
2259         find_symbols_map findsymbols;
2260         findsymbols(poly);
2261         if ( findsymbols.syms.size() == 0 ) {
2262                 return poly;
2263         }
2264         lst syms;
2265         exset::const_iterator i=findsymbols.syms.begin(), end=findsymbols.syms.end();
2266         for ( ; i!=end; ++i ) {
2267                 syms.append(*i);
2268         }
2269
2270         // make poly square free
2271         ex sfpoly = sqrfree(poly.expand(), syms);
2272
2273         // factorize the square free components
2274         if ( is_a<power>(sfpoly) ) {
2275                 // case: (polynomial)^exponent
2276                 const ex& base = sfpoly.op(0);
2277                 if ( !is_a<add>(base) ) {
2278                         // simple case: (monomial)^exponent
2279                         return sfpoly;
2280                 }
2281                 ex f = factor_sqrfree(base);
2282                 return pow(f, sfpoly.op(1));
2283         }
2284         if ( is_a<mul>(sfpoly) ) {
2285                 // case: multiple factors
2286                 ex res = 1;
2287                 for ( size_t i=0; i<sfpoly.nops(); ++i ) {
2288                         const ex& t = sfpoly.op(i);
2289                         if ( is_a<power>(t) ) {
2290                                 const ex& base = t.op(0);
2291                                 if ( !is_a<add>(base) ) {
2292                                         res *= t;
2293                                 }
2294                                 else {
2295                                         ex f = factor_sqrfree(base);
2296                                         res *= pow(f, t.op(1));
2297                                 }
2298                         }
2299                         else if ( is_a<add>(t) ) {
2300                                 ex f = factor_sqrfree(t);
2301                                 res *= f;
2302                         }
2303                         else {
2304                                 res *= t;
2305                         }
2306                 }
2307                 return res;
2308         }
2309         if ( is_a<symbol>(sfpoly) ) {
2310                 return poly;
2311         }
2312         // case: (polynomial)
2313         ex f = factor_sqrfree(sfpoly);
2314         return f;
2315 }
2316
2317 } // namespace GiNaC
2318
2319 #ifdef DEBUGFACTOR
2320 #include "test.h"
2321 #endif