when there are multiple variables with the same maximum degree, the one with
[ginac.git] / 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