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