ginac/factor.cpp

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