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