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