ginac/factor.cpp

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