ginac/normal.cpp

   1 /** @file normal.cpp
   2  *
   3  *  This file implements several functions that work on univariate and
   4  *  multivariate polynomials and rational functions.
   5  *  These functions include polynomial quotient and remainder, GCD and LCM
   6  *  computation, square-free factorization and rational function normalization. */
   7
   8 /*
   9  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
  10  *
  11  *  This program is free software; you can redistribute it and/or modify
  12  *  it under the terms of the GNU General Public License as published by
  13  *  the Free Software Foundation; either version 2 of the License, or
  14  *  (at your option) any later version.
  15  *
  16  *  This program is distributed in the hope that it will be useful,
  17  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19  *  GNU General Public License for more details.
  20  *
  21  *  You should have received a copy of the GNU General Public License
  22  *  along with this program; if not, write to the Free Software
  23  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  24  */
  25
  26 #include <stdexcept>
  27 #include <algorithm>
  28 #include <map>
  29
  30 #include "normal.h"
  31 #include "basic.h"
  32 #include "ex.h"
  33 #include "add.h"
  34 #include "constant.h"
  35 #include "expairseq.h"
  36 #include "fail.h"
  37 #include "inifcns.h"
  38 #include "lst.h"
  39 #include "mul.h"
  40 #include "ncmul.h"
  41 #include "numeric.h"
  42 #include "power.h"
  43 #include "relational.h"
  44 #include "pseries.h"
  45 #include "symbol.h"
  46 #include "utils.h"
  47
  48 #ifndef NO_NAMESPACE_GINAC
  49 namespace GiNaC {
  50 #endif // ndef NO_NAMESPACE_GINAC
  51
  52 // If comparing expressions (ex::compare()) is fast, you can set this to 1.
  53 // Some routines like quo(), rem() and gcd() will then return a quick answer
  54 // when they are called with two identical arguments.
  55 #define FAST_COMPARE 1
  56
  57 // Set this if you want divide_in_z() to use remembering
  58 #define USE_REMEMBER 0
  59
  60 // Set this if you want divide_in_z() to use trial division followed by
  61 // polynomial interpolation (always slower except for completely dense
  62 // polynomials)
  63 #define USE_TRIAL_DIVISION 0
  64
  65 // Set this to enable some statistical output for the GCD routines
  66 #define STATISTICS 0
  67
  68
  69 #if STATISTICS
  70 // Statistics variables
  71 static int gcd_called = 0;
  72 static int sr_gcd_called = 0;
  73 static int heur_gcd_called = 0;
  74 static int heur_gcd_failed = 0;
  75
  76 // Print statistics at end of program
  77 static struct _stat_print {
  78         _stat_print() {}
  79         ~_stat_print() {
  80                 cout << "gcd() called " << gcd_called << " times\n";
  81                 cout << "sr_gcd() called " << sr_gcd_called << " times\n";
  82                 cout << "heur_gcd() called " << heur_gcd_called << " times\n";
  83                 cout << "heur_gcd() failed " << heur_gcd_failed << " times\n";
  84         }
  85 } stat_print;
  86 #endif
  87
  88
  89 /** Return pointer to first symbol found in expression.  Due to GiNaC´s
  90  *  internal ordering of terms, it may not be obvious which symbol this
  91  *  function returns for a given expression.
  92  *
  93  *  @param e  expression to search
  94  *  @param x  pointer to first symbol found (returned)
  95  *  @return "false" if no symbol was found, "true" otherwise */
  96 static bool get_first_symbol(const ex &e, const symbol *&x)
  97 {
  98         if (is_ex_exactly_of_type(e, symbol)) {
  99                 x = static_cast<symbol *>(e.bp);
 100                 return true;
 101         } else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) {
 102                 for (unsigned i=0; i<e.nops(); i++)
 103                         if (get_first_symbol(e.op(i), x))
 104                                 return true;
 105         } else if (is_ex_exactly_of_type(e, power)) {
 106                 if (get_first_symbol(e.op(0), x))
 107                         return true;
 108         }
 109         return false;
 110 }
 111
 112
 113 /*
 114  *  Statistical information about symbols in polynomials
 115  */
 116
 117 /** This structure holds information about the highest and lowest degrees
 118  *  in which a symbol appears in two multivariate polynomials "a" and "b".
 119  *  A vector of these structures with information about all symbols in
 120  *  two polynomials can be created with the function get_symbol_stats().
 121  *
 122  *  @see get_symbol_stats */
 123 struct sym_desc {
 124         /** Pointer to symbol */
 125         const symbol *sym;
 126
 127         /** Highest degree of symbol in polynomial "a" */
 128         int deg_a;
 129
 130         /** Highest degree of symbol in polynomial "b" */
 131         int deg_b;
 132
 133         /** Lowest degree of symbol in polynomial "a" */
 134         int ldeg_a;
 135
 136         /** Lowest degree of symbol in polynomial "b" */
 137         int ldeg_b;
 138
 139         /** Maximum of deg_a and deg_b (Used for sorting) */
 140         int max_deg;
 141
 142         /** Maximum number of terms of leading coefficient of symbol in both polynomials */
 143         int max_lcnops;
 144
 145         /** Commparison operator for sorting */
 146         bool operator<(const sym_desc &x) const
 147         {
 148                 if (max_deg == x.max_deg)
 149                         return max_lcnops < x.max_lcnops;
 150                 else
 151                         return max_deg < x.max_deg;
 152         }
 153 };
 154
 155 // Vector of sym_desc structures
 156 typedef std::vector<sym_desc> sym_desc_vec;
 157
 158 // Add symbol the sym_desc_vec (used internally by get_symbol_stats())
 159 static void add_symbol(const symbol *s, sym_desc_vec &v)
 160 {
 161         sym_desc_vec::iterator it = v.begin(), itend = v.end();
 162         while (it != itend) {
 163                 if (it->sym->compare(*s) == 0)  // If it's already in there, don't add it a second time
 164                         return;
 165                 it++;
 166         }
 167         sym_desc d;
 168         d.sym = s;
 169         v.push_back(d);
 170 }
 171
 172 // Collect all symbols of an expression (used internally by get_symbol_stats())
 173 static void collect_symbols(const ex &e, sym_desc_vec &v)
 174 {
 175         if (is_ex_exactly_of_type(e, symbol)) {
 176                 add_symbol(static_cast<symbol *>(e.bp), v);
 177         } else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) {
 178                 for (unsigned i=0; i<e.nops(); i++)
 179                         collect_symbols(e.op(i), v);
 180         } else if (is_ex_exactly_of_type(e, power)) {
 181                 collect_symbols(e.op(0), v);
 182         }
 183 }
 184
 185 /** Collect statistical information about symbols in polynomials.
 186  *  This function fills in a vector of "sym_desc" structs which contain
 187  *  information about the highest and lowest degrees of all symbols that
 188  *  appear in two polynomials. The vector is then sorted by minimum
 189  *  degree (lowest to highest). The information gathered by this
 190  *  function is used by the GCD routines to identify trivial factors
 191  *  and to determine which variable to choose as the main variable
 192  *  for GCD computation.
 193  *
 194  *  @param a  first multivariate polynomial
 195  *  @param b  second multivariate polynomial
 196  *  @param v  vector of sym_desc structs (filled in) */
 197 static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
 198 {
 199         collect_symbols(a.eval(), v);   // eval() to expand assigned symbols
 200         collect_symbols(b.eval(), v);
 201         sym_desc_vec::iterator it = v.begin(), itend = v.end();
 202         while (it != itend) {
 203                 int deg_a = a.degree(*(it->sym));
 204                 int deg_b = b.degree(*(it->sym));
 205                 it->deg_a = deg_a;
 206                 it->deg_b = deg_b;
 207                 it->max_deg = std::max(deg_a, deg_b);
 208                 it->max_lcnops = std::max(a.lcoeff(*(it->sym)).nops(), b.lcoeff(*(it->sym)).nops());
 209                 it->ldeg_a = a.ldegree(*(it->sym));
 210                 it->ldeg_b = b.ldegree(*(it->sym));
 211                 it++;
 212         }
 213         sort(v.begin(), v.end());
 214 #if 0
 215         std::clog << "Symbols:\n";
 216         it = v.begin(); itend = v.end();
 217         while (it != itend) {
 218                 std::clog << " " << *it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << endl;
 219                 std::clog << "  lcoeff_a=" << a.lcoeff(*(it->sym)) << ", lcoeff_b=" << b.lcoeff(*(it->sym)) << endl;
 220                 it++;
 221         }
 222 #endif
 223 }
 224
 225
 226 /*
 227  *  Computation of LCM of denominators of coefficients of a polynomial
 228  */
 229
 230 // Compute LCM of denominators of coefficients by going through the
 231 // expression recursively (used internally by lcm_of_coefficients_denominators())
 232 static numeric lcmcoeff(const ex &e, const numeric &l)
 233 {
 234         if (e.info(info_flags::rational))
 235                 return lcm(ex_to_numeric(e).denom(), l);
 236         else if (is_ex_exactly_of_type(e, add)) {
 237                 numeric c = _num1();
 238                 for (unsigned i=0; i<e.nops(); i++)
 239                         c = lcmcoeff(e.op(i), c);
 240                 return lcm(c, l);
 241         } else if (is_ex_exactly_of_type(e, mul)) {
 242                 numeric c = _num1();
 243                 for (unsigned i=0; i<e.nops(); i++)
 244                         c *= lcmcoeff(e.op(i), _num1());
 245                 return lcm(c, l);
 246         } else if (is_ex_exactly_of_type(e, power))
 247                 return pow(lcmcoeff(e.op(0), l), ex_to_numeric(e.op(1)));
 248         return l;
 249 }
 250
 251 /** Compute LCM of denominators of coefficients of a polynomial.
 252  *  Given a polynomial with rational coefficients, this function computes
 253  *  the LCM of the denominators of all coefficients. This can be used
 254  *  to bring a polynomial from Q[X] to Z[X].
 255  *
 256  *  @param e  multivariate polynomial (need not be expanded)
 257  *  @return LCM of denominators of coefficients */
 258 static numeric lcm_of_coefficients_denominators(const ex &e)
 259 {
 260         return lcmcoeff(e, _num1());
 261 }
 262
 263 /** Bring polynomial from Q[X] to Z[X] by multiplying in the previously
 264  *  determined LCM of the coefficient's denominators.
 265  *
 266  *  @param e  multivariate polynomial (need not be expanded)
 267  *  @param lcm  LCM to multiply in */
 268 static ex multiply_lcm(const ex &e, const numeric &lcm)
 269 {
 270         if (is_ex_exactly_of_type(e, mul)) {
 271                 ex c = _ex1();
 272                 numeric lcm_accum = _num1();
 273                 for (unsigned i=0; i<e.nops(); i++) {
 274                         numeric op_lcm = lcmcoeff(e.op(i), _num1());
 275                         c *= multiply_lcm(e.op(i), op_lcm);
 276                         lcm_accum *= op_lcm;
 277                 }
 278                 c *= lcm / lcm_accum;
 279                 return c;
 280         } else if (is_ex_exactly_of_type(e, add)) {
 281                 ex c = _ex0();
 282                 for (unsigned i=0; i<e.nops(); i++)
 283                         c += multiply_lcm(e.op(i), lcm);
 284                 return c;
 285         } else if (is_ex_exactly_of_type(e, power)) {
 286                 return pow(multiply_lcm(e.op(0), lcm.power(ex_to_numeric(e.op(1)).inverse())), e.op(1));
 287         } else
 288                 return e * lcm;
 289 }
 290
 291
 292 /** Compute the integer content (= GCD of all numeric coefficients) of an
 293  *  expanded polynomial.
 294  *
 295  *  @param e  expanded polynomial
 296  *  @return integer content */
 297 numeric ex::integer_content(void) const
 298 {
 299         GINAC_ASSERT(bp!=0);
 300         return bp->integer_content();
 301 }
 302
 303 numeric basic::integer_content(void) const
 304 {
 305         return _num1();
 306 }
 307
 308 numeric numeric::integer_content(void) const
 309 {
 310         return abs(*this);
 311 }
 312
 313 numeric add::integer_content(void) const
 314 {
 315         epvector::const_iterator it = seq.begin();
 316         epvector::const_iterator itend = seq.end();
 317         numeric c = _num0();
 318         while (it != itend) {
 319                 GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric));
 320                 GINAC_ASSERT(is_ex_exactly_of_type(it->coeff,numeric));
 321                 c = gcd(ex_to_numeric(it->coeff), c);
 322                 it++;
 323         }
 324         GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
 325         c = gcd(ex_to_numeric(overall_coeff),c);
 326         return c;
 327 }
 328
 329 numeric mul::integer_content(void) const
 330 {
 331 #ifdef DO_GINAC_ASSERT
 332         epvector::const_iterator it = seq.begin();
 333         epvector::const_iterator itend = seq.end();
 334         while (it != itend) {
 335                 GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric));
 336                 ++it;
 337         }
 338 #endif // def DO_GINAC_ASSERT
 339         GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
 340         return abs(ex_to_numeric(overall_coeff));
 341 }
 342
 343
 344 /*
 345  *  Polynomial quotients and remainders
 346  */
 347
 348 /** Quotient q(x) of polynomials a(x) and b(x) in Q[x].
 349  *  It satisfies a(x)=b(x)*q(x)+r(x).
 350  *
 351  *  @param a  first polynomial in x (dividend)
 352  *  @param b  second polynomial in x (divisor)
 353  *  @param x  a and b are polynomials in x
 354  *  @param check_args  check whether a and b are polynomials with rational
 355  *         coefficients (defaults to "true")
 356  *  @return quotient of a and b in Q[x] */
 357 ex quo(const ex &a, const ex &b, const symbol &x, bool check_args)
 358 {
 359         if (b.is_zero())
 360                 throw(std::overflow_error("quo: division by zero"));
 361         if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric))
 362                 return a / b;
 363 #if FAST_COMPARE
 364         if (a.is_equal(b))
 365                 return _ex1();
 366 #endif
 367         if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)))
 368                 throw(std::invalid_argument("quo: arguments must be polynomials over the rationals"));
 369
 370         // Polynomial long division
 371         ex q = _ex0();
 372         ex r = a.expand();
 373         if (r.is_zero())
 374                 return r;
 375         int bdeg = b.degree(x);
 376         int rdeg = r.degree(x);
 377         ex blcoeff = b.expand().coeff(x, bdeg);
 378         bool blcoeff_is_numeric = is_ex_exactly_of_type(blcoeff, numeric);
 379         while (rdeg >= bdeg) {
 380                 ex term, rcoeff = r.coeff(x, rdeg);
 381                 if (blcoeff_is_numeric)
 382                         term = rcoeff / blcoeff;
 383                 else {
 384                         if (!divide(rcoeff, blcoeff, term, false))
 385                                 return *new ex(fail());
 386                 }
 387                 term *= power(x, rdeg - bdeg);
 388                 q += term;
 389                 r -= (term * b).expand();
 390                 if (r.is_zero())
 391                         break;
 392                 rdeg = r.degree(x);
 393         }
 394         return q;
 395 }
 396
 397
 398 /** Remainder r(x) of polynomials a(x) and b(x) in Q[x].
 399  *  It satisfies a(x)=b(x)*q(x)+r(x).
 400  *
 401  *  @param a  first polynomial in x (dividend)
 402  *  @param b  second polynomial in x (divisor)
 403  *  @param x  a and b are polynomials in x
 404  *  @param check_args  check whether a and b are polynomials with rational
 405  *         coefficients (defaults to "true")
 406  *  @return remainder of a(x) and b(x) in Q[x] */
 407 ex rem(const ex &a, const ex &b, const symbol &x, bool check_args)
 408 {
 409         if (b.is_zero())
 410                 throw(std::overflow_error("rem: division by zero"));
 411         if (is_ex_exactly_of_type(a, numeric)) {
 412                 if  (is_ex_exactly_of_type(b, numeric))
 413                         return _ex0();
 414                 else
 415                         return b;
 416         }
 417 #if FAST_COMPARE
 418         if (a.is_equal(b))
 419                 return _ex0();
 420 #endif
 421         if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)))
 422                 throw(std::invalid_argument("rem: arguments must be polynomials over the rationals"));
 423
 424         // Polynomial long division
 425         ex r = a.expand();
 426         if (r.is_zero())
 427                 return r;
 428         int bdeg = b.degree(x);
 429         int rdeg = r.degree(x);
 430         ex blcoeff = b.expand().coeff(x, bdeg);
 431         bool blcoeff_is_numeric = is_ex_exactly_of_type(blcoeff, numeric);
 432         while (rdeg >= bdeg) {
 433                 ex term, rcoeff = r.coeff(x, rdeg);
 434                 if (blcoeff_is_numeric)
 435                         term = rcoeff / blcoeff;
 436                 else {
 437                         if (!divide(rcoeff, blcoeff, term, false))
 438                                 return *new ex(fail());
 439                 }
 440                 term *= power(x, rdeg - bdeg);
 441                 r -= (term * b).expand();
 442                 if (r.is_zero())
 443                         break;
 444                 rdeg = r.degree(x);
 445         }
 446         return r;
 447 }
 448
 449
 450 /** Pseudo-remainder of polynomials a(x) and b(x) in Z[x].
 451  *
 452  *  @param a  first polynomial in x (dividend)
 453  *  @param b  second polynomial in x (divisor)
 454  *  @param x  a and b are polynomials in x
 455  *  @param check_args  check whether a and b are polynomials with rational
 456  *         coefficients (defaults to "true")
 457  *  @return pseudo-remainder of a(x) and b(x) in Z[x] */
 458 ex prem(const ex &a, const ex &b, const symbol &x, bool check_args)
 459 {
 460         if (b.is_zero())
 461                 throw(std::overflow_error("prem: division by zero"));
 462         if (is_ex_exactly_of_type(a, numeric)) {
 463                 if (is_ex_exactly_of_type(b, numeric))
 464                         return _ex0();
 465                 else
 466                         return b;
 467         }
 468         if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)))
 469                 throw(std::invalid_argument("prem: arguments must be polynomials over the rationals"));
 470
 471         // Polynomial long division
 472         ex r = a.expand();
 473         ex eb = b.expand();
 474         int rdeg = r.degree(x);
 475         int bdeg = eb.degree(x);
 476         ex blcoeff;
 477         if (bdeg <= rdeg) {
 478                 blcoeff = eb.coeff(x, bdeg);
 479                 if (bdeg == 0)
 480                         eb = _ex0();
 481                 else
 482                         eb -= blcoeff * power(x, bdeg);
 483         } else
 484                 blcoeff = _ex1();
 485
 486         int delta = rdeg - bdeg + 1, i = 0;
 487         while (rdeg >= bdeg && !r.is_zero()) {
 488                 ex rlcoeff = r.coeff(x, rdeg);
 489                 ex term = (power(x, rdeg - bdeg) * eb * rlcoeff).expand();
 490                 if (rdeg == 0)
 491                         r = _ex0();
 492                 else
 493                         r -= rlcoeff * power(x, rdeg);
 494                 r = (blcoeff * r).expand() - term;
 495                 rdeg = r.degree(x);
 496                 i++;
 497         }
 498         return power(blcoeff, delta - i) * r;
 499 }
 500
 501
 502 /** Sparse pseudo-remainder of polynomials a(x) and b(x) in Z[x].
 503  *
 504  *  @param a  first polynomial in x (dividend)
 505  *  @param b  second polynomial in x (divisor)
 506  *  @param x  a and b are polynomials in x
 507  *  @param check_args  check whether a and b are polynomials with rational
 508  *         coefficients (defaults to "true")
 509  *  @return sparse pseudo-remainder of a(x) and b(x) in Z[x] */
 510
 511 ex sprem(const ex &a, const ex &b, const symbol &x, bool check_args)
 512 {
 513         if (b.is_zero())
 514                 throw(std::overflow_error("prem: division by zero"));
 515         if (is_ex_exactly_of_type(a, numeric)) {
 516                 if (is_ex_exactly_of_type(b, numeric))
 517                         return _ex0();
 518                 else
 519                         return b;
 520         }
 521         if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)))
 522                 throw(std::invalid_argument("prem: arguments must be polynomials over the rationals"));
 523
 524         // Polynomial long division
 525         ex r = a.expand();
 526         ex eb = b.expand();
 527         int rdeg = r.degree(x);
 528         int bdeg = eb.degree(x);
 529         ex blcoeff;
 530         if (bdeg <= rdeg) {
 531                 blcoeff = eb.coeff(x, bdeg);
 532                 if (bdeg == 0)
 533                         eb = _ex0();
 534                 else
 535                         eb -= blcoeff * power(x, bdeg);
 536         } else
 537                 blcoeff = _ex1();
 538
 539         while (rdeg >= bdeg && !r.is_zero()) {
 540                 ex rlcoeff = r.coeff(x, rdeg);
 541                 ex term = (power(x, rdeg - bdeg) * eb * rlcoeff).expand();
 542                 if (rdeg == 0)
 543                         r = _ex0();
 544                 else
 545                         r -= rlcoeff * power(x, rdeg);
 546                 r = (blcoeff * r).expand() - term;
 547                 rdeg = r.degree(x);
 548         }
 549         return r;
 550 }
 551
 552
 553 /** Exact polynomial division of a(X) by b(X) in Q[X].
 554  *
 555  *  @param a  first multivariate polynomial (dividend)
 556  *  @param b  second multivariate polynomial (divisor)
 557  *  @param q  quotient (returned)
 558  *  @param check_args  check whether a and b are polynomials with rational
 559  *         coefficients (defaults to "true")
 560  *  @return "true" when exact division succeeds (quotient returned in q),
 561  *          "false" otherwise */
 562 bool divide(const ex &a, const ex &b, ex &q, bool check_args)
 563 {
 564         q = _ex0();
 565         if (b.is_zero())
 566                 throw(std::overflow_error("divide: division by zero"));
 567         if (a.is_zero())
 568                 return true;
 569         if (is_ex_exactly_of_type(b, numeric)) {
 570                 q = a / b;
 571                 return true;
 572         } else if (is_ex_exactly_of_type(a, numeric))
 573                 return false;
 574 #if FAST_COMPARE
 575         if (a.is_equal(b)) {
 576                 q = _ex1();
 577                 return true;
 578         }
 579 #endif
 580         if (check_args && (!a.info(info_flags::rational_polynomial) ||
 581                            !b.info(info_flags::rational_polynomial)))
 582                 throw(std::invalid_argument("divide: arguments must be polynomials over the rationals"));
 583
 584         // Find first symbol
 585         const symbol *x;
 586         if (!get_first_symbol(a, x) && !get_first_symbol(b, x))
 587                 throw(std::invalid_argument("invalid expression in divide()"));
 588
 589         // Polynomial long division (recursive)
 590         ex r = a.expand();
 591         if (r.is_zero())
 592                 return true;
 593         int bdeg = b.degree(*x);
 594         int rdeg = r.degree(*x);
 595         ex blcoeff = b.expand().coeff(*x, bdeg);
 596         bool blcoeff_is_numeric = is_ex_exactly_of_type(blcoeff, numeric);
 597         while (rdeg >= bdeg) {
 598                 ex term, rcoeff = r.coeff(*x, rdeg);
 599                 if (blcoeff_is_numeric)
 600                         term = rcoeff / blcoeff;
 601                 else
 602                         if (!divide(rcoeff, blcoeff, term, false))
 603                                 return false;
 604                 term *= power(*x, rdeg - bdeg);
 605                 q += term;
 606                 r -= (term * b).expand();
 607                 if (r.is_zero())
 608                         return true;
 609                 rdeg = r.degree(*x);
 610         }
 611         return false;
 612 }
 613
 614
 615 #if USE_REMEMBER
 616 /*
 617  *  Remembering
 618  */
 619
 620 typedef std::pair<ex, ex> ex2;
 621 typedef std::pair<ex, bool> exbool;
 622
 623 struct ex2_less {
 624         bool operator() (const ex2 p, const ex2 q) const
 625         {
 626                 return p.first.compare(q.first) < 0 || (!(q.first.compare(p.first) < 0) && p.second.compare(q.second) < 0);
 627         }
 628 };
 629
 630 typedef std::map<ex2, exbool, ex2_less> ex2_exbool_remember;
 631 #endif
 632
 633
 634 /** Exact polynomial division of a(X) by b(X) in Z[X].
 635  *  This functions works like divide() but the input and output polynomials are
 636  *  in Z[X] instead of Q[X] (i.e. they have integer coefficients). Unlike
 637  *  divide(), it doesn´t check whether the input polynomials really are integer
 638  *  polynomials, so be careful of what you pass in. Also, you have to run
 639  *  get_symbol_stats() over the input polynomials before calling this function
 640  *  and pass an iterator to the first element of the sym_desc vector. This
 641  *  function is used internally by the heur_gcd().
 642  *
 643  *  @param a  first multivariate polynomial (dividend)
 644  *  @param b  second multivariate polynomial (divisor)
 645  *  @param q  quotient (returned)
 646  *  @param var  iterator to first element of vector of sym_desc structs
 647  *  @return "true" when exact division succeeds (the quotient is returned in
 648  *          q), "false" otherwise.
 649  *  @see get_symbol_stats, heur_gcd */
 650 static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_iterator var)
 651 {
 652         q = _ex0();
 653         if (b.is_zero())
 654                 throw(std::overflow_error("divide_in_z: division by zero"));
 655         if (b.is_equal(_ex1())) {
 656                 q = a;
 657                 return true;
 658         }
 659         if (is_ex_exactly_of_type(a, numeric)) {
 660                 if (is_ex_exactly_of_type(b, numeric)) {
 661                         q = a / b;
 662                         return q.info(info_flags::integer);
 663                 } else
 664                         return false;
 665         }
 666 #if FAST_COMPARE
 667         if (a.is_equal(b)) {
 668                 q = _ex1();
 669                 return true;
 670         }
 671 #endif
 672
 673 #if USE_REMEMBER
 674         // Remembering
 675         static ex2_exbool_remember dr_remember;
 676         ex2_exbool_remember::const_iterator remembered = dr_remember.find(ex2(a, b));
 677         if (remembered != dr_remember.end()) {
 678                 q = remembered->second.first;
 679                 return remembered->second.second;
 680         }
 681 #endif
 682
 683         // Main symbol
 684         const symbol *x = var->sym;
 685
 686         // Compare degrees
 687         int adeg = a.degree(*x), bdeg = b.degree(*x);
 688         if (bdeg > adeg)
 689                 return false;
 690
 691 #if USE_TRIAL_DIVISION
 692
 693         // Trial division with polynomial interpolation
 694         int i, k;
 695
 696         // Compute values at evaluation points 0..adeg
 697         vector<numeric> alpha; alpha.reserve(adeg + 1);
 698         exvector u; u.reserve(adeg + 1);
 699         numeric point = _num0();
 700         ex c;
 701         for (i=0; i<=adeg; i++) {
 702                 ex bs = b.subs(*x == point);
 703                 while (bs.is_zero()) {
 704                         point += _num1();
 705                         bs = b.subs(*x == point);
 706                 }
 707                 if (!divide_in_z(a.subs(*x == point), bs, c, var+1))
 708                         return false;
 709                 alpha.push_back(point);
 710                 u.push_back(c);
 711                 point += _num1();
 712         }
 713
 714         // Compute inverses
 715         vector<numeric> rcp; rcp.reserve(adeg + 1);
 716         rcp.push_back(_num0());
 717         for (k=1; k<=adeg; k++) {
 718                 numeric product = alpha[k] - alpha[0];
 719                 for (i=1; i<k; i++)
 720                         product *= alpha[k] - alpha[i];
 721                 rcp.push_back(product.inverse());
 722         }
 723
 724         // Compute Newton coefficients
 725         exvector v; v.reserve(adeg + 1);
 726         v.push_back(u[0]);
 727         for (k=1; k<=adeg; k++) {
 728                 ex temp = v[k - 1];
 729                 for (i=k-2; i>=0; i--)
 730                         temp = temp * (alpha[k] - alpha[i]) + v[i];
 731                 v.push_back((u[k] - temp) * rcp[k]);
 732         }
 733
 734         // Convert from Newton form to standard form
 735         c = v[adeg];
 736         for (k=adeg-1; k>=0; k--)
 737                 c = c * (*x - alpha[k]) + v[k];
 738
 739         if (c.degree(*x) == (adeg - bdeg)) {
 740                 q = c.expand();
 741                 return true;
 742         } else
 743                 return false;
 744
 745 #else
 746
 747         // Polynomial long division (recursive)
 748         ex r = a.expand();
 749         if (r.is_zero())
 750                 return true;
 751         int rdeg = adeg;
 752         ex eb = b.expand();
 753         ex blcoeff = eb.coeff(*x, bdeg);
 754         while (rdeg >= bdeg) {
 755                 ex term, rcoeff = r.coeff(*x, rdeg);
 756                 if (!divide_in_z(rcoeff, blcoeff, term, var+1))
 757                         break;
 758                 term = (term * power(*x, rdeg - bdeg)).expand();
 759                 q += term;
 760                 r -= (term * eb).expand();
 761                 if (r.is_zero()) {
 762 #if USE_REMEMBER
 763                         dr_remember[ex2(a, b)] = exbool(q, true);
 764 #endif
 765                         return true;
 766                 }
 767                 rdeg = r.degree(*x);
 768         }
 769 #if USE_REMEMBER
 770         dr_remember[ex2(a, b)] = exbool(q, false);
 771 #endif
 772         return false;
 773
 774 #endif
 775 }
 776
 777
 778 /*
 779  *  Separation of unit part, content part and primitive part of polynomials
 780  */
 781
 782 /** Compute unit part (= sign of leading coefficient) of a multivariate
 783  *  polynomial in Z[x]. The product of unit part, content part, and primitive
 784  *  part is the polynomial itself.
 785  *
 786  *  @param x  variable in which to compute the unit part
 787  *  @return unit part
 788  *  @see ex::content, ex::primpart */
 789 ex ex::unit(const symbol &x) const
 790 {
 791         ex c = expand().lcoeff(x);
 792         if (is_ex_exactly_of_type(c, numeric))
 793                 return c < _ex0() ? _ex_1() : _ex1();
 794         else {
 795                 const symbol *y;
 796                 if (get_first_symbol(c, y))
 797                         return c.unit(*y);
 798                 else
 799                         throw(std::invalid_argument("invalid expression in unit()"));
 800         }
 801 }
 802
 803
 804 /** Compute content part (= unit normal GCD of all coefficients) of a
 805  *  multivariate polynomial in Z[x].  The product of unit part, content part,
 806  *  and primitive part is the polynomial itself.
 807  *
 808  *  @param x  variable in which to compute the content part
 809  *  @return content part
 810  *  @see ex::unit, ex::primpart */
 811 ex ex::content(const symbol &x) const
 812 {
 813         if (is_zero())
 814                 return _ex0();
 815         if (is_ex_exactly_of_type(*this, numeric))
 816                 return info(info_flags::negative) ? -*this : *this;
 817         ex e = expand();
 818         if (e.is_zero())
 819                 return _ex0();
 820
 821         // First, try the integer content
 822         ex c = e.integer_content();
 823         ex r = e / c;
 824         ex lcoeff = r.lcoeff(x);
 825         if (lcoeff.info(info_flags::integer))
 826                 return c;
 827
 828         // GCD of all coefficients
 829         int deg = e.degree(x);
 830         int ldeg = e.ldegree(x);
 831         if (deg == ldeg)
 832                 return e.lcoeff(x) / e.unit(x);
 833         c = _ex0();
 834         for (int i=ldeg; i<=deg; i++)
 835                 c = gcd(e.coeff(x, i), c, NULL, NULL, false);
 836         return c;
 837 }
 838
 839
 840 /** Compute primitive part of a multivariate polynomial in Z[x].
 841  *  The product of unit part, content part, and primitive part is the
 842  *  polynomial itself.
 843  *
 844  *  @param x  variable in which to compute the primitive part
 845  *  @return primitive part
 846  *  @see ex::unit, ex::content */
 847 ex ex::primpart(const symbol &x) const
 848 {
 849         if (is_zero())
 850                 return _ex0();
 851         if (is_ex_exactly_of_type(*this, numeric))
 852                 return _ex1();
 853
 854         ex c = content(x);
 855         if (c.is_zero())
 856                 return _ex0();
 857         ex u = unit(x);
 858         if (is_ex_exactly_of_type(c, numeric))
 859                 return *this / (c * u);
 860         else
 861                 return quo(*this, c * u, x, false);
 862 }
 863
 864
 865 /** Compute primitive part of a multivariate polynomial in Z[x] when the
 866  *  content part is already known. This function is faster in computing the
 867  *  primitive part than the previous function.
 868  *
 869  *  @param x  variable in which to compute the primitive part
 870  *  @param c  previously computed content part
 871  *  @return primitive part */
 872 ex ex::primpart(const symbol &x, const ex &c) const
 873 {
 874         if (is_zero())
 875                 return _ex0();
 876         if (c.is_zero())
 877                 return _ex0();
 878         if (is_ex_exactly_of_type(*this, numeric))
 879                 return _ex1();
 880
 881         ex u = unit(x);
 882         if (is_ex_exactly_of_type(c, numeric))
 883                 return *this / (c * u);
 884         else
 885                 return quo(*this, c * u, x, false);
 886 }
 887
 888
 889 /*
 890  *  GCD of multivariate polynomials
 891  */
 892
 893 /** Compute GCD of polynomials in Q[X] using the Euclidean algorithm (not
 894  *  really suited for multivariate GCDs). This function is only provided for
 895  *  testing purposes.
 896  *
 897  *  @param a  first multivariate polynomial
 898  *  @param b  second multivariate polynomial
 899  *  @param x  pointer to symbol (main variable) in which to compute the GCD in
 900  *  @return the GCD as a new expression
 901  *  @see gcd */
 902
 903 static ex eu_gcd(const ex &a, const ex &b, const symbol *x)
 904 {
 905 //std::clog << "eu_gcd(" << a << "," << b << ")\n";
 906
 907         // Sort c and d so that c has higher degree
 908         ex c, d;
 909         int adeg = a.degree(*x), bdeg = b.degree(*x);
 910         if (adeg >= bdeg) {
 911                 c = a;
 912                 d = b;
 913         } else {
 914                 c = b;
 915                 d = a;
 916         }
 917
 918         // Normalize in Q[x]
 919         c = c / c.lcoeff(*x);
 920         d = d / d.lcoeff(*x);
 921
 922         // Euclidean algorithm
 923         ex r;
 924         for (;;) {
 925 //std::clog << " d = " << d << endl;
 926                 r = rem(c, d, *x, false);
 927                 if (r.is_zero())
 928                         return d / d.lcoeff(*x);
 929                 c = d;
 930                 d = r;
 931         }
 932 }
 933
 934
 935 /** Compute GCD of multivariate polynomials using the Euclidean PRS algorithm
 936  *  with pseudo-remainders ("World's Worst GCD Algorithm", staying in Z[X]).
 937  *  This function is only provided for testing purposes.
 938  *
 939  *  @param a  first multivariate polynomial
 940  *  @param b  second multivariate polynomial
 941  *  @param x  pointer to symbol (main variable) in which to compute the GCD in
 942  *  @return the GCD as a new expression
 943  *  @see gcd */
 944
 945 static ex euprem_gcd(const ex &a, const ex &b, const symbol *x)
 946 {
 947 //std::clog << "euprem_gcd(" << a << "," << b << ")\n";
 948
 949         // Sort c and d so that c has higher degree
 950         ex c, d;
 951         int adeg = a.degree(*x), bdeg = b.degree(*x);
 952         if (adeg >= bdeg) {
 953                 c = a;
 954                 d = b;
 955         } else {
 956                 c = b;
 957                 d = a;
 958         }
 959
 960         // Calculate GCD of contents
 961         ex gamma = gcd(c.content(*x), d.content(*x), NULL, NULL, false);
 962
 963         // Euclidean algorithm with pseudo-remainders
 964         ex r;
 965         for (;;) {
 966 //std::clog << " d = " << d << endl;
 967                 r = prem(c, d, *x, false);
 968                 if (r.is_zero())
 969                         return d.primpart(*x) * gamma;
 970                 c = d;
 971                 d = r;
 972         }
 973 }
 974
 975
 976 /** Compute GCD of multivariate polynomials using the primitive Euclidean
 977  *  PRS algorithm (complete content removal at each step). This function is
 978  *  only provided for testing purposes.
 979  *
 980  *  @param a  first multivariate polynomial
 981  *  @param b  second multivariate polynomial
 982  *  @param x  pointer to symbol (main variable) in which to compute the GCD in
 983  *  @return the GCD as a new expression
 984  *  @see gcd */
 985
 986 static ex peu_gcd(const ex &a, const ex &b, const symbol *x)
 987 {
 988 //std::clog << "peu_gcd(" << a << "," << b << ")\n";
 989
 990         // Sort c and d so that c has higher degree
 991         ex c, d;
 992         int adeg = a.degree(*x), bdeg = b.degree(*x);
 993         int ddeg;
 994         if (adeg >= bdeg) {
 995                 c = a;
 996                 d = b;
 997                 ddeg = bdeg;
 998         } else {
 999                 c = b;
1000                 d = a;
1001                 ddeg = adeg;
1002         }
1003
1004         // Remove content from c and d, to be attached to GCD later
1005         ex cont_c = c.content(*x);
1006         ex cont_d = d.content(*x);
1007         ex gamma = gcd(cont_c, cont_d, NULL, NULL, false);
1008         if (ddeg == 0)
1009                 return gamma;
1010         c = c.primpart(*x, cont_c);
1011         d = d.primpart(*x, cont_d);
1012
1013         // Euclidean algorithm with content removal
1014         ex r;
1015         for (;;) {
1016 //std::clog << " d = " << d << endl;
1017                 r = prem(c, d, *x, false);
1018                 if (r.is_zero())
1019                         return gamma * d;
1020                 c = d;
1021                 d = r.primpart(*x);
1022         }
1023 }
1024
1025
1026 /** Compute GCD of multivariate polynomials using the reduced PRS algorithm.
1027  *  This function is only provided for testing purposes.
1028  *
1029  *  @param a  first multivariate polynomial
1030  *  @param b  second multivariate polynomial
1031  *  @param x  pointer to symbol (main variable) in which to compute the GCD in
1032  *  @return the GCD as a new expression
1033  *  @see gcd */
1034
1035 static ex red_gcd(const ex &a, const ex &b, const symbol *x)
1036 {
1037 //std::clog << "red_gcd(" << a << "," << b << ")\n";
1038
1039         // Sort c and d so that c has higher degree
1040         ex c, d;
1041         int adeg = a.degree(*x), bdeg = b.degree(*x);
1042         int cdeg, ddeg;
1043         if (adeg >= bdeg) {
1044                 c = a;
1045                 d = b;
1046                 cdeg = adeg;
1047                 ddeg = bdeg;
1048         } else {
1049                 c = b;
1050                 d = a;
1051                 cdeg = bdeg;
1052                 ddeg = adeg;
1053         }
1054
1055         // Remove content from c and d, to be attached to GCD later
1056         ex cont_c = c.content(*x);
1057         ex cont_d = d.content(*x);
1058         ex gamma = gcd(cont_c, cont_d, NULL, NULL, false);
1059         if (ddeg == 0)
1060                 return gamma;
1061         c = c.primpart(*x, cont_c);
1062         d = d.primpart(*x, cont_d);
1063
1064         // First element of divisor sequence
1065         ex r, ri = _ex1();
1066         int delta = cdeg - ddeg;
1067
1068         for (;;) {
1069                 // Calculate polynomial pseudo-remainder
1070 //std::clog << " d = " << d << endl;
1071                 r = prem(c, d, *x, false);
1072                 if (r.is_zero())
1073                         return gamma * d.primpart(*x);
1074                 c = d;
1075                 cdeg = ddeg;
1076
1077                 if (!divide(r, pow(ri, delta), d, false))
1078                         throw(std::runtime_error("invalid expression in red_gcd(), division failed"));
1079                 ddeg = d.degree(*x);
1080                 if (ddeg == 0) {
1081                         if (is_ex_exactly_of_type(r, numeric))
1082                                 return gamma;
1083                         else
1084                                 return gamma * r.primpart(*x);
1085                 }
1086
1087                 ri = c.expand().lcoeff(*x);
1088                 delta = cdeg - ddeg;
1089         }
1090 }
1091
1092
1093 /** Compute GCD of multivariate polynomials using the subresultant PRS
1094  *  algorithm. This function is used internally by gcd().
1095  *
1096  *  @param a   first multivariate polynomial
1097  *  @param b   second multivariate polynomial
1098  *  @param var iterator to first element of vector of sym_desc structs
1099  *  @return the GCD as a new expression
1100  *  @see gcd */
1101
1102 static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
1103 {
1104 //std::clog << "sr_gcd(" << a << "," << b << ")\n";
1105 #if STATISTICS
1106         sr_gcd_called++;
1107 #endif
1108
1109         // The first symbol is our main variable
1110         const symbol &x = *(var->sym);
1111
1112         // Sort c and d so that c has higher degree
1113         ex c, d;
1114         int adeg = a.degree(x), bdeg = b.degree(x);
1115         int cdeg, ddeg;
1116         if (adeg >= bdeg) {
1117                 c = a;
1118                 d = b;
1119                 cdeg = adeg;
1120                 ddeg = bdeg;
1121         } else {
1122                 c = b;
1123                 d = a;
1124                 cdeg = bdeg;
1125                 ddeg = adeg;
1126         }
1127
1128         // Remove content from c and d, to be attached to GCD later
1129         ex cont_c = c.content(x);
1130         ex cont_d = d.content(x);
1131         ex gamma = gcd(cont_c, cont_d, NULL, NULL, false);
1132         if (ddeg == 0)
1133                 return gamma;
1134         c = c.primpart(x, cont_c);
1135         d = d.primpart(x, cont_d);
1136 //std::clog << " content " << gamma << " removed, continuing with sr_gcd(" << c << "," << d << ")\n";
1137
1138         // First element of subresultant sequence
1139         ex r = _ex0(), ri = _ex1(), psi = _ex1();
1140         int delta = cdeg - ddeg;
1141
1142         for (;;) {
1143                 // Calculate polynomial pseudo-remainder
1144 //std::clog << " start of loop, psi = " << psi << ", calculating pseudo-remainder...\n";
1145 //std::clog << " d = " << d << endl;
1146                 r = prem(c, d, x, false);
1147                 if (r.is_zero())
1148                         return gamma * d.primpart(x);
1149                 c = d;
1150                 cdeg = ddeg;
1151 //std::clog << " dividing...\n";
1152                 if (!divide_in_z(r, ri * pow(psi, delta), d, var))
1153                         throw(std::runtime_error("invalid expression in sr_gcd(), division failed"));
1154                 ddeg = d.degree(x);
1155                 if (ddeg == 0) {
1156                         if (is_ex_exactly_of_type(r, numeric))
1157                                 return gamma;
1158                         else
1159                                 return gamma * r.primpart(x);
1160                 }
1161
1162                 // Next element of subresultant sequence
1163 //std::clog << " calculating next subresultant...\n";
1164                 ri = c.expand().lcoeff(x);
1165                 if (delta == 1)
1166                         psi = ri;
1167                 else if (delta)
1168                         divide_in_z(pow(ri, delta), pow(psi, delta-1), psi, var+1);
1169                 delta = cdeg - ddeg;
1170         }
1171 }
1172
1173
1174 /** Return maximum (absolute value) coefficient of a polynomial.
1175  *  This function is used internally by heur_gcd().
1176  *
1177  *  @param e  expanded multivariate polynomial
1178  *  @return maximum coefficient
1179  *  @see heur_gcd */
1180 numeric ex::max_coefficient(void) const
1181 {
1182         GINAC_ASSERT(bp!=0);
1183         return bp->max_coefficient();
1184 }
1185
1186 numeric basic::max_coefficient(void) const
1187 {
1188         return _num1();
1189 }
1190
1191 numeric numeric::max_coefficient(void) const
1192 {
1193         return abs(*this);
1194 }
1195
1196 numeric add::max_coefficient(void) const
1197 {
1198         epvector::const_iterator it = seq.begin();
1199         epvector::const_iterator itend = seq.end();
1200         GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
1201         numeric cur_max = abs(ex_to_numeric(overall_coeff));
1202         while (it != itend) {
1203                 numeric a;
1204                 GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric));
1205                 a = abs(ex_to_numeric(it->coeff));
1206                 if (a > cur_max)
1207                         cur_max = a;
1208                 it++;
1209         }
1210         return cur_max;
1211 }
1212
1213 numeric mul::max_coefficient(void) const
1214 {
1215 #ifdef DO_GINAC_ASSERT
1216         epvector::const_iterator it = seq.begin();
1217         epvector::const_iterator itend = seq.end();
1218         while (it != itend) {
1219                 GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric));
1220                 it++;
1221         }
1222 #endif // def DO_GINAC_ASSERT
1223         GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
1224         return abs(ex_to_numeric(overall_coeff));
1225 }
1226
1227
1228 /** Apply symmetric modular homomorphism to a multivariate polynomial.
1229  *  This function is used internally by heur_gcd().
1230  *
1231  *  @param e  expanded multivariate polynomial
1232  *  @param xi  modulus
1233  *  @return mapped polynomial
1234  *  @see heur_gcd */
1235 ex ex::smod(const numeric &xi) const
1236 {
1237         GINAC_ASSERT(bp!=0);
1238         return bp->smod(xi);
1239 }
1240
1241 ex basic::smod(const numeric &xi) const
1242 {
1243         return *this;
1244 }
1245
1246 ex numeric::smod(const numeric &xi) const
1247 {
1248 #ifndef NO_NAMESPACE_GINAC
1249         return GiNaC::smod(*this, xi);
1250 #else // ndef NO_NAMESPACE_GINAC
1251         return ::smod(*this, xi);
1252 #endif // ndef NO_NAMESPACE_GINAC
1253 }
1254
1255 ex add::smod(const numeric &xi) const
1256 {
1257         epvector newseq;
1258         newseq.reserve(seq.size()+1);
1259         epvector::const_iterator it = seq.begin();
1260         epvector::const_iterator itend = seq.end();
1261         while (it != itend) {
1262                 GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric));
1263 #ifndef NO_NAMESPACE_GINAC
1264                 numeric coeff = GiNaC::smod(ex_to_numeric(it->coeff), xi);
1265 #else // ndef NO_NAMESPACE_GINAC
1266                 numeric coeff = ::smod(ex_to_numeric(it->coeff), xi);
1267 #endif // ndef NO_NAMESPACE_GINAC
1268                 if (!coeff.is_zero())
1269                         newseq.push_back(expair(it->rest, coeff));
1270                 it++;
1271         }
1272         GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
1273 #ifndef NO_NAMESPACE_GINAC
1274         numeric coeff = GiNaC::smod(ex_to_numeric(overall_coeff), xi);
1275 #else // ndef NO_NAMESPACE_GINAC
1276         numeric coeff = ::smod(ex_to_numeric(overall_coeff), xi);
1277 #endif // ndef NO_NAMESPACE_GINAC
1278         return (new add(newseq,coeff))->setflag(status_flags::dynallocated);
1279 }
1280
1281 ex mul::smod(const numeric &xi) const
1282 {
1283 #ifdef DO_GINAC_ASSERT
1284         epvector::const_iterator it = seq.begin();
1285         epvector::const_iterator itend = seq.end();
1286         while (it != itend) {
1287                 GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric));
1288                 it++;
1289         }
1290 #endif // def DO_GINAC_ASSERT
1291         mul * mulcopyp=new mul(*this);
1292         GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
1293 #ifndef NO_NAMESPACE_GINAC
1294         mulcopyp->overall_coeff = GiNaC::smod(ex_to_numeric(overall_coeff),xi);
1295 #else // ndef NO_NAMESPACE_GINAC
1296         mulcopyp->overall_coeff = ::smod(ex_to_numeric(overall_coeff),xi);
1297 #endif // ndef NO_NAMESPACE_GINAC
1298         mulcopyp->clearflag(status_flags::evaluated);
1299         mulcopyp->clearflag(status_flags::hash_calculated);
1300         return mulcopyp->setflag(status_flags::dynallocated);
1301 }
1302
1303
1304 /** xi-adic polynomial interpolation */
1305 static ex interpolate(const ex &gamma, const numeric &xi, const symbol &x)
1306 {
1307         ex g = _ex0();
1308         ex e = gamma;
1309         numeric rxi = xi.inverse();
1310         for (int i=0; !e.is_zero(); i++) {
1311                 ex gi = e.smod(xi);
1312                 g += gi * power(x, i);
1313                 e = (e - gi) * rxi;
1314         }
1315         return g;
1316 }
1317
1318 /** Exception thrown by heur_gcd() to signal failure. */
1319 class gcdheu_failed {};
1320
1321 /** Compute GCD of multivariate polynomials using the heuristic GCD algorithm.
1322  *  get_symbol_stats() must have been called previously with the input
1323  *  polynomials and an iterator to the first element of the sym_desc vector
1324  *  passed in. This function is used internally by gcd().
1325  *
1326  *  @param a  first multivariate polynomial (expanded)
1327  *  @param b  second multivariate polynomial (expanded)
1328  *  @param ca  cofactor of polynomial a (returned), NULL to suppress
1329  *             calculation of cofactor
1330  *  @param cb  cofactor of polynomial b (returned), NULL to suppress
1331  *             calculation of cofactor
1332  *  @param var iterator to first element of vector of sym_desc structs
1333  *  @return the GCD as a new expression
1334  *  @see gcd
1335  *  @exception gcdheu_failed() */
1336 static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const_iterator var)
1337 {
1338 //std::clog << "heur_gcd(" << a << "," << b << ")\n";
1339 #if STATISTICS
1340         heur_gcd_called++;
1341 #endif
1342
1343         // Algorithms only works for non-vanishing input polynomials
1344         if (a.is_zero() || b.is_zero())
1345                 return *new ex(fail());
1346
1347         // GCD of two numeric values -> CLN
1348         if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) {
1349                 numeric g = gcd(ex_to_numeric(a), ex_to_numeric(b));
1350                 if (ca)
1351                         *ca = ex_to_numeric(a) / g;
1352                 if (cb)
1353                         *cb = ex_to_numeric(b) / g;
1354                 return g;
1355         }
1356
1357         // The first symbol is our main variable
1358         const symbol &x = *(var->sym);
1359
1360         // Remove integer content
1361         numeric gc = gcd(a.integer_content(), b.integer_content());
1362         numeric rgc = gc.inverse();
1363         ex p = a * rgc;
1364         ex q = b * rgc;
1365         int maxdeg =  std::max(p.degree(x),q.degree(x));
1366
1367         // Find evaluation point
1368         numeric mp = p.max_coefficient();
1369         numeric mq = q.max_coefficient();
1370         numeric xi;
1371         if (mp > mq)
1372                 xi = mq * _num2() + _num2();
1373         else
1374                 xi = mp * _num2() + _num2();
1375
1376         // 6 tries maximum
1377         for (int t=0; t<6; t++) {
1378                 if (xi.int_length() * maxdeg > 100000) {
1379 //std::clog << "giving up heur_gcd, xi.int_length = " << xi.int_length() << ", maxdeg = " << maxdeg << endl;
1380                         throw gcdheu_failed();
1381                 }
1382
1383                 // Apply evaluation homomorphism and calculate GCD
1384                 ex cp, cq;
1385                 ex gamma = heur_gcd(p.subs(x == xi), q.subs(x == xi), &cp, &cq, var+1).expand();
1386                 if (!is_ex_exactly_of_type(gamma, fail)) {
1387
1388                         // Reconstruct polynomial from GCD of mapped polynomials
1389                         ex g = interpolate(gamma, xi, x);
1390
1391                         // Remove integer content
1392                         g /= g.integer_content();
1393
1394                         // If the calculated polynomial divides both p and q, this is the GCD
1395                         ex dummy;
1396                         if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) {
1397                                 g *= gc;
1398                                 ex lc = g.lcoeff(x);
1399                                 if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative())
1400                                         return -g;
1401                                 else
1402                                         return g;
1403                         }
1404 #if 0
1405                         cp = interpolate(cp, xi, x);
1406                         if (divide_in_z(cp, p, g, var)) {
1407                                 if (divide_in_z(g, q, cb ? *cb : dummy, var)) {
1408                                         g *= gc;
1409                                         if (ca)
1410                                                 *ca = cp;
1411                                         ex lc = g.lcoeff(x);
1412                                         if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative())
1413                                                 return -g;
1414                                         else
1415                                                 return g;
1416                                 }
1417                         }
1418                         cq = interpolate(cq, xi, x);
1419                         if (divide_in_z(cq, q, g, var)) {
1420                                 if (divide_in_z(g, p, ca ? *ca : dummy, var)) {
1421                                         g *= gc;
1422                                         if (cb)
1423                                                 *cb = cq;
1424                                         ex lc = g.lcoeff(x);
1425                                         if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative())
1426                                                 return -g;
1427                                         else
1428                                                 return g;
1429                                 }
1430                         }
1431 #endif
1432                 }
1433
1434                 // Next evaluation point
1435                 xi = iquo(xi * isqrt(isqrt(xi)) * numeric(73794), numeric(27011));
1436         }
1437         return *new ex(fail());
1438 }
1439
1440
1441 /** Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X)
1442  *  and b(X) in Z[X].
1443  *
1444  *  @param a  first multivariate polynomial
1445  *  @param b  second multivariate polynomial
1446  *  @param check_args  check whether a and b are polynomials with rational
1447  *         coefficients (defaults to "true")
1448  *  @return the GCD as a new expression */
1449 ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args)
1450 {
1451 //std::clog << "gcd(" << a << "," << b << ")\n";
1452 #if STATISTICS
1453         gcd_called++;
1454 #endif
1455
1456         // GCD of numerics -> CLN
1457         if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) {
1458                 numeric g = gcd(ex_to_numeric(a), ex_to_numeric(b));
1459                 if (ca || cb) {
1460                         if (g.is_zero()) {
1461                                 if (ca)
1462                                         *ca = _ex0();
1463                                 if (cb)
1464                                         *cb = _ex0();
1465                         } else {
1466                                 if (ca)
1467                                         *ca = ex_to_numeric(a) / g;
1468                                 if (cb)
1469                                         *cb = ex_to_numeric(b) / g;
1470                         }
1471                 }
1472                 return g;
1473         }
1474
1475         // Check arguments
1476         if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) {
1477                 throw(std::invalid_argument("gcd: arguments must be polynomials over the rationals"));
1478         }
1479
1480         // Partially factored cases (to avoid expanding large expressions)
1481         if (is_ex_exactly_of_type(a, mul)) {
1482                 if (is_ex_exactly_of_type(b, mul) && b.nops() > a.nops())
1483                         goto factored_b;
1484 factored_a:
1485                 ex g = _ex1();
1486                 ex acc_ca = _ex1();
1487                 ex part_b = b;
1488                 for (unsigned i=0; i<a.nops(); i++) {
1489                         ex part_ca, part_cb;
1490                         g *= gcd(a.op(i), part_b, &part_ca, &part_cb, check_args);
1491                         acc_ca *= part_ca;
1492                         part_b = part_cb;
1493                 }
1494                 if (ca)
1495                         *ca = acc_ca;
1496                 if (cb)
1497                         *cb = part_b;
1498                 return g;
1499         } else if (is_ex_exactly_of_type(b, mul)) {
1500                 if (is_ex_exactly_of_type(a, mul) && a.nops() > b.nops())
1501                         goto factored_a;
1502 factored_b:
1503                 ex g = _ex1();
1504                 ex acc_cb = _ex1();
1505                 ex part_a = a;
1506                 for (unsigned i=0; i<b.nops(); i++) {
1507                         ex part_ca, part_cb;
1508                         g *= gcd(part_a, b.op(i), &part_ca, &part_cb, check_args);
1509                         acc_cb *= part_cb;
1510                         part_a = part_ca;
1511                 }
1512                 if (ca)
1513                         *ca = part_a;
1514                 if (cb)
1515                         *cb = acc_cb;
1516                 return g;
1517         }
1518
1519 #if FAST_COMPARE
1520         // Input polynomials of the form poly^n are sometimes also trivial
1521         if (is_ex_exactly_of_type(a, power)) {
1522                 ex p = a.op(0);
1523                 if (is_ex_exactly_of_type(b, power)) {
1524                         if (p.is_equal(b.op(0))) {
1525                                 // a = p^n, b = p^m, gcd = p^min(n, m)
1526                                 ex exp_a = a.op(1), exp_b = b.op(1);
1527                                 if (exp_a < exp_b) {
1528                                         if (ca)
1529                                                 *ca = _ex1();
1530                                         if (cb)
1531                                                 *cb = power(p, exp_b - exp_a);
1532                                         return power(p, exp_a);
1533                                 } else {
1534                                         if (ca)
1535                                                 *ca = power(p, exp_a - exp_b);
1536                                         if (cb)
1537                                                 *cb = _ex1();
1538                                         return power(p, exp_b);
1539                                 }
1540                         }
1541                 } else {
1542                         if (p.is_equal(b)) {
1543                                 // a = p^n, b = p, gcd = p
1544                                 if (ca)
1545                                         *ca = power(p, a.op(1) - 1);
1546                                 if (cb)
1547                                         *cb = _ex1();
1548                                 return p;
1549                         }
1550                 }
1551         } else if (is_ex_exactly_of_type(b, power)) {
1552                 ex p = b.op(0);
1553                 if (p.is_equal(a)) {
1554                         // a = p, b = p^n, gcd = p
1555                         if (ca)
1556                                 *ca = _ex1();
1557                         if (cb)
1558                                 *cb = power(p, b.op(1) - 1);
1559                         return p;
1560                 }
1561         }
1562 #endif
1563
1564         // Some trivial cases
1565         ex aex = a.expand(), bex = b.expand();
1566         if (aex.is_zero()) {
1567                 if (ca)
1568                         *ca = _ex0();
1569                 if (cb)
1570                         *cb = _ex1();
1571                 return b;
1572         }
1573         if (bex.is_zero()) {
1574                 if (ca)
1575                         *ca = _ex1();
1576                 if (cb)
1577                         *cb = _ex0();
1578                 return a;
1579         }
1580         if (aex.is_equal(_ex1()) || bex.is_equal(_ex1())) {
1581                 if (ca)
1582                         *ca = a;
1583                 if (cb)
1584                         *cb = b;
1585                 return _ex1();
1586         }
1587 #if FAST_COMPARE
1588         if (a.is_equal(b)) {
1589                 if (ca)
1590                         *ca = _ex1();
1591                 if (cb)
1592                         *cb = _ex1();
1593                 return a;
1594         }
1595 #endif
1596
1597         // Gather symbol statistics
1598         sym_desc_vec sym_stats;
1599         get_symbol_stats(a, b, sym_stats);
1600
1601         // The symbol with least degree is our main variable
1602         sym_desc_vec::const_iterator var = sym_stats.begin();
1603         const symbol &x = *(var->sym);
1604
1605         // Cancel trivial common factor
1606         int ldeg_a = var->ldeg_a;
1607         int ldeg_b = var->ldeg_b;
1608         int min_ldeg = std::min(ldeg_a,ldeg_b);
1609         if (min_ldeg > 0) {
1610                 ex common = power(x, min_ldeg);
1611 //std::clog << "trivial common factor " << common << endl;
1612                 return gcd((aex / common).expand(), (bex / common).expand(), ca, cb, false) * common;
1613         }
1614
1615         // Try to eliminate variables
1616         if (var->deg_a == 0) {
1617 //std::clog << "eliminating variable " << x << " from b" << endl;
1618                 ex c = bex.content(x);
1619                 ex g = gcd(aex, c, ca, cb, false);
1620                 if (cb)
1621                         *cb *= bex.unit(x) * bex.primpart(x, c);
1622                 return g;
1623         } else if (var->deg_b == 0) {
1624 //std::clog << "eliminating variable " << x << " from a" << endl;
1625                 ex c = aex.content(x);
1626                 ex g = gcd(c, bex, ca, cb, false);
1627                 if (ca)
1628                         *ca *= aex.unit(x) * aex.primpart(x, c);
1629                 return g;
1630         }
1631
1632         ex g;
1633 #if 1
1634         // Try heuristic algorithm first, fall back to PRS if that failed
1635         try {
1636                 g = heur_gcd(aex, bex, ca, cb, var);
1637         } catch (gcdheu_failed) {
1638                 g = *new ex(fail());
1639         }
1640         if (is_ex_exactly_of_type(g, fail)) {
1641 //std::clog << "heuristics failed" << endl;
1642 #if STATISTICS
1643                 heur_gcd_failed++;
1644 #endif
1645 #endif
1646 //              g = heur_gcd(aex, bex, ca, cb, var);
1647 //              g = eu_gcd(aex, bex, &x);
1648 //              g = euprem_gcd(aex, bex, &x);
1649 //              g = peu_gcd(aex, bex, &x);
1650 //              g = red_gcd(aex, bex, &x);
1651                 g = sr_gcd(aex, bex, var);
1652                 if (g.is_equal(_ex1())) {
1653                         // Keep cofactors factored if possible
1654                         if (ca)
1655                                 *ca = a;
1656                         if (cb)
1657                                 *cb = b;
1658                 } else {
1659                         if (ca)
1660                                 divide(aex, g, *ca, false);
1661                         if (cb)
1662                                 divide(bex, g, *cb, false);
1663                 }
1664 #if 1
1665         } else {
1666                 if (g.is_equal(_ex1())) {
1667                         // Keep cofactors factored if possible
1668                         if (ca)
1669                                 *ca = a;
1670                         if (cb)
1671                                 *cb = b;
1672                 }
1673         }
1674 #endif
1675         return g;
1676 }
1677
1678
1679 /** Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X].
1680  *
1681  *  @param a  first multivariate polynomial
1682  *  @param b  second multivariate polynomial
1683  *  @param check_args  check whether a and b are polynomials with rational
1684  *         coefficients (defaults to "true")
1685  *  @return the LCM as a new expression */
1686 ex lcm(const ex &a, const ex &b, bool check_args)
1687 {
1688         if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric))
1689                 return lcm(ex_to_numeric(a), ex_to_numeric(b));
1690         if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)))
1691                 throw(std::invalid_argument("lcm: arguments must be polynomials over the rationals"));
1692
1693         ex ca, cb;
1694         ex g = gcd(a, b, &ca, &cb, false);
1695         return ca * cb * g;
1696 }
1697
1698
1699 /*
1700  *  Square-free factorization
1701  */
1702
1703 // Univariate GCD of polynomials in Q[x] (used internally by sqrfree()).
1704 // a and b can be multivariate polynomials but they are treated as univariate polynomials in x.
1705 static ex univariate_gcd(const ex &a, const ex &b, const symbol &x)
1706 {
1707         if (a.is_zero())
1708                 return b;
1709         if (b.is_zero())
1710                 return a;
1711         if (a.is_equal(_ex1()) || b.is_equal(_ex1()))
1712                 return _ex1();
1713         if (is_ex_of_type(a, numeric) && is_ex_of_type(b, numeric))
1714                 return gcd(ex_to_numeric(a), ex_to_numeric(b));
1715         if (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))
1716                 throw(std::invalid_argument("univariate_gcd: arguments must be polynomials over the rationals"));
1717
1718         // Euclidean algorithm
1719         ex c, d, r;
1720         if (a.degree(x) >= b.degree(x)) {
1721                 c = a;
1722                 d = b;
1723         } else {
1724                 c = b;
1725                 d = a;
1726         }
1727         for (;;) {
1728                 r = rem(c, d, x, false);
1729                 if (r.is_zero())
1730                         break;
1731                 c = d;
1732                 d = r;
1733         }
1734         return d / d.lcoeff(x);
1735 }
1736
1737
1738 /** Compute square-free factorization of multivariate polynomial a(x) using
1739  *  Yun´s algorithm.
1740  *
1741  * @param a  multivariate polynomial
1742  * @param x  variable to factor in
1743  * @return factored polynomial */
1744 ex sqrfree(const ex &a, const symbol &x)
1745 {
1746         int i = 1;
1747         ex res = _ex1();
1748         ex b = a.diff(x);
1749         ex c = univariate_gcd(a, b, x);
1750         ex w;
1751         if (c.is_equal(_ex1())) {
1752                 w = a;
1753         } else {
1754                 w = quo(a, c, x);
1755                 ex y = quo(b, c, x);
1756                 ex z = y - w.diff(x);
1757                 while (!z.is_zero()) {
1758                         ex g = univariate_gcd(w, z, x);
1759                         res *= power(g, i);
1760                         w = quo(w, g, x);
1761                         y = quo(z, g, x);
1762                         z = y - w.diff(x);
1763                         i++;
1764                 }
1765         }
1766         return res * power(w, i);
1767 }
1768
1769
1770 /*
1771  *  Normal form of rational functions
1772  */
1773
1774 /*
1775  *  Note: The internal normal() functions (= basic::normal() and overloaded
1776  *  functions) all return lists of the form {numerator, denominator}. This
1777  *  is to get around mul::eval()'s automatic expansion of numeric coefficients.
1778  *  E.g. (a+b)/3 is automatically converted to a/3+b/3 but we want to keep
1779  *  the information that (a+b) is the numerator and 3 is the denominator.
1780  */
1781
1782 /** Create a symbol for replacing the expression "e" (or return a previously
1783  *  assigned symbol). The symbol is appended to sym_lst and returned, the
1784  *  expression is appended to repl_lst.
1785  *  @see ex::normal */
1786 static ex replace_with_symbol(const ex &e, lst &sym_lst, lst &repl_lst)
1787 {
1788         // Expression already in repl_lst? Then return the assigned symbol
1789         for (unsigned i=0; i<repl_lst.nops(); i++)
1790                 if (repl_lst.op(i).is_equal(e))
1791                         return sym_lst.op(i);
1792
1793         // Otherwise create new symbol and add to list, taking care that the
1794         // replacement expression doesn't contain symbols from the sym_lst
1795         // because subs() is not recursive
1796         symbol s;
1797         ex es(s);
1798         ex e_replaced = e.subs(sym_lst, repl_lst);
1799         sym_lst.append(es);
1800         repl_lst.append(e_replaced);
1801         return es;
1802 }
1803
1804 /** Create a symbol for replacing the expression "e" (or return a previously
1805  *  assigned symbol). An expression of the form "symbol == expression" is added
1806  *  to repl_lst and the symbol is returned.
1807  *  @see ex::to_rational */
1808 static ex replace_with_symbol(const ex &e, lst &repl_lst)
1809 {
1810         // Expression already in repl_lst? Then return the assigned symbol
1811         for (unsigned i=0; i<repl_lst.nops(); i++)
1812                 if (repl_lst.op(i).op(1).is_equal(e))
1813                         return repl_lst.op(i).op(0);
1814
1815         // Otherwise create new symbol and add to list, taking care that the
1816         // replacement expression doesn't contain symbols from the sym_lst
1817         // because subs() is not recursive
1818         symbol s;
1819         ex es(s);
1820         ex e_replaced = e.subs(repl_lst);
1821         repl_lst.append(es == e_replaced);
1822         return es;
1823 }
1824
1825 /** Default implementation of ex::normal(). It replaces the object with a
1826  *  temporary symbol.
1827  *  @see ex::normal */
1828 ex basic::normal(lst &sym_lst, lst &repl_lst, int level) const
1829 {
1830         return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
1831 }
1832
1833
1834 /** Implementation of ex::normal() for symbols. This returns the unmodified symbol.
1835  *  @see ex::normal */
1836 ex symbol::normal(lst &sym_lst, lst &repl_lst, int level) const
1837 {
1838         return (new lst(*this, _ex1()))->setflag(status_flags::dynallocated);
1839 }
1840
1841
1842 /** Implementation of ex::normal() for a numeric. It splits complex numbers
1843  *  into re+I*im and replaces I and non-rational real numbers with a temporary
1844  *  symbol.
1845  *  @see ex::normal */
1846 ex numeric::normal(lst &sym_lst, lst &repl_lst, int level) const
1847 {
1848         numeric num = numer();
1849         ex numex = num;
1850
1851         if (num.is_real()) {
1852                 if (!num.is_integer())
1853                         numex = replace_with_symbol(numex, sym_lst, repl_lst);
1854         } else { // complex
1855                 numeric re = num.real(), im = num.imag();
1856                 ex re_ex = re.is_rational() ? re : replace_with_symbol(re, sym_lst, repl_lst);
1857                 ex im_ex = im.is_rational() ? im : replace_with_symbol(im, sym_lst, repl_lst);
1858                 numex = re_ex + im_ex * replace_with_symbol(I, sym_lst, repl_lst);
1859         }
1860
1861         // Denominator is always a real integer (see numeric::denom())
1862         return (new lst(numex, denom()))->setflag(status_flags::dynallocated);
1863 }
1864
1865
1866 /** Fraction cancellation.
1867  *  @param n  numerator
1868  *  @param d  denominator
1869  *  @return cancelled fraction {n, d} as a list */
1870 static ex frac_cancel(const ex &n, const ex &d)
1871 {
1872         ex num = n;
1873         ex den = d;
1874         numeric pre_factor = _num1();
1875
1876 //std::clog << "frac_cancel num = " << num << ", den = " << den << endl;
1877
1878         // Handle trivial case where denominator is 1
1879         if (den.is_equal(_ex1()))
1880                 return (new lst(num, den))->setflag(status_flags::dynallocated);
1881
1882         // Handle special cases where numerator or denominator is 0
1883         if (num.is_zero())
1884                 return (new lst(num, _ex1()))->setflag(status_flags::dynallocated);
1885         if (den.expand().is_zero())
1886                 throw(std::overflow_error("frac_cancel: division by zero in frac_cancel"));
1887
1888         // Bring numerator and denominator to Z[X] by multiplying with
1889         // LCM of all coefficients' denominators
1890         numeric num_lcm = lcm_of_coefficients_denominators(num);
1891         numeric den_lcm = lcm_of_coefficients_denominators(den);
1892         num = multiply_lcm(num, num_lcm);
1893         den = multiply_lcm(den, den_lcm);
1894         pre_factor = den_lcm / num_lcm;
1895
1896         // Cancel GCD from numerator and denominator
1897         ex cnum, cden;
1898         if (gcd(num, den, &cnum, &cden, false) != _ex1()) {
1899                 num = cnum;
1900                 den = cden;
1901         }
1902
1903         // Make denominator unit normal (i.e. coefficient of first symbol
1904         // as defined by get_first_symbol() is made positive)
1905         const symbol *x;
1906         if (get_first_symbol(den, x)) {
1907                 GINAC_ASSERT(is_ex_exactly_of_type(den.unit(*x),numeric));
1908                 if (ex_to_numeric(den.unit(*x)).is_negative()) {
1909                         num *= _ex_1();
1910                         den *= _ex_1();
1911                 }
1912         }
1913
1914         // Return result as list
1915 //std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << endl;
1916         return (new lst(num * pre_factor.numer(), den * pre_factor.denom()))->setflag(status_flags::dynallocated);
1917 }
1918
1919
1920 /** Implementation of ex::normal() for a sum. It expands terms and performs
1921  *  fractional addition.
1922  *  @see ex::normal */
1923 ex add::normal(lst &sym_lst, lst &repl_lst, int level) const
1924 {
1925         if (level == 1)
1926                 return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
1927         else if (level == -max_recursion_level)
1928                 throw(std::runtime_error("max recursion level reached"));
1929
1930         // Normalize children and split each one into numerator and denominator
1931         exvector nums, dens;
1932         nums.reserve(seq.size()+1);
1933         dens.reserve(seq.size()+1);
1934         epvector::const_iterator it = seq.begin(), itend = seq.end();
1935         while (it != itend) {
1936                 ex n = recombine_pair_to_ex(*it).bp->normal(sym_lst, repl_lst, level-1);
1937                 nums.push_back(n.op(0));
1938                 dens.push_back(n.op(1));
1939                 it++;
1940         }
1941         ex n = overall_coeff.bp->normal(sym_lst, repl_lst, level-1);
1942         nums.push_back(n.op(0));
1943         dens.push_back(n.op(1));
1944         GINAC_ASSERT(nums.size() == dens.size());
1945
1946         // Now, nums is a vector of all numerators and dens is a vector of
1947         // all denominators
1948 //std::clog << "add::normal uses " << nums.size() << " summands:\n";
1949
1950         // Add fractions sequentially
1951         exvector::const_iterator num_it = nums.begin(), num_itend = nums.end();
1952         exvector::const_iterator den_it = dens.begin(), den_itend = dens.end();
1953 //std::clog << " num = " << *num_it << ", den = " << *den_it << endl;
1954         ex num = *num_it++, den = *den_it++;
1955         while (num_it != num_itend) {
1956 //std::clog << " num = " << *num_it << ", den = " << *den_it << endl;
1957                 ex next_num = *num_it++, next_den = *den_it++;
1958
1959                 // Trivially add sequences of fractions with identical denominators
1960                 while ((den_it != den_itend) && next_den.is_equal(*den_it)) {
1961                         next_num += *num_it;
1962                         num_it++; den_it++;
1963                 }
1964
1965                 // Additiion of two fractions, taking advantage of the fact that
1966                 // the heuristic GCD algorithm computes the cofactors at no extra cost
1967                 ex co_den1, co_den2;
1968                 ex g = gcd(den, next_den, &co_den1, &co_den2, false);
1969                 num = ((num * co_den2) + (next_num * co_den1)).expand();
1970                 den *= co_den2;         // this is the lcm(den, next_den)
1971         }
1972 //std::clog << " common denominator = " << den << endl;
1973
1974         // Cancel common factors from num/den
1975         return frac_cancel(num, den);
1976 }
1977
1978
1979 /** Implementation of ex::normal() for a product. It cancels common factors
1980  *  from fractions.
1981  *  @see ex::normal() */
1982 ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const
1983 {
1984         if (level == 1)
1985                 return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
1986         else if (level == -max_recursion_level)
1987                 throw(std::runtime_error("max recursion level reached"));
1988
1989         // Normalize children, separate into numerator and denominator
1990         ex num = _ex1();
1991         ex den = _ex1();
1992         ex n;
1993         epvector::const_iterator it = seq.begin(), itend = seq.end();
1994         while (it != itend) {
1995                 n = recombine_pair_to_ex(*it).bp->normal(sym_lst, repl_lst, level-1);
1996                 num *= n.op(0);
1997                 den *= n.op(1);
1998                 it++;
1999         }
2000         n = overall_coeff.bp->normal(sym_lst, repl_lst, level-1);
2001         num *= n.op(0);
2002         den *= n.op(1);
2003
2004         // Perform fraction cancellation
2005         return frac_cancel(num, den);
2006 }
2007
2008
2009 /** Implementation of ex::normal() for powers. It normalizes the basis,
2010  *  distributes integer exponents to numerator and denominator, and replaces
2011  *  non-integer powers by temporary symbols.
2012  *  @see ex::normal */
2013 ex power::normal(lst &sym_lst, lst &repl_lst, int level) const
2014 {
2015         if (level == 1)
2016                 return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
2017         else if (level == -max_recursion_level)
2018                 throw(std::runtime_error("max recursion level reached"));
2019
2020         // Normalize basis and exponent (exponent gets reassembled)
2021         ex n_basis = basis.bp->normal(sym_lst, repl_lst, level-1);
2022         ex n_exponent = exponent.bp->normal(sym_lst, repl_lst, level-1);
2023         n_exponent = n_exponent.op(0) / n_exponent.op(1);
2024
2025         if (n_exponent.info(info_flags::integer)) {
2026
2027                 if (n_exponent.info(info_flags::positive)) {
2028
2029                         // (a/b)^n -> {a^n, b^n}
2030                         return (new lst(power(n_basis.op(0), n_exponent), power(n_basis.op(1), n_exponent)))->setflag(status_flags::dynallocated);
2031
2032                 } else if (n_exponent.info(info_flags::negative)) {
2033
2034                         // (a/b)^-n -> {b^n, a^n}
2035                         return (new lst(power(n_basis.op(1), -n_exponent), power(n_basis.op(0), -n_exponent)))->setflag(status_flags::dynallocated);
2036                 }
2037
2038         } else {
2039
2040                 if (n_exponent.info(info_flags::positive)) {
2041
2042                         // (a/b)^x -> {sym((a/b)^x), 1}
2043                         return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
2044
2045                 } else if (n_exponent.info(info_flags::negative)) {
2046
2047                         if (n_basis.op(1).is_equal(_ex1())) {
2048
2049                                 // a^-x -> {1, sym(a^x)}
2050                                 return (new lst(_ex1(), replace_with_symbol(power(n_basis.op(0), -n_exponent), sym_lst, repl_lst)))->setflag(status_flags::dynallocated);
2051
2052                         } else {
2053
2054                                 // (a/b)^-x -> {sym((b/a)^x), 1}
2055                                 return (new lst(replace_with_symbol(power(n_basis.op(1) / n_basis.op(0), -n_exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
2056                         }
2057
2058                 } else {        // n_exponent not numeric
2059
2060                         // (a/b)^x -> {sym((a/b)^x, 1}
2061                         return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
2062                 }
2063         }
2064 }
2065
2066
2067 /** Implementation of ex::normal() for pseries. It normalizes each coefficient and
2068  *  replaces the series by a temporary symbol.
2069  *  @see ex::normal */
2070 ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const
2071 {
2072         epvector newseq;
2073         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
2074                 ex restexp = i->rest.normal();
2075                 if (!restexp.is_zero())
2076                         newseq.push_back(expair(restexp, i->coeff));
2077         }
2078         ex n = pseries(relational(var,point), newseq);
2079         return (new lst(replace_with_symbol(n, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
2080 }
2081
2082
2083 /** Implementation of ex::normal() for relationals. It normalizes both sides.
2084  *  @see ex::normal */
2085 ex relational::normal(lst &sym_lst, lst &repl_lst, int level) const
2086 {
2087         return (new lst(relational(lh.normal(), rh.normal(), o), _ex1()))->setflag(status_flags::dynallocated);
2088 }
2089
2090
2091 /** Normalization of rational functions.
2092  *  This function converts an expression to its normal form
2093  *  "numerator/denominator", where numerator and denominator are (relatively
2094  *  prime) polynomials. Any subexpressions which are not rational functions
2095  *  (like non-rational numbers, non-integer powers or functions like sin(),
2096  *  cos() etc.) are replaced by temporary symbols which are re-substituted by
2097  *  the (normalized) subexpressions before normal() returns (this way, any
2098  *  expression can be treated as a rational function). normal() is applied
2099  *  recursively to arguments of functions etc.
2100  *
2101  *  @param level maximum depth of recursion
2102  *  @return normalized expression */
2103 ex ex::normal(int level) const
2104 {
2105         lst sym_lst, repl_lst;
2106
2107         ex e = bp->normal(sym_lst, repl_lst, level);
2108         GINAC_ASSERT(is_ex_of_type(e, lst));
2109
2110         // Re-insert replaced symbols
2111         if (sym_lst.nops() > 0)
2112                 e = e.subs(sym_lst, repl_lst);
2113
2114         // Convert {numerator, denominator} form back to fraction
2115         return e.op(0) / e.op(1);
2116 }
2117
2118 /** Numerator of an expression. If the expression is not of the normal form
2119  *  "numerator/denominator", it is first converted to this form and then the
2120  *  numerator is returned.
2121  *
2122  *  @see ex::normal
2123  *  @return numerator */
2124 ex ex::numer(void) const
2125 {
2126         lst sym_lst, repl_lst;
2127
2128         ex e = bp->normal(sym_lst, repl_lst, 0);
2129         GINAC_ASSERT(is_ex_of_type(e, lst));
2130
2131         // Re-insert replaced symbols
2132         if (sym_lst.nops() > 0)
2133                 return e.op(0).subs(sym_lst, repl_lst);
2134         else
2135                 return e.op(0);
2136 }
2137
2138 /** Denominator of an expression. If the expression is not of the normal form
2139  *  "numerator/denominator", it is first converted to this form and then the
2140  *  denominator is returned.
2141  *
2142  *  @see ex::normal
2143  *  @return denominator */
2144 ex ex::denom(void) const
2145 {
2146         lst sym_lst, repl_lst;
2147
2148         ex e = bp->normal(sym_lst, repl_lst, 0);
2149         GINAC_ASSERT(is_ex_of_type(e, lst));
2150
2151         // Re-insert replaced symbols
2152         if (sym_lst.nops() > 0)
2153                 return e.op(1).subs(sym_lst, repl_lst);
2154         else
2155                 return e.op(1);
2156 }
2157
2158
2159 /** Default implementation of ex::to_rational(). It replaces the object with a
2160  *  temporary symbol.
2161  *  @see ex::to_rational */
2162 ex basic::to_rational(lst &repl_lst) const
2163 {
2164         return replace_with_symbol(*this, repl_lst);
2165 }
2166
2167
2168 /** Implementation of ex::to_rational() for symbols. This returns the
2169  *  unmodified symbol.
2170  *  @see ex::to_rational */
2171 ex symbol::to_rational(lst &repl_lst) const
2172 {
2173         return *this;
2174 }
2175
2176
2177 /** Implementation of ex::to_rational() for a numeric. It splits complex
2178  *  numbers into re+I*im and replaces I and non-rational real numbers with a
2179  *  temporary symbol.
2180  *  @see ex::to_rational */
2181 ex numeric::to_rational(lst &repl_lst) const
2182 {
2183         if (is_real()) {
2184                 if (!is_rational())
2185                         return replace_with_symbol(*this, repl_lst);
2186         } else { // complex
2187                 numeric re = real();
2188                 numeric im = imag();
2189                 ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl_lst);
2190                 ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl_lst);
2191                 return re_ex + im_ex * replace_with_symbol(I, repl_lst);
2192         }
2193         return *this;
2194 }
2195
2196
2197 /** Implementation of ex::to_rational() for powers. It replaces non-integer
2198  *  powers by temporary symbols.
2199  *  @see ex::to_rational */
2200 ex power::to_rational(lst &repl_lst) const
2201 {
2202         if (exponent.info(info_flags::integer))
2203                 return power(basis.to_rational(repl_lst), exponent);
2204         else
2205                 return replace_with_symbol(*this, repl_lst);
2206 }
2207
2208
2209 /** Implementation of ex::to_rational() for expairseqs.
2210  *  @see ex::to_rational */
2211 ex expairseq::to_rational(lst &repl_lst) const
2212 {
2213         epvector s;
2214         s.reserve(seq.size());
2215         for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
2216                 s.push_back(split_ex_to_pair(recombine_pair_to_ex(*it).to_rational(repl_lst)));
2217                 // s.push_back(combine_ex_with_coeff_to_pair((*it).rest.to_rational(repl_lst),
2218         }
2219         ex oc = overall_coeff.to_rational(repl_lst);
2220         if (oc.info(info_flags::numeric))
2221                 return thisexpairseq(s, overall_coeff);
2222         else s.push_back(combine_ex_with_coeff_to_pair(oc,_ex1()));
2223         return thisexpairseq(s, default_overall_coeff());
2224 }
2225
2226
2227 /** Rationalization of non-rational functions.
2228  *  This function converts a general expression to a rational polynomial
2229  *  by replacing all non-rational subexpressions (like non-rational numbers,
2230  *  non-integer powers or functions like sin(), cos() etc.) to temporary
2231  *  symbols. This makes it possible to use functions like gcd() and divide()
2232  *  on non-rational functions by applying to_rational() on the arguments,
2233  *  calling the desired function and re-substituting the temporary symbols
2234  *  in the result. To make the last step possible, all temporary symbols and
2235  *  their associated expressions are collected in the list specified by the
2236  *  repl_lst parameter in the form {symbol == expression}, ready to be passed
2237  *  as an argument to ex::subs().
2238  *
2239  *  @param repl_lst collects a list of all temporary symbols and their replacements
2240  *  @return rationalized expression */
2241 ex ex::to_rational(lst &repl_lst) const
2242 {
2243         return bp->to_rational(repl_lst);
2244 }
2245
2246
2247 #ifndef NO_NAMESPACE_GINAC
2248 } // namespace GiNaC
2249 #endif // ndef NO_NAMESPACE_GINAC