ginac/factor.cpp

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