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