- See if __GNUC__ < 2.97 before using std::vector<..,malloc_alloc>. Sorry,
[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 "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 special cases where numerator or denominator is 0
1870         if (num.is_zero())
1871                 return (new lst(_ex0(), _ex1()))->setflag(status_flags::dynallocated);
1872         if (den.expand().is_zero())
1873                 throw(std::overflow_error("frac_cancel: division by zero in frac_cancel"));
1874
1875         // Bring numerator and denominator to Z[X] by multiplying with
1876         // LCM of all coefficients' denominators
1877         numeric num_lcm = lcm_of_coefficients_denominators(num);
1878         numeric den_lcm = lcm_of_coefficients_denominators(den);
1879         num = multiply_lcm(num, num_lcm);
1880         den = multiply_lcm(den, den_lcm);
1881         pre_factor = den_lcm / num_lcm;
1882
1883         // Cancel GCD from numerator and denominator
1884         ex cnum, cden;
1885         if (gcd(num, den, &cnum, &cden, false) != _ex1()) {
1886                 num = cnum;
1887                 den = cden;
1888         }
1889
1890         // Make denominator unit normal (i.e. coefficient of first symbol
1891         // as defined by get_first_symbol() is made positive)
1892         const symbol *x;
1893         if (get_first_symbol(den, x)) {
1894                 GINAC_ASSERT(is_ex_exactly_of_type(den.unit(*x),numeric));
1895                 if (ex_to_numeric(den.unit(*x)).is_negative()) {
1896                         num *= _ex_1();
1897                         den *= _ex_1();
1898                 }
1899         }
1900
1901         // Return result as list
1902 //std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << endl;
1903         return (new lst(num * pre_factor.numer(), den * pre_factor.denom()))->setflag(status_flags::dynallocated);
1904 }
1905
1906
1907 /** Implementation of ex::normal() for a sum. It expands terms and performs
1908  *  fractional addition.
1909  *  @see ex::normal */
1910 ex add::normal(lst &sym_lst, lst &repl_lst, int level) const
1911 {
1912         if (level == 1)
1913                 return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
1914         else if (level == -max_recursion_level)
1915                 throw(std::runtime_error("max recursion level reached"));
1916
1917         // Normalize and expand children, chop into summands
1918         exvector o;
1919         o.reserve(seq.size()+1);
1920         epvector::const_iterator it = seq.begin(), itend = seq.end();
1921         while (it != itend) {
1922
1923                 // Normalize and expand child
1924                 ex n = recombine_pair_to_ex(*it).bp->normal(sym_lst, repl_lst, level-1).expand();
1925
1926                 // If numerator is a sum, chop into summands
1927                 if (is_ex_exactly_of_type(n.op(0), add)) {
1928                         epvector::const_iterator bit = ex_to_add(n.op(0)).seq.begin(), bitend = ex_to_add(n.op(0)).seq.end();
1929                         while (bit != bitend) {
1930                                 o.push_back((new lst(recombine_pair_to_ex(*bit), n.op(1)))->setflag(status_flags::dynallocated));
1931                                 bit++;
1932                         }
1933
1934                         // The overall_coeff is already normalized (== rational), we just
1935                         // split it into numerator and denominator
1936                         GINAC_ASSERT(ex_to_numeric(ex_to_add(n.op(0)).overall_coeff).is_rational());
1937                         numeric overall = ex_to_numeric(ex_to_add(n.op(0)).overall_coeff);
1938                         o.push_back((new lst(overall.numer(), overall.denom() * n.op(1)))->setflag(status_flags::dynallocated));
1939                 } else
1940                         o.push_back(n);
1941                 it++;
1942         }
1943         o.push_back(overall_coeff.bp->normal(sym_lst, repl_lst, level-1));
1944
1945         // o is now a vector of {numerator, denominator} lists
1946
1947         // Determine common denominator
1948         ex den = _ex1();
1949         exvector::const_iterator ait = o.begin(), aitend = o.end();
1950 //std::clog << "add::normal uses the following summands:\n";
1951         while (ait != aitend) {
1952 //std::clog << " num = " << ait->op(0) << ", den = " << ait->op(1) << endl;
1953                 den = lcm(ait->op(1), den, false);
1954                 ait++;
1955         }
1956 //std::clog << " common denominator = " << den << endl;
1957
1958         // Add fractions
1959         if (den.is_equal(_ex1())) {
1960
1961                 // Common denominator is 1, simply add all fractions
1962                 exvector num_seq;
1963                 for (ait=o.begin(); ait!=aitend; ait++) {
1964                         num_seq.push_back(ait->op(0) / ait->op(1));
1965                 }
1966                 return (new lst((new add(num_seq))->setflag(status_flags::dynallocated), den))->setflag(status_flags::dynallocated);
1967
1968         } else {
1969
1970                 // Perform fractional addition
1971                 exvector num_seq;
1972                 for (ait=o.begin(); ait!=aitend; ait++) {
1973                         ex q;
1974                         if (!divide(den, ait->op(1), q, false)) {
1975                                 // should not happen
1976                                 throw(std::runtime_error("invalid expression in add::normal, division failed"));
1977                         }
1978                         num_seq.push_back((ait->op(0) * q).expand());
1979                 }
1980                 ex num = (new add(num_seq))->setflag(status_flags::dynallocated);
1981
1982                 // Cancel common factors from num/den
1983                 return frac_cancel(num, den);
1984         }
1985 }
1986
1987
1988 /** Implementation of ex::normal() for a product. It cancels common factors
1989  *  from fractions.
1990  *  @see ex::normal() */
1991 ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const
1992 {
1993         if (level == 1)
1994                 return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
1995         else if (level == -max_recursion_level)
1996                 throw(std::runtime_error("max recursion level reached"));
1997
1998         // Normalize children, separate into numerator and denominator
1999         ex num = _ex1();
2000         ex den = _ex1(); 
2001         ex n;
2002         epvector::const_iterator it = seq.begin(), itend = seq.end();
2003         while (it != itend) {
2004                 n = recombine_pair_to_ex(*it).bp->normal(sym_lst, repl_lst, level-1);
2005                 num *= n.op(0);
2006                 den *= n.op(1);
2007                 it++;
2008         }
2009         n = overall_coeff.bp->normal(sym_lst, repl_lst, level-1);
2010         num *= n.op(0);
2011         den *= n.op(1);
2012
2013         // Perform fraction cancellation
2014         return frac_cancel(num, den);
2015 }
2016
2017
2018 /** Implementation of ex::normal() for powers. It normalizes the basis,
2019  *  distributes integer exponents to numerator and denominator, and replaces
2020  *  non-integer powers by temporary symbols.
2021  *  @see ex::normal */
2022 ex power::normal(lst &sym_lst, lst &repl_lst, int level) const
2023 {
2024         if (level == 1)
2025                 return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
2026         else if (level == -max_recursion_level)
2027                 throw(std::runtime_error("max recursion level reached"));
2028
2029         // Normalize basis
2030         ex n = basis.bp->normal(sym_lst, repl_lst, level-1);
2031
2032         if (exponent.info(info_flags::integer)) {
2033
2034                 if (exponent.info(info_flags::positive)) {
2035
2036                         // (a/b)^n -> {a^n, b^n}
2037                         return (new lst(power(n.op(0), exponent), power(n.op(1), exponent)))->setflag(status_flags::dynallocated);
2038
2039                 } else if (exponent.info(info_flags::negative)) {
2040
2041                         // (a/b)^-n -> {b^n, a^n}
2042                         return (new lst(power(n.op(1), -exponent), power(n.op(0), -exponent)))->setflag(status_flags::dynallocated);
2043                 }
2044
2045         } else {
2046
2047                 if (exponent.info(info_flags::positive)) {
2048
2049                         // (a/b)^x -> {sym((a/b)^x), 1}
2050                         return (new lst(replace_with_symbol(power(n.op(0) / n.op(1), exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
2051
2052                 } else if (exponent.info(info_flags::negative)) {
2053
2054                         if (n.op(1).is_equal(_ex1())) {
2055
2056                                 // a^-x -> {1, sym(a^x)}
2057                                 return (new lst(_ex1(), replace_with_symbol(power(n.op(0), -exponent), sym_lst, repl_lst)))->setflag(status_flags::dynallocated);
2058
2059                         } else {
2060
2061                                 // (a/b)^-x -> {sym((b/a)^x), 1}
2062                                 return (new lst(replace_with_symbol(power(n.op(1) / n.op(0), -exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
2063                         }
2064
2065                 } else {        // exponent not numeric
2066
2067                         // (a/b)^x -> {sym((a/b)^x, 1}
2068                         return (new lst(replace_with_symbol(power(n.op(0) / n.op(1), exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
2069                 }
2070         }
2071 }
2072
2073
2074 /** Implementation of ex::normal() for pseries. It normalizes each coefficient and
2075  *  replaces the series by a temporary symbol.
2076  *  @see ex::normal */
2077 ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const
2078 {
2079         epvector new_seq;
2080         new_seq.reserve(seq.size());
2081
2082         epvector::const_iterator it = seq.begin(), itend = seq.end();
2083         while (it != itend) {
2084                 new_seq.push_back(expair(it->rest.normal(), it->coeff));
2085                 it++;
2086         }
2087         ex n = pseries(relational(var,point), new_seq);
2088         return (new lst(replace_with_symbol(n, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
2089 }
2090
2091
2092 /** Implementation of ex::normal() for relationals. It normalizes both sides.
2093  *  @see ex::normal */
2094 ex relational::normal(lst &sym_lst, lst &repl_lst, int level) const
2095 {
2096         return (new lst(relational(lh.normal(), rh.normal(), o), _ex1()))->setflag(status_flags::dynallocated);
2097 }
2098
2099
2100 /** Normalization of rational functions.
2101  *  This function converts an expression to its normal form
2102  *  "numerator/denominator", where numerator and denominator are (relatively
2103  *  prime) polynomials. Any subexpressions which are not rational functions
2104  *  (like non-rational numbers, non-integer powers or functions like sin(),
2105  *  cos() etc.) are replaced by temporary symbols which are re-substituted by
2106  *  the (normalized) subexpressions before normal() returns (this way, any
2107  *  expression can be treated as a rational function). normal() is applied
2108  *  recursively to arguments of functions etc.
2109  *
2110  *  @param level maximum depth of recursion
2111  *  @return normalized expression */
2112 ex ex::normal(int level) const
2113 {
2114         lst sym_lst, repl_lst;
2115
2116         ex e = bp->normal(sym_lst, repl_lst, level);
2117         GINAC_ASSERT(is_ex_of_type(e, lst));
2118
2119         // Re-insert replaced symbols
2120         if (sym_lst.nops() > 0)
2121                 e = e.subs(sym_lst, repl_lst);
2122
2123         // Convert {numerator, denominator} form back to fraction
2124         return e.op(0) / e.op(1);
2125 }
2126
2127 /** Numerator of an expression. If the expression is not of the normal form
2128  *  "numerator/denominator", it is first converted to this form and then the
2129  *  numerator is returned.
2130  *
2131  *  @see ex::normal
2132  *  @return numerator */
2133 ex ex::numer(void) const
2134 {
2135         lst sym_lst, repl_lst;
2136
2137         ex e = bp->normal(sym_lst, repl_lst, 0);
2138         GINAC_ASSERT(is_ex_of_type(e, lst));
2139
2140         // Re-insert replaced symbols
2141         if (sym_lst.nops() > 0)
2142                 return e.op(0).subs(sym_lst, repl_lst);
2143         else
2144                 return e.op(0);
2145 }
2146
2147 /** Denominator of an expression. If the expression is not of the normal form
2148  *  "numerator/denominator", it is first converted to this form and then the
2149  *  denominator is returned.
2150  *
2151  *  @see ex::normal
2152  *  @return denominator */
2153 ex ex::denom(void) const
2154 {
2155         lst sym_lst, repl_lst;
2156
2157         ex e = bp->normal(sym_lst, repl_lst, 0);
2158         GINAC_ASSERT(is_ex_of_type(e, lst));
2159
2160         // Re-insert replaced symbols
2161         if (sym_lst.nops() > 0)
2162                 return e.op(1).subs(sym_lst, repl_lst);
2163         else
2164                 return e.op(1);
2165 }
2166
2167
2168 /** Default implementation of ex::to_rational(). It replaces the object with a
2169  *  temporary symbol.
2170  *  @see ex::to_rational */
2171 ex basic::to_rational(lst &repl_lst) const
2172 {
2173         return replace_with_symbol(*this, repl_lst);
2174 }
2175
2176
2177 /** Implementation of ex::to_rational() for symbols. This returns the
2178  *  unmodified symbol.
2179  *  @see ex::to_rational */
2180 ex symbol::to_rational(lst &repl_lst) const
2181 {
2182         return *this;
2183 }
2184
2185
2186 /** Implementation of ex::to_rational() for a numeric. It splits complex
2187  *  numbers into re+I*im and replaces I and non-rational real numbers with a
2188  *  temporary symbol.
2189  *  @see ex::to_rational */
2190 ex numeric::to_rational(lst &repl_lst) const
2191 {
2192         if (is_real()) {
2193                 if (!is_rational())
2194                         return replace_with_symbol(*this, repl_lst);
2195         } else { // complex
2196                 numeric re = real();
2197                 numeric im = imag();
2198                 ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl_lst);
2199                 ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl_lst);
2200                 return re_ex + im_ex * replace_with_symbol(I, repl_lst);
2201         }
2202         return *this;
2203 }
2204
2205
2206 /** Implementation of ex::to_rational() for powers. It replaces non-integer
2207  *  powers by temporary symbols.
2208  *  @see ex::to_rational */
2209 ex power::to_rational(lst &repl_lst) const
2210 {
2211         if (exponent.info(info_flags::integer))
2212                 return power(basis.to_rational(repl_lst), exponent);
2213         else
2214                 return replace_with_symbol(*this, repl_lst);
2215 }
2216
2217
2218 /** Implementation of ex::to_rational() for expairseqs.
2219  *  @see ex::to_rational */
2220 ex expairseq::to_rational(lst &repl_lst) const
2221 {
2222         epvector s;
2223         s.reserve(seq.size());
2224         for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
2225                 s.push_back(split_ex_to_pair(recombine_pair_to_ex(*it).to_rational(repl_lst)));
2226                 // s.push_back(combine_ex_with_coeff_to_pair((*it).rest.to_rational(repl_lst),
2227         }
2228         ex oc = overall_coeff.to_rational(repl_lst);
2229         if (oc.info(info_flags::numeric))
2230                 return thisexpairseq(s, overall_coeff);
2231         else s.push_back(combine_ex_with_coeff_to_pair(oc,_ex1()));
2232         return thisexpairseq(s, default_overall_coeff());
2233 }
2234
2235
2236 /** Rationalization of non-rational functions.
2237  *  This function converts a general expression to a rational polynomial
2238  *  by replacing all non-rational subexpressions (like non-rational numbers,
2239  *  non-integer powers or functions like sin(), cos() etc.) to temporary
2240  *  symbols. This makes it possible to use functions like gcd() and divide()
2241  *  on non-rational functions by applying to_rational() on the arguments,
2242  *  calling the desired function and re-substituting the temporary symbols
2243  *  in the result. To make the last step possible, all temporary symbols and
2244  *  their associated expressions are collected in the list specified by the
2245  *  repl_lst parameter in the form {symbol == expression}, ready to be passed
2246  *  as an argument to ex::subs().
2247  *
2248  *  @param repl_lst collects a list of all temporary symbols and their replacements
2249  *  @return rationalized expression */
2250 ex ex::to_rational(lst &repl_lst) const
2251 {
2252         return bp->to_rational(repl_lst);
2253 }
2254
2255
2256 #ifndef NO_NAMESPACE_GINAC
2257 } // namespace GiNaC
2258 #endif // ndef NO_NAMESPACE_GINAC