ginac/factor.cpp

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